Analytical Three-dimensional Magnetohydrostatic Equilibrium Solutions for Magnetic Field Extrapolation Allowing a Transition from Non-force-free to Force-free Magnetic Fields

For the extrapolation of magnetic fields into the solar corona from measurements taken in the photosphere (or chromosphere) force-free magnetic fields are typically used. This does not take into account that the lower layers of the solar atmosphere are not force-free. While some numerical extrapolation methods using magnetohydrostatic magnetic fields have been suggested, a complementary and numerically comparatively cheap method is to use analytical magnetohydrostatic equilibria to extrapolate the magnetic field. In this paper, we present a new family of solutions for a special class of analytical three-dimensional magnetohydrostatic equilibria, which can be of use for such magnetic field extrapolation. The new solutions allow for the more flexible modelling of a transition from non-force-free to (linear) force-free magnetic fields. In particular, the height and width of the region where this transition takes place can be specified by choosing appropriate model parameters.


Introduction
Modelling the magnetic field in the solar atmosphere is of great importance for our interpretation of many of solar observations, in particular in the solar corona T. Neukirch tn3@st-andrews.ac.uk T. Wiegelmann wiegelmann@mps.mpg.de 1 (e.g. Wiegelmann, Thalmann, and Solanki, 2014). Because coronal magnetic fields cannot be measured routinely with the required accuracy, extrapolation methods with photospheric magnetic field measurements as boundary conditions are normally used, usually assuming that the magnetic field is force-free (see e.g. recent reviews by Wiegelmann and Sakurai, 2012;Régnier, 2013). In recent years measurements of the magnetic field in the chromosphere have also become available (e.g. Harvey, 2012). An overview of measurements of photospheric and chromospheric magnetic fields can, for example, be found in the paper by Lagg et al. (2017).
While the assumption of force-free magnetic fields is well satisfied in large parts of the solar corona due to the low plasma β, the lower parts of the solar atmosphere (chromosphere and photosphere) can in general not be considered to be force-free (e.g. Metcalf et al., 1995;De Rosa et al., 2009;Wiegelmann, Thalmann, and Solanki, 2014) and hence magnetohydrostatic (MHS) models, including pressure and gravity, would be more appropriate for these regions. Developing numerical extrapolation methods for the magnetostatic case has, for example, been attempted in Wiegelmann and Neukirch (2006) (including pressure only), Gilchrist and Wheatland (2013), and Zhu and Wiegelmann (2018).
These numerical approaches are usually computationally expensive. Hence, a complementary approach for including pressure and gravity forces in magnetic extrapolation, which is computationally relatively cheap, would be to use analytical three-dimensional MHS equilibrium solutions. Obviously, just as for 3D force-free fields only a limited number of analytical 3D MHS solutions suitable for magnetic field extrapolation are known, with the known 3D MHS solutions useful for extrapolation purposed being comparable in status to 3D linear forcefree solutions. We emphasize that in order to find analytical solutions to the MHS equations in 3D one has to make a number of assumptions, which may limit the applicability of the method to some extent. Hence, this approach using analytical MHS equilibria has to be regarded as an alternative method which allows one to get a reasonably fast extrapolation method including a non-forcefree part of the solar atmosphere, but not as a replacement for the numerical approaches mentioned above.
Subsets of these three-dimensional MHS solutions have been used for modelling both global solar magnetic field models, using spherical coordinates (e.g. Bagenal and Gibson, 1991;Gibson and Bagenal, 1995;Gibson, Bagenal, and Low, 1996;Hoeksema, 1993, 1994;Zhao, Hoeksema, and Scherrer, 2000;Ruan et al., 2008), and local solar coronal structures, using Cartesian or cylindrical coordinates (e.g. Aulanier et al., 1999Aulanier et al., , 1998Petrie, 2000;Gent et al., 2013;Gent, Fedun, and Erdélyi, 2014;Wiegelmann et al., 2015Wiegelmann et al., , 2017. Other applications include, for example, models of the magnetic fields of stars (e.g. MacTaggart et al., 2016) and their interaction with exoplanets (e.g. Lanza, 2008Lanza, , 2009. In this paper, we shall focus on 3D MHS solutions in Cartesian geometry, i.e. assuming a constant gravitational force, which we take to point in the negative z-direction (hence the coordinate z has the meaning of height above the photosphere from now on).
If we follow the argument that small plasma-β should imply nearly forcefree magnetic fields (and vice versa), there should be a marked transition from non-force-free to force-free fields with increasing height when moving from the photosphere through the chromosphere into the corona. Of the currently known solutions the one which comes closest to showing this feature has been suggested by Low (1991Low ( , 1992 and has an exponential height profile of the non-force-free current density. This solution has been used repeatedly for modelling purposes (see e.g. Aulanier et al., 1998Aulanier et al., , 1999Wiegelmann et al., 2015Wiegelmann et al., , 2017. It is the aim of this paper to provide another set of 3D MHS solutions whose non-force-free current density have a more flexible dependence on height. A general problem which arises in the use of this class of 3D MHS solutions for extrapolation purposes is that one has to be careful to avoid generating regions of negative plasma pressure or density (see e.g. Petrie and Neukirch, 2000;Petrie, 2000;Gent et al., 2013). This problem is caused by the mathematical structure of the expressions for the plasma pressure and density, both of which are written as the difference between a positive background pressure and density, and potentially negative terms, depending on the a priori unknown magnetic field solution. In theory, this problem can always be solved by either increasing the background plasma pressure and density or by decreasing the amplitude of the negative terms. In practice, the former often leads to an unrealistically high value of the plasma-β throughout the model domain, whereas the latter can cause the loss of meaningful spatial structures in plasma pressure and density. Being able to control the solution structure better should be of advantage for modelling purposes and may also have a (positive) bearing on the problem of keeping plasma pressure and density positive everywhere.
The structure of the paper is as follows. In Section 2 we briefly summarise the basic theory of the particular class of 3D MHS solutions we use in this paper and then present the calculation leading to the the new set of solutions in Section 3. In Section 4, we present some example solutions and in Section 5 a discussion and conclusions.

Basic Theory
The MHS equations are given by Here B denotes the magnetic field, j the current density, p the plasma pressure, ρ the mass density, Ψ the gravitational potential and µ 0 is the permeability of free space. In this paper we use Cartesian geometry with a constant gravitational force pointing in the negative z-direction, i.e. Ψ = gz with g being the constant gravitational acceleration. The general theory for this case was first developed by Low (1991Low ( , 1992. Later, Neukirch and Rastätter (1999) presented a somewhat simpler, albeit equivalent, formulation and we will follow their approach in the following brief summary. The main assumption made is that the current density can be written in the form It can be shown (Low, 1991;Neukirch and Rastätter, 1999) that F has to be a function of B z and z. The form for F suggested by Low (1991Low ( , 1992 was i.e. linear in B z , but with an arbitrary function f (z). Clearly, the function F is responsible for the non-force-free, i.e. perpendicular, part of the current density in this approach, although it should be noted that it will generally also contribute to the parallel part of the current density. Hence, the choice of the free function f (z) influences the dependence of the non-force-free part of the current density with height. Choosing F to be a linear in B z leads to a linear partial differential equation for B z (e.g. Low, 1991). Alternatively, representing the magnetic field B in the form (e.g. Nakagawa and Raadu, 1972) one can show that (Neukirch and Rastätter, 1999) where is the Laplace operator, and is the two-dimensional Laplace Operator in x and y. Additionally, one finds that and hence that The plasma pressure and the plasma density are given by the expressions where g is the constant gravitational acceleration. The plasma temperature can be determine from the plasma density and pressure via an equation of state, for example, the ideal gas law withμ the mean atomic weight of the plasma and k B the Boltzmann constant.

A New Family of Solutions
Solutions to Equation 7 have been found for several different choices of the function f (z), using either separation of variables (see e.g. Low, 1991) or the Green's function method (e.g. Petrie and Neukirch, 2000). We will use separation of variables in this paper. Leaving f (z) unspecified for the time being, separation of variables leads to in Cartesian coordinates x, y, z with the equation forΦ being where k 2 = k 2 x + k 2 y . In Equation 18 we have suppressed the parametric dependence ofΦ on k x and k y in the Cartesian coordinate case (or on k and m in the cylindrical coordinate case). At this point it is also convenient to normalise all coordinates by a typical length scale L, as well as regarding α as the normalised αL and k as kL from now on.
In practical applications, for example using potential or linear force-free magnetic fields to extrapolate photospheric field measurements, different versions of periodic boundary conditions are often imposed (e.g. Seehafer, 1978;Alissandrakis, 1981;Otto, Büchner, and Nikutowski, 2007). If periodic boundary conditions in x and y are imposed k x , k y take on discrete values and the integrals in Equation 17 are replaced by infinite sums (which are then truncated in applications to observational data). With such boundary conditions Fast Fourier Transformation (FFT) techniques can be used to implement the boundary conditions very effectively.
As pointed out by Neukirch (1995) Equation 18 is analogous to a one-dimensional Schrödinger equation with the potential −k 2 f (z) with energy eigenvalue −(α 2 − k 2 ) (or alternatively −k 2 [f (z) − 1] with energy eigenvalue −α 2 ). Thus solutions to Equation 18 are known for a large class of functions f (z), but not all of these will be of use for modelling solar magnetic fields. Examples of functions f (z) which have been used for modelling solar magnetic fields are f (z) = f 0 =constant, for which the solutions forΦ are exponential functions (e.g. Neukirch, 1997a;Petrie, 2000;Petrie and Neukirch, 2000), and f (z) = a exp(−κz), with Bessel functions with an argument that is proportional to exp(−κz/2) as solutions (e.g. Low, 1991Low, , 1992Aulanier et al., 1998Aulanier et al., , 1999. In this paper we aim to find solutions to Equations 7 or 18, respectively, which show a transition from non-force-free behaviour to (linear) force-free behaviour as the height z increases from the photosphere through the chromosphere into the corona. Hence, we would like f (z) to possibly decrease from non-zero values to a value close to zero when reaching the corona. While it is also in principle possible to model such a behaviour with the exponential f (z) mentioned above, the transition is only controlled by the parameter κ. Because we would like to be able to control both the height at which the transition from non-force-free to force-free fields happens as well as the range in height over which this happens we here suggest to use the function with a, b, z 0 and ∆z real parameters. The parameters a and b are dimensionless and we assume a > 0 and 1 ≥ b > 0 in this paper. At z = 0 we have The parameters a and b control the overall magnitude of f in the regimes The parameter z 0 controls where this transition takes place and the parameter ∆z is the length scale over which the transition happens. We see that f → 0 for (z − z 0 )/∆z 0 if b = 1 is chosen. This means that any perpendicular electric currents go to zero above the region in which the transition happens and the magnetic field becomes a (linear) force-free field.
We remark that the form of f (z) in Equation 19 includes the possibility of f (z) = f 0 = constant, by choosing b = 0. As discussed above, for this special case the z-dependence of the solutions is given by exponential functions (e.g. . Also shown are the two limiting cases described in the main text. Neukirch, 1997a;Petrie, 2000;Petrie and Neukirch, 2000). In principle, Equation 19 also includes the case of an exponentially decaying function of z (e.g. Low, 1991) by choosing b = 1 and letting As stated above, in this case the z-dependence is given by Bessel functions with exponential arguments of the form exp(−κz/2), with κ/2 = 1/(∆z) in our limit (see e.g. Low, 1991). However, as one can see from the expression inside the limit, this case combines two of the parameters of f into one parameterā and hence has less flexibility in applications. Although, as discussed before, the model with an exponential f (z) also allows for a transition from non-force-free to force-free fields one can only control the length scale over which this happens, but one does not have an additional specific height at which this happens such as z 0 in the f considered in this paper. We show example plots of the three different forms of f (z) discussed above in Figure 1.
Substituting the function f (z) given in Equation 19 into Equation 18 we obtain: Making use of the coordinate transformation with 0 < η < 1, we can transform Equation 22 into the differential equation where withk = k∆z andᾱ = α∆z. Equation 24 can be solved using the hypergeometric function 2 F 1 (a, b, c; z) (see e.g. Abramowitz and Stegun, 1965;Olver et al., 2010) in the formΦ with γ = √ C 2 , δ = √ C 1 , and A, B constants to be determined by boundary conditions at z = z min (here we take z min = 0) and z = z max . For example, if we want the solution to tend to zero as z → ∞ (η → 0), we have to choose B = 0. The other coefficient A will be determined by the boundary condition at z = z min = 0 in the same way as in the potential or linear force-free magnetic field case. The only difference is the functional dependence of the individual modes on z.
We notice that for given values of a, b andᾱ the parameters C 1 and C 2 become negative for small k (in particular C 1 < 0 ifk 2 <ᾱ 2 /[1 − a(1 − b)] and C 2 < 0 whenk 2 <ᾱ 2 /[1 − a(1 + b)]). This implies that γ, δ or both become imaginary. The change in nature of the solution is similar to the change from exponential to trigonometric in the linear force-free case and while this can in principle be incorporated into a Green's function approach (for a discussion see e.g. Chiu and Hilton, 1977;Wheatland, 1999;Petrie and Lothian, 2003;Priest, 2014) one usually has to impose restrictions on the values of a, b andᾱ in the case with periodic boundary conditions and discrete Fourier modes. As in the linear force-free case the specific bounds are determined by the smallest value of the wave vector k 2 . We will discuss an example in section 4.

Illustrative Examples
Instead of using an arbitrary magnetogram we investigate the effects of the solution parameters on the structure of the magnetic field and the plasma pressure, density and temperature, by using a relatively simple, doubly periodic example which allows us to control the shape of the bottom boundary conditions and hence enables us to undertake a study of the solutions that highlights the features of the solutions more clearly. To achieve this we will use the following boundary condition for B z at z = 0: This choice is based on a special case of the bivariate von Mises distribution which is used, for example, in directional statistics (see e.g. Mardia and Jupp, 1999) and we combine two of these functions, one with a positive sign and one with a negative sign, to balance the magnetic flux through the lower boundary. Herex = πx/L andỹ = πy/L, with the domain size in the x-and y-directions given by −1 ≤x,ỹ ≤ 1. All other quantities with a tilde that are directly related to length scales are also normalised by L. Furthermore, B 0 is a reference magnetic field value, theκ ij (> 0) values determine the width of the flux distribution, with largerκ ij values making the width to smaller, theμ ij values specify the positions of the maximum or minimum of the magnetic flux distribution. The function I 0 (x) is a modified Bessel function of the first kind (see e.g. Abramowitz and Stegun, 1965;Olver et al., 2010). The denominator is included to normalise the integral of the bivariate von Mises distribution over the x-y box to unity, so that the positive and the negative magnetic flux through the bottom boundary cancel exactly and the total magnetic flux through the bottom boundary vanishes, independently of the values chosen forκ ij andμ ij . The remaining boundary condition imposed is B z → 0 as z → ∞. The function 28 is periodic in the x-and the y-direction and can easily be expanded into a Fourier series, which allows us to find the general solution for these boundary conditions without problems (for mathematical details see Appendix A). We show surface plots of the exact expression 28 and an expansion based on the first ten Fourier modes in both x and y in Figure 2. The parameter values used in this plot areμ x1 =μ y1 = −μ x2 = −μ y2 = 1.2/π ≈ 0.382 and κ x1 = κ x2 = κ y1 = κ y2 = 10. The maximum difference between the exact plot and the ten-mode approximation is of the order of 10 −5 and hence the tenmode approximation is considered to be sufficiently accurate for this choice of parameter values.
For the MHS examples we will be showing we have used the values z 0 = 0.2L, ∆z = 0.1z 0 and b = 1.0, which means that for z z 0 the magnetic field tends towards a potential state (in the case α = 0) or a linear force-free state (if α = 0). By choosing b < 1.0 one could retain a controlled level of MHS behaviour of the magnetic field above z 0 .
We will present solutions for different values of the parameter a below the maximum value for a, which is given by As discussed in Section 3 the condition 29 follows directly from the definition of γ and for given b, α and minimum value of k 2 corresponds to the value of a at which γ becomes imaginary. One can also see from Equation 29 that to have a max > 0 one has to have k 2 min > α 2 which is the condition for linear force-free field modes to drop off exponentially with height.
We present magnetic field line plots for four different parameter combinations in Figures 3 and 4. For reference we present the potential magnetic field for the given boundary conditions (a = 0.0, α = 0.0) in panel (a) and the linear forcefree magnetic field with α = 0.5 (a = 0) in panel (b) of Figures 3 and 4. We compare these two cases with two field line plots for non-zero values of a, one roughly half of a max and the other 0.99 a max . To highlight the region where the field lines change the most between the different cases, we also present the same field line plots on a domain with z ≤ 2z 0 in Figures 5 and 6.
As expected, the main difference between the linear force-free case and the two MHS cases can be seen for z ≤ z 0 . This is particularly obvious in Figure 4, in which a considerable steepening of the field lines is noticeable in the region below z 0 = 0.2. This change in the field line behaviour can be attributed to the value of the smallest γ in the series expansion approaching zero which leads to a less rapid decrease of the lowest order mode with height for z ≤ z 0 . Due to our choice of a relatively small value for ∆z this behaviour changes relatively sharply around z ≈ z 0 and hence no major changes in the magnetic field can be seen at heights above z 0 . The magnetic field above z 0 is therefore basically identical to a linear force-free field.
In Figures 7 and 8 we show the spatial variation of and ∆ρ = 1 g df dz In these two figures ∆p is normalised by B 2 max /µ 0 , where B max is the maximum value of B z at z = 0, and ∆ρ is normalised by B 2 max /(µ 0 gL) Figure 7 shows the variation of ∆p and ∆ρ with height for 0 ≤ z ≤ 2z 0 at the position of the maximum of |B z | on the lower boundary, i.e. x = y =μ x1 ≈ −0.382. We have chosen these x and y-coordinates because one expects the largest variations of ∆p and ∆ρ to happen there. As one can see both ∆p and ∆ρ have negative values and they generally increase with z from their lowest value at z = 0 until approaching zero around z = z 0 , as expected. The increase of the amplitude of ∆p and ∆ρ and their slower increase with z with increasing value of a is also obvious. The latter is caused by the Fourier modes, in particular the lowest order mode, to decrease less fast with height when a increases. Although a detailed analysis is a bit more complicated, one can roughly regard the z dependence of each mode as being given by exp(−2γz) below and exp(−2δz) above z ≈ z 0 . As a → a max we have γ → 0 for the lowest order mode and hence the mode decreases more slowly with z. Figure 4. The same field line plots as in Figure 3, but viewed along the y-direction (i.e. projected onto the x-z-plane). As in Figure 3 we have −1 ≤ x/L ≤ 1 (horizontal direction) and 0 ≤ z/L ≤ 2 (vertical direction). The field lines in the MHS cases (panels (c) and (d)) show clear signs of steepening below z/L = z 0 /L = 0.2 (second tickmark on the vertical z-axis), as expected.
We also note that ∆ρ has a local maximum and minimum just below z 0 . This is caused by the df dz term in ∆ρ, which is largest at z 0 . This term becomes the larger the sharper the gradient of f at z = z 0 is made, i.e. the smaller ∆z becomes (the derivative actually tends to a δ-function in the limit ∆z → 0 in which case the tanh-function tends to a step function). For this reason one should not attempt to make ∆z too small.
To keep the pressure and density positive everywhere we need to make the background pressure p 0 (z) and the corresponding background density ρ 0 (z) = − 1 g dp0 dz positive and larger at every height z than the minimum values of ∆p and ∆ρ for all x, y values at the same height (see e.g. Wiegelmann et al., 2015). As one can see from (c) (d) Figure 5. The same field line plots as in Figure 3, but with 0 ≤ z ≤ 2z 0 to make the difference between the various cases more obvious.
(a) (b) (c) (d) Figure 6. The same field line plots as in Figure 5, but viewed along the y-direction (i.e. projected onto the x-z-plane). needed to be satisfied for z < z 0 in the case b = 1.0 because the contributions to the pressure and the density by ∆p and ∆ρ become arbitrarily small very fast for z > z 0 . In Figure 9 we show a comparison of variation with z of the absolute values of ∆p and ∆ρ (again at x = y =μ x1 ≈ −0.382), for a solution of the family presented in this paper and a solution for f (z) =ā exp(−κz) (Low, 1991). We remark that in order to start with the same pressure at z = 0 one has to set

SOLA
(e) (f ) Figure 8. Surface plots of ∆p and ∆ρ at the heights z = 0 ((a) and (b)), z = z 0 /2 ((c) and (d))and z = z 0 ((e) and (f)), showing the variation of pressure and density in the x-and y-directions. a = 2a. We have also chosen the inverse length scale κ = 1/z 0 , corresponding to the height of the transition of the solutions presented in this paper from non-force-free to force-free. These plots clearly show that while for both sets of solutions the pressure and density deviations decrease rapidly for z > z 0 , the decrease is much more rapid for our solutions, which, as already stated above, (b) Figure 9. Comparison of the variation with height z of the absolute values of the pressure (panel (a)) and density (panel(b)) deviation from a stratified background atmosphere for the solutions presented in this paper (solid lines) and for a solution with an exponential f (z) (Low, 1991). The plots have a logarithmic y-axis to emphasize the difference of the quantities of the solutions, especially the much faster drop-off of the solutions presented in this paper above the transition height z 0 . The parameters chosen for this plot are α = 0.5, a = 0.48, z 0 /L = 0.2 and ∆z = 0.1z 0 for the solid line, whereas we have chosenā = 0.96 and κ = 1/z 0 (see main text for definitions). The x-and y-coordinates of these plots are the position of the maximum value of |Bz| on the lower boundary, x/L = µx/π ≈ −0.38, y/L = µy/π ≈ −0.38.
reduces the need for artificially adjusting the stratified background atmosphere in this region to avoid negative density or pressure values. One could of course argue that by increasing κ, i.e. reducing the length scale over which the exponential f (z) and hence the pressure and density deviations decay, one could in principle achieve a similar effect for the exponential f (z). However, for this f (z) this would also decrease the effect of the non-force-free part of the current density in the regions below z 0 and hence reduce the effect of the non-force-free part of the current density, which somewhat defeats the purpose of using such a magnetic field model for extrapolation in the first place. Because the solutions presented in this paper can guarantee a very rapid transition to a linear force-free magnetic field, but still allow us to control the region z < z 0 separately this problem does not arise in the new family of solutions.

Discussion and Conclusions
We have presented a new family of three-dimensional MHS solutions based on the general theory originally developed by Low and co-workers (e.g. Low, 1985;Bogdan and Low, 1986;Low, 1991Low, , 1992, although we have used the alternative mathematical formulation by Neukirch and Rastätter (1999). The main motivation for trying to find these solutions was that they are of potential importance for analytical non-force-free magnetic field extrapolation methods. The new family of solutions allows for the MHS nature of the equilibria to be limited to a domain below a specific pre-determined height with a smooth transition to a potential or linear force-free solution possible above height. This is achieved by choosing one of the fundamental free functions of the theory in the form of a hyperbolic tangent. It is important to emphasize that this is just one possibility of choosing the parameters for this equilibrium family and that MHS solutions in the full domain are also possible.
While such a transition is also possible with the free function mentioned above chosen in the form a decaying exponential function (e.g. Low, 1991Low, , 1992Aulanier et al., 1998Aulanier et al., , 1999Wiegelmann et al., 2015Wiegelmann et al., , 2017, the new solution family allows much more control over the transition from the MHS domain to the non-MHS domain, in particular regarding the departures of the plasma pressure and density from a stratified atmosphere. This property can be of importance for keeping the plasma pressure and density of the solution positive and should make this equilibrium family interesting for magnetic field extrapolation purposes when a non-force-free layer has to be included. We emphasize that we do not regard these analytical extrapolation methods as replacements for numerical methods for non-force-free magnetic field extrapolation, but as a numerically relatively cheap complementary method, which could be used as an initial "quick look" tool. Obviously, the limitations that arise from having to make a number of strong assumptions to allow analytical progress have to born in mind when applying these methods. withΦ nm (z) the solution of Equation 22 for the wave vector k 2 nm = π 2 L 2 (n 2 + m 2 ). (40) We remark that n or m are allowed to take on the value 0, but not both at the same time. Comparing Equations 33 and 39 one can easily see thatā nm = a nm /(k 2 nmΦnm (0)) etc.