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 2 β = 8πp∕B (i.e., the ratio of gas and magnetic pressures). Assume the magnetic field is sufficiently weak such that 2 B ∕8π ≪ p (i.e., high plasma-β). In the case of total pressure equilibrium between plasma within the magnetic structure and the ambient plasma,

B2 pi + ---= pamb, (13 ) 8π
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
− 1 Δ ϱ∕ϱ ≈ − β . (14 )
If instead of thermal equilibrium, the interior and the external fluids possessed the same specific entropy s, the density deficit would be
( ) Δϱ ∕ϱ ≈ − (γ1β)−1, where γ1 = ∂-lnp- (15 ) ∂ lnϱ s
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 (τs ∼ L ∕cs, where L is a typical length-scale) or Alfvèn crossing time (τA ∼ L ∕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 dz = gϱ, 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).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 Dvz Dt = gΔ ϱ∕ϱ. 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

( ) −1 Hp = dlnp- , (16 ) dz
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*).
View Image
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.
View Image
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.

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).

View Image
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.

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.

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 ∝ ϱ. (17 )
View Image
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 B = B ˆx. For simplicity, assume that the upflow is axisymmetric with respect to the vertical direction ˆz. Let the gradients of the local velocity field be

∂vx- ∂vy- ∂x = ∂y = α, (18 ) ∂vz-= 𝜖α, (19 ) ∂z
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-ϱ Dt = − (2 + 𝜖)α, (20 ) D lnB -------= − (1 + 𝜖)α. (21 ) Dt
This results in the scaling relation between field strength and mass density as
κ 1-+-𝜖 B ∝ ϱ , where κ = 2 + 𝜖. (22 )
When the expansion is isotropic (𝜖 = 1), the field strength scales as 2∕3 B ∝ ϱ. 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 r ≪ Hp 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 r ∼ Hp, 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 1∕2 B ∝ ϱ. 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.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ϱ = ϱγ dpdz, (23 ) ad p 1 dz
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ϱ- dϱ = dz dz. (24 )
If ∂ϱad- ∂ϱ dz < 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 ∂ϱad ∂ϱ dz > 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 ∂ϱad-= ∂ϱ- dz dz, the stratification is neutral to adiabatic perturbations and is marginally stable. This implies the buoyancy force does not diminish nor amplify vertical perturbations.
View Image
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.

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

| dln-T-|| ∇ad = d ln p | , (25 ) s ∇ = dlnT-. (26 ) dln p
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.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).

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 m s–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 8 10 from the photosphere to the corona.

View Image
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.5

View Image
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 16 10 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.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.,

( 2) d-- B-- dz p + 8 π = − gϱ. (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.

View Image
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.

Figure 14* illustrates the different modes of the magnetic buoyancy instability. The k ⊥ B 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 14*b). 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-- B- C22 N-2 dz ln ρ < − g V2 , (28 ) A
where z is height, VA = B (4πϱ)−1∕2 is the Alfvén speed, N is the Brunt–Väisälä frequency, and ∘ ------ Cs = γ1p∕ϱ is the adiabatic sound speed. The Brunt–Väisälä frequency is defined by
g d N 2 = ------ln(pϱ− γ1). (29 ) γ1 dz
From the above expression, it is easy to see that N 2 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- B- dz ln ρ < 0. (30 )
When the layer is subadiabatically stratified (entropy increases with z), 2 N > 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*).

The k ∥ B mode is called the undular mode (Figure 14*c). 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 14*d and 14*e) is given by (Acheson, 1979*)

2 ( 2 2 ) d--ln B < − Cs- k2(1 + kz-) + N- , (31 ) dz g || k2⊥ V 2A
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,
2 2 d--ln B < − Csk||(k2 + k2) ≤ 0. (32 ) dz gk2⊥ ⊥ z
Namely, there exist unstable modes when
dB --- < 0. (33 ) dz
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 k y is that in the transverse direction, and k = ∘k2--+-k2- x y. The interchange mode (10 ky∕kx = 10 ∼ ∞) 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 k ∕k = 1010 y x mode is reduced, while that of the k ∕k = 0 y x mode is increased.

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*).

View Image
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.
View Image
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).

View Image
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.

View Image
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.

View Image
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.

View Image
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.

View Image
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.

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.

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.

View Image
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.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) = qrB (r) 𝜃 l, 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*).

View Image
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 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.

View Image
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.

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.

View Image
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).

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.01c s, where c s 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.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).

View Image
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).

View Image
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*).

View Image
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.

View Image
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.

View Image
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.

View Image
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.

View Image
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.

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.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.

View Image
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.

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.

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.

View Image
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.

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 Å 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 6 7 −2 − 1 10 –10 erg cm s. 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.

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.

View Image
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.

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 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.

View Image
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.

View Image
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.

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.

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 𝒞 of radius R (for the following discussion, consider 𝒞 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 -- B (r,𝜃,z) = B (r,z) + B ′(r,𝜃,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

-- ∂B--= ∇ × (v × B-) + ∇ × E, (34 ) ∂t
where E = v′ ×-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 𝒞 is then
∫ -- ∮ ∂B-- ( -- -- ) Φ˙𝒞 = ∂t ⋅ dS = v × B + E ⋅ dl, (35 ) 𝒞 ∂𝒞 = ˙Φm + ˙Φf, (36 )
˙ -- -- ---- -- -- Φm = 2πR [v × B ]𝜃 =-2πR-(vzBr-− Bz vr), (37 ) Φ˙f = 2πR ℰ𝜃 = 2πR v′zB ′r − B ′zv′r. (38 )
Equation (36*) indicates that the magnetic flux enclosed within 𝒞 changes as a result of (1) advection of the mean magnetic field -- B by the mean flow -- 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 ˙ Φ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.

View Image
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 Bl(r)ˆl and B𝜃(r)ˆ𝜃 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

c ∫ c∫ 2π cr I(r) = --- (∇ × B )l ⋅ dS =-- B 𝜃rd𝜃 = --⟨B 𝜃(r)⟩. (39 ) 4π 2 0 2
In numerical experiments modeling the dynamics of twist flux tubes, the magnetic profiles B (r ) l 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

2 2 Bl(r) = B0 exp(− r ∕a ), (40 ) B 𝜃(r) = qrBl(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.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

-v2ϱ--- W e = B2∕4 π, (42 ) t
where v is the rise velocity of the tube. The criterion for maintaining coherence of the flux tube is W e ≲ 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

−1 qcr = a , (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 40Watch/download Movie 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 Alfvén speeds comparable to the local sound speed.
View Image
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: gif-Movie (736 KB) Volumetric rendering of a rising flux tube experiencing the developing of the helical kink instability (Fan et al., 1998).

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 𝒫. The line integral of the gas and magnetic pressure gradient forces about this curve vanish. That is,

∮ ∮ ( 2 ) r × ∇ (− p ) ⋅ dl = r × ∇ − B-- ⋅ dl = 0, (44 ) gas 8π
where r is the displacement vector of a point in the curve from 𝒫. 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
∮ [ ] tz = r × 1-(B ⋅ ∇)B ⋅ dl, (45 ) 4π
which cannot be balanced by pressure gradient forces. This unbalanced torque drives rotational motions at the photosphere.10
View Image
Figure 41: Distribution of the vertical component of vorticity (ωz) 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 ωz. Image reproduced with permission from Fan (2009a*), copyright by AAS.

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 ωz 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*).

View Image
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.

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.

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 2 α = (∇ × B ) ⋅ B ∕B. 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

1--4 2∂-α τl = 16a B l ∂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 h = j ⋅ B c. Since vector magnetograms at a single height only give the contribution to h c 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 −1 c j × B vanishes everywhere in the volume. This implies that the current density j is proportional to B, or

∇ × B = α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 −8 −1 1 × 10 m to −8 −1 5 × 10 m (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 = V A ⋅ BdV , (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 ∫ 2 V B ∕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 B n 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

∫ ∫ HR = V (A ⋅ B )dV − V(A0 ⋅ B0)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 B0 ⋅ ˆn|∂V = B ⋅ ˆn|∂V.

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

∫ -1- 2 2 EF (B, P ) = 8π |B| − |P | dV , (50 ) V
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

∮ dHR--= 2 (Ap × E ) ⋅ ˆndS, (51 ) dt ∂V
where A p is the vector potential field of the potential field satisfying ∇ ⋅ A = 0 p and A ⋅ ˆn | = 0 p ∂V (Berger and Field, 1984). Under ideal evolution, E = − v × B so that
∮ dHR-- dt = − 2 [(Ap ⋅ B )v − (Ap ⋅ v )B ] ⋅ ˆndS. (52 ) ∂V
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 ∫ dHRdt 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
dHR ∫ -----= 2 [(Ap ⋅ B)vn − (Ap ⋅ v)Bn ]dS, (53 ) dt
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.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.

View Image
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. (2013), 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.

View Image
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--- − 16----T----- σ ≈ 10 6 −3∕2(e.s.u). (54 ) Spitzer 10 K

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

( ) ( ) ( ) S ≈ 10−13 ----T----- --L---- ----VA----- . (55 ) 106 K −3∕2 109 cm 108 cm s−1

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.,

η = η0 (vd∕vcri − 1)α + ηc for vd ≥ vcri, (56 ) = ηc for vd < vcri. (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 −1∕2 Vinflow ∕VA = S (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.

View Image
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 (Karlický 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

− 5 5∕2 κ = 1.84-×-10--T----≈ 10−6 erg s−1 K− 1 cm −1, (58 ) lnΛ
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 (5∕2 κ ∝ T), thermal conduction is very effective in the high temperature (6 T > 10 K) corona, whereas it is almost negligible in the chromosphere and below.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):

[ 2 ] ∂B-- 4π-η -1-- --D--- ∂t = ∇ × v × B − c j − enej × B + (j × B ) × B cϱiν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

[ ] ∂B--= ∇ × u × B − 4π-ηj , (60 ) ∂t c
u = v + v + v , (61 ) Hall Amb vHall = − H ∇ × B, (62 ) v = M (∇ × B ) × B, (63 ) Amb c H = ------, (64 ) 4πene --D2--- M = 4π ϱ ν . (65 ) i in
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πηj) = ∇ × (− η∇ × B ) c 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.

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

2 [ ] ∂(B--∕8π)-+ ∇ ⋅ -c-E × B = E ⋅ j. (66 ) ∂t 4π
In terms of the Poynting flux c-E × B 4π, 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 2 2ϱv, 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 2 EAmb ⋅ j = − 4πM FL, 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*).

Table 1: Summary of the impact of various effects present in MHD of partially ionized plasmas.
Plasma Flow Ohmic Hall Ambipolar
Maintains topology (i.e., E perpendicular to B)? Yes No Yes Yes
Contributes to Poynting flux? Yes Yes Yes Yes
Dissipative? No Yes No Yes

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 ω 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

E = − v × B + ηj∥ + ηcj⊥, (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 = η.
View Image
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.

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.

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 ϱ, temperature T, and field strength B. The energy equation includes a Newton cooling term given by

d-𝜀 𝜀 −-𝜀0(ϱ) dt = − τ(ϱ) , (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 (2 3 10 – 10 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.

  Go to previous page Scroll to top Go to next page