Preprocessing of vector magnetograms for magnetohydrostatic extrapolations

Context. Understanding the 3D magnetic field as well as the plasma in the chromosphere and transition region is important. One way is to extrapolate the magnetic field and plasma from the routinely measured vector magnetogram on the photosphere based on the assumption of the magnetohydrostatic (MHS) state. However, photospheric data may be inconsistent with the MHS assumption. Therefore, we must study the restriction on the photospheric magnetic field, which is required by the MHS system. Moreover, the data should be transformed accordingly before MHS extrapolations can be applied. Aims. We aim to obtain a set of surface integrals as criteria for the MHS system and use this set of integrals to preprocess a vector magnetogram. Methods. By applying Gauss' theorem and assuming an isolated active region on the Sun, we related the magnetic energy and forces in the volume to the surface integral on the photosphere. The same method was applied to obtain restrictions on the photospheric magnetic field as necessary criteria for a MHS system. We used an optimization method to preprocess the data to minimize the deviation from the criteria as well as the measured value. Results. By applying the virial theorem to the active region, we find the boundary integral that is used to compute the energy of a force-free field usually underestimates the magnetic energy of a large active region. We also find that the MHS assumption only requires the x-, y-component of net Lorentz force and the z-component of net torque to be zero. These zero components are part of Aly's criteria for a force-free field. However, other components of net force and torque can be non-zero values. According to new criteria, we preprocess the magnetogram to make it more consistent with the MHS system and, at the same time close, to the original data.


Introduction
The force-free assumption ∇×B = α(r)B has long been the basis for a magnetic field extrapolation from the solar photosphere to the corona (Wiegelmann & Sakurai 2012). As a result, a number of force-free field extrapolation techniques have been developed over half a century. The simplest case is the potential field or current-free field extrapolation in which α = 0 (Schmidt 1964;Schatten et al. 1969). A next step with a spatially constant α is the linear force-free field extrapolation (Chiu & Hilton 1977;Seehafer 1978;Alissandrakis 1981). At last, a model with non-constant α is called the nonlinear force-free field (NLFFF) extrapolation. The NLFFF extrapolations include the (1) upward integration method (Nakagawa 1974;Wu et al. 1985Wu et al. , 1990Cuperman et al. 1991;Demoulin & Priest 1992;Song et al. 2006), (2) Grad-Rubin method (Grad & Rubbin 1958;Sakurai 1981;Amari et al. 1997Amari et al. , 1999Amari et al. , 2006Wheatland 2004Wheatland , 2006, (3) relaxation method (Mikic et al. 1988;Roumeliotis 1996;Valori et al. 2005;Jiang & Feng 2012;Fan et al. 2011;Guo et al. 2016), (4) optimization method (Wheatland et al. 2000;Wiegelmann 2004;Wiegelmann & Inhester 2010), (5) boundary-element method (Yan & Sakurai 2000;Yan & Li 2006;He & Wang 2006) and (6) forward-fitting method (Malanushenko et al. 2012;Aschwanden 2013Aschwanden , 2016. In a force-free model, some necessary conditions of the magnetic field on the boundary have to be fulfilled. Molodenskii (1969), Molodensky (1974), and Aly (1984Aly ( , 1989 obtained several integral relations on the boundary corresponding to the vanishing net magnetic force and vanishing torque, respectively. If these integral relations are not fulfilled then the boundary is not consistent with a force-free field. Based on this,  proposed a preprocessing algorithm to modify the magnetogram within the error margins of the measurements to minimize the net magnetic force and torque. The resulting vector magnetogram is more suitable for a force-free extrapolation. Further developments include using the method of simulated annealing (Fuhrmann et al. 2007), extending to the spherical geometry (Tadesse et al. 2009), adding a new term concerning chromospheric longitudinal fields (Yamamoto & Kusano 2012), and dealing with the potential and non-potential components separately (Jiang & Feng 2014).
With increasing spatial resolution of the measured vector magnetogram, we can study the magnetic field in the lower solar atmosphere in detail. However, the force-free assumption is not valid anymore in this layer as the plasma β is much larger than it is in the corona (Gary 2001). A straightforward approach is to take into account the plasma in the extrapolation. Zhu et al. (2013), Article number, page 1 of 11 arXiv:2010.06174v1 [astro-ph.SR] 13 Oct 2020 Zhu et al. (2016), andMiyoshi et al. (2020) proposed the magnetohydrodynamic relaxation method to obtain a magnetohydrostatic (MHS) solution with a non-force-free layer near the bottom boundary. Gilchrist & Wheatland (2013) and Gilchrist et al. (2016) extended the Grad & Rubbin (1958) method to compute the MHS equilibria. Wiegelmann & Neukirch (2006) used the optimization method to solve the MHS equations without gravitational force in the Cartesian coordinate system. Wiegelmann et al. (2007) extended the code with gravitational force in spherical coordinate system. Recently, we extended the optimization method by introducing the gravitational force (Zhu & Wiegelmann 2018) in the Cartesian coordinate system, tested the code with an realistic radiative MHD simulation (Zhu & Wiegelmann 2019), and also applied the code to a Sunrise/IMaX dataset (Zhu et al. 2020). It is worth noting that the new algorithm ensures the positive definiteness of gas pressure and mass density.
As a result of the MHS assumption, we can also define several integral relations just as we did for the force-free field. A preprocessing algorithm (extended from the NLFFF case) is proposed to modify the vector magnetogram within the error margins of the measurement. The resulting magnetogram is expected to be more consistent with the assumption of a MHS extrapolation.
The remainder of the paper is organized as follows: In Sect. 2 we define the integral relations for a MHS system and also apply the virial theorem to an active region. In Sect. 3 we describe the optimization algorithm to derive consistent boundary conditions for a MHS extrapolation. In Sect. 4 we apply the algorithm to an example of the observed vector magnetogram. In Sect. 5 we draw conclusions.
As discussed in many papers (e.g., Chandrasekhar 1961;Molodenskii 1969;Molodensky 1974;Aly 1984), the above force balance equation may be written in the form analogous to equations of elasticity as follows: where T i j = 1 2 B 2 δ i j − B i B j is the Maxwellian tensor and f includes forces such as pressure gradient (−∇p) and gravitational force (−ρẑ).

Equations of momenta
Assuming a Cartesian coordinate, integrating Eq. (3) over the volume V with surface S, we get We apply these three relations to an isolated active region with the photosphere as the bottom boundary. As to a closed current system, the magnetic field falls off as 1/r 3 . Therefore as r → ∞, the lateral and top surface integrals with magnetic terms approach zero. We note that in −∇pdv = − S pds the transverse surfaces make zero net contribution since, without magnetic fields, the forces acting on the opposite side boundaries by the external plasma pressure are equal and opposite. Then we have net Lorentz force restrictions written as where S 1 is the bottom boundary. We note that Eqs. (7) and (8) are exactly the same as those for a force-free field. However, that is not the case vertically. Eq. (9) implies that the difference between the pressure at the bottom boundary and the weight of the plasma above is compensated by the Lorentz force.
The net Lorentz force torque restrictions can also be obtained in a similar way with the cross product of Eq.
(3) and r as follows: We see Eq. (12) is exactly the same as that for a force-free field, while that is not the case in two transverse directions. Eqs. (10) and (11) imply that the plasma may induce rotational moments relative to the x or y axes.

Virial theorem
Multiplying Eq. (3) by r = xx + yŷ + zẑ and integrate over volume V with a surface S, we derive the virial theorem which relates the magnetic energy and forces inside the volume to the integral on the surface. If the electric currents are closed, the magnetic field falls off as 1/r 3 , and as r → ∞ the surface integral approaches zero. Therefore the magnetic energy vanishes if f = 0. That means a force-free field does not exist in an isolated current system (Chandrasekhar 1961).
where E t = 1 2 B 2 + p γ−1 − ρz is the total energy including magnetic energy 1 2 B 2 , internal energy p γ−1 and gravitational potential energy −ρz. If γ = 4 3 then the left-hand side (LHS) of Eq. (14) becomes exactly the total energy. Eq. (14) relates the total energy and pressure in the volume to the integral on the surface.

Apply to an active region
Let us focus on an active region. Assume a semi-infinite half-space with the bottom boundary on the photopshere, the magnetic field of an isolated active region decreases with distance as 1/r 3 , thus the surface integral in Eq. (13) on the side and top boundaries will approach zero. Then we get: where S 1 is the bottom boundary.

Effect of translation on the integration
In the force-free case, the integral at the right-hand side (RHS) of Eq. (15) is invariant under translation. However, it seems from the LHS of Eq. (15) that this nice property is not fulfilled any more in a MHS system. Suppose a translation: (x, y, z) −→ (x + ∆x, y + ∆y, z + ∆z). Then B (r ) = B(r), dx = dx, and dy = dy. Thus we have the RHS of Eq. (15) in the new coordinate system as follows: where ∆r = (∆x, ∆y, ∆z). We note that Eqs. (7) and (8) are used to infer formula (18) from (17). Thus the integration on the bottom boundary is invariant under the translation in the same horizontal plane with ∆z = 0.

Estimate the magnetic energy
Suppose r = 0 on the bottom boundary, Eq. (15) becomes Article number, page 3 of 11 As with a force-free field, the integral at the RHS of Eq. (19) is often used to estimate the magnetic energy in the semi-infinite volume. However, the integral is not equal to the magnetic energy in a MHS equilibrium. In the following analysis, we try to judge, in a MHS equilibrium, whether the integral is an underestimate or an overestimate of the magnetic energy. Consider the case of an axisymmetric monopole sunspot (see Fig. 1). The field line approaches a radial direction since the gas pressure becomes small. Near the spot, however, the field line has to bend to compensate the pressure gradient. According to Sec. 2.3.1, the boundary integration of Eq. (19) is invariant under the translation on the bottom plane. For the sake of convenience, we chose r = 0 at the center of the spot. As gas pressure is depleted in the spot, the negative pressure gradient usually points to the interior of the spot. To compensate this force, a reversed Lorentz force is generated transversely. Then a vertically downward component of the Lorentz force is required to make sure that the Lorentz force is perpendicular to the magnetic field. The negative net Lorentz force has been confirmed in most active regions by statistical studies (Moon et al. 2002;Tiwari 2012;Liu et al. 2013;Liu & Hao 2015). According to Fig. 1, far away from the spot, the contribution of f · r vanishes as the field line becomes radial. However, near the spot, the contribution of f · r is always negative. Thus we get, in a monopole spot case, The above inequality also holds in an active region with multiple spots. See Fig. 2 with two spots, the f · r of the left spot is equal to that of the monopole spot case. As to the right spot, with the coordinate transformation r = r 1 + R, we have where R is transverse while V fdv is vertical. Term R· V fdv vanishes is because the horizontal vector R is perpendicular to the vector V fdv, which is vertical due to the axisymmetry assumption. It is worth noting that, in the multiple spots case, the magnetic field line deviates from the radial direction because of the attraction from other spots. This effect, however, is minor near the spot. In regions far away from spots where field lines are not radial any more, the magnitude of f · r decreases with altitude rapidly because f decreases exponentially (depending on temperature) while r can only grow linearly. Thus the analysis expression above still works. Therefore, we conclude that, as in an active region with sunspots, the surface integral of Eq. (19) is usually an underestimate of the magnetic energy of the active region.

Preprocessing method
To see if a vector magnetogram can be served as the boundary condition for a MHS system, three parameters are proposed (similar to those for a force-free field ) as follows: A vector magnetogram is suitable for a MHS extrapolation at least if f lux 1, f orce 1 and torque 1. Since the vector magnetogram sometimes do not fulfill the aforementioned criteria, a preprocessing procedure is required to modify the data within the freedom of the noise. The algorithm is based on the preprocessing method that was developed by . As the restrictions on the boundary values for a MHS equilibrium are weaker compared with those for a force-free field, we need to change the preprocessing procedure accordingly. Non-zero values of z-component of net Lorentz force and x-, y-component of net torque are allowed to exist in a MHS equilibrium. Therefore these three numbers should be retained during the preprocessing process.
To do so, we define the functional where are the z-component of net Lorentz force and the x-, y-component of net Lorentz torque, respectively. The summation is over all p grid nodes on the photosphere. The weighting factors µ n are as yet undetermined. The terms L 1 and L 2 correspond to the net force and net torque constraints. The term L 3 measures the deviation between the original data and the preprocessed data. The term L 4 controls the smoothing. For computational reasons, sufficiently smooth data are necessary for an optimization to obtain a good solution. The new algorithm differs from that for the force-free field ) by introducing three quantities: a 0 , a 1 , and a 2 , which ensure that the preprocessing does not change the corresponding integration values. The strategy of preprocessing is to use the gradient descent method to minimize L, and meanwhile make all L n small as well. The magnetic field is optimized as follows: with a gradient descent method The three functional derivatives at node (q) are defined as Smoothing was performed for all three components. Effects from terms that have mixed products of vertical and transverse magnetic field components in functional Eqs. (30-33) are not included when evaluating B z . This is designed for the fact that B z is measured with much higher accuracy than B x and B y (Martínez Pillet et al. 2011).

Application to Sunrise/IMaX data
For the test we used a combined vector magnetogram in which the Sunrise/IMaX data (Martínez Pillet et al. 2011;Solanki et al. 2017) is embedded in the HMI data (Scherrer et al. 2012). We have used this dataset to extrapolate the magnetic field as well as the plasma using two approaches (Wiegelmann et al. 2017;Zhu et al. 2020). The top panels of Fig. 3 show the original vector magnetogram within the IMaX field of view (FOV). The combined data have a flux imbalance of f lux = 0.013 (almost balanced). The MHS criteria are not fulfilled to some extent with f orce = 0.091 and torque = 0.066. For comparison, Aly's criteria for force-free field are largely violated with f orce = 0.29 and torque = 0.32. Before the magnetogram is preprocessed, we first need to choose the appropriate µ n . There are four parameters of µ n . Only three of them are independent since only the ratio of the parameters really counts. We further assume µ 1 = µ 2 D 2 ≡ µ 12 , where D is the average edge length of the magnetogram. Thus we could define L 12 = L 1 + D 2 L 2 . This assumption gives the same weight to the momentum and torque constraints. As only the ratio of the parameters counts, we specify µ 12 = 1/B ave , where B ave is the average magnetic field strength in the magnetogram. Then only two independent parameters remain (µ 3 and µ 4 ). A survey of the two parameters for the combined magnetogram shows that, as was also found in , log(L 3 ) and log(L 4 ) are almost determined by the ratio of µ 3 to µ 4 , while log(L 12 ) depends on the magnitude of µ 3 and µ 4 .
With a deviation of the magnetic field value (i.e., a finite L 3 ), we obtain a smoothed solution that satisfies the criteria of a MHS system. The deviation is tolerable as long as L 3 does not exceed noise level of the magnetogram, that is, L 3 = 2.2 × 10 −9 in this case. The noise is retrieved from the HMI dataset (Hoeksema et al. 2014). Fig. 4 shows optimal µ 3 and µ 4 combinations at which log(L 3 ) equals to noise level of the magnetogram. As µ 3 and µ 4 decrease, L 12 converges. Tab. 1 shows L-values as well as three summations of the initial data and the preprocessed data. The three summations are not changed during the preprocessing owing to the special design of the algorithm for fixing them. An improvement on L 12 is obvious in that L 12 is reduced by five orders of magnitude, which enforce a good compliance with the MHS criteria. The quantity L 3 has to be finite because we allow the field values to deviate from the observed values. Even though the smoothness is hardly visible in Fig. 3, L 4 after preprocessing is 1/20 smaller than it was before preprocessing. Now we get f orce = 7.0 × 10 −5 and torque = 8.7 × 10 −5 with an optimal values of µ 3 = 1.0 × 10 −3 and µ 4 = 8.9 × 10 −4 , which are more suitable for a MHS extrapolation. Fig. 5 shows field lines of different models. The original line-of-sight magnetogram was used for the potential field modeling, the force-free preprocessing was applied for the NLFFF modeling, and the MHS preprocessing was applied for the MHS modeling. Fig. 1. Cartoon shows, in a MHS state, a typical interaction of the magnetic force and plasma forces in an active region with single compact polarity. Table 1. L-values and other indexes before and after preprocessing.

Sunspot
Original 8.5 × 10 −3 0.0 1.1 × 10 −7 0.20 0.11 0.14 Preprocessed 2.4 × 10 −8 2.2 × 10 −9 6.4 × 10 −9 0.20 0.11 0.14 All extrapolations were done in a 2336 × 1824 × 128 box. A central box with 936 × 936 × 128 dimensions above the IMaX FOV was cut to display the result. It is clear from fig. 5 that field lines that share the same footpoints have different structures. A relatively large difference (e.g., the connectivity) can be found between panel (a) and (b) due to the force-free currents in the NLFFF. The difference between panel (b) and (c) as a consequence of the perpendicular component of the current are much smaller. Field lines in both panel (b) and (c) have quite a similar pattern, but still have a different connectivity as well as length of field lines. In Fig. 6 we compare the simultaneous chromospheric observation by Sunrise/SuFI 3968 Å (Gandorfer et al. 2011;Solanki et al. 2017) with field lines within subvolumes spanning the 600-1400 km height range. We note that the extrapolation data were cut according to the SuFI FOV. The observed slender fibrils are generally believed to outline the magnetic fields in this layer. We find that most field lines trace the fibrils nicely. However, deviations can also be observed in some regions.

Conclusions
We have investigated the virial theorem for a MHS system. By applying it to a large active region, we find that the surface integral that was often used to compute the energy of the force-free magnetic field is a lower bound of the magnetic energy in the semiinfinite volume. We also find a set of surface integrals as necessary criteria for a MHS system. These integrals are equivalent to Aly's criteria for a force-free field. According to the new set of criteria, we proposed an optimization algorithm to preprocess the vector magnetogram with the aim to use the result as a suitable input for a MHS extrapolation. The optimization strategy is similar to that designed in  for the force-free modeling, which is to force the data to fulfill the criteria for the MHS system and be sufficiently smooth within the freedom of the measurement noise. We also show an application of the preprocessing method to a combined vector magnetogram in which the IMaX data are embedded in the HMI data.