To demonstrate how the aforementioned three steps are implemented, we first apply this new 3D ER model to an MHD system that only contains point-wise measurements of (B, J) fields together with MHD Eqs 1a, b as has been extensively investigated by the traditional least-squares method (e.g., Denton et al., 2020; Torbert et al., 2020). This reconstruction model can be considered an ER model for an ideal MHD system because the effect of resistivity (η) is not included. Extension to a more comprehensive nonlinear ER model that uses Eqs. 1–3 with point-wise measurements of (B, J) and (U, E) fields and incorporates the effects of plasma resistivity contained in Eqs. 2, 3 near the EDRs will be presented in our future investigations.

Here, we follow Torbert et al. (2020) and prescribe the form of the reconstructed field by expressing the time-independent magnetic field B as a quadratic function of the spatial coordinate r by a second-order Taylor expansion of a vector field B ( r ) ≈ B ( r 0 ) + [ D r B ( r 0 ) ] ( r − r 0 ) + 1 2 ( r − r 0 ) T [ D r 2 B ( r 0 ) ] ( r − r 0 ) , (4) where r 0 = ( 1 / 4 ) ∑ α = 1 : 4 r α is the barycenter of the tetrahedron defined at its four vertices by the locations of the four MMS spacecraft ( r α ( α = 1,2,3,4 ) ) , with D r B ( r 0 ) and D r 2 B ( r 0 ) being the first- and second-order derivatives of B at r 0 , respectively. The new ER model presented here is independent of the coordinate system though we have chosen to employ Geocentric Solar Ecliptic (GSE) coordinates. The terminology, notations and various manipulations of the tetrahedron geometry formed by a four-point satellite configuration have been discussed previously (e.g., Chanteur, 1998; Harvey, 1998; Robert et al., 1998; Dunlop et al., 2002). In addition to the barycenter, we may also define four face-centers ( r F α = ( 1 / 3 ) ∑ β ≠ α r β ) and six edge-centers ( r α β = ( r α + r β ) / 2 ) of the tetrahedron that can be easily calculated from the coordinates of the vertices. In practice, the coefficients of the derivatives in Eq. 4 will be determined by the reconstruction model based on the measurements. Hence, we may define the reconstruction model by rewriting Eq. 4 into the following explicit form for the ith component of the magnetic field B i ( r ) = B 0 i + ∑ j = 1 3 C 0 i , j Δ x j + 1 2 ∑ j , k = 1 3 D 0 i , j k Δ x j Δ x k ,     i = 1,2,3 , (5) where r = ( x 1 , x 2 , x 3 ) and Δ x j = x j − x 0 , j . The resulting smooth 3D magnetic field will be determined by thirty model parameters { B 0 i , C 0 i , j , D 0 i , j k } constrained by the MMS measurements. Given these model parameters, the spatial derivatives of the B field and the associated divergence ( ∇ ⋅ B ) and vorticity ( ∇ × B ) fields can be evaluated analytically and thus their valuations are available at any spatial point, though the measurements ( B ^ , J ^ ) are only available at the four vertices. Note that, physically, δ ( r ) ≡ ∇ ⋅ B ( r ) = ∑ i = 1 3 ∂ B i ( r ) / ∂ x i ≡ 0 for any value of r. Specifically, δ ( r 0 ) = 0 leads to ∑ i = 1 3 C 0 i , i = 0 and ∑ i = 1 3 ∑ j = 1 3 D 0 i , i j Δ x j = 0 for the quadratic expression of B given in Eq. 5. Likewise, the plasma current density J can also be evaluated analytically from the modeled B field by Eq. 1a. When using these analytic expressions, the field values and constraints evaluated at barycenter, four vertices and four face centers, such as J ( r F α ) = J F α and δ ( r 0 ) = δ ( r α ) = δ ( r F α ) = 0 , are of particular importance.

Given MMS measurements at the vertices ( r α ) of the magnetic field ( B ^ ) from the MMS Fluxgate Magnetometer (FGM) instruments (Russell et al., 2016) and particle current density ( J ^ ) from the Fast Plasma Investigation (FPI) sensors (Pollock et al., 2016), the model parameters { B 0 i , C 0 i , j , D 0 i , j k } in Eq. 5 can often be derived by minimizing a loss function as defined below. Here, the loss function characterizes 1) the model-measurement differences between the modeled ( B , J ) and measured ( B ^ , J ^ ) parameters and 2) the model departures corresponding to the violation of the MHD Eqs. 1a, b. For a linear system, the minimization procedure can also be reduced to solving a set of linear algebraic equations (e.g., Menke, 1989). Depending on whether the number of the adopted constraints is smaller than, equal to, or greater than the number of model parameters, the solution derived from the least-squares method could be under-, even-, or over-determined, respectively. Previous reconstruction models have focused on the even-determined solutions of a quadratic loss function (e.g., Dunlop et al., 1988; Torbert et al., 2020), for which the measurement errors were not explicitly considered. The new ER model presented here adopts a new method that derives the model parameters by directly minimizing a generalized loss function using a stochastic optimization method that contains random measurement errors and consists of a flexible number of constraints. As a result, the solution is always programmatically feasible regardless of whether the physical constraints defined by the MHD equations are linear or nonlinear and whether the system is under-, even-, or over-determined.

The generalized loss function (L) has the following form L = L O + w A ε A L A + w B ε B L B + w C ε C L C , (6) where the individual components of the loss function ( L O , L A , L B , L C ) are given by L O = 1 12 ∑ α = 1 4 ∑ i = 1 3 [ B i ( r α ) − B ^ α , i ] 2 , (7a) L A = 1 12 ∑ α = 1 4 ∑ i = 1 3 [ J i ( r α ) − J ^ α , i ] 2 , (7b) L B = 1 9 [ δ 2 ( r 0 ) + ∑ α = 1 4 δ 2 ( r α ) + ∑ α = 1 4 δ 2 ( r F α ) ] or   L B ∗ = 1 5 [ δ 2 ( r 0 ) + ∑ α = 1 4 δ 2 ( r α ) ] , and (7c) L C = 1 4 ∑ α = 1 4 [ μ 0 J ( r F α ) ⋅ ( Δ r β γ × Δ r β δ ) − ( B ¯ β γ ⋅ Δ r β γ + B ¯ γ δ ⋅ Δ r γ δ + B ¯ δ β ⋅ Δ r δ β ) ] 2 (7d) with Δ r β γ = ( r γ − r β ) being the edge vector connecting the vertices r β and r γ and B ¯ β γ = ( B ^ β + B ^ γ ) / 2 being the mean magnetic field on the edge Δ r β γ calculated by the measured B ^ field by using a linear approximation to obtain the field along an edge, between two spacecraft measurements. In Eq. 7, we use i to denote the dimensional index ranging 1-3 and use Greek letters to denote tetrahedron points or faces ranging 1–4. The components L O and L A each consist of twelve terms or twelve constraints and represent the differences of the modeled and measured fields at the vertices r α . Thus, L O and L A correspond to the model-measurement differences in the loss function. The component L B consists of nine physical constraints, which requires minimization of δ 2 ( r ) = ( ∇ ⋅ B ) 2 at nine particular spatial points (i.e., one barycenter r 0 , four vertices r α and four face centers r F α ). Because the measurements do not directly enter the expression, L B corresponds to the model departures or violations from the above MHD equations. L B can be replaced by L B ∗ , which neglects the face-center constraints. The component L C consists of four approximate physical constraints derived from the generic MHD equation obtained by applying Stokes’ Theorem to Ampere’s Law ( μ 0 ∬ S J ˜ ⋅ d S = ∮ C B ^ ⋅ d l ) on the four tetrahedron faces, which derives the current density components normal to the tetrahedron faces ( J ˜ ) by using the linear curlometer technique from the measured B ^ (Dunlop et al., 1988). A minimization between J ˜ and J projecting onto the normal directions of four tetrahedron faces yields L C . Thus, L C also possesses the nature of the model-measurement differences. Note that, as previously denoted, each face-center r F α in Eq. 7d is defined by other three vertices ( r β , r γ , r δ ) . Specification of the weighting factors ( w A , w B , w C ) in Eq. 6 determines the selection of the loss function components to be included in the reconstruction model. The scaling parameters ( ε A , ε B , ε C ) in Eq. 6 depend on the characteristic length scale of the tetrahedron and the dimensional factors of the loss functions. We will discuss the settings of these parameters in more detail below.

We first note the similarities and differences between L A and L C in Eqs. 6, 7. Both loss function components adopt the differences in current densities as constraints. L A is the difference between the modeled J and the measured particle current density J ^ at four vertices whereas L C is the difference between the modeled J components and the current density J ˜ components derived from the curlometer technique (i.e., using B ^ ) on the four tetrahedron face-centers. When both B ^ and J ^ are available and include direct measurement errors of the same order, L A is more accurate to be included in the generalized loss function L than L C because the J ˜ value used in L C contains additional errors due to the linear approximation assumed in the curlometer technique. On the other hand, if only B ^ , but not J ^ , is available (in which case L A will not be available) or if the errors in J ^ are far greater than those in B ^ , then L C is preferred to L A for inclusion in L. In Denton et al. (2020), J ˜ derived from the curlometer technique is used to modify the particle current density J ^ to produce a composite current density, which together with the measured B ^ is used to build the reconstruction model. Our approach of introducing different constraints L A and L C for different current densities J ^ (measured directly by FPI) and J ˜ (derived from the curlometer technique) evaluated at different spatial locations provides a clear physical significance and algorithmic flexibility.

In this new 3D ER model, the model parameters in Eq. 5 are solved by directly minimizing the loss function L defined in Eq. 6 using a stochastic optimization algorithm called the SPSA method (Spall, 1998a; Spall, 1998b; Spall, 2003) through an iterative procedure that also naturally incorporates the errors for the measured fields ( B ^ , J ^ ) . A comprehensive introduction to the algorithm with detailed procedures of implementation to the current problem is presented in Supplementary Appendix A. Note that the generalized loss function L defined by Eqs. 6, 7 is in a quadratic form with respect to the model parameters { B 0 i , C 0 i , j , D 0 i , j k } because the MHD Eqs. 1a, b are linear. Minimization of a quadratic loss function is equivalent to solving a set of linear algebraic equations for the model parameters (e.g., Menke, 1989; Axelsson, 1996). When the model parameters are obtained by directly minimizing the loss function L the corresponding MHD system could be either linear or nonlinear. The nonlinearity occurs in our new 3D ER model when Eqs. 2, 3 are also included as additional constraints. Nonlinear systems are not unusual in various empirical models. For example, in Roelof et al. (1993), the loss function for reconstructing global magnetospheric images based on the extreme ultraviolet (EUV) and energetic neutral atom (ENA) measurements is highly nonlinear, for which the loss function can only be directly minimized. Furthermore, when the loss function contains measurements, it also contains random measurement errors. The SPSA method effectively solves problems containing random errors by including the errors in the solutions. In addition, we will show later through examples that the SPSA method can solve slightly under-determined problems that could not be solved directly by the traditional least-squares approach.

In practice, the SPSA method solves for the model parameters that minimize the following dimensionless loss function (y) with a random perturbation that characterizes the measurement errors (Supplementary Eqs. A4a, b in Supplementary Appendix A) L θ = L / B 00   and (8a) y = L θ + ε σ , (8b) where L is given by Eq. 6, B 00 is the measured mean magnetic field (defined by Supplementary Eq. A1 in Supplementary Appendix A) used to normalize the general loss funciton L, ε σ = N ( 0 , σ 2 ) represents a random variable having a normal distribution with zero mean and σ 2 variance that characterizes the random measurement errors. The first-order SPSA algorithm is adopted to solve for the model parameters in this paper. The specifications of various model parameters including the weighting coefficients and scaling parameters in Eq. 6 and the algorithmic procedures of the recursive formulations are presented in Supplementary Appendix A. In Section 3, we will detail the application of this SPSA-based ER model to a specific EDR case using MMS measurements and discuss the relationship between the SPSA model variance σ 2 in Eq. 8b and the variances of the random errors of the measured B ^ and J ^ fields, σ B 2 and σ J 2 , respectively. Note that L θ in Eq. 8a is a deterministic variable whereas y in Eq. 8b is a randome variable. A stochastic optimation method such as SPSA algorithm optimizes loss functons associated with randome variables.

When all the loss function components in Eq. 6 ( L O , L A , L B , L C ) are included, then, the total number of constraints is thirty-seven (37). This number is greater than the number of model parameters (30) and the problem is significantly over-determined. Because L C adopts a linear approximation in the curlometer technique it is expected to introduce additional errors in the modeled fields near the reconnection regions where the field curvature is large. As a result, our default setting for the reconstruction model is to set w C = 0 , i.e., to not include L C in the generalized loss function L. This reduces the total number of constraints for the default setting to thirty-three (33) and thus renders the problem, i.e., solving thirty model parameters, only slightly over-determined. The model departures in the loss function component L B shown in Eq. 7c are the application of the MHD equation ∇ ⋅ B = 0 to nine particular points on the tetrahedron (the four vertices, the four face-centers, and the barycenter). When L B is replaced by L B ∗ that only applies ∇ ⋅ B = 0 to the barycenter plus four vertices, the total number of the constraints is reduced to twenty-nine (29) and the problem becomes slightly under-determined. Our numerical experiments show that the model parameters resulting from the SPSA method yield only slight and negligible (∼1–3%) differences when the problem is changed between slightly over-determined and slightly under-determined. On the other hand, the numerical solution to a set of under-determined linear algebraic equations no longer exists or cannot be calculated directly if the problem were solved by the traditional least-squares method (e.g., Menke, 1989).

To explain why using L B and L B ∗ does not lead to significantly different solutions, we first note that for an even-determined or an over-determined problem with a quadratic loss function, a unique solution can be derived either by directly solving an optimization problem or by solving a set of linear algebraic equations (e.g., Axelsson, 1996; Chong and Zak, 2001). It is also noted that for an over-determined problem, the inclusion of additional measurements or constraints may not change noticeably the existing solution if the newly added constraints are redundant (e.g., Menke, 1989). For an under-determined problem where the number of constraints is less than that of the model parameters, however, the set of linear algebraic equations becomes undetermined and one is no longer able to uniquely solve for the model parameters. Returning to the expressions of the loss functions in Eqs. 6–8, we note that the roles of model parameters and constraints (e.g., B vs. B ^ , or J vs. J ^ ) do not show preference to one or the other. A minimized or a least-squares solution is always formally available for given numbers of model parameters and constraints regardless of their relative magnitudes. Adding four constraints of ∇ ⋅ B = 0 to the four tetrahedron faces is expected to be largely redundant to the already existing constraints of ∇ ⋅ B = 0 at the barycenter and four vertices, thus leading to only slight modifications to the model parameters. Again, it is noted that unlike ( B ^ , J ^ ) that are only available on the four vertices, the analytic B-field as expressed by Eq. 5 and all its derived fields such as J and ∇ ⋅ B are available on any spatial point. Furthermore, in terms of the uniqueness of the solution, either the random noise term or the under-determined constraints in the loss function y in Eq. 8 could lead to the non-uniqueness of the solution. Note that the stochastic optimization algorithm minimizes the random variable y defined by Eq. 8b rather than the deterministic physical loss function L θ defined by Eq. 8a. We will discuss this issue in more detail in the next section. From the perspective of constraint redundancy, it is also noted that given the analytic expression in Eq. 5 for B, the relation ∇ ⋅ J = ∇ ⋅ ( ∇ × B / μ 0 ) ≡ 0 will be automatically satisfied regardless of what the model parameters are. As a result, one cannot introduce a constraint component for J similar to L B based on the redundant relation of ∇ ⋅ J = 0 .

Note that the measurement errors here include uncertainties in both the instrumentation and subsequent processing of the data. However, the random errors in Eq. 8 are associated with the unbiased instrument noise. Here, we estimate the errors by directly calculating the parameter variability included in the data series. In Figure 2, we show both the means ( B 0 , J 0 ) and the normalized standard deviations ( σ B , σ J ) of the magnitudes for the measured B ^ and J ^ fields. The averages are taken over the four spacecraft and over moving windows with widths of 7, 11, and 15 time steps, respectively. The calculated B 0 is approximately equal to the characteristic value of the temporal mean of the magnetic field B 00 defined in Supplementary Eq. A1. It is noted that the mean fields are not noticeably sensitive to the width of the moving window. This implies that the sampling rate of the measurements is high enough to resolve the temporal variability of the fields. It is also noted from Figure 2 that there is no systematic variation of σ B with respect to B 0 , whereas σ J is inversely proportional to J 0 . The weighting factor w A in Eq. 6 is proportional to the σ B 2 / σ J 2 parameter that can be calculated from the values shown in the figure. The weighting factors ( w B , w C ) are prescribed to ( 1,0 ) for the default setting of the reconstruction model. Note that setting w B = 1 here also means that we give no preference between the model-measurement differences L 0 and the model departures L B . To set the final model parameter σ used in Eq. 8, we note that σ B directly derived from the measured B ^ contains both the unbiased random errors required for the construction of ε σ in Eq. 8 and possibly also the biased errors associated with the parameter retrieval and data processing issues. In addition, a smaller σ in the random loss function y will yield a more numerically accurate solution, though its usefulness may be limited by the measurement errors; any numerical accuracy achieved that is higher than the measurement errors after setting σ → 0 does not contain additional information as the results are ultimately limited by the uncertainty in the measurements. As a result, our default setting for σ in the algorithm as shown in Eq. 8b takes a conservative value of σ = 0.1 σ ¯ B , where σ ¯ B is the time-averaged standard deviation of B ^ as shown in Figure 2A.

Given model parameters { B 0 i , C 0 i , j , D 0 i , j k } , a smooth 3D solution for ( B , J ) can be plotted to be visulized. But before addressing these visualizations, we begin our discussions here with evaluations of the quality factors associated with these results. In Figure 3, we show the relative differences of the fields ( B , J ) reconstructed at every time step based on the MMS-measured fields ( B ^ , J ^ ) shown in Figure 1. The indices ( γ B , γ J ) can be considered as the normalized loss function components ( L O , L A ) corresponding to the model-measurement differences, which can be used as a set of accuracy indicators of the reconstruction model and are defined as: γ B = ∑ α = 1 4 ∑ i = 1 3 [ B i ( r α ) − B ^ α , i ] 2 ∑ α = 1 4 ∑ i = 1 3 B ^ α , i 2   and (9a) γ J = ∑ α = 1 4 ∑ i = 1 3 [ J i ( r α ) − J ^ α , i ] 2 ∑ α = 1 4 ∑ i = 1 3 J ^ α , i 2 . (9b)

The results from a pair of reconstructions with σ = 0.1 σ ¯ B and σ = σ ¯ B , respectively, are presented in Figure 3. The default setting, which has a smaller measurement noise of σ = 0.1 σ ¯ B , yields a more accurate reconstruction field as characterized by smaller indices ( γ B , γ J ) . On the other hand, if the measurement noise in the loss function y amounts to σ ¯ B , such that σ ∼ σ ¯ B ∼ 0.1 as shown in Figure 2, then, a numerical solution of B with γ B < σ ¯ B can be considered to be an acceptable or valid solution. Our default setting of σ = 0.1 σ ¯ B leads to a numerical solution of B with γ B ≪ σ ¯ B , which can be considered an accurate solution. It should also be noted that because of the existence of measurement errors in B ^ (i.e., σ > 0 ), a deterministic and idealized solution with γ B ≡ 0 is considered to be as accurate as one with γ B < σ . Figure 3 shows that far greater errors exist in the modeled current density γ J than those in the magnetic field γ B . This is largely expected since the modeled current density J is a quantity derived from the prescribed B field and contains fewer free parameters and therefore is expected to lead to greater errors in J than in B. This is another reason for us to set σ so that it is much smaller than σ ¯ B in Eq. 8, which yields a solution also with an acceptable error in the reconstructed J field. Comparison between the two panels in Figure 3 shows that the magnitude of the errors in the reconstruction model is sensitive to the measurement errors. This feature can be further confirmed by examining γ J variation with the time-dependent measurement errors. We note from Figure 2 that σ J changes more significantly with time than σ B , which leads to a significant variation in the weighting factor w A . Figure 4 shows both w A ( = σ B 2 / σ J 2 ) and γ J for a standard setting of σ = 0.1 σ ¯ B . The figure shows a negative correlation between w A and γ J . Since errors in B ^ are nearly constant, Figure 4 shows a strong positive correlation between the errors in the measured J ^ and the modeled J. Overall, Figures 2–4 show that the accuracy of the reconstructed fields from the stochastic optimization algorithm is limited by the measurement uncertainties, with more accurate measurements unsurprisingly resulting in a more accurate solution for the reconstruction based upon those measurements.

Unlike reconstruction models based on the GSR technique, where the reconstructed fields are mainly derived by various physical relations, the new ER model presented here is mainly data-driven, directly fitting the modeled fields to the measured fields. Since the design of the general loss function highlighted in Eqs. 6, 7 also contains a component of the model departures characterizing a few physical constraints, the validity of those constraints can be used as a measure of the quality of the ER model in addition to the two indices ( γ B , γ J ) for measuring the accuracy of the solution. Here, one important constraint is the vanishing of the divergence of the magnetic field ( ∇ ⋅ B = 0 ), which is also used as a constraint of the loss function component L B in Eq. 7. Dunlop et al. (1988) introduced an index of the ratio of the divergence to the vorticity of the magnetic field as a quality indicator to measure the robustness of the reconstructed current density J field. In Figure 5, we show the following two quality indicators Q model and Q curl representing the ratio of the divergence to the vorticity of the magnetic field based on the reconstructed B field and the measured B ^ field by curlometer technique, respectively. These are defined as Q model = ∑ α = 1 4 ( ∇ ⋅ B ) α 2 ∑ α = 1 4 ∑ i = 1 3 ( ∇ × B | α , i ) 2   and (10a) Q curl = [ | ∇ ⋅ B ^ | | ∇ × B ^ | ] curlometer . (10b)

In the above, Q model is calculated by evaluating ∇ ⋅ B and ∇ × B analytically based on the modeled B field at the four vertices, whereas Q curl is calculated by evaluating the volume-averaged | ∇ ⋅ B ^ | and | ∇ × B ^ | based on the measured B ^ field following the schemes shown in Dunlop et al. (2002) and Middleton and Masson (2016). Note that Q curl has also been used as an objective index that measures the quality of the J fields reconstructed from the curlometer technique (Dunlop et al., 2002). On the other hand, the index Q model defined in Eq. 10a measures the quality or robustness of both B and J fields derived from the ER model. Figure 5 shows that typically Q model ≤ 0.01 for a standard setting of the model parameter σ = 0.1 σ ¯ B , whereas the typical values of Q curl are much greater, Q curl ≥ 0.1 . Since Q curl is calculated directly by a linear approximation from the measured B ^ field, it contains errors from both the linear approximation and measurement errors (Dunlop et al., 1988; Dunlop et al., 2002). Comparison between the two panels in Figure 5 also shows that the increase in the measurement errors by changing σ = 0.1 σ ¯ B to σ = σ ¯ B only increases the quality indicator Q model by a factor of ∼3 (note the changed scale on the ordinate between the two panels for just the “model” result, not the “curl” result). The relation of Q model ≪ Q curl is still valid for σ = σ ¯ B . As a result, we can conclude based on Figure 5 that the uncertainties in the reconstructed fields based on the MMS-measured B ^ near the EDR using the curlometer technique are mostly contributed by the linear approximation used in the technique. It should be pointed out that when an ER model is formulated and solved as an even-determined problem based on the traditional least-squares method, one may impose a condition of vanishing Q model everywhere ( Q model ≡ 0 ) . In this case, an alternative constraint corresponding to model departures, say, a vanishing variance of the modeled B field in a particular diretion M, i.e., ∂ 2 B / ∂ M 2 = 0 , needs to be introduced into the reconstruciton model (e.g., Torbert et al., 2020).

We now turn to the loss function component L C . The default setting of the weighting factors is ( w B , w C ) = ( 1,0 ) . This means that the constraint of the precise physical relation of ∇ ⋅ B = 0 is fully utilized whereas the constraint of matching the modeled current components to the ones derived from the curlometer technique on the tetrahedron faces is neglected. Again, setting w B = 1 here also means that we give no preference between two sets of constraints of the model-measurement differences and the model departures. Our analysis of the quality indicators ( Q model , Q curl ) derived from the runs without L C shown in Figure 5 can be considered as a rationale for the default setting of w C = 0 . It is noted that the curlometer technique developed in Dunlop et al. (1988) and Middleton and Masson (2016) applies a linear approximation to the entire volume of the tetrahedron, whereas in L C the linear approximation applies only to the four individual tetrahedron faces. Hence, it is worthwhile examining quantitatively the effect of the loss function component L C on the performance of the reconstruction model. In Figure 6, we show the indices ( r B , r J ) and the quality indicators ( Q model , Q curl ) for a sensitivity run of the reconstruction model with all parameters in default settings, except w C that is set to 1. Comparing Figure 6A with Figure 3A, we find that the inclusion of L C significantly reduces the accuracy of the reconstructed fields. This is expected because the inclusion of L C not only introduces a linear approximation in the calculation of the current density J from the magnetic field B, but also enhances the degree of over-determination of the model. Both of these are expected to increase the errors of a least-squares solution. Comparing Figure 6B with Figure 5A on the modeled quality indicators Q model derived from different runs underscores the same conclusion—i.e., that the inclusion of L C makes the model performance worse. However, Figure 6B also shows that Q model is still significantly smaller than Q curl (note the different scales), indicating that a linear approximation in a loss function component only partially affects the model performance. This sensitivity investigation of setting w C = 1 also demostrates the flexibility of the new ER model that directly minimizes the general loss function with its components being able to be included or excluded without changing the model framework.

At this stage, it is also interesting to examine a largely under-determined setting of excluding both L B and L C in the generalized loss function L by setting ( w B , w C ) = ( 0,0 ) in Eq. 6. There are only twenty-four (24) constraints in L O and L A , all given by the MMS measurements, whereas the reconstruction model contains thirty (30) model parameters that need to be determined. Hence, the problem is largely under-determined and the solution cannot be uniquely solved. For a stochastic optimization, such as the one based on the SPSA method, there is no fundamental difference in the non-uniqueness of the solution either due to the lack of constraints or due to random errors in the loss function. In other words, the unknown parameters for the reconstruction model can always be formally solved by minimizing the generalized loss function y in Eq. 8. Figure 7 shows the indices ( γ B , γ J ) and the quality indicators ( Q model , Q curl ) for a sensitivity run of the reconstruction model that sets ( w B , w C ) = ( 0,0 ) . We note from Figure 7A that the relative differences between the modeled and measured fields ( γ B , γ J ) are much less than those shown in Figures 3A, 6A. However, the quality indicator Q model shown in Figure 7B is much greater than those derived by any approach shown above including Q curl derived by the curlometer technique. Figure 7 shows that even though one can construct an empirical model that leads to a very good fit between the modeled and the measured fields at the prescribed spatial points, the fields may not necessarily satisfy some physical relations, such as ∇ ⋅ B = 0 . This is due to the following two facts: 1) the fields contain errors, either in the measured field or in the modeled field derived from the measurements and 2) the numerical evaluation of the physical relation based on the discrete measurements involves a small difference between two large quantities. The divergence of a vector field contains two components of variation corresponding to variations in the magnitude and direction of the vector. For a deformation vector field that is mainly confluent-diffluent—i.e., divergence is mainly caused by the change in direction—the calculation of the divergence of the vector field generally involves a small difference of two large quantities (e.g., Holton, 2004). In this case, small errors in the B field will be greatly amplified in calculating ∇ ⋅ B unless an additional constraint or assumption of ∇ ⋅ B = 0 , or | ∇ ⋅ B | being small, is explicitly included in the model or algorithm development. A similar assumption of “charge neutrality” in plasma physics is also used as an explicitly imposed constraint in developing various MHD models (e.g., Gurnett and Bhattacharjee, 2005).

The other more important implication of this test run for a largely under-determined setting with only 24 constraints for a 30-parameter reconstruction model is that the current ER model can be directly applied to reconstructing fields with a set of incomplete measurements. The stochastic optimization algorithms can solve for model parameters under the same algorithmic framework regardless whether the problem is over- or under-determined. For example, for a default setting of the current ER model with 33 constraints, the algorithm can be directly applied to an incomplete set of MMS measurements if the ( B ^ , J ^ ) measurements from one spacecraft are not available. Under such a circumstance, the same algorithm with 27 (= 33−6) constraints will produce a reconstruction field ( B , J ) that fits the measured ( B ^ , J ^ ) at three vertices having available measurements plus ∇ ⋅ B = 0 being satisfied at all four vertices, all within the measurement errors. The results shown in Figure 7 also suggest that the deterioration of the reconstructed fields due to lack of the needed constraints is gradual. On the other hand, the algorithm based on the traditional least-squares method that solves a set of linear algebraic equations (e.g., Torbert et al., 2020) becomes inapplicable once the problem changes from an even- to under-determined one.

We now present the reconstructed fields based on the MMS measurements shown in Figure 1. A reconstruction model can be developed in either an L-M-N coordinate system derived from the minimum variance analysis or in a fixed system, such as GSE that is used in the present reconstruction model. One purpose of adopting the L-M-N coordinate system to develop a reconstruction model is to take advantage of the ability to neglect changes in the minimum variance direction to convert a slightly under-determined problem into an even-determined one (Denton et al., 2020; Torbert et al., 2020). When the reconstructed field varies rapidly with time, the constructed L-M-N coordinate may also change accordingly. Under such a circumstance, the reconstructed fields at different time instances cannot be directly compared with each other. Our reconstruction model based on the SPSA stochastic optimization method can automatically accommodate an over-determined or under-determined setting of the model as discussed above. As a result, the fields reconstructed at different temporal instances but on the common, fixed GSE coordinate system can be directly compared.

We present the reconstructed fields in a local GSE coordinate X-Y-Z such that ( X , Y , Z ) = ( X ′ , Y ′ , Z ′ ) − ( X 0 , Y 0 , Z 0 ) , (11) where X ′ - Y ′ - Z ′ define the generic GSE coordinate system and ( X 0 , Y 0 , Z 0 ) = ( − 1.373 × 10 5 , 2.70 × 10 4 , 2.32 × 10 4 )  km is determined by the satellite constellation, which corresponds to the GSE coordinate of the mean barycenter averaged over the measurement time shown in Figure 1. In Figure 8, we show the reconstructed B fields projected into and its magnitude |B| ( = B 1 2 + B 2 2 + B 3 2 , in nT) evaluated on the X-Z plane of Y = 0 at six time instances of (a) t = 1.172  s, (b) t = 1.904  s, (c) t = 2.636  s, (d) t = 3.368  s, (e) t = 4.100  s, and (f) t = 4.833  s after 22:34 UT. The figure shows that both the magnetic configuration and the intensity of the magnetic field change noticeably with time. The reconnection region is characterized by a weak |B| and a reversal of the orientation (or a near anti-parallization) of the B vectors across the region. It is noted that a weak |B| also means a weak confinement to the motions of energetic electrons. This will lead to localized very fine-scale energy spectra and angular distributions that could be correlated with the remote magnetic topologies through the gyro-sounding process as revealed by the data from the Fly’s Eye Energetic Particle Spectrometer (FEEPS) onboard the MMS spacecrafts (Cohen et al., 2021; Turner et al., 2021). The development and evolvement of these two features can be easily identified in this figure. To provide a better view on the development of the reconnection region, we show in Figure 9A the superimposed B fields at two neighboring time instances of t = 2.636  s and t = 3.368  s on the same plot. The figure shows the development of an anti-parallel B field having a nearly opposite direction and an equal magnitude with a significantly weak B field sandwiched between the two regions at t = 3.368  s and especially in the region of X > 0 . Though the reconstruction in the present model is under the X-Y-Z coordinate whereas the reconstruction in Torbert et al. (2020) was presented in the L-M-N coordinate, the configuration of the reconstructed B-field shown in Figure 9A is qualitatively similar to that shown in Torbert et al. (2020). Figure 9B shows the corresponding cross tail current J on the Y-Z plane of X = 0 that shows a significant intensification in its magnitude due to the development of the reconnection event.