Magnetohydrodynamic Equilibria

Thomas Wiegelmann

doi:10.1093/acrefore/9780190871994.013.3

Summary

Magnetohydrodynamic equilibria are time-independent solutions of the full magnetohydrodynamic (MHD) equations. An important class are static equilibria without plasma flow. They are described by the magnetohydrostatic equations

$\begin{array}{l} j \times B = \nabla p + ρ \nabla Ψ, \\ \nabla \times B = μ_{0} j, \\ \nabla \cdot B = 0. \end{array}$

$B$ is the magnetic field, $j$ the electric current density, $p$ the plasma pressure, $ρ$ the mass density, $Ψ$ the gravitational potential, and $µ_{0}$ the permeability of free space. Under equilibrium conditions, the Lorentz force $j \times B$ is compensated by the plasma pressure gradient force and the gravity force.

Despite the apparent simplicity of these equations, it is extremely difficult to find exact solutions due to their intrinsic nonlinearity. The problem is greatly simplified for effectively two-dimensional configurations with a translational or axial symmetry. The magnetohydrostatic (MHS) equations can then be transformed into a single nonlinear partial differential equation, the Grad–Shafranov equation. This approach is popular as a first approximation to model, for example, planetary magnetospheres, solar and stellar coronae, and astrophysical and fusion plasmas.

For systems without symmetry, one has to solve the full equations in three dimensions, which requires numerically expensive computer programs. Boundary conditions for these systems can often be deduced from measurements. In several astrophysical plasmas (e.g., the solar corona), the magnetic pressure is orders of magnitudes higher than the plasma pressure, which allows a neglect of the plasma pressure in lowest order. If gravity is also negligible, Equation 1 then implies a force-free equilibrium in which the Lorentz force vanishes.

Generalizations of MHS equilibria are stationary equilibria including a stationary plasma flow (e.g., stellar winds in astrophysics). It is also possible to compute MHD equilibria in rotating systems (e.g., rotating magnetospheres, rotating stellar coronae) by incorporating the centrifugal force. MHD equilibrium theory is useful for studying physical systems that slowly evolve in time. In this case, while one has an equilibrium at each time step, the configuration changes, often in response to temporal changes of the measured boundary conditions (e.g., the magnetic field of the Sun for modeling the corona) or of external sources (e.g., mass loading in planetary magnetospheres). Finally, MHD equilibria can be used as initial conditions for time-dependent MHD simulations. This article reviews the various analytical solutions and numerical techniques to compute MHD equilibria, as well as applications to the Sun, planetary magnetospheres, space, and laboratory plasmas.

1. Introduction

The evolution of many complex physical systems such as planetary magnetospheres and stellar coronae can be distinguished by their dynamic phases (magnetic storms, solar eruptions, flares) and slowly evolving quasi-static phases. During the slow phases, the plasma and magnetic field are in equilibrium and the time-dependent dynamic terms in full magnetohydrodynamics (MHD) can be neglected. Even during these quiet times, the systems can slowly evolve. One example for a quiet phase is the mass loading of the Earth’s magnetotail: The magnetosphere is formed by interaction of the solar wind with the Earth’s magnetic field. The solar wind compresses the Earth’s magnetic field on the dayside and stretches it to a magnetotail on the nightside. While the magnetic and plasma forces compensate each other in the magnetotail, energy and plasma transport from the solar wind into the magnetotail takes place, a process that can be well described as a sequence of static equilibria. An analysis of these equilibria can tell whether the system can become unstable and transit toward a dynamic phase (e.g., a magnetospheric substorm). A similar distinction of quasi-static and dynamic processes can be applied to coronae of the Sun and other stars. During quiet times, the photosphere and corona are magnetically connected and the coronal plasma is in equilibrium and can be modeled using magnetic field measurements in the photosphere (visible surface) of the Sun and stars as boundary conditions. Plasma convection below the surface causes slow changes in the boundary conditions and thereby coronal models based on these measurements evolve in time. They thus show, for example, the accumulation of magnetic energy and the formation of strong electric currents. Again, an analysis of these coronal equilibria can predict the possibility and likelihood of a transition to a dynamic phase with flares and coronal mass ejections. While equilibrium theory does not allow the investigation of the dynamic phase itself, it is very useful to study the slow evolution of systems and to investigate the physical conditions leading to instabilities, which are, in fact, associated with a fast dynamic evolution.

j \times B = \nabla p + ρ \nabla Ψ,

(1)

\nabla \times B = μ_{0} j,

(2)

\nabla \cdot B = 0.

(3)

Plasma equilibria are described by the magnetohydrostatic (MHS) equations (Equations 1–3), where Equation 1 describes an equilibrium between the Lorentz force, the plasma pressure gradient force, and the gravitational force. The electric current density is computed via Ampère’s law (Equation 2). Inserting this expression in the Lorentz force makes the equations nonlinear and, consequently, it is extremely difficult to find analytical solutions in three dimensions (3D). Even fully numerical approaches are difficult to implement because of the nonlinearity. More easily accessible are configurations that are either only two dimensional (2D) (requiring some symmetry assumption) or cases for which the nonlinearity can be transformed away. While it is certainly possible to combine these approaches, there has historically been a clear distinction between problems that can be computed analytically and those where powerful computers are needed. Before powerful computers and complex simulation codes became common tools, many works concentrated on basic physical understanding, the analytical investigation of nonlinear effects, and the development of self-consistent models. One example for physical systems where nonlinear and self-consistent MHD equilibrium models have been studied for about four decades is the Earth’s magnetosphere, which may be modeled by MHS models without gravity. In this case, the MHS equations in 2D can be reduced to a single partial differential equation, the Grad–Shafranov equation. This is a nonlinear partial differential equation. The Grad–Shafranov equation, which is here derived from the MHS equations, is one of the key equations for plasma physics because it is also highly relevant for a kinetic plasma description. A one-dimensional (1D) solution, the so-called Harris sheet, was originally derived from kinetic plasma physics but is a MHS equilibrium as well. An extension of the Grad–Shafranov theory allows additional terms in 2D MHS equilibria (gravity, magnetic shear, or guide field flow). Even investigating time-dependent effects as a sequence of equilibria (under certain constraints) is possible. In some cases, analytical solutions are still possible, but often at least part of the problem has to be treated numerically by, for example, solving an ordinary differential equation, which can be easily done by solvers available in numerical packages. The user of such packages hardly notices any difference from analytical solutions, and these methods are therefore often dubbed semi-analytical. Two-dimensional MHS equilibrium theory with a Grad–Shafranov of how this theory has been used to model the Earth’s magnetosphere is given in section “Earth’s Magnetosphere.”

Nonetheless, since the early era of computers, there has been a demand for solutions in full 3D to model, for example, the solar coronal magnetic field, which is difficult to measure. Historically, a number of solutions have been found using linearized equations such as linear force-free fields, which are still popular. In these models, the Lorentz force vanishes, which means that electric currents either vanish (potential fields) or they are co-aligned with the magnetic field. An advantage of linear equations is that solutions of the equations can be superposed to create new solutions. Well-known examples are solutions in the form of orthogonal function systems like spherical harmonics and Bessel functions. While these classes of solutions are in principle analytical (for each base function), computers are needed to calculate the large number of coefficients (e.g., for the spherical harmonics) in order to match these models to input data (e.g., magnetograms in the solar photosphere). The section “Magnetohydrostatics in 3D” outlines some basic theory of equilibria in 3D, including the theory of Euler potentials as well as methods of how MHS equilibria can be linearized. As an add-on, similarities of MHS and stationary hydrodynamics (SHD) are investigated. Applications of the 3D theory to the Sun are discussed in “Solar Magnetic Fields and Plasma” and to rotating planetary magnetospheres and magnetic stars in “Rotating Planetary Magnetospheres.” Applications to fusion plasmas and, in particular, tokamaks are discussed briefly in the section “Laboratory Plasmas.”

Nowadays it is also possible to compute nonlinear equilibria in full 3D, which require sophisticated numerical methods and considerable computational resources. Often simpler, linear models are used as an initial guess for the computation of complex nonlinear equilibria. These kinds of nonlinear 3D models are not discussed in the form of a general theory, but the specification (e.g., which quantities can be measured?) of particular physical systems needs to be taken into account, which is demonstrated in “Solar Magnetic Fields and Plasma.”

2. Theory

2.1 From Magnetohydrodynamics to Equilibrium Theory

The time-dependent ideal magnetohydrodynamics (MHD) equations are (see any textbook on plasma physics; e.g., Priest, 1982, 2014, 2019; Biskamp, 1993; Goedbloed & Poedts, 2004; Schindler, 2006) as follows:

\frac{\partial ρ}{\partial t} + \nabla \cdot (ρ v) = 0,

(4)

ρ \frac{\partial v}{\partial t} + ρ (v \cdot \nabla) v = j \times B - \nabla p - ρ \nabla Ψ,

(5)

\frac{\partial B}{\partial t} = \nabla \times E,

(6)

\nabla \times B = μ_{0} j,

(7)

\nabla \cdot B = 0,

(8)

E + v \times B = 0,

(9)

+ Equation of state, Energy equation,

(10)

where $B$ is the magnetic field, $E$ is the electric field, $j$ is the electric current density, $p$ is the plasma pressure, $ρ$ is the mass density, $v$ is the flow velocity, $Ψ$ is the gravitational potential, $t$ is the time, and $µ_{0}$ is the permeability of free space.

For a given gravitational potential $Ψ$ , seven evolution equations (Equations 4, 5, and 6) for the seven variables $B, v, ρ$ have to be solved. The electric current $j$ is linked to $B$ by Ampère’s law (Equation 7). The MHD equations are closed by an equation of state or an energy equation. A common simplification is the assumption of an isothermal plasma with constant $T$ or an isentropic plasma that directly links $ρ$ and $p$ . $E$ is linked to $B$ by ideal Ohm’s law (Equation 9). Finally, $B$ is subject to the solenoidal condition (Equation 8). Strictly, this is only a condition for the initial $B$ because once $\nabla \cdot B = 0$ , this holds forever. However, ensuring that this also holds numerically is much more demanding.

Usually the time-dependent MHD equations (Equations 4–10) can be solved only numerically. Under certain physical conditions, three kinds of simplification may become acceptable: neglecting physical terms, reducing the dimensions, and linearizing the equations.

•

For neglecting terms in the MHD equations, in particular, the following ensues:

–

No time-dependence (i.e., time derivatives are $\partial / \partial t = 0$ ), which leads to stationary MHD (investigated in “Stationary MHD With Plasma Flow”).

–

Additional neglect of plasma flow $v = 0$ leads to magnetohydrostatics (MHS) (see “Magnetohydrostatics in 2D” and “Magnetohydrostatics in 3D”). Often MHS is considered to be the standard approach for an MHD equilibrium.

–

$\partial / \partial t \approx O (ϵ) ≪ 1$ All quantities change so slowly in time that the dynamic term in MHD $ρ \partial v / \partial t$ can be neglected and the plasma is always near equilibrium. This leads to sequences of magnetostatic equilibria, which are studied in “Earth’s Magnetosphere.”

–

The gravity term $ρ \nabla Ψ$ is neglected. This assumption is well justified, for example, for modeling the Earth’s magnetosphere (see “Magnetohydrostatics in 2D”).

–

All forces except the Lorentz force $j \times B$ are neglected. This approach, called the (nonlinear) force-free field approach, is valid in the solar corona (up to a distance of about 2.5 solar radii). Force-free fields (including the special case of current-free potential or vacuum fields) are also relevant in some laboratory plasmas (e.g., they surround the high plasma $β$ area in tokamaks. Applications to the solar corona are studied in “Solar Magnetic Fields and Plasma” and laboratory plasmas in “Laboratory Plasmas.”

•

If symmetries of the system are exploited, then the number of independent variables are reduced from 3D to 2D, or even 1D. This approach is often used for studying magnetospheres and laboratory devices (axisymmetric).

•

Linearizing the original nonlinear equations can be done by special assumptions regarding the nature of electric currents and the dominance of a stationary background magnetic field, small v, or density gradients.

–

A popular assumption is that the electric current density is parallel (with a global constant $α$ ) to the magnetic field, which is the linear force-free approach (also called constant- $α$ force-free) (see the section “Linear Force-Free Fields”).

–

Another approach is that the currents are perpendicular to the gravity force, which allows a linearization of the MHS equations and is studied in “Linear MHS.”

–

The current density is perpendicular to any other external force (e.g., the centrifugal force). This allows under certain circumstances a linearization of the equations and a set of analytical solutions by separation of variables. This kind of equilibrium is studied in “Rotating Planetary Magnetospheres.”

Additionally, reducing the number of dimensions is possible but usually not necessary because linear equations can be solved effectively also in full 3D (e.g., with the help of a fast Fourier transform). These linear approaches are therefore popular due to their easier mathematical treatment, but less realistic, because nonlinear effects are often important to reveal the true nature of complex physical systems like solar and space plasmas.

2.2 Magnetohydrostatics in 2D

In the following, equilibria with one ignorable spatial coordinate are investigated. In this case, all physical quantities depend only on two spatial variables, but the vector fields can have non-vanishing components in all three coordinates. It has been shown by Edenstrasser (1980) that only three classes of symmetric equilibria with one ignorable coordinate exist: straight, helical, and circular symmetries. The MHS equations can be reduced in these cases to a single partial differential equation, the (generalized) Grad–Shafranov equation. The corresponding coordinate systems are orthogonal. Widely used are Cartesian, cylindrical, and spherical coordinates, a formulation in other geometries (e.g., toroidal is possible). In the following, Grad–Shafranov theory is described in Cartesian geometry. Spherical and cylindrical geometries are applied in the sections “Rotating Planetary Magnetospheres and Rotating Magnetic Stars” and “Laboratory Plasmas.” The most general case is a helical symmetry (see Biskamp, 1993, section 3.2). A 2D Cartesian coordinate system where $x$ and $z$ are the non-invariant coordinates and the configuration is translationally invariant in $y$ is used here. The magnetic field is given in terms of a flux function $A$ by

B = \nabla A \times e_{y},

(11)

where $e_{y}$ is the unit vector in the $y$ direction. This approach automatically fulfills the solenoidal condition (Equation 14). For simplicity, gravity and the magnetic field component in the invariant direction $B_{y}$ are neglected.

j \times B = \nabla p,

(12)

\nabla \times B = μ_{0} j,

(13)

\nabla \cdot B = 0.

(14)

In this approximation, surfaces $A = constant$ defined as so-called magnetic surfaces (e.g., L shells in magnetospheres). The plasma pressure is constant along magnetic field lines, which are contours of $A$ . If the plasma pressure becomes a function of the flux function, $p = p (A)$ , it is therefore automatically constant along field lines (a sufficient condition). To ensure that $p$ is constant along field lines, it is not necessary, however, that $p$ is a function of $A$ in the strict sense. On topologically distinct field lines with the same value of $A$ , the pressure can have different values on the distinct lines, but it is only required to be constant along each field line. By inserting Equation 11 in Equations 12 and 13, the magnetostatic equations reduce to a Grad–Shafranov equation.

2.2.1 Grad–Shafranov Equation

In the following, the Grad–Shafranov equation in Cartesian geometry is discussed. For a formulation in cylindrical geometry, see “Laboratory Plasmas.”

- Δ A = j_{y} (A) = \frac{\partial p}{\partial A} .

(15)

The functional form of $A \to p (A)$ or equivalent $A \to j_{y} (A)$ is still undefined. It is this dependence as well as the boundary conditions which shape the field lines and magnetic surfaces. For any choice, one obtains solutions of the MHS equations. A linear dependence in $j_{y} (A)$ corresponds to a quadratic form in $p (A) \propto A^{2}$ . More sophisticated than arbitrarily choosing such a functional form is to derive the relation $A \to p (A)$ from first principles by using kinetic plasma theory or statistical mechanics. Under the condition of a local thermodynamic equilibrium and a quasi-neutral plasma. one obtains (see Chapter 6.2 in Schindler, 2006):

p (A) \propto \exp (- A) .

(16)

2.2.2 Liouville Equation

With this approach, the Grad–Shafranov equation (Equation 15) takes the mathematical form of a Liouville equation (Liouville, 1853):

- Δ A = λ \exp (- A) .

(17)

This equation has applications not only in plasma physics, but also in other areas of physics (e.g., hydrodynamics). For examples see Walker (1915) for an application to a self-gravitating gas layer for nebular theory and Schmid-Burgk (1967) for astrophysical applications. As is seen in “Stationary MHD With Plasma Flow,” magnetostatic plasma equilibria and stationary flows of a fluid can be described with a similar mathematical structure. This nonlinear partial differential equation in 2D has the nice property that it may be solved analytically with the help of a generating function. In this approach, the 2D partial differential equation is reduced to an equation with one complex variable in the form

4 \frac{\partial^{2} A}{\partial u \partial \bar{u}} = λ \exp (- c A),

(18)

where $u = x + i y$ and $\bar{u} = x - i y$ for configurations invariant in $z$ and where $x, y, z$ are the Cartesian coordinates. Alternatively used is $u = x + i z$ to model a 2D plane $x, z$ for configurations invariant in $y$ . This equation has the general solution

A (u, \bar{u}) = \frac{2}{c} \ln \frac{1 + \frac{c λ}{8} | ξ {(u)}^{2} |}{| \frac{\partial ξ}{\partial u} |},

(19)

where any analytical function of $ξ (u)$ generates an exact solution of Equation 18 in the form given in Equation 19, and thereby a magnetostatic equilibrium. An analysis of the mathematical procedure, existence, proofs, and so on is outside the scope of this article (but can be found, e.g., in Bandle, 1975).

2.2.3 Harris Sheet

In 1D (with $z$ , say, as the independent variable), the Cartesian Grad–Shafranov equation reduces to

- \frac{d^{2} A (z)}{d z^{2}} = \exp (- 2 A),

(20)

which can be integrated directly analytically or by using the generating function $ξ (u) = \exp (i u) = \exp [i (x + i z)]$ in the Liouville approach in Equation 19. It has the Harris sheet (see Harris, 1962) as solution. The Harris sheet is a well-known equilibrium solution in 1D, where plasma and magnetic pressure are in equilibrium and contain a smooth current. This equilibrium is not only a solution of MHS, but remains valid in kinetic theory and has been originally found by Harris for kinetic equilibria. The Harris sheet is often a good local approximation for the plasma sheet in the magnetotail. As presented in “Earth’s Magnetosphere,” even involved sequences of magnetospheric MHS equilibria are based on modifications of the Harris sheet. For the flux function, one derives after integration:

Open in new tab

Figure 1. The Harris sheet. All quantities depend only on $z$ . Shown is the magnetic field $B_{x}$ (solid), the plasma pressure $p$ (dashed), the magnetic pressure $B_{x}^{2} / 2$ (dash-dotted), and the electric current density $j_{y}$ (dash-double-dotted). The dotted line shows the sum of magnetic and plasma pressure, which is constant in a Harris sheet.

A (z) = \ln \cosh (z),

(21)

from which all other quantities [magnetic field $B_{x} (z)$ , plasma pressure $p (z)$ , current density $j_{y} (z)$ ] can be computed:

B_{x} = - d A / d z = - \tanh (z),

(22)

p = \frac{1}{2} \exp (- 2 A) = \frac{1}{2 \cosh^{2} (z)},

(23)

j_{y} = - d^{2} A / d z^{2} = - \frac{1}{\cosh^{2} (z)} .

(24)

Please note that all other components of the magnetic field and current density vanish. Figure 1 shows a plot of the quantities as a function of $z$ as well as (dotted) the total pressure $p + B^{2} / 2$ , which is constant in the Harris sheet.

2.2.4 Including Magnetic Shear, 2.5D Configurations

A natural extension of the 2D theory is to take into account a magnetic field component $B_{y}$ in the invariant $y$ direction. The magnetic field is represented as

B = \nabla A \times e_{y} + B_{y} e_{y},

(25)

and the Grad–Shafranov equation becomes

- Δ A = \frac{\partial}{\partial A} (p (A) + \frac{B_{y}^{2} (A)}{2}) .

(26)

This approach is often called 2.5D because the equilibria depend on only two independent coordinates $(x, z)$ but contain all three components of the magnetic field and electric current. A particular case is the choice $p (A) = 0$ , which leads to force-free equilibria in 2.5D. Force-free equilibria are particularly popular to model the solar corona. The magnetic field component $B_{y} (A)$ in the invariant $y$ direction is a function of the flux function. Applications to the solar corona are discussed in “Solar Magnetic Fields and Plasma.”

2.2.5 Including Gravity

For equilibria with an external gravity field, the pressure can be written as

P = P (A, Ψ) = p (A) \cdot f (Ψ),

(27)

where $Ψ$ is the gravitational potential. For a constant gravitational acceleration, one finds:

P = P (A, Ψ) = p (A) \cdot \exp (- \frac{z}{H}),

(28)

where $H$ is the scale height. Configurations including both the plasma pressure and gravity are important for modeling regions of solar and stellar atmospheres, where the force-free assumption is not valid (e.g., in the chromosphere). Particular 2D solutions with gravity are not discussed here, but configurations with gravity are considered in more detail in the sections “Linear MHS” and “Nonlinear MHS,” where linear and nonlinear magnetostatic equilibria in 3D are addressed.

2.3 Magnetohydrostatics in 3D

To find solutions of MHS equations in full 3D without making any symmetry assumptions is extremely difficult. Some basic theory on how this can be done is outlined here. Numerically computed equilibria with application to the Sun are presented in “Solar Magnetic Fields and Plasma.”

2.3.1 Euler Potentials

A generalization of the flux function to 3D are Euler potentials

B = \nabla α \times \nabla β,

(29)

which reduce to the flux function representation defined in Equation 11 when $α = A$ and $β = y + const$ . Euler potentials are not necessarily unique and might not exist globally in a given physical domain (see Chapter 5.1.2 in Schindler, 2006). Euler potentials can also be used to derive a vector potential for $B = \nabla \times A$ in the form $A = α \nabla β$ . It is obvious that by representing the magnetic field with Euler potentials, the solenoidal condition (Equation 14) is automatically fulfilled. Euler potentials have the useful property that magnetic field lines are defined as intersections of surfaces on which the Euler potentials $α$ and $β$ are constant. For MHS equilibria without gravity, the plasma pressure is constant on magnetic field lines and $p = p (α, β)$ . For configurations with gravity, one has $p = p (α, β, Ψ)$ , where $Ψ$ is the gravitational potential. By inserting the ansatz (Equation 29) into the MHS Equations 1–3, a system of nonlinear coupled partial differential equations can be derived (see, e.g., Zwingmann, 1987):

\nabla β \cdot \nabla \times (\nabla α \times \nabla β) = μ_{0} \frac{\partial p (α, β, Ψ)}{\partial α},

(30)

- \nabla α \cdot \nabla \times (\nabla α \times \nabla β) = μ_{0} \frac{\partial p (α, β, Ψ)}{\partial β},

(31)

- ρ = \frac{\partial p (α, β, Ψ)}{\partial Ψ} .

(32)

At first sight, hardly anything is gained by reformulating the MHS equations in this form, because Equations 30–32 look even more complicated than the original MHS equations and they contain higher-order derivatives. Furthermore, the use of the definition (Equation 29) leads intrinsically to nonlinear equations and a linearization, as is explained in “Linear Force-Free Fields,” and “Linear MHS” is unsuitable here. As a consequence, the Euler potential approach for computing 3D equilibria never became as popular as the flux function representation and the associated Grad–Shafranov equation in 2D. Nevertheless, for specific cases, 3D solutions of Equations 30–32 have been found. The non-uniqueness of Euler potentials does have an advantage, however. If an MHS solution in the form of Euler potentials is found, this non-uniqueness can be used to construct solutions of the stationary MHD equations, with field-aligned incompressible plasma flow, with the help of a transformation procedure introduced by Gebhardt and Kiessling (1992). This procedure is briefly discussed in the section “Stationary MHD With Plasma Flow.”

2.3.2 Potentials for the Current Density

It is interesting to note (see Low, 1991) that it is possible to represent also the electric current density with potentials in the form

μ_{0} j = \nabla U \times \nabla F,

(33)

where $U (x, y, z)$ and $F (x, y, z)$ are functions of the coordinates. A representation in spherical and cylindrical geometry is also popular. This representation is in particular helpful when studying a system containing external forces, which can also be computed as the gradient of a potential. Examples are an external gravity field (see the section “Linear MHS”) and rotating systems with a centrifugal force (see “Rotating Planetary Magnetospheres”), or both (in the section “Rotating Magnetic Stars”). Choosing one of the potentials in Equation 33 proportional to the potential of the external force allows, under certain circumstances, a linearization of the force–balance equation with analytical solutions by using separation of variables.

2.4 Stationary MHD With Plasma Flow (SMHD)

The equations for stationary MHD equilibria with incompressible plasma flow are given as:

ρ (v \cdot \nabla) v = j \times B - \nabla p,

(34)

\nabla \times B = μ_{0} j,

(35)

\nabla \cdot B = 0,

(36)

\nabla \cdot ρ v = 0,

(37)

\nabla \times (v \times B) = 0,

(38)

\nabla \cdot v = 0,

(39)

where $v$ is the plasma velocity. Note that gravity has been neglected because it would require also the inclusion of compressibility. For field-aligned flows, Equation 38 is fulfilled automatically. After applying the vector identity $(v \cdot \nabla) v = \nabla v^{2} / 2 - v \times (\nabla \times v)$ to the left-hand side of Equation 34, one can see that the set of Equations 34–39 has two clearly distinguishable subsets: MHS and stationary incompressible hydrodynamics (SHD).

2.4.1.1 MHS as SMHD Subclass

If all terms containing the plasma flow $v$ are neglected, Equations 34–39 reduce to MHS:

\frac{1}{μ_{0}} (\nabla \times B) \times B - \nabla p = 0,

(40)

\nabla \cdot B = 0.

(41)

2.4.1.2 SHD as SMHD Subclass

However, if all terms containing the magnetic field $B$ are neglected, Equations 34–39 reduce to the equations of SHD:

- ρ v \times (\nabla \times v) - \nabla (\frac{1}{2} ρ v^{2} + p) = 0,

(42)

\nabla \cdot v = 0.

(43)

2.4.1.3 Mathematical Similarities Between SHD and MHS

It is interesting to note that MHS (Equations 40 and 41) and SHD (Equations 42 and 43) have a similar mathematical structure of the form

(\nabla \times c) \times c = \nabla Π,

(44)

\nabla \cdot c = 0,

(45)

where for MHS

c = B,

(46)

Π = μ_{0} p,

(47)

and for SHD

c = \sqrt{ρ} v,

(48)

ρ_{2}

Π = - \frac{ρ}{2} v^{2} - p .

(49)

This similarity has been pointed out in several publications (e.g., Grad, 1960; Tsinganos, 1981; Tasso & Throumoulopoulos, 1998). Gebhardt and Kiessling (1992) provided a transformation technique to compute stationary MHD equilibria if MHS equilibria in the form of Euler potentials (see the section “Euler Potentials”) have been found. Similarly, the velocity field can also be represented as

v = \nabla f \times \nabla g,

(50)

which automatically fulfills the incompressibility condition (Equation 39). As the Euler potentials are not unique, a suitable transformation between the Euler potentials of the magnetic field and velocity field may be found. For the theory of these transformations in 3D, see the original publication by Gebhardt and Kiessling (1992). The following section contains the main features of these transformation techniques for systems in 2D, where the Euler functions reduce to a flux function, as shown for MHS in “Magnetohydrostatics in 2D.”

2.4.1 SMHD Equilibria in 2D

Here one aims to solve the stationary MHD equations (34–39) in 2D and assumes that one has already found a solution of the corresponding MHS equations (40 and 41). Equation 38 is solved by assuming that the velocity is parallel (or antiparallel) to the magnetic field. It is convenient to introduce the Alfvén speed

v_{A} = \frac{B}{\sqrt{μ_{0} ρ}},

(51)

and the Alfvén Mach number

M_{A} = \frac{v}{v_{A}},

(52)

which reduces the stationary MHD equations (34–39) to

B \cdot \nabla M_{A} = 0

(53)

\nabla Π = \frac{(1 - M_{A}^{2}) (\nabla \times B) \times B}{μ_{0}} - \frac{| B |^{2}}{2 μ_{0}} \nabla (1 - M_{A}^{2}),

(54)

\nabla \cdot B = 0.

(55)

For static equilibria without plasma flow $(M_{A} = 0)$ , these equations reduce to the well-known MHS equations (40 and 41). Just as in “Magnetohydrostatics in 2D,” the magnetic field in 2D is represented by a flux function,

B = \nabla α \times e_{y} .

(56)

Since flux functions (like Euler potentials) are not unique, it is possible to transform to another flux function $A$ . From Equation 53, $M_{A}$ is constant on magnetic field lines and therefore is a function of the flux functions $[M_{A} = M_{A} (α) = M_{A} (A)]$ . This property is used to simplify Equation 54 by the choice

(1 - M_{A}^{2}) {(\frac{\partial α}{\partial A})}^{2} \equiv 1,

(57)

which eliminates all terms containing the Mach number in Equation 54 and the stationary MHD equations (53–55) become

- \frac{1}{μ_{0}} Δ A = \frac{\partial Π}{\partial A},

(58)

α = \pm \int \frac{d A}{\sqrt{1 - M_{A} {(A)}^{2}}} .

(59)

Equation 58 obviously has the same mathematical form as the Grad–Shafranov Equation 15 derived for MHS equilibria. Equation 59 has been derived by integrating Equation 57. With help of an appropriate Mach number profile $M (A)$ , Equation 59 can be used to deduce a stationary MHD equation from the solution of the Grad–Shafranov equation. The transformation procedure can become quite complicated, however, and special care has to be taken if trans-Alfvénic flows are considered because for $v = v_{A}, M_{A} = 1$ and the denominator in Equation 59 vanishes. An explicit example of such a transformation (limited to sub-Alfvénic flows) is presented in the section “Magnetotail-like Solutions of the Liouville Equation.”

3. Applications

3.1 Earth’s Magnetosphere

The Earth’s magnetosphere is formed by the interaction of the solar wind with the Earth’s magnetic field. In a vacuum, the field would be roughly dipolar, but due to the interaction with the solar wind, the Earth’s magnetic field is compressed at the dayside and stretched out to a long tail (more than about 100 R_E ) at the nightside, known as the “magnetotail.” Locally the current sheet in the tail is often modeled using the Harris sheet model (see section “Harris Sheet”), but it is possible to construct more sophisticated equilibrium models analytically, as is described.

3.2 Magnetotail-Like Solutions of the Liouville Equation

Several authors have applied the Liouville equation (Equation 17) to model magnetospheric plasmas, but the Liouville equation can also be derived from a Maxwell–Vlasov system for kinetic equilibria, as done, for example, in Kan (1973) to model a current sheet in the magnetotail. The true 2D character of the Liouville equations allows a derivation of an analytical magnetospheric configuration, including the transition region from dipolar to stretched configurations, as done in Manankova (2003). The tail approximation (see section “Magnetotail Convection Model”) cannot be applied here because the approximation of stretched configurations is not valid in the dipolar region. Schindler and Birn (2004) found magnetospheric solutions of the Liouville equation (Equation 17) using the generating function

ξ (u) = 2 \hat{l} \exp (i (u / \hat{l} + \sqrt{\frac{u / \hat{l}}{ϵ}})),

(60)

Open in new tab

Figure 2. (a) Field lines [contour plot of flux function $A (x, y)$ ] for the magnetotail-like solution equation (Equation 61) of the Liouville equation as found by Schindler and Birn (2004). (b) Corresponding electric current density, $- j$ , is plotted for better visualization. (c) and (d) Stationary MHD equilibria, created with the help of a transformation procedure.

Adapted from Nickeler and Wiegelmann (2010), Figure 6. The original was published under Creative Commons License 3.0.

where $\hat{l}$ is used for normalization and $ϵ$ is a free parameter. Inserting Equation 60 in the general solution (Equation 19) leads to the solution:

A (x, y) / \hat{A} = \ln (\frac{\cosh (\frac{y}{\sqrt{2 ϵ} \sqrt{r + x}} + y)}{\sqrt{\frac{1}{r} (\frac{1}{4 ϵ} + \sqrt{\frac{r + x}{2 ϵ}}) + 1}}),

(61)

with $r = \sqrt{x^{2} + y^{2}}$ , where the coordinates here are normalized with respect to $\hat{l}$ . For $ϵ \to \infty$ , this Ansatz for the generating function $ξ$ produces the Harris sheet solution (see “Harris Sheet”). Figure 2(a) shows the magnetic field lines for this solution and panel (b) the current density.

3.2.1 Including Field-Aligned Plasma Flow

As explained in “Stationary MHD With Plasma Flow,” stationary equilibria with field-aligned incompressible plasma flow can be constructed with the help of a transformation. This was used in Nickeler and Wiegelmann (2010) by adding a plasma flow in the lobe of the magnetotail solution discussed in “Magnetotail-like Solutions of the Liouville Equation.” Figure 2(c) shows the Alfvén Mach number profile $M_{A} (x, y)$ (ratio of plasma flow velocity to Alfvén Mach velocity), which was chosen in the form

M_{A} (A) = M_{1} + (M_{2} - M_{1}) \tanh (d (A - A_{c})),

(62)

where $M_{1}, M_{2}, d$ and $A_{c}$ are free parameters (see Nickeler & Wiegelmann, 2010). The variable $d \propto 1 / l_{flow}$ is an inverse length and controls the spatial scales of the plasma flow (Mach number on field lines). Note that $M_{A}$ in Equation 62 was limited to sub-Alfvénic flows and the transformation (Equation 59) was used to compute a stationary MHD solution from Equation 61. As a consequence, two additional current sheets form where the gradient of $M_{A}$ is large [see Figure 2(d)].

3.2.2 Sequences of Magnetostatic Equilibria

Physical processes in the magnetosphere can be separated into the slow evolution loading phase (energy in the magnetotail increases) and a dynamic eruptive phase (magnetospheric substorms), where energy and matter (plasma) are released. The slowly evolving loading phase can be well modeled with a sequence of magnetostatic equilibria (e.g., in the form of the model developed by Schindler & Birn, 1982), which is summarized in section ” Magnetotail Convection Model.”

3.2.3 Magnetotail Convection Model

The magnetotail convection model developed by Schindler and Birn (1982) starts from the time-dependent ideal MHD equations (4–10) without gravity, where the system is closed by an isentropic equation of state:

\frac{\partial}{\partial t} (\frac{p}{ρ^{γ}}) + v \cdot \nabla (\frac{p}{ρ^{γ}}) = 0

(63)

where $γ$ is the polytropic index. The original model by Schindler and Birn (1982) used a 2D geometry $(x, z)$ . as in Figure 3; the theory was extended later to 3D by Birn and Schindler (1983). The physical conditions in the Earth’s magnetotail allowed them to introduce a smallness parameter $ϵ$ to distinguish different scales (e.g., the length scale of variation in the $x$ coordinate changes much more slowly than in $z$ , leading to $ϵ = L_{z} / L_{x} ≪ 1$ ). There is a general theory of asymptotic expansions to find approximate solutions for systems with multiple scales (Schindler, 2006). Here specific application to the magnetotail is studied. The ratio of the speed of convection flows $v_{c}$ to the speed of MHD waves $v_{M H D}$ is of the order of $ϵ = v_{c} / v_{M H D} ≪ 1$ and the typical convection time-scale $t_{c}$ is long compared with the travel-time of MHD waves, leading to $ϵ = t_{M H D} / t_{c} ≪ 1$ . Rewriting Equations 4–10 with variables O(1) times powers of $ϵ$ yields a set of equations for every power of $ϵ$ The authors solved up to $ϵ^{1}$ , neglecting higher-order terms. Consequently, the flow and time-dependent terms vanish in the equation of motion (Equation 5), which reduces to the MHS approach where the Lorentz force is compensated by the plasma pressure gradient. In 2D, this equation reduces to the Grad–Shafranov equation (Equation 15). In other equations such as the ideal Ohm’s law (Equation 9), terms of the order of $ϵ$ remain. Using the definition of the magnetic flux function (Equation 11), ideal Ohm’s law is written as

Open in new tab

Figure 3. A 2D model of the Earth’s magnetosphere. The original model, developed by Schindler and Birn (1982), was restricted to the magnetotail marked here with a rectangular domain. In a subsequent work by Becker et al. (2001), the near-Earth dipolar magnetic field was included into the convection model. The shaded area marked with $V (A, t)$ is the volume between two field lines. In the mathematical limit of a vanishing distance between the two field lines, $V (A, t)$ becomes the differential volume.

Adapted from Schindler and Birn (1982), Figure 1. Reproduced with permission from John Wiley & Sons.

Open in new tab

Figure 4. A numerical solution of the magnetotail convection model (limited to the rectangular area in Figure 3). (a) A few example field lines [contour plot of flux function $A (x, z)$ ) for the initial state]. (b) The initial smooth current-density profile $j_{y} (x, z)$ . (c) and (d) The same quantities at the end of the loading phase. The tail-like magnetic field lines have been compressed and, as a consequence, a thin current sheet has formed.

Adapted from Wiegelmann and Schindler (1995), Figure 1 and 2. Reproduced with permission from John Wiley & Sons.

\frac{\partial A}{\partial t} + v \cdot \nabla A = \frac{d A}{d t} = 0.

(64)

Just as in MHS theory, the flux function is constant on magnetic field lines. Unlike a static equilibrium, flux function and pressure are time-dependent, with $A = A (x, z, t)$ and $p = p (A, t)$ . While at every time the configuration is in equilibrium (with the Lorentz force compensated by pressure gradient), both quantities change slowly and the temporal evolution has to follow Ohm’s law as well as the entropy equation (Equation 63). In the entropy equation, the mass density (constant on field lines) can be related to the differential volume $V$ (see Figure 3) of each field line in the form $V (A, t) = m_{i} N / ρ$ , where $N$ is the particle number density and $m_{i}$ the average particle mass. Taking the convective derivative $d / d t = \partial / \partial t + v \cdot \nabla$ into account, Equation 63 expresses entropy conservation for each flux tube:

\frac{d}{d t} (p {(A, t)}^{1 / γ} V (A, t)) = 0.

(65)

Using further $\partial A / \partial x = ϵ ≪ \partial / \partial z$ and neglecting terms of order $ϵ^{2}$ , the Grad–Shafranov equation becomes, after integration in $z$ ,

\frac{1}{2 μ_{0}} {(\frac{\partial A}{\partial z})}^{2} + p (A, t) = \frac{1}{2 μ_{0}} B_{0}^{2} (x, t) .

(66)

The integration constant $B_{0} (x, t)$ corresponds to the magnetic field in the lobe. This is a quantity that can be deduced from measurements and, consequently, $B_{0} (x, t)$ specifies the boundary conditions for the quasi-static model. Under the quasi-static equilibrium conditions defined earlier, the magnetic pressure in the lobe corresponds to the plasma pressure in the center of the magnetotail $P_{0} (x, t) = B_{0} {(x, t)}^{2} / 2$ and Equation 66 describes a local Harris sheet-like equilibrium (see “Harris Sheet”) in $z$ . Note that terms of the order of $E$ are considered and lead to a parametric dependence of Equation 66 on $x$ and $t$ via the $z$ integration constant $B_{0}$ . This constant is used to model the flux transfer from the front side of the magnetosphere to the magnetotail and is given by

B_{0} (x, t) = \sqrt{2} {(1 + t)}^{k_{1} / 2} {(1 + Λ {(1 + t)}^{k_{2}} x)}^{- n / 2},

(67)

where $k_{1}$ and $k_{2}$ are free parameters. They are deduced from measurements, as specified later. The measurements show $Λ = 0.036$ and $n = 1.1$ . The equation of state (Equation 63) reduces in this approximation to entropy conservation on each flux tube (Equation 65).

Schindler and Birn (1982) provided analytical solutions of the system (Equations 65–67) under the assumption that the plasma pressure evolves self-similarly in time when $p (A, t) = p_{1} (t) p_{2} (A)$ and $j_{y} (A, t) = j_{1} (t) j_{2} (A)$ . These analytical solutions are only possible for the choice $k_{2} = k_{1} (1 / γ - 1 / 2)$ in Equation 67. For the values deduced from measurements, namely $k_{1} = k_{2} = 2$ , the system of Equations 65–67 has to be solved numerically (Wiegelmann & Schindler, 1995). The numerical solution is not self-similar but shows a strong increase of the electric current density in the near-Earth magnetotail, as shown in Figure 4. While the formation of current sheets can be modeled by the sequence of MHS equilibria, plasma instabilities on microscopic scales (well below the MHD scale) can occur and lead to dynamic processes (e.g., magnetic substorms), which cannot be modeled by equilibrium theory. In this sense, a sequence of equilibria can explain why the slow loading phase of the magnetotail ends and the eruptive phase is initiated. While the Schindler and Birn (1982) model is limited to 2D configurations in the magnetotail, the theory has been extended to 3D by Birn and Schindler (1983). The approximation $L_{z} / L_{x} ≪ 1$ is valid in the magnetotail but breaks down in the dipolar region close to Earth. Including this area requires a full numerical approach (Becker et al., 2001). It was found that current sheets form in the transition region between the dipolar and tail-like field at the night-side magnetosphere.

3.3 Solar Magnetic Fields and Plasma

In the solar corona, the magnetic pressure is about two to four orders of magnitude higher than the plasma pressure (the dimensionless number $β = p / p_{mag} ≪ 1$ ), non-magnetic forces can be neglected in lowest order (Neukirch, 2005), and the Lorentz force approximately vanishes. Under this assumption, the MHS equations (1–3) reduce to

(\nabla \times B) \times B = 0,

(68)

\nabla \cdot B = 0.

(69)

This approach gives a force-free field model.

3.3.1 Force-Free Magnetic Fields

For a review article dedicated entirely to solar force-free fields, see Wiegelmann and Sakurai (2012). Within this article, only a few basic properties of this approach are addressed. The cross product in Equation 68 vanishes only if the two terms, $µ_{0} j = \nabla \times B$ and $B$ , are parallel (or antiparallel) to each other (or if one term vanishes). This leads to

\nabla \times B = α B,

(70)

B \cdot \nabla α = 0,

(71)

where $α$ is a function of space. Taking the divergence of Equation 70 and using the solenoidal condition (Equation 69) gives Equation 71, which constrains $α$ to be constant along magnetic field lines. One has to distinguish three cases:

•

$α = 0$ : The electric current density vanishes. This assumption gives rise to a potential field model and is described in the section “Potential Fields.”

•

$α$ is globally constant. Electric currents are parallel or antiparallel to the magnetic field. This approach is called a Linear Force-Free Field model (LFFF) and is discussed in section “Linear Force-Free Fields.”

•

$α$ varies in space, but is constrained by Equation 71. This is the generic case and called a Nonlinear Force-Free Field model (NLFFF). This approach is discussed in “Nonlinear Force-Free Fields for 2D” and in “Nonlinear Force-Free Fields in 3D.”

Open in new tab

Figure 5. A potential field source-surface model based on a synoptic magnetogram observed with SDO/HMI in March 2015.

3.3.2 Potential Fields

In current-free configurations, the magnetic field can be derived from a scalar potential

B = - \nabla ϕ,

(72)

And, after inserting in the solenoidal Equation 69,

Δ ϕ = 0

(73)

must be solved together with appropriate boundary conditions. For this model, the inner boundary conditions are the measured magnetic field in the solar photosphere (radial or line-of-sight component of the photospheric magnetic field). These measured photospheric magnetic fields provide the inner boundary conditions (at a radius $r = r_{sun}$ , which is the radius of the Sun) for the model. At the outer boundary of the model, which is called the source surface (often placed at $r = 2.5 r_{sun}$ ), the magnetic field is assumed to become radial. This is a simple way to mimic the effect of the solar wind flow in stretching and opening magnetic field lines. This approach is called the Potential Field Source Surface model; an example is shown in Figure 5. Equation 73 has been solved by a spherical harmonic decomposition, and a synoptic magnetogram observed with the Helioseismic and Magnetic Imager (HMI) on board the Solar Dynamics Observatory (SDO) was used as a photospheric boundary condition. It is also possible to solve Equation 73 in Cartesian geometry with different methods. This is not addressed here [see Priest (2014), section 3.3 for explicit potential field solutions in Cartesian and spherical geometry], but Cartesian potential fields as a special case with $α = 0$ are investigated in section “Linear Force-Free Fields.” Potential fields are easy to compute and require as a boundary condition only the line-of-sight (LOS) photospheric field. Thus, they have been used frequently in the past and are computed also nowadays as a lowest order approximation. One reason is that for a given vertical flux distribution on the boundaries, potential fields have a minimum energy [see Priest (2014), section 3.3 for a mathematical proof]. More sophisticated models (e.g., nonlinear force-free fields) have a higher energy. The amount of excess energy in these models is an upper limit for the free magnetic energy of the configuration. Only configurations with free magnetic energy can become unstable (a necessary but not sufficient condition) and power solar eruptions like flares and coronal mass ejections. During these dynamic processes, the free magnetic energy or parts of it are converted to thermal and kinetic energy. Another reason for computing potential fields is that they can be used as initial states for numerical computation of nonlinear models (see “Nonlinear Force-Free Fields in 3D”).

3.3.3 Linear Force-Free Fields

For linear force-free fields, the force-free function $α (r)$ in Equations 70 and 71 reduces to a constant. This automatically satisfies Equation 71, leading to

\nabla \times B = α B,

(74)

where now $α$ is a parameter. Like potential fields, linear force-free fields require only the radial or vertical magnetic field in the photosphere as a photospheric boundary condition. While it is possible to formulate the linear force-free equations in any geometry, this approach never became popular for spherical geometry because a constant value of $α$ in the entire global corona is unrealistic. Linear force-free fields are (and have been more in the past) popular for modeling solar active regions for two reasons: First, the equations are much easier (and faster) to solve, and second, they need only the vertical photospheric field as a boundary condition. Such measurements have been available for decades, whereas horizontal fields are more difficult to measure.

For a constant $α$ , a Helmholtz equation can be derived by taking the curl of Equation 70 and applying vector identities, leading to

Δ B + α^{2} B = 0,

(75)

which can be solved by a number of numerical techniques. The fast Fourier transform method (Alissandrakis, 1981) is widely used in the solar community [explicit analytical linear force-free solutions in Cartesian and spherical geometry can be found in Priest (2014), section 3.4]. While it is known that $α$ has to be constant, its exact value is a priori unknown. If vector magnetograms are available, the value of $α$ can be deduced from the ratio of the vertical components of $\nabla \times B$ and $B$ on the surface (Hagino & Sakurai, 2004):

α = \sum μ_{0} j_{z} sign (B_{z}) / \sum | B_{z} |,

(76)

where the vertical electric current density $j_{z} (x, y)$ is

j_{z} = \frac{1}{μ_{0}} (\frac{\partial B_{y}}{\partial x} - \frac{\partial B_{x}}{\partial y}) .

(77)

By applying these formulas, the photospheric magnetic field measurements are used as the boundary condition for the model and also to specify the free model parameter $α$ . Before the launch of SDO in 2010, vector magnetograms had been only occasionally available, whereas LOS magnetograms had been available for decades. Consequently, alternative methods to derive the constant $α$ have been developed, mainly by computing a number of LFFF models with the same boundary conditions but with several different values of $α$ . The magnetic field lines traced from the solutions of the LFFF models have been compared with loops seen in coronal EUV images. A number of quantitative criteria have been developed to optimize the value of $α$ which best fits the coronal observations (see, e.g., Carcedo et al., 2003; Conlon & Gallagher, 2010).

A single value of $α$ often does not fit well to structures (coronal loops) seen in different parts of an active region. The optimal value of $α$ can have a rather large scatter and even change its sign. Variations of $α$ , however, contradict the assumption of a linear force-free field, and this discrepancy can only be resolved if nonlinear effects are taken into account, as is done in the section “Nonlinear Force-Free Fields in 2D” and in the section “Nonlinear Force-Free Fields in 3D”. Furthermore, there is another problem with linear force-free fields, namely, that the amount of magnetic energy of these configurations is not bounded if the computations extend in a vertical direction from the photosphere to infinity. In such cases, the computed magnetic energy becomes unrealistic at large distances from the Sun. This property is not a problem, however, for computations of a solar-active region with a vertical extend comparable with the lateral extension of the active region. There are some limitations on the allowed values of $α$ (which depend on the size of the region) and non-uniqueness of linear force-free fields. It is well outside the scope of this article to discuss these details (but see Priest, 2014, section 3.4].

3.3.4 Nonlinear Force-Free Fields in 2D

Due to the difficulty of taking nonlinearities into account, they are first investigated in 2D by using again the method of prescribing the magnetic field with the help of a flux function, similar to Cartesian MHS equilibria in the section “Grad–Shafranov Equation.” Using the flux-function approach and deriving a Grad–Shafranov equation is possible also in spherical and cylindrical geometry. For curvilinear coordinate systems, the mathematics can become quite involved. The interested reader can find a derivation of the Grad–Shafranov equation for these cases in Schindler (2006, section 5.2). A class of rotationally symmetric nonlinear force-free equilibria has been developed in Low and Lou (1990). In spherical geometry, the magnetic field is represented as

B = \frac{l}{r \sin θ} (\frac{1}{r} \frac{\partial A}{\partial θ} e_{r} - \frac{\partial A}{\partial r} e_{θ} + Q e_{φ}),

(78)

Open in new tab

Figure 6. Several 2D-axisymmetric configurations. The color-coded magnetogram shows the vertical field in the photosphere. The white curves are closed magnetic field lines (also called magnetic loops), which have been traced starting from the negative magnetic polarity area (shown in blue).

Open in new tab

Figure 7. The axisymmetric NLFFF solution found by Low and Lou (1990) is very flexible in the sense that by rotation and shifting, 3D-like solutions in 3D-Cartesian geometry can be created. The NLFFF solution shown in Figure 6(e) was rotated by an angle of $φ = 0, π / 8, π / 4, π / 2$ in Figure 7(a–d), respectively. Figure 7(e) and (f) show 3D views of a potential field and MHS equilibrium, respectively. The black lines are closed magnetic field lines.

where $A$ is the flux function, and $Q$ is related to (but not identical with) the $φ$ component of the magnetic field $B$ . Note that in spherical geometry, the flux function is not a component of the vector potential. With this approach, one obtains a Grad–Shafranov equation with the independent spherical variables $r$ and $θ$ , whereas the configuration is axisymmetric in $ϕ$ by construction:

\frac{\partial^{2} A}{\partial r^{2}} + \frac{1 - μ^{2}}{r^{2}} \frac{\partial^{2} A}{\partial μ^{2}} + Q \frac{d Q}{d A} = 0,

(79)

where $µ = c o s θ$ . Low and Lou (1990) found solutions for

Q (A) = λ A^{1 + 1 / n} (α = \frac{d Q}{d A} ~ A^{1 / n}),

(80)

with a separation ansatz

A (r, θ) = \frac{P (μ)}{r^{n}} .

(81)

A few example field lines for this kind of nonlinear force-free field (NLFFF) solution are shown in Figure 6(e). This kind of solution is very popular to mimic 3D NLFFF equilibria in Cartesian geometry, as shown in Figure 7. In panel (a), the symmetry in $ϕ$ is still obvious. If the whole configuration is rotated by an angle of $φ = π / 8, π / 4, π / 2$ , as shown in Figure 7(b–d), this symmetry becomes less obvious after a transformation to Cartesian geometry and after cutting off regions that are located below the photosphere after rotation. These equilibria are often used to test the performance of numerical 3D nonlinear force-free extrapolation codes.

3.3.5 Including Plasma Pressure: 2D MHS With Shear

The ansatz (Equation 78) used by Low and Lou (1990) for force-free fields has been extended to magnetostatic equilibria (without gravity) by Wiegelmann and Neukirch (2006):

- \frac{\partial^{2} A}{\partial r^{2}} - \frac{1 - μ^{2}}{r^{n}} \frac{\partial^{2} A}{\partial μ^{2}} = μ_{0} r^{2} (1 - μ^{2}) \frac{d Λ}{d A} + b_{ϕ} \frac{d b_{ϕ}}{d A} .

(82)

For the choice $d Λ / d A = 0$ , this reduces to the force-free ansatz found by Low and Lou (1990); $Λ$ contains the non-magnetic forces due to a plasma pressure gradient. Self-similar (and thereby separable) solutions are only possible for the choice $Λ (A) = k A^{s}$ and the additional constraint that the power of the radial component r has to be the same in all terms (see Wiegelmann & Neukirch, 2006). Figure 6(f) shows an example solution.

3.3.6 Comparison of Different Models

All models have exactly the same vertical photospheric magnetic flux distribution (color-code in Figure 6) and the magnetic field lines have been traced starting from the negative magnetic polarity region. Note that the photospheric horizontal magnetic field components are different for different models. The different model assumptions are related to the electric current density:

•

No electric currents, which naturally is the simplest approach and is shown in Figure 6(a) and Figure 7(e). The loops connect the negative and positive polarity areas on the photosphere and perpendicularly cross the polarity inversion line.

•

The current is strictly parallel to the magnetic field lines. This is the force-free approach because the Lorentz force $j \times B$ vanishes. First shown are the effects of currents which are proportional with a global constant $α$ . As shown in Figure 6(b), for $α = 1$ , the field lines become slightly twisted. If the current is increased to $α = - 4$ in Figure 6(c), the shearing significantly grows and leads to a large twisting of the field lines. The amount of shearing is the same for currents that are antiparallel to the field [ $α = - 4$ in Figure 6(d)], but the field lines become sheared and are spiraling in the opposite direction. The structure becomes somewhat more complex for nonlinear force-free configurations (where $α$ changes in space), shown in Figure 6(e) and Figure 7(a).

•

The current is not parallel to the magnetic field and the Lorentz force $j \times B$ does not vanish. In this case, the Lorentz force is compensated by the plasma pressure gradient force and the configuration is magnetostatic [MHS, see Figure 6(f) and Figure 7(f)]. The effects are more visible in the 3D view in Figure 7. The Lorentz force leads to more vertical loops in the strong negative footpoint area (compared to potential and force-free fields), but the shearing and horizontal spiraling of the field lines are significantly lower than for the force-free field shown in Figure 7(a). The loops extend significantly higher into the atmosphere and, in the weak positive polarity footpoint area, they bend backward.

3.3.7 Linear MHS

A specific ansatz for the electric current density (see Low, 1991) is

\nabla \times B = α B + f (z) \nabla B_{z} \times e_{z},

(83)

which is a generalization of linear force-free fields in Equation 74 and the $α B$ term is the same as in linear force-free models. The second term on the right-hand side of Equation 83 is a special form of Equation 33; this part of the electric current density is perpendicular to $e_{z}$ . A popular choice is

f (z) = a \exp (- κ z),

(84)

where $a$ is a free parameter and $κ$ an inverse length. As pointed out by Aulanier et al. (1998), the form (Equation 83) causes a nonlinearity in the field-aligned current density because this component is caused by the $α B$ and additional parts of the horizontal current. Nevertheless, it is preferable to call the ansatz (Equation 74) linear magnetostatic because the corresponding MHS equations can be solved by Fourier transforms, similar to the linear force-free field equations in the section “Linear Force-Free Fields.” The shortcomings of linear models, as outlined in that section, apply here too. Both free parameters $α$ and $α$ have to be globally constant in the entire computational domain. This excludes strongly localized concentrations of electric current and Lorentz forces. The approach (Equation 83) has been applied also in spherical geometry (Bogdan & Low, 1986). In this case, the electric current density flows on spherical shells. An extension of the theory allows the addition of a field-aligned electric current (see Neukirch, 1995).

3.3.8 Nonlinear Force-Free Fields in 3D

In the generic case, there are no symmetries and nonlinear effects have to be taken into account. The system of equations (Equations 68 and 69 or, equivalently 70 and 71) must be solved. In principle, solutions using Euler potentials are possible but, due to the shortcomings and difficulties of Euler potentials, it is more common to avoid them in practical computations of solar coronal magnetic fields. The boundary conditions are deduced from measurements of the photospheric magnetic field vector in the solar photosphere. Either the three components of this vector $[B_{x} (x, y), B_{y} (x, y), B_{z} (x, y)]$ are used directly in Equations 68 and 69 or the distribution of the force-free function $α (x, y)$ is computed with the help of the vertical electric-current distribution $j$ (see Equation 77), giving

α (x, y) = \frac{μ_{0} j_{z} (x, y)}{B_{z} (x, y)} .

(85)

3.3.8.1 Solving a Well-Posed Problem: Grad–Rubin

It has been shown by Bineau (1972) that well-posed boundary conditions for nonlinear force-free fields are provided by the vertical magnetic field $B_{z} (x, y)$ and $α (x, y) [or j_{z} (x, y)]$ , which is required only on the surface region with either positive or negative magnetic polarity. The distribution of $α (x, y)$ needs to be prescribed only for one polarity because $α$ is constant along magnetic field lines. Therefore, a field line connecting positive and negative polarity regions cannot have different values of $α$ in both polarities. Existence and uniqueness of solutions based on these conditions have been proved under the condition of small and smoothly distributed $α$ . Solving Equations 70 and 71 is called the Grad–Rubin method (Grad & Rubin, 1958). Solving a mathematically well-posed problem is certainly tempting for many mathematicians and mathematical physicists, but the principal problem is that the measured horizontal magnetic field vector provides two distinct well-posed problems, namely the distribution of $α$ for the positive and negative polarity regions. For ideal measurements, both well-posed problems should lead to the same solution, but in reality, the horizontal magnetic field measurements contain measurement errors and inconsistencies and the two solutions can differ significantly (see Schrijver et al., 2008). For solar applications, the Grad–Rubin method is implemented as follows:

Compute a potential field with $B_{z} (x, y)$ as boundary condition, see the section “Potential Fields.”

Inject a current into the system by distributing $α$ along magnetic field lines via Equation 71.

These additional currents change the magnetic field, which is computed via Equation 70.

Naturally the field lines change, too, and the field-aligned currents are redirected in step 2.

Steps 2–4 are repeated until the solution converges. Convergence is not guaranteed, however.

3.3.8.2 Solving Force-Free Equations by Optimization or MHD Relaxation

Another possibility to solve the force-free equations iteratively is to use Equations 68 and 69 and the measured photospheric magnetic field vector $(B_{x} (x, y), B_{y} (x, y), B_{z} (x, y))$ as boundary conditions for the model. For this method, it is not necessary to compute the $α$ distribution. This is an advantage because $α$ can contain large errors, in particular in regions with low vertical field $B$ in the denominator of Equation 85. Providing all three components of the photospheric magnetic-field vector as boundary conditions leads to an overdetermined problem, which is a disadvantage of the method. For arbitrary distributions of the horizontal magnetic field vector, a solution of the force-free equations (Equations 68 and 69) is not guaranteed. It is therefore necessary to check if the observed boundary conditions are consistent with the assumption of a force-free field. For inconsistent boundary conditions, a preprocessing procedure can be applied which corrects the boundary field within the measurement errors of the horizontal fields in the photosphere. Codes that use this approach are the MHD relaxation and the optimization methods. They iterate for a force-free field as follows:

Compute a potential field with $B_{z} (x, y)$ as boundary condition, see the section “Potential Fields.”

Replace the horizontal fields in the bottom layer of the potential field by the measured (or preprocessed) horizontal field $B_{x} (x, y), B_{y} (x, y)$ .

This produces discontinuities at the lower boundary and a violation of the force-free Equations 68 and 69.

Iterate for a force-free field (solution of Equations 68 and 69) by solving

\frac{\partial B}{\partial t} = μ F

where $µ$ is an iteration parameter and $F$ is a function of the magnetic field and its derivatives (see Wiegelmann & Sakurai, 2012, section 6).

Repeat step 4 until the solution converges and $F$ vanishes.

Convergence is not guaranteed.

In the MHD relaxation method, $F$ is derived from simplified time-dependent MHD equations and in the optimization approach by minimizing a functional (see Wiegelmann, 2008; Wiegelmann & Sakurai, 2012, for review papers).

3.3.9 Application of Nonlinear Force-Free Models in 3D

Several studies use (nonlinear) force-free models to investigate structures on the Sun (e.g., using the 3D magnetic field line structure to guide the interpretation of coronal images). Also popular is to compute time series of force-free equilibria whereby the configuration is in equilibrium at each instance while the measured boundary conditions (photospheric magnetic field vector) change in time. Figure 8(a) shows one example of such a series of equilibria for a solar-active region. In Figure 8(b), the integrated free magnetic energy of the active region in time is compared with the monitored flaring activity. While the dynamic phase of a flare can obviously not be modeled by static equilibria, the free-energy content of these equilibria can tell whether there is energy available to power a flare (Jing et al., 2009). This process of accumulating free magnetic energy can be well modeled from a sequence of equilibria.

Open in new tab

Figure 8. (a) A nonlinear force-free field model for an active region observed on August 25, 2001, where the green lines are magnetic field lines. (b) Free magnetic energy computed from a sequence of nonlinear force-free models. The free energy (black diamonds) has been computed from the difference in energy of a nonlinear force-free model and a potential model for the same vertical magnetic flux in the photosphere. The integrated magnetic flux in the photosphere is shown in orange. The gray curves show the flaring activity, which resulted in an X-class flare (very large flare) at the peak of this curve between 16:23 and 16:45. During flares, free magnetic energy is converted into flaring energy. Interestingly, the energy already dropped about 15 min before the flare was detected (the drop is highlighted by the green line).

Adapted from Jing et al. (2009), part of Figures 4 and 5. Copyright AAS. Reproduced with permission.

3.3.10 Conclusions on Nonlinear Force-Free Models in 3D

Unfortunately, none of the methods (Grad–Rubin, MHD relaxation, optimization) can be guaranteed to converge toward a unique solution. If the codes are fed with consistent and complete boundary conditions, all methods can reconstruct the correct unique solution with rather low errors (as shown in Schrijver et al., 2006). For observational vector magnetograms, which contain measurement errors, inconsistencies, and even missing horizontal fields in some regions, the output of the codes can differ significantly, as seen in another code-comparison paper by De Rosa et al. (2009). Different implementations of these methods have incorporated certain algorithms to deal with noise and inconsistencies in the photospheric magnetic-field measurements. Since the launch of SDO in 2010, high-quality vector magnetograms have been routinely available and the quality of the 3D model fields have significantly improved. Furthermore, nonlinear force-free fields are routinely compared with coronal images. Some algorithms to compute nonlinear force-free fields incorporate information deduced from coronal images; for example, the projection of coronal loops or the 3D structure of stereoscopic reconstructed loops (see, e.g., Malanushenko et al., 2012; Aschwanden, 2013; Chifu et al., 2017). In these models, information from coronal images is already used during the iterative computation toward a force-free configuration. Consequently, nonlinear force-free field models have become a standard tool to guide the analysis of solar coronal observations. Despite the difficulties and challenges in computing nonlinear force-free fields, they provide more powerful and realistic field models than the potential and linear force-free models discussed in the sections “Potential Fields” and “Linear Force-Free Fields.”

3.3.11 Nonlinear MHS

As pointed out in “Nonlinear Force-Free Fields in 3D,” there can be inconsistencies between the measurements in the photosphere and the force-free model in the solar corona. These inconsistencies are not just due to noise, but are also attributed to the fact that in the lowest layers of the solar atmosphere (photosphere and chromosphere), magnetic and non-magnetic forces are equally important. Consequently, at a rather thin (about 2 Mm) layer above the height where the measurements are taken, the Lorentz force does not necessarily vanish and the force-free equations (68 and 69) are not fulfilled. Consequently, a step forward is to take non-magnetic forces such as pressure gradients and gravity into account and solve the full nonlinear MHS equations (Equations 1–3). Corresponding codes (Grad–Rubin, MHD relaxation, optimization) for solar applications are based on the corresponding nonlinear force-free models discussed in “Nonlinear Force-Free Fields in 3D.” These codes iterate self-consistently the 3D magnetic field, plasma pressure, and densities. While the corresponding implementations and first tests for such nonlinear MHS codes have been performed, they are in their infancy and still cannot be used as routinely as the corresponding nonlinear force-free models. Nonlinear MHS codes are an active area of research, and it is expected that they will be further developed and used more frequently beyond the present (2019) and into the future.

3.4 Rotating Planetary Magnetospheres and Rotating Magnetic Stars

3.4.1 Rotating Planetary Magnetospheres

For planetary magnetospheres, the gravitational potential of the planet can often be neglected, giving for such rigidly rotating systems without external gravity

j \times B = \nabla p + ρ \nabla V,

(87)

\nabla \times B = μ_{0} j,

(88)

\nabla \cdot B = 0,

(89)

where $V$ is the centrifugal potential for a magnetosphere in strict co-rotation. The potential is defined as

V = - \frac{Ω^{2} r^{2}}{2} .

(90)

A method to solve these equations analytically in cylindrical geometry $(r, ϕ, z)$ under special assumptions (see Neukirch, 2009) follows. As in the section “Linear MHS,” the electric current density is assumed to be perpendicular to the external force, namely the centrifugal force:

μ_{0} j = \nabla F \times \nabla V .

(91)

Open in new tab

Figure 9. MHD equilibria in rotating systems. (a) Dipolar solution for a rotating magnetosphere. Adapted from Neukirch (2009), Figure 1. Reproduced with permission from Taylor & Francis. (b) A superposition of dipolar and quadrupolar fields. Adapted from Wilson and Neukirch (2018), part of Figure 12. Reproduced with permission from Taylor & Francis. (c) Coronal magnetic field (displaced dipole) of a rotating star. Adapted from Al-Salti and Neukirch (2010), part of Figure 1. In all panels, the colored curves are magnetic field lines. Reproduced with permission from Quality Control ESO.

This kind of representation is reminiscent of the Euler potential representation (for the magnetic field $B$ ) introduced in the section “Euler Potentials” and is an example of the general formulation given in the section Potentials for the current density. Here, one of these potentials coincides with the centrifugal potential $V (r, ϕ, z)$ , while the potential $F (r, ϕ, z)$ remains a free function. Inserting Equation 91 into the force–balance Equation 87 gives

\frac{1}{μ_{0}} (B \cdot \nabla F) \nabla V - \frac{1}{μ_{0}} (B \cdot \nabla V) \nabla F = \nabla p + ρ \nabla V .

(92)

The pressure is a function of the potentials $p = p (F, V)$ and consequently $\nabla p = \partial p / \partial V \nabla V + \partial p / \partial F \nabla F$ . Because the potentials $F$ and $V$ are linearly independent, the terms containing $\nabla F$ and $\nabla V$ in Equation 92 have to be fulfilled independently, leading to

\frac{\partial p}{\partial F} = - \frac{1}{μ_{0}} (B \cdot \nabla V),

(93)

ρ = - \frac{\partial p}{\partial V} + \frac{1}{μ_{0}} (B \cdot \nabla F) .

(94)

This set of equations (93 and 94) is still difficult to solve for an arbitrary free potential $F (r, ϕ, z)$ . It can be linearized using a specific form (as suggested by Low, 1991; Neukirch, 2009),

F (r, ϕ, z) = κ (V) B \cdot \nabla V,

(95)

with a free scalar function $κ$ . Because of the linearity, new solutions can be created by the superposition of known solutions. Formulas for explicit analytical solutions can be found in Neukirch (2009). Figure 9(a) shows a dipolar solution. Subsequently, Wilson and Neukirch (2018) made use of the linearity of the solution and computed a configuration as the superposition of a dipolar and quadrupolar configuration, which is shown in Figure 9(b).

3.4.2 Rotating Magnetic Stars

Stellar applications usually include the gravitational field of the central body. For rotating configurations, including an external gravity field, one still has to solve Equations 87–89, but the potential $V$ as defined in Equation 90 includes now centrifugal and gravitational potentials. In the following, a spherical coordinate system $(r, θ, ϕ)$ is used and the combined potential (Al-Salti & Neukirch, 2010) is given as

V = - \frac{Ω^{2} r^{2} \sin^{2} θ}{2} - \frac{G M_{0}}{r},

(96)

Open in new tab

Figure 10. The cylindrical coordinate system ( $R, Φ, Z$ ) used for computing tokamak equilibria. $Φ$ is the invariant direction. The curves indicate a torus, which traps the plasma (indicated with dots).

Reproduced from Takeda and Tokuda (1991), Figure 2.1, with permission from Elsevier.

where $G$ is the gravitational constant and $M_{0}$ is the mass of the star.

Like Equation 91, the electric current density is described by two potentials. The equations are more complicated because $V$ is a function of $r$ and $θ$ and contains the gravitational potential (see Equation 96). Nevertheless, the same linearization as in Equation 95 is possible. It is not possible, however, to find analytical solutions by separation of variables, and a numerical approach is necessary. Figure 9(c) shows a configuration based on a displaced dipole found by Al-Salti and Neukirch (2010).

3.5 Laboratory Plasmas

An important application of MHD equilibrium theory for laboratory plasmas is the modeling of fusion devices, in particular of tokamaks and stellarators. Here only some basic fusion-relevant MHD equilibrium theory is reviewed.

Open in new tab

Figure 11. Different classes of topology for tokamak MHD equilibria. (a) Separatrix surface is outside the plasma. (b) Separatrix surface and surface of plasma coincide. (c) Separatrix surface is inside the plasma. The shaded areas indicate regions containing plasma and the thick black lines are the separatrix surfaces.

Reproduced from Takeda and Tokuda (1991), Figure 2.3, with permission from Elsevier.

For more details on MHD theory for a tokamak plasma, see Freidberg (1982). An overview of analytical and numerical approaches to compute MHD equilibria in 2D and 3D is given in Biskamp (1993, Chapter 3). The equilibrium state of fusion plasmas is modeled by the equations of magnetostatics (Equations 1–3). Fusion-relevant computations of MHD equilibria also have engineering applications. A realistic model of a fusion device can help in the design of the corresponding machines (e.g., tokamaks or stellarators) before they are built. Moreover, in fusion theory the gravitational term can be neglected and

j \times B = \nabla p,

(97)

\nabla \times B = μ_{0} j,

(98)

\nabla \cdot B = 0

(99)

must be solved. For a tokamak equilibrium, it is sufficient to solve the magnetostatic equations (Equations 97–99) in 2D due to its toroidal symmetry. This is not the case for the complex geometry of stellarators, which require the magnetostatic equations in 3D to be solved. The aim of this research is to trap plasma in the device. Therefore fusion devices with stable plasma equilibria are of interest. Stable means that plasma instabilities should not occur. The operation of tokamaks in a pulsed mode is not considered here.

3.5.1 2D Grad–Shafranov Equilibria With Application to Tokamaks

In a toroidally symmetric cylindrical geometry $(R, Φ, Z)$ (see Figure 10), the magnetostatic equations (Equations 97–99) can be reduced to a Grad–Shafranov equation (see, e.g., McCarthy, 1999):

- (\frac{\partial^{2} ψ}{\partial R^{2}} - \frac{1}{R} \frac{\partial ψ}{\partial R} + \frac{\partial^{2} ψ}{\partial Z^{2}}) = μ_{0} R j_{Φ} = μ_{0} R^{2} \frac{\partial p}{\partial ψ} + \frac{1}{2} \frac{\partial F^{2}}{\partial ψ},

(100)

where the configuration is 2D and invariant in $Φ$ . Here $ψ (R, Z)$ is the flux function and $j_{Φ}$ is the electric current density in the invariant toroidal direction $Φ$ . $p (ψ)$ and $F (ψ)$ are a priori unknown free source functions. Equation (100) is the cylindrical version of Equation 26, which contains both plasma pressure and a magnetic field component in the invariant direction. An overview of simple and sophisticated analytical solutions for computing tokamak-relevant MHD equilibria can be found in McCarthy (1999).

Key quantities to characterize the tokamak geometry are the inverse aspect ratio $ϵ$ and the safety factor $q$ , which are defined as follows:

ϵ = \frac{a}{R_{0}},

(101)

where $R_{0}$ and $a$ are the major and minor radius of the tokamak torus. The safety factor depends on the toroidal current or equivalently on the magnetic field. In the earliest tokamak designs, the safety factor $q$ was defined as the ratio of the numbers of times a (helical) field line travels around the (long) toroidal to the (short) poloidal direction. For a circular cross section, the Kruskal–Shafranov limit,

q \geq 1,

(102)

has to be fulfilled to ensure stability of the plasma configuration, but fusion machines are usually operated at higher values of $q$ . Here explicit expressions for $q$ are not given because they depend on the particular tokamak design (e.g., often D-shaped cross sections are used (see Takeda & Tokuda, 1991, section 2). An example of a modern tokamak currently under construction is ITER with the specification $R_{0} = 6.2 m$ , $a = 2.0 m$ , and $q = 3.0$ .¹

3.5.2 Analytical Solutions

3.5.2.1 Solov’ev Solution

Quite a simple analytical solution of Equation 100 has been found by Solov’ev (1968) using the assumption that the source functions $p = A / µ_{0} ψ$ and $F^{2} = 2 B ψ + F_{0}^{2}$ are linear in $ψ$ . Here $A, B, F_{0}$ are constants, which can be used to specify $j_{Φ}$ . With this approach, the right-hand side of Equation 100 does not depend explicitly on $ψ$ and the equation can be integrated:

ψ = ψ_{vac} - \frac{A R^{2}}{8} - \frac{B Z^{2}}{2},

(103)

where $ψ_{vac}$ is any current-free or vacuum solution, which can be superposed due to the linearity of the Grad–Shafranov equation under this assumption.

3.5.2.2 Quadratic and Dissimilar Source Functions

Analytical solutions of Equation 100 have also been found for the linear Grad–Shafranov equation. In this case, the source functions $p$ and $F^{2}$ are quadratic in $ψ$ (see McCarthy, 1999, section IIB). A wide and flexible class of analytical solutions can be constructed by choosing dissimilar source functions (e.g., with $p (ψ)$ linear and $F^{2} (ψ)$ a second-order polynomial in $ψ$ ). The corresponding solutions can be expressed in the form of Bessel functions and have free parameters that can be used to fit to experiments. This is usually achieved numerically by least-square fitting algorithms (see McCarthy, 1999, sections II and III).

3.6 Numerical Computations in 2D

For solving the cylindrical Grad–Shafranov Equation 100 in the generic nonlinear case (when the source functions are neither linear nor quadratic, nor a combination of both), one has to apply numerical methods. For details, see the review paper by Takeda and Tokuda (1991), “Computation of MHD equilibrium of tokamak plasma.” The paper also addresses briefly the computation of non-axisymmetric 3D equilibria, which are relevant for stellarators. The description of a numerical code for solving the Grad–Shafranov equation in toroidal (tokamak) geometry, with the help of a variational algorithm, can be found in Lütjens et al. (1996). The code was dubbed CHEASE (Cubic Hermite Element Axisymmetric Static Equilibrium) and is frequently used by several laboratories to model their tokamak devices.

In 2D configurations, magnetic field lines (magnetic surfaces in 3D) can be visualized as contours of the flux function, as shown in Figure 11. The three panels show different topologies of the magnetic surfaces, which are used for classifying the equilibria. In Figure 11(a), the separatrix surface is located outside the plasma region (gray). In Figure 11(b), the boundary of the plasma and the separatrix surface coincide, and in Figure 11(c), the separatrix surfaces lies inside the plasma region. For a tokamak equilibrium, one usually aims for equilibria of type (a), whereas cases (b) and (c) occur in specific devices (e.g., in divertor tokamaks and doublet tokamaks, respectively) (see Takeda & Tokuda, 1991, section 2, and references therein).

3.7 3D Numerical Solutions With Application to Stellarators

Due to its complex geometry, it is not possible to model stellarators in 2D and consequently one cannot reduce the magnetostatic equations (Equations 97–99) to a single partial differential equation. Suggestions on how fusion-relevant equilibria can be computed numerically have already been made by Grad and Rubin (1958). The suggested method, dubbed the Grad–Rubin method, has become popular to compute nonlinear force-free equilibria for astrophysical applications (see “Solar Magnetic Fields and Plasma”).

Chodura and Schlueter (1981) developed a code for computing 3D equilibria with the help of an energy principle in the form

W = \int (\frac{| B |^{2}}{2 μ_{0}} + \frac{p}{γ - 1}) d V,

(104)

where $γ \geq 0$ is the adiabatic index. As pointed out by Chodura and Schlueter (1981), minimizing the energy functional (Equation 104) should in principle lead to a stable equilibrium. Hirshman and Whitson (1983) developed a code using a steepest decent algorithm to minimize the functional (Equation 104) and thereby solved the magnetostatic equations (Equations 97–99) in 3D. The code VMEC (Variational Moments Equilibrium Code) has become a standard tool for computing stellarator equilibria and has also been applied to tokamaks. An overview of the VMEC code and a list of fusion devices for which it has been applied can be found in a fusion-wiki.²

4. Conclusion

This article has given an overview of the theory and application of magnetohydrodynamic (MHD) systems in equilibrium. Exact solutions are usually only possible if the number of dimensions is reduced to 2D or for special cases in 3D, where these special cases often include a linearization. While both approaches have been successful in studying physical systems such as magnetospheres and stellar atmospheres, they are not fully satisfying. The real world is three-dimensional and nonlinear. Within the past two decades, since the 1990s, the study of these generic cases (nonlinear and 3D) has become possible due to powerful computers and the development of sophisticated numerical techniques. In certain areas such as the solar corona, subsets of equilibrium theory (e.g., nonlinear force-free field models) are applicable because in the solar corona magnetic forces dominate over plasma forces. Despite a further simplification compared with magnetohydrostatic (MHS) models, developing optimal numerical codes to solve these equations is very challenging. Major problems have been and still are noise and inconsistencies of the measured boundary conditions (the magnetic field vector in the photosphere of the Sun and stars) and the additional complication that the Lorentz force vanishes only in the corona, but not necessarily in the lower solar (stellar) atmosphere (photosphere and chromosphere) and also not in the heliosphere and astrospheres (interaction regions of stellar winds with the interstellar medium), where the pressure gradient force, the external gravitational force, and the inertial forces caused by plasma flows (stellar winds) become important. For rapidly rotating stars and planetary magnetospheres, it is furthermore necessary to include the centrifugal force as an additional inertial force. Taking all these additional physical effects into account is a challenge. Popular approaches in 3D include a special class of solutions where the electric currents are restricted to flow perpendicular to external forces (like gravity, centrifugal force, or both). Such assumptions about the currents are mathematically convenient because they often allow a linearization of the system of equations, but unfortunately, this approach is often not physically justified. Furthermore, it has been demonstrated that this approach does not allow the inclusion of the centrifugal force and plasma flow on open field lines. To do so (without any assumptions on the currents), one has to use nonlinear, fully numerical approaches. Computer codes for nonlinear MHS equilibria in 3D are under development. The difficulties are the more complex physics in these approaches, but also the problem that additional measurements are needed to prescribe the boundary conditions. Challenges in the near future are also the inclusion of plasma flows (stellar winds) and centrifugal forces within a consistent model. Certainly, MHD equilibrium theory is not an isolated research topic, but embedded in the larger area of time-dependent MHD and non-MHD processes in plasma physics. Taking equilibrium solutions as initial states for time-dependent MHD simulations has been a standard approach in magnetospheric as well as in solar physics. The first approaches were in 2D, whereas now 3D simulations are state of the art. Care has to be taken that the initial MHD equilibria are not too simple because that can lead to misleading results. This can happen, for example, if finding numerical stability, when using a linear equilibrium as input, does not exclude instabilities (eruptions, dynamics) for nonlinear equilibria. Even more involved are sequences of MHD equilibria that form strong current concentrations and strongly localized Lorentz forces, because the assumptions of MHD are not satisfied within these structures and more sophisticated plasma models (e.g., kinetic theory) may need to be adopted.

Browse All

Share Link

Article contents

Magnetohydrodynamic Equilibria

Magnetohydrodynamic Equilibria

Summary

Keywords

Subjects

1. Introduction

2. Theory

2.1 From Magnetohydrodynamics to Equilibrium Theory

2.2 Magnetohydrostatics in 2D

2.2.1 Grad–Shafranov Equation

2.2.2 Liouville Equation

2.2.3 Harris Sheet

2.2.4 Including Magnetic Shear, 2.5D Configurations

2.2.5 Including Gravity

2.3 Magnetohydrostatics in 3D

2.3.1 Euler Potentials

2.3.2 Potentials for the Current Density

2.4 Stationary MHD With Plasma Flow (SMHD)

2.4.1.1 MHS as SMHD Subclass

2.4.1.2 SHD as SMHD Subclass

2.4.1.3 Mathematical Similarities Between SHD and MHS

2.4.1 SMHD Equilibria in 2D

3. Applications

3.1 Earth’s Magnetosphere

3.2 Magnetotail-Like Solutions of the Liouville Equation

3.2.1 Including Field-Aligned Plasma Flow

3.2.2 Sequences of Magnetostatic Equilibria

3.2.3 Magnetotail Convection Model

3.3 Solar Magnetic Fields and Plasma

3.3.1 Force-Free Magnetic Fields

3.3.2 Potential Fields

3.3.3 Linear Force-Free Fields

3.3.4 Nonlinear Force-Free Fields in 2D

3.3.5 Including Plasma Pressure: 2D MHS With Shear

3.3.6 Comparison of Different Models

3.3.7 Linear MHS

3.3.8 Nonlinear Force-Free Fields in 3D

3.3.8.1 Solving a Well-Posed Problem: Grad–Rubin

3.3.8.2 Solving Force-Free Equations by Optimization or MHD Relaxation

3.3.9 Application of Nonlinear Force-Free Models in 3D

3.3.10 Conclusions on Nonlinear Force-Free Models in 3D

3.3.11 Nonlinear MHS

3.4 Rotating Planetary Magnetospheres and Rotating Magnetic Stars

3.4.1 Rotating Planetary Magnetospheres

3.4.2 Rotating Magnetic Stars

3.5 Laboratory Plasmas

3.5.1 2D Grad–Shafranov Equilibria With Application to Tokamaks

3.5.2 Analytical Solutions

3.5.2.1 Solov’ev Solution

3.5.2.2 Quadratic and Dissimilar Source Functions

3.6 Numerical Computations in 2D

3.7 3D Numerical Solutions With Application to Stellarators

4. Conclusion

Further Reading

References

Notes

Related Articles