Estimating gradients of physical fields in space

This study focuses on the development of a multi-point 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 non-simultaneous 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.

1 Introduction

Multi-point 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 four-point 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 multi-spacecraft missions (Ogilvie et al., 1977; Escoubet et al., 2001; Liu et al., 2005; Friis-Christensen 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 non-zero, 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 four-point 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.

2 The problem

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 higher-order 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):

$$\begin{array}{}\text{(1)}& j\left(\mathit{x}\right)=f\left(\mathit{x}\right)+w\left(\mathit{x}\right).\end{array}$$

Here, x is an r-component 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 multi-index notation (see Appendix A) as a sum of multi-variate polynomials:

$$\begin{array}{}\text{(2)}& f\left(\mathit{x}\right)=\sum _{\left|\mathit{\alpha }\right|=\mathrm{0}}^{\mathrm{\infty }}{g}_{\mathit{\alpha }}{\mathit{x}}^{\mathit{\alpha }}.\end{array}$$

In this equation, g_α is the coefficient of the polynomial x^α, and we employ the properties of multi-index notation (Properties 3 and 4 in Appendix A). By comparing Eq. (2) with the Taylor expansion of f(x) around $\stackrel{\mathrm{‾}}{\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{‾}}{\mathrm{0}}\right)$ as follows:

$$\begin{array}{}\text{(3)}& g_{\mathit{\alpha }}={\displaystyle \frac{f_{,\mathit{\alpha }}\left(\stackrel{\mathrm{‾}}{\mathrm{0}}\right)}{\mathit{\alpha }\mathrm{!}}},\end{array}$$

where we employ Property 8 of the multi-index notation. Suppose we aim to approximate f(x) using polynomials up to degree d. We define

$$\begin{array}{}\text{(4)}& {\displaystyle }{p}_d\left(\mathit{x}\right)\equiv \sum _{\left|\mathit{\alpha }\right|=\mathrm{0}}^d{g}_{\mathit{\alpha }}{\mathit{x}}^{\mathit{\alpha }},\text{(5)}& {\displaystyle }{p}_d^{+}\left(\mathit{x}\right)\equiv \sum _{\left|\mathit{\alpha }\right|=d+\mathrm{1}}^{\mathrm{\infty }}{g}_{\mathit{\alpha }}{\mathit{x}}^{\mathit{\alpha }}.\end{array}$$

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(\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

$$\begin{array}{}\text{(6)}& j\left(\mathit{x}\right)=p_d\left(\mathit{x}\right)+p_d^{+}\left(\mathit{x}\right)+w\left(\mathit{x}\right).\end{array}$$

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:

$$\begin{array}{}\text{(7)}& j_m=p_d\left({\mathit{x}}_m+\mathit{\delta }\mathit{x}\right)+p_d^{+}\left({\mathit{x}}_m+\mathit{\delta }\mathit{x}\right)+w\left({\mathit{x}}_m+\mathit{\delta }\mathit{x}\right)+\mathit{\delta }j.\end{array}$$

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.

Figure 1 Multi-spacecraft measurement in time and space. (a) The trajectories of multiple probes or spacecraft that fly in formation in time and space are represented by arrows. (b) Measurements of a scalar field f are made along the trajectories.

3 The solution

We define the error between the measured field and the approximating polynomial as s_m, given by

$$\begin{array}{}\text{(8)}& s_m=j_m-p_d\left({\mathit{x}}_m\right).\end{array}$$

To quantify the total error, we employ the weighted least-squares method, which constructs the total error as a weighted sum of all individual errors:

$$\begin{array}{}\text{(9)}& S=\sum _{mn}^M{s}_m{W}_{mn}{s}_n.\end{array}$$

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

$$\begin{array}{}\text{(10)}& W_{mn}={\mathit{\delta }}_{mn}/M.\end{array}$$

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{̃}}{g}}_{\mathit{\beta }}$ results in a set of $\left(\begin{array}{c}d+r\\ r\end{array}\right)$ equations for ${\stackrel{\mathrm{̃}}{g}}_{\mathit{\beta }}$:

$$\begin{array}{}\text{(11)}& {\left.{\displaystyle \frac{\partial S}{\partial g_{\mathit{\beta }}}}\right|}_{{\stackrel{\mathrm{̃}}{g}}_{\mathit{\beta }}}=\mathrm{0}.\end{array}$$

We define the matrix R with elements:

$$\begin{array}{}\text{(12)}& R_{\mathit{\beta }\mathit{\alpha }}\equiv \sum _{mn}^M{\mathit{x}}_m^{\mathit{\beta }}{W}_{mn}{\mathit{x}}_n^{\mathit{\alpha }}.\end{array}$$

Additionally, we define

$$\begin{array}{}\text{(13)}& J_{\mathit{\beta }}\equiv \sum _{mn}^M{\mathit{x}}_m^{\mathit{\beta }}{W}_{mn}{j}_n.\end{array}$$

With these notations, taking into account Eqs. (9), (8), (7), and (4), Eq. (11) can be explicitly expressed as a system of equations:

$$\begin{array}{}\text{(14)}& J_{\mathit{\beta }}=\sum _{\left|\mathit{\alpha }\right|\le d}{R}_{\mathit{\beta }\mathit{\alpha }}{\stackrel{\mathrm{̃}}{g}}_{\mathit{\alpha }}.\end{array}$$

This linear system of equations consists of $\left(\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{̃}}{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{‾}}{\mathrm{0}}$ can be determined. The approximation of the field f(x) is then given by

$$\begin{array}{}\text{(15)}& {\stackrel{\mathrm{̃}}{p}}_d\left(\mathit{x}\right)=\sum _{\left|\mathit{\alpha }\right|\le d}{\stackrel{\mathrm{̃}}{g}}_{\mathit{\alpha }}{\mathit{x}}^{\mathit{\alpha }}.\end{array}$$

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 Existence and uniqueness of the solution and implications for multi-spacecraft mission design

4.1 The requirement of a unique solution

From Eq. (14) there exists a unique set of solutions for ${\stackrel{\mathrm{̃}}{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

$$\begin{array}{}\text{(16)}& R_{\mathit{\beta }\mathit{\alpha }}=\sum _{mnlsk}^M{\mathit{x}}_m^{\mathit{\beta }}(P^T{)}_{ml}{O}_{ls}{O}_{sk}{P}_{kn}{\mathit{x}}_n^{\mathit{\alpha }}.\end{array}$$

Considering the relation rank(A^TA)=rank(A) and the invertibility of OP, we have

$$\begin{array}{}\text{(17)}& \text{rank}\left(\mathbf{R}\right)=\text{rank}\left(\mathbf{OPX}\right)=\text{rank}\left(\mathbf{X}\right),\end{array}$$

where the matrix X is defined by

$$\begin{array}{}\text{(18)}& X_{n\mathit{\alpha }}\equiv {\mathit{x}}_n^{\mathit{\alpha }}.\end{array}$$

Therefore, the uniqueness of the solution in Eq. (14) is equivalent to the rank of X being $\left(\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(\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(\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

$$\begin{array}{}\text{(19)}& \sum _{\left|\mathit{\alpha }\right|\le d}{a}_{\mathit{\alpha }}{\mathit{x}}_m^{\mathit{\alpha }}=\mathrm{0}.\end{array}$$

A similar result has also been obtained for multi-variate 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 second-order surface.

To illustrate the second condition, we take d=2, r=3 as an example. The matrix X in this case is given by

$$\begin{array}{}\text{(20)}& \mathbf{X}=\left[\begin{array}{cccccccccc}\mathrm{1}& x_{\mathrm{1}}& y_{\mathrm{1}}& z_{\mathrm{1}}& x_{\mathrm{1}}^{\mathrm{2}}& x_{\mathrm{1}}{y}_{\mathrm{1}}& x_{\mathrm{1}}{z}_{\mathrm{1}}& y_{\mathrm{1}}^{\mathrm{2}}& y_{\mathrm{1}}{z}_{\mathrm{1}}& z_{\mathrm{1}}^{\mathrm{2}}\\ \mathrm{1}& x_{\mathrm{2}}& y_{\mathrm{2}}& z_{\mathrm{2}}& x_{\mathrm{2}}^{\mathrm{2}}& x_{\mathrm{2}}{y}_{\mathrm{2}}& x_{\mathrm{2}}{z}_{\mathrm{2}}& y_{\mathrm{2}}^{\mathrm{2}}& y_{\mathrm{2}}{z}_{\mathrm{2}}& z_{\mathrm{2}}^{\mathrm{2}}\\ \mathrm{1}& x_{\mathrm{3}}& y_{\mathrm{3}}& z_{\mathrm{3}}& x_{\mathrm{3}}^{\mathrm{2}}& x_{\mathrm{3}}{y}_{\mathrm{3}}& x_{\mathrm{3}}{z}_{\mathrm{3}}& y_{\mathrm{3}}^{\mathrm{2}}& y_{\mathrm{3}}{z}_{\mathrm{3}}& z_{\mathrm{3}}^{\mathrm{2}}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}\\ \mathrm{1}& x_M& y_M& z_M& x_M^{\mathrm{2}}& x_M{y}_M& x_M{z}_M& y_M^{\mathrm{2}}& y_M{z}_M& z_M^{\mathrm{2}}\end{array}\right].\end{array}$$

If all the points lie on a second-order algebraic surface, we can express the surface formally with appropriately chosen coefficients a_α as

$$\begin{array}{}\text{(21)}& \begin{array}{rl}{a}_{\left(\mathrm{0},\mathrm{0},\mathrm{0}\right)}& +a_{\left(\mathrm{1},\mathrm{0},\mathrm{0}\right)}x+a_{\left(\mathrm{0},\mathrm{1},\mathrm{0}\right)}y+a_{\left(\mathrm{0},\mathrm{0},\mathrm{1}\right)}z+a_{\left(\mathrm{2},\mathrm{0},\mathrm{0}\right)}{x}^{\mathrm{2}}\\ & +a_{\left(\mathrm{1},\mathrm{1},\mathrm{0}\right)}xy+a_{\left(\mathrm{1},\mathrm{0},\mathrm{1}\right)}xz+a_{\left(\mathrm{0},\mathrm{2},\mathrm{0}\right)}{y}^{\mathrm{2}}\\ & +a_{\left(\mathrm{0},\mathrm{1},\mathrm{1}\right)}yz+a_{\left(\mathrm{0},\mathrm{0},\mathrm{2}\right)}{z}^{\mathrm{2}}=\mathrm{0},\end{array}\end{array}$$

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(\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 second-order 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(\begin{array}{c}\mathrm{2}+\mathrm{3}\\ \mathrm{3}\end{array}\right)$.

These two conditions have important implications for the orbit design of future multi-spacecraft 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 well-known restriction that at least four measurement points are needed, and these points should not lie on a first-order algebraic surface, or in other words a plane (Fig. 2c), such that with appropriately chosen coefficients a_α they satisfy

$$\begin{array}{}\text{(22)}& a_{\left(\mathrm{0},\mathrm{0},\mathrm{0}\right)}+a_{\left(\mathrm{1},\mathrm{0},\mathrm{0}\right)}x+a_{\left(\mathrm{0},\mathrm{1},\mathrm{0}\right)}y+a_{\left(\mathrm{0},\mathrm{0},\mathrm{1}\right)}z=\mathrm{0}.\end{array}$$

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 second-order algebraic surface, which can be defined by Eq. (21) (Fig. 2f). Typical examples of a second-order 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.

Figure 2 Distributional requirement of measurement points in estimating gradients up to first order (d=1) in (a) one-dimensional (r=1), (b) two-dimensional (r=2), and (c) three-dimensional (r=3) space and gradients up to second order in (d) one-dimensional, (e) two-dimensional, and (f) three-dimensional space. In each panel, the equation defines the hypersurface from which the distribution of measurement points should deviate. The black dots denote the measurement points that lie on some hypersurface defined by the equation with a particular set of coefficients. Such hypersurfaces are represented by the lines in panels (b) and (e), the plane in panel (c), the sphere in panel (f), and the dots in panels (a) and (d). The blue dots represent the additional measurement points not lying on the hypersurfaces.

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 non-linearly 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 non-linear 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 four-spacecraft mission or a lack of spacecraft to form a tetrahedron, resulting in only three spacecraft providing data that lie on a plane. Even in well-functioning four-spacecraft 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 four-spacecraft 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(\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 right-multiplying X by a full-rank 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

$$\begin{array}{}\text{(23)}& {\mathbf{X}}_{nh}^{\prime }=\sum _l^{\left(\begin{array}{c}d+r\\ r\end{array}\right)}{X}_{nl}{G}_{lh}\end{array}$$

such that

$$\begin{array}{}\text{(24)}& X^{\prime }=\left[\begin{array}{cc}{X}_{\mathrm{11}}^{\prime }& \mathrm{\cdots }\\ X_{\mathrm{21}}^{\prime }& \mathrm{\cdots }\\ \mathrm{⋮}& \mathrm{\ddots }\\ X_{M\mathrm{1}}^{\prime }& \mathrm{\cdots }\end{array}\left|\begin{array}{c}\stackrel{N}{\overbrace{\begin{array}{ccc}\mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\end{array}}}\end{array}\right.\right]\end{array}$$

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(\begin{array}{c}d+r\\ r\end{array}\right)$ diagonal matrix and A and B are unitary matrices of order M and $\left(\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 left-multiplying 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

$$\begin{array}{}\text{(25)}& ({\mathbf{X}}^{\prime }{)}^{\mathrm{T}}\mathbf{W}\mathit{j}=({\mathbf{X}}^{\prime }{)}^{\mathrm{T}}{\mathbf{WX}}^{\prime }{\mathbf{G}}^{-\mathrm{1}}\stackrel{\mathrm{̃}}{\mathit{g}},\end{array}$$

where j and $\stackrel{\mathrm{̃}}{\mathit{g}}$ are column vectors. ${\mathbf{G}}^{-\mathrm{1}}\stackrel{\mathrm{̃}}{\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

$$\begin{array}{}\text{(26)}& \left[\begin{array}{c}{\mathit{J}}^{\prime }\\ \\ \mathrm{0}\\ \mathrm{⋮}\\ \mathrm{0}\end{array}\right]=\left[\begin{array}{cc}\begin{array}{c}{\mathbf{R}}^{\prime }\\ \\ \mathrm{0}& \mathrm{\cdots }\\ \mathrm{⋮}& \mathrm{\ddots }\\ \mathrm{0}& \mathrm{\cdots }\end{array}& \underset{N}{\underbrace{\begin{array}{ccc}\mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\end{array}}}\end{array}\right]\left[\begin{array}{c}{\mathit{g}}^{\prime }\\ \\ ({\mathbf{G}}^{-\mathrm{1}}\stackrel{\mathrm{̃}}{\mathit{g}}{)}_{\left(\begin{array}{c}d+r\\ r\end{array}\right)-N+\mathrm{1}}\\ \mathrm{⋮}\\ ({\mathbf{G}}^{-\mathrm{1}}\stackrel{\mathrm{̃}}{\mathit{g}}{)}_{\left(\begin{array}{c}d+r\\ r\end{array}\right)}\end{array}\right],\end{array}$$

where J^′ and R^′ contain non-vanishing components and ${\stackrel{\mathrm{̃}}{\mathit{g}}}^{\prime }$ includes the first $\left(\begin{array}{c}d+r\\ r\end{array}\right)-N$ components of ${\mathbf{G}}^{-\mathrm{1}}\stackrel{\mathrm{̃}}{\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{̃}}{\mathit{g}}$ the soluble ${\stackrel{\mathrm{̃}}{\mathit{g}}}^{\prime }$. By solving for them from

$$\begin{array}{}\text{(27)}& {\mathit{J}}^{\prime }={\mathbf{R}}^{\prime }{\stackrel{\mathrm{̃}}{\mathit{g}}}^{\prime },\end{array}$$

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{̃}}{\mathit{g}}=\left[f\right(\stackrel{\mathrm{‾}}{\mathrm{0}}),{\partial }_{x}f(\stackrel{\mathrm{‾}}{\mathrm{0}}),{\partial }_{y}f(\stackrel{\mathrm{‾}}{\mathrm{0}}),{\partial }_{z}f(\stackrel{\mathrm{‾}}{\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(\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 right-multiplication 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(\begin{array}{c}d+r\\ r\end{array}\right)$, and $\left(\genfrac{}{}{0ex}{}{d+r}{r}\right)-N<h\le \left(\begin{array}{c}d+r\\ r\end{array}\right)$. In the first $\left(\begin{array}{c}d+r\\ r\end{array}\right)-N$ columns, in addition to the rest being zeros, $\left(\begin{array}{c}d+r\\ r\end{array}\right)-N$ unities are so placed that the following two conditions are met.

1. G^′ has full rank.

2. Let the row (column) index of unity be i (j). If $\left(\begin{array}{c}u-\mathrm{1}+r\\ r\end{array}\right)<i\le \left(\begin{array}{c}u+r\\ r\end{array}\right)$ for some u and $\left(\begin{array}{c}v-\mathrm{1}+r\\ r\end{array}\right)<j\le \left(\begin{array}{c}v+r\\ r\end{array}\right)$ for some v, then we should have u=v.

5 Analytical error analysis

While a unique solution can be obtained for estimating ${\stackrel{\mathrm{̃}}{g}}_{\mathit{\alpha }}$ and ${\stackrel{\mathrm{̃}}{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⁻¹ be the inverse of R. We multiply Eq. (14) by ${\left(R^{-\mathrm{1}}\right)}_{\mathit{\gamma }\mathit{\beta }}$, sum over β, and obtain

$$\begin{array}{}\text{(28)}& \begin{array}{rl}& \sum _{\left|\mathit{\beta }\right|\le d}{\left(R^{-\mathrm{1}}\right)}_{\mathit{\gamma }\mathit{\beta }}\sum _{mn}{\mathit{x}}_m^{\mathit{\beta }}{W}_{mn}[p_d\left({\mathit{x}}_n+\mathit{\delta }\mathit{x}\right)\\ & +p_d^{+}\left({\mathit{x}}_n+\mathit{\delta }\mathit{x}\right)+w\left({\mathit{x}}_n+\mathit{\delta }\mathit{x}\right)+\mathit{\delta }j]={\stackrel{\mathrm{̃}}{g}}_{\mathit{\gamma }},\end{array}\end{array}$$

where use was made of Eqs. (13) and (7). According to the binomial theorem for multi-variate polynomials (see Eq. A1), we have the decomposition

$$\begin{array}{}\text{(29)}& p_d\left(\mathit{x}+\mathit{\delta }\mathit{x}\right)=\sum _{\left|\mathit{\alpha }\right|\le d}{g}_{\mathit{\alpha }}\sum _{\stackrel{\mathrm{‾}}{\mathrm{0}}\le \mathit{\lambda }<\mathit{\alpha }}\left(\begin{array}{c}\mathit{\alpha }\\ \mathit{\lambda }\end{array}\right){\mathit{x}}^{\mathit{\lambda }}\mathit{\delta }{\mathit{x}}^{\mathit{\alpha }-\mathit{\lambda }}+\sum _{\left|\mathit{\alpha }\right|\le d}{g}_{\mathit{\alpha }}{\mathit{x}}^{\mathit{\alpha }}.\end{array}$$

Substituting this into Eq. (28), subtracting g_γ from both sides, and defining the error in estimating g_γ as

$$\begin{array}{}\text{(30)}& \mathit{\delta }{g}_{\mathit{\gamma }}\equiv {\stackrel{\mathrm{̃}}{g}}_{\mathit{\gamma }}-g_{\mathit{\gamma }},\end{array}$$

we obtain the complete expression for the error:

$$\begin{array}{}\text{(31)}& \begin{array}{rl}& \sum _{\left|\mathit{\beta }\right|\le d}{\left(R^{-\mathrm{1}}\right)}_{\mathit{\gamma }\mathit{\beta }}\sum _{mn}^M{\mathit{x}}_m^{\mathit{\beta }}{W}_{mn}\left[\sum _{\left|\mathit{\alpha }\right|\le d}{g}_{\mathit{\alpha }}\sum _{\stackrel{\mathrm{‾}}{\mathrm{0}}\le \mathit{\lambda }<\mathit{\alpha }}\left(\begin{array}{c}\mathit{\alpha }\\ \mathit{\lambda }\end{array}\right){\mathit{x}}_n^{\mathit{\lambda }}\mathit{\delta }{\mathit{x}}^{\mathit{\alpha }-\mathit{\lambda }}\right.\\ & \left.+p_d^{+}\left({\mathit{x}}_n+\mathit{\delta }\mathit{x}\right)+w\left({\mathit{x}}_n+\mathit{\delta }\mathit{x}\right)+\mathit{\delta }j\phantom{\sum _{\stackrel{\mathrm{‾}}{\mathrm{0}}\le \mathit{\lambda }<\mathit{\alpha }}\left(\begin{array}{c}\mathit{\alpha }\\ \mathit{\lambda }\end{array}\right)}\right]=\mathit{\delta }{g}_{\mathit{\gamma }}.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(32)}& {\displaystyle \frac{{\left(\mathit{\delta }{g}_{\mathit{\gamma }}\right)}_{\mathrm{t}}}{g_{\mathit{\gamma }}}}=\sum _{\left|\mathit{\alpha }\right|>d}{\displaystyle \frac{g_{\mathit{\alpha }}}{g_{\mathit{\gamma }}}}\sum _{\left|\mathit{\beta }\right|\le d}{\left(R^{-\mathrm{1}}\right)}_{\mathit{\gamma }\mathit{\beta }}\sum _{mn}{\mathit{x}}_m^{\mathit{\beta }}{W}_{mn}\left({\mathit{x}}_n^{\mathit{\alpha }}+\mathit{\delta }\mathit{x}\right).\end{array}$$

It is obvious that three factors combine to create this error. The first is the ratio of higher-order 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

$$\begin{array}{}\text{(33)}& {\displaystyle \frac{{\left(\mathit{\delta }{g}_{\mathit{\gamma }}\right)}_{\mathrm{t}}}{g_{\mathit{\gamma }}}}=\sum _{\left|\mathit{\alpha }\right|>d}{\displaystyle \frac{g_{\mathit{\alpha }}}{g_{\mathit{\gamma }}}}{q}_{\mathit{\alpha }\mathit{\gamma }}^{\mathit{\#}}\underset{m}{max}|{\mathit{x}}_m{|}^{\left|\mathit{\alpha }\right|-\left|\mathit{\gamma }\right|},\end{array}$$

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

$$\begin{array}{}\text{(34)}& \sum _m{\mathit{x}}_m=\stackrel{\mathrm{‾}}{\mathrm{0}}.\end{array}$$

Thus, Eq. (33) can be re-expressed as

$$\begin{array}{}\text{(35)}& {\displaystyle \frac{{\left(\mathit{\delta }{g}_{\mathit{\gamma }}\right)}_{\mathrm{t}}}{g_{\mathit{\gamma }}}}=\sum _{\left|\mathit{\alpha }\right|>d}{\displaystyle \frac{g_{\mathit{\alpha }}}{g_{\mathit{\gamma }}}}{\displaystyle \frac{\mathrm{1}}{q_{\mathit{\alpha }\mathit{\gamma }}}}{L}^{\left|\mathit{\alpha }\right|-\left|\mathit{\gamma }\right|},\end{array}$$

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):

$$\begin{array}{}\text{(36)}& q_{\mathit{\alpha }\mathit{\gamma }}={\displaystyle \frac{L^{\left|\mathit{\alpha }\right|-\left|\mathit{\gamma }\right|}}{{\sum }_{\left|\mathit{\beta }\right|\le d}{\left(R^{-\mathrm{1}}\right)}_{\mathit{\gamma }\mathit{\beta }}{\sum }_{mn}{\mathit{x}}_m^{\mathit{\beta }}{W}_{mn}{\mathit{x}}_n^{\mathit{\alpha }}}}.\end{array}$$

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{̃}}{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 low-pass 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{̃}}{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

$$\begin{array}{}\text{(37)}& f\left(\mathit{x}\right)-{\stackrel{\mathrm{̃}}{p}}_d\left(\mathit{x}\right)=p_d^{+}\left(\mathit{x}\right)-\sum _{\left|\mathit{\gamma }\right|\le d}\mathit{\delta }{g}_{\mathit{\gamma }}{\mathit{x}}^{\mathit{\gamma }}.\end{array}$$

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(\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 left-multiply Eq. (27) by ${\mathbf{R}}^{\prime -\mathrm{1}}$ to obtain

$$\begin{array}{}\text{(38)}& \begin{array}{rl}& \sum _h^{\left(\begin{array}{c}d+r\\ r\end{array}\right)-N}{\left(R^{\prime -\mathrm{1}}\right)}_{lh}\sum _{mn}^M{\left(X^{\prime \mathrm{T}}\right)}_{hm}{W}_{mn}[p_d\left({\mathit{x}}_n+\mathit{\delta }\mathit{x}\right)\\ & +p_d^{+}\left({\mathit{x}}_n+\mathit{\delta }\mathit{x}\right)+w\left({\mathit{x}}_n+\mathit{\delta }\mathit{x}\right)+\mathit{\delta }j]={\stackrel{\mathrm{̃}}{g}}_l^{\prime }.\end{array}\end{array}$$

By defining

$$\begin{array}{}\text{(39)}& \mathit{\delta }{g}_l^{\prime }\equiv {\stackrel{\mathrm{̃}}{g}}_l^{\prime }-g_l^{\prime },\end{array}$$

we have the following expression for it:

$$\begin{array}{}\text{(40)}& \begin{array}{rl}& \sum _h^{\left(\begin{array}{c}d+r\\ r\end{array}\right)-N}{\left(R^{\prime -\mathrm{1}}\right)}_{lh}\sum _{mn}^M{\left(X^{\prime \mathrm{T}}\right)}_{hm}{W}_{mn}\\ & \left[\sum _{\left|\mathit{\alpha }\right|\le d}{g}_{\mathit{\alpha }}\sum _{\stackrel{\mathrm{‾}}{\mathrm{0}}\le \mathit{\lambda }<\mathit{\alpha }}\right(\begin{array}{c}\mathit{\alpha }\\ \mathit{\lambda }\end{array}){\mathit{x}}_n^{\mathit{\lambda }}\mathit{\delta }{\mathit{x}}^{\mathit{\alpha }-\mathit{\lambda }}+p_d^{+}\left({\mathit{x}}_n+\mathit{\delta }\mathit{x}\right)\\ & +w\left({\mathit{x}}_n+\mathit{\delta }\mathit{x}\right)+\mathit{\delta }j\phantom{\sum _{\stackrel{\mathrm{‾}}{\mathrm{0}}\le \mathit{\lambda }<\mathit{\alpha }}\left(\begin{array}{c}\mathit{\alpha }\\ \mathit{\lambda }\end{array}\right)}]=\mathit{\delta }{g}_l^{\prime }.\end{array}\end{array}$$

The relative error caused by truncation is given by

$$\begin{array}{}\text{(41)}& {\displaystyle \frac{{\left(\mathit{\delta }{g}_l^{\prime }\right)}_{\mathrm{t}}}{g_l^{\prime }}}=\sum _{\left|\mathit{\alpha }\right|>d}{\displaystyle \frac{g_{\mathit{\alpha }}}{g_l}}\sum _h^{\left(\begin{array}{c}d+r\\ r\end{array}\right)-N}{\left(R^{\prime -\mathrm{1}}\right)}_{lh}\sum _{mn}^M{\left(X^{\prime \mathrm{T}}\right)}_{hm}{W}_{mn}\left({\mathit{x}}_n^{\mathit{\alpha }}+\mathit{\delta }\mathit{x}\right).\end{array}$$

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

$$\begin{array}{}\text{(42)}& {\displaystyle \frac{{\left(\mathit{\delta }{g}_l^{\prime }\right)}_{\mathrm{t}}}{g_l^{\prime }}}=\sum _{\left|\mathit{\alpha }\right|>d}{\displaystyle \frac{g_{\mathit{\alpha }}}{g_l}}{\displaystyle \frac{\mathrm{1}}{q_{\mathit{\alpha }l}^{\prime }}}{L}^{\left|\mathit{\alpha }\right|-u},\text{if }\left(\begin{array}{c}u-\mathrm{1}+r\\ r\end{array}\right)<l\le \left(\begin{array}{c}u+r\\ r\end{array}\right),\end{array}$$

where the quality factor is given by

$$\begin{array}{}\text{(43)}& \begin{array}{rl}{q}_{\mathit{\alpha }l}^{\prime }=& {\displaystyle \frac{L^{\left|\mathit{\alpha }\right|-u}}{{\sum }_h^{\left(\begin{array}{c}d+r\\ r\end{array}\right)-N}{\left(R^{-\mathrm{1}}\right)}_{lh}{\sum }_{mn}^M{\left(X^{\prime \mathrm{T}}\right)}_{hm}{W}_{mn}{\mathit{x}}_n^{\mathit{\alpha }}}}.\\ & \text{if }\left(\begin{array}{c}u-\mathrm{1}+r\\ r\end{array}\right)<l\le \left(\begin{array}{c}u+r\\ r\end{array}\right)\end{array}\end{array}$$

Therefore, this method designed for cases when measurement points are not well distributed has good accuracy.

6 Summary and discussion

The techniques for calculating linear gradients of general physical fields and quadratic gradients of magnetic fields using four-point measurements have been widely applied in the context of multi-spacecraft 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 multi-point techniques that rely on more points to estimate quadratic and higher-order gradients.

In summary, we have analytically explored the general method to estimate gradients of fields in space based on multi-point 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(\begin{array}{c}d+\mathrm{3}\\ \mathrm{3}\end{array}\right)$ spacecraft are needed, and these spacecraft should not lie on a dth-order surface in space. In particular, at least 10 points that are not on a second-order 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(\begin{array}{c}d+\mathrm{4}\\ \mathrm{4}\end{array}\right)$ measurement points that do not lie on a dth-order 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 co-moving 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 ath-order 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 multi-point techniques that rely on gradients of physical fields. Additionally, they hold significance for the future design of multi-spacecraft missions aimed at studying physics associated with quadratic or higher-order 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 divergence-free 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 zeroth-order gradients of electric currents and first-order gradients of magnetic fields are correlated by $\mathrm{\nabla }\times \mathit{B}={\mathit{\mu }}_{\mathrm{0}}\mathit{j}$. Beyond these linear constraints, non-linear 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 divergence-free 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 higher-order 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).

Appendix A: Multi-index notation

Here we list the properties of a multi-index notation tailored for multi-variate functions (Riachy et al., 2011).

Let $\mathit{\alpha }=\left({\mathit{\alpha }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\alpha }}_r\right)$ be an r-tuple of non-negative integers ${\mathit{\alpha }}_i,i=\mathrm{1},\mathrm{\dots },r;i,r\in \mathbb{N}$. α is called a multi-index. The symbol in bold x denotes a vector in ℝ^r. As for time and space, r=4.

For multi-indices $\mathit{\alpha },\mathit{\beta }\in {\mathbb{N}}^r$, the following properties are either defined or deduced.

1. 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)$.

2. 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{\}}$.

3. 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}$.

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

5. Factorial: α!=α₁!⋯α_r!.

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

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

8. Higher-order 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 }}$.

9. Denote by 1_i∈ℕ^r the multi-index 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)$.

10. 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}}}$.

11. Binomial theorem:

    $$\begin{array}{}\text{(A1)}& {\left(\mathit{x}+\mathit{y}\right)}^{\mathit{\alpha }}=\sum _{\stackrel{\mathrm{‾}}{\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}$$

Appendix B: Scaling coordinates

In addition to the scaling of a temporal coordinate by a characteristic speed to obtain x₀=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).