An Optimization Principle for Computing Stationary MHD Equilibria With Solar Wind Flow

In this work we describe a numerical optimization method for computing stationary MHD-equilibria. The newly developed code is based on a nonlinear force-free optimization principle. We apply our code to model the solar corona using synoptic vector magnetograms as boundary condition. Below about two solar radii the plasma $\beta$ and Alfv\'en Mach number $M_A$ are small and the magnetic field configuration of stationary MHD is basically identical to a nonlinear force-free field, whereas higher up in the corona (where $\beta$ and $M_A$ are above unity) plasma and flow effects become important and stationary MHD and force-free configuration deviate significantly. The new method allows the reconstruction of the coronal magnetic field further outwards than with potential field, nonlinear force-free or magneto-static models. This way the model might help to provide the magnetic connectivity for joint observations of remote sensing and in-situ instruments on Solar Orbiter and Parker Solar Probe.


Introduction
Traditionally the global structure of the coronal magnetic field is modelled with the help of source surface potential field models (PFSS, see Schatten, Wilcox, and Ness, 1969). In potential field models electric currents are neglected and they require as photospheric boundary condition only line-of-sight magnetic field measurements. Nevertheless, the effect of a solar wind flow is considered by the introduction of a so called 'source surface' as outer boundary condition. At this artificial surface, usually at 2.5R s , all field lines are assumed to become radial, thereby mimicking the effect of the solar wind. Going beyond this simple potential field approach with non-potential global coronal magnetic field models has been an active research topic for several years with various physical models (see, e.g., Mackay and Yeates, 2012;Wiegelmann, Petrie, and Riley, 2017, for review articles on global coronal magnetic field modelling.). Recently several of the non-potential methods (seven different codes) have been compared in the framework of an ISSI-meeting and the results have been published (Yeates et al., 2018). While all of the models are non-potential and incorporate electric currents, the physical assumptions and computational implementations are different and the models used also different input data. One group of codes are nonlinear force-free extrapolations (3 different implementations have been compared), which do not consider time dependence and plasma effects. As boundary condition these codes require measurements of the photospheric magnetic field vector. Codes for solving the nonlinear force-free equations have been first applied to active regions and different numerical methods and implementations have been intensively compared and evaluated in a number of studies, e.g., Schrijver et al. (2006); Metcalf et al. (2008); Schrijver et al. (2008);De Rosa et al. (2009);DeRosa et al. (2015). Active region and global nonlinear force-free codes finally solve the same nonlinear force-free equations, but do so with different numerical implementations. Well known codes for global nonlinear force-free computations are based on the Grad-Rubin method (see Amari et al., 2013Amari et al., , 2014, an optimization principle (see Wiegelmann, 2007;Tadesse et al., 2014) and force-free electrodynamics (see Contopoulos, Kalapotharakos, and Georgoulis, 2011;Contopoulos, 2013).
Other approaches go beyond the force-free assumption and take effects like plasma forces, flows and time-dependence into account. In the linear magnetostatic approaches (see Bogdan and Low, 1986;Neukirch, 1995) the Lorentz force is compensated by plasma pressure and gravity forces. Coronal equilibria with plasma pressure and steady 3D nonlinear flows have been found in Nickeler et al. (2017). The evolving magnetofrictional method (see Mackay and van Ballegooijen, 2006;Yeates, 2014) solves also the force-free equations, but takes time-dependence into account and requires therefore a time sequence of photospheric magnetograms (radial component only) as boundary condition and incorporates flux transport. MHD codes from different research groups (see, e.g., Mikic and Linker (1994); Mikić et al. (1999) and Feng et al. (2012);Feng (2020)) are capable to solve the full MHD equations.
For global corona models these codes use the radial photospheric field as boundary conditions. It is possible, however, to limit the MHD simulations to the assumption of a zero-β plasma, and to incorporate additional observations, e.g., use results from magnetofrictional simulations (see Yeates et al., 2018, for details).
An interesting point is that only the (static) nonlinear force-free codes make use of the photospheric vector magnetograms, whereas all other methods compared in Yeates et al. (2018) use the radial photospheric field only. As a consequence the study of Yeates et al. (2018) revealed that the nonlinear force-free extrapolations are superior in active regions (where the vector magnetic field measurements are most accurate) while quiet Sun features like filament channels are better modelled by other approaches. These findings stimulated us using vector magnetograms as boundary condition and go beyond the force-free assumption. To do so we develop a new stationary MHD code with field aligned plasma flow, which uses synoptic vector magnetograms as boundary condition.
We organize the paper as follows: In Section 2 we present the stationary MHD-equations and how we aim to solve them with the help of an optimization principle. Section 3 contains an application and evaluation of the newly developed code to a synoptic vector magnetogram observed with SDO/HMI. Finally we draw conclusions and give an outlook for future work in Section 4.

Stationary Compressible MHD
To model the coronal magnetic field and plasma environment we use the equations of stationary compressible ideal MHD.
where B is the magnetic field, p the plasma pressure, ρ the mass density, v the flow velocity, Ψ the gravity potential, µ 0 the permeability of free space, T the temperature, R the gas constant and E the electric field. The stationary MHD equations are given by the force balance Equation 1, the solenoidal condition (Equation 2), mass continuity Equation 3, an energy equation or an equation of state (Equation 4) and ideal Ohm's law (Equation 5). We are aware that nonideal effects as, e.g., caused by turbulence, can lead to a violation of ideal Ohm's law and play an important role for coronal heating and solar wind expansion (see Cranmer et al., 2015, for a review article). For special cases (2.5D, incompressible flow and no gravity) stationary solutions of resistive MHD have been found (see, e.g., Throumoulopoulos, 1998;Throumoulopoulos and Tasso, 2000;Nickeler et al., 2014). Studying turbulence and resistivity is well outside the scope of this work, however. For simplicity we assume an isothermal equation of state in Equation 4, which leads to a linear relation of plasma pressure and density. We replace p by ρRT in Equation 1 to reduce the number of independent quantities. For the highly conducting coronal plasma we have to consider ideal Ohms law, which is for a vanishing electric field satisfied if v × B = ∇f , with an arbitrary scalar function f , which is constant on magnetic field lines. For the particular choice f = 0, which we use here, the plasma flows and magnetic fields are parallel.

Optimization Principle for Stationary MHD
Minimizing a functional L of quadratic terms to compute 3D coronal magnetic field models has been introduced by Wheatland, Sturrock, and Roumeliotis (2000) for computing nonlinear force-free fields. First implementations of this optimization approach have been done in Cartesian geometry to model active regions. A spherical implementation for global force-free optimization has been done in Wiegelmann (2007) and was first applied with synoptic vector magnetograms as boundary condition in Tadesse et al. (2014). A magneto-hydro-static optimization code has been developed in Wiegelmann et al. (2007) for global computations in spherical geometry. Within this work we try to extend this optimization approach towards stationary MHD. We aim to solve the stationary MHD Equations 1-5 by minimizing a functional L which we define as where all terms in the functional have a quadratic form as defined below. This means that the stationary MHD equations are solved when L is zero.
If the term L force vanishes, Equation 1 is satisfied. We applied vector identities to the flow terms to bring them into a similar form as the magnetic terms. The linearity in the isothermal equation of state (Equation 4) has been used to substitute p by ρRT . The other parts of the functional are linked to the solenoidal condition, the steady-state continuity equation, and the condition that the flow is field aligned: In the L angle(B,v) term we use a weighting with the Alfvén Mach number M A to give a stronger weight to regions with strong flow. The B 2 in the nominator of Equation 7 originates from the nonlinear force-free and magnetohydro-static optimization codes. Dividing by this term ensures that sufficient weight is given to the equilibrium in weak field regions and gives the terms L force and L divB the same dimensionality. While in the functional all terms are of quadratic form, it might be convenient for humans to monitor also the angle between magnetic fields and flows in degree.
This formula is similar to the definition of the weighted average angle between magnetic field and electric currents used for evaluating the quality of force-free computations (see Schrijver et al., 2006, for details.).

Testing and Application of the Method
Previous versions of optimization codes have been tested first with known analytical or semi-analytical equilibria, e.g. Low and Lou (1990) for nonlinear force-free active region models and Neukirch (1995) for global magnetostatic equilibria. Unfortunately we are not aware of analytic 3D solutions of stationary MHDequilibria with compressible plasma flow and gravity. We therefore use the code to construct a numerical equilibrium by minimizing the functional (Equation 6) and monitoring the individual terms (Equations 7 -10).

Boundary Conditions and Initial Force-Free Model
If all non-magnetic terms in the stationary MHD Equations 1-5 are neglected, we get as a subclass the nonlinear force-free field equations, which are given as As boundary data we use a synoptic vector magnetogram from SDO/HMI for Carrington rotation 2099 as shown in Figure 1, which has been observed between 13/07/2010 and 09/08/2010. Because of the grid convergence problem at the poles (see Wiegelmann (2007)) and because accurate vector field measurements at poles are still not available, we cut out polar regions and limit our computation to latitudes θ = 20 . . . 160. The spatial grid resolution is 2 degrees. As initial state we compute a nonlinear force-free field up to r = 10R s , which solves the force-free Equations 12 -13 with the global force-free optimization approach as described in (Wiegelmann, 2007;Tadesse et al., 2014). The force-free iteration itself uses a potential field as initial state and the potential field solution is also kept on the theta boundaries. We would like to point out that the initial state for the force-free iteration is not a PFSS model, because that is radially limited to the source surface located at 2.5R s . We use just the decaying mode of a spherical harmonic representation and for the tests done here limit them to l = 12. This corresponds to a special case of the global linear magnetostatic model developed in Bogdan and Low (1986) with α = 0 and a = 0.

Initial Parker Solar Wind Solution
If all magnetic terms in the stationary MHD Equations 1-5 are neglected, we get as a subclass the stationary hydrodynamic equations where the gravitational potential is Ψ = −GM s /r, where G is the gravity constant, M S the solar mass and r the distance from the center of the sun. To initialize the plasma and flow variables we use the spherical symmetric solution of the stationary hydrodynamic Equations 14 -16, which was found by Parker (1958) and describes the solar wind. We use an isothermal solution which can be solved analytically by using the Lambert W function (see Cranmer, 2004, for details). In the isothermal case there are only two free parameters, T and p (or alternatively T and ρ). We use T = 3MK and p is so specified that the average plasma beta at r = 1R s is β = 0.01, whereas the magnetic pressure was computed (and averaged over the sphere) from the initial force-free magnetic field. The sound velocity is c s = 157 km/s and the corresponding critical radius is at r c = 3.84R s . Figure 2 shows several quantities of this spherically symmetric solution as a function of the radius r. Panel a) contains the plasma pressure gradient force (black), gravity force (green) and the flow force (red). At low coronal heights the flow force is small and gravity and pressure forces compensate each other. With increasing distance from the sun the flow becomes more and more important. The flow velocity is shown in panel b). In panel c) we compare the magnetic pressure (green), plasma pressure (black) and solar wind ram pressure (red). At low heights in the low β corona (black line in panel d) the magnetic pressure dominates. This is also the region which is usually computed using the force-free assumption. With increasing height the plasma pressure and kinematic pressure become more important and finally dominate over the magnetic pressure. The red line in panel d) shows the Alfvén Machnumber M A . In the lower corona M A is very small, but increases to values above unity with increasing distance from the sun. In regions with β 1 and M A 1 it is justified to neglect the influence of plasma effects and flows and this is a justification why these regions are usually modelled under the force-free assumption. Higher up in the corona, however, β and M A exceed unity and the force-free assumption is not valid. Therefore plasma and flow effects have to be taken into account, which is done here in the framework of stationary MHD.

Evaluation of the Initial State and Optimization
In the initial state the force-free Equations 12-13 and the hydrodynamic Equations 14-16 are fulfilled separately. Table 1 shows in the line 'NLFFF, Parker' the discretisation errors of the terms L force , L divB and L cont in the initial state. The term L angle(B,v) is not a discretisation error, it indicates how well the magnetic field and the plasma flow are aligned. For a perfect equilibrium with field aligned Table 1. Value of the different terms of the functional L (see Equations 6-10) for the initial state and the final stationary MHD-equilibrium. The evolution of these quantities during the iteration is shown in Figure 3. We monitor also the weighted angle between magnetic field and plasma flow as defined in Equation 11. flow all terms would be zero. While this initial state contains a force-free magnetic field configuration and a spherically symmetric hydrodynamic solar wind, the magnetic field and plasma flow are not aligned. The averaged weighted angle between v and B is 31.1 • and the corresponding term of the functional L angle(B,v) is more than two orders of magnitude larger than the discretisation errors in final state after relaxation. We minimize the functional L iteratively with a steepest descent method, which has been successfully used in the global nonlinear force-free and magnetohydro-static optimization codes (see Wiegelmann, 2007;Wiegelmann et al., 2007, for details). Figure 3 shows the evolution of L in black and its individual terms: L force in green,L divB in blue, L cont in orange and L angle (B,v) in red. While the L force and L divB term stay approximately constant, the term L angle(B,v) decreases rapidly until it is of the order of the discretisation errors of the initial equilibrium, which corresponds to an average weighted angle between flow and magnetic field of 2 • . An exception is the term L cont which increases somewhat. To understand the evolution of the L cont term we must consider that while the initial force-free magnetic field is already complex, the plasma flow is smooth and strictly radial in the initial state, which results in a low initial value of L cont . This is not the case in the final equilibrium, when the flow becomes field aligned and thereby more complex. In the final state all terms of L are roughly of the same order and of the level of the discretisation error of the initial NLFFF-solution (see also Table 1, line 'Stationary MHD'). It can therefore be considered that a solution of the stationary MHD equations with field aligned plasma flow has been reached to a good approximation. Figure 4 shows in the top panel 'NLFFF' selected field lines for the initial force-free equilibrium. In the bottom panel 'FlowMHD' field lines (same starting points) of the final stationary MHD equilibrium are shown. It seems that low lying field lines are hardly affected, whereas further away from the sun the field lines become considerably more radial in the stationary MHD equilibrium. Physically this can be understood in the sense that the solar wind flow stretches and opens the field lines. In Figure 5 we investigate the differences more quantitatively. Panel a) shows how radial the magnetic field is as a function of the radius, where |B r | and |B| have been averaged over the whole sphere. The black line 'NLFFF' shows the initial force-free field and the red line 'FlowMHD' the stationary MHD solution. A value of unity means that the field is fully radial while a value of zero means the field is horizontal. The plot corroborates that the magnetic field becomes more radial on average with increasing distance from the Sun. At low heights both solutions agree, but above about 2R s the stationary MHD solution becomes significantly more radial then the initial field, as was already qualitatively seen in the field line plots in Figure 4. Figure 5 panel b) shows the relative quadratic difference (quadratic difference of the magnetic field vectors averaged over the entire sphere at each radial distance and divided by the averaged magnetic field strength). This plot confirms that below about 2R s the initial force-free and final stationary MHD equilibria are almost identical and deviate higher up in the corona. Figure 6 shows some quantities at a height of r = 2.5R s , which is usually the distance the source surface in located in PFSS-models. Panel a) shows the radial magnetic field component B r . Between the positive (red) and negative (green) magnetic field areas is the polarity inversion line. Around this line the plasma pressure increases (panel b) and the kinematic pressure p kin = ρv 2 /2 of the solar wind decreases (panel c). A reduced kinematic pressure around the polarity inversion line happens naturally, because here the magnetic field and consequently the aligned plasma flow are not dominated by the radial component. Table 2. Comparison of initial and final magnetic field. The first four columns contain information on the initial state and the investigated area. Column one names the investigated area, either the entire sphere 'Global', or an active region 'AR', or a quiet Sun area 'QS' (see the white and black boxes in Figure 1. Column two contains the initial magnetic field model, either a potential field 'Potential' or a nonlinear force-free field 'NLFFF'. The third, fourth and fifth column contain the temperature on the initial hydrodynamic equilibrium, the sound velocity cs and the related critical radius rc where the solar wind passes the sound velocity. The remaining columns six, seven and eight contain the results. Computed is the relative quadratic differences of the initial and final magnetic field as a function of the radius and subsequently it is checked at which radial distance from the sun the difference exceeds 10 −4 , 10 −3 and 10 −2 , respectively.

Area
Initial B Initial T cs rc 10 −4 10 −3 10 −2 An important question is to which extent simpler force-free extrapolations are justified in some areas. From the equilibrium investigated so far, it seems that non-force-free effects become important above about 2R s . In the following we want to quantify this further and investigate how the radii where the transition from force-free to non-force-free occur depend on the initial conditions. We also investigate if the results based on global averaging (marked 'Global' in Table 2) are different from a selected active region (marked 'AR' in Table 2 and with a white box in Figure 1 top panel) and a selected quiet sun area (marked 'QS' in Table 2 and with a black box in Figure 1). We investigate the effect of the initial magnetic field (either a nonlinear force-free field (marked with NLFFF) or a potential field (marked Potential). We also investigate the effect of the initial hydrodynamic Parker solar wind equilibrium by varying the temperature. Similar as in Figure 5 panel b) we compute the relative quadratic difference of the initial magnetic field and the final stationary MHD-equilibrium. In Table 2 we check at which radii the relative quadratic differences exceed thresholds of 10 −4 , 10 −3 , and 10 −2 , respectively.
The areas with r below the radius defined by a relative quadratic difference of 10 −4 can be considered as force-free. For areas with r larger than the radius defined by a relative quadratic difference of 10 −2 the effects of plasma flow and plasma forces are important. Between the two radii the transition between forcefreeness and stationary MHD take place. For the equilibrium investigated earlier (first row in Table 2) this transition area is between 1.80R s and 2.54R s for the global case. Restricting the computations on a selected active region or quiet sun area hardly changes anything. The corresponding radii change only by less than 0.1R s (upward for the active region and downward for the quiet sun). The influence of using a potential field as initial state instead of a force-free one is very small, too and leads to a very small (well below 0.1R s ) shift upwards. Changing the initial solar wind equilibrium has a significantly stronger effect. For a lower temperature the sound velocity c s decreases and the critical radius r c where the flow transits the sound velocity increases. The opposite happens for an increased temperature. This results in a higher sound velocity and the radius r c , where the wind passes the sound velocity is lower. Consequently one has significant faster flows at lower heights. It is found that the plasma flow starts having a significant effect already at a radial solar distance of approximately r c /2.

Conclusions and Outlook
In this paper we developed an optimization principle for the computation of stationary MHD equilibria, using synoptic vector magnetograms as boundary condition. As a first step we used a rather simple spherically symmetric and isothermal solar wind model as initial condition for the plasma and flow variables. We found that the newly developed code converged towards a stationary MHD equilibrium because the residual errors of the equilibrium are comparable with discretisation errors of nonlinear force-free solution. Below about 2R s the MHDsolution is very similar to a nonlinear force-free field. The reason is that because of the low plasma β and small Alfvén Mach number M A in these regions, plasma and flow effects can be neglected. Higher up in the corona the solar wind flow stretches out and opens up the magnetic field line. Such an effect of the solar wind is considered already in PFSS-model, but with an artificial source surface, whereas in stationary MHD the influence of the wind flow on the magnetic field is computed self-consistently. While the choice of the initial magnetic field configuration hardly influences this result, the profile of the outflow velocity is important. Fast flows in lower heights do obviously require that the effects of these flows have to be taken into account.
In the application of the method there is certainly room for improvement. Naturally one could relax the isothermal condition used here just for reasons of simplicity and use a more involved energy equation. In principle it is also possible to take measurements into account for the solar wind profile, which does not necessarily have to be spherically symmetric in the initial state. Including the observations of coronal loops to constrain the magnetic field as done in the forcefree models of e.g. Aschwanden (2013); Chifu, Wiegelmann, and Inhester (2017) can in principle also be incorporated in stationary MHD models. It is a useful feature of optimization principles that additional constraints are straightforward to implement by additional terms in the functional L. Figure 1. As boundary data we use a synoptic vector magnetogram from SDO/HMI for Carrington rotation 2099, which has been observed between 13/07/2010 and 09/08/2010. The polar regions have been cut out. One active region (white box in top panel) and a quiet sun area of the same size (black box) are investigated in Section 3.5. Figure 2. The variation of the initial conditions, based on a combination of a force-free magnetic field and a spherically symmetric Parker wind solution. a) Shows the different force terms b) Solar wind velocity as function of r c) Plasma and kinematic pressure from Parker's model, magnetic pressure from force-free field model. d) Plasma Beta and Alfvén Mach number: If both quantities are small, plasma and flow forces can be neglected for computing the coronal