GALS – Gradient Analysis by Least Squares

We present a method, GALS (Gradient Analysis by Least Squares) for estimating the gradient of a physical field from multi-spacecraft observations. To obtain the best possible spatial resolution, the gradient is estimated in the frame of reference where structures in the field are essentially locally stationary. The estimates are refined iteratively by a least squares method. We show that GALS is not very sensitive to the spacecraft configuration and resolves structures much smaller than the characteristic size of the spacecraft distribution. Furthermore, GALS requires little user input. GALS has been tested on synthetic magnetic field data and data from the Cluster FGM instrument. GALS will also be useful for other types of data. The results indicate that GALS is robust and superior to the curlometer method for estimating the current from magnetic field measurements.


Introduction
The fundamental equations of space physics, such as the MHD equations, Maxwell's equations, and the Vlasov equation, are all first order differential equations relating the temporal evolution to spatial gradients of physical fields. To compare in situ observations with theory we must hence be able to calculate gradients from measurements. In space physics the calculation of space and time gradients from satellite measurements is not a trivial problem. Using a single spacecraft we can determine how a measured field varies Correspondence to: M. Hamrin (hamrin@space.umu.se) along the satellite orbit, but without complementary information we cannot tell whether these variations are spatial or temporal. Furthermore, we have no information about gradients in directions perpendicular to the spacecraft orbit. The unique capability of resolving three-dimensional spatial variations was an important motivation for the Cluster mission, comprising four identical spacecraft launched in 2000. Descriptions of methods for analyzing multi-spacecraft data, for example the curlometer method, the wave telescope technique, and the discontinuity analyzer, are collected in Paschmann and Daly (1998). In this study we will compare GALS to the curlometer and to a single spacecraft method (Lühr et al., 1996).
From four simultaneous measurements in space we can obtain reasonable estimates of the spatial derivatives. For example, the curlometer method (Robert et al., 1998) estimates the rotation of the magnetic field, ∇×B, and (neglecting the displacement current) the corresponding current density j, can be calculated. However, the curlometer cannot estimate spatial variations on length scales smaller than the spacecraft configuration. De Keyser et al. (2007) recently described a method based on least squares fitting for calculating gradients from multi-point observations. This method is demonstrated to be very robust, and it can provide reliable error estimates.
Gradients with scalelengths significantly shorter than the distance between spacecraft can be resolved by combining a discontinuity analysis (Dunlop and Woodward, 1998) to determine the orientation and velocity of the boundary with single spacecraft techniques (Lühr et al., 1996) to compute gradients. The spatial resolution of these methods is determined by the data sampling frequency and the velocity of the spacecraft relative to the discontinuity. These methods can produce very good results under favorable conditions, but their application requires significant effort and skill. This paper presents a method called GALS (Gradient Analysis by Least Squares), which in many cases can    estimate spatial gradients on length scales shorter than the characteristic size of the spacecraft configuration. Since GALS is based on weighted least squares similar to those used by De Keyser et al. (2007), it is possible to produce error estimates of comparable quality. However, in this first study we will not discuss error estimates. Instead, we will focus on describing how GALS combines the best features of the curlometer and the De Keyser et al. (2007) methods with the high resolution of the discontinuity analysis/single spacecraft techniques. The gradients of scalar physical fields are invariant under Galilei transformations. A simple consequence of Galilei invariance is that there almost always exist a frame of reference where the local time derivative is zero. The only exception is the case when the spatial gradient vanishes. In that particular case the time derivative cannot be transformed away.
The quality of the gradient estimates that can be obtained from spacecraft observations of the field depends crucially on the choice of frame of reference. The best possible estimate is obtained in the frame of reference where the field locally appears stationary. Notice that "stationary" here means that the partial time derivative of the physical field is zero at the point in space-time where the gradient is computed. In this reference frame, time variations at nearby points are also minimized, although they cannot be completely eliminated. As will be explained below, the usually rather high time resolution obtained from in situ satellite measurements may in this optimal frame of reference be converted into a spatial sampling along the gradient at distances much shorter than the satellite separation.

Method
Let the position of satellite s at time t be given by r s =R s (t). In the following, the satellite subscript s will be omitted unless it is needed for clarity. Assuming a cluster of S satellites, we will estimate the space and time gradients of a physical field F along the trajectory R cm (t)= 1 S S s=1 R s (t) of the center of mass of the satellites. At the reconstruction time t p the gradients are determined at the point r p =R cm (t p ). Transforming to a frame of reference moving with the velocity P r p =d t R cm (t p ) of the center of mass, we introduce the new coordinates Here, τ is the time relative to the reconstruction time and x 0 (τ ) is the corresponding spatial coordinates relative to the reconstruction point r p moving with velocity P r p . Using these substitutions, the function F (r, t) describing the field in the original coordinate system is in the new system replaced by F 0 (x 0 , τ )=F (x 0 (τ )+r p +τP r p , t p +τ )=F (r, t).
The first derivatives of F 0 form a linear approximant F 0 L (x 0 , τ )=F p +x 0 ·∂ 0 x 0 F 0 p +τ ∂ τ F 0 p to the exact F 0 , valid for small x 0 and τ . Considering a reference frame moving with velocity u relative to the satellites and introducing the new coordinates x u =x 0 −uτ , we find that the approximant is transformed into where F u p =F 0 p =F p and ∂ x u F u p =∂ x 0 F 0 p =∂ x F p are independent of the reference frame, but the time derivatives are related by ∂ τ F u p =∂ τ F 0 p +u·∂ x 0 F 0 p . From this we see that by choosing a velocity we can always find a reference frame where ∂ τ F u p =0. As illustrated in Fig. 1, we can convert the temporal resolution of the observations into a spatial resolution by transforming to a reference frame where ∂ τ F u p =0. Figure 1a shows the situation in the center of mass frame of reference. The satellites sample the field at the discrete times τ n =n t, n=0, ±1, . . . , ±N. At the time τ n satellite s is at x 0 sn =x 0 s (τ n ) and observes the field F sn =F 0 (x 0 sn , τ n ). The spatial coordinates x 0 sn of the satellites are essentially independent of τ n , and take on only four distinct values. While the satellites are at rest, a point of constant F 0 L moves to the right with velocity u, that is, F 0 L (τ u, τ )≈F p . Figure 1b illustrates the same scene in the reference frame moving with velocity u. At x u =0 the value of F u L is now time independent and F u L (0, τ )≈F p to first order in τ . Since the field is as stationary as possible, the distortion caused by time averaging is minimized in this reference frame. The satellites moving through the structure now provide measurements at the coordinates x u sn =x 0 sn −uτ n , with a spatial separation u t along the gradient direction. In many cases, this will allow us to resolve the gradient with higher resolution and better statistics than using only four points with separations determined by the satellite configuration.
The field and its derivatives at the reconstruction point are estimated by minimizing the weighted least squares problem In this first study, this problem is solved by a simple iterative method, described in the Appendix. The Maxwell equation ∇·B=0 is used as constraint in our least squares routines. The weight function W is initially an isotropic Gaussian whose width is determined by the characteristic spacecraft separation L, but its width along the gradient will be adjusted in the algorithm to a value <L whenever possible.
The GALS method is designed to optimize the resolution in the direction of the gradient by transforming to a reference frame where ∂ τ F u p =0 and reducing . The resolution in the plane perpendicular to the gradient is still determined by the satellite separation, but since this is a plane of minimal field variations the resolution is not critical.
Structures in the observed field will in practice have a finite lifetime. By specifying a parameter, which we call the coherence time T c , the user can to some extent tune GALS to the problem and the data at hand. In the algorithm, a time window is created by restricting the τ -values to τ n ≤T c /2. The optimal choice will depend on the data at hand, but also on the physical phenomenon of interest. If T c is chosen longer that the actual coherence time of relevant structures in the field, it may be impossible to find a meaningful velocity. On the other hand, if T c is chosen too short, the GALS algorithm may not be able to interpret data from separate spacecraft as a single, coherent, structure. For a particular feature in the data, it is often possible to estimate a suitable T c by inspecting when and how well this feature is reproduced in the data from the other spacecraft. In practice, it often makes sense to use the rule of thumb that T c is chosen to match the "cluster transition time", i.e., the approximate time from when the first spacecraft enters the structure until the last spacecraft leaves it.

Results
In this section we present the capacity of GALS as an estimator for the current density. The performance of GALS is investigated by using both synthetic and real magnetic field data, and the GALS results are compared with what can be obtained with two other tools for estimating the current: the curlometer (Robert et al., 1998;Dunlop et al., 2002) and the single-spacecraft technique (Lühr et al., 1996).
The curlometer uses multi-spacecraft, e.g. Cluster, magnetic field data to calculate the current density from ∇×B (assuming that the displacement current is negligible in Ampere's law). The current is estimated from four simultaneous measurements of the magnetic field, and the curlometer therefore captures the instantaneous current density as a snapshot in time. Since the curlometer technique is based on the assumption that the current is constant over the tetrahedral volume defined by the four spacecraft, the spatial resolution of the curlometer current density depends directly on the satellite configuration.
The single spacecraft technique can very efficiently resolve, e.g., thin current sheets when it is combined with a method such as the discontinuity analyzer for obtaining the velocity of the spacecraft relative to the sheet. Considering the sheet to be stationary, the observed total derivative is interpreted as a spatial gradient and the resolution is determined by the data sampling rate rather than by the size of the satellite tetrahedron. The quality of the velocity estimate is of crucial importance when determining the full current density vector with the single-spacecraft technique. When analyzing a convecting, planar, current sheet passing all four spacecraft this is generally not a problem, but investigations of more general current density structures with the single-spacecraft technique are limited by the lack of good velocity estimates.
The least squares method for gradient calculation developed by De Keyser et al. (2007) is more robust than the curlometer, but they did not explicitly exploit the advantages of a moving reference system. In the absence of significant correlations between the spacecraft, GALS will calculate gradients by a method very similar to that of De Keyser et al. (2007). However, whenever a coherent structure can be found, GALS will automatically attempt to determine its velocity and and estimate the gradients by a method related to the single spacecraft technique. We will show that GALS combines the robustness of the least squares method with the high spatial resolution of the single spacecraft technique. We will also show that it is possible to tune the behaviour of GALS to the problem at hand by changing the coherence time, T c .
The setup for a test run of GALS on synthetic magnetic field data is shown in Fig. 2. Four identical spacecraft observe a Gaussian current sheet that is infinite in the y and z directions. The spacecraft and the current sheet approach each other along the x direction with a relative velocity of 20 km/s. The current sheet width is W =20 km, which is up for a test run on synthetic magnetic field data. Relaed reference system, the spacecraft and the structure are h +5 km/s and −15 km/s, respectively. the curlometer (dark green), yellow, blue, and green). Th dotted line. (d) GALS estim (red, green and blue) for the B ter Λ(B y ) obtained by GALS lighted time window indicate Fig. 2. Setup for a test run on synthetic magnetic field data. Relative to a fixed reference system, the spacecraft and the structure are moving with +5 km/s and −15 km/s, respectively. small compared to the characteristic size of the satellite tetrahedron, L=200 km. The tetrahedron has elongation and planarity E=P =0.75 (Paschmann and Daly, 1998).
Panel (a) of Fig. 3 shows the B y component of the magnetic field sampled by the four spacecraft. The other field components are identically zero, B x ≡B z ≡0. The solid magenta line and the dotted black line show the field reconstructed by GALS and the true magnetic field at the center of mass of the satellites. It is evident that GALS reconstructs B y correctly although the current sheet is narrow and could be difficult to resolve. The coherence time T c is chosen as 12 s, as indicated by the width of the highlighted yellow region.
Panels (b) and (c) show the two non-zero components of the current density estimates, J x and J z , obtained from GALS, the curlometer and the single-spacecraft method, respectively. As indicated in Fig. 2, the true current is along the z-direction. However, the curlometer (dark green line) produces an artificial J x , which is bi-polar and comparable to the true current in magnitude. The corresponding J x from GALS (magenta) is three orders of magnitude smaller and not visible in Fig. 3b. In panel (c) we see that GALS (magenta) correctly identifies the magnitude and position of the true current sheet (black dotted line). The single-spacecraft method, utilizing the known relative velocity of 20 km/s, also clearly identifies the narrow current sheet for each satellite crossing. In contrast, the curlometer cannot properly resolve this narrow current sheet, and underestimates the peak current density by more than a factor of two.
In Fig. 3d we see that GALS identifies the structure in the B y data and obtains the correct velocity (−20, 0, 0) km/s with respect to the satellites in the neighborhood of the current sheet. Further away the magnetic gradient is weak, and the velocity is not well defined. Figure 3e shows the resolution in the reconstruction of B y and how it adapts to the data. In the region of the strong magnetic field gradient (B y )≈4 km, indicating that GALS resolves the structure on scale lengths much shorter than L. The exact value of the resolution is determined by details in the GALS algorithm. Since the results presented here are based on a prototype code, the detailed behaviour of the resolution will not be discussed further in this article.
Next we will illustrate GALS capacity to resolve structures on various scale lengths by using real magnetic field data from the FGM instrument (Balogh et al., 1997) on board the Cluster satellites. Figures 4 and 5 show a Cluster magnetopause crossing with a thin current sheet on 30 March 2002. This magnetopause crossing has already been investigated in some detail in the literature (De Keyser et al., 2005;Panov et al., 2006). In this paper we present this sample event just to show the capabilities of GALS, and not to investigate any specific physical details of the event.
In panels (a-c) of Figs. 4 and 5, we show the x, y and z magnetic field components as observed by the four spacecraft (red, yellow, blue, green). We see that the main contribution to this current sheet comes from the B z component. The magenta curve in the same panels correspond to the GALS reconstructed magnetic field along the x, y and z directions. Panels (d-f) contain the current estimated with three different methods: 1) The GALS current (magenta); 2) The curlometer current (dark green); 3) The current from the single-spacecraft method (red, yellow, blue, green). The single-spacecraft current is calculated using a current sheet velocity of (−23, −11, −4) km/s, derived by a timing analysis of the encounter of the gradient in B z . The results from each satellite are time-shifted to make the peaks coincide. In the last three panels of Figs. 4 and 5 we show the GALS estimated velocity components u x , u y and u z (red, green and blue) for each magnetic field component (B x , B y and B z ). No pre-processing of the data was performed but the output from all methods were smoothed by a 0.25 s sliding window.
In the GALS reconstruction in Fig. 4 we have used a rather small value of the coherence time, T c =0.6 s, indicated by the highlighted yellow region. The result is that the current from GALS is smooth and similar to the curlometer result. Choosing a small value of T c , only very few observations from a short time interval are included in the GALS least squares problem. GALS will then not be able to determine a meaningful velocity, and the resulting current will be similar to the curlometer result, which can be viewed as a snapshot in time of the current along the satellites' center of mass.
In Fig. 5, the same event is analyzed with the coherence time T c =6 s according to our rule of thumb. The coherence time is indicated by the highlighted region. Using a larger   value of T c , more space-time observation points are included in the GALS least squares problem. The extra information from the included time domain is used to improve the spatial resolution. Of course, this requires that the shape of the investigated convective structure is approximately stable on these time scales. In this case, the GALS current will therefore be more similar to the result from the single-spacecraft method, which is highly capable to resolve small scale structures. In panels (d-f) of Fig. 5 we clearly see that the GALS result rather closely follows the single-spacecraft currents. The curlometer, on the other hand, fails to resolve the small scale variations. Notice that the estimated thickness of the current sheet is about 50 km, which is of the order of the proton gyroradius (in the magnetosheath plasma the gyroradius is about 50 km and in the magnetospheric plasma about 25 km).
Panels (i) of Figs. 4 and 5 show the structure velocity obtained in the GALS reconstruction of the B z data which causes the dominant part of the current sheet. Using a large coherence time, T c =6 s as in Fig. 5, we see in panel (i) that the velocity around the time of the current sheet crossing (about 13:11:46.50) is approximately consistent with the result from the timing analysis, (−23, −11, −4) km/s. On the other hand, using a smaller T c =0.6 s as in Fig. 4, the velocity estimate in the B z reconstruction does not stabilize on a specific value during the magnetopause crossing. However, it should be noted that the velocity estimated by GALS is not necessarily closely related to the velocity of a current sheet or some other physical structure in the magnetic field. In many cases it is completely unrelated to the velocity of large scale field structures. As discussed in Sect. 2, the velocity obtained by GALS describes the motion of the frame of reference where the field signal is locally stationary. In that sense, the GALS velocity should rather be regarded as a local phase velocity than as the velocity of a physical structure.     Fig. 3. (a-c) Magnetic field B x , B y , and B z components in the GSE coordinate system observed by the four spacecraft (red, yellow, blue, green) together with the reconstructed magnetic field from GALS (magenta).
(d-f) Current density J x , J y , and J z components from GALS (magenta), the curlometer (dark green) and the single-spacecraft method (red, yellow, blue, green). The highlighted time window indicates the used coherence time T c =0.6 s used in GALS. (g-i) The GALS velocity u x , u y and u z (red, green and blue) obtained for the magnetic field components B x , B y and B z .

Summary and discussion
Above we have shown that it is possible to tune GALS to the problem at hand by changing the coherence time, T c . Choosing a small value of T c , the current resulting from GALS becomes similar to that from the curlometer while a larger T c improves the spatial resolution and makes the current more similar to the result from the single-spacecraft method. This demonstrates that GALS is a tool that easily can be applied to large amounts of data and yet provides a resolution comparable to the single-spacecraft method.
We have applied the GALS method to synthetic magnetic field data and data from the FGM instrument on board the Cluster satellites. The results indicate that GALS is a robust method to obtain high-resolution spatial gradients from multi-spacecraft observations. Furthermore, GALS is simple to use since the user needs to specify only a single parameter, the coherence time T c . This makes GALS complementary to the discontinuity analyzer/single-spacecraft technique, which demands substantial user input.
Comparing GALS with the curlometer and the singlespacecraft method, we find that GALS resolves thin current sheets better than the curlometer. GALS can resolve structures on scales considerably smaller than the characteristic size of the satellite configuration, and our tests indicate that it is less sensitive to the spacecraft configuration than the curlometer. The single-spacecraft method and GALS provide similar resolutions, but the single-spacecraft method requires a separate calculation of the velocity of the current sheet.
GALS is related to the method of De Keyser et al. (2007) in the sense that they are both based on a weighted least squares fit to the data. However, De Keyser et al. (2007) does not take advantage of the freedom to choose a frame of reference moving with convecting structures in the field. In the absence of coherent convective structures, GALS will compute gradients on scale lengths determined by the spacecraft separation, similar to the De Keyser et al. (2007) method. However, when a coherent structure is found, GALS will automatically determine its velocity and produce gradients with resolution similar to the single-spacecraft technique.
Current sheets are important in space physics, since they are ubiquitous at boundaries between plasma populations. Here, we have illustrated the principles behind GALS using current sheets, but it should be clear from the design that neither a current nor a sheet is essential. However, the method requires that the surface of constant F is not too strongly curved within the volume spanned by the spacecraft configuration. Without this assumption, which is implicit in the curlometer as well as the single-spacecraft methods, it seems impossible to derive meaningful information about gradients. Moreover, GALS requires that the life time of relevant structures in the field is not too short.
We consider our results promising and we expect that the GALS method will be useful when reliable high resolution estimates of spatial gradients are needed. In the future, we intend to develop reliable error estimates, investigate the capacity to separate spatial and temporal variations, and apply the method to N =4 spacecraft.

Appendix A
The GALS algorithm A1 Initial setup Let τ n =n t, n=0, ±1, . . . , ±N, be the discrete times at which the satellites sample the field. At the time τ n satellite s is at x 0 sn =x 0 s (τ n ) and observes the field F sn =F 0 (x 0 sn , τ n ). The field sampled by the satellites can be expressed as a Taylor expansion where F 0 p and its derivatives at x 0 =0 and τ =0 are five unknown parameters. The field and its derivatives at the reconstruction point are estimated by minimizing the weighted least squares problem where the weight function is spatially isotropic. The matrix κ 0 is where L is the characteristic size of the spacecraft configuration (Paschmann and Daly, 1998) and =1/T c . The coherence time T c is a user-defined parameter that designates the time scale during which the local structure of the field is expected to retain its approximate shape. Figure 1a illustrates the sampling of the field in the satellites frame of reference. The spatial coordinates x 0 sn of the satellites are essentially independent of n (i.e., time), and take on only four distinct values. Although the spatial resolution will be as low as if we used only one observation per spacecraft, minimizing Eq. (A2) will give us a first estimate of F 0 p , ∂ x 0 F 0 p , and ∂ τ F 0 p . These estimates can then be improved by iteration.

A2 Iterative improvements
The initial estimates of the field and its gradients can be improved by choosing a new frame of reference that follows the field structure rather than the satellites. As long as the gradient does not vanish, we can always find a velocity, parallel to the gradient, such that the field at a point moving with that velocity has a constant value. To obtain an estimate of this velocity we introduce for i≥1 a unit vectorû i parallel to the estimated spatial gradient, Assuming u 0 =0, the velocity, relative to the satellites, of the field structure at the reconstruction point is then determined as The velocity u i describes how a plane, tangent to a surface of constant F at (r p , t p ), moves along the gradient. In some simple cases, such as the magnetic field generated by a convecting plane current sheet, u i will be closely related to the velocity of the current sheet. However, in more complicated situations u i may be completely unrelated to the motion of large scale field structures. Notice that u i should not be interpreted as the velocity of a physical object, but is more similar to a locally defined phase velocity. To improve the resolution, an anisotropic weight function is used in the iterations. Introduce two orthogonal unit vectorsv i andŵ i in the plane perpendicular toû i . In this plane the field will be slowly varying, and a resolution ∼L is acceptable. To resolve the gradient as well as possible we focus on measurements taken close to this plane through the origin and within the coherence time (|τ n | T c /2). Hence, we define the weight function W i as The resolution i is mainly determined by the magnitude of the velocity u i =|u i | of the satellites in the structure's frame of reference. As illustrated in Fig. 1b, consecutive observations from each satellite will be separated by a distance t>u i in the direction of the gradient. Thus we can typically choose i ∼ t>u i L and still have measurements with significant weight from all spacecraft. However, if the velocity is low or the coherence time short, so that T c >u i <L, the algorithm will increase i in order to include data from all satellites.

A3 The algorithm
After obtaining the initial estimates according to Eq. (A2-A4), the GALS algorithm can be summarized as follows: For i=1, 2, . . . I. Transform to the frame of reference moving with velocity u i by defining II. Solve the weighted least squares problem n,s F i p + x i sn · ∂ x i F i p + τ n ∂ τ F i p − F sn 2 W i (x i sn , τ n ) to obtain estimates of ∂ x i F i p and ∂ τ F i p in this reference frame.
III. Calculate a refined velocity estimate according to Eq. (A6).
These calculations are repeated until the velocity update falls below a threshold. The final results are then transformed back to the initial frame of reference according to The GALS algorithm is applied to each timestep independently. For the case of magnetic field observations, GALS treats each component, B x , B y , and B z , and it obtains estimates of the velocity u and resolution along the spacecraft orbit for structures in each field component. However, the three field components are coupled together with an additional constraint in our least squares routines by the use the Maxwell equation ∇·B=0.