Testing magnetohydrostatic extrapolation with radiative MHD simulation of a solar ﬂare

Context. On the sun, the magnetic ﬁeld vector is measured routinely solely in the photosphere. By using these photospheric measure- ments as a boundary condition, we developed magnetohydrostatic (MHS) extrapolation to model the solar atmosphere. The model makes assumptions about the relative importance of magnetic and non-magnetic forces. While the solar corona is force-free, this is not the case with regard to the photosphere and chromosphere. Aims. The model has previously been tested with an exact equilibria. Here we present a more challenging and more realistic test of our model with the radiative magnetohydrodynamic simulation of a solar ﬂare. Methods. By using the optimization method, the MHS model computes the magnetic ﬁeld, plasma pressure and density self- consistently. The nonlinear force-free ﬁeld (NLFFF) and gravity-stratiﬁed atmosphere along the ﬁeld line are assumed as the initial conditions for optimization. Results. Compared with the NLFFF, the MHS model provides an improved magnetic ﬁeld not only in magnitude and direction, but also in magnetic connectivity. In addition, the MHS model is capable of recovering the main structure of plasma in the photosphere and chromosphere.


Introduction
The nonlinear force-free field (NLFFF; Wiegelmann & Sakurai 2012;Guo et al. 2017) is considered a state-of-the-art model of the solar coronal magnetic field thanks to the low plasma β (Gary 2001) in the corona. This, however, is not the case in the photosphere and chromosphere, where plasma β is close to or even larger than unity. By calculating the net Lorentz force of AR7216, Metcalf et al. (1995) find that the magnetic field is not force-free in the photosphere. Zhu et al. (2016Zhu et al. ( , 2017 derive the magnetic field configuration of active regions using the magnetohydrodynamic (MHD) relaxation approach. They find the forcefree assumption fails beneath a height of 1.8 Mm. In order to study the magnetic field in the non-force-free region, a magnetohydrostatic (MHS) extrapolation should be used.
The MHS extrapolation is introduced on the basis of the MHS equilibria assumption, which is supposed to be a better approximation in the low layers of the Sun than the force-free assumption. The governing equations of the MHS equilibria can be written as where B, p and ρ are the magnetic field, plasma pressure and plasma density, respectively. We note that the equations above have been normalized using the following constants: ρ 0 = 2.7 × 10 −1 g cm −3 (density), T 0 = 6 × 10 3 K (temperature), g = 2.7 × 10 4 cm s −2 (gravitational acceleration), L 0 = RT 0 µg = 1.8×10 7 cm (length), p 0 = ρ 0 RT 0 µ = 1.3×10 5 dyn cm −2 (plasma pressure), and B 0 = 4πp 0 = 1.3 × 10 3 G (magnetic field). Low (1985Low ( , 1991Low ( , 1992 and Neukirch & Rastätter (1999) found the MHS equations can be linearized when the current consists of a linear component parallel to the magnetic field and another component perpendicular to gravity. Using this so-called linear MHS model, Aulanier et al. (1998Aulanier et al. ( , 1999 derived the magnetic field and plasma distribution of AR7722 and AR7986. More recently, Wiegelmann et al. (2015Wiegelmann et al. ( , 2017 use this model to extrapolate the MHS equilibria by SUNRISE/IMaX magnetogram. The unprecedented resolution (40 km pixel −1 ) of IMaX enabled the model to resolve the thin non-force-free layer with several tens of grids. However, the linear MHS model excludes the high concentration of electric currents and Lorentz forces. This is somewhat similar to the well-known limitations of linear force-free field in modeling the corona.
A few methods have been developed to solve the MHS equations in the general case numerically. Hu & Dasgupta (2006 propose an approach for deriving the non-force-free field by superposing one potential field and two linear forcefree fields. The Grad-Rubin iteration procedure is extended by Gilchrist & Wheatland (2013), Gilchrist et al. (2016) to compute the MHS equilibria. The MHD relaxation method, which applies an "evolution technique", is capable of yielding the MHS solution (McClymont & Mikic 1994;Jiao et al. 1997;Zhu et al. 2013). Wiegelmann & Inhester (2003), Wiegelmann & Neukirch (2006) and Wiegelmann et al. (2007) use the optimization method to treat the MHS equations without gravity in different coordinates systems. We extend the optimization method to model the system with a gravity force in Zhu & Wiegelmann (2018, hereafter Paper I). In that work, we test our model using a perfect MHS equilibria (Low 1985(Low , 1991. However, the actual Sun is very dynamic and a lot more complex. Here we aim to investigate whether the MHS extrapolation functions in a realistic situation. In this paper, we apply our numerical code to reconstruct a snapshot of the solar flare simulation. The simulation includes 3D radiative transfer in the convection zone and photosphere, optically thin radiation and field-aligned heat conduction in the corona, etc. These features set the simulation into a realistic environment. The organization of the paper is as follows. In Sect. 2, we briefly describe the simulation and assess the physical state of the reference snapshot. In Sect. 3, we introduce the optimization method for the MHS model and the numerical setup to extrapolate the reference snapshot. In Sect. 4, we evaluate the results through several metrics. In Sect. 5, we discuss a few factors that can impact the ability of our model to reconstruct the magnetic field and plasma. Concluding remarks are presented in Sect. 6. Cheung et al. (2019) carry out a radiative MHD simulation of a solar flare. They use the Max-Planck-Institute for Aeronomy/University of Chicago Radiation Magneto-hydrodynamics (MURaM) code, which allows for simulations spanning from the upper convection zone into the solar corona (Vögler et al. 2005;Rempel 2017). The initial setup consists of bipolar sunspots. After the thermally relaxed equilibra are reached, they impose the emergence of another twisted bipolar flux system at the bottom boundary. The accumulation of the magnetic energy finally leads to an abrupt relaxation which powers coronal mass ejection and a flare.

The active region in the radiative MHD simulation
Among all the snapshots that are available for download 1 , we choose the last one (at about 8 minutes after the flare) as the reference model for our test. The original data spans ±49.152 and ±24, 576 Mm in the x-and y-axes, respectively. In the z-axis, the data spans from 7.5 Mm beneath the photosphere to 41.6 Mm above. The grid spacing is 192 (64) km in the horizontal (vertical) direction. In this study, we focus on the MHS equilibrium in the lower atmosphere. For this purpose, the data from photosphere to 8.192 Mm above are extracted. The photosphere is cut at the average geometrical height corresponding to optical depth unity. This is done by using "hgcr_hslice.pro" with input parameter "k = 116" which is suggested in "readme.txt". Both "hgcr_hslice.pro" and "readme.txt" can be downloaded from the above data link. Photospheric magnetogram and the magnetic field lines in the sub-volume are illustrated in Fig. 1.
The system continues to develop the entire time over the course of the simulation. However, Fig. 2a shows, at the moment we are considering, that the time-varying term is smaller than other terms in the momentum equation. The reference snapshot is not strictly static (see panel b). Part of the velocity comes from the violence of the flare eruption. To see how the moving fluid affects the magnetic field and plasma, we look at the momentum equation in a steady state.
where ρ, p, u and B are mass density, pressure, velocity and magnetic field, respectively. The inertial term (P v = ρv 2 ), plasma pressure (P g = p) and magnetic pressure(P b = B 2 /2) in the above equation compensate for one another. Figure 3 shows the influence of these three terms. We see from panel a that: (1) plasma pressure affects the magnetic field significantly near the 1 http://purl.stanford.edu/dv883vb9686 photosphere; (2) the magnetic field dominates the three terms above 300 km. Panel b shows that the strength of inertial term is stronger than that of plasma pressure in the layer 1 Mm < z < 2 Mm. We deduce that in the extrapolation: (1) plasma pressure should be taken into account near the photosphere; (2) excluding the inertial term has limited effects on the pattern of the magnetic field; (3) excluding the inertial term has significant impact on the plasma distribution above 1 Mm. We conclude that this snapshot is suitable, albeit imperfect, to test the MHS extrapolation.

Method
The optimization method used to solve the MHS equations in this study is described in detail in Paper I. Here we give a short outline and its recent improvements. The method involves the A162, page 2 of 8 minimization of the functional with where ω a and ω b the weighting functions and ν the Lagrangian multiplier. W is a diagonal matrix to incorporate the measurement error of the magnetic field and gas pressure. [B opt p opt ρ opt ] and [B obs p obs ρ obs ] are the optimized and observed quantities in the photosphere respectively. The weighting functions and slow boundary injection are found to be capable of improving the results (Wiegelmann 2004;). Here, for simplicity, we set the weighting function to unity over the whole region and use a fixed lower boundary. We note that in Eqs. (5) and (6), the denominator (B 2 + p) is different than that of the previous code (B 2 ). We make the change in order to avoid nonphysically large contribution to the functional L from the weak magnetic field region. Then, the solution of the MHS equations can be optimized by with p = Q 2 and ρ = R 2 which ensure the positive p and ρ.
The gradient descent method is used to find the minimum: where δL B , δL Q and δL R (see Appendix A) are the functional derivatives respect to function B, Q and R, respectively. µ B , µ Q and µ R control the step length. In Paper I, pressure distribution on the photospheric pressure was derived by solving the following Poisson's equation where f ph is the 2D horizontal Lorentz force on the photosphere, In this study, however, the Poisson method does not work effectively because of the dynamics in the simulation. Instead, a simple vertical magnetic flux tube is assumed to derive the pressure distribution from the force balance condition. That is: where p quiet is the plasma pressure in the quiet Sun where the magnetic field is extremely weak. We set p = 10 −3 p quiet when negative pressure appeared from Eq. (12).

Apply to a MURaM snapshot
Based on a photospheric vector magnetogram from the MURaM simulation, we extrapolate the magnetic field, plasma pressure and density on 512 × 256 × 128 grid points. The grid spacing (192/64 km in horizontal/vertical direction) is the same as the one in the simulation. The procedure shown below to perform the MHS extrapolation is slightly different from the one presented in Paper I: 1. Calculate a NLFFF model (Wiegelmann 2004;. 2. Distribute the pressure in the photosphere by using p + B 2 z /2 = p quiet . Calculate the pressure along the magnetic field line with the gravity-stratified assumption (see 1D temperature profile in Fig. 4). Calculate the density in the computational box by using ideal gas with a 1D temperature model.
3. Iterate for B, Q and R by Eqs. (8)-(10). This step is repeated until L reaches its minimum.

Plasma solution
The MHS model computes the plasma in the computational box, which is an important advantage over the NLFFF model. Figures 5 and 6 compare the pressure and density of the two models at different levels. At the bottom boundary, the plasma of the MHS model is very similar to that of the reference model, which denotes that p + B 2 z /2 = p quiet does indeed make sense. Above the bottom boundary, the MHS model is able to recover the main structures of the plasma in the photosphere and lower chromosphere (below 1 Mm). For example: the depletion of pressure and density in the strong field region, the spiral structure (indicated by arrows in Figs. 5 and 6) in the new emerging region, and the "tentacle" around the spots as well. From another perspective, Fig. 7 illustrates pressure and density in a vertical plane. We clearly see the main structures below 1 Mm are recovered well. Above 1 Mm, however, a large deviation appears. As we did not use the temperature data in the modeling, the plasma cannot be determined precisely. Another reason for the inaccurate plasma result comes from the dynamic nature of the atmosphere. Figure 8 shows the differences in the flux tube, which is widely studied in solar activity events. The field lines of different models start from the same seeds where the vertical electric current is strong. The comparison signifies that the MHS model is superior to the NLFFF model in reconstructing long twisted field lines. This leads not only to a different amount of twist, but also to the connectivity of the magnetic field. As illustrated in Fig. 9 the left footpoint is connected to another footpoint far away in the twin emerging spot in both the reference and MHS model. In the NLFFF, however, it connects to the closer side of the untwisted spot. The NLFFF line simply fails to bend up before it touches the photosphere, which results in a large difference in connectivity. It is worth noting that the vector magnetogram used in NLFFF extrapolation is preprocessed to smooth and fulfill the force-free conditions. Without preprocessing the photospheric boundary, the data are inconsistent with the forcefree assumption and consequent force-free computations do not converge (which is not shown here). There have been some studies reporting that the preprocessing improves the extrapolation considerably Metcalf et al. 2008).

3D magnetic structure in the models
To quantify how well each is extrapolation performed, we use 6 metrics as noted in Schrijver et al. (2006) and Barnes et al. (2006) to compare vector field b (our results) and B (the reference MHD model): where C vec , C CS , E N , E M and E e are vector correlation, Cauchy-Schwarz inequality, normalized vector error, mean vector error A162, page 4 of 8 and magnetic energy metric, respectively. N is the number of grid points in the computation box. In addition, the field line divergence (FLD) metric measures how well the magnetic connectivity is recovered. For both the simulation and extrapolations, a score can be given to any point, where the field line is closed, by the distance between the endpoints divided by the length of the field line. A single score is assigned by the fraction of the area at the lower boundary that has an FLD below 20% (10%, as adopted in Barnes et al. 2006 andMetcalf et al. 2008). The domain of comparison (within a grid size of 210 × 180×64) has the same field-of-view (FOV) as Fig. 1b, extending upward to the height of 4.1 Mm.
A quantitative evaluation of the two models is given in Table 1. The MHS model scores better than NLFFF in all metrics. It is not clear, based on the minute difference between C vec and C CS , MHS extrapolation is a significant improvement over the NLFFF extrapolation. However, for the metric 1 − E m , the MHS extrapolation scores 0.65, which is an improvement over the NLFFF's score of 0.57. It should be noted that C vec and C CS are sensitive to differences in angle between the two vectors that are being compared, whereas E n and E m are sensitive to both angle and norm differences. The FLD is the most rigorous metric among the six we used. It also indicates that the MHS model performs better than the NLFFF model.
We also present a comparison of the first four metrics in every plane along the height (see Fig. 10). Again, the MHS extrapolation scores better. It is worth noting that the weak field region can contribute significantly to the metrics of C CS and E M (see Eqs. (14) and (16)). That is the reason why the two metrics drop dramatically in the layers near the photosphere. The scores of two different extrapolations get closer as height increases, which indicates that the influence of the plasma effect decreases as height increases.

Magnetic twist and quasi-separatrix layers (QSLs) in the models
Finally, we compare the magnetic twist numbers T w , T g and squashing factor Q s for different models. T w and T g are defined by Berger & Prior (2006) as the twist number of one magnetic field line about the neighboring line and the axis, respectively. They are different in most situations. The self-helicity which could be measured by local twist T w is negligible for a large number of flux tube (Demoulin et al. 2006). In this study, T g is applicable to evaluate the twist since only one flux rope with finite size exists. T w and T g are given as: where T(s) is the unit tangent vector to the axis curve, V(s) is a unit vector normal to T(s) and points to the secondary curve, J is the parallel component of the electric current. The integration is carried out along the specific field line. QSLs are the generalized topological structures (Demoulin et al. 1996) which are the regions with high squashing factor Q s . Q s is defined by mapping the field line (Titov et al. 2002): where (x 1 , y 1 ) and (x 2 , y 2 ) are the two footpoints of one field line. Both indexes are important in studying the 3D magnetic structure. We used the code developed by Liu et al. (2016) to calculate T w and Q s . To compute T g , we first have to locate the axis. We assume as the axis any one of a number of magnetic field lines that penetrate the vertical square (4×4 Mm, thick black line from the top view shown in the left panel of Fig. 11) and a score can be assigned to this line based on the average T g of the surrounding lines. Then we define the axis as the field line whose score is the greatest (see the axis and sample surrounding field lines in Fig. 11). A similar approach is used in Guo et al. (2010). For each model, the axis is determined by the same way. As a result, the average T g about the axes of the reference, MHS and NLFFF model is 0.95, 0.94 and 0.88, respectively. The twists of the field lines close to the axis are sufficiently reconstructed by the MHS model. The twist number T w is illustrated in Fig. 12a-f. The white ellipses represent the central area that is mainly recovered, whereas the black ellipses show the mismatched areas. We can see a similar mismatch in the two extrapolations since the MHS model uses the NLFFF solution as the initial condition. If we look to some T w structures with small size, the MHS model includes more of them that exist in the reference model. That is because the MHS model uses the unpreprocessed magnetogram while the NLFFF model uses the preprocessed one. Without preprocessing, NLFFF does not converge because the data are inconsistent with the force-free assumption. Preprocessing makes the data force-free consistent, but it also naturally removes structures related to finite Lorentz-forces. The QSLs represent the layers which separate magnetic systems with different topologies. We can see the QSLs, which are indicated by the white arrow (Fig. 12j), mainly divide the central area into two topological parts. However, there are additional QSLs in the two extrapolations which divide the same area into three parts. So we should be very careful to use the QSLs result from extrapolation.

Discussion
The inability to generate an MHS model that sufficiently recovers the reference model is rather disappointing. This has led us to examine the new extrapolation in order to distinguish certain factors that impact out ability to produce a robust model.  = 5 a represents case with potential field as the initial condition.

Choice of free parameters and initial conditions
The three parameters µ B , µ p and µ ρ control the step length when the Eq. (4) is minimizing. µ B is computed at every step by the line search method. We assume µ p = µ B because Q has the same dimension as B. We further assume µ ρ = µ B . Table 2 shows the case study with varied between 0.1 and 10. We find the functional L decreases with increasing . However, the most accurate magnetic field generated by the code when = 5. Figure 13 shows how the free parameter affect the pressure and density result. The density is much more affected than the pressure by .
A larger leads to a stronger contrast of the density because of the higher weight introduced by . From the above case studies, we set = 5 as the optimized value for the extrapolation. The choice of the initial condition of magnetic field has significant influence on the results (Wiegelmann 2004;Schrijver et al. 2006;Zhu & Wiegelmann 2018). We ran an additional test of = 5 with potential field as the initial condition. The test is named as 5 a in Table 2 and Fig. 13. It shows that the MHS extrapolation with the initial potential field has lower scores for the magnetic field evaluation and less accurate plasma distribution. Table 1 also compares the magnetic field of reference model with the result for case "MHS (a) " for which the pressure and density in the bottom boundary is provided. We observe a limited increase in the metrics for the magnetic field comparison. Figures 14 and 15 show the plasma distribution at the extremely low height (≤0.32 Mm) benefits from the application of the accurate bottom boundary condition. Above 1 Mm, however, the plasma result shows no improvement. The increasing dynamics at the higher atmosphere make the result resistant to the choice of the plasma boundary condition at the bottom. Nevertheless, the simplest boundary condition with the plasma uniformly distributed is not able to yield a satisfactory result. So, the plasma boundary described in Sect. 3.1 remains an adequate choice.

Conclusions
In this work, we applied the MHS extrapolation code to model an active region that serves as a snapshot of the flare simulation. The magnetic field, plasma pressure, and density are computed consistently by the model.
The MHS model is capable of reconstructing the main structures of pressure and density in the photosphere and lower choromosphere (below 1 Mm). As we did not use the temperature data to constrain the plasma, the plasma solution above 1 Mm cannot be relied on. The deviation gets larger as the height increases.
Generally, the magnetic field solution is improved compared with the NLFFF not only in terms of magnitude and direction, but magnetic connectivity as well. As to the flux tube, which is of vital importance for eruptive events, the twists of the field lines are constructed well. However, some twist structures derived outside the active region are incorrect. It is worth noting that even right inside the active region, the extrapolated QSLs could be unrealistic. The Q s value is computed by tracing the field lines from which the high nonlinearity is introduced. A tiny deviation shown in Fig. 9 near the bald patch area results in a large difference in connectivity.
Paper I developed the MHS extrapolation and tested the model using a semi-analytic MHS solution. In this work, we present a more challenging test of our model with a radiative MHD simulation of a solar flare. The next activity planned is an application for the IMaX magnetogram (Martínez Pillet et al. 2011) on board the SUNRISE balloon-borne solar observatory (Barthol et al. 2011;Berkefeld et al. 2011) during its second flight (Solanki et al. 2017).
The variables in Eqs. (8) Ω a = λΩ a , (A.8) wheren is the inward unit vector on the surface S.