Nonlinear Magnetic Gradients and Complete Magnetic Geometry From Multispacecraft Measurements

Topological configurations of the magnetic field play key roles in the dynamics of space plasmas. This study presents a novel algorithm that can estimate the quadratic magnetic gradient as well as the complete geometrical features of magnetic field lines, based on magnetic field and current density measurements by a multiple spacecraft constellation at four or more points. The explicit estimators for the linear and quadratic gradients, the apparent velocity of the magnetic structure and the curvature and torsion of the magnetic field lines can be obtained with well‐predicted accuracies. The feasibility and accuracy of the method have been verified with two demonstrations. The algorithm has been successfully applied to exhibit the geometrical structure of a flux rope. This algorithm has wide applications for studying a variety of magnetic configurations in space plasmas.

Attempts to partially solve this problem have used physical constraints to assist the complete determination of the magnetic gradients (Vogt et al., 2009). The symmetries in plasma structures and the electromagnetic field laws can also be useful. It has been found by Shen, Rong, Dunlop et al. (2012) that, for a force-free magnetic structure in which the current is field-aligned, the 3-D magnetic gradient can be completely obtained with 3 spacecraft magnetic measurements. In their derivation, Ampere's law     0 B j and the solenoidal condition of the magnetic field   0 B are used to reduce the equations. Furthermore, if the force-free magnetic structure is steady and moving with a known relative velocity, only two S/C magnetic observations are needed to gain the complete nine components of the linear magnetic gradient (Shen, Rong, & Dunlop, 2012). Y. Y. Liu et al. (2019) have suggested a method to get the nonlinear distribution of the magnetic field in a stationary plasma structure by fitting the second-order Taylor expansion based on 4 S/C magnetic measurements and one S/C current density observation, where Ampere's law and the solenoidal condition of the magnetic field was also used. Torbert et al. (2020) have successfully obtained the 3D distribution of the magnetic field by using the four-point magnetic and particle/current density measurements of MMS. In their exploration, they have applied a fitting method to the magnetic field to the third order in magnetic gradient (( )    m l k j B 0 ), called the "25-parameter fit." In addition, the second order derivative of the magnetic field at the minimum variation direction in the local MVA characteristic coordinates was regarded as zero. Denton et al. (2020) have reconstructed the magnetic structures with a fitting method in the local characteristic coordinates system  and explored the contributions of various linear and quadratic terms in detail. However, there still exists no explicit solution to the determination of the quadratic magnetic gradient based on multiple spacecraft measurements.
With multiple S/C magnetic observations, geometrical features of the magnetic field lines can be obtained (Rong et al., 2011;Shen et al., 2003Shen et al., , 2011Shen et al., , 2014Shen, Liu, et al., 2008;Shen, Rong, et al., 2008;Lavraud et al., 2016;Xiao et al., 2018). The geometry of the magnetic field lines (MFLs) so obtained includes the tangential direction (just the direction of the magnetic vector), principal direction (along the curvature vector), binormal vector (the normal of the osculating plane of one MFL), curvature, and torsion. However, the torsion of the MFLs has not been obtained in these previous methods. The reason for this is that the torsion of the MFLs depends on the quadratic magnetic gradient, which needs 10-point S/C magnetic measurements (Chanteur, 1998) to be deduced. Therefore, it is necessary to explore the calculation of the torsion of MFLs based on observations of a limited number of S/C, in order to gain this more complete picture of MFLs in space.
This problem is addressed herein, where an explicit algorithm has been derived to estimate the quadratic magnetic gradient as well as the complete geometrical parameters of the MFLs based on measurements with a limited number of spacecraft. This approach has a wide range of applications for analyzing the magnetic structure in space plasmas.
In Section 2 and Appendices A and B, the algorithm for the quadratic magnetic gradient is derived. In Section 3 and Appendix C, the general differential geometric properties of magnetic field lines are addressed. In Sections 4, the new approach has been demonstrated on the Harris current sheets and cylindrical flux ropes in comparison with the analytical features of these two kinds of magnetic structures illustrated in Appendices D and E. Section 5 presents an actual application of this algorithm for analyzing the magnetic structure of one flux rope observed by MMS. Section 6 gives the summary and discussions.

The Estimators for the Linear and Quadratic Gradients of Magnetic Field
It is very important to obtain the quadratic gradient of the magnetic field. With it, we can grasp more accurately the structure of the magnetic field and, uncover the complete geometrical structure of the MFLs, including the Frenet coordinates and curvature, as well as the torsion. In this section, we obtain the explicit estimator of the quadratic magnetic gradient based on magnetic field and current density measurements from a multi-S/C constellation.
We present the derivation of this algorithm as follows.
The configuration of the four-spacecraft constellation (Cluster or MMS) is illustrated in Figure 1. (1) In this study, the Greek subscripts or superscripts apply to spacecraft and      , , , 1, 2, 3, 4; while the Latin subscript c indicates the barycenter.
The apparent motional velocity of the magnetic field structure relative to the S/C constellation reference frame is denoted as V , which may vary from point to point (Hamrin, et al., 2008). The velocity of the S/C constellation relative to the proper reference frame of the magnetic field structure is    V V. We establish the Cartesian coordinates ( 1 2 3 , , x x x ) in the S/C constellation reference, and choose the 3 x axis along the direction of    V V with its basis   3 / V x V . The configuration of the S/C constellation is characterized by the volume tensor, which is defined (Harvey, 1998;Shen et al., 2003)  . (2) Here, and in the text below, we apply some Latin subscripts or superscripts (other than c) to denote Cartesian coordinates with  , , , , , 1, 2, 3 i j k e m n and  , , , 1, 2 p q s r .
(i) The linear gradients of the magnetic field and current density at the barycenter As the MMS S/C cross a magnetic structure, the four S/C measure the magnetic field with high accuracy and time resolution Russell et al., 2016). The magnetic field observed by the th S/C at The MMS S/C can measure the distributions of ions and electrons with sufficient accuracy to yield the local current density (Pollock et al., 2016;Torbert et al., 2015Torbert et al., , 2020 as To obtain the magnetic field and its first order gradient at the barycenter of the MMS constellation, we first neglect the second order magnetic gradient under the linear approximation. With four S/C, simultaneous magnetic observations, the magnetic field and its linear gradient at the barycenter of the S/C constellation can be obtained with the previous methods established by Harvey (1998) and Chanteur (1998). In order to suppress the fluctuating components in the magnetic field and obtain the magnetic gradient at higher accuracy, we make use of the time series of the magnetic observations by the four S/C to get the magnetic gradient with the method first put forward by De Keyser, et al. (2007). In their approach, the time series SHEN ET AL.

10.1029/2020JA028846
3 of 32 , , x x x ) are the Cartesian coordinates in the S/C constellation reference frame. The S/C constellation is composed of four spacecraft (the number of spacecraft can be more 4), whose barycenter is at the point C. C has been set as the origin of the coordinates in the constellation reference frame. The apparent motional velocity of the magnetic field structure relative to the S/C constellation reference is V . Conversely, the velocity of the S/C constellation relative to the proper reference of the magnetic field structure is    V V.
data of the four S/C do not need to be synchronized. Appendix A gives the explicit estimator of the linear gradient of magnetic field in space and time with corrections from the quadratic gradients.
Based on Equations A14 and A15 in Appendix A, the magnetic field and its first order derivatives at the barycenter of the MMS constellation under the linear approximation are as follows. . The barycenter coordinates (Harvey, 1998)   In the above Equations 5-8, the accuracy is found to first order due to omission of the second order gradients. We will correct the magnetic field and its first order derivatives at the barycenter with the second order derivatives of the magnetic field according to Appendix A and will further obtain the corrected quadratic magnetic gradient by iteration (see (vii) later). The corrected magnetic field and its first order gradient at the barycenter will then have second order accuracy.
In this investigation, we have neglected the magnetic gradients with orders higher than two, so that the current density can be regarded as linearly varying. According to Equations A14 and A15 in Appendix A, the current density at the barycenter is and the linear gradient of the current density at the barycenter is of which the component form is Generally, the electron and ion measurements have different time resolutions. So that the electron and ion current densities and their linear gradients at the barycenter can be first calculated separately with Equations 9 and 10, and finally added to obtain the total current density and its linear gradient at the barycenter.
The gradient of the time derivative of the magnetic field is equivalent to the time derivative of the magnetic gradient, that is, (iii) The transformations between the temporal and spatial gradients of the magnetic field in different reference frames This approach will make use of the proper reference frame of the magnetic structure, so as to determine the second order gradient in the direction of the apparent motion of the magnetic structure, that is, the longitudinal quadratic gradient of the magnetic field. To do this, we need to find the apparent velocity V of the magnetic structure relative to the spacecraft constellation. For space plasmas, this relative velocity is much less than the speed of the light in vacuum, that is,  c V . Shi et al. (2006) have first obtained the velocity of the magnetic structure relative to the spacecraft with the temporal and spatial variation rates of the magnetic field under the assumption of stationarity. Hamrin et al. (2008) have obtained the apparent velocity of the magnetic structure using a proper reference frame. Here we give a concise discussion on the transformations between the temporal and spatial gradients of the magnetic field in different reference frames.
The time and space coordinates (t, r) in the S/C constellation reference frame and the corresponding time and space coordinates (  t ,  r ) in the proper reference frame of the magnetic structure obey the Galilean transformations, that is,   t t,    r r Vt (see also Figure 1) (The Eulerian description is applied in each reference frame.). The magnetic fields observed in the S/C constellation frame and the proper frame of the It is obvious that the magnetic gradient in these two reference frames are also identical, that is, The relationship between the time derivative of the magnetic field in the S/C constellation, , and time derivative of the magnetic field in the proper reference frame of the magnetic structure, which is the same formula as given by Song and Russell (1999) and Shi et al. (2006).
In the proper reference frame of the magnetic structure, The component form of the above formula is The above equation has a unique solution of the apparent velocity and a proper reference frame can be . Thus the apparent velocity of the magnetic structure relative to the S/C constellation is (Hamrin et al., 2008;Shi et al., 2006) It is noted that the apparent velocity of the magnetic structure can vary with time.  / V V is a characteristic, directional vector, so that we can define  / V V as the directional vector of the 3 x axis in the S/C constellation reference frame, that is,   We can further investigate the transformation between the time derivatives of the magnetic gradients in the two different reference frames. Similarly to the linear magnetic gradients in Equation 13, the quadratic magnetic gradients in the S/C constellation frame and the proper frame of the magnetic structure are identical, that is, The relationship between the time derivative of the magnetic gradient in the S/C constellation frame,     , t t B r , and the time derivative of the magnetic gradient in the proper frame of the magnetic structure, in the proper reference frame and Equation 19, this reduces to 10.1029/2020JA028846 6 of 32 which is the formula relating the time derivative of the linear magnetic gradient to the quadratic magnetic gradient in the S/C constellation reference frame. With this general formula the gradient of the linear magnetic gradient in the direction of apparent velocity is readily obtained as shown below in (iv).
(iv) The longitudinal gradient of Based on Equation 20, the gradient of the linear magnetic gradient along the 3 x direction at the barycenter c r satisfies The right hand side of the above equation can be obtained from Equation 12, so that nine components of the quadratic magnetic gradient can be obtained. Formula 22 is applicable for both steady and unsteady magnetic structures.
Furthermore, due to the symmetry of the quadratic gradient, of which the right hand side is given by Equation 20, so that 6 more components of the quadratic magnetic gradient can be obtained. Now only       p q l , 1,2, 1, 2, 3 B p q l are to be found, which involve   4 3 12 components. Considering the symmetry of the quadratic magnetic gradient,      p q l q p l B B, only   3 3 9 of these components are independent.
The gradient of the current density will be needed for the estimation of the remaining components of the quadratic magnetic gradient.
(v) Three components and two constraints for the quadratic magnetic gradient using the gradient of current density From Ampere's law, we get the constraints that with which we can obtain some components of the quadratic magnetic gradient if j is known (for simplicity, we replace  0 j by j.). If the electromagnetic fields are strongly varying, , with the electric displacement current included. However, in this investigation we only consider the slow-varying electromagnetic fields with the limitation j is not an independent constraint due to Equation 22. It is a surplus condition, which we have not used because Equation 22 can yield the longitudinal gradient directly already. Furthermore, , so that the gradient of the current density only provides 9 − 3 − 1 = 5 independent constraints.
The transverse quadratic gradient of the longitudinal magnetic field, that is, the quadratic gradient of the magnetic component 3 B in the plane orthogonal to the direction of motion (or where again Ampere's law    B j has been used. Thus, Equation 24 leads to where  p j q is used. The above formula yields the transverse quadratic magnetic gradient of the longitudinal magnetic field and contains three independent components of the quadratic magnetic gradient at the barycenter.
SHEN ET AL.
10.1029/2020JA028846 7 of 32 There are still 6 components of the quadratic magnetic gradient remaining to be determined, that is, q s c B t r , which are the transverse quadratic gradients of the transverse magnetic field.
Two additional constraints can be obtained from which is at the barycenter.
Based on Ampere's law, therefore, three more components of the quadratic magnetic gradient and two constraints on it can be obtained with the gradient of current density as shown in Equations 25-27. Now four constraints are to be found for the complete determination of the quadratic magnetic gradient.

(vi) The last four constraints
The magnetic field is divergence-free, that is, It is noted that the sum over k is made in the above formula.
, there are only two independent constraints, that is, There are therefore only two constraints left to be found.
Using magnetic rotation analysis (MRA) , see also Appendix B), the remaining two constraints can be obtained from the properties of the magnetic field. As shown in Appendix B, based on MRA, the magnetic rotation tensor has three characteristic directions ( 1 X , 2 X , 3 X ), as illustrated here in Figure 2. The coordinate line 3 X is along 3 X . In the third characteristic direction 3 X , the magnetic unit vector  B B b has no rotation, and the square of the magnetic rotation rate is at each point of the coordinate line 3 X (as indicated in Figure 2), we have Since the magnetic unit vector b obeys  ˆ1 b b , the above constraint contains only two independent component equations, which can be chosen as The three characteristic directions ( 1 X , 2 X , 3 X ) have a relationship with the base vectors ( x ) of the S/C coordinates ( 1 2 3 , , x x x ), as follows: The first order partial derivative obeys: and the second order partial derivative obeys: is omitted in the above equations. Therefore, Equation 34 reduces to x is orthogonal to both 3 x and 3 X , that is, (and as illustrated in Figure 2). Thus All the terms in the right hand side of the above equation are known. With Equation 59 developed in the next section, we can express the second order gradients of the components of the magnetic unit vector on the two sides of Equation 38 in terms of the magnetic gradients. With Equation 59, we get where p = 1, 2. All the terms in the right hand side of the above equation can be determined with (59), (8), and The above two equations are valid at the barycenter.
In addition, from Equations 29 and 26, we can obtain and The above two equations are also valid at the barycenter.
So far, we have obtained all the components of the quadratic gradient    c B at the barycenter. The accuracy of the quadratic gradient is to first order, just as that for the magnetic gradient.

(vii) Recalculating the magnetic gradients by iteration
In order to enhance the accuracy of the magnetic quantities, we can correct the estimate of the field and its linear gradient at the barycenter with the quadratic magnetic gradient obtained above (based on Equations A8 and A13 in Appendix A). Subsequently, we can further go through the above steps (ii-vi) to get the corrected quadratic magnetic gradient with better accuracy.
The procedure is as follows: By Taylor expansion as in Equation A3, the magnetic fields measured by the four spacecraft are Based on Equation A8 in Appendix A, we obtain the magnetic field at the barycenter, corrected by the quadratic magnetic gradient, as: where, the general volume tensor  R is as defined in Equation A9.
in Appendix A, we get the first order magnetic gradient at the barycenter corrected from the quadratic magnetic gradient as 10.1029/2020JA028846 10 of 32 Furthermore, we can perform the above steps (ii-vi) to obtain the corrected quadratic magnetic gradient using these updated estimates. The quadratic magnetic gradient obtained in this iterative sense has a higher accuracy, while errors in the magnetic field, its linear gradient and the apparent velocity of the magnetic structure at the barycenter, are of second order in L/D, where L is the size of the S/C constellation and D is the characteristic scale of the magnetic structure.
We can call this new method as Nonlinear Magnetic Gradient (NMG) algorithm for convenience.
To summarize the NMG algorithm, we proceed as follows determine the three characteristic directions ( 1 2 3,, X X X ) using MRA, and define the coordinate base , such as to fix the Cartesian coordinates    It should be noted that the magnetic field, the linear magnetic gradient and the quadratic magnetic gradient are all identical in different reference frames. We will test all these estimators in detail in Section 4.
Given the magnetic field c B , the first order magnetic gradient    c B and the quadratic magnetic gradient    c B , the complete geometry of the magnetic field lines of the magnetic structure can be determined. We will find the estimators for the geometrical parameters of the MFLs in the next section.

Determining the Complete Geometry of Magnetic Field Lines Based on Multiple S/C Measurements
The geometry of the MFLs plays a critical role in the evolution of the space plasmas. In this section, we will extract the estimators for the complete geometry of the MFLs, from the linear and quadratic gradients of the magnetic field estimated in Section 2. SHEN ET AL.

The Natural Coordinates and Curvature of the MFLs
The geometry of the magnetic field lines is illustrated in Figure 3. The directional magnetic unit vector is , which is also the tangential vector of the MFLs. The MFLs are usually turning, and the bending of MFLs is characterized by the curvature vector, that is, where "s" is the arc length along the MFLs. Shen et al. (2003Shen et al. ( , 2011 first presented the estimator of the curvature of MFLs, which has found many applications in multi-point data analysis.
Here, a brief description of it is given and we will then investigate further the complete geometry of the MFLs as well as the explicit estimators.
The gradient of the magnetic field    c B at the barycenter from multi-spacecraft measurements has already been expressed in Section 2.
The gradient of the magnetic strength B = |B| is while at the barycenter of the S/C constellation, Similarly, the gradient of the unit magnetic vector b is with Equation 49, the above Equation 51 becomes Hence, the gradient of the unit magnetic vector b at the barycenter is All the coefficients on the right hand side of the above formula involve values at the barycenter (Shen, et al., 2003): B . Equation 53 obeys the condition that:   The curvature of the MFLs at the barycenter is therefore All the coefficients on the right hand side of the above formula involve values at the barycenter.
, indicating that the obtained curvature vector is orthogonal to the magnetic field.
The radius of the curvature of the MFLs is The principal normal vector of the MFLs is is the magnetic unit vector; κ is the curvature vector of the magnetic field line, K and N are the principal normal and binormal, respectively. The magnetic field line is also twisting with torsion.

The binormal vector of the MFLs is
The above expressions collectively describe the estimators of the magnetic curvature analysis approach (Shen et al., 2003, where   , ,b K N constitute the natural coordinates, or the Frenet frame (trihedron).
The unit magnetic vector b , principal normal K , and binormal N are orthogonal to each other.
Usually, the MFLs not only bend, but also twist, such as the helical MFLs manifested in a flux rope. The twisting of the magnetic field is closely related to the magnetic helicity of the global magnetic structure, which is generally conserved in space (Berger, 1999). The local twist of the MFLs can be described quantitatively by the torsion. In order to get the complete geometry of the MFLs, therefore, the torsion should be known. The torsion of the MFLs is defined as Therefore, the quadratic gradient of the magnetic field B is essential for the calculation of the torsion of the MFLs.
We now investigate the relationship between the torsion of the MFLs and the quadratic gradient of the unit magnetic vector b; as well as with the quadratic magnetic gradient B.
To do this, we need to first deduce the expression of the quadratic gradient of the unit magnetic vector in terms of the linear and quadratic magnetic gradients.
The quadratic gradient of the unit magnetic vector b is Thus, the estimator of the quadratic gradient of b at the barycenter is expressed as Based on this definition, the torsion of the MFLs is So that the torsion of the MFLs at the barycenter of the S/C constellation is The above formula is one of the estimators of the torsion of the MFLs that is dependent on the linear and quadratic gradients of the unit magnetic vector b .
SHEN ET AL.
where the condition  b 0 j j N has been used. Appendix C presents another verification of the expression (63) for clarity. It seems that Equation 63 is invalid as B = 0 or   0. However, there is no field line as B = 0, while for   0, the field line is a straight line and its torsion has no fixed value, and thus is meaningless.
Therefore, the torsion of the MFLs at the barycenter can be written as All the coefficients on the right hand side of the above formula involve values at the barycenter. Equation 64 is another estimator of the torsion of the MFLs, expressed in terms of the linear and quadratic gradients of the magnetic field. The two different estimators of the torsion of the MFLs (62) and (64) are obviously equivalent.

Demonstrations
In this section, the estimators put forward in Sections 2 and 3 will be demonstrated for model current sheets and flux ropes, which can occur in the magnetosphere, in order to verify the validity and accuracy of this approach. A one-dimensional Harris current sheet model (Harris, 1962) and a Lundquist-Lepping cylindrical force-free flux rope model (Lundquist, 1950) are used for these two typical structures, respectively. Appendices D and E present, analytically, the geometrical features of these two kinds of magnetic structures. The tests below have shown that the estimators of the quadratic magnetic gradients and the complete geometry of the MFLs are obtained with good accuracy compared to the models, so we expect they can find wide applications in investigating the magnetic structures and configurations in space plasmas with multi-S/C measurements.
The operative calculating steps can follow the steps (a-g) summarized at the end part of Section 2. The curvature and torsion of magnetic field lines can be obtained by Equations 54 and 64 (or 62) in Section 3.

One-Dimensional Harris Current Sheet
For the one-dimensional Harris current sheet, the magnetic field can be formulated as Equation (D1)  . As shown in Figure 4a, we set an arbitrary S/C constellation trajectory from (2, 2, 2) R E to (−2, −2, −2) R E during 100 s. The S/C constellation is assumed to be a regular tetrahedron with a separation of L = 100 km. The analytic values of the magnetic field and the current density at the barycenter of the four S/Cs are shown in panels (b) and (c) in Figure 4, respectively.
In this test, we have set n = 10, and make N = 4n = 40 points to calculate the spatial and temporal gradient of the magnetic field at the barycenter of the S/C constellation with the method in Appendix A. Therefore, we can get the spatial gradient of the vector field within the interval 5-95 s. Furthermore, the temporal and spatial scale corresponds to the time resolution of the field sampling (i.e., T = 1s) and the characteristic size of tetrahedron (L = 100 km). The magnetic field and the nonzero component   x B z of the linear magnetic gradient at the barycenter are derived with Equations 5 and 6 and shown in Figures 4b and 4d, respectively, which are in good agreement with their analytic values as given by Appendix A (the circles represent the results derived by the method, while the black solid line denotes the analytic results derived by theoretical formula). The current density at the barycenter can also be derived with   After completing the above steps, iterative operation and error analysis are necessary and we will discuss these later.

Two-Dimensional Force-Free Flux Ropes
In this subsection, we attempt to investigate the complete geometry of magnetic field lines for a classic force-free flux rope model. In this model, the three components of the magnetic vector in cylindrical coordinates can be expressed (Lundquist, 1950) as: where r is the distance from the central axis,  1 is the characteristic scale of the flux rope, and J is the Bessel function. In this test, we adopt  0 60 B nT and   1 / E R . Then the characteristic scale of the flux rope is The trajectory of the 4 SCs with a separation of L = 100 km is set to be from (−2, 0, 0) R E to (2, 0, 0) R E during 100 s (indicating the actual structure velocity is 257 km/s [0.0404 Re/s]), which is shown in Figure 5a. The average magnetic field measured by four S/C is illustrated in Figure 5b, the bipolar signature of By and the enhancement of Bz around the flux rope's center is apparent.
By repeating the same procedures as in Section 4.1, the quadratic magnetic gradient can be readily acquired (Figures 5c-5g). One can find that the results derived by the method are in good agreement with the analytic results obtained in Appendix E. The variations of curvature and torsion of the MFLs confirm that the magnetic topological structure is different from those of the current sheet (Figures 5h and 5i). It can also be seen from Figures 5h and 5i that the curvature of the MFLs contains a minimum, and the torsion of the MFLs contains a maximum, at the center. This indicates that the straighter and more twisted the MFLs, the nearer to the center of flux rope, implying the nonplanar and helical structure of the flux rope. This test shows that the results obtained by the approach are in good agreement with the analytical results, indicating that the estimators obtained in Sections 2 and 3 are reliable and applicable. SHEN ET AL.

Errors
The errors of the estimators put forward in this study may arise from two types of sources: the underlying measurement errors and the truncation errors. The key measurement errors include the error in the measured magnetic field B and that of the current density j derived from the plasma moment data (which will be seen in the application in Section 5). The truncation errors arise from terms beyond the differential order considered here and represent neglected behavior of the magnetic structure and plasmas.
The spatial truncation errors can be approximately measured by L/D, where D is the typical spatial size of magnetic structure and L is the size of tetrahedron of four SC. When L/D is very small, the truncation errors are generally small. However, as L/D grows large, the truncation errors may become significant. The iterative operation allows us to attempt to get more accurate and reliable results. Figure 6 compares the results of the calculations made with no iteration; with the first and second iterations, and theoretical calculation with L/D = 0.3. It can be seen that the iterations yield more accurate results. However, the second iteration in these examples did not produce better results than the first iteration.  (Figures 7f-7j), although after the first or second iterations they are improved.
Through the above analysis, one can conclude that the most accurate results are those derived by at least one iteration, especially when L/D is larger than 0.5. Thus, it is necessary to perform the first iteration when L/D is larger than 0.5.
In addition, we have performed a calculation on the modeled geomagnetosphere with no symmetry and presented the obtained results in the supplementary information, where the structure of the cusp regions has been demonstrated reasonably as in Figure S1 and S2.

Application: Magnetic Flux Rope
In this section, we have applied the approach developed in Sections 2 and 3 to investigate the magnetic structure and geometry of a magnetic flux rope at magnetopause, observed by MMS during 2015-10-16 13:04:33-13:04:35, which is the second of two sequential flux ropes reported by Eastwood et al., (2016), and has been analyzed by many researchers (e.g., C. Zhang et al., 2020). Here, we have used the high-resolution magnetic field data measured by the fluxgate magnetometer, operating at 128 vectors per second in burstmode Russell et al., 2016), and the plasma data provided by FPI (fast plasma investigation, measuring electrons at cadence of 30 ms and ions at cadence of 150 ms) (Pollock et al., 2016;Torbert et al., 2015). To calculate the quadratic magnetic gradient, the plasma moments are interpolated to obtain a 1/128 s time resolution to match that of the magnetic field data and to derive the current density. Note that the MMS constellation is often nearly a regular tetrahedron with its separation scale of L  20 km during this time interval.
Typically, there are many waves affecting the magnetic field at various frequencies in space plasmas. If we calculated the time variation rates of the magnetic field and the linear and quadratic magnetic gradients directly, the errors caused by these waves would be so large that we would miss the underlying global features of the magnetic structure. To get rid of the influence of the waves, the magnetic field ( Figure 8a) and current density (Figures 8b-8d) data have been filtered by a low-pass filter to eliminate disturbances with frequencies higher than 1 Hz from the data. In order to apply the method in Appendix A to calculate the temporal and spatial gradients of the magnetic field and current density, we have adopted n = 10 time points on each spacecraft to form a set of data. Thus, there are in total N = 4n = 40 points in a group of data. With SHEN ET AL.

10.1029/2020JA028846
18 of 32 this approach, the calculated temporal and spatial gradients of the magnetic field and current density have rather high accuracy.
We have derived the magnetic rotation features of the flux rope by using the MRA method illustrated in Appendix B . Figure 8e shows the time series of the magnetic minimum rotation direction 3 X , which is approximately stable and nearly parallel to GSE + Y direction. Assuming the flux rope is cylindrically symmetric, 3 X could be approximately regarded as the orientation n of the flux rope axis, that is,  3 n X . The helical angle of the MFLs can be defined as As shown in Figure 8f, the SHEN ET AL.

10.1029/2020JA028846
19 of 32 , and is illustrated in Figure 8g. One can find that, the apparent velocity at the leading edge of flux rope is larger than that at the trailing edge, suggesting that the flux rope is decelerating and not stable during this interval. Assuming that the flux rope is steady and has a force-free magnetic field, Eastwood et al., (2016)  , when the flux rope is nearly steady. Considering the complicated motion and structure of flux rope and the different data processing approaches applied, the small discrepancy among the results is not surprising.
By using the estimators in Sections 2 and 3, the magnetic gradients and geometry of the flux rope can be obtained and these are demonstrated in Figure 9. The total 27 components of the quadratic gradient of magnetic field have been obtained with the estimators in Section 2, which are illustrated in panels (a-i) of Figure 9. It can be found that the order of the quadratic gradient of the magnetic field is generally less than 10 −2 nT/km 2 , while that of the first-order magnetic gradient is ∼10 −1 nT/km. The complete geometry of the SHEN ET AL.
10.1029/2020JA028846 20 of 32 MFLs in the flux rope can be derived by the estimators in Section 3, which is illustrated in Figures 9j-9l. It can be seen that the curvature of MFLs reaches its minimum value of ∼0.80 × 10 −3 /km (Figure 9j) and the torsion reaches its maximum value of ∼0.012/km 2 (Figure 9l) at ∼34.1 s, when the helical angle is the largest (Figure 8f). These features indicate that this flux rope is a typical one and is consistent with the 2-D flux rope model in Appendix E. The maximum curvature of the MFLs is about ∼3.0 × 10 −3 /km, while accordingly the minimum radius of the curvature of the MFLs is ∼330 km. We can choose this as the characteristic scale of the flux rope, that is, D = 330 km. Furthermore, assuming the flux rope has a cylindrical helical structure, the torsion of MFLs can also be obtained directly from the curvature and helical angle from Equation E9

Summary and Discussions
The quadratic magnetic gradient is a key parameter of the magnetic field, with which the fine structure of a magnetic structure can be revealed, as well as the twisting property of the magnetic field. However, up to now, the quadratic magnetic gradient from multi-S/C constellation measurements has not been explicitly calculated. Chanteur (1998) showed that in order to get the quadratic magnetic gradient from multipoint magnetic observations, in general, the number of S/C in the constellation has to be equal to or larger than 10, which is difficult to realize in present space exploration. Fortunately, the MMS constellation can not only provide rather accurate four-point magnetic field, but can also produce very good four-point current density estimates from particle measurements, such as to allow the quadratic magnetic gradient problem to be solved in the manner discussed here.
SHEN ET AL.
10.1029/2020JA028846 22 of 32 This paper provides a method to obtain the linear and quadratic magnetic gradients as well as the apparent velocity of the magnetic structure, based on the four-point magnetic field and current density observations and give their explicit estimators. Furthermore, the complete geometry of the magnetic field lines is revealed on the bases of these linear and quadratic magnetic gradients, and the estimator for the torsion of the MFLs is given. Simple, but relevant, tests on this novel algorithm have been made for a Harris current sheet and a force-free flux rope model, and the effectiveness and accuracy of these estimators have been verified.
In this approach (NMG algorithm), the physical quantities to be determined are as follows: the magnetic field c B (three parameters); the linear magnetic gradient   , as deduced from MRA (two independent constraints); resulting in a total of 12 + 3+9 + 5+2 + 2 = 33 independent parameters or constraints.
We note that the contribution of the current density measurements in this approach is the first order gradient of the current density, which is related to the quadratic magnetic gradient by Ampere's law. Considering the conservation of the current density   0 j and   3 B already obtained from the constraint equation j yields only 2 × 3 − 1 = 5 independent con- provides only 3 − 1 = 2 independent constraints. Therefore, the linear and quadratic magnetic gradients, and the apparent velocity of the magnetic structure, can be completely determined based on the four-point magnetic field and current density measured by the MMS constellation.
The calculations have been expressed as being carried out in the S/C constellation frame. The NMG algorithm proceeds as follows. First, under the linear approximation, the temporal and spatial gradients of the ) and of the current density ( Therefore, all the 18 independent components of the quadratic magnetic gradient can be calculated.
The quadratic magnetic gradient, obtained with no iteration, has a truncation error of the first order in L/D because the linear approximation has been made. To find a more accurate quadratic magnetic gradient, an iterative procedure can be performed. In this procedure, the magnetic field, the linear magnetic gradient, and the time derivative of the linear magnetic gradient are corrected by using the quadratic magnetic gradient calculated initially and the above steps are then repeated, so as to achieve the components of the SHEN ET AL.

10.1029/2020JA028846
23 of 32 corrected quadratic magnetic gradient. After this first iteration, the magnetic field, linear magnetic gradient, the apparent velocity of the magnetic structure at the barycenter of the S/C tetrahedron all have their accuracies improved significantly and have truncation errors in the second order of L/D, while the accuracy of the quadratic magnetic gradient obtained is also enhanced.
This algorithm is valid for both stationary and nonstationary structures, whether the magnetic structures are moving at a constant velocities or accelerating/decelerating. It is noted that the magnetic field, linear and quadratic magnetic gradients are identical for different inertial frames of reference.
With the magnetic field, linear and quadratic magnetic gradients found, the complete geometry of the MFLs can be determined, including the natural coordinates or Frenet coordinates (tangential unit vector, principal normal and binormal), curvature, and torsion. The corresponding estimators for the geometrical features have been given.
The algorithm for estimating the quadratic magnetic gradient and the geometry of the MFLs have been demonstrated with the Harris current sheet and cylindrical flux rope, and its correctness has been verified. It is found that, the errors of the linear quadratic magnetic gradients, apparent velocity of the magnetic structure, and the geometrical parameters are of first order in L/D when no iteration is made. If one iteration is performed, the accuracies of the linear magnetic gradient, apparent velocity of the magnetic structure, curvature of the MFLs are improved significantly and their errors appear at the second order in L/D, while the accuracies of the quadratic magnetic gradient and the torsion of the MFLs are also enhanced.
To determine the first order magnetic gradient and apparent relative velocity of the magnetic structure, this algorithm is more accurate than the previous approaches based on the linearity approximation (Chanteur, 1998;Harvey, 1998;Shi et al., 2006).
We have also applied the algorithm developed in this research to investigate the magnetic structure of one flux rope measured by MMS (Eastwood et al., 2016), showing good results. The applicability of this approach is therefore verified.
If the magnetic gradients with orders higher than two are neglected the magnetic field can be expressed as with the MMS magnetic field and current density measurements, the linear and quadratic magnetic gradients at the barycenter are obtained, such that the local spatial distribution of the magnetic field, as well as the MFLs, can be obtained.

Appendix A: The Explicit Estimators for the Linear Gradients of Field in Space and Time
De Keyser et al. (2007) has put forward an algorithm for calculating the gradients in space and time of a field, which they called Least-Squares Gradient Calculation (LSGC). Here we will find the explicit estimator of the four-dimensional linear gradients of a scalar field or one component of the vector field with corrections from the quadratic gradients.
Considering the 4 S/C of the constellation obtained time series of measurements on a certain physical quantity investigated, as illustrated in Figure 1 in Section 2. Here the S/C constellation reference frame is used. Assuming each S/C makes observations at n times, in total N = 4n measurements are made by the constellation, which form a set of data. In the previous work of De Keyser, et al. (2007), it is supposed that, in this area of spacetime, the physical quantity concerned is approximately varying linearly, and the linear gradients of field in space and time are about homogeneous. Here, however, the homogeneity assumption is not necessary because the quadratic gradients have been included. In the S/C

Journal of Geophysical Research: Space Physics
Obviously, in the S/C constellation reference frame, the four S/C are nearly motionless and their space co- x y z do not change with time during typical structure crossing events.
In the S/C constellation reference frame, at the space time          , , ; 1, 2, 3, 4 a a a a a x x y z t , the physical quantity measured is , , , . The spacetime coordinates at the central point satisfy gives the value of physical quantity f at the barycenter of the constellation with corrections by the quadratic gradients. Furthermore,

Appendix D: Geometry of the MFLs in One-Dimensional Current Sheets
Here, the properties of the one-dimensional current sheets are the update to the results in Appendix D in the previous study (Shen et al., 2003). It is assumed that the magnetic field in the 1 dimensional current sheets is   ˆx . Let the z axis to be along the normal to the 1 dimensional current sheets.
SHEN ET AL.

10.1029/2020JA028846
27 of 32 The components of the magnetic field in the x and y directions are invariants, that is,     x 0, 0 y . Therefore the components of the magnetic field in the Cartesian coordinates are Const. However, as shown in actual observations, the component y B is not constant, which maximizes at the center of neutral sheets and is decreasing away from the center of the current sheets (Rong, et al., 2012). The MFLs in the magnetotail current sheets often have a shape of helix in the neutral sheets (Shen, Liu, et al., 2008).
As shown in Figure E1, polar coordinates are used. The central axis is along the z axis, the arc length is s, the distance from the central axis is r, and the azimuthal angle is .The radial unit vector is ˆr e , and the azimuthal unit vector is  e . The tangent vector of the MFLs is On the contrary, if the curvature  and torsion  of the cylindrical spiral MFLs have been measured, the helix angle, the distance from the central axis, the spiral pitch and the rotation frequency can be expressed as Any arbitrary magnetic field line can locally be fitted by a cylindrical spiral arc with the same curvature and torsion. The curvatures of the magnetic field lines are always nonnegative. However, the torsion of one MFL can be either positive or negative. When   0, the helix angle   0, the magnetic field line is locally a right-hand cylindrical spiral; while   0,   0, it is a left-hand one.