Magnetohydrostatic modeling of AR11768 based on a SUNRISE/IMaX vector magnetogram

Context. High resolution magnetic field measurements are routinely done only in the solar photosphere. Higher layers like the chromosphere and corona can be modeled by extrapolating the photospheric magnetic field upward. In the solar corona, plasma forces can be neglected and the Lorentz force vanishes. This is not the case in the upper photosphere and chromosphere where magnetic and non-magnetic forces are equally important. One way to deal with this problem is to compute the plasma and magnetic field self-consistently with a magnetohydrostatic (MHS) model. Aims. We aim to derive the magnetic field, plasma pressure and density of AR11768 by applying the newly developed extrapolation technique to the SUNRISE/IMaX data. Methods. An optimization method is used for the MHS modeling. The initial conditions consist of a nonlinear force-free field (NLFFF) and a gravity-stratified atmosphere. Results. In the non-force-free layer, which is spatially resolved by the new code, Lorentz forces are effectively balanced by the gas pressure gradient force and the gravity force. The pressure and density are depleted in strong field regions, which is consistent with observations. Denser plasma, however, is also observed at some parts of the active region edges. In the chromosphere, the fibril-like plasma structures trace the magnetic field nicely. Bright points in SUNRISE/SuFI 3000 {$\AA$} images are often accompanied by the plasma pressure and electric current concentrations. In addition, the average of angle between MHS field lines and the selected chromospheric fibrils is $11.8^\circ$, which is smaller than those computed from the NLFFF model ($15.7^\circ$) and linear MHS model ($20.9^\circ$). This indicates that the MHS solution provides a better representation of the magnetic field in the chromosphere.


Introduction
Due to the high plasma β in the photosphere and chromosphere, non-magnetic forces should be taken into account when we study these layers. If we neglect dynamics and plasma flow, then the resulting static state of the system can be described by the so-called magnetohydrostatic (MHS) assumption, which is determined by the balance of Lorentz force, pressure gradient and gravity force, together with the solenoidal condition of the magnetic field.
A more challenging problem is to solve the MHS equations in the generic case, which can be done only numerically. A few methods have been developed. The magnetohydrodynamic (MHD) relaxation method was introduced to derive the MHS solution by using an "evolution technique" (McClymont & Mikic 1994;Jiao et al. 1997;Zhu et al. 2013Zhu et al. , 2016Miyoshi et al. 2019). The Grad-Rubin iteration, which is well known in the calculation of the nonlinear force-free field (NLFFF), was extended by Gilchrist & Wheatland (2013) and Gilchrist et al. (2016) to solve the MHS equations. Wiegelmann & Inhester (2003) and Wiegelmann & Neukirch (2006); Wiegelmann et al. (2007) develop the optimization method to treat the MHS equations. We recently extended this method by introducing the gravity force (Zhu & Wiegelmann 2018). In addition, our new algorithm ensures that the resulting plasma pressure and density are positive definite. More recently, we further test our code with the radiative MHD simulation of a solar flare (Zhu & Wiegelmann 2019). This challenging test (solar flares are intrinsically non-static and even non-stationary) provides a solid foundation for the application of the code to real observations. It is worth noting that, besides the traditional NLFFF models, a new type of NLFFF model (Malanushenko et al. 2012;Aschwanden 2013Aschwanden , 2016 has been introduced recently. This alternative approach uses the line-of-sight magnetogram to constrain the potential field. With the help of a forward-fitting method the angle between the magnetic field and coronal or chromospheric loops is minimzed. In this model, the assumption of a force-free photosphere is not used either. This works well in the corona, but may not be the ideal way to describe the non-force-free chromospheric field. In this paper, we apply our code, for the first time, to a vector magnetogram obtained by the IMaX instrument on the Sunrise balloon-borne solar observatory. Since the field-of-view (FOV) of IMaX is limited to part of the active region, an SDO/HMI vector magnetogram is used to cover the unobserved parts. The organization of the paper is as follows. In Section 2, we describe the dataset used in this paper. In Section 3, we analyze the results and compare them with other models. Conclusions and perspectives are presented in Section 4.

Magnetohydrostatic Extrapolation and its application to AR11768
The MHS extrapolation computes the magnetic field, plasma pressure and density consistently with the help of an optimization principle. The algorithm was described and tested in detail in Zhu & Wiegelmann (2018. The primary dataset used in this work was recorded by the vector magnetograph IMaX (Martínez Pillet et al. 2011) onboard the Sunrise balloon-borne solar observatory Berkefeld et al. 2011;Gandorfer et al. 2011) during its second flight, refer to as Sunrise II . The IMaX data have a pixel size of 40 km a FOV that contains 936 × 936 pixels (50 × 50 ). This FOV is limited to part of AR11768. The data were inverted by Kahil et al. (2017) using the SPINOR inversion code (Frutiger et al. 2000) that builds on STOPRO routines (Solanki 1987). A one-component atmospheric model with three optical depth nodes (at logτ=-2.5, -0.9, 0) for the temperature and height-independent magnetic field is applied. The affect of the inversion (not the same inversion used here) on the extrapolation result was studied by Wiegelmann et al. (2010) and found to be minor. The 180 • ambiguity is removed with an acute angle method  which minimizes the angle with, in this case, the corresponding HMI vector magnetogram. We note that there is a square pattern in the transverse field of IMaX ( Fig.1). It originates from disambiguating IMaX with HMI data in which a square patten also exists (seen in the transverse field). Essentially, the square pattern is the IMaX noise that is modulated to match the HMI spatial resolution by disambiguating.
Although this noise appears in big mosaics (big compared to the IMaX resolution), the size of every single one is still very small compared to the active region. Therefore the effect they have on the extrapolation is expected to be similar to the effect of the normal noise which has been studied by Zhu & Wiegelmann (2018). In that paper we found, based on the quantitative assessment for the extrapolation, the influence of the random noise (20% level) on the magnetic field to be less than 4%. As pointed out by De Rosa et al. (2009), it is necessary to have flux-balance within the FOV and to catch the magnetic connectivity in order to find unique solutions. So we have embedded the IMaX data into the SDO/HMI vector magnetogram (Pesnell et al. 2012;Scherrer et al. 2012) taken closest in time to the analyzed IMaX magnetogram. Fig. 1 shows both IMaX and combined vector magnetogram. Fig. 2 depicts the flow chart of our code applied to the combined magnetogram. The following steps are plotted. First, compute a NLFFF (Wiegelmann 2004) with a preprocessed magnetogram , the net Lorentz force and torque are removed within the error margin of the measurement by using a minimization principle to make the data compatible with the force-free assumption) and a gravity-stratified atmosphere along field lines with a 1D temperature profile (Fig. 3). The pressure on the bottom boundary is computed using p = p quiet − 1 2 B z 2 , where p quiet is the pressure in the quiet region. The bottom density is determined by assuming a uniform temperature of 6000 Kelvin on the photosphere. Second, carry out a further optimization to achieve an MHS equilibrium with the original magnetogram. Note that although a temperature profile is given at the initial state to relate the gas pressure and density, the initial temperature profile is no longer a restriction on the pressure and density optimization. There are two options to input the bottom magnetic boundary. One is to replace the preprocessed magnetogram with the original magnetogram directly, which is used in this study. The other is to change the magnetic field gradually from the preprocessed value to the original value. The side and top boundaries are fixed to their initial values which are potential fields. The computation is performed in a 2336 × 1824 × 128 box with a 40 km grid spacing both in horizontal and vertical directions. All of the following analyses are restricted to a 936 × 936 × 128 box above the IMaX FOV (unless otherwise stated).

Solution consistency
Both the magnetogram and the disequilibrium of the initial atmosphere (due to the nonuniform gas pressure and density in the photosphere) drive the evolution of the system to an MHS state. Fig. 4 illustrates the compensation of forces in the initial state and in the final equilibrium. The horizontal component of the forces are shown in panels (a) and (c), while the vertical components are illustrated in panels (b) and (d). In the initial state, we find the residual force is greater than the Lorentz force (see panels (a) and (b)), so that this state is obviously far from an MHS equilibrium. As mentioned before, this promotes further optimization. We see, in the final solution, that the Lorentz force is balanced effectively by the pressure gradient and vertically also by gravity (see panels (c) and (d)). On average, the residual force is 43% of the Lorentz force in the transverse direction while the percentage is only 24% in the vertical direction below 2 Mm. These residual forces are not close to 0. The main contributions come from the photosphere and regions with very high or low β. The inadequate boundary condition of plasma as well as the noise in the measured magnetic field prevent the MHS extrapolation from balancing forces well on the photosphere. In a very high plasma region (β > 10) Lorentz force are too weak to act against plasma forces, which could result in a ratio of residual force over Lorentz force that is much larger than 1. In a very low plasma region (β < 0.1) plasma forces are too weak to against Lorentz force, which could result in a ratio that is close to 1. We recompute the ratios excluding the bottom boundary and only in regions where 0.1 ≥ β ≤ 10. The numbers are 23% and 7.6% for the transverse and vertical direction, respectively, which means that the major part of the Lorentz force is balanced. Although those two numbers are still not very small, they are acceptable considering that we have embedded two data sets in each other (recorded by different instruments and obtained by   Similar results were also found both in the LMHS modelings (Aulanier et al. 1998(Aulanier et al. , 1999Wiegelmann et al. 2015Wiegelmann et al. , 2017 and the previous tests of our MHS code (Zhu & Wiegelmann 2018. However, we also find that, at some parts of the active region edges, plasma pressure and density are enhanced. This has never been reported in previous extrapolations. Fig. 6 depicts selected Lorentz force vectors in a FOV that is outlined by the red square in Fig. 5 (c). We see that regions with an enhanced gas pressure which are encircled by blue ellipses are surrounded with a net inward flux of the Lorentz force. As a result, the plasma is squeezed together. It is also worth noting that the fibril-like plasma pattern traces the magnetic field in Fig. 5 (c) and Fig. 7 (a). Such localized concentrations reflect the nonlinear nature of the MHS system, which cannot be observed in the LMHS model (see Fig. 7). According to the model of plasma β over an active region developed by Gary (2001), the magnetic field dominates plasma above a height of 1 ∼ 2 Mm (see Fig. 3 in Gary (2001)). Fig. 8 (a) shows that the β (z) of the MHS equilibrium is right inside the β range illustrated in Gary (2001). A high plasma β is a necessary, but not a sufficient condition for the magnetic field to be nonforce-free. To see whether the magnetic field is non-force-free in the high β region, the current-weighted sine of the angle between the magnetic field and the electrical current density (Schrijver et al. 2006

Plasma
As shown in Fig. 8 (b), σ decreases fast from 0.7 at the photosphere to less than 0.1 above 2.0 Mm. Since the effective vertical resolution of the solution is roughly the same as the horizontal resolution of the magnetogram, the high resolution IMaX data allow us to study this non-force-free layer in detail. Coarser data (e.g. HMI) with a few grid points to resolve this layer meet with difficulties when focussing on the lower atmopshere. Fig. 8 (c) shows Lorentz force distributions at a height of 0.4 Mm. We see strong Lorentz forces are mainly located at edges of strong- field features, where plasma β drops precipitous. This is a natural result caused by the strong plasma forces at edges. The great plasma differences at edges of magnetic features in the lower atmosphere are routinely observed.
3.3. Relation between plasma pressure, current density and photospheric brightness Fig. 9 shows the mapping of plasma and current density onto an image acquired by Sunrise/SuFI. Note that the extrapolation data are cut according to the SuFI FOV of 15 × 38 . Bright points are clearly seen in the inter-granular lanes in panel (c). They are typically regarded as nearly vertical slender flux tubes with kG magnetic fields (Solanki 1993;Nagata et al. 2008;Requerey et al. 2014). The lateral radiative inflow makes the tube hot and bright (Spruit 1976), making them visible as bright points at wavelengths sensitive to temperature (Schüssler et al. 2003;Riethmüller et al. 2010). Fig. 9 (a) (b) show that regions of high plasma pressure and strong electric current density coincide. Most of them are located near the edges of magnetic flux tubes. These flux tubes which appear as photospheric bright points at SuFI 3000 Å (see Fig. 9 (c)) are often accompanied by high plasma pressure (below 400 km) and electric current around them. At the edges of flux tubes, the plasma and magnetic field interact, which leads to the high current and co-spatial high plasma pressure. Such a high current density around magnetic bright points is also reminiscent of the electric current sheets expected to bound flux tubes. It is worth noting that many bright points do not have corresponding enhancements in the electric current. This may be related to the local dynamics. It must also be kept in mind that the MHS model does not take into account the radiation.

Comparisons of magnetic fields and chromospheric fibrils
The SuFI instrument provides diffraction-limited images at 3968 Å with contributions from both the photosphere and low-to-midchromosphere (Jafarzadeh et al. 2017). The observed slender fibrils at this wavelength are the dominant structures (see Fig. 10 (a)) in the SuFI FOV. It is generally believed that long fibrils in the chromosphere outline magnetic fields in this layer (de la Cruz Rodríguez & Socas-Navarro 2011; Jing et al. 2011;Schad et al. 2013;Leenaarts et al. 2015;Zhu et al. 2016). We plot field lines within the sub-volume spanning the 600-1400 km height range, which implies that low field lines (≤ 600 km) are ignored. Seed points of the field lines are uniformly selected in the photosphere. Fig. 10 shows that most fibrils trace field lines nicely. The similar field line patterns of different extrapolation models imply that plasma forces have limited impact on the large scale magnetic field of this potential-like active region, at least in the height range from 600-1400 km. In order to quantitatively show the degree of agreement between fibrils and field lines we compute the angle θ between fibrils and the plane-of-the-sky component of the magnetic field. The computation is not carried out on every pixel on the image. Instead, to show more clearly the discrepancy among different models, we focus on the regions where large differences between magnetic vectors of the MHS model and the NLFFF model appear (see contours in Fig. 10). We have tried to identify all fibrils in the regions of interest. For each fibril, one point is picked for the statistics. Then the method introduced by Inhester et al. (2008) is used to estimate the orientation of the fibril. The key point of the method is to use the gradient and Hessian matrix of the image intensity to determine the orientation. Then θ is computed vertically within the 600-1400 km height range. The smallest one is specified as the discrepancy between the fibril and the magnetic vector. For totally 26 points picked (hence 26 fibrils identified), the average of θ for the MHS model, NLFFF model and LMHS model are 11.8 • , 15.7 • and 20.9 • , respectively. There are 18 points (nearly 70%) at which the MHS model's θ is less than those of the other two models.

Conclusion and perspectives
In this work we apply, for the first time, a nonlinear MHS code to model the solar lower atmosphere starting from a real magnetogram. A combined vector magnetogram to cover the whole active region is used as the boundary input. MHS equilibrium is constructed in which Lorentz forces are effectively balanced by plasma forces. The pressure and density depletion take place in the strong field regions together with enhancements in the active region edges. In the low β layers, the fibril-like plasma structures clearly outline the magnetic field lines. A thin non-force-free layer is resolved within about 50 grid points in the z direction. The high plasma pressure and co-spatial high electric current appear around the bright points prominent in SuFI 3000 Å images. Although similar general patterns of