Gradient estimation using configurations of two or three spacecraft

The forthcoming three-satellite mission Swarm will allow us to investigate plasma processes and phenomena in the upper ionosphere from an in-situ multi-spacecraft perspective. Since with less than four points in space the spatiotemporal ambiguity cannot be resolved fully, analysis tools for estimating spatial gradients, wave vectors, or boundary parameters need to utilise additional information such as geometrical or dynamical constraints. This report deals with gradient estimation where the planar component is constructed using instantaneous three-point observations or, for quasi-static structures, by means of measurements along the orbits of two close spacecraft. A new least squares (LS) gradient estimator for the latter case is compared with existing finite difference (FD) schemes and also with a three-point LS technique. All available techniques are presented in a common framework to facilitate error analyses and consistency checks, and to show how arbitrary combinations of planar gradient estimators and constraints can be formed. The accuracy of LS and FD planar gradient estimators is assessed in terms of prescribed and adjustable discretization parameters to optimise their performance along the satellite orbits. Furthermore, we discuss the implications of imperfect constraint equations for error propagation, and address the effects of sub-scale structures. The two-spacecraft LS scheme is demonstrated using Cluster FGM measurements at a planar and essentially force-free plasma boundary in the solar wind where all three different types of constraints to construct out-of-plane derivatives can be applied.


Introduction
Over the past decade, magnetic measurements made by Low Earth Orbiting (LEO) satellites such as Ørsted (Neubert et al., 2001), CHAMP (Reigber et al., 2002), and also SAC-C have significantly advanced our understanding of the geomagnetic field and its dynamics. Models like POMME (e.g., Maus et al., 2006) and CHAOS (Olsen et al., 2006) resolve variability on small spatial and temporal scales and allow to study, e.g., crustal magnetization and secular variations in unprecedented detail. Further improvements in geomagnetic main (internal) field modelling will require identifying the internal field signal at LEO altitudes with even higher accuracy, and here the contributions from electrical currents in space impose severe limitations. Separating these external contributions from the internal field in magnetic measurements is thus a key challenge in geomagnetism (Lühr et al., 2009).
A current system that arguably causes the most severe problems in geomagnetic field modelling from LEO satellite data is associated with magnetic field-aligned currents (FAC) that flow along magnetic flux tubes in the auroral zone. These FAC connect the ionospheric end of auroral flux tubes where electromagnetic energy is dissipated through conduction (Pedersen) currents with the equatorial magnetosphere. Auroral FAC are highly variable in both time and space, in particular during magnetospheric substorms. At LEO altitudes they affect the measurements of geomagnetic field missions directly. In the transient phase when the auroral current circuit is established, FAC are communicated by shear Alfvén waves that travel along geomagnetic field lines and do not suffer from geometrical attenuation. Auroral phenomena are closely linked with FAC. The so-called auroral electrojet is an ionospheric current system fed by FAC that gives rise to magnetic disturbances of up to several 1000 nT on the ground. For a discussion of FAC and Alfvén waves in the auroral context, see Paschmann et al. (2002), Vogt (2002), and Keiling (2009). A FAC model based on ten years of CHAMP measurements was recently presented by He et al. (2012).
The upcoming Swarm mission consists of three LEO spacecraft on polar orbits at altitudes around 400 km and 530 km. The lower two will fly side-by-side with a longitudinal separation of about 1.5 • . In addition to scalar and vector magnetometers, the Swarm satellites will also be equipped with an electric field instrument (EFI) based on ion drift measurements. Since processes that vary in time and also in several independent spatial dimensions cannot be resolved by single-spacecraft measurements, the multi-spacecraft mission Swarm will introduce a new quality in geomagnetic observations from LEO satellites.
To tap the full potential of Swarm multi-point observations, analysis techniques need to be tailored to the specifics of the mission. The problem of FAC estimation based on insitu measurements from two Swarm satellites was addressed by Ritter and Lühr (2006) who used simulated FAC structures to test their method. Vogt et al. (2009) introduced an analysis framework for three-spacecraft configurations based on the set of planar reciprocal vectors as a generic tool to estimate the components of parameter vectors (spatial gradients, wave vectors, boundary normals) in the plane spanned by the three position vectors. To estimate the out-of-plane component of parameters vectors, additional information has to be considered, e.g., in the form of geometrical constraints, stationarity assumptions, or other conditions. Shen et al. (2012) carried this approach further and constructed a scheme to estimate the full magnetic gradient matrix for stationary fieldaligned current structures from measurements along the orbits of two spacecraft. This paper brings together the various approaches to gradient and FAC estimation for two-spacecraft and threespacecraft configurations such that measurements and constraints can be freely combined. The resulting data analysis framework is supposed to facilitate consistency checks and error analyses, thus enhancing both the significance and the accuracy of parameter estimates based on Swarm measurements. Using the least squares (LS) principle, we derive a planar gradient estimator for quasi-static structures observed along the orbits of two spacecraft. Existing finite difference (FD) schemes are presented in a canonical form and then compared with the new LS technique. The force-free constraint formulated by Shen et al. (2012) is converted into the same algebraic form as the geometrical constraints of Vogt et al. (2009). We study statistical errors and truncations errors of the discretization schemes, and discuss possibilities to improve their performance. The two-spacecraft LS gradient estimator is applied to magnetic field measurements of two Cluster satellites in the solar wind during their transit through a directional discontinuity. Geometrical constraints and the force-free condition are used to produce curl estimates that can be compared with the boundary curl profile obtained from four-point crossing time analysis.
The case of S = 4 spacecraft received special attention because of the Cluster mission. Here all instantaneous parameters of a general linear field in space (offset and all components of the gradient) can be uniquely determined from the data if no further constraints are employed (Vogt et al., 2008). Hence, in this case all linear gradient estimators give identical results. Explicitly, the identity of the homogeneous least squares and the linear spatial interpolation methods was shown by Chanteur and Harvey (1998). The equivalence of the latter with the finite difference (FD) approach follows directly from the algebra of reciprocal vectors (see Eq. 14.6 in ) that ensures the general FD condition where f σ and f τ denotes the data at positions r σ and r τ , respectively, and is the gradient estimator based on reciprocal vectors k σ , see . The multi-spacecraft gradient estimation problem is overdetermined when the number of spacecraft S > 4. A LS solution still exists (Harvey, 1998;De Keyser et al., 2007;De Keyser, 2008) and can be written in a form similar to Eq. (2), namely, see Vogt et al. (2008). The generalised reciprocal vectors are defined through the position tensor in a mesocentric frame of reference characterised by σ r σ = 0.
Ann. Geophys., 31, 1913-1927, 2013 www.ann-geophys.net/31/1913/2013/ With less than four spacecraft, the gradient estimation problem is underdetermined and its solution no longer unique. Non-collinear configurations of S = 3 spacecraft define a plane, and measurements give information only on the planar component ∇ p f of the gradient. An estimator for ∇ p f is the minimum norm solution of the LS problem (Vogt et al., 2009): Here the q σ 's denote planar reciprocal vectors given by where r τ ν = r ν − r τ are relative position vectors, (σ, τ, ν) is the cyclic permutation of (1, 2, 3) with σ in the first position, and n is the normal vector defined through n = r 12 × r 13 .
The corresponding normal unit vector isn = n/|n|. The planar component ∇ p and the normal component ∇ n of the del (nabla) operator ∇ are formally related through wheren·∇ ≡ ∂/∂n is the directional derivative along the unit vectorn, i.e., the normal derivative. The full gradient can be reconstructed as follows: To obtain an estimate for the normal derivative, Vogt et al. (2009) considered two types of geometrical constraints. When the spatial gradient ∇f is parallel to a given unit vectorê, then When the spatial gradient ∇f is perpendicular to a given unit vectorê, then These formulas for the normal derivatives were derived in the three-spacecraft context, but they are more general as they do not depend on the specifics of the planar gradient estimator, and thus can also be used for two-spacecraft estimation of quasi-static spatial gradients as described in Sect. 3. Any linear estimator for both the full gradient or its planar component can be written in a form analogous to Eqs. (2), (3), and (6) with the reciprocal vectors possibly replaced by others. We are going to refer to this representation as the canonical form of a linear gradient estimator, and call the associated set of vectors like k σ or q σ its canonical base vectors. The canonical form facilitates comparison with other estimators, and allows to use elements of the error analysis framework developed by  for Eq. (2), see also Vogt et al. (2008Vogt et al. ( , 2009).

Gradient estimation schemes for Swarm
This section deals with gradient estimation schemes for the Swarm mission like those proposed by Ritter and Lühr (2006) and Shen et al. (2012). The two Swarm satellites a and b orbiting at the same altitude collect data in the plane spanned by the separation vector and the average velocity vector. Spatial structures that vary only weakly during the transit of the satellites can be considered quasi-static in this context. For such structures an estimate of the planar gradient can be obtained by combining the along-track variations with the measured differences between the two satellites. Since the configuration does not probe the normal gradient, additional conditions have to be taken into account. The general procedure of first estimating the planar gradient and then inferring its normal component is thus similar to gradient estimation from three-spacecraft data as described by Vogt et al. (2009), with the limitation that time-varying gradients have to be excluded from the two-spacecraft case.

Geometry of planar four-point configurations
When the local trajectories of two spacecraft a and b are in one plane, and measurements at two (or N) points along each orbit are taken into account, one effectively obtains planar four-point (or planar 2N-point) configurations. Spacecraft velocities are V a and V b . The mesocenter is moving at the velocity The difference vector allows to express the spacecraft positions at a centre time t * as follows: Since spatial gradients are independent of constant offsets, we may assume the frame of reference to be mesocentric: r * = 0. Mesocentric position vectors are r a = −r and r b = r . Throughout this paper, we assume coordinate systems to be mesocentric unless the more general case is explicitly addressed. Figure 1 shows the main parameters a planar four-point configuration for the estimation of stationary spatial gradients. We need to choose a time interval t that gives suitable separations of measurement points along the orbits. The resulting configuration consists of the following four points: For the two Swarm satellites that orbit side-by-side, the directions of their spacecraft velocities V a and V b differ slightly, but the magnitudes can be considered identical, otherwise the orbits would soon be out-of-phase: |V a | = |V b |. This implies that the difference vector must be perpendicular to the mesocenter velocity V * . Following Ritter and Lühr (2006) and Shen et al. (2012), we now establish a local planar coordinate system where one axis (v) is aligned with the mesocenter velocity. The coordinate unit vector is thusv = V * /|V * |. The second axis (u) completes the planar orthogonal frame, see Fig. 1. With the definitions the (u, v) coordinates of the measurement points are The set of discretization parameters (L, M, , m) represents four degrees of freedom of which two are prescribed by the orbit geometry, namely, the (local) values of L and m/M = ε. L is the spatial separation between the satellite orbits, and ε indicates how much they deviate from an ideal parallel case, with ε = 1 corresponding to perpendicular orbits. In the case of Swarm, L ≈ 80 km at the equator (the total separation is 160 km) and about less than half this value in the auroral zone. The parameter m/M = ε is not larger than 0.015 and can thus be treated as a small parameter.
The remaining two degrees of freedom are captured in the parameters M and or, alternatively, the scaled variables M/L = µ and /L = λ. We may call them adjustable parameters because M or µ can be changed through the time interval t, and or λ can be varied by introducing a time shift τ between the reference times of the two spacecraft. Such a time shift was employed by Ritter and Lühr (2006)  to enforce a symmetric configuration with 0. Geometrically, µ = M/L measures the along-track separation relative to the across-track distance, and λ = /L indicates how close the polygon is to an equal-sided trapezoid. To obtain error formulas and quality measures in terms of the adjustable parameters, we are going to write the formulas also in terms of the set (L, m/M, M/L, /L) = (L, ε, µ, λ). For the measurement points r ± a,b we obtain Note that the parameters t and τ should be increased only within certain bounds to ensure that the quasi-static assumption is not violated. Both the time interval and the time shift effectively increase the time scale on which the structures are required to vary only weakly.

LS approach to planar gradient estimation
A least squares (LS) estimator for planar four-point configurations can be constructed by minimising the cost function with respect to the planar gradient estimator∇ p f and the constant value f * . In a non-mesocentric frame, we need to replace r σ by (r σ −r * ). The problem is overdetermined, and the algebra is analogous to the full (three-dimensional) gradient estimation problem when measurements from S > 4 satellites are available (Vogt et al., 2008). The solution can be written in the form Ann. Geophys., 31, 1913-1927, 2013 www.ann-geophys.net/31/1913/2013/ where the canonical base vectors q σ are (minimum-norm) solutions of the equations Rq σ = r σ . Here R = 4 σ =1 r σ r T σ is the position tensor. As long as the four measurement points are located in a plane but not on a line, the position tensor has rank 2, with two positive eigenvalues λ 1 , λ 2 (note that R is positive definite), and corresponding eigenvectorsê 1 ,ê 2 . The third eigenvalue is zero, and the eigenvector is normal to the plane spanned by the four measurement points. The pseudo-inverse of R is then given by and the canonical base vectors are In (u, v) coordinates, the position tensor R is a (2 × 2) matrix with a regular inverse R −1 as long as the four measurement points are not collinear. Representations of the tensor R −1 and the vectors q ± a,b in (u, v) coordinates are given in Appendix A.
The accuracy of the solutions q σ , and thus also of the gradient estimates, is controlled by the condition number c = λ 1 /λ 2 of R. The condition number is conveniently analysed in (u, v)-coordinates where R is a (2 × 2) matrix so that c may be expressed as see Appendix A for a proof. Here T and D are the trace and the determinant of the matrix, respectively. The condition number is a function of T 2 /D. Using the expressions in Appendix A we obtain The logarithm of the condition number c = c(λ, µ) is displayed in Fig. 2 for ε = 0.01. The ideal case c = 1 corresponds to λ = 0 and µ = 1. As long as λ is not much larger and µ not much smaller than unity, the condition number assumes tolerable values. Very large values of λ have a negative effect on the stability of the matrix inversion. In such a case the orbit phases should be adjusted through an appropriate time shift τ as in Ritter and Lühr (2006).

FD approach to planar gradient estimation
With |r | 2 = L 2 + 2 , |r | sin β = L, and L cos β = sin β, Eq. (10) of Shen et al. (2012) can be rearranged as follows: when the values f b and f a are replaced by arithmetic means: The expressions can be combined to write the FD planar gradient estimator of Shen et al. (2012) in (u, v) coordinates in the canonical form where the canonical base vectors are Ritter and Lühr (2006) also employed a FD approach to approximate partial derivatives for estimating the FAC density. They used x and y to denote the along-track and the cross-track coordinate, respectively. Furthermore, they applied a time shift to adjust the orbit phases of satellites a and b such that the separation vector 2 r is perpendicular to the mesocenter velocity, thus effectively enforcing 0. The estimation scheme for partial derivatives is found from Eqs. (3) and (5) in Ritter and Lühr (2006): where the equivalences with the present notation are Using the appropriate finite differences for dB y and dB x , we obtain Comparison with Eqs. (36) and (37) shows that the ∂B y /∂x and ∂B x /∂y equations correspond to the → 0 limit of the FD formulas for ∂B u /∂v and ∂B v /∂u of Shen et al. (2012). Hence there is no need for a separate assessment of the Ritter and Lühr (2006) FD method.
In Appendix A it is found that the FD approach to planar four-point gradient estimation differs from the LS approach only in terms proportional to ε:

Boundary integral approach to curl estimation
A second approach to FAC estimation mentioned by Ritter and Lühr (2006) is the boundary integral (BI) method based on Ampère's integral law where path integration is along the boundary of the area A, see also Dunlop et al. (1988). For the polygon associated with the planar four-point configuration of Fig. 1, the boundary integral is most naturally discretized using the trapezoidal rule on each of the four legs. The algebra is summarised in Appendix B. A salient conclusion from the analysis concerns the equivalence with curl estimation based on the FD schemes discussed in Sect. 3.3, namely, the discrete form of the boundary integral gives exactly the same result as the FD curl estimator formed by combining the appropriate discrete partial derivatives. Another useful result of the algebra in Appendix B is a particularly compact representation of the current estimator: where the area A is the modulus of the oriented area

Force-free constraint for normal derivatives
The estimators discussed in Sects. 3.2, 3.3, and 3.4 utilise all the information from measurements of planar four-point configurations. To construct estimators for derivatives in the direction normal to the spacecraft plane, one needs to employ additional information, e.g., in the form of geometrical or physical constraints. Two types of geometrical constraints that can be used for this purpose are presented in Sect. 2, see also Vogt et al. (2009). In the auroral context, an appropriate condition follows from the observation that electrical currents j are typically aligned with the ambient magnetic field B. Since the associated force term j × B in the magnetohydrodynamic equation of motion vanishes, we may refer to this condition as the force-free constraint. The FAC density j || can then be inferred from the normal current density j n through see also Ritter and Lühr (2006). In the case of circular orbits, n is the radial direction, and thus j || = j n / sin I where I is the inclination of the magnetic field. Shen et al. (2012) pursued this approach further and showed that all normal derivatives can be constructed from the planar gradient if the forcefree constraint and the divergence-free nature of the magnetic field are taken into account. In Appendix C these conditions are rearranged into the form that allows direct comparison with constraint equations (12) and (13). LS estimators for ∇ p ×B and ∇ p ·B are found from the corresponding gradient formulas: FD estimators are obtained by replacing the canonical base vectors: q σ → p σ .

Accuracy of gradient estimation schemes
In assessing the quality of gradient estimators and their robustness against undesirable effects like measurement noise or deviations from underlying assumptions, we consider separately the planar gradient (Sect. 4.1), the normal derivatives (Sect. 4.2), and deviations from linearity (Sect. 4.3).

Accuracy of planar gradient estimation
The main sources of error considered in the analysis of four-point Vogt and Paschmann, 1998) and three-point (Vogt et al., 2009) gradient estimators are instrumental noise (physical errors) and inaccuracies in the spacecraft positions (geometrical errors). To assess their relative importance for planar gradient estimation from Swarm measurements, we follow Vogt and Paschmann (1998) and compare the instrumental error δB ∼ 1 nT with |∇B| δr where δr denotes the inaccuracy in spacecraft position. The magnetic field variation associated with a typical auroral FAC density of 1 µA m −2 is |∇B| ∼ 10 −3 nT m −1 . For LEO satellites, δr is in the range of several ten meters, hence |∇B| δr is expected to be one to even two orders of magnitude smaller than δB which means that errors due to positional inaccuracies should be negligible compared to instrumental errors. This observation simplifies the error analysis considerably. We first consider the three-spacecraft estimator of the instantaneous (not necessarily quasi-static) planar gradient as introduced by Vogt et al. (2009). The accuracy of the method was studied in the same paper (Sect. 5.2 and Appendix B) using an error analysis framework introduced by Chanteur (1998) that essentially rests on the algebraic structure of the gradient estimator. Assuming the instrumental errors to be mutually uncorrelated and isotropic, the covariance of the planar gradient vector of a component B j of the magnetic field is given by where Q = 3 σ =1 q σ q T σ is the reciprocal tensor. In the expression for the square magnitude error instrumental errors are amplified by the factor For details, see Vogt et al. (2009). Now we look at the FD and LS four-point estimators of quasi-static planar gradients presented in Sects. 3.3 and 3.2. Since each estimator is available in canonical form and thus in the same algebraic structure as the three-spacecraft scheme discussed in the previous paragraph, both FD and LS estimators can be studied in the same way as before. One just needs to insert the appropriate canonical base vectors to obtain the error amplification factors for the FD scheme: and for the LS scheme: where D is the determinant of the position tensor, see Eq. (A3). The FD and LS sets of canonical base vectors differ only by terms of order O(ε), hence their amplification factors are identical up to O ε 2 . For the Swarm satellites a and b, ε 10 −2 which implies that the LS and FD estimators should give almost identical results in practice. Figure 3 displays the product i.e., the error amplification factor of the LS estimator normalised by the mean square inter-spacecraft distance (1/4) 4 σ =1 |r σ | 2 . The normalised error amplification factor of the LS estimator turns out to be identical with T 2 /4D where T and D are the trace and the determinant, respectively, of the position tensor R, and thus shows the same behaviour as the condition number c, see Sect. 3.2. Configurational error amplification occurs most prominently for large values of λ = /L and hence can be reduced through orbit phase adjustment.
Significant differences between the error amplification factors of the LS and FD estimators are noticed only for larger values of ε. To illustrate the effect, Fig. 4 shows for ε = 0.5 the error amplification factors of the FD (left panel) and the LS (right panel) planar gradient estimators, both normalised by the mean square inter-spacecraft distance. Here the LS estimator turns out to be more robust than the FD estimate with respect to error propagation.
Another possible source of error concerns the assumption of quasi-static structures. With three spacecraft, planar gradients can be obtained also for spatial structures that undergo at least slow temporal changes. In the case of planar four-point configurations realised with only two spacecraft a and b, however, the velocities V a and V b of the spacecraft relative to the spatial phenomenon of interest must be known. Using data from the three-spacecraft Space Technology 5 mission, Wang et al. (2009)  along the ambient magnetic field with speeds in the range of 100 km s −1 . Attributing the resulting change in the measured magnetic field to the spatial structure could lead to much larger errors, therefore, efforts should be made to identify Alfvénic structures using both magnetic and electric field data.

Accuracy of normal derivative estimation
The use of constraint equations in the reconstruction of normal derivatives means that an additional source of error comes into play, namely, when the constraints are not satisfied exactly. For a detailed discussion of the two geometrical constraintsê ⊥ ∇f andê ± ∇f , see Vogt et al. (2009). In brief, the angle γ between the normal directionn and the true gradient ∇f is found to be of key importance. Furthermore, error indicators that should not become too small are |ê ×n| for the parallel constraintê ± ∇f , and |ê ·n| for the perpendicular constraintê ⊥ ∇f . For the force-free constraint, the error introduced by a nonzero force j × B is given in Appendix C. The error magnitude assumes its maximum value µ 0 j ⊥ /n ·B when the force vector is in the spacecraft plane. Here j ⊥ is the current component perpendicular to the magnetic field. For Swarm, the denominatorn ·B will not be far from unity in the auroral zone which means that here the method should be reasonably robust against errors. At low latitudes, normal derivative estimates should be more cautiously interpreted.

Nonlinearity and sub-scale structures
The first three terms of the Taylor expansion of a twodimensional scalar field can be written as follows: where the first two terms constitute the corresponding linear field, G = ∇ p f is the gradient, and H is the Hessian matrix of second-order derivatives. Without loss of generality, we may identify the reference position vector r • with the origin of the mesocentric coordinate frame, thus r • = r * = 0. Linear gradient estimators are called consistent if they reproduce the true gradient G = ∇f when applied to a general linear field. In canonical form with base vectors q σ , the condition is equivalent to σ q σ = 0 (null vector) and σ q σ r T σ = 1 (identity matrix). In Appendix A it is shown that linear LS estimators are consistent by construction. For the FD gradient estimators discussed in Sect. 3.3, the conditions are easily verified in (u, v) coordinates. Consistency implies at least first-order accuracy in the sense that the estimated gradient differs from the true gradient only in terms that are linear in the set of discretization parameters such as the separation lengths L and M. Second-order accuracy is achieved when a gradient estimator eliminates also the quadratic term in the Taylor expansion. For the four-point FD and LS estimators of Sect. 3.3, the algebra is straightforward but very lengthy. The final result for the FD estimator can be written as with f uu = ∂ 2 f/∂u 2 , f uv = ∂ 2 f/∂u∂v, and f 2 (r) = 1 2 r T Hr. Using a curvature length scale L c to approximate the second derivatives f uu , f uv ∼ f * /L 2 c , we find that, although second-order terms do not cancel exactly, they are proportional to O(ε) and thus small if the discretization parameters (L, M, , m) do not exceed the curvature length scale L c . For Swarm the discretization scheme is thus effectively secondorder as long as we are not concerned with inhomogeneity scales that are smaller than the spacecraft separations, i.e., with sub-scale structures.

Example
The gradient analysis framework proposed in this paper is tailored to the two Swarm satellites orbiting the upper ionosphere and auroral zone in an east-west side-by-side configuration. To demonstrate the analysis methods, one would ideally use multi-spacecraft FAC observations from such a configuration in the LEO environment with its high ambient magnetic field. Unfortunately, past and ongoing multispacecraft missions like Cluster (Escoubet et al., 1997) or Space Technology 5 (e.g., Slavin et al., 2008) produced crossings of the auroral zone in string-of-pearls configurations. To test the different types of constraints for computing normal derivatives, we need to consider current structures in other geospace regions that share the main characteristics of auroral current sheets, i.e., the structures should be planar, force-free, and stationary in a well-defined frame of reference. As shown in Appendix D, the force-free condition is satisfied at a planar plasma boundary characterised by the profile P , normal unit vectorŝ, and speed U , if both the magnitude change |P | and the normal component P s are sufficiently small. Directional discontinuities in the interplanetary magnetic field typically meet these requirements, see Knetter et al. (2004). Cluster FGM measurements at such a plasma boundary are displayed in Fig. 5. The transits of S/C 2 (red) and S/C 3 (green) happen at approximately the same time and (in the reference frame moving with the structure) in a side-by-side configuration, so that the event is well suited for the demonstration of two-spacecraft gradient estimators.
Information from all four Cluster spacecraft are used to determine the key boundary parameters and their errors through crossing time analysis as described by Vogt et al. (2011). We obtain the following values for boundary speed and orientation: The accuracy of these estimates is controlled by the geometry of the Cluster tetrahedron as given by its inter-spacecraft distance L = 4190 km, planarity P = 0.51, and elongation E = 0.45, as well as by the timing accuracy δt ∼ 0.1 s. The resulting uncertainties in boundary orientation and speed can be obtained using the error formulas and Figs. 1-3 in Vogt et al. (2011), and they turn out to be 2 • and 3 %, respectively. Hence the reference frame in which the plasma structure appears stationary is indeed well defined. Furthermore, the similarity of the four boundary profiles as displayed in Fig. 5 strongly suggests that the structure is planar on the inter-spacecraft separation scale. Figure 6 shows that the field magnitudes vary only little during the whole time interval, and the normal components B s are small, in particular until about 01:38:55 UTC, so that we can assume the overall structure to be force-free, see Appendix D. A more refined analysis based on the decomposition of the boundary curl profile into components parallel and perpendicular to the local magnetic field vector is presented below, see Fig. 8. It turns out that a small portion of the boundary is not force-free, with significant effects on the corresponding curl estimator. The boundary profile P = P (t) superposed in Fig. 5 was obtained through averaging the shifted time series B σ . The derivativeṖ = dP /dt (computed through finite differencing) and the boundary slowness vector m =ŝ/U yield the following estimator for the curl of the magnetic field: see Vogt (2014) and Appendix D. The three components of the resulting boundary curl profile are the magenta lines in Fig. 7, as they should be observed at the mesocenter of the two spacecraft S/C 2 and S/C 3. After entering the boundary at 01:38:52 UTC, the curl increases very rapidly for about 0.3 s, followed by a much slower decrease. Displayed in the same Fig. 7 are curl estimates produced by the two-spacecraft LS scheme using FGM measurements from Cluster S/C 2 and S/C 3. Three different constraints were applied to compute the derivatives normal to the plane spanned by the spacecraft separation vector and the velocity vector, namely, (a) the gradient is assumed to be parallel to the boundary normal unit vectorŝ (red), (b) the gradient is assumed to be perpendicular to the (local) magnetic field vector obtained as the average of the measurements at the two spacecraft (green), and (c) the force-free condition is assumed (blue). The three constraints yield almost identical results for the GSE y component of the estimated curl because this component is determined by derivatives in the four-point plane, i.e., by in-plane gradient estimates that are not affected by the conditions to construct the out-of-plane derivatives. Note that the spatial separation of S/C 2 and S/C 3 was mostly in GSE z direction, and the boundary was moving essentially along the x axis, so the unit vector normal to the four-spacecraft configuration was mainly along the y axis.
The time interval chosen for the analysis was t = 0.4 s (equivalent to nine samples), so the differences between leading and trailing positions in the quasi-four-point configuration were 2 t = 0.8 s. Assuming an instrumental error of the order δB 0.5 nT, the configuration yields geometrical errors for the planar gradient of about δB · σ |q σ | 2 1 nT Mm −1 . Shortly before and after the rapid increase at 01:38:52 UTC, actual deviations of the two-spacecraft curl estimates from the boundary profile are larger. The y components of the estimated curl profiles are clearly smeared out because the rapid increase is not resolved by the finite differencing window of 0.8 s. For the geometrical (parallel and perpendicular) conditions, this effect is also present in the z component of the curl. The overall shape is reproduced but the resulting peak values turn out to be lower than in the curl estimate derived from the four-spacecraft averaged boundary profile.
Unlike the geometric conditions, the curl estimate produced with the force-free constraint yields a z-component at 01:38:52 UTC that is far too large to be explained by differencing scale effects. Figure 8 shows that the force-free condition is not satisfied at that time. The curl component parallel to the local magnetic field (magenta) is much larger than the perpendicular component (black) only about ±0.7 s around 01:38:53 UTC, and during this time interval the z component of the two-spacecraft estimate is indeed close to the boundary curl profile.

Summary
In preparation of the three-satellite mission Swarm where two spacecraft move on neighbouring orbits, this paper provides a generic framework for gradient and FAC estimation. Measurements in the plane spanned by either three spacecraft or by two spacecraft and the orbital velocity are utilised to estimate the planar gradient component. The normal component is constructed in a second step from a suitable geometrical or physical constraint in algebraic form. Both steps are independent in the sense that planar gradient estimators and constraints for the construction of normal derivatives can be freely combined.
For probing quasi-static structures by means of four-point configurations generated by two spacecraft on neighbouring orbits, available FD estimators of the planar gradient were brought into a canonical form, and compared with a new LS planar gradient estimator. Comparison of the canonical base vectors revealed that the FD estimators proposed by Ritter and Lühr (2006) and Shen et al. (2012) are algebraically identical. Furthermore, they yield the same numerical approximation of electrical currents as the boundary integral approach. The new LS planar gradient estimator was shown to differ from the FD estimators only in terms proportional to the parameter combination ε = m/M which is of the order 10 −2 for configurations generated by the Swarm satellites a and b. We studied statistical errors associated with measurement noise, and discretization errors induced by deviations of the actual field from a linear model. The results suggest that the adjustable time interval t should be chosen to make the resulting along-track separation M of measurements comparable with the across-track discretization length L given by the separation of the two satellites. For the Swarm mission where L changes substantially along the satellite orbits, the time interval t should thus be continuously varied in accordance with L. Furthermore, orbit phase adjustment of the two satellites through a time shift τ can significantly reduce the errors of the planar gradient estimators studied here. Note that such an adjustment requires the quasi-static assumption to be valid over a longer time scale.
For the construction of normal derivatives we considered three types of constraint equations, namely, the two geometrical conditions of Vogt et al. (2009), and the combination of force-free and divergence-free conditions as introduced by Shen et al. (2012). The latter was rewritten to ease comparison with the geometrical constraints, and to analyse its robustness against errors when the underlying assumptions are not exactly met. In practice, the predictions of different constraints may be compared among each other and with ground-based instruments, thus facilitating consistency checks and error analyses. This approach was illustrated for a transition of the Cluster tetrahedron through a planar plasma boundary in the solar wind. A reference curl profile was constructed from four-point crossing time analysis. Two-spacecraft curl estimates were produced for three different types of constraints. As expected, large errors can result when constraint equations are not satisfied. Furthermore, it was demonstrated how the finite spatial resolution of the satellite array affects curl estimates when sub-scale structures are present.

Least-squares planar gradient estimator for four-point configurations
In mesocentric (u, v) coordinates, the position tensor is the (2 × 2) matrix: The inverse of the matrix can be written as follows One compares with Eqs. (39) and (40) to find that the LS approach to planar four-point gradient estimation differs from the FD approach discussed in Sect. 3.3 only in terms proportional to ε: Note that M is the main requirement for this ordering to be meaningful. The condition may hold even if L ∼ M is not met, e.g., when for the study of thin elongated structures one chooses a small time interval t to improve the along-track resolution.
The algebra of LS gradient estimators further yields which implies that linear LS gradient estimators are consistent by construction, see also Sect. 4.3.
To prove the formula for the condition number of a (2 × 2) matrix given in Eq. (34), we denote their eigenvalues as λ 1 and λ 2 and make use of the fact that the determinant D = λ 1 · λ 2 and the trace T = λ 1 + λ 2 are invariant under coordinate transformations. We then form where c = λ 1 /λ 2 is the condition number. With the (natural) logarithm x = ln c we may write which finally gives and thus Eq. (34).

Equivalence of FD and BI curl estimation schemes using planar four-point configurations
The planar four-point configuration sketched in Fig. 1 naturally suggests to evaluate the boundary integral B · ds by means of the trapezoidal rule on each of the four legs: This yields sixteen dot products of which eight cancel. The remaining eight terms can be rearranged to yield Inserting Eqs. (19) and (20), we obtain and thus for the discrete version of the boundary integral The oriented area A = An is composed of contributions from two triangles: Using again Eqs. (19) and (20), the expressions can be rearranged to yield thus A = A 1 + A 2 = 4 r × V * t = 4 LMn. The modulus is A = 4 LM. The resulting curl estimator is identical with the combination of partial derivatives Appendix C

Force-free and divergence-free constraints
As demonstrated by Shen et al. (2012), the conditions allow to construct the normal derivatives ∂/∂n =n · ∇ of the magnetic field B from the planar gradient ∇ p B. We now show that the resulting constraint equations can be written in the form given by Eq. (51) that facilitates comparison with other constraint equations.
Like any other vector, the normal derivative ∂B/∂n can be decomposed into its planar component (perpendicular ton) and its normal component (parallel to ±n) as follows:
The divergence-free condition 0 = ∇ · B = ∇ p · B +n · ∂B ∂n (C4) immediately gives the normal component in terms of ∇ p · B. The planar component is found from the force-free condition using ∇ = ∇ p + ∇ n = ∇ p +n (∂/∂n) and ∂/∂n =n · ∇: Cross-multiplication withn yields where the identity a × (b × c) = b(a · c) − c(a · b) was used.
Another cross-multiplication withn gives the final form of the planar component of ∂B/∂n: The derivation further demonstrates which error results when the observed structure is not exactly force-free. We note µ 0 j = ∇ ×B, measure the mismatch through t = ×B, and obtain for the resulting error in ∂B/∂n: where j = |j |. The error is zero only in the exceptional case of the force vector j × b and thus also t being parallel to the normal direction, i.e., when both the electrical current and the magnetic field are in the spacecraft plane. The error is maximum in magnitude when the force vector itself is in the spacecraft plane, then where |t| = | ×B| = sin (j , B), and j ⊥ = j |t| is the current component perpendicular to the magnetic field. and the normal component P s are small. Measurements at such a boundary can be modelled as B σ (t) = M(t, r σ ) + δB σ (t) = P (t − m · r σ ) + δB σ (t) , see Vogt (2014). Here m =ŝ/U is the boundary slowness vector, and the terms m · r σ are identical with the lags in crossing time analysis. The spatiotemporal model M(t, r) = P (t − m · r) yields the following estimator for the curl of the field: Note that thus the first term is proportional to the change in intensity. The second term is proportional to the normal component because m·P = mP s . Therefore, planar structures with B s close to zero and sufficiently small variations in magnetic intensity can be considered approximately force-free (j × B ≈ 0).