1 Introduction

The study of how magnetic flux emerges onto the Sun is motivated by a number of fundamental questions about solar magnetic activity. First of all, a complete picture of the solar dynamo requires not only an understanding of how magnetic fields are generated and amplified in the solar interior, it also requires an understanding of the transport processes that bring magnetic fields to the solar atmosphere. Observations of how active regions and smaller scale flux emergence events appear give telltale signs of the structure of magnetic fields below the surface. However, the interpretation of what the observed patterns of magnetic activity mean requires us to know how they got there in the first place. Are buoyant flux bundles that form active regions twisted, and what does that mean for the solar dynamo? How deep are sunspots anchored, and do their formation and decay influence the activity cycle?

In addition to its relation to the solar dynamo problem, the study of magnetic flux emergence helps us understand the basic physical mechanisms that are responsible for a myriad of dynamical atmospheric phenomena. For instance, emerging magnetic flux can help energize the solar atmosphere by building up free magnetic energy, or it lowers it by triggering eruptive events such as flares, coronal mass ejections (CMEs), and jets. The interaction of emerging flux with pre-existing magnetic field in the atmosphere can lead to the creation of current sheets, and magnetic reconnection in these regions can lead to impulsive releases of energy. The basic physical processes that control the rate of magnetic reconnection and the means by which magnetic energy is released into other forms of energy (kinetic, thermal, and radiative) are of interest for the larger astrophysical and plasma physics community.

Magnetic flux emergence is a vibrant and active area of research in solar physics driven forward partly by advances in observational capabilities of the solar atmosphere (and interior) and partly by the availability of supercomputers that facilitate large-scale numerical simulations. For the uninitiated, an attempt to digest the diverse range of studies (both theoretical and observational) in the flux emergence literature can be a daunting task. We hope to make this task easier. This review article aims to provide

  1. 1.

    a primer on the physics that govern the behavior of emerging magnetic flux,

  2. 2.

    an introduction of how magnetic flux emergence plays a key role in many aspects of solar physics, and

  3. 3.

    a broad overview of both established and recent developments in this very active field of research.

The main scope of this review is to cover theoretical aspects of magnetic flux emergence. We hope that, upon exposure to the information contained in this article, the reader will be sufficiently familiar with the key physical concepts that they can become discerning readers of the flux emergence literature. Readers interested in observational studies of flux emergence will find relevant references. However, there are simply too many studies to be included in a review of this limited scope (which is focused on theory).

1.1 Structure of this review article

The article is structured as follows. Section 1.2 expands on the introduction and discusses how flux emergence fits in with a broad range of problems in solar physics.

Section 2 provides a brief discussion of magnetohydrodynamics (MHD), which is the most useful theoretical framework for studying magnetic flux emergence. Section 2.1 introduces the basic MHD equations capturing conservation principles and Faraday’s law of electromagnetic induction. Section 2.2 discusses how constitutive relations describing the material properties of the plasma are used to supplement the MHD equations. The science questions that are being addressed and the choice of constitutive relations largely determines the differences between different models.

Section 3 details how certain physical mechanisms and processes play key roles in the progression of magnetic flux emergence.

Section 4 discusses the relation between flux emergence and eruptive events such as jets and CMEs.

Section 5 discusses how data-driven models are starting to be used in the study of flux emergence.

Finally, Section 6 will give some thoughts on open questions.

1.2 Science questions and challenges in flux emergence

Emerging flux plays a key role in many aspects of solar magnetic activity. The role of theories and numerical modeling of emerging flux is to clarify the various physical processes involved in flux emergence and therefore contribute to a better understanding of the solar magnetic activity and ultimately improve our ability to forecast hazardous space weather events.

Specific science questions in the modeling of emerging flux include:

  • What is the mechanism that brings the magnetic flux from the interior to the atmosphere?

  • How does emerging flux transport the magnetic energy and helicity?

  • How much flux is trapped below the surface during flux emergence, and what is the contribution of this trapped flux to the solar dynamo?

  • What are the roles of emerging flux in the free energy accumulation and triggering of the transient events such as jets, flares, CMEs?

  • What are the physical properties of subsurface magnetic structures that rise and eventually emerge onto the surface?

  • Does flux emergence occur as the rise of coherent bundles or as smaller elementary units? (Zwaan, 1978, 1985)?

  • What is the difference between flare-productive sunspots and quiet sunspots?

  • How do convective flows impact the morphology and physical character of emerging flux (e.g., Fan et al., 2003; Cheung et al., 2007a)?

  • Is magnetic twist necessary for magnetic flux to emerge (e.g., Moreno-Insertis and Emonet, 1996; Murray et al., 2006)?

  • How do the individual and statistical properties (e.g., Hagenaar, 2001; Hagenaar et al., 2003; Iida et al., 2012; Otsuji et al., 2011) of emerging flux relate to the solar activity cycle?

  • What are the observational consequences of emerging magnetic flux (Bruzek, 1967; Zwaan, 1978), and what are the physical mechanisms responsible?

  • Can we predict the appearance of new emerging flux region?

  • What are the physical ingredients necessary for a realistic model of emerging flux?

Figures 13 show a selection of observations of emerging flux regions covering the photosphere to corona. They illustrate the range in temperature, density, and plasma-β (ratio of gas to magnetic pressures) found in emerging flux regions. By definition, flux emergence consists of magnetic field reaching the atmospheric layers from the solar interior. Models that describe this evolution need to take into accounts physical effects that are important for the convection zone, the photosphere, the chromosphere, transition region, and eventually the corona. Different models attempt to include the relevant effects to varying degrees of comprehensiveness. The choice depends on the science question, but all models face the challenge of capturing the wide range of length and time scales relevant for flux emergence. Whereas active regions are born over the course of days and may spawn eruptions throughout their lifetimes (e.g., see Figure 3), Alfvénic crossing times over active region loops can be of the order of minutes or less. While G-band bright points can be just a few tens of km across in diameter (see Figure 1), sunspots can have diameters exceeding 10 Mm. The abrupt changes in physical regimes also pose challenges. For instance, whereas radiative cooling can be approximated by optically thin losses in the corona, the full radiative transfer equations need to be solved for the photosphere and chromosphere (in order to have thermodynamic structures that comparable to observations). The plasma-β is greater than unity in the convection zone, but tiny in the corona (where the magnetic field is dominant).

Figure 1:
figure 1

Observation of a small-scale emerging flux event within an active region. The left two panels show the G-band intensity and Stokes V signal from the photosphere. The right two panels show the chromospheric intensity in Hα and Ca II H. Image reproduced with permission from Guglielmino et al. (2010), copyright by AAS.

Figure 2:
figure 2

Observations of an emerging flux region by SDO/HMI and SDO/AIA. Images from AIA show the response of the corona to the birth of an active region at the solar surface. Image reproduced with permission from Centeno (2012), copyright by AAS.

Figure 3:
figure 3

Hinode/SOT CaII H observation of NOAA AR 11158, which produced a series of M and X-class flares as it emerged into the solar atmosphere. This particular image shows the flare ribbons and post-flare loops after the X-class flare on Feb. 15, 2011 (Image credit: Joten Okamoto).

Whether during solar maximum or minimum, the Sun’s surface is pervaded by magnetic fields at all scales. High-sensitivity measurements of the photospheric magnetic field show this is true even in so-called ‘quiet-Sun’ regions (Lites et al., 1996; Domínguez Cerdeña et al., 2003; Harvey et al., 2007; Orozco Suárez et al., 2007; Lites et al., 2008; Pietarila Graham et al., 2009; Danilovic et al., 2010b,a). Turbulent motions of the plasma from scales beyond supergranulation (L ∼ 30 Mm) down to scales below granulation (L ∼ 1 Mm) weave the magnetic field into complex configurations that lead to the formation of current sheets that facilitate magnetic reconnection. When Parker (1988) put forth his nanoflare model of coronal heating based on the dissipation of current sheets in the corona, he considered continuous horizontal photospheric motions of line-tied magnetic fields as the source of energy build-up in the corona. This by itself is considered sufficient for the production of current sheets, which then dissipate via magnetic reconnection. When emerging flux is introduced into this picture (even if the problem of interest is not coronal heating), one can expect magnetic reconnection to be even more pervasive. This is due to the fact that magnetic fields emerging through the convection zone (a low plasma-β environment) are not constrained to emerge into the atmosphere with an orientation that is aligned with pre-existing field (e.g., see Figure 4). As demonstrated in several models of flux emergence (see Section 4), this misalignment will create current sheets that are amenable to reconnection. So studies of flux emergence are also relevant for magnetic reconnection research, and the fact that reconnection can occur in different plasma regimes in the solar atmosphere (see Section 3.7.1) means flux emergence events are natural experiments for magnetic reconnection. One of the challenges of flux emergence research is to create models with reconnection behavior that is justified by the underlying microphysics and matches observed behavior at large (observable) scales.

Figure 4:
figure 4

SDO/AIA observations of NOAA AR 11112 in UV and EUV channels. Tarr et al. (2014) estimate that, over the course of two days, the quantity of (steadily) reconnected flux between the emerging flux region and preexisting field is comparable to a large M- or a small X-class flare. Image courtesy of L. A. Tarr.

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\varrho } \over {Dt}} + \varrho \nabla \cdot{\bf{v}} = 0,$$
(1)

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

$$\varrho {{D{\bf{v}}} \over {Dt}} = \nabla \cdot\sigma + {\rm{f}},$$
(2)

where σ is the stress tensor and f the body (volumetric) force. In most applications of solar MHD,

$${\bf{f}} = \varrho {\bf{g}},$$
(3)

where g is the gravitational acceleration and

$${\sigma _{ij}} = - p{\delta _{ij}} + {M_{ij}},$$
(4)
$${M_{ij}} = - {{{B^2}} \over {8\pi }}{\delta _{ij}} + {{{B_i}{B_j}} \over {4\pi }}.$$
(5)

δ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. \({\bf{M = }}{1 \over {4\pi }}(\;\;{\bf{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 \over {4\pi }}{\partial \over {\partial {x_j}}}({B_i}{B_j}) = {1 \over {4\pi }}{B_j}{{\partial {B_i}} \over {\partial {x_j}}}\) (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 \({\bf{M = }}{1 \over c}{\bf{J}}\;\;{\bf{B}}\), where \({\bf{j}} = {c \over {4\pi }}\) 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

$$\varrho {{Ds} \over {Dt}} = {Q \over T},$$
(6)

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.

$$\varrho {{D\varepsilon } \over {Dt}} = Q - p\nabla \cdot{\bf{v}},$$
(7)

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

$${{\partial {\bf{B}}} \over {\partial t}} = - c\nabla \times {\bf{E}},$$
(8)

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

$${\bf{E}} = - {c^{ - 1}}{\bf{v}} \times {\bf{B}}.$$
(9)

This yields the familiar Eulerian form of the ideal MHD induction equation

$${{\partial {\bf{B}}} \over {\partial t}} = \nabla \times ({\bf{v}} \times {\bf{B}}).$$
(10)

In Lagrangian form, this equation becomes

$${{D{\bf{B}}} \over {Dt}} = - {\bf{B}}(\nabla \cdot{\bf{v}}) + ({\bf{B}}\cdot\nabla){\bf{v}}.$$
(11)

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 \over 2}\varrho {v^2}\). Similarly, the dot product of B with the induction equation leads to an evolution equation for the magnetic energy density B2/8π. Combining these two equations with the equation for specific internal energy leads to an equation for the total energy density, \((\varrho \varepsilon + {B^2}/8\pi + {1 \over 2}\varrho {v^2})\), 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 −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 e and g (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 (e, 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

$$c{\bf{E}} = - {\bf{v}} \times {\bf{B}} + \eta \nabla \times {\bf{B}},$$
(12)

where η = c2/4πσc is the magnetic diffusivity of the plasma.Footnote 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.

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.

Figure 5:
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.

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.

Figure 6:
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.

3 The Physics of Magnetic Flux Emergence

3.1 Magnetic buoyancy

Magnetic flux emergence studies begin with the premise that the appearance of bipolar magnetic flux structures on the solar surface results from the buoyant rise of magnetic plasma from the convection zone into the overlying atmosphere. This mechanism was first invoked by Parker (1955) to explain the formation of active regions and sunspots. For an detailed and broad account of studies that examine the origin and transport of buoyant magnetic structures in the solar interior, we refer the reader to the articles by Moreno-Insertis (1997) and Fan (2009b).

As an introduction to the concept of magnetic buoyancy, consider a magnetic structure embedded in the solar convection zone with magnetic field strength B. Define β = 8πp/B2 (i.e., the ratio of gas and magnetic pressures). Assume the magnetic field is sufficiently weak such that B2/8πp (i.e., high plasma-β). In the case of total pressure equilibrium between plasma within the magnetic structure and the ambient plasma,

$${p_i} + {{{B^2}} \over {8\pi }} = {p_{{\rm{amb}}}},$$
(13)

where pi and pamb denote the internal and ambient gas pressures, respectively. If the structure were in thermal equilibrium with its surroundings (i.e., Te = Tamb), the internal mass density would be lower than the ambient value by a relative amount

$$\Delta \varrho / \varrho \approx - {\beta ^{ - 1}}.$$
(14)

If instead of thermal equilibrium, the interior and the external fluids possessed the same specific entropy s, the density deficit would be

$$\Delta \varrho / \varrho \approx - {({\gamma _1}\beta)^{ - 1}},{\rm{where}}\,{\gamma _1} = {\left({{{\partial \ln p} \over {\partial \ln \varrho }}} \right)_s}$$
(15)

is Chandrasekhar’s first adiabatic index (Chandrasekhar, 1967). In either case, the density deficit leads to an upward directed buoyancy force of ∇ϱg, which would lead to the ascent of the magnetic structure.

3.1.1 Introducing buoyant magnetic structures in simulations of emerging flux

In numerical experiments simulating the rise of buoyant magnetic structures (tubes, sheets, etc), it is often desirable to begin with an initial condition such that the forces acting on the magnetized plasma vanish with the exception of buoyancy. This is desirable since any unbalanced pressure gradient or Lorentz force will lead to an evolutionary adjustment (i.e., compression, expansion or distortion) of said structure over a sound-crossing (τSL/cs, where L is a typical length-scale) or Alfvèn crossing time (τAL/CA). Depending on the choice of parameters, τS and τA can be much smaller than the rise time due to buoyancy. So an out-of-equilibrium initial condition would be dramatically affected before the dynamics of interest take place. This is undesirable, especially for numerical experiments that attempt to isolate various physical effects by means of parameter studies. These considerations are also important in studies that examine the development of magnetic buoyancy instabilities, which require equilibrium background states. Since such setups are rather important, this section will provide some details on how equilibria are typically constructed.

Consider a rather common setup in the flux emergence literature, in which an initial twisted magnetic flux tube is embedded in a plane-parallel stratification. A plane-parallel stratification is one such that the MHD quantities are only functions of height z. The initial stratification is chosen such that the right hand side of Eq. (2) is zero (i.e., ∇ · σ + gϱ = 0). For plane-parallel, non-magnetic stratifications, this reduces to \({{dp} \over {dz}} = g\varrho\), the hydrostatic condition. In general, the introduction of a magnetic flux tube into the atmosphere will lead to a deviation from mechanical equilibrium since the Lorentz force associated with the newly introduced magnetic field manifests as both a magnetic pressure gradient force and an magnetic tension force (see Section 2.1).Footnote 2 To restore mechanical equilibrium, the gas pressure distribution within the tube must be decreased appropriately. The decrease in gas pressure can be obtained via a combination of changes in gas temperature and density (or any other pair of independent state variables). In the context of flux emergence studies, the choice is usually made that results in a density deficit. This means hydrostatic equilibrium is no longer preserved and the tube will experience a buoyant acceleration of \({{D{v_z}} \over {Dt}} = g\Delta \varrho / \varrho\) It is common in a number of 3D simulations to find initially horizontal tubes with a density profile that varies along the tube axis. This is often done to induce the development of an Ω-loop.

While a large fraction of numerical experiments of magnetic flux emergence begin with an initially buoyant magnetic structure embedded within a stratified layer, there are also many exceptions. For instance, while the simulation of Magara and Longcope (2001) began with an initially horizontal, twisted flux tube, the initial rise of the tube was not due to magnetic buoyancy (plasma in the tube had the same density as the exterior). Instead, vertical momentum was artificially imparted to a segment of the tube to initiate the development of an Ω-loop.

In some simulations, the magnetic structure is initially not even embedded in the computational domain. Rather than inserting a buoyant magnetic structure in the model convection zone, some models introduce rising magnetic flux by using time-dependent boundary conditions at the bottom boundary to inject flux into the computational domain. For instance, in a series of realistic magnetoconvection simulations (Stein and Nordlund, 2006; Stein et al., 2011; Stein and Nordlund, 2012), Stein and collaborators injected magnetic field into the computational domain (which captures the photosphere and the top 20 Mm of the convection zone) by imposing horizontal magnetic field of uniform orientation and strength in convective upflows crossing the bottom boundary. Similarly, the flux emergence and active region formation simulations of Cheung et al. (2010) introduce a semi-toroidal twisted flux tube by injecting the structure through the bottom boundary. While the injected magnetic structure already has vertical momentum as a result of being embedded in upflows, they are still introduced in a way that results in a density deficit relative to the mean stratification (even non-magnetic upflows have this property), which means that such structures still experience buoyancy forces.

3.2 Effect of stratification in the convection zone

The solar convection zone is strongly stratified and the density contrast between the bottom of the convection zone (at a depth of about z = 200 Mm) and the surface is roughly 106 (with ϱbottom ∼ 10−1 g cm−3 and ϱτ =1 ∼ 10−7 g cm−3, see Spruit, 1974; Christensen-Dalsgaard et al., 1991). Let p be the mean gas pressure as a function of depth in the convection zone. At the top of the convection zone, the pressure scale height, defined as

$${H_p} = {\left({{{d\ln p} \over {dz}}} \right)^{ - 1}},$$
(16)

is roughly 150 km and increases with depth. With this in mind, consider a fluid element originally at pressure equilibrium with its surroundings at z = 0 and rising one pressure scale height. If it maintained its initial pressure and density throughout the ascent, it would have a pressure enhancement of order unity relative to its new surroundings. Such a pressure enhancement cannot be maintained longer than a sound-crossing time. For a granule with radius 500 km and assuming a sound speed of 7 km s−1, this is roughly 1 min. Since they are subsonic, with speeds on the order of 0.1–1 km s−1, granular upflows must expand horizontally to maintain rough pressure equilibrium with its surroundings. Based on this simple concept of lateral pressure equilibrium, Spruit et al. (1987) predicted that a buoyant flux tube rising through the convection zone will flatten into a sheet-like structure as it approaches the base of the photosphere (see Figure 7).

Figure 7:
figure 7

Figure from Spruit et al. (1987) illustrating the conjectured behavior of a buoyant flux tube as it approaches the photospheric base (indicated by the flat horizontal line) from below. The arrows indicate the direction of plasma motion. The severe change in aspect ratio of the tube’s cross-section is due to the diminishing pressure scale height at the top of the convection zone. This predicted behavior is a robust result in numerical MHD simulations of emerging flux. Image reproduced with permission, copyright by Springer.

Figure 8 shows a sequence from a 2D simulation of the rise of a buoyant magnetic flux tube over multiple pressure scale heights. In this model by Toriumi and Yokoyama (2011), the convection zone is modeled by an adiabatically stratified polytrope spanning more than 20 pressure scale heights. The initial depth of the flux tube is at a depth of 20 Mm. In the initial rise phase (t < 400 in dimensionless units), the evolution of the flux tube is consistent with results from earlier work modeling the rise of 2D flux tubes, the likes of which demonstrate how the horizontal pressure gradient results in a tendency for the tube to split into counter-rotating vortex rolls (Schüssler, 1979; Longcope et al., 1996). The flux tube in this case retains coherence because of its transverse component, which imparts a magnetic tension which counteracts the horizontal pressure gradient (Moreno-Insertis and Emonet, 1996; Emonet and Moreno-Insertis, 1998, see Section 3.6 for an extended discussion about the role of magnetic twist). As the buoyant rise proceeds, the flux tube impinges into the near-surface layers. Due to the diminishing scale heights, the flux tube flattens into a pancake-like structure, as predicted by Spruit et al. (1987).

Figure 8:
figure 8

Sequence of snapshots from a 2D MHD simulation of the buoyant of a twisted buoyant flux tube through a polytropic model of the convection zone. The left and right panels show the magnetic field strength and out-of-plane component of vorticity, respectively. The top of the convection zone in this model is at z = 0. Image reproduced with permission from Toriumi and Yokoyama (2011), copyright by AAS.

Radiative MHD simulations by Cheung et al. (2010) show that a similar flattening and horizontal expansion occurs even in the presence of convective flows driven by photospheric cooling (see Figure 9). In their simulation of the rise of a twisted semi-torus, the drop in ambient pressure lead to strong horizontal expansion of the apex of the rising torus. This expansion filled the near-surface layers of the convection zone with dispersed magnetic flux.

Figure 9:
figure 9

The top panel shows the photospheric distribution of the vertical component of the magnetic field resulting the simulated rise of a semi-toroidal flux tube from at depth of 7.5 Mm to the surface. The bottom panel shows a vertical cross-section (at y = 0) of the magnetic field strength in the photosphere and in the subsurface layers captured by the simulation. Due to the density and pressure drop over this height range, the emerged flux is dispersed and covers a much larger area that the subsurface roots of the torus. Image reproduced with permission from Cheung et al. (2010), copyright by AAS.

3.2.1 Scaling relation between density and magnetic field strength

It is worthwhile to examine possible scaling laws for how the field strength of a rising magnetic structure changes with decreasing depth in the convection zone. For this purpose, one is interested in following the change of field strength B along the trajectory of a fluid element. The simplest case is illustrated in panel (a) of Figure 10, which shows a horizontal flux tube that has risen and expanded. The tube rises as a whole and does not have motion (more precisely, flow gradients) along its axis. By virtue of mass and magnetic flux conservation, it is easy to show that the magnetic field strength in the tube will scale as

$$B \propto \varrho.$$
(17)
Figure 10:
figure 10

Two different scenarios considered for deriving the scaling relation between density and magnetic field strength in Section 3.2.1. Panel (a) illustrates the case where a horizontal flux tube rises uniformly and expands only in the directions perpendicular to the tube axis. Panel (b) shows a scenario where a localized segment of a flux tube rises and expands (‘magnetic bubble’ case). Due to the diminishing scale heights near the surface, it has a flattened structure (cf. Figures 7 and 8).

Panel (b) shows a case where a localized segment of the tube is approaching the near-surface layers of the convection zone and experienced strong expansion due to the diminishing pressure scale height. To derive the corresponding scaling relation, consider Eqs. (1) and (11), which are the continuity and ideal induction equations in Lagrangian form. Following Cheung et al. (2010), consider a Lagrangian magnetic element in the middle of this ‘magnetic bubble’, which is threaded by (predominantly) horizontal magnetic field. Assume that the field is uniform and aligned in the x-direction, so that \({\bf{B = }}B\hat x\). For simplicity, assume that the upflow is axisymmetric with respect to the vertical direction \(\hat z\). Let the gradients of the local velocity field be

$${{\partial {v_x}} \over {\partial x}} = {{\partial {v_y}} \over {\partial y}} = \alpha,$$
(18)
$${{\partial {v_z}} \over {\partial z}} = \epsilon\alpha,$$
(19)

where α and ϵα can be considered horizontal and vertical expansions rates, respectively. ϵ is a measure of the anisotropic of the expansion rates in the vertical and horizontal directions. For ϵ = 1, expansion in the vertical and horizontal directions occur at equal rates. For ϵ ≈ 0, the expansion of the fluid element occurs predominantly in the horizontal directions. Substituting α and ϵα for the velocity gradients in Eqs. (1) and (11), respectively, one obtains

$${{D\ln \varrho } \over {Dt}} = - (2 + \epsilon)\alpha,$$
(20)
$${{D\ln B} \over {Dt}} = - (1 + \epsilon)\alpha.$$
(21)

This results in the scaling relation between field strength and mass density as

$$B \propto {\varrho ^\kappa },\,{\rm{where}}\,\kappa = {{1 + \epsilon } \over {2 + \epsilon }}.$$
(22)

When the expansion is isotropic (ϵ = 1), the field strength scales as Bϱ2/3. When the expansion is predominantly horizontal (ϵ = 0), the scaling relation becomes Bϱ1/2.

In the context of magnetic flux emergence, we are interested in how magnetic fields evolve on their journey from the solar interior to reach the surface. The evolution of a flux bundle in the deep solar interior with radius rHp with axial perturbations on scales λ ≫ Hp (i.e., a thin flux tube, cf. Spruit, 1981b) should be adequately described by Bϱ. Even in the case when the rHp, this scaling relation may hold since the slow variation in the axial direction means that the effect of flow gradients along the field (which leads to field intensification by stretching) is negligible. Pinto and Brun (2013) carried out global, 3D anelastic simulations of the rise of a toroidal flux tube embedded in a turbulent convection zone with a sustained dynamo field. The flux tube was introduced such that it has uniform buoyancy along its length. As a result, the entire tube tends to rise to the surface (except for interactions with convective flows). They found that the scaling relation from such a setup followed Bϱα, where α ≲ 1. In the AR region formation simulation of Cheung et al. (2010), the flattening of the rising loop as it approaches the top of the convection was shown to be more suitably described by Bϱ1/2. A localized Ω-loop that begins its journey from the deep solar interior and reaches the surface will likely experience a transition between the two scaling relations (i.e., α transitioning from 1 to 1/2).

3.2.2 Stability of the stratification to convective perturbations

The density and pressure contrast alone do not determine the (in)ability of material to emerge into the solar atmosphere. The stability of the stratification with respect to perturbations is more important. Consider a hydrostatic layer with pressure p(z), density ϱ(z), and temperature T(z) as functions of height. What would happen if a fluid element at a certain height in this layer were displaced by height dz? In the present discussion, we restrict our attention to purely hydrodynamic (i.e., non-magnetic) perturbations.Footnote 3 Consider the case where the fluid displacement is adiabatic (i.e., the fluid maintains its specific entropy s). As a result of the displacement to higher (or lower) ambient pressure, the fluid element would compress (or expand), resulting in a change in the density given by

$$d{\varrho _{{\rm{ad}}}} = {\varrho \over p}{\gamma _1}{{dp} \over {dz}}dz,$$
(23)

where γ1 was introduced as Chandrasekhar’s first adiabatic constant in Eq. (15). For the same change in height dz, the change in ambient density of the surrounding stratification is

$$d\varrho = {{d\varrho } \over {dz}}dz.$$
(24)

If \({{\partial {\varrho _{{\rm{ad}}}}} \over {dz}}\; < \;{{\partial \varrho } \over {dz}}\), a fluid element adiabatically displaced downward would have higher density than its new surroundings and if it were displaced upward, it would be less dense. In both cases, the buoyancy force will accelerate the fluid element to further the displacement from its original position. In this scenario, the stratification is considered to be convectively unstable to adiabatic perturbations. If \({{\partial {\varrho _{{\rm{ad}}}}} \over {dz}}\; < \;{{\partial \varrho } \over {dz}}\), the buoyancy force would be a restoring force that accelerates the displaced parcel back toward its original position. In this case the stratification is convectively stable. When \({{\partial {\varrho _{{\rm{ad}}}}} \over {dz}}\; < \;{{\partial \varrho } \over {dz}}\), the stratification is neutral to adiabatic perturbations and is marginally stable. This implies the buoyancy force does not diminish nor amplify vertical perturbations.

The classifications of the convective (in)stability of a stratification can be phrased in terms of logarithmic temperature gradients. The Schwarzschild criterion states that a stratification is convectively unstable if ∇ > ∇ad, where

$${\nabla _{{\rm{ad}}}} = {\left. {{{d\ln T} \over {d\ln p}}} \right\vert_s},$$
(25)
$$\nabla = {{d\ln T} \over {d\ln p}}.$$
(26)

The derivative in ∇ad can be calculated given an EOS (e.g., T = T(p, s)) and does not depend on the height profiles of T and p for the stratified layer of interest. The derivative in ∇ refers specifically to the height distributions of pressure and temperature in the chosen stratification. The cases ∇ = ∇ad and ∇ < ∇ad indicate marginal stability and stability of the stratification, respectively. For a detailed derivation of how one may compute ∇ad for general equations of state, we refer the reader to Chapter 4 of Kippenhahn and Weigert (1994). Chapter 6 of the same book also provides a detailed description of convective stability criteria for stellar interiors.

Figure 11 shows a plot of ∇, ∇ad, and the superadiabaticity δ = ∇ − ∇ad as functions of height in a 3D radiative hydrodynamic simulation of the near-surface layers of the convection zone and photosphere.Footnote 4 The plots shows how δ is positive in the convection zone (i.e., it is superadiabatic) and negative in the photosphere. So when a fluid element emerges from the convection zone into the photosphere, it is entering a stably stratified region that is difficult to penetrate. Within a distance of about one pressure scale height, the fluid element would become so dense relative to its surroundings that the (anti-)buoyancy force it to overturn and return to the convection zone. For this reason, a robust feature of many numerical simulations of flux emergence from the convection zone into the solar atmosphere show (at least a temporary) halt of the rising structure just beneath the photosphere. While subadiabaticity acts to suppress the further passage of magnetic flux into the atmosphere, other agents act to overcome this barrier. They include magnetic buoyancy instabilities (Section 3.3), the action of turbulent convective flows (Section 3.4).

Figure 11:
figure 11

Horizontally-averaged logarithmic temperature gradients in a radiative hydrodynamic simulation of the near-surface layers of the convection zone and overlying photosphere. The solid line shows ∇ad, the dotted line shows ∇, and the dashed line shows δ = ∇ − ∇ad. The photosphere has δ < 0 and is convectively stable. Image reproduced with permission from Cheung et al. (2007b), copyright by ESO.

3.2.3 Detection and inference of the subsurface structure of magnetic structures before and during emergence

From the results described above, one would expect horizontal flows in the near-surface layers of the convection zone to be much more pronounced than corresponding flows in the deeper layers as a buoyant flux rope rises toward the surface. Birch et al. (2013) applied helioseismic holography to two samples of Doppler flow measurements from SoHO/MDI. One sample contains pre-emergence ARs while the other sample contains observations without emergence. Based on results from their samples (more than 100 members per sample), they rule out the existence of large-scale coherent flows exceeding 15 ms−1 in the top 20 Mm of the convection zone over the course of one day before the arrival of emerging flux at the surface.

In contrast, photospheric horizontal flows coherent over AR-scales of up to 1–2 km s−1 have been reported (Toriumi et al., 2012) during episodes of emerging flux. The flow speeds from the AR formation simulations of Cheung et al. (2010) are consistent with the latter result. However, they are not directly comparable to helioseismic work of Birch et al. (2013) since the simulations only cover the top 7.5 Mm of the convection zone and do not take more than one day to traverse that depth to reach the photosphere.

Recently, Ilonidis et al. (2011) reported the helioseismic detection of emerging flux regions at depths of ∼ 65 Mm prior to their appearance at the photosphere. Their reported acoustic time anomalies of 12 to 16 s are one to two orders of magnitude larger than the travel time perturbations of pre-emerging flux regions reported by Leka et al. (2013). It is still an open question how the two results can be reconciled. Perhaps the effective acoustic travel time shifts from the two helioseismic schemes have different physical origin. In any case, flux emergence models that capture the photospheric layers and the top tens of Mms of the convection zone (e.g., the model of Stein and Nordlund, 2012, captures ‘only’ the top 20 Mm of the convection zone) are sorely needed to answer such questions. These models are needed for synthetic observables that can be used to to study what the helioseismic signals means in terms of subsurface conditions.

3.2.4 The extreme density contrast between the solar photosphere and corona

The previous discussion focused on the convection zone and photosphere. In fact the density contrast between the photospheric base and the corona can be even more drastic. Figure 12 shows a plane-parallel stratification used by Archontis et al. (2004) for the background atmosphere of his numerical simulations of flux emergence from the convection zone into the corona. The model stratification consists of a series of plane-parallel layers stacked on top of each other in mechanical equilibrium. In this model, the mass density drops by a factor of 108 from the photosphere to the corona.

Figure 12:
figure 12

A typical example of a model stratification used for numerical simulations of flux emergence from the convection zone into the corona. In this figure the height scale z is expressed in units of 170 km. The solid line (labeled P0) indicates the height profile of gas pressure in the initial plane-parallel stratification in units of p0 = 1.4 × 105 dyn cm−2. The dashed line indicates the mass density profile in units of ϱ0 =3 × 10−7 g cm−3. The dash-dotted line indicates the temperature profile in units of T0 = 5600 K. The solid line labeled Pm indicates the magnetic pressure (in units of p0) associated with a vertical cut through the mid-plane of the initially horizontal tube. Image reproduced with permission from Archontis et al. (2004), copyright by AAS.

It is not immediately obvious how one can intuitively ‘comprehend’ such a large density drop so an illustrative example is in order. Consider the case of coronal mass ejections (CMEs). The largest CMEs are global-scale events that occupy significant fractions of the coronal volume of the Sun (see Figure 13). It is estimated that the largest CMEs have mass contents of ∼ 1016 g (see Webb and Howard, 2012; Chen, 2011, for reviews of CME modeling and observations). What volume of material at photospheric densities would give an equivalent mass content? For a photospheric density of 10−7 g cm−3, it suffices to consider the mass contained with a box of 1 Mm2 extent in the horizontal directions and Hp ∼ 100 km extent in the vertical direction to enclose the same mass content. In other words, a single granule contains as much mass as the largest CMEs.Footnote 5

Figure 13:
figure 13

The solar atmosphere is strongly stratified: the mass contained in a single granule (left, Hinode SOT image of a sunspot surrounded by granulation) is comparable to the mass content of the largest CMEs (right, composite of running difference images from the SOHO/LASCO C2 and SDO/AIA instruments. Both quantities are of order 1016 g. (Credit for right image: NASA CDAW Data Center.)

With the previous example in mind, consider the mass content associated with an entire emerging AR. If all the mass threaded by the emerging field lines were brought up into the corona, there would be sufficient mass for hundreds, if not thousands of CMEs. The absence of such extreme events suggests that most of the mass originally attached to the emerging field lines must somehow never reach the higher solar atmosphere. Yet, since the magnetic field does indeed reach the corona (at least in ARs), there must be physical mechanisms to remove mass from the emerging flux structure. The nonlinear evolution of magnetic buoyancy instabilities and interaction of emerging flux with convective flows lead to undulated magnetic field lines that are amenable to reconnection, which is an efficient way to remove the mass burden from emerging flux (see Sections 3.4.33.4.6).

3.3 Magnetic buoyancy instabilities

3.3.1 Interchange and undular modes and their stability criteria

The concept of magnetic buoyancy instabilities is central to the physics of magnetic flux emergence. As mentioned in Section 3.2.2, the strong subadiabaticity of the photosphere acts like a barrier for passage of flux from the solar interior to the chromosphere and beyond. Magnetic buoyancy instabilities has been shown to be an effective physical mechanism for the aiding the further rise of magnetic flux from the photospheric layers into the chromosphere and corona.

To begin, consider a horizontal magnetic flux sheet of finite vertical extent embedded in a stratified layer. The flux sheet is embedded so that it is initially in mechanical equilibrim with the non-magnetic layers above and below, i.e.,

$${d \over {dz}}\left({p + {{{B^2}} \over {8\pi }}} \right) = - g\varrho.$$
(27)

The magnetic pressure of the flux sheet acts to support the gas above, which is a source of free gravitational potential energy. Magnetic buoyancy instabilities are pathways permitted by MHD evolution to liberate this free potential energy.

Magnetic buoyancy instabilities have close correspondence with their purely hydrodynamic counterpart. However, the presence of the magnetic field (mediated by the Lorentz force) introduces features which are distinct to the magnetic cases. The magnetic field introduces a preferred direction in the plane, so that different modes of instability arise depending on the angle between B and the wavevector of the perturbation k (Kruskal and Schwarzschild, 1954; Newcomb, 1961; Acheson, 1979; Hughes and Proctor, 1988). Since we are mainly interested in how magnetic buoyancy instabilities launch magnetic fields into the solar atmosphere, the following discussion will ignore the effects of rotation.

Figure 14 illustrates the different modes of the magnetic buoyancy instability. The kB mode it is called the interchange mode. It is similar to usual hydrodynamic Rayleigh-Taylor instability or convective instability, with the field lines remaining unbent so the resulting structure is two-dimensional (Figure 14b). In a narrow sense the so-called “magnetic Rayleigh-Taylor instability” refers to the interchange mode. Assuming adiabatic perturbations (i.e., δp/p = γ1δϱ/ϱ), the instability criterion of the interchange mode is given by (Acheson, 1979)

$${d \over {dz}}\ln \left({{B \over \rho }} \right) < - {{C_2^2} \over g}{{{N^2}} \over {V_A^2}},$$
(28)

where z is height, VA = B(4πϱ)−1/2 is the Alfvén speed, N is the Brunt-Väisälä frequency, and \({C_s} = \sqrt {{\gamma _1}p/\varrho}\) is the adiabatic sound speed. The Brunt-Väisala frequency is defined by

$${N^2} = {g \over {{\gamma _1}}}{d \over {dz}}\ln (p{\varrho ^{ - {\gamma _1}}}).$$
(29)

From the above expression, it is easy to see that N2 is proportional to the entropy gradient of the layer. In the case of N = 0, i.e., neutral for convective instability, Eq. (28) reduces to

$${d \over {dz}}\ln \left({{B \over \rho }} \right) < 0.$$
(30)

When the layer is subadiabatically stratified (entropy increases with z), N2 > 0 and criterion (28) is more stringent and requires a steeper gradient in B/ϱ to overcome the stabilizing tendency of the stratification. The interchange mode mode has been considered as a mechanism of breaking up the large-scale troidal flux at the base of the convection zone into the thin, buoyant flux tubes that are eventually observed as sunspots at the surface (Parker, 1977; Cattaneo and Hughes, 1988; Cattaneo et al., 1990; Matthews et al., 1995; Wissink et al., 2000a).

Figure 14:
figure 14

Schematic illustration of the different modes of the magnetic buoyancy instability of a horizontal flux sheet. Image reproduced with permission from Matsumoto et al. (1993), copyright by AAS.

The kB mode is called the undular mode (Figure 14c). In this mode the field lines undulate and the plasma slides down along the magnetic fields from the crests to the troughs, thus the rising crests get lighter and the troughs get heavier, facilitating the amplification of the perturbation. The instability criterion of the undular mode (or the mixed mode, see Figure 14d and 14e) is given by (Acheson, 1979)

$${d \over {dz}}\ln B < - {{C_s^2} \over g}\left({k_\parallel ^2(1 + {{k_z^2} \over {k_ \bot ^2}}) + {{{N^2}} \over {V_A^2}}} \right),$$
(31)

where kz, k, and k are the wavenumbers in the vertical direction and in directions parallel and perpendicular to magnetic field, respectively. By comparing with Eq. (30), one can see that the instability criterion is broader in the undular mode. For instance, when N = 0,

(32)

Namely, there exist unstable modes when

$${{dB} \over {dz}} < 0.$$
(33)

The instability criterion Eq. (33) can be also obtained by the energy principle (Newcomb, 1961). For the undular mode to grow in an adiabatic stratification, the instability criterion is satisfied as long as the field strength decreases with height (as opposed to B/ϱ for the interchange mode). The undular mode can therefore be unstable even when the equilibrium is stable for the interchange mode.

The undular mode is also called Parker instability because it was first considered by Parker (1966) as the mechanism for interstellar cloud formation (see also Shu, 1974; Zweibel and Kulsrud, 1975; Matsumoto and Shibata, 1992; Kim and Hong, 1998). The same concept was later applied to flux emergence in the solar interior (Acheson, 1979; Spruit and van Ballegooijen, 1982). See Fan (2009b) for a review of flux emergence in the deep convection zone.

3.3.2 Linear growth rate and eigenfunctions

Which mode grows faster depends on the magnetic profile as well as the background stratification. However, in general, the interchange modes grow faster than the undular modes when the magnetic field is isolated and there is a discontinuity (or sharp gradient). If the direction of the magnetic field is constant, i.e., with no magnetic shear, the behaviour of the interchange mode is similar to the hydrodynamic Rayleigh-Taylor instability. Namely, the wavenumber of the largest growth rate is limited by diffusion and hence in the ideal limit, the growth rate is larger for larger wavenumber. In the presense of magnetic shear, the high wavenumber modes are suppressed (Cattaneo and Hughes, 1988; Cattaneo et al., 1990; Kusano et al., 1998; Nozawa, 2005).

In contrast to the interchange mode, the undular model exhibits a characteristic length-scale even in the absence of diffusion. At sufficiently small length-scales, the restoring effect of magnetic tension (stemming from bent field lines) dominates over buoyancy. The undular mode therefore has a critical wavenumber above which the instability is stabilized by magnetic tension. The growth rate of the undular mode has a maximum value at wavelength λ = 2π/k ∼ 10Hp. Pariat et al. (2004) noted that the typical separation between Ellerman bombs in an observed emerging flux region is comparable to this length-scale.

Figure 15 shows the linear growth rates of modes for the magnetic buoyancy instability of an isolated magnetic sheet embedded in a super-adiabatically stratified layer mimicking the solar convection zone (Nozawa, 2005). The left panel shows the case witout magnetic shear: kx is wavenumber in the direction of the magnetic field and ky is that in the transverse direction, and \(k = \sqrt {k_x^2 + k_y^2}\). The interchange mode (ky/kx = 1010 ∼ ∞) has growth rate that monotonically increases with k, while the pure undular mode (ky/kx = 0) and the mixed modes have the maximum growth rate at finite k. The right panel shows the case with magnetic shear (i.e., variation of background field orientation with height). In comparison to the case without magnetic shear, the growth rate of the ky/kx = 1010 mode is reduced, while that of the ky/kx = 0 mode is increased.

Figure 15:
figure 15

Growth rates of different modes for the magnetic buoyancy instability of an isolated flux sheet. Here kx and ky are the wavenumbers parallel and perpendicular to the magnetic field. The case ky/kx = 0 represents the pure undular mode while ky/kx → ∞ represents the interchange mode. The left panel shows the growth rates for the case without magnetic shear (in the vertical direction) while the right panel shows the corresponding growth rates with magnetic shear. Image reproduced with permission from Nozawa (2005), copyright by PASJ.

The undular mode has eigenfunctions with velocity amplitudes that are larger in the horizontal directions than in the vertical direction. This corresponds to the fact that the dominant motion of the undular mode is the plasma draining along field lines. Figure 16 shows the shape of the eigenfunction calculated by Horiuchi et al. (1988). This paper considers the case of a galaxy where the temperature is constant and the gravity varies as a function of the height from the galactic plane, but the basic characteristics are the same for the solar atmosphere. The eigenfunction and the growth rate of the undular mode in a super adiabatically stratified layer is also shown in Nozawa et al. (1992). It is also worth noting that the eigenfunction of the interchange mode tends to concentrate near the boundary layer of the magnetic flux where there is a sharp magnetic gradient (Cattaneo and Hughes, 1988).

Figure 16:
figure 16

The eigenfunction of the undular mode. Image reproduced with permission from Horiuchi et al. (1988), copyright by PASJ.

3.3.3 Nonlinear evolution of magnetic buoyancy instabilities in emerging flux

Shibata et al. (1989a,b) first applied the undular mode to explain the emerging flux in the solar atmosphere and performed 2D MHD simulations. Figure 17 shows the result of the 2D simulation by Shibata et al. (1989b). The initial condition is a horizontal magnetic flux sheet embedded in a stratified layer consisting of an isothermal photosphere/chromosphere and an overlying hot corona. A initial flux sheet is embedded in the photosphere/chromosphere, and a small perturbation is given in horizontal velocity (in order to be close to the eigenfunction, see the dicussion above) locally at the center of the flux sheet. The perturbation then grows and the crest of the undulating magnetic field rises into the upper atmopshere, forming so-called Ω-loops. Since the height of Ω-loops is larger than many scale heights, field-aligned downflows become supersonic and standing shocks are formed at loop footpoints. The rise velocity of the Ω-loops (10–15 km s−1) and the velocity of the field-aligned downflow (30–50 km s−1) are consistent with the observations (e.g., Chou and Zirin, 1988; Otsuji et al., 2007).

Figure 17:
figure 17

2D MHD simulation of undular mode instability of a flux sheet in the photosphere. Image reproduced with permission from Shibata et al. (1989b), copyright by AAS.

An extension of the work on the nonlinear evolution on the undular mode in the flux emergence context was carried out by Kaisig et al. (1990), Shibata et al. (1990a), and Nozawa et al. (1992). Because of the superadiabatic stratification, convective intensification of magnetic fields (Parker, 1978; Spruit, 1979; Spruit and Zweibel, 1979) was found to occur at the footpoints of the emerging Ω-loops, but otherwise the nonlinear evolution in the atmosphere looks similar to Shibata et al. (1989a). The reason is that, as the magnetic field rises through the photosphere and the chromosphere over many scale-heights, the gas pressure rapidly decreases, and hence above a certain height the magnetic pressure becomes dominant over the gas pressure. The ensuing evolution is basically the expansion into the free space driven by the magnetic pressure of emerging flux itself. Shibata et al. (1990b) found a group of self-similar solutions that describe the nonlinear evolution of the undular mode reasonably well.

The first 3D MHD simulations of solar flux emergence were carried out by Matsumoto and Shibata (1992) and Matsumoto et al. (1993). The initial condition is essentially the same as previous 2D simuation, but the extra degree of freedom available in 3D allows the interchange mode to also play a role. It was found that Ω-shaped emerging loops similar to 2D simulation were also formed, as shown in Figure 18. Thus, the overall evolution in the nonlinear stage was determined by the undular mode, though the presence of the interchange mode led to fine structure with variations perpendicular to the background field.

Figure 18:
figure 18

Magnetic field lines in the 3D MHD simulation of emerging flux. Image reproduced with permission from Matsumoto et al. (1993), copyright by AAS.

The effect of magnetic shear on the nonlinear evolution of the magnetic buoyancy instability in the context of emerging flux was also examined in 2D (Kusano et al., 1998) and in 3D (Nozawa, 2005). In 2D, the general finding is that if the magnetic field at the upper interface between the magnetic sheet overlying plasma is mostly perpendicular to k (i.e., the component out of the simulation plane is dominant), the interchange mode grows at the interface and only small scale structure is produced. On the other hand, if the magnetic field at the interface is mostly parallel to k (i.e., dominated by horizontal component in the simulation plane) the resulting structure is similar to the pure undular mode. In 3D, the interchange model initially dominates since it has larger growth rate at small-scales (see Figure 15). However, in the non-linear phase the undular mode dominates the overall evolution. This is consistent with the result without magnetic shear.

Similar nonlinear simulations of magnetic buoyancy instabilities in the context of the formation of flux ropes at the base of the convection zone have been extensively carried out (e.g., Cattaneo and Hughes, 1988; Matthews et al., 1995; Wissink et al., 2000a; Fan, 2001a; Kersalé et al., 2007; Favier et al., 2012).

3.3.4 Suppression of horizontal expansion

As discussed in Section 3.2, an isolated magnetic flux tube tends to preferentially expand in the horizontal direction upon reaching the sub-adiabatically stratified photosphere (see Figures 7 and 8). The horizontal expansion weakens the magnetic field and, hence, causes the loss of coherence and prevents the further rise. This effect was already found in the 3D simulation by Matsumoto et al. (1993). Moreover, Matsumoto et al. (1993) performed a simulation of the rise of isolated flux tube with periodic boundary conditions to mimic the rise of multiple flux tubes. Figure 19 shows a snapshot in the simulation during the nonlinear phase. It was found that, even though the flux tubes are isolated in the photosphere, they rise and expand horizontally and eventually collide and merge each other. Thus, the horizontal expansion is suppressed and the magnetic flux can continue to rise vertically. Krall et al. (1998) showed by 2D simulation that the presense of the ambient magnetic field in the atomosphere can also suppress the horizontal expansion and, thus, helps the rise of the magnetic flux in the atmosphere.

Figure 19:
figure 19

Surface of constant magnetic field strength from the result of 3D simulation of the rise of isolated flux tubes. Image reproduced with permission from Matsumoto et al. (1993), copyright by AAS.

Another way to suppress the horizontal expansion is the magnetic tension by the azimuthal component of a twisted flux tube. As discussed in Section 3.6 certain amount of twist is also required to keep the coherence of a flux tube as it rises through the convection zone. Therefore, many authors carried out simulations of the buoyant rise of a twisted flux tube embedded in the convection zone in the way described in Section 3.1.1.

Figure 20 shows one typical example by Fan (2001b). Initially a twisted flux tube is embedded just below the photosphere. The density is decreased at the center so that only a part of the flux tube becomes buoyant, leading to the formation of an Ω-shaped flux tube. Note that the individual field lines do not necessarily have a simple Ω-shape except for those near the central axis of the tube: see panels (a) and (b) of Figure 20. As the flux tube expands the azimuthal component becomes more dominant. So the outer most field lines that emerged at fisrt are oriented almost perpendicular to the axis of the tube, as seen in panel (f) of Figure 20. It can be also seen in the distribution of the vertical magnetic field at the photosphere. These features are commonly observed in many simulation studies (e.g., Magara and Longcope, 2003; Manchester IV et al., 2004; Archontis et al., 2004; Toriumi et al., 2011) though there are significant and interesting differences in the emergence dynamics that depends on the parameters, as discussed in the following sections.

Figure 20:
figure 20

Result of 3D simulation of the emergence of a twisted flux tube. Panels (a) and (b) show magnetic field lines at different times. In both panels, the red line indicates locus of the tube axis, which is seen to be trapped in the photosphere while the peripheral field lines are able to reach the corona. Image reproduced with permission from Fan (2001b), copyright by AAS.

When the bulk of the rising part of the flux tube emerged in the atmosphere and their magnetic pressure dominates the surrounding gas pressure, the tube starts to expand by its own magnetic pressure gradient similarly to the case of 2D undular mode discussed above. Magara (2004) presented a detailed analysis of the forces acting on the individual field lines and found that the dynamics of the periphery field lines is akin to the free expansion of outer field lines in 2D undular mode simulations (Shibata et al., 1989a,b). In contrast, inner field lines are subject to strong confinement by adjacent twisted field lines.

3.3.5 Instability in the photosphere and two-step emergence

The works reviewed in the previous sections consisted of initial conditions in which horizontal flux sheets were embedded into the near-surface layers (or even in the photospheric) layers. This type of setup facilitated the detailed study of magnetic buoyancy instabilities but does not address how the magnetic flux reached those layers. The discussion in this section focuses addresses this connection.

Figure 21 shows the result of a 2.5D simulation of the emergence of a twisted flux tube by Magara (2001). In this case, the axis of the tube is perpendicular to the simulation plane. Due to invariance in this third direction, there is no undulation along the tube axis. The tube is initially located at a depth of 1 Mm in the convection zone and rises due to an imposed density deficit. As the tube impinges upon the interface between the convection zone and photosphere, it expands preferentially in the horizontal direction (see also Section 3.2) so that the outermost field lines (projected onto the plane) become nearly horizontal (e.g., see panel (c) of figure). This is similar to the setup for magnetic buoyancy instabilities and the mode with the faster growth rate can take over. In this instance the operation of the Parker instability (undular mode) was identified as the mechanism for the further passage of flux into the corona. Strictly speaking, however, the presence of the longitudinal component of the field (i.e., out-of-plane component) means the instability becomes a mixed mode in 3D. In any case, the primary physical mechanism for the formation of Ω-loops is the draining of plasma along field lines, which accentuates the buoyancy of the crest of the loop.

Figure 21:
figure 21

2.5D simulation of the cross section of an emerging twisted flux tube by Magara (2001). In all panels the color coding indicates the density stratification and the arrows indicate the flow velocity in the plane of the simulation domain. The development of a magnetic buoyancy instability at the top of the tube (see panels d to f) leads to the dramatic expansion of the emerging flux into the coronal layers. Image reproduced with permission, copyright by AAS.

In models where a rising flux tube impinges into the atmosphere, the horizontally expanded fields at the base of the photosphere never reach a true state of mechanical equilibirum. Instead, they continue to expand horizontally while continuously gaining flux due to pile-up from below. As reported by Archontis et al. (2004) from their 3D emergence simulations of individual twisted flux tubes, the pile-up of flux can increase the field strength to the point where the criterion for the magnetic buoyancy instability is satisfied. They reported that at the moment when the instability criterion was satisfied, Ω-loops began to be launched into the coronal layers.

The arrival of magnetic flux below the photosphere, the suppression of its emergence, and its subsequent rapid rise into corona due to magnetic buoyancy instabilities is a salient feature of a number of numerical models of flux emergence (see also Murray et al., 2006; Hood et al., 2012) and is termed two-step emergence (Toriumi et al., 2011; Toriumi and Yokoyama, 2011). Figure 22 shows another example. Unlike many models which begin with horizontal flux tubes, this model has an initial condition consisting of a semi-toroidal, twisted flux tube embedded in the convection zone (Hood et al., 2012). The drainage of plasma along the legs of the flux tube enhances the buoyancy at the crest, which rises toward the photosphere and flattens like a pancake. Despite the different initial condition (cf. Archontis et al., 2004), when the pancake becomes strong enough so that the instability criteria (Eq. 31) is satisfied, the buoyancy instability starts to grow. In the case shown in Figure 22 two loops (for flux bundles) grow because of the mixture of undular and interchange modes.

Figure 22:
figure 22

Isosurface of the magnitude of the magnetic field in the 3D simulation of the emergence of twisted and toloidal flux tube. The times are t = 0 (top left), 40 (top right), 60 (bottom left), and 80 (bottom right). Image reproduced with permission from Hood et al. (2012), copyright by Springer.

Whether and when the instability in the photosphere occurs depends on the properties of the underlying structure. In the case of flux tubes, the degree of magnetic twist in the tube plays an important role in determining whether the instability may occur.Footnote 6 A common choice of magnetic profiles for twisted flux tubes is such that the longitudinal (Bl) and azimuthal field (Bθ) components are related by Bθ(r) = qrBl(r), where r is the distance from the tube axis. In this case q is a measure of the twist inherent in the tube (see Eqs. 40 and 41).

Figure 23 shows magnetic field lines of emerging flux tubes in simulations with different values of q (Murray et al., 2006). In the strongly twisted case (q = 0.3, top), the horizontal expansion is suppressed by the magnetic tension of the azimuthal field, so that the magnetic dome smoothly expands into the upper atomosphere. In the weakly twisted case (q = 0.1, bottom), the horizontal expansion dominates and the magnetic field fails to emerge into the upper atomosphere (at least within the time of the simulation). In the intermidiate case, the buoyancy instability starts to grow and the loops are rising at two location in the magnetic pancake. The parameter study of 2.5D simulations by Toriumi et al. (2011) adds support to this point. Figure 24 (reproduced from their paper) plots the height of the tube apex as a function of time for simulations with different values of q. In the strongly twisted regime (q > 0.25), the deceleration of the emerging flux in the photosphere is not significant, while in the weak twist regime (q ≤ 0.05) the magnetic flux does not emerge into the upper atomosphere. In the intermediate regime (0.1 ≥ q ≥ 0.2) the magnetic flux stays is temporarily trapped in the photosphere, and after some delay is seen to rapidly launch into the coronal.

Figure 23:
figure 23

Traces of magnetic field lines in three separate numerical experiments of twisted flux tube emergence. Panels (a)–(c) show simulations with an initial twist parameters of q = 0.3, q = 0.2, and q = 0.1, respectively. The radius initial radius of the tube in these experiments is a = 2.5 (see Eqs. (40) and (41)). Image reproduced with permission from Murray et al. (2006), copyright by ESO.

Figure 24:
figure 24

From the parameter study of twisted flux tube emergence, a plot showing the height of the tube apex as functions of time and magnetic twist (q). With sufficiently weak twist (q = 0.05), the tube is trapped in the photosphere. For twist q = 0.1 and higher, the magnetic field is launched into the corona due to magnetic buoyancy instabilities. The time delay between arrival at photosphere and passage into the corona decreases with increasing twist. Image reproduced with permission from Toriumi et al. (2011), copyright by PASJ.

3.3.6 Rayleigh-Taylor instability at the tops of emerged loops

Isobe et al. (2005a) found in the 3D simulation of the undular mode of a magnetic sheet that the top of the emerged loops become denser than the inner part, forming a top-heavy configuration that is unstable to the Rayleigh-Taylor instability. Figure 25 shows the temporal evolution of the field lines and density distributions. One can see that as the Ω-loops rise the outer most part becames denser than the inner part.

Figure 25:
figure 25

Magnetic Rayleigh-Taylor instability. The color shows mass density distribution in panels (a)–(c) and current density distribution in panel (d). Image reproduced with permission from Isobe et al. (2007a).

In this Lorentz-force dominated regime, only the interchange modes grow, generating filamentary structures aligned along the magnetic field (similar to arch filament systems, see Bruzek, 1967). In addition, the sinking spikes and rising plumes produce magnetic shear and hence electric current between them, as shown in panel (d) of Figure 25. Isobe et al. (2005a) conjectured that dissipation of the small scale currents by magnetic reconnection may contribute to the coronal heating in the emerging flux structure.

The mechanism of the formation of the top-heavy structure in the emerging loops was further investigated by Isobe et al. (2006). Two effects contribute to the top-heavy structure. The first is the compression of the outermost part of the emerging flux between the corona above and the rising inner loop below. The second, which was more efficient in the particular case they analyzed, is due to the deviation from self-similar evolution (Shibata et al., 1990b) so that the outermost field lines have larger radius of curvature (i.e., more horizontal) and hence the plasma draining along these field lines is less effective.

Although it has not attracted much attention in the literature, the top-heavy configuration is also commonly seen in numerical simulations of the emergence of twisted flux tubes, e.g., Figure 3 of Toriumi et al. (2011). However, it was pointed out by Arber et al. (2007) that if the effect of partial ionization in the chromosphere were taken into account, the neutral drift reduces the amount of lifted mass by the emerging flux and the top-heavy configuration does not form. The effect of partial ionization will be discussed in detail in Section 3.8.

The operation of the magnetic Rayleigh-Taylor instability in the upper solar atmosphere is also important for understanding prominence dynamics. In a series of papers, Hillier et al. (2011, 2012a,b) performed numerical MHD simulations of the nonlinear evolution of the instability in idealized models of quiescent prominences. In their simulations, the background atmosphere and magnetic field follow the Kippenhahn-Schlüter equilibrium (Kippenhahn and Schlüter, 1957; Priest, 1982). This magnetostatic model consists of upward concave field lines with density enhancements at the dips. The initial background field is invariant in the vertical direction. To trigger the Rayleigh-Taylor instability, a buoyant flux tube (in terms of localized density deficit) is inserted into background model at some height. The top-heavy contact discontinuity between the buoyant tube the overlying background state is then given random velocity perturbations with an amplitude of order 0.01cs, where cs is the local sound speed. The nonlinear evolution of the system leads to buoyantly rising cavities and dense downflows with morphologies that are qualitatively similar to those discovered in Hinode/SOT observations of quiescent prominences (Berger et al., 2008, 2010). This finding suggests that on the Sun, emerging flux rising into quiescent prominences maybe a important driver of the local dynamics.

3.4 Magnetoconvection

Outside of sunspots, the Sun’s photosphere is pervaded by an ever evolving tessellation of granules (see left panel of Figure 13). Individual granules on the solar surface have horizontal length-scales of ∼ 1 Mm and lifetimes ∼ 10–20 min. They are the most conspicuous manifestation of convection and is evidence for convective energy transport as the dominant mechanism for the outward transport of energy from the solar interior into the solar atmosphere. By treating radiative heating and cooling in a consistent and realistic manner, radiative convection simulations of solar granulation have been remarkably successful at reproducing observed/observable properties of solar granulation (Nordlund, 1985; Steffen et al., 1989; Stein and Nordlund, 1989, 1998; Wedemeyer et al., 2004; Leenaarts and Wedemeyer-Böhm, 2005; Cheung et al., 2007b; Beeck et al., 2012). In the remainder of this section we will briefly describe the underlying physics of solar surface convection and focus on their relevance for magnetic flux emergence. For details concerning observations and models of solar granulation, we encourage the reader to consult the Living Reviews article by Nordlund et al. (2009).

The dominant mechanism of the vertical transport of energy in the solar convection zone is by means of the enthalpy flux Fenthalpy = vz(ϱε + p). A net positive-Fenthalpy results from correlations between upflows and enhancements in internal energy. In other words, upflows in convective cells have relatively higher temperature and internal energy. At the solar surface, this property is responsible for granules (upflows) being relatively bright and hot compared to downflow lanes in between granules. The transition from the convection zone to the photosphere occurs in layer of a few tens of km where the opacity of the gas drops precipitously (Stein and Nordlund, 1998). Over this distance, the gas goes from being optically thick to optically thin and radiation takes over as the dominant mechanism for vertical energy transport. The radiative cooling experienced by the surface plasma leads to compression of the gas, which creates an overdense layer near optical depth unity with a diminished specific entropy. This gradient in the specific entropy drives the development of the convective instability and the creation of downflows.

The following sections highlight a few physical concepts concerning the near-surface layers of the convection zone that are particularly pertinent to the topic of magnetic flux emergence.

3.4.1 Treatment of the convection zone in various models of flux emergence

Different groups in the flux emergence modeling community have taken various approaches to the treatment of the convection zone in their models. The choice of how the convection zone is treated depends on the science question being addressed. For some questions, it suffices to treat the convection zone as a plane-parallel stratified layer. For others, 3D radiative transfer calculations are need to yield realistic cooling rates that drive convective flows.

Beginning from the very first numerical simulations of magnetic flux emergence from the convection zone to the atmosphere (Shibata et al., 1989a), there has been a steadily growing body of work that model the rise of buoyant magnetic fields through a model convection zone that is adiabatically (or superadiabatically) stratified but without convective flows. By intention, many these studies set out to study the magnetic flux emergence process in idealized setups in order to isolate or to highlight how other physical mechanisms drive the emergence process. For instance, the use of a series of plane-parallel stratified layers (e.g., see Figure 12) allows study of the initiation and growth of magnetic buoyancy instabilities (see Section 3.3). Other studies systematically examine the dynamics that occur when emerging magnetic flux reconnect with pre-existing coronal field (Galsgaard et al., 2005) and for this purpose, it is desirable to have simple setups that allow one to control the angle of alignment between the two flux systems. These types of parameter studies are much harder to control when convective flows are present in the system.

To further simplify the setup, some models ignore the convection zone entirely by choosing computational domains that are entirely above it. For instance, the numerical simulations such as those by Forbes and Priest (1984) and Fan and Gibson (2003, 2004) set the lower boundary of the domain to be at the base of the photosphere. The emergence of a flux tube into model corona was driven by imposing a time-dependent −v × B electric field consistent with the vertical passage of a twisted semi-torus through a horizontal plane. Such a setup allowed the authors to precisely control the amount of twist they wished to inject into the model corona associated with the emerging torus in order to study the onset and growth of the non-linear kink instability in the low plasma-β regime.

Some numerical simulations that address the influence of convective flows on emerging flux do so without resorting to performing full 3D radiative transfer calculations to determine the radiative cooling term in near the photosphere. There exists different simplifying approaches. In one approach, one begins with a stack of idealized hydrostatic stratifications on top on one another (similar to Figure 12), making sure that the layer representing the convection zone is superadiabatically stratified (δ > 0). The convective instability would then amplify randomly imposed perturbations into convective plumes.Footnote 7 With an appropriately chosen bottom boundary condition for energy input and a specified cooling function near z = 0, the system would evolve to a dynamic convective pattern that is statistically stationary (i.e., in terms of horizontal energy fluxes, RMS speed of flows, etc.). This is the approach taken by Isobe et al. (2008), who examined how a uniform magnetic field embedded in a computational domain with fully-developed convection to study how magnetoconvection can lead to flux emergence. We will continue discussing this work in Section 3.4.5.

In parallel with developments in numerical simulations of emerging flux, realistic modeling of photospheric granular convection (Nordlund, 1985; Steffen et al., 1989; Stein and Nordlund, 1989, 1998; Wedemeyer et al., 2004; Leenaarts and Wedemeyer-Böhm, 2005; Cheung et al., 2007b; Beeck et al., 2012) and magnetoconvection (Vögler et al., 2005; Khomenko et al., 2005; Schüssler and Vögler, 2006; Stein and Nordlund, 2006; Vögler and Schüssler, 2007; Cameron et al., 2007; Schüssler and Vögler, 2008; Steiner et al., 2008; Kitiashvili et al., 2010a,b; Leenaarts et al., 2010; Sainz Dalda et al., 2012) have matured to a state where simulations could reproduce a wide range of observable from high-resolution observations of the photosphere. In order to produce granular convection simulations with granular flow speeds, brightness contrast, and thermodynamic structure with sufficiently high fidelity to observations, their simulations must use an EOS that includes partial ionization of elements in the gas mixture at solar abundances. This is essential since the latent heat released in recombining atoms (mainly hydrogen and helium) has a direct impact on important thermodynamic variables such as the adiabatic index γ1 and the adiabatic temperature gradient ∇ad. A simple way to intuitively understand the effects of latent heat is to consider a fluid element that rising adiabatically in an atmosphere. As it rises, the fluid element expands and the temperature of the element decreases. If ions are allowed to recombine with free electrons, the latent heat released will suppress the temperature decrease. This effect was studied in detailed (Rast and Toomre, 1993), who found that ionization enhances buoyancy driving and makes convection more vigorous. This picture provides a simple reason why partial ionization changes ∇ad and thus the stability of a stratification. Usually, changes in ionization degree is included in the EOS in order to realistically capture the partitioning of the net energy flux through the system between radiative, enthalpy, and kinetic fluxes. Of course, the resulting mean stratification has δ > 0. However, even with an equation of state that includes ionization changes, one may choose to construct a adiabatic background stratification. This approach was taken by Leake and Linton (2013), who wanted to focus their study on other effects (in this case, Pedersen current dissipation) and neglect convective flows.

3.4.2 Asymmetry between upflows and downflows

Convection in a stratified layer with multiple pressure scale heights leads to an asymmetry between upflows and downflows (Hurlburt et al., 1984). Whereas material in upflows expand, material in downflows compress. Together with mass conservation, this leads to broad, gentle (i.e., lower speeds) upwellings separated by narrower and faster downflow lanes. When a magnetic flux structure (say, a flux tube) is rising through the convection zone, different portions of it will encounter downflows. The downflows exert an effective aerodynamic drag that tends to pull portions of the flux tube downward. Figure 26 shows results from a numerical simulation of buoyant flux tubes interacting with flows in a convective layer (Fan et al., 2003). In this particular simulation the flux tube is insufficiently buoyant so that some segments of the tube that encounter convective flows flows are pinned down to the bottom of the computational domain. Given sufficient magnetic buoyancy (i.e., higher field strengths), the authors found that is was possible for the flux tube to rise effectively unimpeded by convective downflows. Abbett et al. (2004) extended this study by performing the simulations over multiple turnover times to examine the long term behavior of the magnetic field in a convecting layer (more on this later).

Figure 26:
figure 26

Vertical distribution of the buoyancy force (i.e., relative density perturbation) in a simulation of buoyant flux tubes interacting with ambient convective flows. Dark regions show antibuoyant downflows. The white horizontal structure in the top panel shows the initially horizontal flux tube. In a later snapshot (lower panel), two portions of the tube are seen to be pinned down by downflows. Image reproduced with permission from Fan et al. (2003), copyright by AAS.

The simulations performed by Fan et al. (2003) were carried out in the anelastic approximation, which is able to capture high density contrasts in stratified layers but assumes that flow speeds are negligible compared to the speed of sound. While this is suitable for deep layers in the convection zone, it is not suitable for the near-surface layers where downflows can attain speeds with Mach numbers M ∼ 0.1–1. However, the general conclusions drawn from the work of Fan et al. (2003) and Abbett et al. (2004) still hold in compressible MHD (e.g., see Bushby and Archontis, 2012, for idealized simulations of flux tube interacting with compressible convection). As an extension of this type of study, Cheung et al. (2007a) performed radiative MHD simulations of the rise of twisted flux tubes in granular convection. They found that, due to the high Mach numbers found in convective downflows, it is virtually impossible for buoyant flux tubes to rise unimpeded by the convection. As a result of the very vigorous convection near the solar surface, the morphology of emerging flux regions (EFRs) will be highly structured by the surface granulation pattern. On the other hand, emerging magnetic structures with sufficient field strengths can influence the granulation pattern in a way that temporarily yields anomalously elongated granules and dark lanes that are coincident with upflows (see Figure 27 for an example). Subsequent numerical simulations with similar setups have since confirmed the robustness of these results (see, e.g., Tortosa-Andreu and Moreno-Insertis, 2009; Martínez-Sykora et al., 2008). Anomalous, elongated granules and darkenings have also been reported in observations of regions of emerging flux (Zwaan, 1985; Strous and Zwaan, 1999; Cheung et al., 2008; Schlichenmaier et al., 2010; Lim et al., 2011).

Figure 27:
figure 27

Continuum intensity images (top, at 5000 Å) and (bottom) synthetic magnetograms of a simulated emerging flux region. The snapshot at t = 66 min shows the appearance of anomalously large granules and intergranular bright points in the emerging flux region. The magnetogram shows mixed polarity structures down to the granulation scale. Image reproduced with permission from Cheung et al. (2008), copyright by AAS.

3.4.3 Undulation of emerging magnetic field lines by granular flows

Consider a magnetic structure with predominantly horizontal field initially embedded immediately below the photosphere. Its interaction with the convective flows would introduce undulations in the magnetic field lines such that the field lines have crests coincident with upflows and troughs coincident with downflows. Figure 28 shows such a scenario from the flux emergence simulations of Tortosa-Andreu and Moreno-Insertis (2009). The time sequence of magnetic field lines in the neighborhood of a granule shows how the initially (mostly) horizontal field is undulated to form Ω-loops that poke through the centers of granules as well as U-loops that submerge at intergranular lanes. The formation of such U-loops is a robust feature of MHD simulations of flux emergence through convective flows (Tortosa-Andreu and Moreno-Insertis, 2009; Cheung et al., 2010; Fang et al., 2010) and of simulations of near-surface magnetoconvection (Stein and Nordlund, 2006; Abbett, 2007).

Figure 28:
figure 28

3D visualization of the interaction of emerging horizontal magnetic field with granular convection (from Tortosa-Andreu and Moreno-Insertis, 2009). The grey scale maps show the distribution of vertical velocity at constant geometrical height in the photosphere. This sequence of magnetic field line visualizations shows how the convective flows undulate emerging field lines at the scale of granulation and leads to the formation of submerging U-loops aligned with downflows. Image reproduced with permission, copyright by ESO.

In the last (bottom) snapshot of Figure 28, the set of U-loops along an intergranular lane are beginning to be pinched due to colliding outflows from the adjacent granules. This brings together oppositely directed flux. Magnetic reconnection between the opposite polarities results in the creation of twisted magnetic structures (or plasmoids) with lengths comparable to the intergranular lane (∼ 1 Mm). If there reconnecting field had an initial horizontal component aligned along the intergranular lane, the plasmoid would have a guide field along its axis and structure would be akin to a granular-scale twisted magnetic flux rope. Figure 29 (also taken from Tortosa-Andreu and Moreno-Insertis, 2009) illustrates an example of such as structure forming as a direct consequence of convective flows acting on the emerging magnetic field.

Figure 29:
figure 29

3D visualization of the creation of small-scale twisted magnetic flux tube at an intergranular lane immediately below the photosphere (from Tortosa-Andreu and Moreno-Insertis, 2009). Images with yellow color-coding show the distribution of vertical velocity at constant geometrical height in the photosphere. Image reproduced with permission, copyright by ESO.

3.4.4 Mixed polarity patterns and serpentine field lines

An observable manifestation of magnetic field undulation by granular flows (discussed in the previous section) is the pattern of mixed polarities fields in photospheric magnetograms of emerging flux regions. Figure 31 shows an observation of an emerging flux region (which eventually becomes NOAA AR 10978) taken with the Spectropolarimeter instrument (Lites et al., 2001; Tsuneta et al., 2008) onboard the Hinode spacecraft (Kosugi et al., 2007). The existence of mixed polarity field down to the granulation scale is striking, but what is more important is that the orientation of mixed polarity pairs is not random. Rather, the mixed polarity fields are, on average, roughly oriented along the direction connecting the preceding and leading polarity spots (see also Otsuji et al., 2011; Centeno, 2012).

The emergence of serpentine field lines into the solar atmosphere have observational consequences beyond the structured mixed polarity patterns in photospheric magnetograms. Pariat et al. (2004) analyzed observations of an emerging flux region from the Flare Genesis Experiment and reported that locations of Ellerman bombs (EBs, transient brightenings in the wings of the Hα line) are correlated with the locations of Bald Patches (BPs), which they define as locations where Bz =0 and B · ∇BZ > 0. BPs are dips in magnetic field field lines. As illustrated in Figure 30, magnetic dips are found in between adjacent loops. In this picture, magnetic reconnection between adjacent loops along a serpentine field line lead to plasma heating in the chromosphere and the appearance of Ellerman bombs. A data-driven 3D MHD model of this process for a particular AR was investigated by Pariat et al. (2009) and results from that work will be discussed in Section 5. The plausibility of magnetic reconnection for the Ellerman bombs was also investigated by means of MHD simulations by Isobe et al. (2007b) and Archontis and Hood (2009). While the former induced the rise of loops with the Parker instability, the latter considered the evolution of an initially uniform, horizontal flux tube embedded in a convecting layer. Both studies support the scenario posed by Pariat et al. (2004). In addition, Bello González et al. (2013) reports that Ellerman bombs are found near locations where opposite polarities at the photosphere are close to each other. By constructive parameterized models of the plasma heating following reconnection events, they showed via non-LTE radiative transfer calculations heating by reconnection can indeed lead to signatures in Hα that would be detected as Ellerman bombs.

Figure 30:
figure 30

Schematic illustration of serpentine field lines in an emerging flux region. In this picture, development of the undular (Parker) magnetic buoyancy instability leads to loops with a characteristic separation of a few Mm. Magnetic dips occur between adjacent opposite polarities and reconnection between field lines from the polarities lead to plasma heating, which show up as Ellerman bombs in Hα observations. Image reproduced with permission from Pariat et al. (2004), copyright by AAS.

Pariat et al. (2004) found that BPs are usually aligned in series along the orientation of the developing active region and that the typical distance between consecutive BPs is a few Mm. This happens to be consistent with the criterion of the Parker (undular mode, see Section 3.3) instability, which yields λ ≳ 2 Mm for instability under photospheric conditions (see also Shibata et al., 1989a). As already discussed in this section, the action of the convective flow on emerging horizontal field naturally results in serpentine field lines. Nevertheless, the essential physics captured by the Parker instability (namely that mass draining along field lines enhanced buoyancy at loop apices) must be at work.

Figure 31:
figure 31

Hinode Spectropolarimeter (SP) observations of an emerging flux in NOAA AR 10978. The upper and lower right panels show the continuum intensity and line-of-sight magnetic flux densities, respectively. The left panels show the Stokes V signal at ± 260 mÅ from the center of the Fe I 6302.5 Å line. The mixed polarity pattern in the simulated magnetograms of Figure 27 shows some resemblance the observed pattern. Image reproduced with permission from Lites (2009), copyright by Springer.

3.4.5 Convection-driven flux emergence

In the majority of numerical experiments performed for studying magnetic flux emergence, the simulation setup involves the introduction of a buoyant magnetic structure into a stratified layer, either by embedding the buoyant structure near the bottom of the computation domain or by introducing it via time-dependent boundary conditions. In some experiments, such as those that study the non-linear development of magnetic buoyancy instabilities (see Section 3.3), the initial condition may consist of a horizontal field embedded in the stratification in mechanical equilibrium. In the latter case, the development of an instability from a perturbation would naturally lead to some portions of initial magnetic structure to become buoyant and emerge.

An alternative approach to the flux emergence problem is to carry out magnetoconvection simulations and let the convective flows self-consistently form buoyant magnetic structures that eventually emerge into the solar atmosphere. The interactions between the magnetic field and the convective flows lead to magnetic flux that is advected up in upflow regions to create granular scale emergence events. Depending on the presence or absence of a large-scale underlying structure, these emergence events may be organized in a coherent fashion (i.e., they can have similiar orientation so that discrete flux emergence events lead to accumulation of flux into pores and spots) or be independent from one another.

Abbett (2007, see Figure 32) performed magnetoconvection simulations using a realistic EOS from the OPAL project (Rogers et al., 1996) for the near-surface layers of the Sun while imposing a radiative cooling term that was tabulated from radiative magnetoconvection simulations of Bercik (2002). The use of a tabulated cooling function meant that they could bypass solving the radiative transfer equation. The computational domain extended from a 2.5 Mm below optical depth unity in the photosphere to 5 Mm above. The magnetoconvection simulation was initiated with the introduction of a uniform horizontal seed field with a strength of 0.01 G. The simulation showed that convective flows acting on original seed field amplified fields to unsigned flux densities beyond 30 G.Footnote 8 The convective flows also create localized concentrations of kG fields as well as small-scale flux loops that emerged from the convection zone into the atmosphere.

Figure 32:
figure 32

Field line rendering from a magnetoconvection simulation by Abbett (2007). The rendering shows collections of field lines that poke through the surface as granular-scale Ω-loops as well as some U-loops that are trapped by downflow lanes. Image reproduced with permission, copyright by AAS.

Isobe et al. (2008, see Figure 33) carried out magnetoconvection simulations in a dimensionless setup using a ideal gas EOS with constant ratio of specific heats (γ = 5/3). The cooling function at the interface between their model convection zone and chromosphere/corona (measured in terms of density contrast from the top of the convecting layer) was chosen to allow the purely hydrodynamic (i.e., non-magnetic) convection run to develop horizontal flow speeds with Mach numbers M ∼ 0.1, which is not atypical for the quiet Sun. The magnetoconvection experiment was initiated by introducing a uniform vertical field into a snapshot with fully-developed convection. The simulation result showed that the convective flows can wind up the existing field, creating small-scale loops underneath the surface which are then brought up through granular upflows. When this occurs, the emerged flux can interact with the pre-existing, predominantly vertical field in layers a few scale heights above the surface. Magnetic reconnection between the emerging and pre-existing field ‘open’ field leads to jets that emit material up into the higher atmosphere acts as a source of magneto-acoustic waves within the atmosphere. Applying this lesson to the Sun, one may conclude that not all waves in the chromosphere and corona need to originate from the photospheric base but may also come from episodic in-situ reconnection events.

Figure 33:
figure 33

In the simulation of Isobe et al. (2008), the convectively-driven emergence of a granular-scale flux loop reconnects with ambient vertical field. Such reconnection events eject mass into the height atmosphere and provides an in-situ source of magneto-acoustic power (as opposed to waves propagating from photospheric base). The central and left panels show an example of such an event in a 3D magnetoconvection simulation. The right panel shows a 2D counterpart. Image reproduced with permission, copyright by AAS.

The first observational detection of granulation scale emergence events was reported by De Pontieu (2002), who noted that emerged flux elements are buffeted by ambient granular flows. They suggested that such events may be related to the so-called horizontal internetwork fields (HIFs) reported by Lites et al. (1996). The connection between the type of event that De Pontieu (2002) reported and HIFs has been strengthened since the launch of the Hinode satellite. The first reports of detection of small-scale flux emergence from Hinode/SP were by Centeno et al. (2007). They performed inversions on the full Stokes parameters and found the appearance of horizontal field in granular interiors preceding the appearance of vertical flux elements. Thereafter, Ishikawa et al. (2008) and Ishikawa and Tsuneta (2009) studied the so-called transient horizontal magnetic fields (THMFs), which are granulation-scale, transient linear polarization signals with essentially no preferred orientation. THMFs are likely flux emergence events driven by magnetoconvection as opposed to the passage of a large-scale (L ≫ 1 Mm) magnetic structure from beneath the photosphere.

Martínez González et al. (2007) used infrared spectropolarimetric data and inferred the height reached by such small-scale emerging loops. In that paper, they concluded that such granular-sized loops are unlikely to reach the photosphere. Shortly after, from coordinated observations between Hinode/SP and the Dutch Open Telescope, Martínez González and Bellot Rubio (2009) reported that more than 20% of small-scale loops identified in Hinode/SP data manage to reach the chromospheric layers (as detected in Mg I b 5173 Å). Furthermore, they reported a typical time delay of 5 min for photospheric loops to reach the chromosphere. Thereafter, Martínez González et al. (2010) continued this line of investigation from Hinode observations of Fe I 6130 Å doublet, Mg I b 5173 Å and Ca II 3969 A lines to follow the passage of a magnetic loops into the chromosphere (see Figure 34). They estimated that the energy flux provided to the chromosphere by these rising loops is of order 106–107 erg cm−2 s−1. As such, they concluded that small-scale emergence of loops can be an important contributor to the energy budget of the chromosphere and overlying corona.

Figure 34:
figure 34

The passage of a magnetic loop from the photosphere to chromosphere as revealed by Hinode/SOT observations. Image reproduced with permission from Martínez González et al. (2010), copyright by AAS.

It remains uncertain whether they result from local dynamo action or from the re-processing of decayed fields from active regions and/or ephemeral regions, since the latter process can introduce scatter to the orientation of the original field. Furthermore, there is the possibility that the global dynamo and small-scale convective dynamo are intricately linked.

Whereas the simulations of Abbett (2007) and Isobe et al. (2008) focused on the response of the solar atmosphere to convection-driven, granular-scale flux emergence events, the work of Stein and collaborators (Stein and Nordlund, 2006; Stein et al., 2011; Stein and Nordlund, 2012) focused on the evolution of magnetic fields rising over multiple pressure-scale heights in the top layers of the convection zone. In their radiative magnetoconvection simulations, the computational domain encompassed the top 20 Mm of the convection zone, which corresponds to a density contrast of more than 104. Instead of imposing an volume filling field as an initial condition, they began with a 3D model of a fully-developed, radiatively driven convecting flow. Since the simulation (1) uses a realistic equation of state for a gas mixture with solar abundance, and (2) solves the radiative transfer equation in 3D using realistic source functions and opacities (assuming local thermodynamic equilibrium), the model contains solar-like convective flow patterns and a realistic mean stratification. With such a model as a initial condition, they proceeded to inject horizontal magnetic field of uniform strength and orientation in upflow regions that cross the bottom boundary. In this series of papers, they performed simulations varying the field strength of the injected magnetic field from as high as 40 kG to as low as 1 kG. As the horizontal field is injected, it has a typical length scale that follows the size of convective upflows at 20 Mm depth (see Figure 35). However, as the loops eventually reach shallower depths where the pressure-scale height (and associated convective cell size) diminishes, the large-scale loops become locally undulated by the convective downflows. This creates a hierarchy of loops, with smaller loops ‘piggy-backing’ on larger loops that are anchored deeper down.

Figure 35:
figure 35

Vertical velocity patterns from a solar-like model of convection. Blue and green indicate upflows whereas yellow and red indicate downflows. The different panels show the velocity distribution at various depths of the model. The horizontal extent of the computation domain is 48 × 48 Mm2 and the bottom boundary is at 20 Mm below optical depth unity. Image reproduced with permission from Stein et al. (2011), copyright the authors.

Figure 36 shows two synthetic brightness images from the simulation reported in Stein and Nordlund (2012). In this case, the feeding of horizontal magnetic field was limited to a 25 × 25 Mm2 square patch at the bottom boundary (the horizontal extent of the computational domain is 48 Mm2). The horizontal field fed in through the boundary had an orientation of 30 degree relative to the direction of the x-axis and a field strength of 1 kG. They found that, after about two days of continuously injecting flux with this configuration, sufficient magnetic flux had emerged at the photosphere for pores to develop.

Figure 36:
figure 36

Synthetic brightness images of a simulated emerging flux region. Overlaid are traces along horizontal magnetic field vectors at the surface. Image reproduced with permission from Stein and Nordlund (2012), copyright by AAS.

3.4.6 Turbulent transport of magnetic flux

Section 3.2 showed that a buoyant magnetic structure reaching the photosphere must undergo drastic expansion. The horizontal outflows associated with this expansion distribute the emerging field over a wide area (e.g., see Figure 9). Given that the emerged flux is dispersed, one question to ask is how the dispersed flux is transported to give spot-sized coherent flux concentrations. Another question is how does the emerging field get rid of its mass burden to reach the tenuous atmosphere. Observational evidence and numerical simulations suggest that serpentine field lines discussed in Section 3.4.4 may play a central role for both processes.

Consider the scenario depicted in Figure 37. In this simplified 2D picture, a sheet of a magnetic field is predominantly horizontal and is in the process of emerging through the surface. It is undulated due to interactions with convective flows. This field line has a footpoint anchored in the subsurface layers. The other end of the field line connects to the opposite polarity (not shown in the figure). The undulated segment of the field line manifests itself at the photosphere as an alternating sequence of positive and negative polarities. From the expansion of granular upwellings, the positive polarities will be advected to the left side while the negative correlations will be advected to the right side. A consequence of this flux expulsion to the intergranular lanes is flux cancellation by opposite polarities. In the 2D picture, this type of reconnection will create O-loops that carry mass away from the emerging field. In 3D these structures will generally be plasmoids. The removal of mass is an important step for flux to continue its rise into the higher atmosphere and this mechanism provides an efficient way to do so (see Lites, 2009; Centeno, 2012, for observational evidence of this process occurring in EFRs). Other mechanisms include magnetic buoyancy instabilities (see Section 3.3) and ambipolar diffusion (see section 3.8.3). However, the latter is only expected to play a significant role in the chromosphere and not in the near-surface layers of the convection zone and photosphere.

Figure 37:
figure 37

Schematic drawing illustrating the effect of granular flows on magnetic polarities from a serpentine field line. Image reproduced with permission from Cheung et al. (2010), copyright by AAS.

The systematic advection of positive and negative polarities in opposite directions by granular flows provides a systematic correlation that is central to turbulent flux transport. To quantify how correlations between plasma velocity and magnetic field lead to, on average, a transport of flux toward growing spots, Cheung et al. (2010) presented the following mean field approach to examine the numerical simulation of AR formation. Consider a cylindrical coordinate system the z-axis (the axis of symmetry) located at the center of a growing flux bundle that eventually becomes a spot. Consider the rate of change of magnetic flux crossing a circular surface \({\cal C}\) of radius R (for the following discussion, consider \({\cal C}\) to be lying at the photospheric base). Now decompose the magnetic and velocity fields into the sum of azimuthal averages and corresponding fluctuating components, namely \({\bf{B}}(r, \theta, z) = \bar B(r,z) + {\bf{B}}(r,\theta,z)\) etc, where the bar and prime denote an azimuthal average and the fluctuation about the average, respectively. For ideal MHD, the corresponding mean-field induction equation is

$${{\partial \overline {\bf{B}} } \over {\partial t}} = \nabla \times ({\bf{\bar v}} \times {\bf{\bar B}}) + \nabla \times E,$$
(34)

where \(E = \overline {{\bf{v}}\;\;{\bf{B}}}\) is the (azimuthally-averaged) mean-field electromotive force resulting from correlations between the magnetic field and velocity field fluctuations (Krause and Rädler, 1980; Rädler, 1980). The time rate of change of magnetic flux crossing the circular surface \({\cal C}\) is then

$${\dot \Phi _{\cal C}} = \int_{\cal C} {} {{\partial \overline {\bf{B}} } \over {\partial t}}\cdot{\bf{d}}S = \oint_{\partial {\cal C}} {(\overline {\bf{v}} \times \overline {\bf{B}} + E)} \cdot{\bf{d}}l,$$
(35)
$$ = {\dot \Phi _{\rm{m}}} + {\dot \Phi _{\rm{f}}},$$
(36)

where

$${\dot \Phi _{\rm{m}}} = 2\pi R{[\overline {\bf{v}} \times \overline {\bf{B}} ]_\theta } = 2\pi R({\overline v _z}{\overline B _r} - {\overline B _z}{\overline v _r}),$$
(37)
$${\dot \Phi _{\rm{f}}} = 2\pi R{{\cal E}_\theta } = 2\pi R\overline {v_z^{\prime}B_r^{\prime} - B_z^{\prime}v_r^{\prime}}.$$
(38)

Equation (36) indicates that the magnetic flux enclosed within \({\cal C}\) changes as a result of (1) advection of the mean magnetic field \(\overline {\bf{B}}\) by the mean flow \(\overline {\bf{v}}\) (first term on the r.h.s.), and (2) correlations between their fluctuating components (e.g., from the serpentine field line process discussed above).

In their simulation of AR formation, Cheung et al. (2010) found that \({\dot \Phi _{\rm{f}}}\) played a dominant role in transporting magnetic flux inward to facilitate flux growth. They also examine the Lorentz force experienced by the positive and negative polarities and found that the former tended to be accelerated radially inward (for a positive spot) and the latter outward. This finding is consistent with the tethered-balloon model of Spruit (1981a). In his model, emerging magnetic loops anchored below the photosphere behave like tethered balloons and their equilibrium state is reached when buoyant elements hover vertically above their footpoints. The evolution toward this equilibrium is driven by the tension force of the magnetic field.

Rempel and Cheung (2014) extended the work of Cheung et al. (2010) by carrying out radiative MHD simulations of the birth and decay of ARs following flux emergence. They report that the dispersal of vertical magnetic flux from a modeled spot during the first two days of the decay phase is well-described by a self-similar solution to the 2D magnetic diffusion equation (Meyer et al., 1974; Mosher, 1977). While the analyses of Cheung et al. (2010) and Rempel and Cheung (2014) shed some light on how flux is transported by turbulent motions, they do not provide self-consistent mean-field models that take into account the non-linear feedback between the Lorentz force and the fluid velocity.

An attempt at a self-consistent mean-field model of the collapse of weakly distributed vertical flux into large-scale coherent spots was presented by Brandenburg et al. (2013). They performed numerical simulations of the evolution of an initially uniform vertical field (i.e., no large-scale emergence) in a stratified layer with forced turbulence. In their simulations, turbulent motions are imposed by a forcing term in the momentum equation describing random, unpolarized plane waves. They found the spontaneous formation of a large-scale flux bundle with maximum field strengths comparable to equipartition values but less than the 3 kG strengths found in sunspots (Kitiashvili et al., 2010b, also report for formation of pores from their magnetoconvection simulations). The spontaneous formation of such a flux concentration is reported to be consistent with the presence of the so-called negative effective magnetic pressure instability (NEMPI). This instability arises in a mean-field description, in which the presence of turbulent motions acting on a magnetic field leads to a suppression of the hydromagnetic pressure and tension. In a related study, Warnecke et al. (2013) carried out direct numerical simulations of stratified forced turbulence acting on an initially horizontal field. Their simulations include an idealized ‘coronal’ envelope in which the turbulent forcing term in the momentum equation is suppressed. From an initially uniform magnetic field, they report that a bipole magnetic pair spontaneously forms within 1–2 turbulent diffusion times (estimated from the size of the computational domain and the effective turbulent diffusivity) and then subsequently decays within half a diffusion time.

3.5 The role of large-scale subsurface structure

Is the role of flux emergence simply to load the near-surface layers of the Sun with magnetic flux, and then to let turbulent effects near the surface take over to cause collapse into a spot? Or is pre-existing large-scale structure in the subsurface layers essential for spot formation? These questions must be answered in the context of observed properties of sunspots and ARs. For instance, observational studies (e.g., Zwaan, 1978, 1985) indicate that spot formation does not result from the collapse of a pre-existing distributed field in the absence of large-scale flux emergence events. Robust statistical trends of ARs may also give clues about the subsurface structure of ARs. Well-known statistical properties of ARs include the equatorward migration of active latitudes during a cycle, Hale’s polarity rule, Joy’s law of AR tilt and the asymmetry between field strengths of leading and trailing polarity spots. These statistical trends have motivated a long line of studies of the rise of buoyant magnetic flux tubes in global spherical domains, either using the thin flux tube approximation (Fan et al., 1993; Moreno-Insertis et al., 1994; Caligari et al., 1995; Weber et al., 2011) or by peforming full 3D MHD simulations (Fan, 2008; Jouve and Brun, 2009; Jouve et al., 2013; Pinto and Brun, 2013; Nelson et al., 2014). Despite the inherent limitations of the thin flux tube approximation (see Section 2.2 of Fan, 2009b, for a detailed discussion), this class of models still holds appeal due to their remarkable success in reproducing (and providing a physical basis for) many of the statistical properties of ARs.

Rempel and Cheung (2014) carried out a number of radiative MHD simulations to test the influence of subsurface flows and large-scale magnetic structure on the spot formation process. In a control experiment, they imposed the kinematic rise of a semi-toroidal magnetic tube through the bottom boundary of a Cartesian computational domain capturing the top 16 Mm of the convection zone. The flux tube has an axial field strength of 21.2 kG at that depth and a total toroidal flux of 1.7 × 1022 Mx. After the semi-torus had passed through the bottom boundary, they switched off the imposed upflows (500 m s−1) at the footpoints of the torus. In this experiment, the associated emerging flux region at the photosphere produced a pair of spots (or large pores that do not have fully developed penumbrae). In a second experiment, a magnetic semi-torus with a lower axial field strength (10.6 kG) but with the same total flux was used. This experiment also produced a pair of spots in the emerging flux region. In a third experiment, the upflows at the torus footpoints were maintained even after the semi-torus had crossed the bottom boundary. Small pores formed in this third experiment but the sustained flux of mass and enthalpy pumped through the footpoints of the torus led to horizontal outflows that prevented the coalescence of these pores into large coherent spots.

A robust feature of simulations modeling the rise of magnetic flux tube under the influence of rotation is the presence of retrograde flows at the apices of the rising loops (see references above as well as Abbett et al., 2001). Such a retrograde flow of a few hundred m s−1 is a consequence of angular momentum conservation. Motivated by this finding and observational evidence for asymmetry between leading and trailing polarities in observed ARs, Rempel and Cheung (2014) repeated their control experiment with one important change. To mimick the presence of a retrograde flow, they superimposed a toroidally-aligned flow of 500 m s−1 on top of the advecting upflow at the bottom boundary. The presence of this flow resulted in an asymmetric model AR with a leading polarity spot that is more coherent than the trailing polarity spot. A comparison between the different experiments carried out in their study suggests that subsurface structure (in particular large-scale flows) has an important influence on the ability of the emerging flux region to yield spots.

Thin flux tube models and flux emergence models that assume the existence of isolated flux tubes as initial or boundary conditions do not addresss the question of how such structures are formed. Recent work on global dynamos is beginning to address this issue. Global anelastic MHD simulations by Nelson et al. (2014) show that large-scale coherent flux bundles can spontaneously form from a dynamo-generated field (albeit in a model rotating at 3 times solar rotation), and that convective flows help advect the magnetic bundles into rising Ω-loops. One limitation of global flux tube (and dynamo) simulations is that they do not capture the top ten or so pressure scale-heights of the convection zone. Complementarily, numerical models of magnetic flux emerging at the surface do not capture the deeper layers treated by global models. Going forward, models capturing the coupling between the deeper and the shallower layers of the convection zone are likely essential for addressing the outstanding questions posed above.

3.6 The role of magnetic twist in flux emergence

A number of observational results suggest that emerging flux regions (especially on AR-scales) often embody magnetic twist to various degrees. On the observational side, vector magnetogram observations of ARs using photospheric spectropolarimetric techniques show that at least some ARs have systematic current distributions (e.g., Leka et al., 1996). Active regions with so-called ‘magnetic tongues’ in the morphology of their magnetic polarities in longitudinal magnetograms is also suggestive of underlying large-scale twist (Canou et al., 2009; Archontis and Hood, 2010; Luoni et al., 2011, see Figure 38 for an example). Furthermore, there is a trend that, within a given cycle, ARs in the northern and southern hemispheres have opposite signs of twist (Pevtsov et al., 1995; Longcope et al., 1998; Hagino and Sakurai, 2004). In the corona, the detection of loop systems with sigmoidal morphologies in X-ray (Rust and Kumar, 1996; Canfield et al., 1999) and EUV (Sterling et al., 2000) observations is suggestive of magnetic twist in ARs.

Figure 38:
figure 38

Evidence for twisted flux tube emergence. The top row shows three SoHO/MDI magnetograms of NOAA AR 10808. The positive and negative polarities have ‘magnetic tongue’ morphology. The bottom row shows three synthetic magnetograms from the simulation of the emergence of a twist flux tube. The striking resemblance in morphology suggests a twisted flux rope structure for NOAA AR 10808. Image reproduced with permission from Archontis and Hood (2010), copyright by ESO.

The examination of magnetic twist in emerging flux and developed ARs is important because the degree of magnetic twist is a proxy for free magnetic energy available to power solar eruptions and flares (Berger and Field, 1984). When twisted magnetic flux emerges onto the solar surface, the transport of current-carrying magnetic field (and its associated magnetic helicity) plays a crucial role in the build-up of this free energy. While some of concepts in the following sections have broader relevance than the immediate topic of this review article, it is nevertheless important to introduce them to reveal their connection with flux emergence.

3.6.1 Idealized, twisted magnetic flux tubes

Magnetic twist is most conveniently defined when considering idealized flux tubes. Consider an axisymmetric, cylindrical magnetic flux tube with \({B_l}(r)\hat l\) and \({B_\theta }(r)\hat \theta\) representing the longitudinal and transverse components of the magnetic field, respectively. By definition, an untwisted flux tube has zero transverse field (Bθ = 0). The presence of a systematic (i.e., azimuthally-averaged) transverse field 〈Bθ〉 about the tube axis indicates the presence of currents in the tube. The net current contained within a circular cross-section of radius r from the tube axis is given by

$$I(r) = {c \over {4\pi }}\int {{{(\nabla \times {\bf{B}})}_l} \cdot } d{\bf{S}} = {c \over 2}\int_0^{2\pi } {{B_\theta }r \, d\theta = {{cr} \over 2}\langle {B_\theta }(r)\rangle.}$$
(39)

In numerical experiments modeling the dynamics of twist flux tubes, the magnetic profiles Bl(r) and Bθ(r) are often truncated (i.e., set to vanish) beyond some radial distance. In such a case, there will be a radius R such that I(r) = 0 for r > R. Nevertheless, the flux tube is still considered to be current-carrying (i.e., have net twist) if I(r) ≠ 0 for r < R.

As pointed out by Murray et al. (2006), a common choice of the magnetic distribution for idealized axisymmetric flux tubes in numerical experiments is

$${B_l}(r) = {B_0}\exp (- {r^2}/{a^2}),$$
(40)
$${B_\theta }(r) = qr{B_l}(r),$$
(41)

where q and a are parameters that specify the twist and thickness of the tube (e.g., Fan et al., 1998, 1999; Cheung et al., 2007a). Other magnetic profiles have also been used (e.g., Moreno-Insertis and Emonet, 1996; Emonet and Moreno-Insertis, 1998; Magara, 2001) in the past. However, as demonstrated by Murray and Hood (2008), the exact distribution is not as important as the strengths of the field components.

The first MHD simulations of the buoyant rise of a (initially stationary) magnetic flux tube by Schüssler (1979) showed that an untwisted tube fragments two counter-rotating vortex rolls after rising a distance comparable to its initial diameter. The horizontal pressure gradient that leads to the horizontal expansion of the initial tube combined with the net vorticity of each of the fragments lead to downward directed aerodynamic lift forces which halt the further rise of the magnetic structure (Longcope et al., 1996). Moreno-Insertis and Emonet (1996) and Emonet and Moreno-Insertis (1998) showed that, given sufficient magnetic twist, magnetic tension can suppress the fragmentation.Footnote 9 Borrowing the terminology used for comparing inertial and tension forces in rising air bubbles, Emonet and Moreno-Insertis (1998) defined the magnetic Weber number as

$$We = {{{v^2}\varrho } \over {B_t^2/4\pi }},$$
(42)

where v is the rise velocity of the tube. The criterion for maintaining coherence of the flux tube is We ≲ 1.

The extension of this type of investigation to 3D MHD simulations indicates that variations in the third dimension can have an impact on the requirement for rising tubes to remain coherent. In their 3D simulations of buoyant magnetic flux tubes Dorch and Nordlund (1998) examined the effect of a randomly varying transverse field on flux tube dynamics. They report that even a weak, random transverse field (with field strength that is an order of magnitude smaller than the longitudinal field) is sufficient to let the tube maintain a more-or-less coherent rise. Abbett et al. (2000) introduced variations in the third dimension in terms of a varying profile of the density (also specific entropy) deficit along the tube. The choice of the profile leads the tube to develop into an Ω-loop structure. They report that, even when the twist of the initial tube is below the critical threshold, the tube can remain coherent provided it has sufficiently large curvature at the apex. In addition to these effects, simulations that include the effects of rotation show that presence of the Coriolis force helps suppress the fragmentation at the tube apex (Wissink et al., 2000b; Abbett et al., 2001).

Krall et al. (1998) carried out the first idealized 2.5D simulations of twisted flux tube emergence from the convection zone into the corona and concluded that the presence of twisted aided the rise of the tube. Murray et al. (2006) carried out a parameter study of 3D flux emergence simulations varying the twist and magnetic field strength of the initial submerged flux tube (see Figure 23). In their experiments, the use of a twist parameter of q = 0.1 (for a = 2.5, see Eqs. (40) and (41)) results in the flux tube being trapped below the photosphere. The field strength resulting from the evolution was found to be too small to initiate magnetic buoyancy instabilities to launch the field into the overlying corona (see Section 3.3.5). When they increased the twist, the rising field retained compactness and successfully emerged. Toriumi et al. (2011) also performed a parameter study of 3D flux emergence experiments and reported results in agreement with Murray et al. (2006). These results, together with considerable observational evidence (see following section), motivates the choice of using twisted flux tubes in numerical simulations of flux emergence.

3.6.2 The helical kink instability in emerging flux tubes

There is an extensive body of research on numerical MHD simulations of the rise and emergence of twisted magnetic flux tubes into the solar atmosphere. The first numerical models of the emergence of twisted flux tubes were motivated by the association of delta sunspots with flares (e.g., Kurokawa et al., 1987; Tanaka, 1991; Ishii et al., 1998; Liu and Zhang, 2006) and by the discovery of sigmoidal structures in X-ray observations (Acton et al., 1992; Canfield et al., 1999). Early studies that attempt to reproduce such structures invoke the nonlinear development of the helical kink instability. Here we do not delve into the detailed derivation of the linear stability analysis of Linton et al. (1996). We simply summarize their finding that when the twist q in Eq. (41) exceeds the critical value

$${q_{{\rm{cr}}}} = {a^{ - 1}},$$
(43)

there exists unstable eigenmodes which grow. Simulations of the nonlinear development of the helical kink instability show the writhing of the tube axis as the instability ensued (Matsumoto et al., 1998; Fan et al., 1998, 1999). Matsumoto et al. (1998) modeled the emergence of a helical kink unstable flux tube from the solar interior into the corona (both modeled as idealized plane-parallel stratifications). Figure 40 shows a visualization of the nonlinear development of the helical kink instability modeled by Fan et al. (1998). These types of studies suggest that kinked structures resulting from the instability are a possible cause for sigmoidal loops observed in X-rays (see Figure 39). To allow the simulation to proceed, they used field strengths of the initial subsurface tube with Alfvien speeds comparable to the local sound speed.

Figure 39:
figure 39

3D rendering of a simulation of a twisted magnetic flux rope undergoing the helical kink instability. While the helical instability is developing, the crests of the flux tube are emerging into the idealized (initially plane-parallel) model atmosphere. Image reproduced with permission from Matsumoto et al. (1998), copyright by AAS.

Figure 40:
figure 40

mpg-Movie (736.40625 kB) Still from a movie showing Volumetric rendering of a rising flux tube experiencing the developing of the helical kink instability (Fan et al., 1998). (For video see appendix)

Fan et al. (1999) simulated the development of the instability in twisted flux tubes buoyantly rising through the solar interior and found that the development of the kink instability increases the buoyancy at the tube apex. Furthermore, horizontal cross-sections of the vertical magnetic field shows a compact bipole pair with an orientation that deviates more than 90 deg from the axial orientation of the initial tube. The horizontal components of the magnetic field near the polarity inversion line of the compact bipole pair is highly sheared. These results led them to conclude that the helical kink instability is a plausible mechanism for the formation of delta spots. Since most observed active regions do not possess delta-type spots, one is led to conclude that the majority of buoyant flux tubes that manage to reach the surface are weakly twisted.

3.6.3 Rotational and shearing motion of photospheric footpoints due to twisted flux rope emergence

In the presence of magnetic twist that is coherent over the scale of sunspots and ARs, it is almost inevitable that an emerging structure exhibits rotational motions at the photosphere. The basic underlying reason is the nature of the Lorentz force. Consider a closed curve lying at a constant height in the atmosphere enclosing some point \({\cal P}\). The line integral of the gas and magnetic pressure gradient forces about this curve vanish. That is,

$$\oint {{\bf{r}} \times \nabla } (- {p_{{\rm{gas}}}})\cdot{\bf{d}}l = \oint {{\bf{r}} \times \nabla } \left({ - {{{B^2}} \over {8\pi }}} \right)\cdot{\bf{d}}l = 0,$$
(44)

where r is the displacement vector of a point in the curve from \({\cal P}\). The same is not true for the line integral of the magnetic tension contribution to the Lorentz force. So, in general, magnetic tension provides a net torque

$${t_z} = \oint {} \left[ {{\bf{r}} \times {1 \over {4\pi }}({\bf{B}}\cdot\nabla){\bf{B}}} \right]\cdot{\bf{d}}l,$$
(45)

which cannot be balanced by pressure gradient forces. This unbalanced torque drives rotational motions at the photosphere.Footnote 10

Consider an idealized twisted flux tube that is initially horizontal below the surface and endowed with a density deficit such that a localized segment of the tube develops into the apex of an Ω-loop. Let the differential density deficit be symmetric about the tube apex. As the Ω-loop emerges through the photosphere, there will be two opposite polarity concentrations. Due to the symmetry of the problem, conjugate points (in a horizontal plane) about the point of symmetry (assumed to be at the origin) are related by Bx(−r0) = Bx(r0), By(−r0) = By(r0), and Bz(−r0) = −Bz(r0). Since magnetic tension is quadratic in B, the torques evaluated about points r0 and −r0 will have equal sign and amplitude.

The aforementioned scenario is a common setup in numerical experiments of flux emergence. Fan (2009a) performed such an experiment and examined the distribution of the vertical component of vorticity. Distributions of wz at the photospheric level at two different times during the emergence event are displayed in Figure 41. The same sense of vorticity in opposite polarity patches of emerging (or emerged) twisted flux tubes has important ramifications for the continued evolution of the magnetic field in the atmosphere. This twists up magnetic field lines in the corona in such a way that certain field lines develop sigmoidal shapes. For a detailed analysis of how photospheric flows during flux emergence can be decomposed into translational, rotation, and shearing components, see the paper by Magara (2006).

Figure 41:
figure 41

Distribution of the vertical component of vorticity (wz) at the photospheric level at two different times in a numerical model of twisted flux tube emergence. In both panels, the two yellow patches are co-incident with the central portions of the twisted tube. Through they have opposite polarity, they have the same sign and amplitude of wz. Image reproduced with permission from Fan (2009a), copyright by AAS.

Depending on the exact configuration of opposite polarity patches in an emerging flux region, vortical flows around opposite polarities will manifest themselves as shear flows. For instance, when opposite polarities are pressed against each other in a compact configuration (e.g., in delta spots, where there is a high gradient of Bz across the polarity inversion line), there tends to be shearing motion across the polarity inversion line. As is in the case of rotation with the same sense by a pair of well-separated sunspots (e.g., Fan, 2009a), shearing flows about a polarity inversion line in a compact AR leads to increased twist in the overlying (emerged) coronal field (e.g., Manchester IV et al., 2004, see Figure 42 of this article). In the case of Fan (2009a), the build-up of twist in the corona is interpreted as the result of propagating torsional Alfvén waves. In the case of Manchester IV et al. (2004), the twisting up of the coronal field is interpreted as the result of the upward propagation of shear Alfvén waves (see also Manchester IV, 2001, 2007). Even in flux emergence simulations that include the vigorous convective flows at the surface, the emergence of sufficiently twisted flux can drive systematic shear shows (Fang et al., 2012). Ultimately, the common agent between these scenarios is the Lorentz force (specifically magnetic tension) providing torques.

Figure 42:
figure 42

Numerical simulation of the evolution of magnetic field from twisted flux rope emergence. Panels (a)–(f) show a vertical plane normal to the axis of the initial submerged tube. The direction of the magnetic field in the plane is depicted with solid white lines. Panels (a)–(c) show the velocity Ux (out of the plane), while panels (d)–(f) show color images of the shear angle measured in degrees. Panels (g)–(i) show 3D visualizations of magnetic field lines and the vertical field strength Bz at the photosphere. Panel (h) shows Ux at the photospheric level along with arrows indicating the flow velocity. Isosurfaces of Ux are shown colored in blue (−19 km s−1) and red (+19 km s−1). Image adapted from Manchester IV et al. (2004), courtesy of W. Manchester.

The winding up of the coronal field from the emergence of a twisted Ω-loop is initiated from the dramatic expansion of field as it reaches near-surface layers of the convection and the overlying atmosphere. Consider an initially horizontal, twisted flux tube with no variation along its axis. Now let a localized section of the tube expand. In ideal evolution, the total number of windings of field lines in the tube is preserved. Parker (1974) showed that the flux tube evolves in such a way that windings migrate from the rest of the flux tube toward the expanded portion. This result can also be described in terms of field line torsion α. In a force-free field (i.e., ∇ × B = αB), it can be shown that B • ∇α = 0. In order words, α is constant along field lines. Even before reaching an equilibrium state, the field line torsion can be locally defined as α = (∇ × B) · B/B2. Longcope and Klapper (1997) showed that in a thin flux tube (i.e., a thin flux bundle), a gradient of α along the tube leads to an axial torque. The torque over a segment of length dl given by

$${\tau _l} = {1 \over {16}}{a^4}B_l^2{{\partial \alpha } \over {\partial l}}dl,$$
(46)

where a is the radius of the thin flux tube and Bl the axial field strength (Longcope and Welsch, 2000; Fan, 2009a).

A magnetic field evolving via the Lorentz force attempts to reach a force-free state (if such a state exists) by equilibrating α along individual field lines. Fan (2009a) investigated the distribution of α along field lines in a flux emergence simulation, and found that the magnitude of α along specific field lines is much smaller in the expanded portions of the emerged coronal field than at points that remain submerged in the convection zone. This is due to the dramatic expansion and lengthening of field lines that expand into the coronal volume. Longcope and Welsch (2000) developed an instructive model of this process of how the equilibration occurs in terms of the shunting of a current (i.e., the current of the submerged tube is shunted at convection zone/corona interface). The shunting of the current launches torsional Alfvén waves (see also appendix F of Parker, 2007) that eventually remove the mismatch between the current/twist in the convection zone and corona. A characteristic of this model is that there is a prolonged equilibration phase (lasting more than a day for AR scales) of sunspot rotation to transfer twist into the corona. This is a plausible explanation for why the shearing term for magnetic energy and helicity energy injection into the corona is more sustained than the contribution by emergence (see also Section 3.6.5).

3.6.4 Quantification of twist in observed ARs

Although the concept of magnetic twist is easily defined for idealized, axisymmetric flux tubes, the same definition cannot be directly applied to quantifying the intrinsic twist or ‘knottedness’ (Moffatt, 1969) in actual non idealized, non-symmetric emerging flux and active regions. There are two reasons for this difficulty. Firstly, it is (as yet) observationally not possible to measure the full 3D distribution of the magnetic field from the photosphere to the corona. Without the full 3D field distribution, there are no direct measurements of volumetric quantities that relate to the twist of an AR. Secondly, in any AR with a morphology more complex than a simple bipole, the concept of a clearly defined axis is ill-defined. For these reasons, indirect methods and proxies has been developed to quantify the ‘twist’ present in observed ARs.

An example of a method that relies only on photospheric field measurements is the use of current helicity hc = j · B. Since vector magnetograms at a single height only give the contribution to hc from jzBz, the latter is often used as a proxy of the total current helicity (Pevtsov et al., 1995; Hagino and Sakurai, 2004). Measurements of this quantity based on ground-based photospheric vector magnetograms indicate that, on average, the ARs in the northern and southern hemispheres have opposite signs of twist.

Some methods for the quantification of twist are based on the assumption that one can reconstruct the 3D coronal field of an using only vector magnetograms (at the photospheric and/or chromospheric levels). The 3D field configuration B(x, y, z) can then be used for calculating global, volumetric quantities that relate to the amount of field line twist or current present in the magnetic configuration. Early attempts considered force-free field extrapolations from magnetograms. By definition, a force-free field is one such that the Lorentz force c−1j × B vanishes everywhere in the volume. This implies that the current density j is proportional to B, or

$$\nabla \times {\bf{B}} = \alpha {\bf{B}},$$
(47)

where α is constant along magnetic field lines and is a measure of field line torsion. When α = 0, this reduces to the current-free condition and the field is potential. When α is uniform but non-zero, the field is a linear force-free field. By fitting linear force-free fields to photospheric magnetograms of ARs, Pevtsov et al. (1995) found that, within the same sunspot cycle, ARs in the northern and southern hemispheres have, on average, opposite senses of magnetic twist, a conclusion consistent studies from surface current helicity measurement (Hagino and Sakurai, 2004). The linear force-field parameter α has amplitudes of order 1 × 10−8 m−1 to 5 × 10−8 m−1 (these amplitudes apply to both positive and negative twist) though the scatter about these values is large.

3.6.5 Quantification of magnetic helicity injection

A rigorous quantity measuring the twist inherent in magnetic fields is magnetic helicity. Simply defined, the magnetic helicity over a certain volume V is

$$H = \int_V {{\bf{A}} \cdot {\bf{B}}dV,}$$
(48)

where A is the vector potential. When Bn = 0 on the boundary ∂V, this quantity is gauge-invariant and has direct correspondence to the Gauss linking number, which measures the pairwise linkage of field lines within the volume V (Moffatt, 1969). Under ideal MHD evolution, H is time-invariant, which means that the linkage of field lines is conserved. Even when magnetic reconnection is present in the system, it has been shown that the helicity decreases at a much lower rate than magnetic energy (Berger, 1984). Such is the appeal of magnetic helicity. Furthermore, the minimization of the magnetic energy v B2/8π dV subject to constant helicity yields the Euler-Lagrange equation ∇ × B = α0B, where α0 is the linear-force free parameter (Woltjer, 1958). When the Euler-Lagrange equation is derived without the constraint of a fixed magnetic helicity, one arrives at the result, for a given distribution of normal field Bn on ∂V, the potential (i.e., current-free) configuration is the one with the lowest energy state.

When there is magnetic flux crossing the boundary ∂V, helicity as defined by Eq. (48) is not invariant under a gauge transformation A′ = A + ∇φ. Since solar applications most often involve volumes with magnetic flux crossing photospheric boundaries, the concept of relative helicity was developed (Berger and Field, 1984; Finn and Antonsen Jr, 1985). The relative helicity defined as

$${H_R} = \int_V {({\bf{A}} \cdot {\bf{B}})dV} - \int_V {({{\bf{A}}_0} \cdot {{\bf{B}}_0})dV}$$
(49)

is gauge-invariant. Here B is the actual magnetic field configuration in V. A0 and B0 are the reference vector potential and magnetic field (∇ × A0 = B0) with \({{\bf{B}}_0}\;\;\hat n{\;_{\partial V}} = {\bf{B}}\;\;\hat n{\;_{\partial V}}\).

Irrespective of HR, the free energy of a magnetic configuration can be defined as

$${E_F}({\bf{B}},{\bf{P}}) = {1 \over {8\pi }}\int_V {} \vert{\bf{B}}{\vert^2} - \vert{\bf{P}}{\vert^2}dV,$$
(50)

where B is the observed or modeled 3D field distribution in volume V and P is the corresponding potential field distribution with the same Bn distribution over ∂V. Since solar flares and eruptions are magnetically driven (see review article by Shibata and Magara, 2011), EF is often taken as the upper bound for energy available to be released during these episodes. For relative helicity calculations, the potential field is often chosen as the reference field, i.e., A0 = Ap and B0 = P.

For both HR and EF, the 3D distribution of the magnetic field in the volume must be known. Since it is not yet possible to measure the full 3D field distribution in the optically thin corona, many have attempted non-linear force-free field (NLFFF) extrapolation to address the problem. It is outside the scope of this article to discuss the applicability and difficulties extrapolation methods and we refer the reader to DeRosa et al. (2009) and references therein for an in-depth review. Without knowledge of the full 3D magnetic field in the corona, an alternative approach is to measure the flux of HR crossing the photosphere during the birth and subsequent evolution of an AR. The rate of relative helicity transport into a volume V is

$${{d{H_R}} \over {dt}} = 2\oint_{\partial V} {} ({{\bf{A}}_p} \times {\bf{E}}) \cdot \hat ndS,$$
(51)

where Ap is the vector potential field of the potential field satisfying ∇ · Ap = 0 and \({{\bf{A}}_p}\;\;\hat n{\;_{\partial V}}\; = 0\) (Berger and Field, 1984). Under ideal evolution, E = −v × B so that

$${{d{H_R}} \over {dt}} = 2\oint_{\partial V} {} [({{\bf{A}}_p} \cdot {\bf{B}}){\bf{v}} - ({{\bf{A}}_p} \cdot {\bf{v}}){\bf{B}}] \cdot \hat ndS.$$
(52)

It is important to note that the integral is over the closed surface ∂V. Applications of Eq. (52) to solar observations are generally limited to evaluation of the integral at the photospheric surface (with V as some volume of interest above photosphere). As such, the interpretation of \(\int {{{d{H_R}} \over {dt}}dt}\) as helicity build-up in the corona (e.g., above an AR) implicitly assumes that helicity flux through the rest of the bounding surface ∂V has zero net contribution. With this assumption the helicity flux equation is simply

$${{d{H_R}} \over {dt}} = 2\int {} [({{\bf{A}}_p} \cdot {\bf{B}}){v_n} - ({{\bf{A}}_p} \cdot {\bf{v}}){B_n}]dS,$$
(53)

where the integral is evaluated over the region of interest at the photosphere. The last equation is the one commonly used in measurements of helicity flux from emerging active regions (e.g., Chae, 2001; Green et al., 2002; Kusano et al., 2002; Nindos et al., 2003; Jeong and Chae, 2007; Liu and Schuck, 2013). In the literature, the first and second terms in Eq. (53) are usually referred to as the emergence and shear terms, respectively.Footnote 11

Magara and Longcope (2003) examined the magnetic energy and helicity fluxes through the photospheric layer in a numerical model of the emergence of a twisted flux tube (see Figure 43). They found that in the early emergence phase, the flux of helicity and magnetic energy are both dominated by the emergence term (i.e., bodily transport in the vertical direction). Thereafter, in the later phases of flux emergence, both quantities are dominated by the shear terms over a longer time period. Further numerical experiments of twisted flux tube emergence (Fan and Gibson, 2004; Fan, 2009a) reinforce this result. Fan (2009a) interpreted this result in terms of the current shunting model of Longcope and Welsch (2000). In short, the stretching of field lines into the corona leads to a decrease of field line tension in the coronal field. The mismatch between field line tension of the field in the corona and convection zone leads to a propagation of twist into the corona (see also Manchester IV, 2007). For AR scales, this process occurs on time scales of days, and is consistent with the slow, sustained injection of helicity by the shear term. See Section 3.6.3 for a further discussion on the Longcope and Welsch (2000) model.

Figure 43:
figure 43

Time profiles of the magnetic energies and relative magnetic helicity as well as their fluxes as computed at the photospheric surface from a 3D numerical simulation of the emergence of a twisted flux tube. Image reproduced with permission from Magara and Longcope (2003), copyright by AAS.

While Eq. (52) for helicity flux through the photosphere is gauge-invariant, it does not by itself offer physical insight on the means by which helicity density is transported within the surface of integration. A straightforward interpretation of the integrand GA = 2[(Ap · B)vn − (Ap · v)Bn] as a helicity flux density yields spurious results since even simple translational motions that do not inject twist will yield apparent helicity flux density maps with large positive and negative contributions (though the surface integral will sum to zero, see Pariat et al., 2005, and references therein). An alternative flux density of magnetic helicity transport was developed by Pariat et al. (2005, called GΦ) and tested by Dalmasse et al. (2014), which can be interpreted as a surface flux density of helicity flux resulting from the twisting of elemental flux tubes and the emergence of such elemental flux tubes. This requires knowledge of the footpoint mapping of field lines, which either requires the 3D magnetic field to be known throughout the coronal volume and/or inference of the footpoint mapping from imaging data of coronal loops. The 3D coronal field reconstruction method of Malanushenko et al. (2012) may be well-suited for this method since the method uses fitted loops from imaging data to constrain the coronal field reconstruction.

3.6.6 Genesis of flux ropes in the solar atmosphere

Do flux ropes in the solar atmosphere result from bodily emergence from the convection zone or are they formed due to post-emergence evolution? Numerous simulations have shown that it is difficult for a coherent flux tube to rise bodily into the corona as a whole due to the heavy plasma that is trapped at the bottom concave (or U-shaped) portions of the winding field lines (e.g., Fan, 2001b; Manchester IV et al., 2004; Archontis et al., 2004; Magara, 2006). While the upper parts of the helical field lines of the twisted flux tube expand into the atmosphere and the corona, the U-shaped parts of the winding field lines remain largely trapped at and below the photospheric layer, and the center of the original tube axis ceases to rise a couple of pressure scale heights above the photosphere (see Figure 20 for an example).

Whether the tube axis can emerge or not is related to important questions such as the fraction of the initial magnetic flux that reaches the upper atmosphere and whether filaments are bodily transported from the convection zone or a product of shearing and flux cancellation at the photosphere. Clarifying the condition for the emergence of axial field is therefore and important topic for simulation studies. Instead of horizontal, cylindrical flux tubes commonly used in simulation studies, Hood et al. (2009) and MacTaggart and Hood (2009b) adopted an semi-toroidal, twisted flux tube as the initial condition. In this case, draining of plasma along the legs of the torus faciliated the emergence of the axial field. In the left panel of Figure 44, it can be seen that the axial field line (in red) has risen above the photosphere. Even without the passage of the axial field line into the solar atmosphere, the emergence of the top portion of a twisted magnetic tube can still result in the formation of one or more flux ropes in the coronal volume. The process by which this occurs is discussed in Section 4.4.

Figure 44:
figure 44

Magnetic fields in the simulations of the emergence of toroidal flux tube from MacTaggart and Hood (2009b), with weaker field (left) and stronger field (right). The red line in each panel indicates the axial field line. The grey horizontal slice shows the distribution of vertical magnetic field in the photosphere. In the stronger field case (right), twisted field lines above the photospheric layer yield a flux rope-like structure. Images reproduced with permission, copyright by ESO.

What do observations say about the relation between the emergence of helical flux tubes (including their axes and bottom halves) and flux ropes in the chromosphere and corona? Based on photospheric magnetogram sequences from the Hinode/SP, Okamoto et al. (2008, 2009) presented a picture of a helically twisted flux tube rising coherently (like a solid rod) from below the photosphere to above optical depth unity. This reported flux emergence event occurred in the neighborhood of an active region with overlying field. They noted the so-called ‘sliding doors’ effect, which is a widening of the channel between opposite polarity field followed by a narrowing of the channel. Another observational feature of this event was the apparent rotation of the transverse field with respect to the polarity inversion line during the reported passage of the twisted tube. The inferred physical scenario was taken to be associated with the subsequent formation of a filament above the region of interest. Lites et al. (2010) used combined Hinode/SOT and TRACE observations to study a different emergence flux region and also concluded that the formation of the filament channel in that AR was due to the emergence of a helical field.

To investigate the plausibility of the twisted flux tube emergence explanation, MacTaggart and Hood (2010) performed numerical simulations of an initially twisted flux tube emerging into an atmosphere with a pre-existing arcade. They were able to reproduce both effects. Based on the simulations, they reported that the sliding doors effect has different phases. First of all, the rise of the twisted tube and its lateral expansion pushes aside plasma threaded by pre-existing arcade field. However, since Bz from the twisted tube begins to contribute to photospheric vertical field, this phase would not be detected as the beginning of the widening phase of the sliding doors effect. As the tube continues to push into the subadibatically stratified photosphere, it continues to expand laterally, which pushes apart the emerged opposite polarities from the polarity inversion line. It is this phase that MacTaggart and Hood (2010) associate with the widening of the observed sliding doors effect. The continued passage of the bottom half of the flux tube into the atmosphere (aided by magnetic buoyancy instabilities, see Section 3.3) leads to the narrowing of the channel. The simulations also reproduced the rotation of the photospheric transverse field with time. MacTaggart and Hood (2010), therefore, concluded that the twisted flux tube emergence scenario has some observational features that are consistent with the scenario presented by Okamoto et al. (2008) and Okamoto et al. (2009).

Vargas Domínguez et al. (2012) disputed the conclusion of MacTaggart and Hood (2010) by pointing out that their numerical simulations predict two additional observable effects which did not fit observations of this supposed flux emergence episode. First of all, they pointed out that the emergence of a twisted flux rope increases the unsigned flux about the polarity inversion line. This occurs in the simulations but not in MDI magnetograms. The second is the lack of evidence of shear flows from the observations, which are a robust feature of the emergence of twisted flux tubes (see Section 3.6.3 for a discussion of the physical mechanism that drives shear flows). Vargas Domínguez et al. (2012) argued that, in the absence of these two observable features, the event reported by Okamoto et al. (2008) and Okamoto et al. (2009) is more consistent with flux cancellation about the polarity inversion line.

Irrespective of the eventual consensus, the purpose of recounting the scientific debate above is to illustrate the potential pitfalls of interpreting time sequences of magnetograms as a 3D representation of the subsurface structure of emerging magnetic fields. This does not mean that temporal evolution of photospheric fields lack valuable information about subsurface structure. Rather, the interpretation should take into account the physics that govern evolution of structures rising in a strongly stratified layer.

3.7 Resistivity and thermal conductivity models

3.7.1 Effect of resistivity models

In Section 2.2.3, it was briefly mentioned that the choice of magnetic resistivity (and, hence, magnetic diffusivity) has an impact on how emerging magnetic flux interacts with surrounding field. The physics of magnetic reconnection itself is a very large topic so covering the entire topic is beyond the scope of this review. Here, we briefly review the effects of resistivity models in MHD modeling of flux emergence. For more comprehensive review of recent progresses of reconnection studies including the kinetic effects, see Yamada et al. (2010).

The classical Spitzer conductivity σSpitzer in fully ionized plasma due to the Coulomb collision of electrons and ions (Spitzer Jr, 1962) is given by

$${1 \over {{\sigma _{{\rm{Spitzer}}}}}} \approx {10^{ - 16}}{T \over {{{10}^6}{{\rm{K}}^{ - 3/2}}}}({\rm{e}}.{\rm{s}}.{\rm{u}}).$$
(54)

Therefore, the Lundquist number S = LVA/λ, where λ = c2/4πσ is magnetic diffusivity, is generally very large in solar and astrophysical plasmas:

$$S \approx {10^{ - 13}}\left({{T \over {{{10}^6}{{\rm{K}}^{ - 3/2}}}}} \right)\left({{L \over {{{10}^9}{\rm{cm}}}}} \right)\left({{{{V_A}} \over {{{10}^8}{\rm{cm}}\,{{\rm{s}}^{ - 1}}}}} \right).$$
(55)

Current and near future computational resources do not allow such an extreme large Lundquist number in numerical simulations. Therefore one should be always careful when interpreting the result of numerical simulation. This is general caution applicable for any numerical simulation, but it is particularly true for reconnection problems where the effect of diffusion (the change of field line connectivity) significantly affects the global topology and dynamics.

The effect of different resistivity models in MHD magnetic reconnection have been studied by many authors. As described in Section 2.2.3, one of the commonly used forms of resistivity (or magnetic diffusivity) is the so-called anomalous resistivity model (Ugai and Tsuda, 1977; Sato and Hayashi, 1979) where resistivity is given by an increasing function of the current density J or the ion-electron drift velocity vd = J/ρ, e.g.,

$$\eta = {\eta _0}{({v_d}/{v_{{\rm{cri}}}} - 1)^\alpha } + {\eta _c}\,{\rm{for}}\,{v_d} \geq {v_{{\rm{dri}}}},$$
(56)
$$ = {\eta _c}\,{\rm{for}}\,{v_d} < {v_{{\rm{cri}}}}.$$
(57)

Here vcri is the threshold above which the anomalous resistivity sets in, and the power index α is usually set to 1 or 2. The constant background resistivity ηc is often set to zero. A similar version uses current density J instead of the drift velocity.

These anomalous resistivity models can be considered heuristic subgrid models based on the theoretical expectation that the kinetic effects such as microscopic instability or wave-particle interaction arise when the current density or drift velocity become large, leading to the increase in effective resistivity. Recent laboratory experiments do show such anomalous increase in resistivity above the Spitzer value when the current sheet width becomes thinner than the ion inertia length (Ji et al., 1998). However, the physical mechanisms of the anomalous resistivity are not yet clarified, and the functional form and the threshold values used in the MHD models are rather ad hoc.

Ugai (2008) and Ugai (2012) present careful studies of the impact of different resistivity models on the properties of magnetic reconnection. The general conclusion is that, in the case of uniform resistivity, reconnection becomes Sweet-Parker type. Namely, a thin elongated current sheet forms and the reconnection rate is scale with the Lundquist number as Vinflow/VA = S−1/2 (Parker, 1957; Sweet, 1958; Biskamp, 1986). On the other hand, if the resistivity is localized due to the anomalous resistivity model, reconnection becomes Petschek type, namely the reconnection rate is almost independent to Lundquist number, Vinflow/VA = log S, and the magnetic energy is converted to the kinetic and thermal energy via standing slow-mode shocks (Petschek, 1964; Ugai, 1986; Scholer, 1989). See Priest and Forbes (2000) for a comprehensive review of the MHD theory of reconnection.

In the context of solar emerging flux, the effect of resistivity model on reconnection was studied by Yokoyama and Shibata (1994). Figure 45 shows the result of their 2D MHD simulation with uniform (left) and anomalous (right) resistivity models. The case with uniform resistivity model shows Sweet-Parker type elongated current sheet and the reconnection is slow. Whereas in the case with anomalous resistivity, Petschek-type fast reconnection with slow shocks develop, plasmoids (magnetic islands, see also Shibata et al., 1992b) are formed in the current sheet and are subsequently ejected.

Figure 45:
figure 45

Temporal evolution of current density (color) and magnetic field lines of 2D simulation of emerging flux with uniform (left) and anomalous (right) resistivity models. Image adapted from Yokoyama and Shibata (1994), courtesy of T. Yokoyama.

Moreover, high-Lundquist number simulations show that multiple plasmoids (magnetic islands) form in the reconnecting current sheet by the tearing instability. Coalescence and ejection of the plasmoids lead to intermittent and bursty reconnection. These processes are commonly observed in the MHD simulations with anomalous resistivity models (e.g., Shibata and Tanuma, 2001; Bárta et al., 2008) as well as in those with uniform but small (S > 104) resistivity (e.g., Loureiro et al., 2005; Ji and Daughton, 2011). Formation of multiple plasmoids and their ejection have been also found in the observation of coronal flares (Takasao et al., 2012) and chromospheric jets (Singh et al., 2012).

As discussed in detail in Section 3.8, the study of flux emergence through the upper photosphere and chromosphere may necessitate the use of a generalized Ohm’s law to capture the effect of neutral-ion interactions. If the Hall effect and ambipolar diffusion were important for the dynamics of flux emergence, a simple localized resistivity model based on anomalous resistivity would not suffice.

How much caution one should pay concerning the details of reconnection physics depends on the problems of interest. Although there are still many unsolved problems in the physics of reconnection, observations have confirmed that the reconnection rate is as large as 0.01–0.1 in the corona (Yokoyama et al., 2001; Isobe et al., 2002, 2005b; Noglik et al., 2005) though it is highly intermittent. It is therefore probable that the change of the global magnetic topology may not be so sensitive to the detail of reconnection process as long as it is sufficiently fast.

On the other hand, how much free energy one can store in an emerging flux region depends on the onset of fast reconnection. If the threshold of current density (or drift velocity vd) is higher, more free magnetic energy cannot be stored before the fast reconnection sets in and the stored energy is released. Yokoyama and Shibata (1994) found that a larger threshold of the anomalous resistivity results in more energetic reconnection.

The structure of the reconnection region also affects the re-distribution of the released energy. For example, the angle between reconnecting field lines affect the fraction of the released energy that goes to the Alfvén waves, which end up as enthalpy flux of the outflow (Takeuchi and Shibata, 2001). Furthermore, the formation of small scale plasmoids will cause the emission of high-frequency waves from the reconnection region (Karlicky et al., 2000, 2010; Isobe et al., 2008). The modes and spectra of the waves determine their dissipation mechanisms, hence the physics of reconnection is crucially important when we are interested in, e.g., the heating mechanism of the atmosphere or the acceleration mechanism of jets.

3.7.2 Field-aligned thermal conduction

The classical Spitzer thermal conductivity (Spitzer Jr, 1962) is given by

$$\kappa = {{1.84 \times {{10}^{ - 5}}{T^5}^{/2}} \over {\ln \Lambda }} \approx {10^{ - 6}}\,{\rm{erg}}\,{{\rm{s}}^{ - 1}}\,{{\rm{K}}^{ - 1}}\,{\rm{c}}{{\rm{m}}^{ - 1}},$$
(58)

where Λ is the Coulomb logarithm. In a fully ionized magnetized plasma such as the solar corona, thermal conduction becomes anisotropic so that the conductivity perpendicular to the magnetic field almost vanishes, whereas the conductivity parallel to the magnetic field remains unchanged.

Because of the nonlinear increase of the thermal conductivity with increasing temperature (κT5/2), thermal conduction is very effective in the high temperature (T > 106 K) corona, whereas it is almost negligible in the chromosphere and below.Footnote 12 Numerical simulations with the realistic thermal conduction in the corona is thus computationally demanding because it requires either an implicit numerical scheme or very small time step for explicit schemes. Many numerical simulation of flux emergence therefore neglect the thermal conduction in the corona. It may be justified, fortunately, if one is interested in the global dynamics of emerging flux. In the corona where the thermal conduction is efficient, the plasma-β is usually low and the dynamics is basically governed by the magnetic field. When the detailed thermodynamics is not important for the science question being address (coronal heating is not in this category), thermal conduction (and, to a certain extend, optically thin radiative losses) may be safely neglected.

There are some problems in which realistic treatment of thermal conduction is essential. One such example is chromospheric evaporation, in which the deposited energy in the corona is transported along the magnetic field to the chromosphere, increasing the gas pressure in the chromosphere and driving an upward evaporating (or ablating) flow. Chromospheric evaporation is established as the mechanism for the plasma supply to post-flare loops (Fisher et al., 1985; Yokoyama and Shibata, 2001) as well as the mechanisms for X-ray jets (Shimojo et al., 2001; Chifor et al., 2008). See Section 4.2 for the acceleration mechanism of jets.

Also, in order to make comparisons between observations and simulation results, one should be careful about how the energy equation is treated in the simulation. This is particularly important for forward modeling, in which the radiation is calculated from the simulation result for a direct comparison with observations. Thermal conduction is the most important factor that determines the temperature and density structure in the corona, whereas in the lower atmosphere radiative transfer is more important. Models that aim to treat the coupling between these various layers in a self-consistent manner must include both effects (e.g., one such code is Gudiksen et al., 2011).

3.8 Partial ionization effects in the lower solar atmosphere

3.8.1 Neutral-ion interaction effects in terms of a generalized Ohm’s law

As mentioned in Section 2.2, the MHD equations capture the principles of mass, momentum, and energy balance, as well as Faraday’s induction equation. These principles are general and apply to astrophysical systems over a wide range of scales. However, they must be supplemented by physical models of the material properties of the plasma. This includes statements about the equation of state, radiative properties of the plasma, as well as the appropriate Ohm’s law that relate the electric field to other MHD quantities.

In the majority of MHD simulations modeling magnetic flux emergence, the Ohm’s law used for the electric field is given by Eq. (12), which includes the contribution from ideal evolution (−v × B) and Ohmic diffusion (η∇ × B). Implicit in this assumption is that contributions due to neutral-ion interactions are negligible. While this may be valid in the fully-ionized corona, some discussion regarding its applicability in the predominantly neutral photosphere and chromosphere is in order.

In the discussion below, we restrict our attention to single-fluid MHD models, which means that the continuity, momentum, and energy equations are still solved for the combined neutral-ion fluids, but the induction equation is modified to take into account a generalized Ohm’s law for the electric field. When interactions between a neutral fluid and an ionized fluid are taken into account, it can be shown that the generalized Ohm’s law includes two additional terms, such that the induction equation becomes (Parker, 2007; Pandey and Wardle, 2008):

$${{\partial {\bf{B}}} \over {\partial t}} = \nabla \times \left[ {{\bf{v}} \times {\bf{B}}} \right] - {{4\pi \eta } \over c}{\bf{j}} - {1 \over {e{n_e}}}{\bf{j}} \times {\bf{B}} + ({\bf{j}} \times {\bf{B}}) \times {\bf{B}}{{{D^2}} \over {c\varrho {\nu _{{\rm{in}}}}}},$$
(59)

where j = c(4π)−1∇ × B is the current density, ne is the electron number density, ϱi is the mass density of the ionized component of the plasma, D = ϱn/(ϱi + ϱn) is the neutral mass fraction, and νin is the rate of ion-neutral collisions. The first and second terms inside the bracket on the r.h.s. of the equation are again the ideal induction and Ohmic induction terms. The third and fourth terms describe the Hall effect and ambipolar diffusion, respectively.

The Hall effect arises when electrons are attached to magnetic field lines while relatively immobile ions are not. This results in a Hall electric field which is proportional to j × B (i.e., the Lorentz force). Ambipolar diffusion, as it is often called in the astrophysical literature (e.g., Parker, 1963; Brandenburg and Zweibel, 1994), is due to collisions of neutrals and ions drifting with respect to each other. This drift results in an effective friction of the neutrals on the ions and leads to an electric field proportional to (j × B) × B.

3.8.2 Consequences of Hall and ambipolar effects

Before reviewing work that address the potential importance of the Hall and ambipolar effects in the context of flux emergence, it may be helpful to gain physical intuition about how the two effects impact the evolution of the magnetic field, assuming they cannot be neglected. The following discussion is summarized in Table 1. First of all, Eq. (59) can be written in the form

$${{\partial {\bf{B}}} \over {\partial t}} = \nabla \times \left[ {{\bf{u}} \times {\bf{B}} - {{4\pi } \over c}\eta {\bf{j}}} \right],$$
(60)

where

$${\bf{u}} = {\bf{v}} + {{\bf{v}}_{{\rm{Hall}}}} + {{\bf{v}}_{{\rm{Amb}}}},$$
(61)
$${{\bf{v}}_{{\rm{Hall}}}} = - H\nabla \times {\bf{B}},$$
(62)
$${{\bf{v}}_{{\rm{Amb}}}} = M(\nabla \times {\bf{B}}) \times {\bf{B}},$$
(63)
$$H = {c \over {4\pi e{n_e}}},$$
(64)
$$M = {{{D^2}} \over {4\pi {\varrho _i}{\nu _{{\rm{in}}}}}}.$$
(65)

As pointed out by Cheung and Cameron (2012), Eq. (60) is of the same form as Eq. (8), which is just the familiar induction equation (for fully-ionized plasmas) including −c−1v × B and Ohmic electric fields. Eq. (60) is different, however, in that the fluid velocity v has been replaced by a generalized velocity u, which is the sum of v and two ‘effective’ velocities vHall and vAmb capturing the contributions due to the Hall and ambipolar terms. The observation that the Hall and ambipolar terms acts like velocities operating on B has an interesting implication, namely that the two effects do not change the topology of the magnetic field. By topology, we mean topological measures such as magnetic helicity that are shown to be constants under ideal MHD evolution (Woltjer, 1958; Moffatt, 1969). Another way of saying this is to observe that both the Hall and ambipolar electric fields are perpendicular to B. Since the parallel electric fields for these effects are both zero, they do not contribute to magnetic reconnection. What this means is that the Ohmic term \((\;\;{{4\pi } \over c}\eta {\bf{j}}) = \;\;\;\;(\;\;\eta \;\;\;{\bf{B}})\) remains the only term that can directly allow for changes in magnetic connectivity (i.e., reconnection). In this restricted sense, the Hall and ambipolar effects are ideal.

Table 1 Summary of the impact of various effects present in MHD of partially ionized plasmas.

To examine the role of the Hall and ambipolar effects on the energy budget of the system, let us consider the evolution equation for B2/8π, which can be obtained by the dot product of Eq. (8) with B/4π, yielding

$${{\partial ({B^2}/8\pi)} \over {\partial t}} + \nabla \cdot \left[ {{c \over {4\pi }}{\bf{E}} \times {\bf{B}}} \right] = {\bf{E}} \cdot {\bf{j}}.$$
(66)

In terms of the Poynting flux \({c \over {4\pi }}{\bf{E}}\;\;\;{\bf{B}}\), both the Hall and ambipolar terms yield non-vanishing contributions (as do the ideal and Ohmic terms). In terms of E · j, the contribution due to the ideal −c−1v × B electric field becomes −FL· v, which is the rate of work done by the Lorentz force on the fluid. This term also appears (but with opposite sign) in the corresponding equation for the kinetic energy density \({1 \over 2}\varrho {v^2}\), and so it not a dissipative effect. The contribution of the Ohmic electric field to EOhm · j = −4πη/c2j2 = −j2/σc. This term is negative-definite and appears (with opposite sign) as a source term in the internal energy equation and specific entropy equation. So (as expected), Ohmic heating is a dissipative effect. The Hall electric field is non-dissipative since the Hall electric field is perpendicular to j. The ambipolar electric field, however, gives a contribution to \({{\bf{E}}_{{\rm{Amb}}}}\;\;{\bf{j}} = \;\;4\pi M{F_L}^2\), where M is the ambipolar coefficient defined in Eq. (65). Like Ohmic dissipation, ambipolar (or often called Pedersen) dissipation appears as a source term in the internal energy equation and specific entropy equation. In summary (see also Table 1), both Ohmic and ambipolar diffusion lead to plasma heating. In fact, since the ambipolar diffusion term is expected to be orders of magnitude larger than Ohmic diffusion under chromospheric conditions, it has been identified as a prime candidate for chromospheric heating (De Pontieu et al., 2001; Goodman, 2000; Khodachenko et al., 2004; Leake et al., 2005; Khomenko and Collados, 2011; Martínez-Sykora et al., 2012).

Although both the Hall and ambipolar effects do not directly change the topology of magnetic fields, their presence may have a major impact on the dynamics of the system, especially in the case of current sheets. First of all, since the magnetic field now evolves with a generalized velocity u instead of the plasma velocity v, the magnetic field is no longer necessarily frozen-in to the plasma even in the absence of Ohmic diffusion. As will be discussed in Section 3.8.3, this has potentially important consequences for the passage of magnetic flux through the solar atmosphere.

As demonstrated by Brandenburg and Zweibel (1994), ambipolar diffusion acts to sharpen current-carrying layers. In the absence of Ohmic diffusion to allow for flux cancellation (and reconnection), ambipolar diffusion would sharpen a current-carrying layer to a tangential discontinuity with singular current density. When Ohmic diffusion exists in tandem with ambipolar diffusion, the latter would sharpen current sheets to the point where Ohmic diffusion increases to balance the local grown in current density (see Cheung and Cameron, 2012, for a steady state solution).

The operation of the Hall electric field may possibly play a role in magnetic reconnection in the chromosphere, though this is a question that remains to be addressed. One of the important results of the Geospace Environmental Modeling (GEM) reconnection challenge (Birn et al., 2001), which compared numerical simulations of reconnecting Harris sheets with codes that ranged from full kinetic, hybrid to Hall-MHD and (non-Hall) MHD, concluded that fast, Petschek-type reconnection is a robust feature of any model as long as the Hall effect is included. Without the presence of the Hall effect, a MHD model with constant resistivity yielded only slow, Sweet-Parker type reconnection rates. However, it should be noted that the simulations carried out in the GEM challenge are for collisionless plasmas (e.g., applicable to the corona). Whether the same results carry over to collisional chromospheric plasmas remains to be tested. On the other hand, the presence of the Hall effect in partially ionized, collisional plasmas yields a dispersion relation which has the same character as that for the Hall effect in fully-ionized plasmas, namely a quadratic dependence of wave frequency w on wavenumber k. This quadratic dependence, which allows for whistler waves, has been invoked by the GEM collaboration as the key physical mechanism for facilitating fast reconnection. Finally, it should be noted that even without the Hall effect, MHD models with localized magnetic resistivity (e.g., anomalous resistivity, see Yokoyama and Shibata, 1994; Uzdensky, 2003) are able to produce fast reconnection.

3.8.3 Influence of ambipolar diffusion on the passage of magnetic flux through the chromosphere

The first simulations examining the potential importance of partial ionization effects on emerging flux was reported by Leake and Arber (2006). In this study, they considered how the use of a generalized Ohm’s law including ambipolar diffusion affected the buoyant rise of a twisted flux tube through an idealized, plane-parallel model of the solar atmosphere. In their paper, they follow Cowling (1957) and Braginskii (1965) and write the electric field (neglecting the Hall term) as

$${\bf{E}} = - {\bf{v}} \times {\bf{B}} + \eta {{\bf{j}}_\parallel } + {\eta _c}{{\bf{j}}_ \bot },$$
(67)

where η is the Ohmic resistivity and ηc is the so-called Cowling resistivity. j and j are components of j parallel and perpendicular to B, respectively. In this formulation, the two diffusivities are components of a resistivity tensor, which is anisotropic with respect to directions perpendicular and parallel to the field. It should be noted that this formation leads to an induction equation that is equivalent to Eq. (60), except the Hall term is assumed to vanish. In Cowling’s formulation, the absence of ambipolar diffusion would mean ηc = η.

Leake and Arber (2006) performed 2.5D simulations (with the tube axis pointing out of the plane) excluding and including the ambipolar effect and found that in latter case, much more magnetic flux was able to rise through chromospheric layers than in the former (see Figure 46). Furthermore, there was substantially less uplift of mass in the latter case. Both results can be explained as per the discussion in Section 3.8.2, where it was mentioned that in the presence of ambipolar diffusion (as well as the Hall term, though this effect is not treated in the current model), magnetic field is no longer frozen into the plasma. So when magnetic field lines are rising through the chromosphere, they do so with less of a mass burden. Furthermore, since ηcη in the chromosphere, there is preferential dissipation of the component of the current perpendicular to B. For this reason, they argue that in the simulation with Cowling resistivity, the emerged field is substantially closer to being force-free (i.e., j × B = 0) than in the simulation without Cowling resistivity.

Figure 46:
figure 46

The left and right panels respectively show field line contours from 2.5D flux emergence simulations through the chromosphere excluding and including Cowling resistivity. Image reproduced with permission from Leake and Arber (2006), copyright by ESO.

The work of Leake and Arber (2006) identified a potentially crucial effect that had not been treated (and remains neglected) in the majority of emerging flux models. Since the implications of this result can be far reaching for our understanding of solar flux emergence, it is worth considering in detail the physics that is treated. First of all, their model treats the effects of partial ionization in the sense that Cowling diffusivity is included as an additional term in the induction equation (i.e., the Pedersen current in a generalized Ohm’s law). The magnitude of this effect depends on the ionization fraction and collision rates between neutrals and ions, as indicated in the expression for the ambipolar coefficient M in Eq. (65). To treat the spatial dependence of this effect, the Cowling resistivity is computed using the local mass density g, temperature T, and field strength B. The energy equation includes a Newton cooling term given by

$${{d\varepsilon } \over {dt}} = - {{\varepsilon - {\varepsilon _0}(\varrho)} \over {\tau (\varrho)}},$$
(68)

where ε0 is a one-dimensional profile of a model of a background solar atmosphere and τ is an adjustment time scale. Both ϵ0 and τ were chosen to mimic the effects of shock heating, chromospheric and coronal/cooling bring about a persistent mean atmosphere. So when a parcel of plasma experiences a decrease in temperature and ϵ due to expansion, it would adjust to the background state (of comparable density) within a typical time scale τ. The equation of state used in the model is actually two relations. The one relating gas pressure to other state variables includes the effect of partial ionization, since the partial pressures of ions, electrons, and neutrals are included in the total pressure. The definition of internal energy density, however, does not include the component due to ionization potential energy. Leake and Linton (2013) addressed this issue by including ionization energy in the equation of state (cf. Section 2.2.1). The main conclusions, namely that Cowling resistivity allows for less mass uplift and evolution toward a more force-free field, remain unchanged. In 3D simulations with and without Cowling resistivity, Arber et al. (2007) reported that the distinctions revealed in 2.5D simulations were harder to ascertain. Robust results which carry over to 3D include the ability of the field to rise without significant mass uplift, the lack of exceedingly low temperature regions due to adiabatic expansion of the plasma and the nearly force-free nature of the emerged field. One big difference they found was that in the 3D case without Cowling diffusivity, the pile up of mass above the emerging magnetic structure lead to the development of magnetic buoyancy instabilities (see Section 3.3 for more on this topic). This complication made it difficult to isolate the true impact of Cowling diffusion.

By itself, ambipolar diffusion is a self-limiting process. This stems from the fact that ambipolar diffusion is dissipative and leads to a heating (in terms of both temperature and specific entropy) of the plasma. Associated with this heating is an increase in the ionization fraction and a decrease in the neutral-ion collision rate. Including this feedback on the thermodynamic state of the plasma leads to a quenching of ambipolar diffusion. An example of this effect was studied by Cheung and Cameron (2012) with simple, 1D numerical experiments of reconnection. In their experiments, heating by ambipolar diffusion in the current layer was very effective in switching off the effect, such that the amount of flux reconnected between simulations with and without ambipolar diffusion were negligible. However, when an energy sink (such as a Newton cooling term) was applied, ambipolar diffusion could proceed to increase the reconnection rate. The presence of such a cooling term in the energy equation of Leake and Arber (2006) and Arber et al. (2007) was likely crucial to keep ambipolar diffusion effective.

Martínez-Sykora et al. (2012) carried out 2D magnetoconvection simulations that include the effect of ambipolar diffusion. Their simulations have a non-LTE treatment of radiative transfer in the chromosphere, so that the effect of scattering of spectral lines is included and the source function both depends on the Planck function and ambient radiation field. This type of treatment may be one of the essential ingredients to quantifying the importance of ambipolar diffusion in the chromospheric layers. Their simulations were not specifically addressing the problem of flux emergence but the results are relevant for discussion. For instance, they found that the magnitude of the ambipolar coefficient can vary by many orders of magnitude a given height. This is due to the dynamical nature of the chromosphere and the steep dependence of the ambipolar coefficient on temperature. Based on this result they caution the use of 1D profiles of the ambipolar coefficient based on mean atmospheric models.

Finally, none of the previously mentioned work take into account the non-LTE nature of ionizing plasma in the chromosphere. Since the recombination time of H at chromospheric conditions (102–103 s, see Kneer, 1980; Leenaarts et al., 2007) is comparable to or longer than the characteristic time between successive passage of acoustic shocks launched from below, the ionization fraction of hydrogen is not in LTE, nor in statistical equilibrium (the latter assumed by Leake and Arber, 2006; Arber et al., 2007). So fully dynamic radiative MHD simulations with non-LTE radiative transfer and time-dependent hydrogen ionization are needed to quantitatively examine the impact of ambipolar diffusion in affecting chromospheric dynamics.

4 Flux Emergence and its Relation to Jets and Eruptions

Emerging flux is the ultimate origin of solar active regions. Therefore, it is no surprise that emerging flux plays a prominent role in many models of solar flares and eruptive events. Most work in the literature addressing the connection between emerging flux and eruptions can be separated into two categories, those in which the emerging flux region is itself a large part of the erupting structure, and those in which emerging flux acts only as a trigger for the eruption of pre-existing magnetic field. The rest of this section focuses only on work that is directly pertinent to flux emergence. For more comprehensive reviews of modeling of flares and eruptions, see Forbes (2000), Klimchuk (2001), Shibata and Magara (2011), and Chen (2011).

4.1 Interaction with pre-existing field; overall picture

Heyvaerts et al. (1977) proposed the well-known model of a solar flare driven by magnetic reconnection of an emerging flux and the pre-existing field. Figure 47 shows a schematic illustration from their study. Since then there have been numerous papers on magnetic reconnection between emerging flux and the ambient field, but the essential picture, namely current sheet formation and reconnection between the emerging flux and the pre-existing field, remains unchanged.

Figure 47:
figure 47

A model of solar flares produced by magnetic reconnection between the emerging flux and the pre-existing field. Image reproduced with permission from Heyvaerts et al. (1977), copyright by AAS.

The first MHD simulation of magnetic reconnection between the emerging flux and pre-existing field was done by Forbes and Priest (1984). In this paper the emerging flux was driven in from the bottom boundary of the computational domain, which initially consisted of uniform, horizontal magnetic field with opposite direction to the emerging flux. They found that, because the emergence is slower than the Alfvén time-scale τA but faster than the reconnection time-scale given by the Sweet-Parker reconnection model τS1/2τA, where S is the Lundquist number, a region of closed loops and a current sheet at the top of the loop are formed. In the later phase reconnection became faster and super-magnetosonic flows (jets) were produced. Thus, the model by Heyvaerts et al. (1977) shown in Figure 47 was demonstrated by MHD simulation.

Later, Shibata et al. (1992b) performed 2D MHD simulations of flux emerging into a stratified atmosphere with a pre-existing horizontal field in the model corona. They self-consistently solved the Parker instability of a flux sheet and reproduce similar morphology and dynamics, namely the formation of Ω-shaped loop, its reconnection with pre-existing field and ejection of jets. They also found that multiple plasmoids were formed in the reconnecting current sheet and ejected intermittently. However, due to the limited resolution the reconnection was likely to be numerical.

Interaction of an emerging twisted flux tube with the pre-existing horizontal coronal field was first studied in 3D by Archontis et al. (2004). In Figure 48 (adopted from their paper), one can recognize the evidence of reconnection between the flux tube and the coronal field in the evolution of the field line connectivities. This probably corresponds to the observation of the formation of coronal loops that connect the newly emerging region and pre-existing flux system (i.e., neighbouring active regions) reported by Longcope et al. (2005).

Figure 48:
figure 48

Time evolution of the magnetic field lines in the simulation by Archontis et al. (2004). The grey isosurfaces show the location of strong fields corresponding to the part of flux tube that remains below the photosphere. The red field lines are traced from the isosurfaces, while the other field lines are traced from the side boundaries. Image reproduced with permission, copyright by AAS.

Following Archontis et al. (2004), the interaction of the emerging twisted flux tube and the horizontal coronal field has been studied in detail (Archontis et al., 2005; Galsgaard et al., 2005; Archontis et al., 2006; Galsgaard et al., 2007). It was found that when the coronal field and the two flux systems are nearly anti-parallel, reconnection between them produces high-temperature and high-velocity outflows. However, when the two flux system have more of a parallel orientation, the amount of reconnected flux is limited and no such energetic features appear. Interestingly, despite having different amounts of reconnected flux in the cases with different orientations, the height-time profiles of the apex of the emerging loops look almost identical (Galsgaard et al., 2007) in different cases. The plasmoids (magnetic islands) in the reconnection region found in 2D simulations (Shibata et al., 1992b; Yokoyama and Shibata, 1994) are also found in 3D simulations by Archontis et al. (2006). Due to the three-dimensionality the plasmoids appear as twisted flux tubes.

2D simulations of flux emerging into pre-existing oblique fieldFootnote 13 was first studied by Yokoyama and Shibata (1995, 1996), in which they found both hot and cool jets along the oblique field lines as shown in Figure 49. The hot and cool jets are adjacent to each other, which is consistent with simultaneous observations of X-ray jets and Hα surges (Canfield et al., 1996).

Figure 49:
figure 49

Temporal evolution of temperature (color), velocity (arrows), and magnetic field lines of 2D simulation of emerging flux with pre-existing oblique fields. Image adapted from Yokoyama and Shibata (1996), courtesy of T. Yokoyama.

3D simulations with an emerging twisted flux tube into oblique ambient field was performed by Moreno-Insertis et al. (2008). Figure 50 shows the 3D topology of the magnetic field. The orange, blue and green lines indicate the field lines in the emerging flux, the pre-existing field, and the reconnected field, respectively. The temperature (red, T = 6.5 MK) and J/B (blue) isosurfaces are also shown to visualize the structure of the jet and the reconnecting current sheet. Many results of 2D simulations (Yokoyama and Shibata, 1995, 1996) are reproduced in 3D simulations. An extension of the work by Moreno-Insertis et al. (2008) was undertaken by Moreno-Insertis and Galsgaard (2013), in which they carried out the flux emergence simulation in a larger domain and for longer duration. They found that, after the main jet phase, a number of eruptive events occurred from the same emerging flux region. This will be discussed further in Section 4.2.

Figure 50:
figure 50

3D visualization of the MHD simulation of an emerging twisted flux tube and pre-existing oblique field by Moreno-Insertis et al. (2008). The orange, blue, and green lines indicate the field lines in the emerging flux, the pre-existing field, and the reconnected field, respectively. The red and blue isosurfaces are that of temperature (T = 6.5 MK) and J/B (blue), respectively. Image reproduced with permission, copyright by AAS.

The overall picture is thus established, namely that emergence flux into a pre-existing flux system causes reconnection between them, produces plasma heating and dynamic events such as jets. However, it seems there are variety of the acceleration mechanisms of jets. In the next section we discuss the different mechanisms of jet acceleration. Moreover, in cases where the pre-existing flux system has more complicated structure, the emerging flux may trigger the large scale eruption of the pre-existing flux system. In some cases the significant fraction of the emerging flux itself erupts. This topic is discussed in Section 4.

Some observations show that jets recur from the same region (Chae et al., 1999; Chifor et al., 2008). The recurrent behavior of reconnection and jets in the MHD simulation of flux emergence was first reported by Murray et al. (2009). As shown in Figure 51, reconnection reversal (or oscillatory reconnection) takes place; reconnection occurs in distinct bursts and the inflow and outflow of one burst become the outflow and inflow in the following burst of reconnection, respectively, and this process repeats. The reversal of reconnection occurs because the reconnection outflow cannot escape from the reconnection site and, hence, the gas pressure gradient in the outflow region increases and eventually reverses the direction of the flows. This phenomena has been further investigated in detail by McLaughlin et al. (2012), who confirmed that the mechanism of the oscillatory reconnection is the local imbalance of forces, primarily the gas pressure gradient, between the neighboring flux systems. Archontis et al. (2010) found that essentially the same process also occurs in 3D simulation with a twisted emerging flux tube. Both in 2D and in 3D, the recurrent, oscillatory reconnection gradually becomes weaker and eventually ceases when the system settle down to a new equilibrium. It may also explain observations of the quasi-periodic propagating transverse waves (Liu et al., 2011b).

Figure 51:
figure 51

Gas pressure (a–c) and magnetic pressure (d–f) and magnetic field lines in 2D MHD simulation of emerging flux and oscillatory reconnection. Image reproduced with permission from Murray et al. (2009), copyright by ESO.

In the case where the sub-photospheric flux system is much longer than the characteristic scale of the undular mode (scale heights ∼ 3000 km), multiple Ω-loops can emerge. Figure 52 shows the result of a 2D MHD simulation of a such case by Isobe et al. (2007b). The neighbouring loops reconnect and cause heating in the lower atmosphere that may be observed as Ellerman bombs. The reconnection events also remove mass from the field lines and allow the emergence of long field lines. A 3D version of this simulation was presented by Archontis and Hood (2009). This process is also caused when realistic granular convective flows are present in the models (see Section 3.4.3).

Figure 52:
figure 52

2D simulation of hierarchical evolution of multiple emerging loops. Image reproduced with permission from Isobe et al. (2007b), copyright by AAS.

4.2 Acceleration mechanisms of jets

Jets, or collimated flows from emerging flux regions have been observed at various wavelengths. Hα observation of the chromosphere often show dark jets, also called surges (Roy, 1973; Brooks et al., 2007). They appear as dark absorbing feature in Hα and represent chromospheric (∼ 104 K) material.

The soft-X-ray observations by Yohkoh revealed numerous X-ray jets (Shibata et al., 1992a; Shimojo et al., 1996), with many of them associated with emerging flux or cancellation events (Shimojo et al., 1998). The typical physical parameters of the jets are: temperature: 3–8 MK, density: 0.7–4.0 × 109 cm−3, and the apparent velocity: 180–530 km s−1, comparable or smaller than the sound velocity. In addition to X-ray observations there are also many imaging and spectroscopic EUV observations of coronal jets (e.g., Chae et al., 1999; Chifor et al., 2008).

Higher resolution observations by Hinode X-ray telescope subsequently revealed two distinct components: a faint and faster component (velocity ∼ VA ∼ 800 km s−1), and a bright and slower one similar to those found in Yohkoh observations (Cirtain et al., 2007). Transverse motion in the jets was also found, suggesting the presence of Alfvén or fast kink-mode waves (Cirtain et al., 2007; Savcheva et al., 2007). The transverse motion has been also found in chromospheric jets observed in spicules and chromospheric jets (De Pontieu et al., 2007; Nishizuka et al., 2008; Liu et al., 2011a; Okamoto and De Pontieu, 2011). High-resolution observations from the Solar Optical Telescope onboard Hinode led to the discovery of numerous small-scale jets in the chromosphere (Shibata et al., 2007; Singh et al., 2012) whose morphology is similar to X-ray and EUV jets (see Figure 53).

Figure 53:
figure 53

Left: Schematic illustration of jets associated with emerging flux observed in different scales. Image reproduced with permission from Shibata et al. (2007), copyright by AAAS. Right: example of observed jets in different scales. From top to bottom, an X-ray jet observed by the X-ray Telescope of Hinode, an EUV jets observed in Fe XII 195 Å filter by TRACE (Nishizuka et al., 2008), and a chromospheric jet observed by the Solar Optical Telescope of Hinode (Singh et al., 2012).

Figure 54:
figure 54

Left: 2D simulations of emerging flux with pre-existing horizontal coronal field including the anisotropic thermal conduction. Top and middle panels show density distribution (color), magnetic field lines, and velocity field. The color in the bottom panes shows temperature distribution. Two evaporation jets are seen in the both side of the emerging loop as dense ejecting structure with coronal temperature. Right: same simulation but without thermal conduction. Image reproduced with permission from Miyagoshi and Yokoyama (2003), copyright by AAS.

It is tempting to invoke reconnection outflow as the cause of jets but the answer is not so simple. First of all, 1D numerical simulations and EUV spectroscopic observations indicate that many, though not all, of the soft X-ray jets are not the reconnection outflow itself, but are so-called evaporation jets driven by the gas pressure gradient (Shimojo et al., 2001; Chifor et al., 2008).

Chromospheric evaporation is a process in which the magnetic energy released via magnetic reconnection in the corona is transported to the chromosphere preferentially along the magnetic field either by anisotropic thermal conduction or high-energy particles. Therefore, in order to study the effect of chromospheric evaporation in MHD simulations coronal thermal conduction is necessary (see discussion in Section 3.7.2).

Miyagoshi and Yokoyama (2003) and Miyagoshi and Yokoyama (2004) performed 2D simulations of emerging flux with pre-existing horizontal coronal field including the anisotropic thermal conduction. Their comparison of simulations with and without thermal conduction demonstrates that thermal conduction is necessary to reproduce high-density evaporation jets. Hot and Alfvénic reconnection outflow was also found in the simulations. The reconnection outflow and the evaporation jet probably correspond to the faint and bright components found by the Hinode/XRT observations (Cirtain et al., 2007).

The acceleration of cool jets is more problematic. It is often naively interpreted that if reconnection occurs in the chromosphere it can launch a jet with chromospheric temperature. Consider the extreme case when all the magnetic energy is used to lift the plasma in the same volume against gravity. The maximum height h is given by

$$h = {{{B^2}} \over {8\pi \rho g}}$$
(69)
$$ = {H \over \beta},$$
(70)

where H = kBT/mg is the scale height and β = 8πρkBT/mB2 is the ratio of gas to magnetic pressures. Considering the typical scale height H ∼ 300 km and plasma beta β ∼ 0.1–10 in the chromosphere, the height of the chromospheric jets would be 30–3000 km, much shorter than observed height of the chromospheric jets (Singh et al., 2012). This is particularly true when reconnection takes place in the lower chromosphere where β is close to unity except in the sunspots. In such a case one needs a mechanism by which the released magnetic energy in the reconnection region is transported upward and subsequently accelerate the plasma where the density is much lower.

Wave propagation and steepening is one such mechanism. Reconnection in the lower atmosphere can generation slow mode wave (Takeuchi and Shibata, 2001) and it is known that as a slow-mode wave propagates upward, its amplitude grows due to the the density contrast, eventually steepening into a shock. The interaction of the shock and the transition region between the chromosphere and the corona can launch the transition region upward, which is observed as a chromospheric jet (Osterbrock, 1961; Shibata and Suematsu, 1982; De Pontieu et al., 2004; Heggland et al., 2009). From their emerging flux simulations with the Bifrost code (which includes radiative transfer, thermal conduction, and changes in ionization state), Martínez-Sykora et al. (2011) reported the occurence of a cool chromospheric jet resulting from the interaction of emerging flux with the overlying arcade field. They identified that reconnection between the emerging flux and ambient field likely occurred in the high-β photosphere. As the reconnected field emerged further, however, the Lorentz force accelerated plasma horizontally. Pressed against a neighbouring wall of stronger (i.e., low-β), predominantly vertical field, the squeezing of the plasma led to a local enhancement of the gas pressure, which accelerated material upwards along the magnetic field lines. While this may be sufficient to account for the penetration of the chromospheric material into the transition region (≈ 4 Mm above the normal height, see also Eq. (70)), the slow-mode wave mechanism likely plays a role in the further evolution of the jet. Takasao et al. (2013) presented a detailed analyses of the 2D simulation result with similar setup with Yokoyama and Shibata (1996) and found there are three types of acceleration mechanisms of cool jets, as summarized in Figure 55.

Figure 55:
figure 55

Schematic illustration of the acceleration mechanisms of chromospheric jets. Image reproduced with permission from Takasao et al. (2013), copyright by PASJ.

Finally, some recent observations suggest there are other types of coronal jets called blowout jets (Moore et al., 2010). This type of jet is thought to be a miniature version of large-scale eruptive events such as filament eruptions and CMEs (see Sakajiri et al., 2004, for similar observations in Hα). The difference from “standard” jets (e.g., Figure 53) seem to be following: while the standard jets can be more or less regarded as a consequence of direct reconnection between emerging flux and ambient field, the blowout jets are more likely to be a store-and-release process. That is, the free energy and twist is at first stored in a quasi-equilibrium magnetic configuration in the atmosphere, which becomes unstable or loses equilibrium and erupts in a way similar to flares and CMEs. See Shibata and Magara (2011) and Chen (2011) for a review of flare/CME models, and see Archontis and Hood (2013) and Moreno-Insertis and Galsgaard (2013) for 3D simulations of emerging flux that lead to blowout-type jets. Indeed this type of model may explain the jets that shows significant twisting (or untwisting) motions (Kurokawa et al., 1987; Patsourakos et al., 2008) as demonstrated by 3D simulation by Pariat et al. (2010).

4.3 Emerging flux as a trigger for eruptions

First of all it should be noted that emerging flux can trigger at least small flares (plasma heating) and associated jets, as proposed by Heyvaerts et al. (1977) and already discussed extensively in Section 4.1. This kind of model can be classified as directly driven models, whereas those models in which the emerging flux trigger the eruption of pre-existing flux system larger than the emerging flux itself can be considered as store and release models (Klimchuk, 2001). In this section we cover only the latter.

Observationally, Feynman and Martin (1995) established the picture of the eruption trigger by emerging flux by studying the 53 quiescent filaments. They found that when the emerging flux is oriented favorably for magnetic reconnection with the magnetic field of the filament channel, it is more likely to trigger the eruption (see also Wang and Sheeley Jr, 1999)).

Chen and Shibata (2000) first demonstrated by 2.5D MHD simulations that emerging flux favorable for reconnection could trigger the eruption of a pre-existing flux rope in the corona, while those not favorable for reconnection could not. Figure 56 schematically illustrates the trigger mechanism of the flux rope in two different cases. In one case (panel a), the emerging flux appears just below the flux rope, which then undergoes reconnection with the pre-existing small loop. The reconnection breaks the mechanical equilibrium of the flux system by reducing the magnetic pressure below the flux rope. The reduced pressure induces the inflows from both sides that carry the magnetic field overlying the flux rope, and reconnection of these field lines brings the flux rope further from mechanical equilibrium. This process is similar to the so-called runaway tether-cutting model (Moore et al., 2001). In the other case (panel b), the emerging flux appears in slightly shifted location and, hence, it undergoes reconnection with the overlying field lines that stabilizes the eruption of the flux rope. As one can see from the figure reconnection weakens the magnetic tension of the overlying field, thus breaking the equilibrium and leading to the eruption of the flux rope. This case is a bit similar to so-called Breakout model (Antiochos et al., 1999) in the sense that the eruption is triggered by the weakening of the tension of the overlying bootstrapping field.

Figure 56:
figure 56

Schematic illustration of the trigger of flux rope eruption by reconnection-favorable emerging flux. Image reproduced with permission from Chen (2008), copyright by the Indian Academy of Sciences.

Lin et al. (2001) studied the similar problem with semi-analytic approach and found that the emerging flux does trigger eruption in some cases but the conditions for the eruption trigger is rather complicated and depends also on parameters other than the field orientation, such as the field strength, distance, and area of the emerging flux region. This conclusion is in agreement with the observational study by Zhang et al. (2008). Dubey et al. (2006) have also carried out a parameter study by 2.5D MHD simulations including the effects of gravity and spherical geometry.

Zuccarello et al. (2008) and MacTaggart and Hood (2009a) performed 2.5D simulations of flux emergence into a quadrapole configuration to study whether the emerging flux triggers breakout-type eruption. Despite differences in the initial setup (see the last part of this subsection) the two simulation studies show good agreement: emerging flux did trigger an eruption but it did not lead to the formation of a flux rope which pinches off the central arcade. Comparison of emerging flux and the shearing motion as the trigger for the breakout-type eruption was carried out by Zuccarello et al. (2009). 2.5D simulations by Leake et al. (2010) also show the formation of multiple flux ropes due to the reconnection of the flux rope and the pre-existing dipole or quadrapole coronal field. However, much of the magnetic shear and its free energy remains in the lower atmosphere and no large-scale eruption was observed.

In order for an emerging flux to trigger the eruption the pre-existing flux system must have some free energy. MacTaggart (2011) modeled by 3D simulation the scenario of a flux tube emerging into a potential field arcade, which as expected did not result in eruptions. In order to study the condition for eruption, Kusano et al. (2012) performed a systematic parameter study of 3D simulations. The models initially have a force-free sheared arcade (as shown in panel (a) of Figure 57) and in each case an emerging bipole magnetic flux inserted from the bottom boundary on the neutral line. They systematically changed the amount of the shear (angle θ0 between the polarity inversion line and the magnetic field) and angle φe of the emerging bipole to study which combinations triggered eruptions. Figure 57 shows one example of their simulations that successfully produced eruption.

Figure 57:
figure 57

Three-dimensional visualization of the case of a successful eruption by Kusano et al. (2012). Each subset represent a birds eye view (a, c, e–h), top view (b), and enlarged side view (d) of the magnetic field at different times. Green tubes represent magnetic field lines with connectivity that differs from the initial state. Selected magnetic fields in the initial state and those retaining the initial connectivity are plotted by blue tubes in (a) and (d), respectively. Red isosurfaces correspond to intensive current layers, and ray scales (white, position) on the bottom plane indicate the distribution of the z component of magnetic field Bz. Image reproduced with permission, copyright by AAS.

Figure 58 shows a diagram that summarize the simulation results. It is found that when φe ∼ 180 degree, the inserted emerging flux has the reconnection-favorable orientation, and in such case the eruption occurs almost independently from the initial condition. On the other hand, even if the shear is strong and a lot of free magnetic energy is stored (θ > 60 degree), the non-reconnection-favorable emerging flux (φe < 100 degree) cannot not trigger the eruption. This diagram is very interesting, though it is not yet clear how universal this behavior. Similar systematic studies should be carried out for different magnetic topologies, in particular with a pre-existing coronal flux rope.

Figure 58:
figure 58

Diagram showing the parameter study of Kusano et al. (2012) with various shear angles θ0 and the emerging flux orientation φe. The orientations of the coronal and emerging fields for different θ0 and φe are also shown. Image reproduced with permission, copyright by AAS.

Figure 59:
figure 59

Snapshots of the field lines of the 3D simulation by Fan and Gibson (2004) as viewed from the side (left panels) and the top (right panels). Image reproduced with permission, copyright by AAS.

Finally, it should be noted that different numerical treatment of emerging flux are adopted in the literature discussed in this section. Many of them use boundary conditions to insert the emerging flux in the domain (Chen and Shibata, 2000; Dubey et al., 2006; Zuccarello et al., 2008, 2009; Kusano et al., 2012), whereas others solve the buoyant rise of the emerging flux from below the photosphere self-consistently (MacTaggart and Hood, 2009a; MacTaggart, 2011; Leake et al., 2010).

4.4 Eruptions of emerging flux ropes

Apart from triggering eruptions, flux emergence can generate self-contained flux ropes which eventually erupt. This approach is more related to the energy storage processes prior to the eruption. The question of whether eruptive flux ropes are formed just prior to eruption or are transported by flux emergence is a question that has attracted considerable attention (Green and Kliem, 2009; Patsourakos et al., 2013).

Even though flux emergence is associated with energy injection, it does not necessary lead to a sufficient build up of free energy nor head down an evolutionary path that leads to an eruptive release of the available free energy. Certain ingredients are needed to induce such behavior. Fan and Gibson (2003, 2004) presented 3D simulations of an kink-unstable flux tube inserted from the lower boundary into the pre-existing potential arcade. It was found that the kink instability developed when a sufficient amount of twist is transported in the corona, and consequently S-shaped field lines and current concentration formed.

S-shaped structures in the corona observed in X-ray are called sigmoids and are known to be a good precursor of eruptions (Rust and Kumar, 1996; Canfield et al., 1999). Therefore, the formation of the sigmoidal structures as a result of the emergence of a twisted flux tube have been studied by many authors (Matsumoto et al., 1998; Magara and Longcope, 2001; Fan, 2009a; Archontis et al., 2009; Archontis and Hood, 2012).

Manchester IV et al. (2004) first reported the self-consistent 3D simulation of the emergence of a twisted flux tube, which drove shear flows that eventually lead to the formation of a detached coronal flux rope which then ejected from the modeled AR (see Figure 42). The difference between Fan and Gibson (2003, 2004) and Manchester IV et al. (2004) is that the scenario studied by the former authors involves the bodily rise of a coherent flux rope (including axial fields and dipped fields beneath the axis) through the photosphere (their bottom boundary) whereas the latter does involve the emergence of the axial fields (nor dipped field lines beneath the axis). The flux rope naturally formed in the latter case from shear flows. The formation of coronal flux ropes by reconnection has been also reported in other simulations of flux emergence (Magara, 2006; Archontis and Török, 2008; Archontis and Hood, 2010; Leake et al., 2013; MacTaggart and Haynes, 2014). Magara (2007) presented a 3D simulation in which three portions of the initial flux tube emerge at once. He found that the resultant magnetic structure is reminiscent to that of models of filaments.

Archontis and Török (2008) reported that if a magnetic flux tube emerges into a non-magnetized atmosphere, the induced shearing motion and reconnection can form a flux rope in the corona. However, the expanding magnetic field surrounding the flux rope acts as a bootstrapping field that suppresses an eruption. In contrast, if the expanding tube is allowed to reconnect with a pre-existing coronal field, the flux rope is able to erupt with an acceleration profile that is similar to those of observed filament eruptions and CMEs (see also Leake et al., 2014). Archontis and Hood (2012) further investigated the condition for eruptions by systematically changing the orientation of the ambient field. As shown in Figure 60, they report that favorable orientations for reconnection increase the likelihood of eruption.

Figure 60:
figure 60

Temporal profile of the height of the emerging (erupting) flux and its vertical velocity for different orientations of the ambient field. Φ = 0, 90, 180 corresponds to the ambient field parallel, perpendicular, and anti-parallel to the axis of the initial flux tube. Image reproduced with permission from Archontis and Hood (2012), copyright by ESO.

Whether the new coronal flux system formed by emerging flux erupts or not depends on the various parameters of the emerging flux itself as well as the interaction of the emerging flux and the pre-existing magnetic fields. Theoretically, we expect that another important factor that determines the stability of the flux system is the upward gradient of the magnetic field (Kliem and Török, 2006; An and Magara, 2013).

Finally, as briefly mentioned in Section 4.2, Moreno-Insertis and Galsgaard (2013) found in the 3D simulation of a twisted emerging flux with an ambient oblique field that multiple eruptions that look like mini-CMEs from the single emerging flux region. They discuss it as a possible mechanism for so-called blowout jets (Moore et al., 2010). From a simulation of the emergence of a highly twisted flux rope across the lower boundary into a pre-existing coronal arcade field, Chatterjee and Fan (2013) reported the occurence of homologous eruptions. These eruptions lead to a series of three CMEs, one of which catches up with a preceding CME and the two merge in a so-called ‘cannibalistic’ fashion. While each of these two studies is interesting in their own right, together they provide insight regarding the physics that drive eruptions at various scales.

5 Data-driven and Data-inspired Models of Emerging Flux

With the ever growing body of work on flux emergence, the field has matured to a point where lessons drawn from observations, theory and numerical experiments are beginning to be applied to modeling of actual emerging flux regions on the Sun. Some of the numerical model efforts in this line of research use observational data to drive their simulations, while others take inspiration from the observations to construct idealized scenarios.

5.1 Flux emergence models with one-way coupling

The first data-driven models of flux emergence did not use observational data. Abbett and Fisher (2003) carried out anelastic simulations of the rise of twisted Ω-loops in the convection zone. Due to limitations of the analestic model, the computational domain of these simulations has a top boundary in the subphotospheric layers. The MHD variables from horizontal planes of the self-contained analestic simulations were then used to set the ghost cell values of a 3D compressible MHD code, which has a computational domain that captures the photosphere and corona (see Figure 61). This way the compressible MHD simulations of flux emergence are driven by evolution in the anelastic simulations. As noted in their paper, an assumption made is that the MHD quantities do not vary significantly in the near-surface layers of the convection zone (i.e., the region skipped between the anelastic and compressible models). Fully compressible simulations of flux emergence in granular convection have since shown otherwise (e.g., Cheung et al., 2007a; Abbett, 2007; Martínez-Sykora et al., 2008, 2009; Tortosa-Andreu and Moreno-Insertis, 2009; Fang et al., 2010, 2012). Nevertheless, the study by Abbett and Fisher (2003) showed that data-driven simulations using different codes to capture different subdomains and physical regimes is possible.

Figure 61:
figure 61

3D rendering illustrating the one-way coupled model of flux emergence by Abbett and Fisher (2003). MHD quantities from anelastic simulations of the rise of Ω-loops in the convection zone are used to drive the bottom boundary of a compressible MHD code. Image reproduced with permission, copyright by AAS.

A more recent example of a model of flux emergence with one-way coupling is the work by Chen et al. (2014). As bottom boundary data for their model of the formation of AR coronal loops, they used MHD quantities sampled at the base of the photospheric layer in a fully-compressible, radiative MHD simulation of AR formation (Rempel and Cheung, 2014, see Sections 3.4.6 and 3.5). The numerical model that provided the bottom boundary data captures the top 16 Mm of the convection zone but only the first few hundred km of the solar atmosphere. In contrast, the computational domain of the Chen et al. (2014) model has a bottom boundary at the photosphere and extends to a height of 73 Mm. In this data-driven simulation, they were able to model the formation of million degree coronal loops above the model AR (see Figure 62).

Figure 62:
figure 62

This model of the formation of AR corona by Chen et al. (2014) is driven at the bottom boundary (z = 0) by a radiative MHD simulation of AR-scale flux emergence (Rempel and Cheung, 2014). The greyscale in panels (a) and (b) shows synthetic magnetograms at the photospheric layer. The yellow shading in all four panels shows EUV images of AR coronal loops synthesized for the 193 A channel of SDO/AIA. Image reproduced with permission from Chen et al. (2014), copyright by ESO.

5.2 Formation of current layers in emerging flux regions

Motivated by the discovery of serpentine field lines, bald patches and Ellerman bombs in emerging flux regions (Pariat et al., 2004, see also Section 3.4.3 in this paper), Pariat et al. (2009) developed an MHD model for driving current build-up in an emerging flux region. To do so, they constructed a potential field extrapolation from a SoHO/MDI magnetogram of an emerging flux region. Using this potential field as an initial condition for an MHD calculation, the introduced Maxwell stresses in the field by applying a expansion flow transverse to the bottom boundary. They did so while imposed a line-tied bottom boundary condition for the magnetic field. While this boundary condition does not truly mimic an emerging flux scenario, it served as a instructive experiment to test where currents would form due to photospheric driving. They reported that normalized current densities (j/B) were found to be especially intense at bald patches, where U-loops are expected to be pinched off following magnetic reconnection and where Ellerman bombs are expected to occur. They note, however, that currents also develop at different locations along separatrix surfaces and conclude that reconnection is not confined to the vicinity of bald patches.

5.3 Active region eruptions following flux emergence

Inspired by the evolution of NOAA AR 10930, which erupted on December 13, 2006 with an X3.4 flare and an Earth-directed CME, Fan (2011) carried out idealized MHD simulations in a spherical subdomain. From observations of the AR 10930, it can be inferred that the formation of the delta spot associated with the eruption was formed following the emergence of a parasitic flux rope into a pre-existing spot. In order to restrict the time-step of the MHD simulation to reasonable values, the initial condition was constructed by smooth an MDI magnetogram so that the maximum field strengths decreased from 3 kG to 200 G (to decrease the Alfvén speed, which limits the simulation time-step). A twisted toroidal flux rope with an east-west orientation was emerged just south of the pre-existing sunspot in such a way that the following polarity of the new region reconnected with the pre-existing spot. This interaction lead to some erosion of the bootstrapping arcade field overlying the new flux rope, which became unstable erupted (see Figure 63). Qualitative comparison with observations (such as the locations of the double ribbon structure) suggests this is a plausible scenario for AR 10930. Since the configuration of the magnetic field was close to what is needed to instigate both the torus and the helical kink instabilities, the ultimate cause for the eruption in the simulation could not be attributed to solely one or the other.

Figure 63:
figure 63

mpg-Movie (8155.41503906 kB) Still from a movie showing 3D rendering from a simulation of an eruption following emergence of a twisted flux tube into a pre-existing sunspot. This is an idealized scenario to mimic the evolution of NOAA AR 10930, which produced an X3.4 flare with an Earth-directed CME. Observations of AR 10930 showed the formation of a delta spot following emergence of a parasitic bipolar region. Reproduced with permission from Fan (2011), copyright by AAS. (For video see appendix)

To model the evolution of an eruptive active region over multiple days, Cheung and DeRosa (2012) performed data-driven simulations using the magnetofriction (hereafter MF) method (Yang et al., 1986; Craig and Sneyd, 1986). This method assumes that plasma velocity v in the induction equation to be proportional to local value of the Lorentz force.Footnote 14 This allows the magnetic field to evolve to a lower energy state while keeping field line connectivity (except at current layers where magnetic diffusion occurs). Cheung and DeRosa (2012) used the MF code to model the evolution of NOAA AR 11158. During the course of its birth (into a relatively quiescent region on the Sun) and passage across the solar disk, AR 11158 launched multiple eruptions and large M- and X-class flares (see Schrijver et al., 2011; Sun et al., 2012).

For their data-driven model, Cheung and DeRosa (2012) used an potential field extrapolation as an initial condition, and imposed Ex and Ey at the bottom boundary to drive the model forward in time. Since they only used longitudinal magnetograms (albeit remapped from SDO/HMI data), the horizontal electric field components are not well-constrained. So they performed numerical experiments under varying assumptions. In particular, they tested how the presence (or absence) of sustained twisting/shearing of the photospheric field affected the structure of the modeled AR. They found that in the case of applied twisting/shearing, multiple flux rope ejections emanated from the sheared arcade above the sharp polarity inversion line (see the animated version of Figure 64). The timing of these flux rope ejections do not coincide with the times of actual observed eruptions since the applied boundary condition was ad hoc. However, the model did give an enhanced horizontal field at the polarity inversion line following flux rope ejections. Such behavior was also reported in an observational study of AR 11158 using SDO/HMI data (Wang et al., 2012).

Figure 64:
figure 64

mpg-Movie (21783.9482422 kB KB) Still from a movie showing Left: Distribution of normalized current density (j/B, green) and magnetic field orientation (arrows) at a vertical cut along the dashed line shown in the right panel. Right: Magnetogram sequence (remapped from SDO/HMI data) used for the data-driven simulation of NOAA AR 11158 by Cheung and DeRosa (2012). Image reproduced with permission, copyright by AAS. (For video see appendix)

The time-evolutionary approach (whether using MF or MHD) should be clearly distinguished from other data-constrained modeling approaches like non-linear force-free field extrapolation (NLFFF; see Mikic and McClymont, 1994; McClymont et al., 1997; DeRosa et al., 2009, for extended discussions), which use individual vector magnetograms to extrapolate the instantaneous coronal field. In the work of Cheung and DeRosa (2012), the evolution of AR 11158 was modeled by driving the simulation with a temporal sequence of magnetograms to capture the AR from emergence to eruption (Gibb et al., 2014, used a similar approached for modeling AR 10977). This was motivated by previous studies using the MF approach to model the evolution of the global corona (Yeates et al., 2007, 2008; Yeates and Mackay, 2009a,b; Yeates et al., 2010). These studies demonstrate that the coronal magnetic configuration at a given time is dependent on prior photospheric evolution. While time-independent NLFFF models may be useful for studying the 3D coronal field in some scenarios (consult McClymont et al., 1997 and DeRosa et al., 2009 for critical assessments regarding NLFFF extrapolation), they do not provide information about how photospheric driving (as manifested in terms of the electric field) transports magnetic flux, energy, and helicity into the atmosphere. Since these central questions of flux emergence are not addressed, a detailed discussion of time-independent data-constrained models is outside the scope of this review.

To advance beyond data-inspired and (ad hoc) data-driven models of emerging and erupting ARs requires robust and reliable measurements of the photospheric vector velocity and magnetic fields. These quantities allow one to compute the −v × B electric field that is required for setting boundary conditions that faithfully describe transport processes associated with flux emergence. The measurement of such quantities is not only useful for MF models, but also for time-dependent MHD models that aim to capture the full dynamical consequences of emerging flux. Recent work on data-driven modeling is preceded by studies that attempt to constrain photospheric flows and their associated electric fields (Chae, 2001; Kusano et al., 2002; Demoulin and Berger, 2003; Longcope, 2004; Welsch et al., 2004, 2007; Schuck, 2005, 2006; Georgoulis and LaBonte, 2006; Ravindra et al., 2008). Advances in this direction of research depend on high-cadence vector magnetograms (e.g., from SDO/HMI or SOLIS) and observations of the atmospheric response to emerging flux. Reliable data sets, together with improvements in methods to constrain photospheric electric fields (Fisher et al., 2010, 2012), will likely stimulate substantial progress in data-driven modeling in the coming years.

6 Open Questions and Final Remarks

The study of flux emergence has received a great push forward from recent advances in observational capability, theory, and numerical modeling. The large number of flux emergence models has provided invaluable insight into the physical mechanisms that are central to the myriad of observable phenomena driven by the passage of magnetic flux from the solar convection zone all the way into the corona. The role of magnetic buoyancy instabilities to bring flux from the photosphere into the corona is well-studied. Furthermore, it is well established that granular convective flows play an important role in the morphology of emerging flux at the photosphere. The formation of sea serpentine field-lines (whether formed from interaction with granular flows or from the Parker instability) is understood to be an important step for removing the mass burden from field lines rising into the corona. Simulations dealing with the emergence of twisted flux ropes show how the Lorentz force drives rotational and shearing motions at the photosphere. The Poynting flux of energy and helicity associated with such motion is shown to be important even beyond the initial emergence phase (i.e., when the vertical unsigned flux at the photosphere is increasing) and are likely key for eruptive phenomena. Models of emerging flux interacting with pre-existing ambient field show how the former can sometimes destablize the pre-existing field to launch jets and CMEs.

A number of challanges and open questions remain. For instance, models of the full life-cyle of active regions from birth to decay remain to be developed. Such models may settle long standing questions about the depth from which active region fields are spawned which, in turn, tells us something about how the dynamo-generated magnetic field is structured in the convection zone. The decay of ARs is also an important topic that needs to be addresssed since the turbulent diffusion of AR-flux over the solar surface and the subduction of flux back into the convection zone likely determine how much remnant flux can be replenished for the large-scale dynamo field.

Helioseismic analysis of pre-emerging flux regions is still a nascent field. Numerical models are needed in order to validate helioseismic techniques which, in turn, are used on real observations to constrain numerical models. Continued work in this regard is likely essential for answering the question of whether we can robustly infer subsurface properties of emerging ARs well before their appearance at the surface.

In the near future, we anticipate substantial developments in data-driven models of emerging flux. Such models will use observed magnetogram sequences for boundary conditions of 3D MHD models. With such models, synthetic diagnostics can be directly compared with observations. Success in reproducing observables that quantatively match observations is crucial for a critical assessment of the maturity of MHD models. Furthermore, advances in this realm are essential for novel approaches to space weather forecasting that go beyond the use of heuristics and probabilistic models.