Reconstructing the Coronal Magnetic Field: The Role of Cross-field Currents in Solution Uniqueness

We present a new 3D magnetohydrostatic (MHS) direct elliptic solver for extrapolating the coronal magnetic field from photospheric boundary conditions in a manner consistent with an assumed plasma distribution. We use it to study the uniqueness of the reconstructed magnetic field as a function of how significant the plasma forcing is on the force balance of the magnetic field. To this end, we consider an analytic MHS model as ground truth. The model uses two free parameters to decompose the current into two parts: a magnetic-field-aligned component and a cross-field component. We perform a comprehensive study of the 2D parameter space to understand under what conditions the ground truth can be reproduced uniquely. We find that current oriented perpendicular to the magnetic field has a smaller solution space than the same amount of current oriented parallel to the magnetic field, and so MHS regimes with larger proportions of plasma-related forcing may be a promising avenue toward finding unique magnetic field reconstructions.


Introduction
Extrapolating the coronal magnetic field from photospheric boundary conditions (BCs) is a challenging problem that has led to a variety of methodologies. These are mostly based on the force-free-field approximation: the lowest-order approximation to the magnetohydrodynamic equations, ignoring all forces except self-balancing magnetic forces (note that the special case of zero currents is known as the potential field). The force-free field is known to be ill posed due to a lack of uniqueness in determining the atmospheric magnetic structure from the boundary data alone (DeRosa et al. 2009;Bineau 1972). This makes direct integration by numerical partial differential equation solvers problematic. Such solvers may or may not numerically converge because there is not a unique mapping between a given set of BCs and the atmospheric magnetic structure in the 3D volume. As a result, the numerical method breaks down.
Iterative techniques, such as the Grad-Rubin and Wheatland-Sturrock-Roumeliotis algorithms, may be used to converge to a solution, but have their own inconsistencies (Grad & Rubin 1958;Wheatland et al. 2000). The former does not converge when strong current densities are present and the latter requires overdetermined boundary information resulting in potential incompatibilities between BCs and governing physics (Inhester & Wiegelmann 2006). Moreover, the lack of uniqueness makes the use of such iterative solvers in determining the likely magnetic field from photospheric observation difficult to justify because there is no guarantee the solution captures the true coronal magnetic field. This is compounded by the possibility that small changes at the photospheric boundary may result in radically different solutions. Beyond this, if the magnetic field extrapolated from the BC is inconsistent with coronal observations, it is difficult to determine how to alter an iterative solver in a manner that could use these coronal observations to shrink the solution space.
Our overarching goal is to build an extrapolation technique that is flexible enough to address nonuniqueness issues by including coronal observations. These might include coronal polarimetric measurements or observations of the coronal plasma distribution. In this paper we focus on the latter, using a finite-difference-based magnetohydrostatic (MHS) computational solver to explore to what extent adding a known inhomogeneous plasma forcing to the magnetic field in the coronal atmosphere reduces the size of the possible solution space. By numerically demonstrating such a collapse of the solution space as a function of the degree of non-force-freeness, we justify future work in which both the lower magnetic boundary and prescription of plasma forcing are used as parameters in an inversion framework; see, e.g., Dalmasse et al. (2016Dalmasse et al. ( , 2019.

Background
The most rudimentary model of the solar coronal magnetic field ignores all electric currents above the photospheric boundary by solving the simple potential field model, ´= B 0. This has the benefit of being computationally efficient, as it is solvable by Fast Fourier Transforms (FFTs). Furthermore, it has been demonstrated that it is competitive with more complicated current-carrying models for reproducing observed magnetic structures such as coronal loops (DeRosa et al. 2009).
Next in complexity are force-free fields, so called because they assume a parallel current that exerts no net Lorentz force on the magnetic field. These may be linear, where the current is everywhere some scalar multiple α of the magnetic field. Or they may be nonlinear, where there is spatial variation of α. A direct vertical integration of nonlinear force-free fields (NLFFFs) from BC data is fundamentally ill posed (Bineau 1972). However, there are many iterative or minimization-based methods that allow the computation of NLFFFs (DeRosa et al. 2009).
Observations of the solar corona argue that it is not generally appropriate to discount all plasma forces in the structured coronal atmosphere (Suess et al. 1996;Gibson et al. 1999). Some solvers account for plasma-derived forces and so can accommodate a non-self-balancing field (e.g., Zhu & Wiegelmann 2018). The solver constructed for the parameter study used in our analysis is of this MHS type.
Finally we note that there are methods which model a time dependence in some magnetohydrodynamical system. These models are very robust and thus invaluable as scientific investigative tools (e.g., Rempel 2017). That said, they can be computationally impractical for quick, iterative calculation.

Analytic Ground-truth Model
In a force-free field, α controls the magnitude of currents parallel to the magnitude field. By introducing an analogous scalar a that controls the magnitude of the current perpendicular to the magnetic field, it is possible to allow control over the forced component of the field in the MHS solution. Larger a and α both drive the field farther from the potential state, representing the presence of more currents and free magnetic energy in the field.
We take the analytic solution from Low (1992), which provides a solution magnetic field B alongside pressure P and density ρ for the following subspace of MHS: where a, α, and κ are scalar parameters and p 0 is the hydrostatic background pressure. Note that α and a behave as described above, where α controls the magnitude of the magnetic field with parallel currents and a impacts the crosscurrent magnetic field. It is important to note that α must be bounded in magnitude by unity (α<1) to maintain a physical system with finite magnetic field at large distances. κ is related to the plasma scale-height, and we follow Low (1992) by setting it to a constant value of 1.7 in this exploratory study. Solutions to the linear subspace of MHS solutions induced by the above restriction on B can be described in terms of Bessel functions of the first kind, denoted J, in cylindrical coordinates as vertical coordinate variable, which allows the wave-mode solution to accommodate gravity, and λ and s are Fourier parameters with restricted domains. For convenience, a detailed description of the model as we employ it is in Appendix A. Importantly, this analytic formulation allows free control of the magnitude of the parallel-and cross-field currents (via α and a, respectively). The parallel-field currents govern how much the magnetic field as a whole tends toward the force-free regime, while the cross-field currents require balance from plasma forces.
By increasing a, we thus effectively increase the magnitude of the plasma forces and are able to explore the impact that such forces may have on the solution space ( Figure 1). Examination of Figures 1(a), (c), and (e) reveals that increasing the parameter α effectively introduces a clockwise rotation of field-line footpoints around the central bipole. Comparing these to Figures 1(b), (d), and (f) shows that increasing a inflates the field, adding additional twist about the polarity inversion lines. Figure 1 illustrates that these effects together may create flux rope structures parallel to the photospheric surface.

The Solver
To inspect the solution space for extrapolations of the BCs of this analytic model, we present an MHS numerical solver to generate solutions from BC and known atmospheric plasma data. The solver is presented as a finite-differences-based method, though generalizable to any sparse matrix differentiation technique. It includes a 3D analytically constructed Jacobian and corresponding least-squares solution. The numerical solver is a direct integration approach to obtaining an MHS solution matching the lower boundary data. As such, it enables the inspection of the parameter space of the analytic model and to what extent the regimes that generate nearly unique magnetic fields rely continuously on the boundary and plasma forcing.
Consider the MHS equations where P and ρ are predetermined in the volume, We assume in Equation (3) that the determination is selfconsistent between P and ρ, or in the case of this work, that m p P=ρk B T for an isothermal temperature T. This is thus equivalent to the following forced nonlinear equation with forcing field F (as defined in Equation (A5)) and mixed BCs: In this analysis, we will use both the BCs and a prescribed 3D forcing F associated with the choice of parameter a. To solve Equation (4) numerically, we first construct an initial guess for the magnetic field based on BC data, for example, the potential field which can be matched up to a given = B z z 0 | with an FFT solution. We choose this because it is numerically efficient, and the literature suggests the potential field may be close to the physically observable system (DeRosa et al. 2009). In these computations, we use N=125,000 (50 3 ) computational points arranged on an equispaced Cartesian lattice in x and y, but with exponentially decaying density in the vertical z direction.
Next we construct a measure of how non-force free the system is by constructing the residual vector for a Newton iteration step n: where  is the Jacobian of the system and ε is a small scalar, which we set to ε=0.1 for the purposes of this study. In this way, we find an update that reduces the MHS imbalance whose divergence is similar to the divergence of the previous vector field. This allows us to ensure that the divergence of this numerical field stays small, without building the divergence equations into the Jacobian itself. We repeat the process until the magnetic field stops changing from iteration to iteration, that is, - || | | is everywhere sufficiently small, e.g., less than 10 −4 , thus demonstrating convergence. Typically for Newton's method this convergence is quick, occurring in less than 10 iterations for reasonable tolerances. The field is modulated by a low-magnitude vector potential to remove any remaining divergence from the field. It is more physically important that the end vector field product be divergence free than that it be perfectly self-balanced, and so the final field is shifted as Finally, we note that in some cases a significant net force on the numeric magnetic field may remain, implying that MHS force balance has not been achieved and the numerical method has broken down. We refer to this as the residual, plotted in Figures 5 and 6, and further described in Appendix B.
We note that the numerically generated magnetic field lines accurately capture many features of the analytic model not present in the initial potential field guess, including the twisted, low-lying inner flux rope tightly bound around the polarity inversion line (compare Figure 1(d) to 2(d)). However, for values of large α, the numerical solution is significantly different from the analytic solution. In particular, the clockwise rotation of field-line footpoints about the central bipole that was noted above is not captured very well by the numerical solutions, although the twist about the polarity inversion lines still appears to be. This effect is most apparent in comparing || || . In Figure 3 we plot this quantity at several horizontal slices for both the analytic and numeric solutions associated with α=0.4, a=4. We see that even though the field lines have some discrepancies, this quantity is well preserved between them.

Parameter Exploration
As discussed above, the fully NLFFF case is not well posed for choices of α greater than some small threshold (Bineau 1972). That is to say, when the forcing field F (see Equation (4)) is zero, we expect different methods to produce different coronal magnetic fields from the same photospheric boundary.
We wish to inspect whether adding a cross-field current reduces the size of the solution space. Adding a purely fieldparallel current to a potential field base, we expect the difference between two valid solutions originating from different methodologies to increase with α, denoting the expanding size of the solution space. We refer to this as the discrepancy between analytic and numerical solutions and interpret it as a measure of the nonuniqueness of the numerical solution.
To this end, we analyze the performance of this numerical scheme for a suite of BCs and plasma-forcing functions corresponding to various parameter choices of the analytic model. The potential field is our initial guess, and we vary α: [0, 1] and a: [0,4], consistent with the range adopted by Low (1992). See Figure 1 for examples of analytic solutions in this range and Figure 2 for the equivalent numerical solutions. In Figure 3. On the left is the spatially varying current helicity of the analytic solution from Low (1992), and on the right the current helicity of the numeric solution with the same BC and plasma forcing. The helicity is predominantly focused along the bundles of field lines following the polarity inversion lines. Figure 4, we plot the discrepancy between numeric and analytic magnetic fields spatially.
The numerical representations obtained are significantly different from the associated analytic solutions in their manifestation of the effects of increasing α. Because they are valid solutions to the BCs and forcing function, this implies that there is nonuniqueness intrinsic to the effect of α for the range of parameters studied here. In contrast, the structures generated by increasing a are well preserved, implying that the solution space is more uniquely determined with regard to that parameter.
It is of particular interest to note that for much of the (α, a) parameter space, we maintain a reasonably small solution space, such that the discrepancy between the numeric field and the analytic one is less than 5%-10%. However, the α=0.8 case has large discrepancies throughout the coronal volume. At first glance it may appear that this is simply correlation between nonpotentiality (i.e., total current magnitude) and lack of uniqueness. However, the total current is much stronger in the (α=0.4, a=4) case than in the (α=0.8, a=0) case (as will be discussed later in Figure 5). In particular, in the analytic solutions the current J is, at its strongest point, This suggests that, far from increasing the size of the discrepancy, a larger current can in fact reduce the size of the solution space provided it has a component perpendicular to the magnetic field. Figure 5 shows the maximum value of the discrepancy, in particular the vector magnetic difference between analytic and numerical solutions within volumes such as those shown in Figure 4, for the full α-a parameter space explored. We see that for a=0, increasing the field-parallel current through α drastically increases the measured discrepancy between the analytic field and the one generated by numerical methods as expected. However, increasing the perpendicular current through the parameter a has a much smaller effect on the discrepancy. Moreover, Figure 5(b) illustrates that for a fixed current strength (shown through the color contouring), the more perpendicular that current is to the magnetic field, the smaller the discrepancy. In other words, the nonunique solution space for a given BC may steadily collapse as the current orientation rotates away from parallel to the magnetic field.
It is our goal now to qualify how confidently we can assert that the discrepancies observed are a result of the properties of MHS itself and not a quirk of the numerical method. To that end, we use the same routine with a different initial guess for the Newton's Method MHS solver. We wish to see if the vertical contours moving outward from the initial guess in Figure 5 are a result of the numerics or more fundamental to the equation.
In Figure 6, we use the analytic solution from Low (1992) with the choice of parameters α=0.4, a=2 as the initial guess for the Newton iterative solver as opposed to the FFT potential field solution. This perturbation of the numerics is found to have little impact on the overall structure of the contours described above; they are still largely vertical, moving out from the island of perfect agreement where the initial guess is located in parameter space. This suggests that this behavior is Figure 5. Contours of maximum discrepancy between the analytic-and numerical-solution magnetic fields with matching boundary conditions, normalized to the greatest magnetic field strength. These computations were completed with a sparse N=27,000 (30 3 ) grid. In panel (a), color denotes the maximum residual or net force on the numeric magnetic field (which should be 0 for MHS balance). That is, the darkest color would denote that 10% or more of the magnetic forces are unbalanced by equivalent plasma forcing. A large net force implies that the numerical method is beginning to break down. In panel (b), color denotes the strength of the current in each simulation. Figure 6. The same domain as in Figure 5, where a different initial guess is used in the Newton solver. The contours and colors are as in Figure 5(a). fundamental to the degeneracy inherent in the equations, not simply an artifact of numerical convergence.

Conclusions
The analytic model and the numerical solver have provided two MHS solutions that share a lower boundary and coronal plasma-forcing profile. Where these two representatives of the population of solutions are close, we may infer that the solution space is small. Where they are far apart, we may infer that there is a large, or highly degenerate, solution space.
We have demonstrated that the solution space for extremely low-current fields is very small. As the parallel current increases, the size of the solution space increases alongside it. However, adding cross-field current has little impact on the solution space. This suggests that rotating current to be more cross-field can have a beneficial impact on the uniqueness of MHS solutions, and so coronal atmospheres with significant cross-field currents may benefit from an MHS direct solver, rather than a back-and-forth iterative NLFFF technique.
In particular, observations of the plasma distribution may be used to remove degeneracy from numerical solutions based on the boundary magnetic field, yielding a reasonably unique magnetic field output. Finally, the incorporation of coronal polarimetric data , as may be obtained by telescopes such as the Coronal Multichannel Polarimeter (Tomczyk et al. 2008)

Appendix A Analytic Model Details
Following Low (1992), for coefficients λ, α, and a, we let κ=1.7, l k = q a 2 , and l a k = s 2 2 2 . We will construct three angularly symmetric solutions and add them together centered at different spatial locations (x n , y n ). Define the following parameters according to the table