the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Estimating gradients of physical fields in space
Yufei Zhou
This study focuses on the development of a multipoint technique for future constellation missions, aiming to measure gradients at various orders, in particular the linear and quadratic gradients, of a general field. It is well established that, in order to estimate linear gradients, the spacecraft must not lie on a plane. Through analytical exploration within the framework of least squares, it is demonstrated that at least 10 spacecraft that do not lie on any quadric surface are required to estimate both linear and quadratic gradients. The spatial arrangement of the spacecraft can be characterized by a set of quality factors. In cases where there is poor temporal synchronization among the spacecraft leading to nonsimultaneous measurements, temporal gradients must be included. If the spacecraft have multiple velocities, by incorporating temporal gradients it is possible to reduce the number of required spacecraft. Furthermore, it is proved that the accuracy of the linear gradient is of second order and that of the quadratic gradient is of first order. Additionally, a method for estimating errors in the calculation is also illustrated.
 Article
(1260 KB)  Fulltext XML
 BibTeX
 EndNote
Multipoint measurements have significantly advanced our understanding of the structures and dynamics of space plasmas. The basic approach involves direct interpretation of the collected data. However, to maximize the potential of these measurements, several techniques have been developed to estimate additional quantities that would otherwise remain inaccessible. One common initial step is to estimate the linear gradients of physical fields, with a particular focus on the magnetic field (Chanteur, 1998; De Keyser, 2008; De Keyser et al., 2007; Dunlop et al., 1988; Harvey, 1998; Hamrin et al., 2008; Vogt et al., 2008, 2009; Shen and Dunlop, 2023). These gradients serve various purposes, such as calculating the electric current density (Dunlop et al., 2015, 2016, 2018), determining the curvature and rotation rate of magnetic field lines (Shen et al., 2003, 2007), locating magnetic nulls crucial for magnetic reconnection (Fu et al., 2015), and determining the dimensionality and velocity of magnetic structures (Shi et al., 2005, 2006; Fadanelli et al., 2019). A recent technique utilized the gradients of normal fields on curved boundary layers to estimate the principal curvatures and directions of boundary layers (Shen et al., 2020; Shao et al., 2023; Zhou et al., 2023).
The recent MMS (Magnetospheric Multiscale) mission has improved particle data measurements with exceptional resolution. With this capability, the electric current can be directly calculated by summing the product of the bulk flux and charge of particles (Burch et al., 2015; Pollock et al., 2016). By leveraging Maxwell's equations and incorporating additional information, such as the electric current measurements from each spacecraft, it becomes possible to estimate not only the linear gradients, but also the quadratic gradients of magnetic fields using fourpoint measurement (Liu et al., 2019; Torbert et al., 2020; Denton et al., 2020; Shen et al., 2021a), though the general estimation of these gradients of an arbitrary field typically requires 10 spacecraft measurements (Chanteur, 1998; Shen et al., 2021c). With both linear and quadratic gradients known, the complete geometry of the magnetic field lines, including their curvature and torsion, can be obtained (Shen et al., 2021a; Torbert et al., 2020). This is of particular use in the reconstruction of key regions such as in reconnection. Unlike other reconstruction methods (see e.g. Sonnerup and Teh, 2008; Hasegawa et al., 2021), the approach utilizing gradients avoids assumptions specific to the reconnection process, thus making it adaptable to a wide range of conditions, though it is still limited by the separation of the spacecraft such that the result may not be as accurate when studying those phenomena of much smaller or larger scales.
At present, there are a growing tendency of enhanced resolution in particle and electric field measurement and an increased number of spacecraft involved in multispacecraft missions (Ogilvie et al., 1977; Escoubet et al., 2001; Liu et al., 2005; FriisChristensen et al., 2006; Angelopoulos, 2008; Burch et al., 2015; Spence et al., 2022; Maruca et al., 2021). An algorithm for the linear and quadratic gradients (ALQG) has been developed that relies on 10 or more measurement points to tackle the general problem of estimating quadratic gradients of physical fields that are not limited to magnetic fields alone (Shen et al., 2021c). In ALQG, the quadratic gradients can be obtained by solving a matrix equation. The characteristic matrix, ℜ^{MN}, which is determined by the positions of the spacecraft within the constellations, has been put forward. If the determinant of the characteristic matrix ℜ^{MN} is nonzero, the full quadratic gradients can be obtained. One application is the measurement of electric charge density using the Poisson equation (Shen et al., 2021c, b). In this approach, the charge density is calculated by summing the diagonal elements of the estimated quadratic gradient matrix of a potential field (Shen et al., 2021b).
However, despite progress in addressing some of the associated challenges, several issues remain unresolved. The first problem revolves around the relationship between the feasibility of estimation and the distribution of measurement points. It is well established that, in fourpoint measurements, linear gradients can be obtained as long as the points do not lie on a plane (Vogt et al., 2009; Shen et al., 2012; Shen and Dunlop, 2023). However, the impact of the point distribution on the estimation of quadratic gradients has not been fully understood. This poses a challenge in determining the optimal distribution that ensures accurate estimation. When four spacecraft are on a plane, it is still possible to obtain the linear gradients on the plane (Vogt et al., 2009; Shen et al., 2012; Shen and Dunlop, 2023). When dealing with quadratic gradients, if a distribution of measurement points is found to be unsuitable for achieving a complete estimation, there is no method available to extract the utmost information regarding the gradients.
The second problem concerns the requirement of simultaneity in measurements, which applies to both the new technique for quadratic gradients and previous techniques for linear gradients (Harvey, 1998; Chanteur, 1998; Hamrin et al., 2008). As the number of spacecraft increases, the issue of temporal synchronization among them becomes more pronounced. One possible approach to mitigate this problem is to incorporate temporal gradients into the analysis (De Keyser et al., 2007; De Keyser, 2008).
The third problem pertains to the accuracy of the estimation process and the associated errors. Although the technique has demonstrated high accuracy when tested on synthetic data, with suggestions that errors in linear gradients are of second order and errors in quadratic gradients are of first order (Shen et al., 2021c), these results have not been deduced analytically. In practical applications, measurements include noise which may also affect estimates of gradients. It is therefore crucial to develop a reliable method for estimating and quantifying errors of various origins.
In this regard, this study presents a further development of ALQG. In addition to calculating quadratic gradients, the results can be applied to reconstruct physical fields and structures in space using polynomials.
We start with the problem of approximation. To approximate a vector field, one approach is to aggregate the approximations of its individual component fields, treating each component as an independent scalar field. This method is useful when there is no additional information available regarding the relationship between the component fields, such as the constraint $\mathrm{\nabla}\cdot \mathit{B}=\mathrm{0}$. For simplicity we consider here the problem of approximating a scalar field in space, and the result can be applied equally well to vector fields.
A field can be seen as a combination of multiple constituent fields originating from different sources. These fields often have distinct temporal and spatial scales. For instance, in the inner magnetosphere during a substorm, the total magnetic field comprises the dipole geomagnetic field (and higherorder moments), disturbances caused by currents (Yang et al., 2012), and other localized and temporary variations. On the bow shock front, various waves superimpose. In most cases, our focus is on specific constituents, such as the disturbance field during a substorm or the shock ramp on a shock front. Therefore, we can express the total field j(x) as the sum of a background field of interest f(x) and wave fields w(x) with smaller scales compared to f(x):
Here, x is an rcomponent vector representing a point. The general case is when r=4 and $\mathit{x}=[vt,x,y,z]=[{x}_{\mathrm{0}},{x}_{\mathrm{1}},{x}_{\mathrm{2}},{x}_{\mathrm{3}}]$, which represents a point in time and space. v is a constant to scale the temporal coordinate with the spatial ones. It can be chosen as the characteristic speed in the system, such as the mean Alfvén speed or flow speed in the region of concern. The choice of v does not impact the general method described below. The scaling of the spatial coordinates is discussed in Appendix B. It is also possible to consider the field in a cut of time and space, that is, at a specific time. Then, r=3 and $\mathit{x}=[x,y,z]=[{x}_{\mathrm{1}},{x}_{\mathrm{2}},{x}_{\mathrm{3}}]$ represents a point in space. Since our objective is to approximate f(x), we can represent it using multiindex notation (see Appendix A) as a sum of multivariate polynomials:
In this equation, g_{α} is the coefficient of the polynomial x^{α}, and we employ the properties of multiindex notation (Properties 3 and 4 in Appendix A). By comparing Eq. (2) with the Taylor expansion of f(x) around $\stackrel{\mathrm{\u203e}}{\mathrm{0}}=[\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{0}]$, we can see that the coefficients g_{α} are related to the gradients ${f}_{,\mathit{\alpha}}\left(\stackrel{\mathrm{\u203e}}{\mathrm{0}}\right)$ as follows:
where we employ Property 8 of the multiindex notation. Suppose we aim to approximate f(x) using polynomials up to degree d. We define
By doing so, we separate the summation in Eq. (2) into a polynomial of degree at most d, denoted as p_{d}(x), and a polynomial in which all the terms have degrees higher than d, denoted as ${p}_{d}^{+}\left(\mathit{x}\right)$. There are $\left({\scriptscriptstyle \begin{array}{c}d+r\\ r\end{array}}\right)=(d+r)\mathrm{!}/r\mathrm{!}d\mathrm{!}$ terms in Eq. (4), resulting in an equal number of coefficients to be determined from measured data. Now we can rewrite Eq. (1) as
When field measurements are conducted using probes, we need to consider the positioning error in time and space, denoted as $\mathit{\delta}\mathit{x}=[v\mathit{\delta}t,\mathit{\delta}x,\mathit{\delta}y,\mathit{\delta}z]$. Suppose we think the total field is measured at x_{m} but due to the positioning error it is actually measured at x_{m}+δx. Taking into account the measurement error in the field, denoted as δj, we can express the sampled data j_{m} as follows:
Note that the scaling factor v for the temporal coordinate x_{m0}=vt_{m} is the same for all the measurement points.
Consider M measurements taken at different points in time and space, yielding data pairs j_{m} and x_{m} for $\mathrm{1}\le m\le M$. The objective is to determine a set of numerical values for g_{α}, where $\left\mathit{\alpha}\right\le d$, that yield the best approximation of f(x) by p_{d}(x) based on these data. It is evident from Eq. (7) that the discarded polynomial ${p}_{d}^{+}\left(\mathit{x}\right)$, the wave field w(x), the measurement error δj, and the positioning error δx all contribute to the final error when solving this problem. The concept of measurement by probes in space is illustrated in Fig. 1.
We define the error between the measured field and the approximating polynomial as s_{m}, given by
To quantify the total error, we employ the weighted leastsquares method, which constructs the total error as a weighted sum of all individual errors:
Here, the weight matrix W_{mn} determines the contribution of each measurement to the approximation. The choice of the weight matrix depends on the specific problem (De Keyser et al., 2007), but in a simple case where all measurements are equally important, it can be expressed as
Generally, it is symmetric and invertible. Minimizing the total error with respect to g_{β} and assuming that this is done when ${g}_{\mathit{\beta}}={\stackrel{\mathrm{\u0303}}{g}}_{\mathit{\beta}}$ results in a set of $\left({\scriptscriptstyle \begin{array}{c}d+r\\ r\end{array}}\right)$ equations for ${\stackrel{\mathrm{\u0303}}{g}}_{\mathit{\beta}}$:
We define the matrix R with elements:
Additionally, we define
With these notations, taking into account Eqs. (9), (8), (7), and (4), Eq. (11) can be explicitly expressed as a system of equations:
This linear system of equations consists of $\left({\scriptscriptstyle \begin{array}{c}d+r\\ r\end{array}}\right)$ equations and unknowns. The tilde notation on g_{α} signifies that it represents an estimated quantity rather than the true value.
The solution to Eq. (14), i.e. the estimation ${\stackrel{\mathrm{\u0303}}{g}}_{\mathit{\alpha}}$, can be obtained directly using common computer programs designed to solve linear systems. By applying the relation in Eq. (3), the gradients up to the dth degree of the field at the origin $\stackrel{\mathrm{\u203e}}{\mathrm{0}}$ can be determined. The approximation of the field f(x) is then given by
It is important to note that, at this stage, the coordinate system, specifically its origin, has not been chosen. In Sect. 5, we will demonstrate that the centre of the measurement points, when chosen as the origin, yields the best reduction of the approximation error resulting from the truncation of the Taylor series.
4.1 The requirement of a unique solution
From Eq. (14) there exists a unique set of solutions for ${\stackrel{\mathrm{\u0303}}{g}}_{\mathit{\alpha}}$ if and only if R has full rank. This requirement has several implications regarding the number, distribution, and velocity of probes in space. To see these, we need to decompose R.
Based on the decomposition of the symmetric and invertible weight matrix as ${W}_{mn}={\sum}_{l,s,k}^{M}({P}^{T}{)}_{ml}{O}_{ls}{O}_{sk}{P}_{kn}$, where O is a diagonal matrix whose elements, when squared, are the eigenvalues of the weight matrix and P is composed of eigenvectors, we can express the matrix R as
Considering the relation rank(A^{T}A)=rank(A) and the invertibility of OP, we have
where the matrix X is defined by
Therefore, the uniqueness of the solution in Eq. (14) is equivalent to the rank of X being $\left({\scriptscriptstyle \begin{array}{c}d+r\\ r\end{array}}\right)$.
The matrix X has rows corresponding to different measurement points and columns corresponding to coefficients g_{α}. To achieve a rank of $\left({\scriptscriptstyle \begin{array}{c}d+r\\ r\end{array}}\right)$ for X, two conditions need to be met. First, the number of measurement points M should be at least $\left({\scriptscriptstyle \begin{array}{c}d+r\\ r\end{array}}\right)$. Second, the points should not all lie on an algebraic surface of degree at most d, ensuring that there is no set of coefficients a_{α} such that
A similar result has also been obtained for multivariate interpolations (Olver, 2006).
Although we present the first condition separately from the second to stress its utility in application, it is contained in the second since a lack of measurement points necessarily makes the existing points lie on a surface prescribed by the second. For example, three points (d=1, r=3) must be on a plane and nine points (d=2, r=3) must be on a secondorder surface.
To illustrate the second condition, we take d=2, r=3 as an example. The matrix X in this case is given by
If all the points lie on a secondorder algebraic surface, we can express the surface formally with appropriately chosen coefficients a_{α} as
and all the points satisfy this equation. This indicates that we can make a linear combination of the columns in Eq. (20) with the coefficients in Eq. 21 and obtain a column of zeros. Thus, the rank of X is lower than $\left({\scriptscriptstyle \begin{array}{c}\mathrm{2}+\mathrm{3}\\ \mathrm{3}\end{array}}\right)$, which is the number of columns it possesses. On the other hand, if the points do not lie on a secondorder algebraic surface, then no set of coefficients exists to linearly combine the columns to reach a column of zeros. In this case the rank is $\left({\scriptscriptstyle \begin{array}{c}\mathrm{2}+\mathrm{3}\\ \mathrm{3}\end{array}}\right)$.
These two conditions have important implications for the orbit design of future multispacecraft missions and for adaptation of this general framework to specific problems in practice, such as measuring electric charges (Shen et al., 2021b) and reconstructing magnetic structures (Liu et al., 2019; Torbert et al., 2020; Shen et al., 2021a). Here we discuss them in detail and illustrate them with Fig. 2.
We first consider simultaneous measurements and r=3. If d=1, which is to estimate the spatial linear gradients, we recover the wellknown restriction that at least four measurement points are needed, and these points should not lie on a firstorder algebraic surface, or in other words a plane (Fig. 2c), such that with appropriately chosen coefficients a_{α} they satisfy
If d=2, in which case both the linear and quadratic gradients are to be estimated, we need at least 10 measurement points. They should not reside on a secondorder algebraic surface, which can be defined by Eq. (21) (Fig. 2f). Typical examples of a secondorder surface include an ellipsoid, an elliptic cone, an elliptic cylinder, or an elliptic paraboloid. Among them the sphere is common, as for the distribution of the probes to date. The geomagnetic stations are on the surface of solid Earth. The Iridium satellite constellation is in the ionosphere.
Next we consider r=4, and those time series data are incorporated to estimate the gradients of fields in time and space. If d=1, at least five points are needed, and they should not lie on a hypersurface in time and space. These five points can be obtained from four spacecraft moving with one velocity, as suggested by previous studies (De Keyser et al., 2007; De Keyser, 2008). If there are only three spacecraft available with identical velocities, the resulting measurement points will inevitably lie in a hyperplane in time and space. Alternatively, if the three spacecraft have at least two kinds of velocities, the measurement points can deviate from a plane and the gradients can be estimated. In the case of d=2, at least 15 points are required, and they should also not belong to a quadratic hypersurface. In this case, 10 spacecraft flying in formation suffice. If there are only nine spacecraft, then at least two velocities are needed.
To prove these necessitates a heavy symbolic computation best suited for a computer to handle effectively. Here we demonstrate the effect brought about by different velocities of spacecraft by considering the y axis in Fig. 2b as a temporal axis. If there is only one spacecraft flying with a constant velocity, its trajectory in time and space forms a straight line, violating the requirement that measurement points should not lie on a linear path. To render the three nonlinearly distributed measurement points, in addition to those by two spacecraft with the same velocity that result in parallel trajectories, it is also possible by a single spacecraft with changing velocity that leads to a nonlinear trajectory. Similarly, treating the z axis in Fig. 2c as a temporal one, the four measurement points can result from three spacecraft with the same velocity, viz. passed by three parallel lines. They can also be generated by two spacecraft with different velocities, each following a trajectory connecting a pair of points.
4.2 When the requirement is not met
In practice, there are situations where the requirement is not met. For d=1, this can occur due to instrument failures in a fourspacecraft mission or a lack of spacecraft to form a tetrahedron, resulting in only three spacecraft providing data that lie on a plane. Even in wellfunctioning fourspacecraft missions, orbital constraints can cause the spacecraft to be nearly coplanar at times. For d=2, many current probes are distributed spherically, such as geomagnetic stations on solid Earth or the Iridium satellite constellation in the ionosphere. The upcoming HelioSwarm mission will consist of only nine spacecraft. In future missions involving 10 or more spacecraft, the same challenges faced by fourspacecraft missions can also arise. Hence, it is crucial to explore whether there exists a method to effectively leverage the available data in such cases.
The direct problem is that Eq. (14) has an infinite number of solutions as the determinant of R becomes zero. One potential approach is to address this problem by excluding certain gradient components from the approximating polynomial (Eq. 4) and relocating them to the truncated one (Eq. 5). By doing so, the degrees of freedom in the problem can be adjusted to match the available measured data. However, it is important to note that randomly dropping components is not suitable, as Sect. 5 will demonstrate that this can lead to overwhelming errors and deterioration of all the estimated g_{α} ($\left\mathit{\alpha}\right=d$). To effectively reduce the degrees of freedom in the approximation, it is necessary to consider the degrees of freedom in the distribution of measurement points, specifically the rank of X, and the surfaces that contain the measurement points.
Suppose that there exist N sets of coefficients a_{α} such that Eq. (19) is satisfied, which indicates that the measurement points lie on the intersection of N distinct surfaces of degrees at most d and that the rank of X is $\left({\scriptscriptstyle \begin{array}{c}d+r\\ r\end{array}}\right)N$. Take $d=\mathrm{1},r=\mathrm{3}$ for example. If N=1 (N=2), then all the points lie on a surface (line) and the rank of X is 3 (2). By rightmultiplying X by a fullrank square matrix G, it is possible to obtain a matrix X^{′} whose last N columns are zeros. To put it more visually, it is possible to have
such that
Each of the last N columns of G is a set of coefficients a_{α} that represents a surface containing the measurement points. Since the process to obtain the G is quite involved and does not affect the scheme to calculate gradients, we will defer the discussion until the scheme is fully revealed. Equation (23) can be interpreted as a singular value decomposition in the form of X=ASB, where S is a $M\times \left({\scriptscriptstyle \begin{array}{c}d+r\\ r\end{array}}\right)$ diagonal matrix and A and B are unitary matrices of order M and $\left({\scriptscriptstyle \begin{array}{c}d+r\\ r\end{array}}\right)$ respectively (Kincaid and Cheney, 2002). If the transformation matrix G in Eq. (23) is chosen to be unitary, we have $\mathbf{G}={\mathbf{B}}^{\mathrm{1}}$ and ${\mathbf{X}}^{\prime}=\mathbf{AS}$.
By leftmultiplying Eq. (14) by G^{T}, making use of Eqs. (12) and (13) and considering the decomposition $\mathbf{I}={\mathbf{GG}}^{\mathrm{1}}$ where I is the identity matrix, we obtain
where j and $\stackrel{\mathrm{\u0303}}{\mathit{g}}$ are column vectors. ${\mathbf{G}}^{\mathrm{1}}\stackrel{\mathrm{\u0303}}{\mathit{g}}$ represents a recombination of the gradient components according to the distribution of the measurement points. In matrix form, this equation is written as
where J^{′} and R^{′} contain nonvanishing components and ${\stackrel{\mathrm{\u0303}}{\mathit{g}}}^{\prime}$ includes the first $\left({\scriptscriptstyle \begin{array}{c}d+r\\ r\end{array}}\right)N$ components of ${\mathbf{G}}^{\mathrm{1}}\stackrel{\mathrm{\u0303}}{\mathit{g}}$. Since Eqs. (23) and (26) can also be obtained by singular value decomposition, the geometrical interpretation of the existence condition of a unique solution given in this study is also equivalent to the one in terms of singular values (De Keyser et al., 2007). With Eq. (26) we have separated from the last N insoluble components of ${\mathbf{G}}^{\mathrm{1}}\stackrel{\mathrm{\u0303}}{\mathit{g}}$ the soluble ${\stackrel{\mathrm{\u0303}}{\mathit{g}}}^{\prime}$. By solving for them from
we can extract the maximum amount of information about gradients from the measurement points of a given distribution.
Let us illustrate this method with a simple example when d=1, r=3, and all the points satisfy z=0. In this case, N=1 and X itself is in the form of X^{′}. Thus, an identity matrix can be used in place of G to give a set of unknowns, ${\mathbf{G}}^{\mathrm{1}}\stackrel{\mathrm{\u0303}}{\mathit{g}}=\left[f\right(\stackrel{\mathrm{\u203e}}{\mathrm{0}}),{\partial}_{x}f(\stackrel{\mathrm{\u203e}}{\mathrm{0}}),{\partial}_{y}f(\stackrel{\mathrm{\u203e}}{\mathrm{0}}),{\partial}_{z}f(\stackrel{\mathrm{\u203e}}{\mathrm{0}}\left)\right]$, which suggests that the gradient along the z direction cannot be estimated, while the rest can still be obtained. This is intuitive in the case of estimating linear gradients. The problem has been addressed previously by the use of reciprocal vectors (Vogt et al., 2009). The benefit of the method here, however, comes from its general applicability in problems of all orders and for future missions that consist of more spacecraft.
Finally we discuss how to obtain G. The possible choices of G are infinite, since the form of X^{′} is invariant upon the linear recombination of the last N columns of G and the random replacement of the first $\left({\scriptscriptstyle \begin{array}{c}d+r\\ r\end{array}}\right)N$ columns as long as the resultant G has full rank. Among all the possible G, the most readily available one is the matrix of Gauss elimination, which we denote by G^{*}. To obtain this matrix, we perform Gauss column elimination on X so that the resultant X^{′} is triangular in its upper left. Each elementary column operation of the elimination is equivalent to rightmultiplication of an elementary matrix. The product of these elementary matrices is G^{*}.
To facilitate the error analysis in Sect. 5, we can also construct from the matrix of Gauss elimination a set of special G, which we denote by G^{′}. The last N columns of G^{′} are those of the G^{*}, ${G}_{lh}^{\prime}={G}_{lh}^{*}$ for $\mathrm{1}\le l\le \left({\scriptscriptstyle \begin{array}{c}d+r\\ r\end{array}}\right)$, and $\left(\genfrac{}{}{0ex}{}{d+r}{r}\right)N<h\le \left({\scriptscriptstyle \begin{array}{c}d+r\\ r\end{array}}\right)$. In the first $\left({\scriptscriptstyle \begin{array}{c}d+r\\ r\end{array}}\right)N$ columns, in addition to the rest being zeros, $\left({\scriptscriptstyle \begin{array}{c}d+r\\ r\end{array}}\right)N$ unities are so placed that the following two conditions are met.

G^{′} has full rank.

Let the row (column) index of unity be i (j). If $\left({\scriptscriptstyle \begin{array}{c}u\mathrm{1}+r\\ r\end{array}}\right)<i\le \left({\scriptscriptstyle \begin{array}{c}u+r\\ r\end{array}}\right)$ for some u and $\left({\scriptscriptstyle \begin{array}{c}v\mathrm{1}+r\\ r\end{array}}\right)<j\le \left({\scriptscriptstyle \begin{array}{c}v+r\\ r\end{array}}\right)$ for some v, then we should have u=v.
While a unique solution can be obtained for estimating ${\stackrel{\mathrm{\u0303}}{g}}_{\mathit{\alpha}}$ and ${\stackrel{\mathrm{\u0303}}{p}}_{d}\left(\mathit{x}\right)$, the accuracy may vary significantly due to various factors. One factor that influences the accuracy is the choice of the weight matrix W_{mn}. If prior information about the background field f and the wave field w is available, it is possible to adapt the weight matrix appropriately to improve the accuracy, as suggested by De Keyser (2008). For general purposes, the plain form of Eq. (10) is sufficient. This form provides a reasonable balance between simplicity and effectiveness in capturing the underlying field characteristics.
Let R^{−1} be the inverse of R. We multiply Eq. (14) by ${\left({R}^{\mathrm{1}}\right)}_{\mathit{\gamma}\mathit{\beta}}$, sum over β, and obtain
where use was made of Eqs. (13) and (7). According to the binomial theorem for multivariate polynomials (see Eq. A1), we have the decomposition
Substituting this into Eq. (28), subtracting g_{γ} from both sides, and defining the error in estimating g_{γ} as
we obtain the complete expression for the error:
The terms in the brackets on the left represent errors of various origins.
Here we consider the error caused by the truncation of Taylor series, i.e. the term containing ${p}_{d}^{+}({\mathit{x}}_{n}+\mathit{\delta}\mathit{x})$. Making use of Eq. (5), we express the relative truncation error in g_{γ} as
It is obvious that three factors combine to create this error. The first is the ratio of higherorder coefficients g_{α} to g_{γ}, which is inherent to the nature of the field being estimated. This ratio can be modelled by ${D}^{\left\mathit{\gamma}\right\left\mathit{\alpha}\right}$, where D is the scale of the field. The second is the values of the measurement points x_{m} which appear in both the inverse of R and the terms after the last summation sign. These values are determined by the choice of the origin and the size and configuration of the measurement points. The third is the positioning error in time and space. The MMS mission consists of four spacecraft that fly in close formation. The separation among them can reach 10 km at times with a relative position error of less than 100 m, i.e. 1 % of the separation. Since as compared to the differences in measurement points x_{m}, δx is usually small, we could ignore it here for the purpose of estimating truncation errors. Then we have the error as a sum of terms at various orders
where ${q}_{\mathit{\alpha}\mathit{\gamma}}^{\mathit{\#}}$ are dimensionless figures that can be calculated by comparing Eq. (33) with Eq. (32). # is used to indicate that ${q}_{\mathit{\alpha}\mathit{\gamma}}^{\mathit{\#}}$ has little physical meaning and will be replaced later.
It is then obvious that, to reduce the error, it is pertinent to choose the centre of measurement points as the origin, and so we have
Thus, Eq. (33) can be reexpressed as
where L is the characteristic dimension of the distribution of measurement points. L can be modelled by the square roots of the eigenvalues of the volumetric tensor (Harvey, 1998). The volumetric tensor R is defined by Eq. (12) when ${W}_{mn}={\mathit{\delta}}_{mn}/M$ and $\left\mathit{\alpha}\right=\left\mathit{\beta}\right=\mathrm{1}$. q_{αγ} are parameters to be calculated by comparing Eq. (35) and Eq. (32):
They are determined by the distribution of measurement points. For a given characteristic scale L of the points, through q_{αγ} the error of estimation can be affected by the distribution of points. Therefore, they can be termed quality factors that indicate whether or not the distribution is sound for the estimation. The absolute values of these quality factors vary from zero to infinity, with larger values representing better quality. In particular, the quality factors of $\left\mathit{\alpha}\right=d+\mathrm{1}$ are the most important since other quality factors correspond to higher orders of L.
In common cases, the accuracy of ${\stackrel{\mathrm{\u0303}}{g}}_{\mathit{\gamma}}$ is of the order of $d+\mathrm{1}\left\mathit{\gamma}\right$. For example, in estimating the gradients up to second order (d=2), the accuracy of the linear gradients is of second order and that of the quadratic gradients is of first order. This conclusion was also suggested by previous tests on synthetic data (Shen et al., 2021c). If the estimation is made up to third order (d=3), the accuracy of the linear gradient could reach third order and that of the quadratic gradient becomes second order.
The total error in Eq. (31) diminishes as the separation between measurement points is reduced, until the effects of errors in position and measurement and the effects of wave fields take hold, though the latter can be mitigated by lowpass filtering. For a specific field, the errors in position and measurement together with the wave field collectively set an upper limit on the accuracy achievable through the manipulation of measurement point configurations. Conversely, for a fixed set of measurement points, the accuracy depends on the comparative magnitudes of background variations versus those of the wave field and measurement errors. In practice, the total relative error can be computed from Eqs. (31) and (32) by modelling the wave field, position error, and measurement error as random variables and by using ${\stackrel{\mathrm{\u0303}}{g}}_{\mathit{\alpha}}$ in place of g_{α} for $\left\mathit{\alpha}\right\le d$ and ${D}^{\left\mathit{\gamma}\right\left\mathit{\alpha}\right}$ in place of ${g}_{\mathit{\alpha}}/{g}_{\mathit{\gamma}}$ for $\left\mathit{\alpha}\right>d$.
We now consider the error involved in the reconstruction of the field, i.e. by Eq. (15). The error is given by
Using Eqs. (5) and (33), we arrive at a trivial conclusion that the degree of the error is at least d+1.
When the degrees of freedom in the measured data are less than the required $\left({\scriptscriptstyle \begin{array}{c}d+r\\ r\end{array}}\right)$, the method described in Sect. 4.2 can be utilized. To analyse the involved error, we apply the foregoing procedures once more. We leftmultiply Eq. (27) by ${\mathbf{R}}^{\prime \mathrm{1}}$ to obtain
By defining
we have the following expression for it:
The relative error caused by truncation is given by
When the G^{′} presented in Sect. 4.2 are used for G, the elements in R^{′} and X^{′} are of the same order of L, as are the elements of R and X. Thus, the error can be expressed as
where the quality factor is given by
Therefore, this method designed for cases when measurement points are not well distributed has good accuracy.
The techniques for calculating linear gradients of general physical fields and quadratic gradients of magnetic fields using fourpoint measurements have been widely applied in the context of multispacecraft missions to advance our understanding of space plasma. However, there are also important quantities and processes associated with the quadratic gradients of other fields that warrant further exploration. For instance, the gradients of velocity play a crucial role in determining fundamental quantities such as viscosity and energy dissipation rate. Overall, the statics and dynamics of physical fields in space are interrelated through their gradients. As the number of spacecraft in a constellation continues to increase, it is helpful to explore and prepare for future missions multipoint techniques that rely on more points to estimate quadratic and higherorder gradients.
In summary, we have analytically explored the general method to estimate gradients of fields in space based on multipoint measurement. Regarding the feasibility of estimation, a general conclusion is that to estimate the complete gradients up to the dth degree using simultaneous measurement, $\left({\scriptscriptstyle \begin{array}{c}d+\mathrm{3}\\ \mathrm{3}\end{array}}\right)$ spacecraft are needed, and these spacecraft should not lie on a dthorder surface in space. In particular, at least 10 points that are not on a secondorder surface are needed to estimate both linear and quadratic gradients. To address the negative effects caused by poor synchronization among spacecraft in a large constellation and to estimate the additional temporal gradients of a field, a time series needs to be taken into account, and it is necessary to have at least $\left({\scriptscriptstyle \begin{array}{c}d+\mathrm{4}\\ \mathrm{4}\end{array}}\right)$ measurement points that do not lie on a dthorder hypersurface in time and space. For linear gradients, these measurement points can be provided by a constellation of four spacecraft with the same velocity or of three spacecraft whose velocities have at least two kinds. For quadratic gradients, 10 comoving spacecraft are sufficient. It is also possible to reduce one spacecraft by adding one more velocity. In situations where the measured data lack degrees of freedom due to ill configuration of spacecraft, which may include a shortage of spacecraft, it becomes necessary to invoke a transformation in order to estimate the gradient components to the best extent possible.
Regarding the accuracy, we have analytically proven that, in an estimation of gradients up to dth order, the order of accuracy of the athorder gradients is at least $d+\mathrm{1}a$. We have also provided quality factors q_{α} to judge the distribution of measurement points and the spacecraft configuration in a constellation. In addition, a method for estimating errors in real time has been presented.
The results obtained offer valuable insights into the development of multipoint techniques that rely on gradients of physical fields. Additionally, they hold significance for the future design of multispacecraft missions aimed at studying physics associated with quadratic or higherorder gradients.
The current study primarily focuses on approximating a single scalar physical field in space. It treats vectors and tensors as aggregations of multiple independent scalar fields. In practice, the constituents of a vector field can be interrelated, as exemplified by the divergencefree condition $\mathrm{\nabla}\cdot \mathit{B}=\mathrm{0}$ for a magnetic field. The gradients of different fields can also be subject to various physical formulas. For instance, the zerothorder gradients of electric currents and firstorder gradients of magnetic fields are correlated by $\mathrm{\nabla}\times \mathit{B}={\mathit{\mu}}_{\mathrm{0}}\mathit{j}$. Beyond these linear constraints, nonlinear constraints exist as well. The gradients of entropy s and velocity v of an isentropic flow satisfy ${\partial}_{\mathrm{t}}s+\mathit{v}\cdot \mathrm{\nabla}s=\mathrm{0}$. Incorporating all conceivable constraints into the current framework is challenging and will be explored in the future. However, when only a limited number of constraints such as the sole divergencefree condition for a magnetic field are taken into account, the existence condition remains unchanged for a complete solution concerning the configuration of measurement points in time and space.
From the numerical point of view, the matrix R (Eq. 12) is likely to be ill conditioned because it is the weighted product of two Vandermonde matrices. This together with the limited resolution of measurement puts a limitation on the practicability of the technique for approximating to higher orders by solving Eq. (14), though the framework is in principle applicable to all orders. Thus, for higherorder approximations, it is necessary to verify in Eq. (31) that the error resulting from the multiplication of δj by the terms outside the square brackets is not substantial. As for quadratic gradients, previous simulations have verified the feasibility, reliability, and accuracy of the technique (Shen et al., 2021c).
Here we list the properties of a multiindex notation tailored for multivariate functions (Riachy et al., 2011).
Let $\mathit{\alpha}=\left({\mathit{\alpha}}_{\mathrm{1}},\mathrm{\dots},{\mathit{\alpha}}_{r}\right)$ be an rtuple of nonnegative integers ${\mathit{\alpha}}_{i},i=\mathrm{1},\mathrm{\dots},r;i,r\in \mathbb{N}$. α is called a multiindex. The symbol in bold x denotes a vector in ℝ^{r}. As for time and space, r=4.
For multiindices $\mathit{\alpha},\mathit{\beta}\in {\mathbb{N}}^{r}$, the following properties are either defined or deduced.

Componentwise sum and difference: $\mathit{\alpha}\pm \mathit{\beta}=({\mathit{\alpha}}_{\mathrm{1}}\pm {\mathit{\beta}}_{\mathrm{1}},\mathrm{\dots},{\mathit{\alpha}}_{r}\pm {\mathit{\beta}}_{r})$.

Partial order $\mathit{\alpha}\le \mathit{\beta}\iff {\mathit{\alpha}}_{i}\le {\mathit{\beta}}_{i},\forall i\in \mathit{\{}\mathrm{1},\mathrm{\dots},r\mathit{\}}$. $\mathit{\alpha}=\mathit{\beta}\iff {\mathit{\alpha}}_{i}={\mathit{\beta}}_{i},\forall i\in \mathit{\{}\mathrm{1},\mathrm{\dots},r\mathit{\}}$.

Given $\mathit{x}=({x}_{\mathrm{1}},\mathrm{\dots},{x}_{r})\in {\mathbb{R}}^{r}$, we have ${\mathit{x}}^{\mathit{\alpha}}={x}_{\mathrm{1}}^{{\mathit{\alpha}}_{\mathrm{1}}}\mathrm{\cdots}{x}_{r}^{{\mathit{\alpha}}_{r}}$.

The total degree of x^{α} is given by $\left\mathit{\alpha}\right={\mathit{\alpha}}_{\mathrm{1}}+\mathrm{\dots}+{\mathit{\alpha}}_{r}$.

Factorial: α!=α_{1}!⋯α_{r}!.

Binomial coefficient: $\left({\scriptscriptstyle \begin{array}{c}\mathit{\alpha}\\ \mathit{\beta}\end{array}}\right)=\left({\scriptscriptstyle \begin{array}{c}{\mathit{\alpha}}_{\mathrm{1}}\\ {\mathit{\beta}}_{\mathrm{1}}\end{array}}\right)\mathrm{\cdots}\left({\scriptscriptstyle \begin{array}{c}{\mathit{\alpha}}_{r}\\ {\mathit{\beta}}_{r}\end{array}}\right)$.

$\stackrel{\mathrm{\u203e}}{b}=(b,\mathrm{\dots},b),b\in \mathbb{N},\stackrel{\mathrm{\u203e}}{b}\in {\mathbb{N}}^{r}$.

Higherorder partial derivative ${\partial}^{\mathit{\alpha}}\equiv {\partial}_{\mathrm{1}}^{{\mathit{\alpha}}_{\mathrm{1}}}\mathrm{\cdots}{\partial}_{r}^{{\mathit{\alpha}}_{r}}$, where ${\partial}_{i}^{{\mathit{\alpha}}_{i}}\equiv \frac{{\partial}^{{\mathit{\alpha}}_{i}}}{\partial {x}_{i}^{{\mathit{\alpha}}_{i}}}$. ${\partial}^{\mathit{\alpha}}f\equiv {f}_{,\mathit{\alpha}}$.

Denote by 1_{i}∈ℕ^{r} the multiindex with zeros for all elements except the ith one, i.e. ${\mathrm{1}}_{i}=\left(\mathrm{0},\mathrm{\dots},\mathrm{0},\mathrm{1},\mathrm{0},\mathrm{\dots},\mathrm{0}\right)$.

The tensor product of two vectors $\mathit{u},\mathit{v}\in {\mathbb{R}}^{r}$ is defined by $\mathit{u}\otimes \mathit{v}=\left({u}_{\mathrm{1}}\mathit{v},\mathrm{\dots},{u}_{r}\mathit{v}\right)\in {\mathbb{R}}^{{r}^{\mathrm{2}}}$.

Binomial theorem:
$$\begin{array}{}\text{(A1)}& {\left(\mathit{x}+\mathit{y}\right)}^{\mathit{\alpha}}=\sum _{\stackrel{\mathrm{\u203e}}{\mathrm{0}}\le \mathit{\beta}\le \mathit{\alpha}}\left(\begin{array}{c}\mathit{\alpha}\\ \mathit{\beta}\end{array}\right){\mathit{x}}^{\mathit{\beta}}{\mathit{y}}^{\mathit{\alpha}\mathit{\beta}}.\end{array}$$
In addition to the scaling of a temporal coordinate by a characteristic speed to obtain x_{0}=vt, it has been suggested that scaling on the three spatial coordinates can be invoked to further improve accuracy when the spatial variations of the physical fields are highly anisotropic (De Keyser et al., 2007). A recent observational study (Liu et al., 2022) showed that the ratios of the characteristic scales parallel to magnetic fields over perpendicular scales are roughly 2:1 for solar wind and magnetosheath plasmas. A corresponding scaling as such can be applied to the spatial coordinates of measurement points before calculating the matrix R and the vector J and solving for the gradients from Eq. (14).
YZ conceived of the study and performed the investigation, with contributions from CS. YZ wrote the paper. CS reviewed and edited the paper.
The contact author has declared that neither of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
This work was supported by the National Natural Science Foundation of China (grant no. 42130202) and the National Key Research and Development Program of China (grant no. 2022YFA1604600).
This research has been supported by the National Natural Science Foundation of China (grant no. 42130202) and the National Key Research and Development Program of China (grant no. 2022YFA1604600).
This paper was edited by Oliver Allanson and reviewed by Johan De Keyser and one anonymous referee.
Angelopoulos, V.: The THEMIS Mission, Space Sci. Rev., 141, 5–34, https://doi.org/10.1007/s1121400893361, 2008. a
Burch, J. L., Moore, T. E., Torbert, R. B., and Giles, B. L.: Magnetospheric Multiscale Overview and Science Objectives, Space Sci. Rev., 199, 5–21, https://doi.org/10.1007/s1121401501649, 2015. a, b
Chanteur, G.: Spatial Interpolation for Four Spacecraft: Theory, in: Analysis Methods for MultiSpacecraft Data, edited by: Paschmann, G. and Daly, P. W., p. 349, ESA Publications Division, Noordwijk, the Netherlands, ISBN: 1608280X, 1998. a, b, c
De Keyser, J.: Leastsquares multispacecraft gradient calculation with automatic error estimation, Ann. Geophys., 26, 3295–3316, https://doi.org/10.5194/angeo2632952008, 2008. a, b, c, d
De Keyser, J., Darrouzet, F., Dunlop, M. W., and Décréau, P. M. E.: Leastsquares gradient calculation from multipoint observations of scalar and vector fields: methodology and applications with Cluster in the plasmasphere, Ann. Geophys., 25, 971–987, https://doi.org/10.5194/angeo259712007, 2007. a, b, c, d, e, f
Denton, R. E., Torbert, R. B., Hasegawa, H., Dors, I., Genestreti, K. J., Argall, M. R., Gershman, D., Contel, O. L., Burch, J. L., Russell, C. T., Strangeway, R. J., Giles, B. L., and Fischer, D.: Polynomial Reconstruction of the Reconnection Magnetic Field Observed by Multiple Spacecraft, J. Geophys. Res.Space, 125, e2019JA027481, https://doi.org/10.1029/2019JA027481, 2020. a
Dunlop, M., Southwood, D., Glassmeier, K.H., and Neubauer, F.: Analysis of multipoint magnetometer data, Adv. Space Res., 8, 273–277, https://doi.org/10.1016/02731177(88)90141X, 1988. a
Dunlop, M. W., Yang, Y.Y., Yang, J.Y., Lühr, H., Shen, C., Olsen, N., Ritter, P., Zhang, Q.H., Cao, J.B., Fu, H.S., and Haagmans, R.: Multispacecraft current estimates at swarm, J. Geophys. Res.Space, 120, 8307–8316, https://doi.org/10.1002/2015JA021707, 2015. a
Dunlop, M. W., Haaland, S., Escoubet, P. C., and Dong, X.C.: Commentary on accessing 3D currents in space: Experiences from Cluster, J. Geophys. Res.Space, 121, 7881–7886, https://doi.org/10.1002/2016JA022668, 2016. a
Dunlop, M. W., Haaland, S., Dong, X.C., Middleton, H. R., Escoubet, C. P., Yang, Y.Y., Zhang, Q.H., Shi, J.K., and Russell, C. T.: Multipoint Analysis of Electric Currents in Geospace Using the Curlometer Technique, in: Electric Currents in Geospace and Beyond, 67–80, John Wiley & Sons, Inc., https://doi.org/10.1002/9781119324522.ch4, 2018. a
Escoubet, C. P., Fehringer, M., and Goldstein, M.: Introduction: The Cluster mission, Ann. Geophys., 19, 1197–1200, https://doi.org/10.5194/angeo1911972001, 2001. a
Fadanelli, S., Lavraud, B., Califano, F., Jacquey, C., Vernisse, Y., Kacem, I., Penou, E., Gershman, D. J., Dorelli, J., Pollock, C., Giles, B. L., Avanov, L. A., Burch, J., Chandler, M. O., Coffey, V. N., Eastwood, J. P., Ergun, R., Farrugia, C. J., Fuselier, S. A., Genot, V. N., Grigorenko, E., Hasegawa, H., Khotyaintsev, Y., Contel, O. L., Marchaudon, A., Moore, T. E., Nakamura, R., Paterson, W. R., Phan, T., Rager, A. C., Russell, C. T., Saito, Y., Sauvaud, J.A., Schiff, C., Smith, S. E., Redondo, S. T., Torbert, R. B., Wang, S., and Yokota, S.: FourSpacecraft Measurements of the Shape and Dimensionality of Magnetic Structures in the NearEarth Plasma Environment, J. Geophys. Res.Space, 124, 6850–6868, https://doi.org/10.1029/2019JA026747, 2019. a
FriisChristensen, E., Lühr, H., and Hulot, G.: Swarm: A constellation to study the Earth's magnetic field, Earth Planet. Space, 58, 351–358, https://doi.org/10.1186/BF03351933, 2006. a
Fu, H. S., Vaivads, A., Khotyaintsev, Y. V., Olshevsky, V., André, M., Cao, J. B., Huang, S. Y., Retinò, A., and Lapenta, G.: How to find magnetic nulls and reconstruct field topology with MMS data?, J. Geophys. Res.Space, 120, 3758–3782, https://doi.org/10.1002/2015JA021082, 2015. a
Hamrin, M., Rönnmark, K., Börlin, N., Vedin, J., and Vaivads, A.: GALSGradient Analysis by Least Squares, Ann. Geophys., 26, 3491–3499, https://doi.org/10.5194/angeo2634912008, 2008. a, b
Harvey, C. C.: Spatial gradients and the volumetric tensor, in: Analysis Methods for MultiSpacecraft Data, edited by: Paschmann, G. and Daly, P. W., p. 307, ESA Publications Division, Noordwijk, the Netherlands, ISBN: 1608280X, 1998. a, b, c
Hasegawa, H., Nakamura, T. K. M., and Denton, R. E.: Reconstruction of the Electron Diffusion Region With Inertia and Compressibility Effects, J. Geophys. Res.Space, 126, e2021JA029841, https://doi.org/10.1029/2021JA029841, 2021. a
Kincaid, D. and Cheney, W.: Numerical Analysis Mathematics of Scientific Computing, American Mathematical Society, 3rd Edn., ISBN: 9780821847886, 2002. a
Liu, Y. Y., Fu, H. S., Olshevsky, V., Pontin, D. I., Liu, C. M., Wang, Z., Chen, G., Dai, L., and Retino, A.: SOTE: A Nonlinear Method for Magnetic Topology Reconstruction in Space Plasmas, Astrophys. J. Suppl. Ser., 244, p. 31, https://doi.org/10.3847/15384365/ab391a, 2019. a, b
Liu, Y. Y., Wang, Z., Chen, G., Yu, Y., Guo, Z. Z., and Xiong, X.: Testing the Linearity of Vector Fields in Cold and Dense Space Plasmas, Astrophys. J., 929, 155, https://doi.org/10.3847/15384357/ac5d4b, 2022. a
Liu, Z. X., Escoubet, C. P., Pu, Z., Laakso, H., Shi, J. K., Shen, C., and Hapgood, M.: The Double Star mission, Ann. Geophys., 23, 2707–2712, https://doi.org/10.5194/angeo2327072005, 2005. a
Maruca, B. A., Rueda, J. A. A., Bandyopadhyay, R., Bianco, F. B., Chasapis, A., Chhiber, R., DeWeese, H., Matthaeus, W. H., Miles, D. M., Qudsi, R. A., Richardson, M. J., Servidio, S., Shay, M. A., Sundkvist, D., Verscharen, D., Vines, S. K., Westlake, J. H., and Wicks, R. T.: MagneToRE: Mapping the 3D Magnetic Structure of the Solar Wind Using a Large Constellation of Nanosatellites, Front. Astron. Space Sci., 8, 665885, https://doi.org/10.3389/fspas.2021.665885, 2021. a
Ogilvie, K. W., von Rosenvinge, T., and Durney, A. C.: International SunEarth Explorer: A ThreeSpacecraft Program, Science, 198, 131–138, https://doi.org/10.1126/science.198.4313.131, 1977. a
Olver, P. J.: On Multivariate Interpolation, Stud. Appl. Math., 116, 201–240, https://doi.org/10.1111/j.14679590.2006.00335.x, 2006. a
Pollock, C., Moore, T., Jacques, A., Burch, J., Gliese, U., Saito, Y., Omoto, T., Avanov, L., Barrie, A., Coffey, V., Dorelli, J., Gershman, D., Giles, B., Rosnack, T., Salo, C., Yokota, S., Adrian, M., Aoustin, C., Auletti, C., Aung, S., Bigio, V., Cao, N., Chandler, M., Chornay, D., Christian, K., Clark, G., Collinson, G., Corris, T., Santos, A. D. L., Devlin, R., Diaz, T., Dickerson, T., Dickson, C., Diekmann, A., Diggs, F., Duncan, C., FigueroaVinas, A., Firman, C., Freeman, M., Galassi, N., Garcia, K., Goodhart, G., Guererro, D., Hageman, J., Hanley, J., Hemminger, E., Holland, M., Hutchins, M., James, T., Jones, W., Kreisler, S., Kujawski, J., Lavu, V., Lobell, J., LeCompte, E., Lukemire, A., MacDonald, E., Mariano, A., Mukai, T., Narayanan, K., Nguyan, Q., Onizuka, M., Paterson, W., Persyn, S., Piepgrass, B., Cheney, F., Rager, A., Raghuram, T., Ramil, A., Reichenthal, L., Rodriguez, H., Rouzaud, J., Rucker, A., Saito, Y., Samara, M., Sauvaud, J.A., Schuster, D., Shappirio, M., Shelton, K., Sher, D., Smith, D., Smith, K., Smith, S., Steinfeld, D., Szymkiewicz, R., Tanimoto, K., Taylor, J., Tucker, C., Tull, K., Uhl, A., Vloet, J., Walpole, P., Weidner, S., White, D., Winkert, G., Yeh, P.S., and Zeuch, M.: Fast Plasma Investigation for Magnetospheric Multiscale, Space Sci. Rev., 199, 331–406, https://doi.org/10.1007/s1121401602454, 2016. a
Riachy, S., Mboup, M., and Richard, J.P.: Multivariate numerical differentiation, J. Comput. Appl. Math., 236, 1069–1089, https://doi.org/10.1016/j.cam.2011.07.031, 2011. a
Shao, P., Shen, C., Ma, Y., Rong, Z., Zhou, Y., Zhang, C., Dunlop, M., and Ji, Y.: Flapping motion configurations of geomagnetotail current sheet, J. Atmos. Sol.Terr. Phys., 243, 106019, https://doi.org/10.1016/j.jastp.2023.106019, 2023. a
Shen, C. and Dunlop, M.: Field gradient analysis based on a geometrical approach, J. Geophys. Res.Space, 128, e2023JA031313, https://doi.org/10.1029/2023JA031313, 2023. a, b, c
Shen, C., Li, X., Dunlop, M., Liu, Z. X., Balogh, A., Baker, D. N., Hapgood, M., and Wang, X.: Analyses on the geometrical structure of magnetic field in the current sheet based on cluster measurements, J. Geophys. Res.Space, 108, 1168, https://doi.org/10.1029/2002JA009612, 2003. a
Shen, C., Li, X., Dunlop, M., Shi, Q. Q., Liu, Z. X., Lucek, E., and Chen, Z. Q.: Magnetic field rotation analysis and the applications, J. Geophys. Res.Space, 112, A06211, https://doi.org/10.1029/2005JA011584, 2007. a
Shen, C., Rong, Z. J., Dunlop, M. W., Ma, Y. H., Li, X., Zeng, G., Yan, G. Q., Wan, W. X., Liu, Z. X., Carr, C. M., and Rème, H.: Spatial gradients from irregular, multiplepoint spacecraft configurations, J. Geophys. Res.Space, 117, A11207, https://doi.org/10.1029/2012JA018075, 2012. a, b
Shen, C., Zeng, G., Zhang, C., Rong, Z., Dunlop, M., Russell, C. T., Escoubet, C. P., and Ren, N.: Determination of the Configurations of Boundaries in Space, J. Geophys. Res.Space, 125, e2020JA028163, https://doi.org/10.1029/2020JA028163, 2020. a
Shen, C., Zhang, C., Rong, Z., Pu, Z., Dunlop, M. W., Escoubet, C. P., Russell, C. T., Zeng, G., Ren, N., Burch, J. L., and Zhou, Y.: Nonlinear Magnetic Gradients and Complete Magnetic Geometry From Multispacecraft Measurements, J. Geophys. Res.Space, 126, e2020JA028846, https://doi.org/10.1029/2020JA028846, 2021a. a, b, c
Shen, C., Zhou, Y., Gao, L., Wang, X., Pu, Z., Escoubet, C. P., and Burch, J. L.: Measurements of the Net Charge Density of Space Plasmas, J. Geophys. Res.Space, 126, e2021JA029511, https://doi.org/10.1029/2021JA029511, 2021b. a, b, c
Shen, C., Zhou, Y., Ma, Y., Wang, X., Pu, Z., and Dunlop, M.: A General Algorithm for the Linear and Quadratic Gradients of Physical Quantities Based on 10 or More Point Measurements, J. Geophys. Res.Space, 126, e2021JA029121, https://doi.org/10.1029/2021JA029121, 2021c. a, b, c, d, e, f
Shi, Q. Q., Shen, C., Pu, Z. Y., Dunlop, M. W., Zong, Q.G., Zhang, H., Xiao, C. J., Liu, Z. X., and Balogh, A.: Dimensional analysis of observed structures using multipoint magnetic field measurements: Application to Cluster, Geophys. Res. Lett., 32, L12105, https://doi.org/10.1029/2005GL022454, 2005. a
Shi, Q. Q., Shen, C., Dunlop, M. W., Pu, Z. Y., Zong, Q.G., Liu, Z. X., Lucek, E., and Balogh, A.: Motion of observed structures calculated from multipoint magnetic field measurements: Application to Cluster, Geophys. Res. Lett., 33, L08109, https://doi.org/10.1029/2005GL025073, 2006. a
Sonnerup, B. U. O. and Teh, W.L.: Reconstruction of twodimensional coherent MHD structures in a space plasma: The theory, J. Geophys. Res.Space, 113, A05202, https://doi.org/10.1029/2007JA012718, 2008. a
Spence, H. E., Alexandrova, O., Arzamasskiy, L., Argall, M. R., Caprioli, D., Case, A. W., Chandran, B. D. G., Chen, L.J., Dors, I., Eastwood, J. P., Forsyth, C., Galvin, A. B., Genot, V. N., Halekas, J. S., Hesse, M., Horbury, T. S., Jian, L., Kasper, J. C., Klein, K. G., Kretzschmar, M., Kunz, M. W., Lavraud, B., Le Contel, O., Mallet, A., Maruca, B., Matthaeus, W. H., Owen, C. J., Retino, A., Reynolds, C., Roberts, O. W., Schekochihin, A. A., Skoug, R. M., Smith, C. W., Smith, S. S., Steinberg, J. T., Stevens, M. L., Szabo, A., TenBarge, J. M., Torbert, R. B., Vasquez, B. J., Verscharen, D., Whittlesey, P. L., Zank, G. P., and Zweibel, E.: An Overview of HelioSwarm: A NASA MIDEX Mission to Reveal the Nature of Turbulence in Space Plasmas, in: AGU Fall Meeting 2022, AGU Fall Meeting 2022, Chicago, IL, United States, 12–16 December 2022, SH12E1485, American Geophysical Union, https://agu.confex.com/agu/fm22/meetingapp.cgi/Paper/1177250 (last access: 1 February 2024), 2022. a
Torbert, R. B., Dors, I., Argall, M. R., Genestreti, K. J., Burch, J. L., Farrugia, C. J., Forbes, T. G., Giles, B. L., and Strangeway, R. J.: A New Method of 3D Magnetic Field Reconstruction, Geophys. Res. Lett., 47, e2019GL085542, https://doi.org/10.1029/2019GL085542, 2020. a, b, c
Vogt, J., Paschmann, G., and Chanteur, G.: Reciprocal Vectors, in: MultiSpacecraft Analysis Methods Revisit, edited by: Paschmann, G. and Daly, P. W., Chap. 4, p. 33, ESA Communications, ISBN: 9879292219376, 2008. a
Vogt, J., Albert, A., and Marghitu, O.: Analysis of threespacecraft data using planar reciprocal vectors: methodological framework and spatial gradient estimation, Ann. Geophys., 27, 3249–3273, https://doi.org/10.5194/angeo2732492009, 2009. a, b, c, d
Yang, J., Toffoletto, F. R., Wolf, R. A., Sazykin, S., Ontiveros, P. A., and Weygand, J. M.: Largescale current systems and ground magnetic disturbance during deep substorm injections, J. Geophys. Res.Space, 117, A0422, https://doi.org/10.1029/2011ja017415, 2012. a
Zhou, Y., Shen, C., and Ji, Y.: Undulated Shock Surface Formed After a Shock–Discontinuity Interaction, Geophys. Res. Lett., 50, e2023GL103848, https://doi.org/10.1029/2023GL103848, 2023. a
 Abstract
 Introduction
 The problem
 The solution
 Existence and uniqueness of the solution and implications for multispacecraft mission design
 Analytical error analysis
 Summary and discussion
 Appendix A: Multiindex notation
 Appendix B: Scaling coordinates
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 The problem
 The solution
 Existence and uniqueness of the solution and implications for multispacecraft mission design
 Analytical error analysis
 Summary and discussion
 Appendix A: Multiindex notation
 Appendix B: Scaling coordinates
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References