2 The Application of MHD Modeling to Magnetic Flux Emergence

The evolution of the buoyant rise of magnetic flux and its interaction with ambient flows is described by the magnetohydrodynamics (MHD) equations. Since the MHD equations are nonlinear, numerical simulations are often relied upon to study the evolution of the system. The derivation of the MHD equations (either from a kinetic or a fluid picture) is beyond the scope of this review and we refer the reader to texts such as Priest (1982*), Choudhuri (1998), Goedbloed and Poedts (2004), Parker (2007*) and Goedbloed et al. (2010).

The rest of this section provides an brief introduction to MHD and how choices of the material properties of the plasma lead to the diverse family of MHD models describing the physics of flux emergence.

2.1 MHD equations

Although we choose not to review the assumptions of MHD nor pursue a full derivation of the equations, it is nevertheless instructive to remind ourselves of the physical principles central to MHD. They are:

  • the Principle of mass conservation,
  • the Principle of momentum conservation,
  • the Principle of energy conservation, and
  • Faraday’s law of electromagnetic induction.

We present the MHD equations in the context of continuum mechanics, which is a general framework for the study of continuous media (i.e., at the macroscopic level) under various conditions. For most theories of continuum mechanics (e.g., Newtonian fluids, the study of elastic materials), the governing equations reflect the first three principles listed above. In the following, D ∕Dt denotes the Lagrangian derivative of a quantity. While the Eulerian version of the equations is often more convenient for implementation in numerical codes, the Lagrangian description has also been used extensively in the context of magnetic flux emergence. For instance, the Lagrangian MHD equations are the basis for the thin flux tube approximation, which model the evolution of discrete, individual flux tubes rising through the bulk of the solar convection zone. For details about the foundations of the thin flux tube approximation and its uses, we refer the reader to reviews by Moreno-Insertis (1997*) and Fan (2009b*). For this review, we prefer to present a Lagrangian description since it serves as the most useful framework for developing a physical picture of what happens to magnetic flux as it rises through the solar interior and eventually emerges into the atmosphere.

In a Lagrangian fluid description, the principle of mass conservation is

D-ϱ-+ ϱ∇ ⋅ v = 0, (1 ) Dt
where ϱ is the mass density and v the material velocity. This equation states that the mass density of a Lagrangian fluid element changes solely because of converging or diverging flow.

The principle of momentum conservation is expressed by Cauchy’s 1st law of motion

ϱ Dv--= ∇ ⋅ σ + f, (2 ) Dt
where σ is the stress tensor and f the body (volumetric) force. In most applications of solar MHD,
f = ϱg, (3 )
where g is the gravitational acceleration and
σ = − pδ + M , (4 ) ij ij2 ij M = − B--δ + BiBj-. (5 ) ij 8π ij 4 π
δ ij is the Kronecker-δ symbol and p denotes the isotropic gas pressure (usually as a function of other quantities). B is the magnetic field and M is the Maxwell stress tensor in cgs units. 1 ∇ ⋅ M = 4π (∇ × B ) × B represents the Lorentz force exerted by the magnetic field on the plasma, with the first term (− ∇B2 ∕8π) representing the magnetic pressure gradient. Together with the gas pressure gradient term, this term in the momentum equation is the reason emerging magnetic field structures expand as they reach tenuous layers of the solar atmosphere. The second term 1--∂-(B B ) = 1-B ∂Bi 4π∂xj i j 4π j ∂xj (since B is solenoidal) represents magnetic tension (or the magnetic curvature force), which tends to straighten magnetic field lines. The Lorentz force can be written in the more familiar looking form ∇ ⋅ M = 1cj × B, where j = c-∇ × B 4π is the current density and c is the speed of light. When viscous effects are important, the total stress tensor σ will also include a contribution from a viscous stress tensor.

The principle of energy balance can be stated in terms of the specific entropy s, so that

Ds Q ϱ--- = --, (6 ) Dt T
where T is the gas temperature and Q represents the volumetric heating rate. Given an appropriate equation of state relating the specific entropy with the specific internal energy 𝜀 and density ϱ, Eq. (6*) can be combined with the continuity equation (1*) to give the evolution equation for internal energy (see, e.g., Section 2.3 of Priest, 1982*), viz.
D 𝜀 ϱ --- = Q − p∇ ⋅ v, (7 ) Dt
where the gas pressure p is given by an appropriate equation of state (usually an ideal gas law for astrophysical plasmas). In the solar atmosphere, Q may include contributions from a non-vanishing divergence of conductive and radiative fluxes as well as local contributions from Joule and viscous dissipation.

Faraday’s induction equation is

∂B ----= − c∇ × E, (8 ) ∂t
where E is the electric field and c is the speed of light. In the regime of ideal MHD where the plasma is a perfect electrical conductor the electric field E′ in the co-moving inertial frame of the plasma vanishes. Assuming the plasma velocity v has speed |v| ≪ c, a Lorentz transformation to the ‘lab’ frame leads to
E = − c− 1v × B. (9 )
This yields the familiar Eulerian form of the ideal MHD induction equation
∂B-- ∂t = ∇ × (v × B). (10 )
In Lagrangian form, this equation becomes
DB ----= − B (∇ ⋅ v) + (B ⋅ ∇ )v. (11 ) Dt
This form of the ideal induction equation provides a simple, fluid-like picture of how the magnetic field evolves with the velocity field. The first term on the right says that a net expansion (compression) of the fluid element leads to a weakening (intensification) of the magnetic field on a Lagrangian fluid element. The second term is the gradient of the velocity field projected onto B, and says that a positive gradient of the flow along magnetic fields lines leads to intensification of the field (and vice versa). If the plasma is not perfectly conducting, additional terms appear in the expression for E. The nature of these terms depends on the model used for the describing the material properties of the plasma (see Sections 2.2.3 and 3.8).

The dot product of v with the momentum equation leads to an evolution equation for the kinetic energy density 1ϱv2 2. Similarly, the dot product of B with the induction equation leads to an evolution equation for the magnetic energy density 2 B ∕8 π. Combining these two equations with the equation for specific internal energy leads to an equation for the total energy density, (ϱ𝜀 + B2 ∕8π + 12ϱv2), as is commonly presented in other texts that introduce the MHD equations.

2.2 Constitutive relations

Equations (1*), (2*), (7*), and (8*) capture the physical principles of mass, momentum and energy balance as well as Faraday’s law of induction. They apply to a wide variety of astrophysical and laboratory conditions but these four equations alone do not comprise a complete set of equations that govern an MHD system. They must be supplemented by so-called constitutive relations, which are statements about the material properties of the plasma.

2.2.1 Equation of state

In a variety of MHD applications (both analytical and numerical), the plasma is assumed to follow adiabatic evolution. This should be recognized as an assumption (justified or otherwise) about the material properties of the plasma. When all of the following, including thermal conduction, radiative cooling/heating, Joule and viscous dissipation are neglected, it can be shown that the initial specific entropy of a fluid element is conserved (i.e., Ds ∕Dt = 0, see Section 24 of Mihalas and Weibel-Mihalas, 1984). In this regime, the plasma temperature and pressure change only due to compression or expansion (i.e., ∇ ⋅ v ⁄= 0).

An example of a constitutive relation is the stress tensor σ. Implicit in the form for σ in Eq. (4*) are assumptions or choices about how the plasma behaves. The existence of the term − pδij is a statement that there exists an isotropic gas pressure. The specific form of the isotropic gas pressure (usually called equation of state, EOS) is yet another constitutive relation. For a cold plasma, p = 0. As will be discussed in more detail in the following sections, a large number of numerical studies of flux emergence from the solar interior into the corona assume an EOS for a perfect, monatomic ideal gas such that the ratio of specific heats cp∕cv = γ = 5∕3. This assumes that the plasma is either completely neutral or completely ionized everywhere and at all times. In the case of plasma in the near-surface layers of the convection zone and overlying photosphere and chromosphere, the degree of ionization ranges from 10−5 to 10− 1 and the gas pressure is a function of both 𝜀 and ϱ (e.g., Stein and Nordlund, 1998*; Vögler et al., 2005*). In the photosphere and below, the assumption of local thermodynamic equilibrium (LTE) is appropriate and thermodynamic variables such as mass density, internal energy density and specific entropy are state variables that are functions of local gas temperature and pressure (and vice versa). In the tenuous chromosphere, the time scales for recombination of hydrogen ions with free electrons are comparable, or longer than the typical time between the passage of shock waves which ionize hydrogen atoms (Kneer, 1980*; Leenaarts et al., 2007*). As such, studies that aim to examine the thermodynamic response of the chromosphere in response to photospheric and coronal evolution will most likely need to solve the rate equations for the hydrogen atom (together with advection terms due to flows) to determine the local pressure and temperature.

2.2.2 Radiative and thermal conductive properties

In addition to the equation of state, the radiative and thermal conductive properties of the plasma are constitutive relations that must be appropriately chosen for the given problem. In the evolution equation for specific internal energy (𝜀, Eq. (7*)), the form of the volumetric heating/cooling term Q depends on these choices. For instance, in the flux emergence model of Shibata et al. (1990a*), a radiative term with a specified cooling time scale was imposed in the photospheric layers to mimic the cooling of photospheric plasma by radiation to the higher atmospheric layers. This simple treatment was sufficient to demonstrate that such cooling would lead to the developing of downflows that convectively intensify emerged flux from a few hundred G strength to beyond 1 kG (see Abbett, 2007*; Fang et al., 2010*, for a more sophisticated parametric treatment of radiatively cooling in the solar atmosphere). For studies that attempt to realistically model the temperature structure of the photosphere, 3D radiative transfer calculations must be carried out to solve for the radiative flux Frad (e.g., Stein and Nordlund, 2006*; Cheung et al., 2007a*; Yelles Chaouche et al., 2009; Martínez-Sykora et al., 2008*, 2009*; Tortosa-Andreu and Moreno-Insertis, 2009*). In such cases, Q is dominated by the term ∇ ⋅ Frad in layers where the plasma transitions from being optically thick to optically thin.

The choice of the thermal conductivity of the plasma is also a constitutive relation. In the solar corona, electrons are strongly magnetized and the transport coefficients are such that thermal conduction is predominantly aligned along magnetic field lines. The presence of such field-aligned thermal conduction leads to efficient energy transport from a reconnection region in the corona to chromospheric footpoints of the field. This can lead to chromospheric evaporation jets which feed mass into the corona (e.g., Yokoyama and Shibata, 2001*). This topic will be discussed in more detail in Section 3.7.2.

2.2.3 Ohm’s law

The form of the electric field given by Eq. (9*) assumes that the plasma is a perfect electrical conductor. When the electrical conductivity of the plasma is finite, non-ideal terms will appear in the expression for E. If a scalar electrical conductivity σc is assumed (i.e., a scalar Ohm’s law), the corresponding electric field is

cE = − v × B + η∇ × B, (12 )
where η = c2∕4π σ c is the magnetic diffusivity of the plasma.1 The additional non-ideal term captures the effect of Ohmic diffusion of the magnetic field, which is the key ingredient for magnetic reconnection. The majority of non-ideal solar MHD models use this form of the electric field in the induction equation. Since the estimated magnetic diffusivity of solar plasma is often much too small for numerical simulations to resolve structures at the diffusive scale, modelers often adopt spatially varying forms of η such that diffusive effects are minimal outside of reconnection layers and sufficient within these layers permit reconnection without introducing spurious numerical effects. One example of such a treatment is the so-call anomalous resistivity (Yokoyama and Shibata, 1994*). This model for the resistivity assumes η ∝ j2, where j is the magnitude of the current density. This dependence was found to permit the onset of Petchek-type fast magnetic reconnection in simulations of emerging flux interacting with pre-existing coronal field. Other choices of the imposed magnetic resistivity such as the so-called hyper-diffusivity scheme (Caunt and Korpi, 2001) are invoked to do a similar job as anomalous resistivity, namely to permit magnetic diffusion where it is most needed. This type of spatially-varying magnetic diffusivity appears to be an important ingredient in models that create fast reconnection jets and plasmoid ejections that accompany flux emergence into the solar corona (see Section 3.7.1 for a discussion how resistivity models affect models of flux emergence).

In weakly ionized plasmas, interactions between neutral and charged species leads to an electric field with additional terms describing ambipolar diffusion (also called the Pedersen effect) and the Hall effect. The inclusion of these effects can lead to evolution of the emerging magnetic field that is distinctly different from more traditional MHD models that simply use Eq. (12*) for the electric field. We will discuss this in more detail in Section 3.8.

2.3 A diversity of MHD models

Not all MHD models are created equal. While the basic principles of physics captured by MHD are general, the freedom for the modeler to choose the constitutive relations that describe the material properties of the plasma is the main reason for the diversity the MHD models that exist for modeling flux emergence.

View Image
Figure 5: Models of magnetic flux emergence can be roughly divided into three categories, though there are large areas of overlap between them. So-called ‘realistic’ models attempt to include all the known important physical ingredients, while idealize models generally focus on studying a more limited set of effects. For case studies of certain observed emerging flux studies, data-driven models are used.

Figure 5* lists some of the choices a modeler makes before performing a theoretical/numerical study of the flux emergence process. Some models are set up for studying certain physical effects in idealized conditions, while other models attempt to treat all (known) relevant effects at once. Some are concerned with interactions at certain length-scales (e.g., scale of granulation) while others try to capture entire active region scales. Some models use magnetic field configurations that mimic specific instances of observed episodes of emerging flux episodes, while others are not concerned with reproducing the observed patterns with high fidelity. As displayed in this figure, MHD models of flux emergence can be roughly divided into three categories: idealized, ‘realistic’, and data-driven models. These categories need not be mutually exclusive. While the division into the three categories is by no means unique, it serves as a guide for us to navigate the literature on MHD modeling of flux emergence.

Idealized models are so-called because they simplify the problem by choosing to neglect only certain effects. There are a number of advantages to idealized models. By choosing to neglect certain effects, idealized models often have simpler setups that allow the modeler to studying other physical effects in isolation. For example, many MHD models use plane parallel polytropes to mimic the average stratification of the solar atmosphere without including the necessary physics that lead to self-consistently generated convective motions and magneto-acoustic waves that channel energy into the chromosphere and corona. Such a choice allows modelers to concentrate on studying effects such as magnetic Rayleigh-Taylor instabilities (see Section 3.3) and the expansion of magnetic flux as it emerges from the convection zone into the tenuous layers of the atmosphere. Another benefit of idealized models is that they are often computationally less demanding and allow for a more extensive exploration of parameter space. In the context of flux emergence, exploration of parameter space includes variations in the flux content or twist of an emerging flux tube, the initial configuration of any pre-existing magnetic field in the atmosphere, the number of emerging flux tubes etc. In addition, simplifications in the model often allow for larger computational domains, finer resolution, and/or longer simulation times (in terms of characteristic timescales of the system).

In contrast to idealized models, so-called realistic models attempt to capture as much as possible all physical processes that are known to be important for dynamics in the solar atmosphere and convection zone, as well as crucial for synthetic diagnostics. The most developed of this class of models in the flux emergence context is from the work of Martínez-Sykora et al. (2008*, 2009*). Their model of flux emergence employs a computational domain that captures the top few Mm of the convection zone as well as the photosphere, chromosphere, and corona. Convective flows in the convection zone are driven by radiatively cooling at the surface that results from solving the 3D radiative transfer problem and a realistic equation of state is used to account for changes in ionization degree in the plasma. These effects are also captured in the flux emergence models of Cheung et al. (2007a*), Tortosa-Andreu and Moreno-Insertis (2009*), and Stein et al. (2011*). However, the Martínez-Sykora et al.* model also includes a chromosphere treated with radiative transfer that includes scattering (i.e., the source function is a linear combination of the Planck function and the local radiation field) and coronal physics such as field-aligned thermal conductive and optically thin radiative losses. This allows them to compute synthetic diagnostics from their flux emergence simulations for the photosphere, chromosphere (e.g., see Figure 6*), and corona, and therefore to compare with observations of all these layers. One drawback of realistic models is that they are often too prohibitively expensive for extensive parameter studies. The feedback between the many physical processes also make it a challenge to distill general physical insight. For these reasons, idealized studies remain essential for providing complementary guidance.

View Image
Figure 6: Synthetic Hinode/SOT Ca ii H image from the flux emergence simulation by Martínez-Sykora et al. (2009*). The line-of-sight is taken to be parallel to the solar surface to mimic observations of such regions at the solar limb. Image reproduced with permission, copyright by AAS.

  Go to previous page Scroll to top Go to next page