Automatic Computation of Linear Magneto-Hydro-Static Equilibria

Deriving the physical parameters of observed phenomena in the solar atmosphere has fundamental importance, as these parameters are then employed to constrain and validate models. Here we report the development of a new computational algorithm based on a magneto-hydro-static model that computes the magnetic field in the solar atmosphere and automatically matches individual magnetic-field lines with observed structures that appear with enhanced emission in extreme-ultraviolet (EUV) images. Presently, for the quiet-Sun regions, we can only measure the vertical photospheric magnetic field Bz\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$B_{z}$\end{document}, as accurate horizontal magnetic field measurements are not available. Thus, vertical photospheric magnetic-field measurements are extrapolated into the upper atmosphere, from the photosphere to the corona, with a magneto-hydro-static model. Free model parameters are then optimized with a downhill-simplex method by comparing quantitatively magnetic-field lines with the enhanced emission of loop structures composing the so-called coronal bright points recorded in EUV images taken with the Atmospheric Imaging Assembly on-board the Solar Dynamics Observatory. The algorithm will be applicable to any solar image data where individual structures with an enhanced emission can be resolved. Most importantly, the algorithm can be employed to obtain the magnetic properties of these structures above the photosphere.

field are mainly obtained in the photosphere. However, to derive the magnetic field above the photosphere, the observed photospheric magnetic field is extrapolated into the higher solar atmosphere under some modeling assumptions. In the low plasma-β corona, the field is modelled under the force-free assumption (see Wiegelmann and Sakurai, 2021, for a recent review on solar force-free fields). For lower heights (upper photosphere and chromosphere) finite Lorentz force has to be considered in the lowest order by magneto-hydro-static (MHS) models (see Zhu, Neukirch, and Wiegelmann, 2022, for a recent review on MHS modelling). In cases of large-scale extrapolations, including the solar wind region, MHD models with plasma flows are required. In the generic case, the related equations are nonlinear and photospheric vector magnetograms alone are not sufficient to constrain an MHD model and one needs either additional observations or make some model assumptions. In principle, coronagraph images can be used together with a tomographic inversion (see Aschwanden, 2011, for a review article on tomography) to derive the 3D plasma density and pressure, and measurements of the solar wind speed can be incorporated into these models. A simpler possibility is to approximate an initial state for the plasma density, pressure and solar wind flow velocity profile from model assumptions. An easy to implement approach is to apply as the initial state for the hydrodynamic MHD-variable spherical symmetric hydrodynamic solutions like the Parker model of the solar wind (see Parker, 1958, for details) and an isothermal coronal plasma. The stationary MHD equations are then solved under the assumption of a field-aligned plasma flow with photospheric vector magnetograms as a boundary condition.
The accuracy of photospheric magnetic-field measurements differs for the different regions of the Sun. While the full magnetic-field vector is measured with high accuracy in active regions, in the quiet Sun, only the vertical magnetic-field component B z is available, because the signal-to-noise ratio of the horizontal magnetic-field component is too poor. The simplest approach to model the magnetic field in the higher solar atmosphere is the current-free potential field, which is uniquely determined from vertical (in spherical geometry radial) magnetic-field measurements. Far away from disk center the line-of-sight and vertical magnetic field differ. Fully neglecting the electric currents, however, is not justified for most of the observed phenomena, and the derived magnetic-field lines will not match the observed coronal, transition-region or chromospheric observations. A possible solution is to apply models which require only vertical measurements as boundary conditions, but have one or several free parameters to model electric currents. A well-known example is a linear force-free (LFF) field model, where the currents are parallel or anti-parallel to the magnetic field and are controlled by the force-free parameter α.
In the time prior to the launch of the Solar Dynamics Observatory (SDO), vector magnetic-field data were rarely available, even in active regions. For those observations, methods were developed to derive the optimum α by comparing the modelled field lines with structures in coronal images to find the optimal value of α. Carcedo et al. (2003) developed a method to compare coronal X-ray images taken with the Solar X-ray Telescope (SXT) on board Yohkoh with extrapolated field lines. Feng et al. (2007b), Feng et al. (2007a), and Conlon and Gallagher (2010) compared LFF model field lines with images from two or three vantage points combining observations from the Extreme ultraviolet Imaging Telescope (EIT) on board the Solar Heliospheric Observatory (SoHO) and Extreme Ultraviolet Imagers (EUVIs) on board the two STEREO spacecraft. Models with the optimum values of α were used to constrain the stereoscopic reconstruction from these data. The main challenge of fitting structures like, for instance, active-region loops is, that the optimal value of the free parameter α might vary for loops in different parts of the active region. Nevertheless, Marsch, Wiegelmann, and Xia (2004) found an example where one value of α was sufficient to match several coronal loops in one active region. Such models have then subsequently been used to guide further investigations. As an example, Xie et al. (2017) used a LFF field model with optimized α to compute the plasma parameters of active-region magnetic loops.
Force-free models are, however, only a good approximation for the low-β corona, but not in the photosphere and chromosphere. In these regions, forces and thereby currents perpendicular to the magnetic field have to be considered. Assumptions regarding the properties of these currents require that the corresponding MHS models have several free parameters, e.g. a parameter α which controls field-aligned currents, a parameter a which controls strictly horizontal currents, and a parameter κ which controls how fast the horizontal currents decrease with height. To find the optimum free parameters, our approach is to generalize the well-known methods for finding the best α for linear force-free fields. In past studies, with only one free parameter, e.g., α, it was possible to scan the full parameter space. Here, we apply a downhill-simplex method in order to minimize a functional containing the three parameters.
The automatic linear magneto-hydro-static (LMHS) computation reported here is tested on observations of small-scale coronal loops in the quiet Sun of the solar atmosphere, known as coronal bright points (CBPs). Loops are the main building blocks of the solar atmosphere (Reale, 2014), with CBPs covering the lower part of the loop-size range. The term CBP describes small-scale multi-loop systems in the low solar corona that appear with enhanced emission in X-ray and EUV emission while rooted in magnetic flux concentrations of opposite polarity (Vaiana et al., 1973). Actually, CBPs are composed of loops at different temperatures, including coronal, transition region (Habbal, Withbroe, and Dowdy, 1990;Kwon, Chae, and Zhang, 2010;Doschek et al., 2010;Mou et al., 2018), as well as chromospheric (e.g. Madjarska et al., 2021). Details on the physical and magnetic properties of CBPs can be found in Madjarska (2019). Figure 1 shows observations of three CBPs (referred to as CBP-1, CBP-2, and CBP-3) observed in emission from plasma at coronal temperatures, recorded in images taken with the 193 Å filter of the Atmospheric Imaging Assembly (AIA) on board the Solar Dynamics Observatory (SDO) (for details see Section 4).
We outline the paper as follows. In Section 2 we introduce a special class of LMHS equations. In Section 3, we revisit a subclass of them, called linear-force-free field models, and discuss which steps are required to upgrade to the more powerful LMHS models. In Section 4, we introduce the observational data used in this study. Section 5 outlines how the free parameters of LMHS can be computed from identified coronal loop structures using a semi-automatic approach. This section generalizes the corresponding methods for linear force-free (LFF) fields. In Section 6, we present a fully automized approach and investigate to which extent the optimum LMHS equilibria can be automatically computed from magnetograms and coronal images without the need to pre-identify individual structures. Finally, we summarize our results, draw conclusions and give an outlook on our further work in Section 7.

Linear Magneto-Hydro-Static Model
The magneto-hydro-static (MHS) equations are given by Here, B is the magnetic field, j is the current density, μ 0 is the permeability of free space, p is the plasma pressure, ρ is the density, and is the gravitational potential. There exist several methods to solve these equations in their generic nonlinear application, e.g., the Grad-Rubin method by Gilchrist and Wheatland (2013) and Gilchrist, Braun, and Barnes (2016), the MHD-relaxation method by Zhu et al. (2013) and Miyoshi, Kusano, and Inoue (2020), and the optimization method by Wiegelmann and Neukirch (2006), Zhu and Wiegelmann (2018), and Zhu and Wiegelmann (2019). These methods are all numerically expensive and require photospheric vector magnetograms as boundary conditions. The heritage of these methods lies in the method-solving nonlinear force-free (NLFF) equations, which is a subset of MHS equations, where the right-hand-side of Equation 1 is equal to zero. The computing time for MHS equilibria is, however, significantly longer (see, e.g., Zhu, Neukirch, and Wiegelmann, 2022) in comparison to an NLFF one. If reliable vector magnetograms are not available, we can still solve the MHS Equations 1-3, by making assumptions considering the electric currents in the form This special class of solutions was found by Low (1991). The first part of the right-hand side αB controls strictly magnetic field-aligned electric currents by the parameter α. If only this part exists, we have a linear force-free model. The second term a exp(−κz)∇B z × e z controls electric currents which are strictly perpendicular to gravity, which means that these currents are horizontal in the x-y plane and do not have a vertical component. The general theory of this approach is described in the original paper (Low, 1991), which is based on a representation of the magnetic field and electric current density by Euler potentials. The special case considered here (Section 2.3 in the Low paper) occurs when one of the Euler potentials coincides with the gravitational potential. The general theory is valid also for more general potentials, e.g., including the centrifugal force to model rotating magnetospheres (see Section 2.6 in the Low paper). For a Cartesian geometry and a constant z-directed gravity force, this leads to a term f (z)∇B z × e z and in the special form of Equation 4 it was chosen f (z) = a exp(−κz), where the parameter a controls the strength of these currents and the parameter κ how fast these currents decrease with height. With increasing height, z the term exp(−κz) decreases exponentially and these currents do almost vanish. In the corona we get exp(−κz) 1 and consequently, the linear magneto-static field reduces to a linear force-free field. Typically, the field-aligned currents do also have a horizontal component, which is by construction parallel to the magnetic field and does not generate a Lorentz force. The currents defined by the last term in Equation 4 are by construction strictly horizontal and contain a field-aligned and a field-perpendicular parts. The latter generates the finite Lorentz force. For a = 0 Equation 4 reduces to the linear force-free model (see, e.g., Seehafer, 1978;Alissandrakis, 1981, for linear force-free methods). An interesting question, as discussed in Aulanier et al. (1999), is whether the MHS equilibria in the form of Equation 4 should be called linear MHS. It was pointed out that the second term in Equation 4 contains the horizontal currents and part of these currents is parallel to the magnetic field while some are perpendicular to it. The perpendicular part of the currents creates the finite Lorentz force, whereas the field-aligned part causes a nonlinearity parallel to the magnetic field. Nevertheless, the resulting equations can be solved with Fast Fourier Transform (FFT) methods and reduced to the linear force-free model for α = 0. Consequently, they are called LMHS in Aulanier et al. (1999). We use this term here as well. It should be noted that a new class of LMHS solutions with a sharper gradient between the forced lower atmosphere (photosphere and chromosphere) and the force-free corona was introduced by Neukirch and . This linear magneto-hydrostatic model is a special case of the generic case of nonlinear magneto-hydrostatics. This is similar to the relation of linear force-free fields and the generic nonlinear force-free fields. If the available observational data and computer resources allow, in both cases nonlinear models are preferable, as they are more general and more realistic, because in the generic cases the solar magnetic field and plasma environment behaves nonlinear. For quiet Sun regions, nonlinear models cannot be applied, because necessary measurements (like the horizontal photospheric magnetic field components) are not available with sufficient accuracy (the noise is higher than the signal). Consequently to model quiet Sun regions, linear models are the best approach we can currently do.
Fitting structures in different locations with a different set of free parameters results in a non-self-consistent model. Nevertheless (when nonlinear models are not available) this approach is helpful because the best-fitting magnetic field lines found by this approach can be used to guide other chromospheric or coronal observations. What these linear models cannot provide is both self-consistency and best-fit parameters in different parts of the field of view (FOV).
An advantage of the approach described with Equation 4 is that only the vertical photospheric magnetic field B z (x, y, z = 0) is used as a boundary condition and that the equations can be solved effectively with the help of FFT, when the parameters α, a, and κ are given. Alternatively to FFT, a Greens' function method was developed by Petrie and Neukirch (2000). The method was applied to data from the Michelson Doppler Imager (MDI) on board SoHO by Aulanier et al. (1999) and to SUNRISE Imaging Magnetograph Experiment (SUNRISE/IMaX) magnetograms by Wiegelmann et al. (2015Wiegelmann et al. ( , 2017. One potential shortcoming of solving equations with FFT is that periodic boundary conditions are required and that the magnetograms used as boundary data have to be flux balanced. This condition is not necessarily fulfilled for arbitrary quiet-Sun areas. We therefore use mirror magnetograms as introduced by Seehafer (1978) for LFF models. These magnetograms mirror the original magnetogram in the range x = 0 . . . L x , y = 0 . . . L y to x = −Lx . . . L x , y = −L y . . . L y by using B z (−x, y, 0) = −B z (x, y, 0) and B z (x, −y, 0) = −B z (x, y, 0). The mirrored magnetograms are flux balanced by construction and are four times larger than the original. This, however, is not a shortcoming for solving the MHS equations with FFT, even on desktop or laptop computers. Solving this set of MHS equations is numerically still somewhat more expensive than solving the LFF equations, even if FFT is used because the high-order Bessel functions have to be evaluated (Low, 1991). This is only necessary for regions where a significant Lorentz force and corresponding electric currents exist. Higher up in the solar atmosphere, the term a exp(−κz) in Equation 4 almost vanishes and the magnetic field becomes linear force-free. A way to monitor the finite Lorentz forces in a particular layer z is the virial theorem (Aly, 1989) The components of the surface integral can be combined to a dimensionless quantity, which tells how force-free a vector magnetogram or any horizontal layer in a computational volume is: This quantity has to be applied to a flux-balanced magnetogram and is frequently used to check whether the vector magnetograms are consistent with the assumption of a force-free model. Metcalf et al. (1995) applied, for instance, the virial theorem to magnetic-field measurements at different heights to monitor the Lorentz force in the photosphere and lower chromosphere. One can naturally apply it also to MHS equilibria. Layers at different heights in the solar atmosphere with force 1 are basically force-free, and the contribution of the second term in Equation 4 is negligible. LMHS models do actually become LFF models above a height z where the term a exp(−κz) in Equation 4 vanishes, which is usually (almost) the case above z = 2 Mm. The effect that the magnetic field energy becomes infinite in the upper half-space for LFF models (see, e.g., Alissandrakis, 1981;Gary, 1989) consequently true also for LMHS. In any finite computational box, the magnetic energy remains final for all models, which is also true for the periodic boundary conditions used in FFT methods and a finite box size in the z-direction. We extrapolated up to a height of z = 13 Mm, which is sufficient to investigate the small-scale loops visible in the EUV images.

Revisiting Linear Force-Free Models
As mentioned above, LMHS equilibria contain three free parameters α, a, and κ. Wiegelmann et al. (2017) already demonstrated how α and a can be constrained from photospheric vector magnetograms in an active region using data taken with SUNRISE/IMaX. As vector magnetograms are not available in the quiet-Sun regions, other methods are then required to optimize these free parameters. Photospheric LOS magnetograms have been routinely available during the SoHO/MDI era and LFF field models were typically employed for active regions data, because they require only B z as boundary conditions. As mentioned above, LFF models also contain a free parameter α which is not constrained by the magnetograms and can be used to match coronal observations. To do so, Carcedo et al. (2003, hereafter C03) developed a quantitative method to compare LFF magnetic-field lines with coronal images from Yohkoh/SXT. The basic steps of the C03 method (see Section 2 in Carcedo et al., 2003, for details) are: i) Align the magnetogram and the coronal image; ii) Identify a coronal loop and its footpoint areas; iii) Compute an LFF model with an arbitrary value of α; iv) Compute a number of field lines starting from locations in the two footpoint areas. Only magnetic-field lines which connect both footpoint areas are further considered. v) Each of these selected field lines is then quantitatively compared with the coronal image.
A Gaussfit is applied to determine how far the centre of the loop (highest brightness in the intensity profile of the coronal loop) and the magnetic-field line are apart. This comparison is done at M positions along the field line. vi) Compute the standard deviation C i (α), which tells how well the LFF model with a certain parameter α agrees with a loop i seen in the image. This can be done by averaging over different magnetic-field lines or by taking the best fitting field line with minimum C i (α). vii) The procedure is repeated for many values of α, thereby scanning the whole parameter space of the LFF model. The value of α with the minimum C(α) in step vi) is the optimum and defines the best-fitting LFF model.
We would like to point out some requirements which have to be taken into account in order to revise the C03 method and apply it to find the best linear MHS model for quiet-Sun regions.
• C03 shows only one coronal loop and its two footpoint areas that could be easily identified. The identification was done by hand, not automatically. Quiet-Sun areas in data from SDO/AIA 193 Å shown in Figure 1 typically show a bundle of small-scale, well-resolved loops thanks to the higher spatial resolution of AIA. For instance, the loop analysed in C03, if recorded in the 171 Å channel of the SDO/AIA, would appear as a set of many loops. The Yohkoh/SXT data used in C03 have a pixel size of 9.81 (∼7.1 Mm) in comparison to SDO/AIA pixel size of 0.6 (0.43 Mm). • It is not always possible to identify both footpoints of loops. As loops' footpoints join in a relatively small area, they create a fan-type appearance, with the bottom part (loop footpoints) of the fan appearing blurred, i.e. the footpoint of each individual loop often cannot be resolved as one requires higher spatial resolution than the one of the present instrumentation. • One constraint of the C03 method is that the entire loop is well visible and not overlapped by other loops crossing along the line of sight.
• As the MHS model has three free parameters (α, a, κ), scanning the entire parameter space is therefore ineffective.
Step 7 in the C03 method has to be replaced by some other method, e.g., a downhill-simplex minimization, instead of scanning the full 3D-parameter space.
In the following, we tackle the associated problems in two steps. First, we develop a semiautomatic method (Section 5), where we identify one loop and its two footpoint areas by hand. Second, as our main goal, we develop a fully automatic algorithm in the sense that for a given magnetogram and a corresponding extreme-ultraviolet (EUV) image, the optimum LMHS solutions are computed automatically (Section 6).

Observational Material
The data used here were taken with SDO/AIA (Lemen et al., 2012) and the Helioseismic Magnetic Imager (HMI)  on board SDO. We used AIA data obtained in the coronal EUV 193 Å channel dominated by emission from Fe XII (log T (K) ∼ 6.2). The temperature response of these channels is given in Mou et al. (2018). The AIA EUV data have a 12 s cadence and 0.6 × 0.6 pixel size. We used HMI line-of-sight magnetograms taken at 45 s cadence. Both datasets were extracted from a 48 hr time series that starts at 00:00 UT on 15 September 2019 and ends on 16 September 2019 at 23:59 UT. This dataset is used in another two forthcoming studies by Madjarska et al. (in prep). The AIA 193 images after the standard data reduction were processed by binning 3 consecutive images in order to increase the signal-to-noise ratio. The HMI data have a 0.5 × 0.5 pixel size but were rescaled to the AIA pixel size of 0.6 × 0.6 , which corresponds to 432 km × 432 km. To achieve a higher signal-to-noise ratio, the HMI images were binned over 8 frames. The choice of 8 frames comes from the cadence of 6 min that needed to be achieved for the purposes of the study by Madjarska et al. (in prep). The AIA 1600 Å images were used to align the HMI data with those from the AIA 193. All images were derogated to 00:00 UT on 16 September 2019. The AIA 193 images after the standard processing were further processed with used the Multi-scale Gaussian Normalization (MGN) code by Morgan and Druckmüller (2014). The resulting images present better enhanced fine structures.

A Semi-Automatic Approach
As a first step, we try to find the best fitting LMHS solution for individual loops and thereby evaluate the effect of the different free parameters. We follow similar steps as prescribed in Carcedo et al. (2003, Section 2) summarized in Section 3. We investigate four loops visible in the AIA 193 images, one in CBP-1, two in CBP-2 (CBP-2a and CBP-2b in Table 1) and one in CBP-3. As a first step, the footpoint areas of each loop are selected by hand. They are marked with white circles in Figures 2-5. For each loop, we investigate six different cases by using the original AIA images (marked 'orig') and the processed ones (marked 'prep').
The counter in the first column of Table 1 is used to label the panels in the Figures 2-5 and 7, 8, 9, 10. Next, we compute a number of 1000 loops from both footpoint areas, which are magnetically connected to the opposite footpoint areas. The starting points are randomly selected within both areas and only closed magnetic-field lines, which magnetically connect these two footpoint areas are considered for further analysis. Each of these selected field lines is then uncurled. The uncurling process represents the extraction of the emissivity    the quiet-Sun AIA 193 images of the CBPs display several close-by loops. Thus, one should not use too many points in the perpendicular direction to avoid interference in the uncurled images by several loops. The Gaussfit provides quantitative information on how far the field lines (which are per construction in the centre of the uncurled image) are apart from the highest emissivity. This quantity is then averaged over all points along the loop and is  defined as a quantitative measure for a field line labelled with i by Optimized linear MHS parameters with a semi-automatic approach. The first column counts the investigated cases (the same counter 01 -24 is used in Figures 2, 3, 4, 5) and the cases' names. For example, CBP-1 we investigated one loop, for CBP-2 two loops (i.e. cases CBP-2a and CBP-2b) and for CBP-3 one loop. We used both the original images, denoted as 'orig' and the processed one as 'prep'. The second column shows the value of n as used in the definition for L in Equation 11. The third column L Pot correspond to the optimum for a potential-field model and column four L LFF for a force-free one, with the optimized value α LFF shown in column five. Finally, column six shows L MHS and the corresponding optimized free parameters α, a, κ are given in the last three columns 7 -9. where N max,j is derived from the Gaussfit at each of the m positions along the field line. We use a Gaussfit in the form where x is the direction perpendicular to the chosen magnetic-field line, A 0 is the height of the Gaussian, A 2 is the width of the Gaussian, and N max is the location of the center of the Gaussian. By construction, the center of the magnetic-field line is at x = 0 and consequently N 2 max is the square of the distance between the magnetic-field line and the center of the EUV loop. Equation 7 generalizes the linear force-free expression introduced in C03 for MHS equilibria, and for a given parameter set α, a, and κ the best fitting field line is defined as C MHS (α, a, κ) = Min(C 2 i (α, a, κ)).
If only one parameter has to be fitted like in C03, one can easily scan the entire parameter space and show how C varies as a function of α. For the three free parameters needed in LMHS, this would be inefficient and we, therefore, use a downhill-simplex minimization procedure (for details see Nelder and Mead, 1965) to find the minimum of C MHS (α, a, κ). Another point is that the C03 method was applied to an image with only one well-visible loop due to the lower resolution of the used image. Our new method is aimed at highresolution images where multiple loops are resolved. The code will automatically find these loops. As coronal loops are brighter than the background, we define a functional where I uncurled is the average brightness of the uncurled image. For convenience, we normalize I uncurled with the average brightness of the image. For a power of n = 0 the average intensity of the uncurled image I uncurled is not considered at all and L MHS = C MHS . For a power of n = 1 and n = 2 the functional L MHS becomes as smaller as larger I uncurled is. This is important for the automatic method, where the coronal loops are not pre-selected, as it gives more weight to fitting the magnetic-field line to brighter loops. For comparison, we also compute L Pot , which corresponds to the best-fitting field line of a potential field, which does of course not have any free parameters. Furthermore, we estimate L LFF for a linear force-free model, which has only α LFF as a free parameter. Table 1 presents the results of this semi-automatic approach for the 24 cases. We also investigated for each case how the MGN processing of the AIA EUV images influences the result.

CBP-1 with Semi-Automatic Algorithm
The example CBP-1 is shown in Table 1 enumerated 01 -06. In all cases, the minimum of L MHS and L LFF is about two orders of magnitude smaller compared with a potential-field line. The MHS results are only slightly better compared with a linear force-free model. The optimum magnetic-field lines for the MHS model, as shown in Figure 2, outline well the CBP-1 loop and hardly differ, as can be seen in panels 01 -06 of Figure 2. The scatter of the optimum α is small and hardly differs between the LFF and MHS models. The values of L are only very little improved from LFF to MHS. We observe a rather large scatter of the force parameter a, which, however, does not lead to differences in the shape of the magnetic-field lines. The observed loop is bright and respectively the values of L decrease for n = 1 and n = 2, when the average brightness of the uncurled images (shown embedded in all panels of the figure) is used in the nominator of L. Using the processed images (examples 02, 04, and 06) lead to lower values for the optimum MHS and LFF field lines, but the algorithm finds about the same loop as well as values of α.

CBP-2a with Semi-Automatic Algorithm
The example CBP-2a is given in Table 1 enumerated 07 -12 and the optimum field lines for the MHS model are shown in Figure 3. Similar to the first example, the optimum magnetic-field line hardly changes for the different parameter sets and the scatter of α is also small. As |α| is relatively small (below unity) the LFF and MHS models do not improve by orders of magnitude over the potential-field approach, but still, they are better by a factor of two. Alike in the first example, the MHS parameters a and κ show a significant scatter. For example 11, the algorithm even finds a = 0, which corresponds to a force-free model.

CBP-2b with Semi-Automatic Algorithm
The example CBP-2b is presented in Table 1 enumerated 13 -18 and the optimum field lines for the MHS model are shown in Figure 4. A visual inspection of the figure directly shows that the optimum field line for the original AIA images (panels 13, 15, and 17) is very different from the optimum field line in the processed images (panels 14, 16, and 18). Not surprisingly this is confirmed in Table 1 and we find optimum values of α of above 4 for the original images, leading to strongly S-shaped field lines in these cases. For the processed images we find an optimum value of α of about 2.5. In both cases, the scatter of α is small for different values of the power n and hardly differs between the MHS and LFF models. Similar to the previous example cases, the MHS parameters a and κ vary significantly.

CBP-3 with Semi-Automatic Algorithm
The example CBP-3 is shown in Table 1 enumerated 19 -24 and the optimum field lines for the MHS model are shown in Figure 5. Again, a visual inspection of the figure directly shows that two very different field lines have been found and the optimum loops correspond to very distinct coronal loops visible in the AIA 193 image. For the processed cases (20, 22 and 24), the same loop was found, and surprisingly, this loop was also found in the original images (case 19 for n = 0 meaning that the brightness of the uncurled image is not considered in L). For these cases the optimum value of α is about −2 and hardly scatters between the LFF and MHS models and also not between the cases. For the original AIA 193 image using n = 1 and n = 2 the algorithm finds a different field line than for n = 0. It is well visible in panels 21 and 23 that a bright loop was chosen. As mentioned earlier for n = 0 the loop brightness is not considered and for the processed images this loop is not brighter than the others. The optimum value of α is about 1.25. This example has another surprise for us. The MHS values a and κ show a low scatter between all six cases (19 -24).

Discussion and Conclusions from the Semi-Automatic Method
Taking all 24 examples, we find that the optimum L value for the MHS model is on average a factor of 168 better than a potential-field model, but only a factor of 1.1 better than a linear force-free model. The optimum value of α is almost the same for MHS and LFF, with an average absolute difference of only 0.06. A possible reason is that close to the footpoints of the loops, the horizontal currents, which are controlled by a have also a field-aligned component. It seems that in most cases the LFF model is already a reasonable approach to finding magnetic loops fitting the structures observed in the coronal images. While on average the absolute difference is only 0.06, for some examples the LMHS model fits the structures significantly better. Any changes (or improvements) of the LMHS model over the LFF model can only happen in the lowest solar atmosphere, where the term a exp(−κz) has a significant contribution, which happens at the height of about z = 2 Mm. Above this, the LMHS model becomes linear force-free. Consequently, improvement of LMHS over LFF can happen only close to the footpoints of the EUV loop (in the upper photosphere and chromosphere) and not at coronal heights. An advantage of the LMHS model is the potential it has to incorporate additional measurements and observations, e.g., horizontal magneticfield measurements and chromospheric observations. The rather large scatter of the MHS parameters a and κ in most cases (except for CBP-3) led us to the question of whether it is possible to fit α from the coronal structures, but constrain a and κ from other measurements, e.g. vector magnetograms as done in Wiegelmann et al. (2017). Vector magnetograms are not available for the quiet Sun investigated here, but generally, such magnetograms are available in active regions.
To investigate the effect of such constraints, we performed MHS extrapolations with fixed values of a = 0.8 and κ = 0.5, which corresponds to quite large forces in the photosphere and a forced layer (drop by 1/e) of 0.9 Mm. On average for this kind of constrained MHS model the L values are a factor of 1.4 worse than the original MHS model and we find a scatter of the absolute difference of the value of α of 0.7. The differences are, however, enormous between the cases. For cases 13, 15, 17, 20, 21 and 23, the difference in α is larger than 2 with a maximum difference of 3.4 for examples 21 and 23. The reason is that the constrained MHS code fitted a different coronal structure than the original MHS and LFF models. We provide images to show this in Appendix 7. In Figures 7-10 we show the obtained magnetic-field lines from the fully optimized MHS model with a thick black line, from the constrained MHS model with a red line, from the potential-field model with a white line and from the linear force-free field model with a green line. If we leave these 6 troublesome cases out, the difference in α for the remaining 18 examples reduces to 0.05. A visual inspection of Figures 7-10 shows also that (except for panels 13, 15, 17, 20, 21, and 23) the loops from the MHS, constrained MHS, and LFF models agree well. In most cases, the potential-field lines (white) are very different from the other models. An exception is CBP-2a, shown in Figure 8, where even a potential field matches the coronal structures quite reasonably. This is because the value of |α| (see Table 1) for this case is low and thereby almost potential.
We conclude that in general, there is a potential for constraining the MHS model by other measurements, but one has to be cautious and further investigations are required. One might also think to do full optimization with all three MHS parameters, but using two or more images, e.g., in addition to the coronal image a chromospheric or transition-region image could be used.

A Fully Automatic Approach
The semi-automatic approach (see Section 5) might be a good choice if one wants to find the optimal parameters for a particular loop as seen in EUV or X-ray images. As mentioned above, in present high-resolution EUV images (SDO/AIA and the Solar Orbiter Extremeultraviolet Imager (SoLO/EUI)) one typically finds several loops in one loop system (CBPs or active regions). Thus, we aim (for co-temporal magnetogram and EUV images) to develop a fully automatic method, i.e. use a line-of-sight magnetogram and a EUV image as input and compute the optimum MHS parameters α, a, and κ, without or with very limited human guidance. We develop the following approach: i) Identify magnetic elements in the magnetogram; ii) Identify pairs of magnetically connected elements; iii) Optimize the free LMHS parameters by minimizing L.
We explain these steps with the three CBP examples in the following. Table 2 Fully automatically investigated cases. The first column gives the CBP name and the magnetic elements that have been employed (first the positive, then the negative element). The row of stars for two cases means that the elements are not connected by a potential-field model. Not shown are results for element pairs having no magnetic connectivity for whatever model. The other columns have the same meaning as in Table 1.

Identify Magnetic Elements in the Magnetogram
Firstly, we automatically identify magnetic elements in the magnetogram. These magnetic elements will then serve as footpoint areas for the field-line integration. To find the magnetic elements we follow the procedure developed in Wiegelmann et al. (2013). Magnetic elements are a number of pixels of the same polarity which are connected. Our goal is to identify strong magnetic elements. As the lower limit for B z we choose 20% of the maximum field strength of the selected quiet-Sun magnetogram. For the three loop systems (CBPs), the estimated limits are 72 G, 160 G and 30 G for CBP-1, CBP-2 and CBP-3, respectively. For the footpoint areas of the loops we should allow that the footpoints of the field lines are also located in the vicinity of the magnetic elements and not only in the strongest field regions. We chose these footpoint areas to have their center at the center of gravity of the magnetic element. The size of the areas is computed by the area of the corresponding magnetic element plus a small region around the element (here 5 pixels). We show the identified magnetic elements for the three examples in Figure 6, panels a, c and e, respectively.

Identify Pairs of Magnetically Connected Elements
After having identified a number of magnetic elements, the next step is to find which elements (of opposite magnetic flux) are magnetically connected with a significant amount of shared magnetic flux. Finding this very precisely is certainly premature because we still do not have a magnetic-field model at this step. One possibility is to apply a potentialfield model and a small number of MHS solutions with different reasonable parameters. The amount of magnetic flux shared by two elements can change depending on the chosen parameters and it is even possible that two magnetic elements are connected with a certain parameter set and not with another. As we will see below, this does indeed hap-pen in two of our example cases. If the total number of element pairs is small, we can skip this step and run our algorithm for all possible element pairs. We apply this approach here because the number of elements, as shown in Figure 6, for our examples is reasonably small.
After finding one or several magnetic-element pairs with shared magnetic flux by fieldline connectivity of these elements, the next step is to evaluate how well these field lines agree with loops seen in the EUV images. This is assessed by the quantities C or L as defined in Section 5. Using C directly, as done in the semi-analytic approach, where wellvisible loops were chosen by hand is, however, not suitable here. The fact that two elements are connected with a magnetic loop does not necessarily mean that this field line outlines a loop visible in EUV. In unfortunate cases, the field lines would be in relatively dark and homogeneous background emission. This would result in low values of C because naturally, the gradient in brightness would also be very small. We want, however, to find parameters for loops, which show an enhanced emission. We, therefore, use the functional L MHS as defined in Equation 11. Here we consider only the case with power n = 1 and the original images.

Optimize the Free LMHS Parameters
In Table 2 we show the results for the automatic method and discuss them below for each of the three CBP examples.

CBP-1 with Fully Automatic Algorithm
The first three rows in Table 2 and Figure 6b show the optimum values of L for the first example. Surprisingly, L is already quite low for a potential-field model for network elements 1 and 2, and becomes much smaller for LFF and MHS. L LFF and L MHS are about a factor of 4 lower compared with what we found in the semi-automatic approach in Section 5 and the value of α is somewhat higher. The best-fitting field line is shown in red in Figure 6b. A visual comparison with the semi-automatic cases in Figure 2 depicts a somewhat shorter loop. For network elements 1 and 3, connected with a green line in Figure 6b, the value of α is about double what we found in the semi-automatic approach and the L values are significantly higher, which can be confirmed visually. The green field line does not match well the visible AIA 193 loop. Finally, a visual inspection of Figure 6b also shows no bright loop structure connecting elements 1 and 4, which is quantitatively confirmed by our algorithm, which outputs a very high value of the functional L MHS > 100.

CBP-2 with Fully Automatic Algorithm
Rows 4 -9 in Table 2 and Figure 6d show the optimum values of L for CBP-2. The lowest value of L for all three models is found for the element pair 4 and 1. The connecting field line is shown with a blue line in Figure 6d. The optimum field is almost potential (very small value of |α| and a is zero). Visually, the blue field line agrees well with a visible loop structure and this case corresponds to the semi-automatic case CBP-2a shown in Figure 3. The semi-automatic approach found somewhat lower values of the functional L and a somewhat stronger deviation from a potential-field model. The red field line connects the two strongest magnetic elements 2 and 1 in Figure 6d and it seems to be almost parallel to the blue field line and it is slightly non-potential with a small negative value of α. The green field line connecting the magnetic elements 2 and 3 and the yellow linking 5 and 1 are almost parallel for a large part and outline nicely the S-shaped structures seen in the AIA 193 CBP-2 image. Both field lines are strongly non-potential with α larger than 4. The magnetic elements 4 and 3 are not magnetically connected within a potential-field model, which is marked by stars in Table 2 and as a black field line in Figure 6d. A white line corresponds to a different structure similar to one of the lines identified with the semi-automatic method for the processed image only, see semi-automatic case CBP-2b in Figure 4. It is surprising that the automatic method finds field lines for all these loops, while the semi-automatic one found some for the original and others for the processed image. In this example, all 6 possible pairs of positive and negative elements are magnetically connected in the LFF and MHS models and the corresponding magnetic-field lines do outline coronal loops seen in the AIA image.

CBP-3 with Fully Automatic Algorithm
Table 2 and Figure 6f show the optimum values of L for the third example, CBP-3. The red field line connects the two strongest element pair 2 and 1. In comparison to other cases, the MHS model depicts a significantly smaller value of L than the LFF model. L MHS is around 90, which is larger than in the other examples. CBP-3 was already our trouble case in the semi-automatic approach (the loops are very faint) in Figure 5 and the same is true here. We find some other element pairs which are magnetically connected, but they do not outline any coronal structure seen in the AIA 193 image, which is also quantitatively found by our algorithm with L values above 100. The element pairs 7 and 1, and 7 and 5, connected with green and blue field lines in Figure 6f outline bright coronal loop structures, but a visual inspection suggests that the footpoints have not been identified correctly. Only part of the loop structure is outlined by the magnetic-field line.

Discussion and Conclusions from the Fully Automatic Algorithm
In the fully automatic code, a magnetogram and an EUV image are given and without any human intervention the code automatically finds the optimum parameters for an MHS model. Neither the number of field lines nor the number of coronal loops which have to be identified has to be predefined as the code finds this automatically. The output of the code is a number of optimized field lines, which match well with coronal loops detected by their enhanced emission in EUV images. These field lines are then automatically overplotted onto the image. For each of the identified bright-loop structures the magnetic elements at both footpoints are given and the optimum set of the free parameters α, a, and κ are calculated. What has to be predefined in the code once is the threshold value for the magnetic element finder (here: 20% of the maximum absolute field strength), which is the same for all examples, the size of the footpoints area (here: the area of the magnetic element plus five pixels in the vicinity), and the power n of I uncurled in Equation 11. Also the number of vertical pixels in the uncurled images (here: seven) was predefined the same for all cases, also the number of test field lines (here: one thousand). The values used here work well for a combination of SDO/HMI magnetograms and SDO/AIA 193 images. It will probably be necessary to adjust them for data from higher resolution instrumentation. For a very large FOV, which includes many magnetic elements, it might be numerically too expensive to let the code run for each possible magnetic element pair, and elements which are magnetically connected should be identified in advance. Finally, there is (as in the semi-automatic approach) also the option to constrain further the free model parameters by other observations, e.g., vector magnetograms in active regions or chromospheric and transition-region images.

Summary, Conclusions, and Outlook
We present here the further development of a semi-automatic method first proposed by Carcedo et al. (2003) for the comparison of magnetic-field lines obtained from a linear force-free model with loop structures visible in Yohkoh/SXT images. First, we adapt the semi-automatic method of Carcedo et al. (2003) for a linear magneto-hydro-static (LMHS) model (Section 5). Second, as our main goal, we develop a fully automatic method which computes optimum LMHS solutions for a given magnetogram and a corresponding extremeultraviolet (EUV) image (Section 6). While the linear-force-free model used by Carcedo et al. (2003) has only one free parameter, α, our LMHS model has three free parameters, α, a, and κ. Thus to make the implementation of the LMHS model and make it computationally feasible while scanning the whole parameter space, we use a downhill-simplex method to minimize a functional. This functional, L (see Section 5), tells us, how well the magnetic field structure agrees with the image. We successfully obtained magnetic-field lines that match the observed loop structure in AIA 193 images. In all cases the LMHS solution (if optimized towards a particular structure) agrees much better than the potential-field solution and in most cases slightly better than with an LFF solution. A shortcoming of all linear magnetic-field models (linear force-free and linear magneto-hydro-static) is that usually, the optimum free model parameters are different for different structures. In active regions, one value of α usually does not fit all visible EUV loops well, but α changes from one field line to another. We find the same here for the small-scale loops in the quiet Sun. The optimum MHS parameter set α, a, and κ changes for different structures. In both cases, this demonstrates the nonlinearity of the Sun. For the solar coronal magnetic fields investigated here, the linear magnetostatic model provided only a small improvement over linear force-free fields. As soon as we reach coronal heights and a low plasma beta, the forces available in MHS vanish and the configuration becomes force-free (see  for nonlinear MHS computations in an active region and a comparison/combination with NLFFF). The additional power MHS models have over force-free models is that they take into account that the solar photosphere and part of the chromosphere have a finite plasma beta and consequently the assumption of a vanishing Lorentz-force is not justified here. MHS-models do have this (rather thin, about 2 Mm) forced layer and become approximately force-free in the corona above. Consequently, the physics that the magnetic field transposes from forced photospheric and chromospheric fields to a force-free corona is also included in linear MHS.
We plan to explore the algorithm on higher resolution imaging and magnetic-field data (e.g. Daniel K. Inouye Solar Telescope, Solar Orbiter Extreme-Ultraviolet Imager and Polarimetric and Helioseismic Imager) and measurements of the chromosphere (from ground and space) and transition region (Interface Region Imaging Spectrograph). In active regions, where vector magnetograms are available, we will employ a nonlinear MHS model. Data Availability The SDO/AIA and HMI data are publicly available from various sources.

Competing interests The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.