Multipoint Field‐Aligned Current Estimates With Swarm

Field‐aligned currents (FACs) are a key component of the magnetosphere‐ionosphere system, providing the transfer of energy and momentum between the distant collision‐free magnetospheric plasma and the collisional ionosphere. The European Space Agency's Swarm mission offers a unique opportunity to explore the FACs low‐altitude end by a broad range of single‐spacecraft and multispacecraft techniques. The present technical report demonstrates the application of dual‐ and three‐satellite FAC estimation methods based on the least squares approach, with the goal to enhance the use of Swarm data in exploring the high‐latitude current structures, beyond the official FAC product presently available to the end user. The dual‐satellite method presents some clear advantages since it provides stabler solutions, can be applied on a more general spacecraft configuration, and offers a robust error estimation scheme. Consequently, the method provides significantly more data near the singularity where the satellites' orbits intersect, allows FAC estimation with configurations that involves the upper Swarm satellite, or to fine tune the constellation geometry to the problem at hand. Similarly, the three‐satellite method, meant to be applied when Swarm forms a close configuration, brings additional valuable information, associated to a different (larger) scale. Moreover, it can estimate the FAC density with high time resolution when instantaneous measurements are used. The performance of the methods are thoroughly analyzed both on Swarm events and on simulated data, and the results are compared with other available methods. Particular emphasis is put on how different FAC estimation methods complement each other and provide consistent results.


Introduction
Auroral displays in the Earth's polar region are generated by energetic particles entering the ionosphere along geomagnetic field lines and offer a vivid impression of the complex plasma dynamics in near-Earth space. Electromagnetic coupling to remote magnetospheric regions such as the magnetotail and the magnetopause is accomplished by field-aligned currents (FACs), that is, electric currents flowing along magnetic field lines (e.g., Lysak, 1990;Paschmann et al., 2002;Vogt, 2002). Magnetosphere-ionosphere coupling through FACs is an essential ingredient in studies of solar-terrestrial relationships and space weather modeling (e.g., Baker et al., 2018;Ohtani et al., 2000).
The distribution and dynamics of FACs is usually derived from the associated magnetic perturbation observed by space-borne magnetometers (see, e.g., Zmuda & Armstrong, 1974;Iijima & Potemra, 1976; and more recently, He et al., 2012He et al., , 2014. Single-satellite measurements (e.g., Freja, FAST, and Champ) cannot reveal the spatiotemporal variability of current structures in geospace and therefore the corresponding FAC density estimation methods rely on implicit assumptions regarding time dependence and geometrical structure (Fung & Hoffman, 1992). In general, stationary infinite planar current sheets, aligned along the geomagnetic east-west direction, are postulated. A more careful analysis takes into account the actual FAC orientation, derived using minimum variance analysis (MVA; Sonnerup & Scheible, 1998) applied to the magnetic measurements. Typical velocities of FAC sheets, in the range of a few hundred meters per second (Wang et al., 2009), are small compared to the spacecraft speed (∼7.6 km/s for a LEO satellite at 500 km altitude) and are usually neglected in such analysis schemes.

10.1029/2018JA026439
orbital phase, to ensure as much as possible the application of Ampère's law on a rectangular quad; the quads are formed by combining virtual positions, that is, two points from each satellite, 5 s separated along the orbit; (iii) the discretized form of the Ampère's law, in the BI formulation, is used to infer the current flowing through the quad and the FAC density is estimated by considering the orientation of local magnetic field. In this article, besides the L2 FAC product, we present results obtained from an in-house implementation of the dual-s/c BI method that follows the same procedure.
Note that by the way it combines the satellite measurements (shifting one sensor data 5-10 s to compensate for orbital lag and using virtual positions 5 s apart along the orbits), the dual-s/c method relies on the assumption of quasi-stationary magnetic field. At mesoscales, the FAC profiles recorded by different satellites with a similar time delay, tend to be well correlated (Gjerloev et al., 2011;, suggesting time stationarity, but this is not a priory guaranteed (see, e.g., Trenchi et al., 2019).
The dual-s/c method from Vogt et al. (2013), to be used in this paper, offers some advantages over the traditional L2 FAC algorithm. One is the least squares approach, which yields more stable solutions for quad configurations that are far from regular (e.g., the extremely elongated configurations close to the geographic poles), when small perturbations in the magnetic field (noise) have a dramatic influence. For such geometries, the advantages of LS technique boils down to the overdetermined character of the problem, that is, three, instead of four, points of observations would suffice to estimate the current density along the normal direction (more details on that are provided in section 3.1). Then, the method allows for a more flexible geometry, where deviations from parallel satellite orbits as well as skewed configurations are tolerated. A robust framework to asses the level of FAC estimation errors as a function of configuration geometry and the level of magnetic noise is also provided. Brief presentations on the dual-s/c method of Vogt et al. (2013) and on the associated error estimation scheme are given in Appendices A and B, respectively.
The three-s/c gradient estimation method (Vogt et al., 2009) has been developed for analyzing data provided by ESA's Cluster four-satellite mission when only three sensors are available. In principle, three (noncollinear) points of observation are enough to infer the gradient components in the plane spanned by the satellites positions. For that, a similar formulation of the problem as in the case of four-point planar configuration can be employed. The method allows to estimate the missing gradient component, that is, along the normal direction, if additional constraints of geometrical (gradient parallel/perpendicular to a certain direction) or physical (e.g., currents flowing only along the magnetic field direction, time stationarity) nature can be used. Note that since curl is a planar operator, three points of observation are enough to fully obtain its component along the normal direction, that is, the normal current density, in case one works with magnetic field data. The method assumes linear spatial variation of the field at the scale of spacecraft constellation and therefore is meant to be applied when the satellites are in close formation. The assumption of time stationarity is not needed, unless one wants to combine satellites virtual positions (e.g., by time-shifting measurements provided by one satellite in order to obtain a more favorable spacecraft configuration).
For some events we provide results from the single-s/c FAC estimation method. The method of Ritter et al. (2013), used to generate the Swarm L2 single-satellite FAC product, applies the discrete version of the Ampère's law to the magnetic field perturbation observed by one satellite, assuming the model of static, infinite planar current sheet oriented perpendicular to the orbit: Along-track variations of magnetic time series are projected onto directions 45 • left and right of spacecraft velocity vector to yield Δb x and Δb y , followed by finite differencing in time (resolution 1 s). The results presented in this article were obtained from an in-house implementation of the method that offers some advantages, like the possibility to use low-or high-resolution magnetic data, filtered or not, as well as to take into account the current sheet inclination with respect to satellite orbit. To obtain the orientation of the current sheet we applied the MVA (see, e.g., Sonnerup & Scheible, 1998) on the magnetic field perturbations perpendicular to the main magnetic field. Figure 1 presents a comparison of FAC density estimates obtained by applying different dual-s/c methods on two auroral crossing events by Swarm that occurred on 4 May 2014 around 19:16 UTC and, respectively, on 22 April 2014 around 05:46 UTC. The events were chosen such that part of the FAC structure is encountered at the highest latitude portion of the orbit, where the quad geometry assumes a highly stretched configuration and the FAC estimations become unreliable (more on this below). In panels (e), different FAC estimates correspond to the LS method (blue line), the BI method implemented in-house (red line, mostly under the blue line), and two dual-s/c L2 estimates available at the ESA server, that is, from the old baseline 0101 (black dashed line) and the more recent baseline 0301 (magenta line).

General Behavior and Comparison With Other Methods
Apart from the regions close to the geographic poles, all the methods provide basically the same current densities. The slightly different amplitude and frequency composition of the L2 estimate, on one hand, and of the LS/BI estimates, on the other hand, can be explained by the different low-pass filter used to preprocess the magnetic data. In this context Vogt et al. (2013) have shown that estimates of current density from the LS method and from other available dual-s/c methods, like the FD method (Ritter & Lühr, 2006;Shen et al., 2012), and BI method (Ritter et al., 2013), are algebraically identical for nearly parallel orbits and small differences in the orbital phases of swA and swC (i.e., up to the first order in the parameters and from Figure A1 in Appendix A).
The baselines 0101 and 0301 of the L2 FAC product were obtained from different versions of the L1b magnetic field files. However, according to our analysis, this change in data calibration explains only partially their large discrepancy at highest latitudes. The fact could indicate a dependence of the L2 algorithm on the magnetic field model used to compute the magnetic perturbations (the 0101 baseline was generated at the beginning of the mission, using outdated magnetic field models, whereas the 0301 baseline relies on updated, Swarm-based, magnetic field model). Yet, since the magnetic field models are potential fields, and because the Swarm lower satellites are close enough one to another to rely on the assumption of linear field variation when the Ampère's law is applied, one should basically obtain the same current estimations regardless of the particular magnetic field model being used or even when the Swarm measurements are directly entering in the calculations. To test that on the dual-s/c LS method, we used magnetic perturbations based on outdated magnetic field models as well as raw Swarm magnetometer readings (no subtraction of model magnetic field) and found that, indeed, the FAC estimations are very close to the solution one obtains when an appropriate, that is, updated, magnetic field model is used. The small discrepancies were typically below In regions close to the geographic poles, the quad configuration becomes singular, that is, the cross-track separation can decrease, in principle, to zero, while the along-track distance remains constant to roughly 40 km. As a consequence, the FAC intensity starts to take artificially high values and a region should be defined, that is, the so-called exclusion zone (EZ), where current density estimates are not provided. For the dual-s/c L2 product, the region above 86 • geographic latitude in both hemispheres has been chosen to serve for this purpose (Ritter et al., 2013). In case of the LS solution, the EZ extension can be automatically tailored to each event based on the estimated level of errors in FAC density (details are provided below and in section 3.3).
It has been mentioned in section 2 that using four sensors in a planar configuration to infer the current density along the normal direction poses an overdetermined problem since three (noncollinear) points are sufficient for that. For the Swarm crossing event that occurred on 22 April 2014, Figure 2a presents a comparison between FAC estimates from dual-s/c LS solution (blue) and two solutions obtained by applying the three-s/c technique on magnetic field perturbation recorded at three apexes of the four-point quad configuration. More precisely, if one designates the leading and the trailing positions of satellites swA and swC by p and m, respectively, the solution shown in green involves the apexes Ap, Cm, and Cp, while the solution shown in magenta involves the apexes Am, Cm, and Ap. The three-and four-point satellite configurations are shown on the bottom part of the figure as projections on the North-East plane of the local NEC frame at three time instances, that is, start time, stop time, and at the middle of the interval. In all projections the spatial scale is the same, indicated in the middle panel.
As can be observed, outside the EZ, the two solutions obtained from the three-s/c method are basically the same as the dual-s/c LS solution, being capable to capture very well both the along-and cross-track variations in magnetic field perturbation (dB). The situation changes a little when no filtering on magnetic field perturbation is applied (not presented): then each of the three-s/c solutions shows a different small-scale structure, with the dual-s/c solution exhibiting an average fluctuation (in the LS sense), but the evolution is the same at greater scales. The small time shifts between the solutions are due to the differences in the position of the mesocenter; indeed, the triangle designated by ApCmCp (corresponding solution in green) has the most advanced mesocenter along the orbit and reaches first a degenerate configuration, reflected in the high values of the predicted FAC current, while the mesocenter of the AmCmAp triangle (corresponding solution in magenta) is trailing along the orbit.
In Figure 2b, the weighted sum J = (J 1 A 1 + J 2 A 2 )∕(A 1 + A 2 ) of the current densities J 1 and J 2 provided by the three-s/c method is shown in black, with A 1,2 the area of the corresponding three-point configuration. Judging in terms of line integral, the quantity at numerator represents the dB line integral for the four-point configuration, since according to our choice of apexes, the contributions along the diagonal Cm-Ap cancel each other being covered in opposite senses. The dual-s/c BI solution from Figure 1, shown here in red line (behind the black line) proves that this is indeed the case. From the same panel one can observe that the Vogt et al. (2013) dual-s/c method, which provides a solution based on the least squares approach to the above mentioned overdetermined problem (blue line), presents a more stable evolution in the critical area near the geographic poles.
In Figure 2c, the single-s/c FAC density estimates based on swA (magenta) and swC (green) data are compared with the dual-s/c LS solution (blue). Panel (d) presents the difference between the magnetic field perturbation seen at swA and swC in the VEC reference frame, defined as follows: unit axis c is pointing to the Earth center, v along ± V m where V m is the quad's mesocenter velocity vector and + or − sign applies for positive and, respectively, negative component of V m along magnetic dipole axis, and e = c × v. Since the single-s/c solutions are capturing only the along-track variation in the east component of dB, one expects major differences with respect to the dual-s/c solution when the cross-track variation of the v component of dB becomes significant. This is indeed the case, see the red component in panel (d). Moreover, the relatively high density FAC structure around 05:46:10, put in evidence by the dual-s/c LS and (in-house implemented) BI methods, but unresolved by the L2 product, is authentic, being caused by a large cross-track difference in the v component of dB.
Panels (f) from Figure 1 show the inaccuracy in current density estimation due to instrumental noise, that is, j B , computed in accordance with the scheme presented in Appendix B. In evaluating the formula (B1), a flat value of 0.2 nT was considered for B (i.e., 0.5 nT instrumental noise, reduced by a factor of 2.5 due to the filtering, see section 3.3). As expected, the error level increases as the satellite approaches the geographic poles due to the decrease in the cross-track distance. By adopting the definition of EZ from section 3.3 (i.e., j B > 0.1 μA/m 2 ; see details there), one obtains the EZ limits shown by vertical black lines in Figures 1 and 2. Figure 2e presents the difference between the dual-s/c LS and BI estimates, while the vertical blue lines mark the region where this difference exceeds the 0.1 μA/m 2 threshold. The difference between dual-s/c LS solution provided by one simulation run and the simulated current (black) is compared with the predicted level of errors in FAC estimation due to instrument noise (red). In (f), many simulation runs are considered, to show that the standard deviation of the difference shown in (e) (black) agrees with the FAC error estimation formula (red). (g) The difference between the dual-s/c LS estimation and the simulated current when an offset in magnetic measurements is present (blue), and a proxy for this difference (red). (h) Similar to (f) but applies on FAC estimates based on filtered magnetic field data. (i) The ratio between the BI and LS FAC error estimations when filtered magnetic field data are used. (j) Differences in current density estimates are investigated for three particular quad configurations shown on the rightmost part of the figure. See text for more details. BI = boundary integral; LS = least squares; FAC = field-aligned current.

Dual-s/c Method Performance on Synthetic Data
The dual-s/c FAC estimation method has been extensively tested on synthetic data to acquire information on its performance in various current configuration. Figure 3 presents the results of one such test, when an array of two Swarm-like sensors is sampling a sequence of upward and downward cylindrical current structures of different scale lengths. The satellite orbits are shown on the right as solid lines that crossed each other, overlaid on the current density contour plot. All the parameters of the simulation (spatial scales, orbits, current intensity, etc.) were chosen to be close to the Swarm context. An artificial noise (randomly distributed values of 0.5-nT standard deviation) has been added to the data to simulate the measurement errors.
Panel (a) presents the vector magnetic field perturbation recorded by swA, that is, dB swA , with x, y, and z components corresponding to the horizontal, vertical, and out of the paper directions for the orbit plot on the right. Panel (b) shows the difference in magnetic field perturbation measured by the two satellites, that is, dB swC −dB swA . In panel (c), the comparison between the LS dual-s/c FAC estimation, that is, j es (red), and the simulated current at the quad's mesocenter, that is, j m (blue), shows a very good agreement, apart from the region where the satellite orbits cross each other (interval between the vertical lines, corresponding to the EZ). The same conclusions apply for the dual-s/c estimation based on the BI method (green), although the differences with respect to simulated current are somewhat larger. Panel (d) compares the single-s/c FAC estimation (red) and the simulated current (blue) at swC position. The single-s/c method does a good Panel (e) shows the difference between the current estimated by LS method and the simulated current, that is, j es − j m (black), together with the errors in FAC estimation provided by equation (B1) from Appendix B, that is, j B (the two lines in red correspond to ± j B ). Within the EZ, the dual-s/c solution is unstable, that is, small fluctuations in dB lead to big changes in j es and the aspect of the black line changes dramatically from one simulation run to the other, depending on the actual values that simulates the magnetic noise. Nevertheless, if one runs the simulation for a sufficiently number of times and compute the standard deviation of j es − j m , one should obtain the behavior predicted by the error formula. This indeed is the case, as shown in panel (f), where the es − m (black), obtained from 1,500 simulation runs, matches very well the j B predicted by equation (B1) (red line, under the black one).
Besides the instrumental noise, the FAC estimations can be affected by the presence of offsets in magnetic field data. To illustrate the effect, we considered noise-free synthetic data and added magnetic offsets in all components of B. In panel (g), the blue line presents the difference between the dual-s/c LS estimation and the simulated current when the offset vector was B off = [0.25, 0.25, 0] nT for swA and − B off for swC. The FAC inaccuracy, that is, j off , takes significant values only closer to the point where the satellites' orbits intersect, with the blue vertical lines indicating the region where the 0.1 μA/m 2 threshold is exceeded. An approximation for j off is presented in Appendix B; the evolution prescribed there by equation (B3) is shown here in red.
Using synthetic data, we further tested how ionospheric electrojet influences the FAC estimations. Figure 4 presents, on the right, the geometry of the simulation: the electrojet is oriented along the X axes, while swA and swC are flying 350 km above, on parallel orbits, inclined with respect to X axes at 30 • . The electrojet has been modeled in two ways, namely, as a straight line current of 800-kA intensity, and as a series of parallel line currents, uniformly distributed over 500 km in latitude, with the same total intensity. Panel (c) on the left of Figure 4 presents the FAC estimated by single-s/c method, when no correction on the orientation of the current sheet is considered: due to the inclined orbit, the satellite observes a cross-track magnetic field perturbation and erroneously attributes it to a FAC. For this rather strong electrojet, the error amounts to ∼0.3 μA/m 2 in the case of the localized electrojet (red) and ∼0.2 μA/m 2 for an electrojet spread in latitude (blue). In the dual-s/c method, where one effectively integrates the magnetic perturbation on a close path, the predicted current density is expected to be vanishingly small, since no current is actually flowing though the quad. This is indeed the case: In panel (d), the dual-s/c FAC density estimations are of the order of 1 nA/m 2 , being actually caused by the breach in the assumption of linear variation for dB. In the analysis above, to clearly show our point, we considered noise-free magnetic measurements recorded by satellites flying on parallel tracks. Nevertheless, we also used nonparallel satellite orbits that have their cross point located above the electrojet to find out that the EZ limits are not affected by this current system.

The Extension of the Exclusion Zone
To delimit the EZ extension, a threshold in current density error, j er , can be used, for example, j er ≤ 0.1 A∕m 2 , or, another possibility, one could work with a certain upper value in the FAC relative errors, for example, j er ∕j es ≤ 10%. Below, the two type of error sources mentioned in section 3.2, that is, instrument noise and offsets, will be discussed separately, to illustrate their specific influence. It is evident that, in reality, both effects are simultaneously present.
Considering first the instrumental noise, one can naturally use the formula for j B provided by equation (B1) from Appendix B. Note that adopting, for example, the first criterion, that is, j B ≤ 0.1 A∕m 2 , does not imply that the FAC density in one particular event is close to the real value by 0.1 μA/m 2 , but rather there is a chance of ∼68.3% to be within this limit, assuming a Gaussian distribution for the errors.
As discussed in section 3.1, within the EZ the four-point satellite configuration becomes very stretched due to the decrease in the cross-track distance. Since the Vogt et al. (2013) method allows for a flexible geometry of the spacecraft array, one could be tempted to optimize the configuration at higher latitude, for example, to reduce the along-track separation (parameter M in Figure A1) in order to probe the EZ region with more regular quads. This would effectively mean to minimize the contribution to FAC estimation error arising from the quad shape, that is, the factor L 2 r ∕L 2 q , in equation (B2). However, reducing the along-track separation would reduce the mean square interspacecraft distance, L 2 r , increasing thus the first term in equation (B2), which captures the inverse dependence of the gradient errors on the array size. Using j B ≤ 0.1 μA/m 2 as definition of the EZ, we tested the net effect of geometry optimization on many Swarm events and arrived at the conclusion that, basically, EZ extension is not changing when the along-track separation is reduced. In connection to this, it is worth mentioning that changing the along-track separation would also imply a scale effect, that is, a change in the spatial resolution of the FAC estimation. We therefore decided to keep the 5-s along-track separation distance along the entire orbit, that is, the same as in the ESA dual-s/c L2 product.
The dual-s/c method relies on the assumption of quasi-stationary magnetic fields; to diminish the influence of small scale fluctuations, always present in real data, a low-pass filter is used in the preprocessing chain (see the discussion in Appendix B). Throughout this paper we worked with a 20-s low-pass Butterworth filter, which provides the flattest frequency response within the pass band. This effectively reduces the influence of instrumental noise on FAC estimation, as shown in Figure 3h. The black line presents the evolution of esf − m , that is the standard deviation of estimated FAC solution obtained from filtered magnetic field data, j esf , with respect to the simulated current at the quad's mesocenter, j m . An overall reduction in the error level is evident, together with a flattening in the region where the satellites' orbits cross (see panel (f) for comparison). The red line designates the error level prescribed by equation (B1) when a random noise of 0.2-nT standard deviation is used. As can be observed, this 2.5 amplitude reduction factor (from B = 0.5 nT in the nonfiltered case) very well describes the level of FAC estimation errors when filtered data are used (a more precise analysis would take into account the dependence of FAC estimation error on frequency). Therefore, in all cases analyzed in the article, the limits of the EZ due to instrumental noise were obtained based on effective magnetic noise of B = 0.2 nT and a threshold in FAC estimation error of 0.1 μA/m 2 .
Using simulated data affected by noise, we evaluate the errors in FAC density estimated through BI method applied to filtered magnetic field data, by computing the standard deviation of BI solution with respect to the simulated current. Panel (i) shows the ratio between what we obtained and the similar standard deviation when LS method is used. We found that, at the EZ limits defined above, the BI method provides estimation errors only slightly higher, by ∼10% (equivalently, the EZ extension in the BI method is only two data points larger). This is a consequence of the almost parallel orbits for the lower satellites ( ∼ 0.012 in this case); the 10.1029/2018JA026439 situation changes if data from swB are combined with swA or swC data (see section 3.5) or when a higher FAC error level is accepted.
In Figure 3 the EZ limits due to instrumental random noise are indicated by vertical lines both for unfiltered (panels c, e, and f) and filtered (panels h and i) dB data. When data are not filtered, the EZ extension is around 7.5 times the standard (i.e., 5 s) along-track separation distance, whereas when filtered data are used, it is around 3 times that distance.
Considering now the presence of offsets in the data, equation (B3) and simulations show that the corresponding EZ extension remains basically the same (around 7 times the standard, 5-s, along-track separation distance), when quad configurations with smaller along-track sizes are used. This type of errors cannot be removed by data prefiltering. However, their impact on FAC estimation decreases as the mission progresses due to the constant data calibration effort. The aspect is easily observed when different baselines of L1b data files are used to analyze EZ intervals where no FAC are expected: for more recent baselines the region with (unreal) current density is smaller (not shown).
For Swarm, the instrumental noise is expected to be around 0.1 ÷ 0.2 nT rms, up to ∼0.5 nT rms for swC after the AMS failure on that satellite . In the present analysis a conservative level of 0.5 nT for all satellites has been chosen. Note, however, that the FAC estimation rests on the assumption of linear magnetic field variations; hence, any random deviation from linearity contributes to the FAC error. In addition to instrumental noise there are also small-scale magnetic fluctuations (e.g., Alfvén waves) that do not show a correlation between the different sensors. The size of these fluctuations depend on the level of smoothing and can easily increase the effective random error to several nanoteslas (unfiltered data). To provide values for the instrument offsets, the error estimates on magnetic field vector, that is, B error , available from the L1b data files, have been considered. For example, in the case of swA, the most recent baseline, that is, 0505, specifies B error values close to ∼0.25 nT in all directions during the events presented in Figure 1 (a decrease from ∼0.75 nT, specified in baseline 0408). To evaluate j off , the most detrimental situation, that is, opposite signs for the along-track components of B off at swA and swC (see equation (B3)), has been chosen.
Using the above criteria to define valid FAC estimates, the dual-s/c LS method allows an EZ systematically smaller than the 4 • spherical sector used in the official L2 product. In order to see how much of the auroral data is gained, an automatic procedure to identify the auroral region location has been used. The procedure has been applied on a number of 25 days, randomly distributed over 4 months of Swarm data. More than 1,500 auroral oval crossings have been identified, with roughly 9.6% of data points falling within the EZ as defined in the L2 product (i.e., missing values in the ESA database). By using the dual-s/c LS method, only 1.6% of the total number of points within the auroral intervals had to be excluded due to instrumental random noise. If one adopts a stricter criterion and require a 2 confidence level (∼95.5% chances of realization) for the 0.1 μA/m 2 threshold in FAC estimation error, the percentage of data points from auroral intervals that had to be excluded is ∼3.4%, roughly the same as when pure instrumental offsets are considered (i.e., ∼3.8%).

Dual-s/c LS Method Applied on Nonrectangular Configurations
In the above discussions, the LS dual-s/c method has been applied to the standard four-point satellite configuration, that is, the configuration described in Ritter et al. (2013), to facilitate comparison with the ESA L2 FAC product. This section illustrates the additional capability of the Vogt et al. (2013) method to work on nonrectangular, skewed configuration.
In the method one can change two parameters that control the quad geometry, that is, the orbital time lag between the satellites and the along-track separation distance (parameters and, respectively, M in the left part of Figure A1). One reason to change the along-track separation distance is to investigate the FAC structures at higher spatial resolution and therefore to put in evidence sharper changes in current density. For example, using an along-track separation of 1-s orbital distance, the FAC estimation errors stay within the tolerable range even at lower latitude, that is, ∼40 • , when the quad is stretched along the east-west direction. As for the higher latitudes, the EZ extension basically does not change, according to the discussion in section 3.3. Similarly, by adjusting the orbital time lag between the satellites one can fine-tune the satellite configuration to better characterize current sheets inclined with respect to spacecraft orbit or one can obtain FAC estimations that rely less on the condition of time stationary current structures (see below).

Figure 5.
Dependence of dual-s/c estimation on current sheet inclination: synthetic data analysis. Left: Trajectories and configuration of two Swarm-like satellites crossing an inclined current sheet. Two quad configurations were used, that is, with no orbital lag between the satellites (blue) and with an orbital lag tuned to the current sheet inclination (magenta). Right: (a and b) magnetic field perturbation recorded by the each satellite. (c) The simulated (green) and predicted (blue and magenta) current densities. (d) The difference between predicted and simulated currents. (e) Field-aligned current density estimation error due to instrumental noise.

Performance of Dual-s/c Method on Inclined Current Sheet
To study how the current sheet inclination impacts on the predictions of dual-s/c method, an analysis based on synthetic data has been performed. Figure 5 presents a case of two Swarm-like satellites moving on parallel orbits that cross a current sheet at 45 • inclination. The parameters of the simulation were chosen to be close to the Swarm context, using a cross-track distance of 48 km and an orbital time lag of 10 s between the satellites. The associated magnetic field perturbation was simulated by a superposition of two tangential hyperbolic profiles that vary only along the sheet's normal direction (see panels a and b on the right). The corresponding current density profile is shown in green on panel (c).
For a quad configuration with the two satellites aligned in orbital phases, that is, the configuration used to generate the L2 FAC product, the dual-s/c method predicts a peak current density around 10% less than the simulated current (see the blue lines in panels c and d), although the configuration was close to a square (see the quad shown in blue on the left of the figure). If one uses a time shift between the satellite positions such that one side of the quad is aligned parallel to the current sheet (configuration shown in magenta), the mismatch between the predicted and simulated current densities decreases to around 2 % (see the magenta lines in c and d). Note that the FAC estimation errors have an accepted level for both configuration, as shown in panel (e) (a flat value of 0.5 nT for B was considered in the computation).
The explanation rests on the fact that a satellite configuration aligned with respect to current sheet samples the structure at a smaller scale along the normal direction, that is, the direction of the magnetic gradient. As a consequence, the linear gradient estimation scheme does a better job and the smoothing effect encounter at the current density peaks is smaller. One can foresee (and the synthetic analysis confirms) that for a perfect square configuration, the most detrimental current sheet inclination is at 45 • . Also, for a fixed along-track separation distance and satellites aligned in orbital phases, that is, as in the case of L2 product, the current sheet estimation for an inclined sheet is less precise at lower latitudes since there the cross-track separation is bigger.

Dual-s/c Method on Configurations That Use Simultaneous Measurements
One quad configuration of potential interest in the dual-s/c method combines pairs of simultaneous measurements, that is, satellite positions that are not compensated for their relative orbital phase. The current density estimations based on such configurations rely less on the time stationarity in dB; for example, if one uses an along-track separation of 1 s, the single-s/c and dual-s/c solutions depend to the same degree on this assumption. This would allow, in principle, to use the dual-s/c method in studies of FAC structures close to the small-scale regime, which exhibit a more dynamic behavior than the current structures in the mesoscale range (Gjerloev et al., 2011;. Certainly, when time nonstationary can be neglected, using a rectangular configurations (or one aligned with the current sheet) with 1-s along-track separation would be more appropriate.
Before illustrating the results on a Swarm event, we discuss the errors in current estimation induced by such configuration. There are two issues related to that: the characteristic scale length along the current sheet normal (roughly along satellite orbit) may increase as compared to the case of standard, 5-s, along-track separation rectangular quad, considering that the relative orbital phase between swA and swC varies between 4 and 10 s during the mission. Similarly, since in most of the cases the FAC sheet is approximately oriented perpendicular to the spacecraft orbit, and relevant quantities are roughly varying along the normal direction, such a configuration is not ideal; indeed, when the along-track gradient of dB east is calculated, dB east is not measured by swA and swC at the same coordinate along the normal axes, as in the standard quad.
To see how important these effects are, we refer again to the analysis on synthetic data presented in section 3.2. We changed the swC orbit such that the same 2-D upward-downward current structure is sampled by an array of two sensors moving on parallel orbits 50 km apart, that is, about the Swarm cross-track separation at 72.5 • latitude. On the right part of Figure 3, the new swC orbit is indicated by red dashed line. Three configurations, illustrated as well on the right in Figure 3, have been considered to estimate the FAC density: a rectangular, 5-s along-track separation quad used in the standard method of Ritter et al. (2013), and two skewed, 1-s along-track separation quads that simulate configurations based on simultaneous satellite positions, with an relative orbital lag of 4 s and, respectively, 7 s.
For each configuration, the difference between the estimated and simulated current at the quad's mesocenter has been calculated; the results are shown in Figure 3j. In case of the rectangular quad (blue line) the difference with respect to simulated current comes from the averaging effect introduced by the scale of the sensor array. The second configuration, corresponding to the red line, has the same along-and cross-track extensions as in the rectangular quad but swA and swC were not aligned with respect to the main direction of variation in dB, which is roughly along y axes in the region of relevance. The difference between red and blue lines, that is, only about few tens of nA/m 2 , quantifies the effect of using a skewed quad instead of a rectangular one. For the third configuration, which has a scale along y axes of 8-s orbital distance, the difference with respect to simulated current values, shown with green line, is higher but still remains below 0.15 μA/m 2 . It is well understood that the above numbers depend on the scale of variation in the current structure; for small-scale FACs, when the assumption of linear field variations starts to be problematic, the discrepancies between estimated and actual current are higher.
In the above simulations, no fluctuations in dB have been considered, for a clear illustration of the analyzed effect. However, considering a level of magnetic noise of B = 0.5 nT and no data filtering, the uncertainty in current density estimation, evaluated with equation (B1), assumes relative small values even for the configurations that use simultaneous satellite positions, that is, 60 and, respectively, 75 nA/m 2 , to be compared with less than 15 nA/m 2 in the case of rectangular quad configuration. Simulations show that data biases of ∼0.25 nT has an even smaller effect on FAC density estimation. One can therefore conclude that the Vogt et al. (2013) method estimates the FAC density pretty accurately if one uses configurations based on simultaneous satellite positions, at least for the typical values of orbital time lag. The error in FAC density estimations remains within a tolerable range and, as an important advantage, the results are less dependent on the condition of time stationarity current sheet structure. Figure 6 presents one Swarm event analyzed with configurations that use simultaneous satellite positions. For this event, the orbital lag between swA and swC was 5 s and the MVA indicated an inclined current sheet, about 35 • with respect to the spacecraft velocity. This fortunate inclination leads to quads roughly aligned with the FAC structure, ideal for probing the filamentary current structures. At the bottom of the plot, the Swarm configuration at three moments are indicated as projections on the NE plane of the local NEC frame, together with the direction along the current sheet (red arrow).
In panel (e), a comparison between two FAC density estimations is presented, one obtained in the standard way, that is, 5-s rectangular configuration applied on filtered dB data (magenta), and one obtained by using The magnetic field perturbation, dB, at swA and swC, respectively. (c) For both satellites, the magnetic field residues (differences between the unfiltered and filtered dB) along satellite velocity direction. (d) The magnetic field residues along east direction. Panel (e) compares the standard dual-s/c solution, that is, the method applied on filtered dB data using a 5-s rectangular quad configuration (magenta), with the solution obtained by using the configuration based on simultaneous satellite positions, separated 1 s in the along-track direction, with no data filtering (blue). (e) The level of FAC estimation error due to instrumental noise for both configurations. The quad that combines pairs of simultaneous satellite position is roughly aligned with the FAC sheet for this event (see the configurations plotted at the bottom). FAC = field-aligned current.
the configuration based on simultaneous satellite positions, separated 1 s in the along-track direction, with no data filtering (blue). In the latter, the fine structure of FAC density is revealed. However, by looking at the small scale magnetic residues observed at swA and swC, that is, the differences between the unfiltered and filtered dB along the spacecraft velocity (panel c) and along the east direction (panel d), it seems that some of these fluctuations are either temporal variations, indicative of some wave activity, or are due to spatial structures with a coherence length smaller that the cross-track distance. This aspect, in line with previous observations (see, e.g., Forsyth et al., 2017;, adversely influences the outcome of the dual-s/c estimation. Therefore, even when using configurations based on simultaneous satellite positions, more appropriate to investigate FAC structures affected by time nonstationarity, the application of dual-s/c method on unfiltered 1-Hz magnetic data should be exercised with care; depending on individual event, a (slightly) filtering of the data might be needed.

Results Involving swB
At least at the beginning of the mission, the orbital planes of swB and of the lower Swarm satellites were close enough one to each other to expect that, occasionally, the standard swA-swC-based FAC density estimate can be complemented with additional estimates based on data provided by the swB-swA and/or swB-swC pair. Figure 7 presents two such events.
For each event, panels (a)-(c) present the magnetic field perturbation in NEC, with swA and swB data shifted in time so that all satellites are aligned in orbital phase (the time shifts are indicated at the top of each plot). Since when using swB in dual-s/c method, the orientation of the quad can be an issue, panel (d) shows, for all two-satellite combination, the angle between the local magnetic field and the quad normal. In the algorithm, a minimum value of 25 • between the local magnetic field and the quad plane has been required for a reliable reconstruction of FAC. In case of the first event (left) this condition is obeyed only for part of the interval when the swB-swA or swB-swC data are combined. Panel (e) compares all three possible dual-s/c FAC estimates based on Swarm data, while in panel (f) the corresponding estimation errors due to instrumental noise are shown (since the satellite tracks are more inclined and at different altitude, the errors caused by instrument offsets are less important, being restricted to smaller regions). On the bottom, the four-point satellite configurations are shown at three instances (beginning, middle and end of the interval), projected on the NE plane of the local NEC coordinate frame. The mesocenters of different quads are shown in the same color as the corresponding FAC estimate.
The first auroral crossing (left) looks stationary, as can be judged from the evolution in dB, and consequently, the different FAC estimates are very similar in the common range of validity. The small misalignment between the current densities can be explained by the differences in the location of the mesocenter. Also, the smaller amplitude predicted by swB-swA and swB-swC pairs at the end of the interval (after ∼ 16:08:45) is linked with the larger extension of the corresponding quads along the north direction, that is, along the magnetic field gradient direction. The fact that swB is orbiting at higher altitude, where the magnetic perturbations are smaller, has only a small influence on the FAC estimation. Indeed, by mapping the swB magnetic field at swA/swC altitude, we arrived at a correction in current density of only ∼1%. This event is an example when swB can be used in combination with one lower satellite to resolve a FAC structure located close to the pole. Indeed, note that the standard swA-swC based FAC estimate is not provided for ∼40 s around 16:08:00, that is, within the corresponding EZ, while the other two estimates show similar evolution, consistent with the changes in dB.
The second event illustrates how using swB data in the dual-s/c method can provide additional information on the FAC structure. Here the current sheet is basically aligned along east-west direction in NEC (see the components of dB), roughly aligned with the sides of the quads as well. Since no time shift between the swA and swB data has been introduced, the differences between swA-swC (green) and swB-swC (blue) FAC estimates can be attributed to the spatial gradient in the current density along the sheet. For example, around 19:46:14 and 19:46:54 (indicated by the vertical dashed lines) the difference in FAC amplitudes are roughly 1 and 0.8 μA/m 2 , respectively, while the distance between the position of the mesocenter along E axis in NEC is around 120 and 140 km, respectively.
For the event presented on the right of Figure 7, the parameter (nonparallelism in the orbital tracks; see Figure A1) of the swB-swC configuration is only ∼ 0.065, and therefore, the FD, BI, and LS methods are expected to provide very similar estimates. However, later into the mission, when starts to take higher values, using LS method is more appropriate since it involves smaller FAC estimation errors. In that respect, see Vogt et al. (2013) where the normalized error amplification factors of the FD and the LS methods are compared for = 0.5 (see Figure 4 there).

Swarm FAC Density Estimation With Three-s/c Method
When all Swarm satellites form a close configuration, the FAC density can be estimated by applying the three-s/c method presented in section 2. The different orbital velocities of the upper and lower satellites makes possible a spacecraft alignment in latitude for roughly four to five auroral crossing events every other 6 days, depending on how strictly we require the satellite to line up (more on this below). However, taking into account the relative drift of swB orbital plane with respect to the orbital planes of the lower satellites, a spatial satellite formation appropriate for three-s/c FAC estimation is more likely to occur at the beginning of the mission. Later on, the separation in longitude between swB and swC/swA might be too large when compared with the longitudinal scale of the FAC sheets, as suggested by the analysis of .
In section 3.1 the three-s/c method has already been applied on virtual positions of swA and swC, the conclusion being that it provides consistent results with the dual-s/c method. Though, what makes the three-s/c method distinct is the possibility to use simultaneous measurements and therefore to drop the assumption of time stationarity in dB. Indeed, this assumption is not needed to infer the current density along the direction normal to the spacecraft plane, unless one wants to combine satellites virtual positions, for example by shifting in time measurements provided by one or two satellites in order to reach a more favorable spacecraft configuration. Therefore, in principle, high-resolution magnetic data, acquired simultaneously by all Swarm satellites, can be used to obtain instant estimations of FAC density. However, this does not imply the method can resolve spatial scale structures smaller than the constellation characteristic lengths. For example, if Swarm encounters a perfectly stationary planar current sheet, what matters is how the scale of current variation compares with the constellation extension along the sheet normal.

Applicability of the Three-s/c Method
Several aspects have to be taken into account when the three-s/c method is applied in the Swarm context. One refers to the orientation of the plane formed by the observation points, which depends on the orbital parameters. Unlike the case of dual-s/c technique, here one has to check whether or not the local ambient magnetic field vector is sufficiently out of the spacecraft plane. For that purpose, in the present analysis a threshold of 25 • has been considered to allow for a reliable reconstruction of FAC from the current density inferred along the normal direction.
Then, right from the start of the operational phase, the characteristic scale lengths of the three-s/c close configurations are larger than in the dual-s/c method and growing as the mission progresses. As a consequence, the aspects of FAC sheets extension in longitude as well as that of deviation from the (assumed) linear field variations are of main concern when evaluating the method validity range. Note that the criterion based solely on the level of FAC estimation errors (the approach used in the dual-s/c method to define its validity range) could be misleading since, according to equation (B2) from B, the errors assume small values for large interspacecraft separation distances even if, for example, the satellites are not sampling the same current structures or the magnetic field is not varying uniformly between the points of observation.
The issues discussed above make clear that for applying the three-s/c method all the satellites should be quasi-simultaneously within the auroral oval and, in addition, some indication should exists that the measurements refer to the same current structure. The latter condition is usually checked by performing a cross-correlation analysis (see, e.g., in Forsyth et al., 2017 a more involved technique has been used that combines cross correlation and linear fitting). To evaluate the occurrence frequently of such circumstances, the first 6 month of data after the mission commissioning, that is, May-October 2014, have been considered. Using an automatic procedure to identify the auroral oval and by requiring small time delays between the corresponding dB profiles, that is, time lag between the leading and trailing satellite less then 30 s, a number of roughly 130 events where individual magnetic field perturbations are reasonably cross-correlated have been identified. The above condition on satellite orbital lag could be relaxed if time shifted, that is, nonsimultaneous, measurements are considered for analysis, relying thus, to some degree, on the assumption of time-stationary FAC current structures.
Another useful parameter to be checked when applying the three-s/c method involves the stability of the solution. Similarly to the case of dual-s/c method (see Appendix A), finding an estimation for the current density along the direction normal to the satellite plane boils down to the inversion of a 3 × 3 matrix with one eigenvalue close to zero. Assuming int and max to be the remaining eigenvalues, the quality of the inversion (i.e., the condition number) is provided by CN = max ∕ int . A high value for CN reflects a stretched configuration, with high impact on the current density error level. Figure 8 presents two auroral crossings by Swarm where the three-s/c method has been applied. In panels (a)-(c), the NEC components of the magnetic field perturbation, dB, indicate the Swarm satellites crossing together the FAC structure. Panel (d) and (e) present the logarithm of the condition number, CN, and, respectively, the angle between the normal to the spacecraft plane and the direction of local ambient magnetic field. Panel (f) shows a comparison between the FAC density estimated by applying the dual-s/c method on filtered dB data (black) and the three-s/c method based on unfiltered magnetic data (magenta). The corresponding FAC estimation errors due to instrumental noise are shown in panel (g). On the bottom, the spacecraft configurations are shown at three instances, that is, beginning, middle, and end of the interval, projected on the NE plane of the local NEC coordinate frame. In the three-s/c method one assumes uniform magnetic field variation between the instantaneous satellite positions (region inside the triangle shown in dashed lines), while the dual-s/c method assumes linear field variation within the area formed by the virtual positions of swA and swC (quad configuration, shown in continuous lines).
For auroral crossings with quasi-identical magnetic perturbation recorded by all sensors, one expects (and the observation confirms) no significant difference between the dual-s/c and three-s/c solutions. The events chosen here partially display differences in dB data recorded by individual satellites, either due to nonstationarity or spatial variations, to illustrate the information added by the three-s/c method. For the 17 June 2014 auroral crossing, a detailed comparison between the dual-s/c and three-s/c solutions is provided in section 4.2.
Simultaneous measurements have been used in the three-s/c method, the fine structure of FAC density reflecting thus changes in the instantaneous (∇×dB) n values, estimated based on unfiltered 1-Hz data. Similarly to the case of dual-s/c method when configurations based on simultaneous measurements are used (see section 3.4.2), the small-scale magnetic field perturbations recorded by each satellite affects the three-s/c solution, although, due to the larger separation of the sensors, less influence on the results is expected. Therefore, depending on each individual event, some (slightly) filtering of the data might be needed. In addition, due to the larger length scale involved, the three-s/c method cannot resolve spatial scales smaller that the dual-s/c method. Note that the two FAC density estimations refer to different points (different latitude and longitude of the corresponding mesocenters), shown with magenta (three-s/c) and black (dual-s/c) crosses at the bottom of Figure 8. Therefore, to facilitate result comparison, the dual-s/c solution has been plotted with a constant time delay that approximately takes into account that effect (note that the distance between the mesocenters is slightly changing in time).
We tested the stability of the three-s/c method with respect to the magnetic field model used to derive dB. Considering the scalar potential associated with the magnetic model, it turned out that the estimation of ∇ × ∇ may assume nonvanishing values, since the approximation of linear field variation is not holding at the (larger) scale of three-s/c formation. Therefore working directly with B, instead of dB, data can induced artificial offsets in the FAC estimation. Nevertheless, when magnetic field perturbations based on several recent magnetic models were used in the analysis, we found only negligible differences between FAC density estimations, up to a few tens of nA/m 2 .

Comparison With Other Methods
As an illustration of how the different FAC estimation methods complement each other and provide consistent results, we present a detailed analysis on the 17 June 2014 auroral crossing shown in section 4.1 (Figure 8, on the right). For this event, the dual-s/c and three-s/c FAC estimations are close one to each other most of the time, except for an interval of ∼20 s centered on 03:54:55, where the difference in current density reached up to ∼2.5 μA/m 2 . This interval will be investigated further below. The MVA indicates a current structure almost perpendicular to the satellite trajectories, that is, only ∼14 • deviation of FAC normal direction with respect to velocity direction of each satellite. The current sheet planarity, measured as MVA eigenvalue ratio, max ∕ min , is well established for swB (∼ 14) and swC (∼ 16), and less so for swA (∼ 7). Figure 9 presents a comparison of the results obtained when the single-s/c, dual-s/c, and three-s/c methods are applied. In order to roughly align the locations where different methods estimate the current density, the dual-s/c and single-s/c solutions have been properly shifted in time with respect to three-s/c solution. To help data interpretation, panel (a) shows the difference between the (filtered) magnetic field perturbation seen at swA and swC in the VEC reference frame, introduced in connection to Figure 2. Panel (b) presents a comparison between three-s/c and different dual-s/c current density estimations based on unfiltered (three-s/c) and, respectively, filtered (dual-s/c) dB data. The color code used to plot FAC density estimations that rely on various satellite configurations has also been adopted in panels (c) and (d). The FAC density estimations based on three-s/c configuration (magenta), and on dual-s/c configurations that involve swA-swC (black), swB-swC (blue), and swB-swA (green). Panel (c) shows, for each satellite configuration, the contributions to BI in the direction aligned to the current sheet. (d) Similar to (c) but refers to the direction perpendicular to the current sheet. Panel (e) compares swA (blue) and swC (red) single-s/c FAC density estimates with the BI contribution along the current sheet in case of swA-swC configuration. At the bottom, the three-s/c and swA-swC configurations are shown projected on the VE plane of the local VEC coordinate frame. The orientation of the current sheet is indicated in the middle inset. See text for details. BI = boundary integral; FAC = field-aligned current.
The quantities plotted in panels (c) and (d) are based on BI concepts. As discussed in section 3.1, in the dual-s/c case the LS and BI methods provide basically the same results outside the EZ region. Similarly, while the (Vogt et al., 2009) method is a complex three-s/c technique with a broad range of applications (e.g., planar or full spatial gradient determination, discontinuity analysis, and wave identification), when it comes to the estimation of field rotation along the direction normal to the array plane, one can also apply Journal of Geophysical Research: Space Physics 10.1029/2018JA026439 the BI method to obtain the same results. The discretized form of (∇ × dB) n in the BI approach is where (i,j) must be a cyclic permutation of 1,2,3 (three-s/c case) or 1,2,3,4 (four-point configuration in the dual-s/c case), r ij = r j − r i is the difference in the position vectors of satellite i and j, and A the area of the triangle/quad formed by the satellite array.
For the analyzed event, all the methods are estimating (∇ × dB) basically along the same direction since, for example, the angle between the plane formed by the instantaneous position of all three satellites and the plane of the swA-swC quad configuration is only around 12 • . If e x and e y are two unit vectors in that plane, one can write r ij = e x (e x ·r ij )+e y (e y ·r ij ) and the sum in equation (1) can be separated in two components, that each quantify the contribution to the BI along e x and e y , respectively. In Figure 9 we have chosen these unit vectors to be the eigenvectors e min and e max associated with the directions of minimum and, respectively, maximum magnetic field variance provided by MVA. Therefore, panel (c) presents contributions from the three-s/c and dual-s/c BI methods associated with the FAC normal direction, which for this event is roughly along v in VEC frame. Likewise, in panel (d) the similar contributions to BI associated with the direction along the current structure are presented.
As can be seen from Figure  The dual-s/c FAC estimations based on swB-swA and swB-swC configurations (green and, respectively, blue lines in panels b-d) provide evidence for a localized (i.e., relatively small scale) spatial variation which primarily affects swA measurements along the current sheet normal direction, as the origin of the difference between swA-swC based dual-s/c and three-s/c FAC estimations. Indeed, if one considers how the satellites are located along the e max direction, that is, swA, followed relatively close by swC and, much far away, swB (see the spacecraft configurations plotted at the bottom of Figure 9), one can observe that for the swB-swC configuration the BI contribution along the current sheet normal is small (panel c, blue line), while in the case of swB-swA, the same quantity (green line) takes slightly higher values but still much smaller then in the case of swA-swC configuration, as if the swA influence is leveled out. As expected, the three-s/c method is also smoothing the local variation encountered by swA. The interpretation of a localized variation which affects swA data is consistent with the result of MVA, where the smaller value for the eigenvalue ratio indicates the sampling of a current structure with a relative higher degree of nonplanarity. It is also worth mentioning that the spatial, instead of temporal, nature of this variation is supported by the dual-s/c FAC estimation that uses simultaneous satellite positions (not shown here but with a similar evolution), which is less dependent on the time stationary assumptions.
In panel (e) the single-s/c FAC density estimates based on filtered swA (blue) and swC (red) data are shown. These solutions are capturing only the along-track variation in the East component of dB. For the swA-swC configuration one expects a contribution to BI along the same direction (approximately along e max for this event) to be an average of the two single-s/c solutions. This is indeed observed, see the black line roughly in between the red and blue lines.

Three-s/c Method Performance on Synthetic Data
The three-s/c FAC estimation method has also been tested on synthetic data to investigate its performance in various current configuration. In Figure 10, similar to Figure 3, the results of one such test are shown: An array of three Swarm-like sensors is sampling a 2-D upward-downward cylindrical current structure of different scale lengths. The satellite orbits are shown on the right, together with the array configuration at three instances, marked by blue dots along the swA trajectory and by blue vertical lines on the left part of the figure. The simulation parameters were close to the Swam context.
Panel (a) presents the three components of the vector magnetic field perturbation recorded by swA, that is, dB swA , while panels (b) and (c) show the difference dB swA − dB swB and dB swA − dB swC , respectively. In panel (d), the comparison between the three-s/c FAC estimation, that is, j es (red), and the simulated current at the mesocenter, that is, j m (blue), shows a very good agreement, apart from the region between the vertical black lines, that is, within the three-s/c EZ, where the three-point configuration is very stretched. For the particular situation presented in Figure 10, this occurs close to the point when the swA and swC orbits are crossing but in general, it occurs whenever the satellites positions are quasi-aligned.
Panel (e) shows (in black line) the difference between j es and j m , together with the level of errors in FAC estimation due to instrumental noise, ±j B (red lines), evaluated through equation (B1) applied to the three-s/c configuration. Within the EZ, the three-s/c estimation j es is unstable, resulting in a dramatic change in the evolution of the black line from one simulation run to the other. In panel (f), the standard deviation of j es − j m (black), obtained from 500 simulation runs, is compared with the j B (red), while panel (g) shows the ratio es − m ∕ B . The latter fluctuates around 1, which fully supports the error estimation formula from equation (B1). Panel (h) presents the average values of j es −j m ; this quantity is vanishingly small provided that the scale of the satellite configuration is small enough compared to the characteristic length of the magnetic field variation. To the end of the interval, that is, around the third moment for which the satellite constellation is shown, the condition is not entirely fulfilled, meaning that the assumption of linear magnetic field variation is broken, causing a bias in the three-s/c FAC estimation of around 0.1 μA/m 2 . The shape of the configuration cannot be the cause of this systematic discrepancy since the logarithm of the three-s/c condition number CN, presented in the last panel, takes reasonably values there. Simulations with noise-free data when data offsets are present (not shown) indicate a corresponding EZ restricted to smaller region when B off components assume values around ±0.25 nT.
We used synthetic data to test how ionospheric electrojet influences the three-s/c FAC estimations. For that purpose, to the setup described in the dual-s/c case, we add an additional sensor that mimics swB. The sensor flies on a parallel orbit, 200 km separated in the cross-track direction and at different altitude, that is, 50 km higher (see Swarm orbit and constellation on the right side of Figure 4.) For the same reasons mentioned in section 3.2 (Ampère's law applied on a close path, through which no current is flowing), the predicted current density is expected to be small, indicating just the degree to which the assumption of linear variation for dB is violated. This is actually the case: Panel (f) on the left part of Figure 4 shows that the three-s/c FAC density estimations are ∼40 nA/m 2 for the localized electrojet (red line) and ∼20 nA/m 2 for an electrojet spread in latitude (blue line). Similarly to the dual-s/c case, by using nonparallel satellite tracks in the simulation such that the constellation becomes singular above the electrojet, we found out that the three-s/c EZ limits are not affected by this current system. Within the limits imposed by the underlying assumptions, the three-s/c method effectiveness to infer the FAC density for complex structures is confirmed by tests on synthetic data. Similarly, the level of FAC estimation errors due to imprecision in magnetic field measurement is also supported by the analysis.

Discussions and Conclusions
The L2 dual-s/c algorithm, designed and implemented before the launch of Swarm, had to deal in a simple manner with the task of FAC density estimation. The statistical location (and orientation) of the auroral oval, imposed processing parameters that ensure a configuration close to a square (the most favorable shape) in that region. Inevitably, the algorithm could not take into account all the aspects of the problem and has not benefited from the expertise accumulated while working with real measurements. The results presented in this technical report substantiate the application of LS dual-s/c and three-s/c FAC density estimation methods on Swarm data, with the goal to enhance the exploration of high-latitude current structures, beyond the official FAC product. The influence of various factors that characterize the current structures (orientation, scale of variation, spatial extent) or the geometry of the spacecraft constellation have been investigated, in order to document the capability of the methods as well as their limitations.
In section 3, the Vogt et al. (2013) dual-s/c method has been thoroughly tested on real and simulated data. Several clear advantages of the method over the standard L2 product, which come either from its LS approach or from the capability to work on flexible geometry, have been established, as follows: -In comparison with other dual-s/c implementations, the LS approach to this overdetermined problem provides more stable FAC estimates close to the geographic poles, where the quad shape become extremely elongated (see Figures 1 and 2). As a consequence, a systematically smaller EZ is obtained. From our evaluation, for a tolerable error in FAC density estimates of 0.1 μA/m 2 , the percentage of FAC data points that have to be excluded due to instrumental noise is ∼1.6% (1 confidence level in FAC density estimation) or ∼3.4% (2 confidence level). When data biases of around 0.25 nT are present, the percentage slightly increases to ∼3.8%. For comparison, the 4 • polar cap sector used in the official Swarm L2 product implies that 9.6% FAC data points are excluded. -The region where the dual-s/c method does not provide reliable results can be naturally defined based on the level of FAC estimation errors, assessed through a corresponding error evaluation scheme. Studies on real and simulated data proved that this criteria is meaningful. -In the method we are using, the results are basically the same irrespective of the magnetic field model used to compute dB, or even when the technique is directly applied on the magnetic measurements, without subtracting a magnetic field model. -Contrary to the case of L2 product, based on tight geometry (i.e., rectangular quads with fixed 5-s along-track separation), in the Vogt et al. (2013) method no simplification on the geometry of the four-point configuration is required. This allows to investigate FAC structures at higher spatial and temporal resolution or to fine-tune the satellite configuration to better characterize inclined current sheets.
In section 3.1 a comprehensive comparison between the dual-s/c LS solution and the solutions provided by other methods was presented. The solution obtained from an in-house implementation of the BI method is basically identical, apart from the region close to the geographic poles, where the LS method better weights the contributions from the four data points, providing thus more stable results. Differences with respect to the single-s/c solutions arise since the latter are capturing only the along-track variation in the east component of dB. At the same time, the three-s/c solutions based on dB recorded at three apexes of the four-point quad configuration are very close to the dual-s/c solution outside the EZ, in line with the observations from 10.1029/2018JA026439 Dunlop et al. (2015). Figure 2 clearly shows that the three-s/c solutions are properly taking into account both the along-and cross-track variations in dB.
Extensive tests on simulated data confirm the methods effectiveness to infer the current density of complex FAC structures. The tests fully support the error analysis scheme provided in the method and were essential in arriving at a meaningful definition of EZ. Other specific analyses investigated (i) the FAC estimation errors in the in-house implemented BI method. Due to the almost parallel tracks of swA and swC, the adopted definition of EZ implies an extension of that region only slightly larger than in the LS method; (ii) the degree to which the electrojet is influencing the FAC estimation: it has been shown that this effect is negligible, the needed FAC density adjustment being of the order of nA/m 2 , that is, 3 orders of magnitude smaller than the typical FAC values; (iii) the effect of current sheet inclination and motion. At optimum latitude, that is, where the quad along-and cross-track distances are roughly the same, the simulations showed an error in FAC estimations up to 10% for a current sheet inclined at 45 • . For the type of quad configurations used in the L2 product (fixed along-track separation and satellites aligned in orbital phase) the FAC estimation of an inclined sheet is less precise at lower latitudes, where the cross-track separation is bigger. The errors due to current sheet motion are less important, that is, below a few percent even for the most dynamical events reported in the literature.
The advantage offered by the dual-s/c method developed in Vogt et al. (2013) to work on a more general spacecraft configuration, opens the possibility to combine simultaneous satellite positions, that is, not compensated for their relative orbital phase, relaxing thus the assumption of time stationarity FAC structure. Simulations show that, for 1-s along-track separation and typical orbital lag between swA and swC, the errors are only slightly higher than in the case of standard, rectangular configuration, provided that the scale of current structure is not too small. On the other hand, when such configuration is used to analyzed unfiltered Swarm data, the estimation might be affected by uncorrelated, local fluctuations measured by each satellite.
Another possibility, derived from the capability to work on flexible geometry configuration, is to combine swB data with data from a lower Swarm satellite, when the former is close enough and roughly aligned in latitude. The LS approach is more appropriate to be used for such nonrectangular/nonparallel configuration, since it involves smaller FAC estimation errors than the implementation based on BI or FD concepts.
In the first event presented in section 3.5, the dual-s/c solutions that use swB measurements are able to resolve FAC structures located close to the pole, that is, within the EZ of the swA-swC solution. The second event illustrates how using swB data in the dual-s/c method can provide additional information on the FAC structure, that is, to estimate the spatial gradient in the current density along the sheet.
As discussed in section 2, the methods used to infer the FAC density rely on certain underlying assumptions. The linear variation in dB is assumed by all methods but starts to be critical once the scale of the sensor array is large, like in the three-s/c method or dual-s/c method that involves swB data. The time stationarity is required in the dual-s/c case and in the three-s/c method if virtual satellite positions are combined, while the single-s/c method heavily relies on the current sheet planarity. In addition, knowing the orientation of the current structure helps us to correct the FAC estimations (single-s/c) or to fine tune the array geometry (dual-s/c and three-s/c) while the information on current sheet motion allows a small correction on FAC estimation. It is therefore important to check these aspects through available techniques like MVA, correlation analysis, and timing analysis, before the application of FAC estimation methods.
The application of three-s/c method to Swarm data (section 4) shows that this technique brings additional valuable information, associated to a different (larger) scale. In the Swarm context, critical aspects refer to the longitudinal extension of the FAC structures, addressable, to some extent, by a cross-correlation analysis of the individual magnetic field perturbations, and the assumption of linear field variation between the points of observations, that is, at a significantly greater distances than in the standard dual-s/c method. When instantaneous observations are used, time stationarity assumption is not needed and data with higher temporal resolution can be used as input. The detailed analysis on the 17 June 2014 auroral crossing (see section 4.2) presents how the different FAC estimation methods complement each other and provide consistent results. The analysis on synthetic data illustrates the performance of three-s/c FAC estimation method and allows to study the aspects of nonlinear field variation (i.e. constellation scale versus gradient characteristic scale). It also fully supports the results provided by the error evaluation scheme (equation (B1) tailored to the case of three-s/c configuration) and proves the negligible effect of the electrojet on FAC estimation (i.e., around 2 orders of magnitude smaller than the typical FAC values).
By analyzing three-satellite conjunction in the auroral region, a number of roughly 130 events with reasonably correlated magnetic field perturbations at all satellites have been identified between May-October 2014 (angular separation of swB orbital plane between approximately 5 • and 11 • ). When time-shifted, that is, nonsimultaneous measurements are considered, relying thus, to some degree, on the assumption of stationary current structures, more events could be analyzed. For large angular separation of swB orbital plane, the application of three-s/c method could be problematic, with the exception of conjunctions that occur near the poles. It is expected the number of conjunction will increase significantly in the future, considering the proposed plan for the evolution of Swarm constellation, that envisages a distinct phase when the orbital planes of swB and swA/swC reaches an angular separation of 180 • , with swB and the lower satellites orbiting the Earth in opposite directions.

Outlook
One could think of several possible directions to continue the present work and improve on the characterization of near-Earth current system using Swarm data: -design and produce quality indicators that evaluate how individual auroral crossing events comply with the underlying assumptions of FAC estimation methods. For example, a low correlation between swA and swC large-scale magnetic field perturbations may indicate events where the dual-s/c method yields unreliable results, see Trenchi et al. (2019). For such events the MVA typically provides small eigenvalue ratio, which implies a low degree of planarity for the current structures. In the same time, for events where the eigenvalue ratio is sufficiently high, the information on current sheet inclination could be incorporated into the method to correct the estimations or to fine-tune the processing parameters (e.g., quad geometry) as discussed in section 3.4. -take advantage on the flexible geometry capability in Vogt et al. (2013) dual-s/c method by using different along-track separations instead of the fix, 5-s separation currently employed. This will reveal information on different FAC scales in latitude. -a study on the spatial coherence between the small-scale fluctuations observed by swA and swC would be highly beneficial for the dual-s/c FAC estimation, particularly for the configurations that rely on simultaneous satellite positions. Ideally, this would provide some criteria/parameters to use unfiltered, or only slightly filtered, 1-Hz resolution data in the method. Related to this subject are the works of Bunescu et al.  Vogt et al. (2009). Currently, we simply compute (∇×dB) along the direction normal to the spacecraft plane but the method allows the computation of full, that is, 3-D spatial gradients based on a range of possible geometrical or physical constraints (gradient parallel or perpendicular to a certain direction, stationarity assumption, etc.). -when swA and swC orbits cross at a point far from the auroral zone, and the magnetic field there has a steady and quiet aspect, the observed FAC solution (not shown) behaves similarly to the curves shown in Figure 3g. The simple form of equation (B3) could then be used, in principle, to acquire information on magnetic field biases present in the data and, provided that the biases are nearly constant, to remove their influence on other FAC intervals that occur near the poles. -explore the possibility to use the three-s/c technique at lower latitude, where the dual-s/c method cannot estimate FACs because the B lines are close to the quad plane. The current system in that region is more complicated since besides the FACs, for example, the interhemispheric FACs and the FAC related to the F region dynamo currents, other currents may exist, for example, gravity driven and diamagnetic currents (see, e.g., Alken, 2016;Lühr et al., 2003Maus & Lühr, 2006;Olsen, 1997;Park et al., 2011, and the references therein). Nevertheless, with the three-s/c method one could perhaps infer an additional component of the local current density vector (besides the radial component, provided by the dual-s/c method) that characterizes this complex current system. -compare the results of the methods presented here with the works of de Keyser et al. (2007) and de Keyser (2008). In these papers, the authors developed a generalized LS method, designed to work with an arbitrary number of satellites, to estimate the gradients of physical quantities by assuming time-stationarity over a certain time scale. The impact of different error contributions (imprecision of measurements, small-scale max ∕ int , where max and min are the eigenvalues of R arranged in descending order. The CN is a factor that enters in the FAC error estimation scheme; its values are presented on the right part of Figure A1 for various geometry configurations. There = ∕L measures the configuration skewness (deviation from an equal-sided trapezoid) and = M∕L measures the along-track separation relative to the across-track distance. As can be seen, the smallest error amplification factor is obtained from a square configuration (i.e., = 0 and = 1).

Appendix B: Accuracy of LS Dual-s/c Method
Following the classification proposed in a number of studies (Chanteur, 1998;Vogt & Paschmann, 1998;Vogt et al., 2009Vogt et al., , 2013 that deal with the accuracy of linear gradient estimators in multipoint mission, the various sources of error can be classified in (a) errors of measurement, related to the instruments intrinsic accuracy; (b) errors in spacecraft positions; and (c) errors due to deviations from linearity. In the case of the four-point quad method, which uses virtual spacecraft positions, an additional source of uncertainty is related to (d) time nonstationarity and/or unsteadiness of the structure that is sampled.
The positional errors are less important in case of Swarm since |∇B| r ≪ B, where r and B are the rms errors in position and, respectively, magnetic field (see Vogt et al., 2013). The effect of time nonstationarity is reduced at small scales by low-pass filtering the magnetic field perturbation, but it can still be present at larger scales. For a moving time stationary structure, the gradient analysis is affected by the spacecraft relative velocity. Considering the Swarm orbital velocity of ∼8 km/s, for a FAC structure velocity of 0.2 km/s, that is, an average value for K p > 4 conditions according to Wang et al. (2009), an error of around ±3% has been found by numeric simulations. Errors due to nonlinear variations are difficult to assess, but a cautious approach requires a quad characteristic dimension sufficiently small, as compared to the spatial scale length of the analyzed structure. In case of Swarm, the latter type of error are expected to affect more the results provided by three-s/c method, where distances between the sensors are larger.
In Vogt et al. (2009Vogt et al. ( , 2013Vogt et al. ( , 2019 an error calculation scheme has been proposed for dual-s/c and three-s/c gradient estimates affected by instrumental random noise. Assuming mutually uncorrelated and isotropic measurement errors, that is, < B B > = ( B) 2 I, with , sensor indexes, the Kronecker symbol, and I the unit matrix, then the mean square error of the radial current j n through the spacecraft plane is provided by with L 2 r = 1∕4 ∑ |r | 2 the mean square inter-spacecraft distance and L −2 q = ∑ |q | 2 the gradient estimation error length. In terms of quad parameters (see Figure A1) one has L 2 r = L 2 [ 1 + 2 + 2 + 2 2 ] = L 2 + M 2 + 2 + 2 M 2 L −2 q = 1 4L 2 2 ( 1 + 2 + 2 + 2 2 1 + 2 ( 2 + 2 ) ) The first term in equation (B2) captures the inverse dependence of the gradient estimation errors on array size. In the dual-s/c case it takes the highest values around the poles where L is smallest. The second term gives the contribution from the quad shape, which varies as well along the orbit, assuming minimum values when the quad takes a square configuration, that is, at latitudes where the along-and cross-track distances are equal. There is a close relation between the second term and the condition number CN, defined in Appendix A, namely, L 2 r ∕L 2 q = 1∕4(1 + CN + CN −1 ). For moderate values of CN, the dependence can be approximated by CN = 4L 2 r ∕L 2 q − 3 (see Vogt et al., 2019). Any measurement offsets B off will produce a current density j off to be evaluated as in the previous section. Equivalently, ifn designates the direction normal to the configuration plane, a more compact expression 10.1029/2018JA026439 (valid also in the three-s/c case) can be used, that is, 0 off =n · ∑ q × B off, (Vogt et al., 2009(Vogt et al., , 2013. A simple formula can be obtained for a nonskewed configuration, when satellites move on nearly parallel tracks (small and ): if one makes the reasonable assumption of constant magnetic bias B off, a and B off, b that affects the vector magnetometers on satellites a and b, respectively, than outside a small interval centered on the orbital cross-point one has that is, the influence is confined to a small region (L not too large) and only the along-track components of B off matter. In terms of boundary integral, (B3) denotes the contributions along the segments parallel to v axes. The contributions along u practically cancel out since the same magnetic value, that is, (B u off, a +B u off, b )∕2, is integrated on segments of roughly the same length. acknowledges financial support from the Deutsche Forschungsgemeinschaft (DFG) through Grant VO 855/4-1 MuSICAL in the context of the DFG Priority Programme SPP 1788 DynamicEarth. The Swarm Level 1b and Level 2 data are publicly accessible through the mission web page provided by ESA (i.e., https://www.esa.int/Swarm). We acknowledge the work of many people involved in data production, calibration, and scientific validation.