SOTE: A Nonlinear Method for Magnetic Topology Reconstruction in Space Plasmas

Assume that the SC tetrahedron crosses a magnetic structure in space. At each sampling point, the SC tetrahedron provides four-point measurements, and between the two adjacent sampling points, the tetrahedron moves a certain distance. If the SC trajectory in each sampling interval could be precisely determined, the continuous trajectory of SC can be obtained from the integration. In the GS reconstruction (Sonnerup & Guo 1996), the negative de Hoffmann–Teller velocity is used as a constant traversing velocity, which may cause considerable errors when crossing a dynamic structure. Here, we use a more accurate method to estimate the SC position during the crossing, based on magnetic-field measurements (Shi et al. 2006).

Assume that at time t = 0, SC1 is located at the origin of the coordinate system, with SC2, SC3, and SC4 near it. The SC1 measurement is expressed as ${{\boldsymbol{B}}}_{1,{\boldsymbol{t}}=0}$, and then Equation (1) gives

$$\begin{eqnarray}&&{\boldsymbol{B}}\left(0,0,0,0\right)={{\boldsymbol{B}}}_{1,{\boldsymbol{t}}=0}.\end{eqnarray} \tag{ 2 }$$

After a sampling interval, i.e., at t = △t, the SC tetrahedron moves a bit. The SC1 leaves the origin and arrives at △r = (△x, △y, △z). We express its measurement now as ${{\boldsymbol{B}}}_{1,{\boldsymbol{t}}=\bigtriangleup {\boldsymbol{t}}}$. Considering no temporal evolution of magnetic fields, Equation (1) gives

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{B}}}_{1,{\boldsymbol{t}}=\bigtriangleup {\boldsymbol{t}}} & = & {\boldsymbol{B}}\left(\bigtriangleup x,\bigtriangleup y,\bigtriangleup z,\bigtriangleup t\right)={\boldsymbol{a}}\bigtriangleup x+{\boldsymbol{b}}\bigtriangleup y+{\boldsymbol{c}}\bigtriangleup z\\ & & +\,{\boldsymbol{d}}\bigtriangleup x\bigtriangleup y+{\boldsymbol{e}}\bigtriangleup x\bigtriangleup z+{\boldsymbol{f}}\bigtriangleup y\bigtriangleup z\\ & & +\,{\boldsymbol{l}}{\left(\bigtriangleup x\right)}^{2}+{\boldsymbol{m}}{\left(\bigtriangleup y\right)}^{2}+{\boldsymbol{n}}{\left(\bigtriangleup z\right)}^{2}+{\boldsymbol{B}}(0,0,0,0).\end{array}\end{eqnarray} \tag{ 3 }$$

Since the sampling interval is usually very small in SC measurements (e.g., 1/128 s for MMS); △x, △y, and △zare small quantities. In this way, the second-order terms in Equation (3) can be neglected, and consequently Equation (3) can be simplified to

$$\begin{eqnarray}&&{{\boldsymbol{B}}}_{1,{\boldsymbol{t}}=\bigtriangleup {\boldsymbol{t}}}={\boldsymbol{a}}\bigtriangleup x+{\boldsymbol{b}}\bigtriangleup y+{\boldsymbol{c}}\bigtriangleup z+{\boldsymbol{B}}(0,0,0,0).\end{eqnarray} \tag{ 4 }$$

Combining Equations (2) and (4), we get

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{B}}}_{1,{\boldsymbol{t}}=\bigtriangleup {\boldsymbol{t}}}-{{\boldsymbol{B}}}_{1,{\boldsymbol{t}}=0} & = & {\boldsymbol{a}}\bigtriangleup x+{\boldsymbol{b}}\bigtriangleup y+{\boldsymbol{c}}\bigtriangleup z\\ & = & \ \left(\begin{array}{ccc}{\boldsymbol{a}} & {\boldsymbol{b}} & {\boldsymbol{c}}\end{array}\right)\cdot \left(\begin{array}{c}\bigtriangleup x\\ \bigtriangleup y\\ \bigtriangleup z\end{array}\right),\end{array}\end{eqnarray} \tag{ 5 }$$

where (a b c) is the Jacobian matrix, which can be obtained from four-point measurements (Paschmann & Daly 1998).

Since the quantities ${{\boldsymbol{B}}}_{1,{\boldsymbol{t}}=0}$, ${{\boldsymbol{B}}}_{1,t=\bigtriangleup {\boldsymbol{t}}}$, and (a b c) are known, the SC displacement △r certainly can be resolved during each sampling interval. Finally, the continuous trajectory of the SC can be obtained from the integration of △r. To minimize the errors, we repeat this procedure using the measurements of SC2, SC3, and SC4 and finally take the average of them.

Having known the continuous trajectory of SC, the time-series data can be transformed to space-series measurements. We combine two sets of four-point measurements, i.e., the measurements at two different time instants (noted as T1 and T2), into one set of eight-point measurements. The time delay from T1 to T2 determines the relative position of the two tetrahedrons and therefore determines the accuracy of the reconstruction results, because an irregular satellite configuration always results in unstable solutions. Taking this issue into account, we propose a method to optimize the displacement between the two SC tetrahedrons, as introduced below.

First of all, we assume that at T1 the SC1 is at the origin of the coordinate system, and the positions of SC2, SC3, and SC4 are expressed as ${{\boldsymbol{R}}}_{\mathrm{SC}2,{\boldsymbol{t}}={\rm{T}}1}$, ${{\boldsymbol{R}}}_{\mathrm{SC}3,{\boldsymbol{t}}={\rm{T}}1}$, and ${{\boldsymbol{R}}}_{\mathrm{SC}4,{\boldsymbol{t}}={\rm{T}}1}$, respectively. During the interval △t (△t = T2–T1), the SC tetrahedron moves a distance along its trajectory, and the displacement is expressed as L, which is also the position of SC1 at T2. Correspondingly, the positions of other three SC are ${\boldsymbol{L}}+{{\boldsymbol{R}}}_{\mathrm{SC}2,{\boldsymbol{t}}={\rm{T}}1}$, ${\boldsymbol{L}}+{{\boldsymbol{R}}}_{\mathrm{SC}3,{\boldsymbol{t}}={\rm{T}}1}$, and ${\boldsymbol{L}}+{{\boldsymbol{R}}}_{\mathrm{SC}4,{\boldsymbol{t}}={\rm{T}}1}$, because the tetrahedron configuration does not change. In total, we have seven nonzero position vectors in the form of (X, Y, Z). According to Equation (1), we expand each vector to the form (X, Y, Z, XY, XZ, YZ, X², Y², Z²) and label them as ${\boldsymbol{X}}1,{\boldsymbol{X}}2,\,...,\,{\boldsymbol{X}}7$, respectively. These vectors constitute a matrix (${\boldsymbol{X}}1^{\prime} ,{\boldsymbol{X}}2^{\prime} ,\,......,\,{\boldsymbol{X}}7^{\prime}$), which partly determines the coefficients in Equation (1). To obtain a stable solution of the coefficients, these vectors should be as linearly independent as possible. To evaluate the independence of these vectors, we define the condition number sequence (CNS). In principle, the 9 × 7 matrix (${\boldsymbol{X}}1^{\prime} ,{\boldsymbol{X}}2^{\prime} ,\,......,\,{\boldsymbol{X}}7^{\prime}$) can be decomposed to 36 7 × 7 submatrices, if we select seven rows from all nine rows. We calculate the condition numbers for all submatrices and sort them from small to large then reserve the first 18 numbers as a CNS. Notice that a CNS is a set containing 18 elements, which could be expressed as $\mathrm{CNS}\left(i\right),i=1,2,\,\ldots ,\,18$. For a fixed T1 time, a different T2 results in a different L and consequently different CNS. The time T2, corresponding to the minimum CNS, will be the best choice.

This theory also explains why we do not add T3 to get another set of four-point measurements. In fact, the time instants T1, T2, and T3 are pretty close, and the SC tetrahedron displacements from T1 to T2 and from T2 to T3 are nearly along the same direction. In such a situation, the additional four vectors at T3 are linearly dependent with the previous vectors, and certainly the condition numbers of the matrix are very large. As a result, we can only include eight-point measurements at most, which provide 24 constraints to resolve Equation (1), as shown below:

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{B}}}_{i} & = & {\boldsymbol{a}}{x}_{i}+{\boldsymbol{b}}{y}_{i}+{\boldsymbol{c}}{z}_{i}+{\boldsymbol{d}}{x}_{i}{y}_{i}+{\boldsymbol{e}}{x}_{i}{z}_{i}\\ & & +\ {\boldsymbol{f}}{y}_{i}{z}_{i}+{\boldsymbol{l}}{x}_{i}^{2}+{\boldsymbol{m}}{y}_{i}^{2}+{\boldsymbol{n}}{z}_{i}^{2}+{{\boldsymbol{B}}}_{0},\end{array}\end{eqnarray} \tag{ 6 }$$

in which i ∈ {1, 2, 3, 4, 5, 6, 7, 8}.

There is a common problem for all reconstruction methods based on multi-time measurements, including the GS method and our SOTE method. Since the time delay △t is necessary for combination of eight-point measurements, the presence of time evolution unavoidable will lead to the time aliasing of the reconstructed maps. Thus, the reconstructed structures are required to be quasi-stationary. In space plasmas, the shock, discontinuity, magnetic hole, and dipolarization front have been thought to be quasi-stationary (during SC crossing) nonlinear structures. Therefore, the SOTE method should be useful for the reconstruction of these structures.

There are 30 unknown coefficients in Equation (1), with 24 constraints provided by the eight-point measurements. Therefore, another six constraints are required. Fortunately, Ampere's theorem,

$$\begin{eqnarray}&&{\rm{\nabla }}\times {\boldsymbol{B}}={\mu }_{0}{\boldsymbol{J}},\end{eqnarray} \tag{ 7 }$$

together with Maxwell's divergence equation,

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{B}}=0,\end{eqnarray} \tag{ 8 }$$

can be used as the constraints, because in SC measurements, the particle moments are usually available. According to Equation (1), the left side of Equation (7) can be expressed as a function containing coordinates and unknown coefficients:

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\left({\rm{\nabla }}\times {\boldsymbol{B}}\right)}_{x}=\displaystyle \frac{\partial {B}_{z}}{\partial y}-\displaystyle \frac{\partial {B}_{y}}{\partial z}={b}_{3}-{c}_{2}+\left({d}_{3}-{e}_{2}\right)x+\left(2{m}_{3}-{f}_{2}\right)y+\left({f}_{3}-2{n}_{2}\right)z={\mu }_{0}{J}_{x}\\ {\left({\rm{\nabla }}\times {\boldsymbol{B}}\right)}_{y}=\displaystyle \frac{\partial {B}_{x}}{\partial z}-\displaystyle \frac{\partial {B}_{z}}{\partial x}={c}_{1}-{a}_{3}+\left({e}_{1}-2{l}_{3}\right)x+\left({f}_{1}-{d}_{3}\right)y+\left(2{n}_{1}-{e}_{3}\right)z={\mu }_{0}{J}_{y}\\ {\left({\rm{\nabla }}\times {\boldsymbol{B}}\right)}_{z}=\displaystyle \frac{\partial {B}_{y}}{\partial x}-\displaystyle \frac{\partial {B}_{x}}{\partial y}={a}_{2}-{b}_{1}+\left(2{l}_{2}-{d}_{1}\right)x+\left({d}_{2}-2{m}_{1}\right)y+\left({e}_{3}-{f}_{1}\right)z={\mu }_{0}{J}_{z},\end{array}\right.\end{eqnarray} \tag{ 9 }$$

where x, y, and z are satellite coordinates, and a, b, c, d, e, f, l, m, and n are elements of the vector constants in Equation (1), with the subscripts indicating their positions in vectors. For simplicity, we consider a satellite at the origin, which means that x, y, and zshould be 0. In this case, Equation (9) can be rewritten as ${b}_{3}-{c}_{2}={\mu }_{0}{J}_{x}$, ${c}_{1}-{a}_{3}={\mu }_{0}{J}_{y}$, and ${a}_{2}-{b}_{1}\,={\mu }_{0}{J}_{z}$, which are three constraints of Equation (1) if the current density of ${\boldsymbol{J}}={ne}({{\boldsymbol{V}}}_{{\boldsymbol{i}}}-{{\boldsymbol{V}}}_{{\boldsymbol{e}}})$ is provided by SC measurements. Similarly, Equation (8) can be expanded to

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\nabla }}\cdot {\boldsymbol{B}} & = & {a}_{1}+{b}_{2}+{c}_{3}+\left(2{l}_{1}+{d}_{2}+{e}_{3}\right)x\\ & & +\,\left({d}_{1}+2{m}_{2}+{f}_{3}\right)y+\left({e}_{1}+{f}_{2}+2{n}_{3}\right)z=0.\end{array}\end{eqnarray} \tag{ 10 }$$

Clearly, to satisfy this equation, we need

$$\begin{eqnarray}&&\left\{\begin{array}{l}{a}_{1}+{b}_{2}+{c}_{3}=0\\ 2{l}_{1}+{d}_{2}+{e}_{3}=0\\ {d}_{1}+2{m}_{2}+{f}_{3}=0\\ {e}_{1}+{f}_{2}+2{n}_{3}=0.\end{array}\right.\end{eqnarray} \tag{ 11 }$$

This provides another four constraints. Since the eight-point measurements and Ampere's theorem have provided 27 constraints, we only need 3 of them. Together, Equations (6)–(8) constitute 30 scalar equations for the 30 unknown coefficients, which can be easily determined by the matrix calculation.

The first three columns of coefficients matrix (i.e., (a b c), see Equation (1)), have a similar format with the Jacobian matrix. Paschmann & Daly (1998) have proposed the method for determining the Jacobian matrix using four-point measurements. However, their method implies the linear assumption (i.e., ${\boldsymbol{B}}={\boldsymbol{a}}x+{\boldsymbol{b}}y+{\boldsymbol{c}}z$), thus, we cannot use this method to determine (a b c). In the second-order frame, the measurements at each point are related to all the unknown coefficients (see Equation (1)), and thereby, we have to consider all 30 equations together.

Using the particle moments measured at each point, we can obtain eight reconstruction results in total. In principle, these results should be same. However, owing to the instrument error or deviation of the two-order assumption from the reality, the eight results may be slightly different. To guarantee the uniqueness of solution, we use the particle moments measured by the SC closest to the geometric center of the heptahedron. In other words, the leading SC (relative to the magnetic structure) at the initial time (T1) is used to provide the constraint of currents.

An artificial SC tetrahedron is designed to move in the simulation box along an asymmetric helical trajectory, as shown in Figure 1(a) (see the red dots). Specifically, this trajectory is defined as $X=0.9\pi \cos (\theta +\pi /4)$, $Y=0.5\pi \sin (\theta +\pi /4)$, $Z=0.9\pi (\theta /\pi -1)$, and $\theta \in [0,2\pi ]$. During this period, four SCs form a constant tetrahedron, with inter-SC separations of ${R}_{21}=0.09\times [0,1,0]$, ${R}_{31}=0.09\times [\tfrac{\sqrt{3}}{2},\tfrac{1}{2},0]$, and ${R}_{41}\,=0.09\times [\tfrac{\sqrt{3}}{6},\tfrac{1}{2},\tfrac{\sqrt{6}}{3}]$, as shown in Figure 1(b). This tetrahedron size (∼0.09, in simulation box) is selected to be comparable with that of the actual MMS, which equals ∼0.5 ion inertial length. A distinct comparison of the scales of SC tetrahedron and the simulation box is shown in Figures 1(a)–(b) by shaded cubes. Eventually, evaluating the magnetic-field components along the trajectory of each SC, we generate a group of quadruple magnetic measurement data. Using a magnetic measurement, we calculate the SC trajectory according to the method introduced in Section 2.1 and then compare the calculation result with the exact trajectory functions to test the accuracy of our method.

Figures 1(c)–(g) show the test results. The X-axis represents the integral length of the helical trajectory (labeled $\sum | \bigtriangleup L|$), while Figures 1(c)–(e) show the coordinates of SC1, in which the dotted black lines are results from our method and the solid red lines are exact trajectory functions. Figure 1(f) presents the error, defined as the deviation of the method result from the exact SC1 position (normalized by the total trajectory length). One can see, at the beginning ($\sum | \bigtriangleup L| =0$), the method result coincides with the exact SC position perfectly, with the error being almost zero. As SC1 moves along its trajectory, the method error increases slightly. To sum up, the error stays small (<4%) during the whole interval, demonstrating that our method is accurate.

In this process, correction is required to eliminate the bad effects resulting from the quasi-2D field or kinetic-scale structures. A quasi-2D field configuration makes the Jacobian matrix incomplete, so the calculation of Equation (5) is meaningless. Also, kinetic-scale structures can cause errors within the SC tetrahedron. To estimate these bad effects, we define a quality index (QI) as the condition number of the Jacobian matrix in Equation (5), labeled as $\mathrm{QI}=\parallel J\parallel \cdot \parallel {J}^{-1}\parallel$, where J is the Jacobian matrix. Such an index is shown in Figure 1(g) (solid black line). We can see that when the QI significantly deviates from its average value, the SC displacement estimated using Equation (5) (blue line in Figure 1(g)) sharply increases and becomes unreliable (the maximum ∣△L∣ can be estimated using the plasma flow velocity). Quantitatively, we treat the SC displacement as unreliable when the QI exceeds two times its median value, QI > 2*median(QI) (see the horizontal dashed line in Figure 1(g)), and correspondingly replace these displacements using the interpolation from adjacent points. For example, in Figure 1(g), there are four intervals, in which the QI exceeds the threshold QI > 2*median(QI) (see the first, third, fourth, and fifth vertical gray shades). We automatically start the correction procedure during these intervals, and consequently, the errors in these intervals do not increase (see Figure 1(f)). However, during the interval at $\sum | \bigtriangleup L| =7$ (see the second vertical gray shade), the QI does not exceed the threshold. In such a situation, the correction procedure is not triggered, and consequently, the errors cannot be removed (see the gray shade in Figure 1(f)). Notice that when applying the SOTE method to SC measurements, this threshold should be carefully defined, in order to obtain an accurate SC trajectory relative to the magnetic structure.

An SC tetrahedron is designed to cross a flux–rope-like structure at [1.67, 1.57, and 1.57] along the −Y direction, as shown in Figure 2(a). The magnetic field measured by SC1 during the crossing is exhibited in Figure 2(b). With this magnetic field, we first calculate the SC trajectory. The result without correction is shown in Figure 2(c) (dotted lines). As can be seen, this result deviates from the exact position (solid lines in Figure 2(c)) during the interval characterized by a large QI (see the vertical shade). Utilizing the technique mentioned above, we correct the SC trajectory and show it in Figure 2(e). We find that the SC trajectory after correction coincides with the exact position well, even during the large-QI interval (see the blue shade in Figure 2(e)). This provides a reliable basis for reconstructing magnetic topology.

Zoom In Zoom Out Reset image size

At $\sum | \bigtriangleup L| =0.375$, the SC tetrahedron provides a set of four-point measurements (see the yellow shade in Figure 2(a)). Then, the tetrahedron moves along its trajectory (see the color lines in Figure 2(a)). Now we need to consider a proper delay to obtain the second set of four-point measurements. We calculate the CNSs with the delay of one to seven sampling intervals of the instrument and show them in Figure 2(f). Clearly, if we use four sampling intervals (△L = 0.218) as the delay, the CNS will be the minimum (see the black line in Figure 2(f)). In this sense, the four-sampling-interval delay is the best choice and certainly provides the second set of four-point measurements. The pink shade in Figure 2(a) marks the SC tetrahedron at that moment.

So far, we obtained eight-point measurements of magnetic field. Taking into account the plasma moments measured by SC1, we can determine all the coefficients in Equation (1) and then reconstructed the magnetic topology. Figures 3(a)–(f) present a comparison between the simulation and reconstruction in the X–Y, X–Z, and Y–Z planes, while Figures 3(m)–(n) show the magnetic topology from simulation and reconstruction. In Figures 3(a)–(f), the colors represent the magnitude of the "out-of-plane" magnetic field, while the black lines represent the "in-plane" magnetic field. Apparently, in the X–Y and X–Z planes, the reconstruction and simulation are similar, and the spiral feature in the X–Z plane is well reproduced by the SOTE method. In the Y–Z plane, however, the reconstruction and simulation are somewhat different. Such a discrepancy may be attributed to the third-order variation of B_x along the Z direction, as shown in Figure 3(c). Figures 3(g)–(i) examine the error of magnitude, $\left|\bigtriangleup {\boldsymbol{B}}\right|=| {{\boldsymbol{B}}}_{\mathrm{SOTE}}-{{\boldsymbol{B}}}_{\mathrm{simu}}|$, during reconstruction, while Figures 3(j)–(l) examine the error of direction, $\theta =\mathrm{arcos}\left(-1\tfrac{| {{\boldsymbol{B}}}_{\mathrm{SOTE}}\,\cdot \,{{\boldsymbol{B}}}_{\mathrm{simu}}| }{| {{\boldsymbol{B}}}_{\mathrm{SOTE}}| \,\cdot \,| {{\boldsymbol{B}}}_{\mathrm{simu}}| }\right)$, during reconstruction. As can be seen, ∣△B∣ is small (<1) within the area about three times the size of the SC tetrahedron but becomes considerably large beyond this area; θ is small (<30°) in most of the area, even in the area where ∣△B∣ is large. This means that the SOTE method is able to reconstruct this structure accurately within an area about three times the size of the SC tetrahedron. In such an area, the 3D magnetic topologies from simulation (Figure 3(m)) and reconstruction (Figure 3(n)) are in good consistency.

Zoom In Zoom Out Reset image size

In this experiment, to match the MMS, we use a small tetrahedron (only six grids) to test the SOTE method. In the following experiments, we will enlarge the SC tetrahedron to comprehensively test the SOTE method.

An asymmetric X-line-like structure, which is the typical topology of magnetopause reconnection, is considered. The specific reconstruction procedure is same as that in Experiment 1, hence it will not be repeated. Figure 4 shows the reconstruction results in the same format as Figure 3. As can be seen, the reconstruction results (Figures 4(d)–(f)) and simulation results (Figures 4(a)–(c)) are quite similar in all the X–Y, X–Z, and Y–Z planes. In particular, the B_x reversal in the Y–Z plane (see the colors in Figure 4(c)) and the asymmetric antiparallel magnetic fields in the X–Y plane (see the black lines in Figure 4(a)) are both reproduced. ∣△B∣ and θ are both small within an area about three times the size of the SC tetrahedron (see Figures 4(g)–(l)). Particularly in the X–Y plane, which is the most important plane for studying magnetic reconnection, ∣△B∣ and θ are considerably small (see Figures 4(g) and (j)). The 3D magnetic topologies from the simulation (Figure 4(m)) and reconstruction (Figure 4(n)) agree well with each other, and both of them exhibit a clear X-line feature.

Zoom In Zoom Out Reset image size

A random magnetic structure (without any specific features) is used for the test to examine whether the SOTE method is applicable to a general case. Figure 5 shows the reconstruction results in the same format as Figure 3. As can be seen, the reconstruction results (Figures 5(d)–(f)) and simulation results (Figures 5(a)–(c)) are very similar in all the X–Y, X–Z, and Y–Z planes. ∣△B∣ and θ are both small within most of the reconstruction area (see Figures 5(g)–(l)). The 3D magnetic topologies from the simulation (Figure 5(m)) and reconstruction (Figure 5(n)) are very consistent, demonstrating that the SOTE method is able to reconstruct a random magnetic structure accurately and therefore can be applicable to a general case.

Zoom In Zoom Out Reset image size