1 Introduction

Looking at a full disk magnetogram (a map showing spatially the line of sight flux density of the magnetic field) of the solar photosphere one sees that the most prominent large scale pattern of magnetic flux concentrations on the solar surface are the bipolar active regions (see Figure 1). When observed in white light (see Figure 2), an active region usually contains sunspots and is sometimes called a sunspot group. Active regions are so named because they are centers of various forms of solar activity (such as solar flares) and sites of X-ray emitting coronal loops (see Figure 3). Despite the turbulent nature of solar convection which is visible from the granulation pattern on the photosphere, the large scale bipolar active regions show remarkable order and organization as can be seen in Figure 1. The active regions are roughly confined into two latitudinal belts which are located nearly symmetrically on the two hemispheres. Over the course of each 11-year solar cycle, the active region belts march progressively from mid-latitude of roughly 35° to the equator on both hemispheres (Maunder, 1922). The polarity orientations of the bipolar active regions are found to obey the well-known Hale polarity law (Hale et al., 1919; Hale and Nicholson, 1925) outlined as follows. The line connecting the centers of the two magnetic polarity areas of each bipolar active region is usually nearly east-west oriented. Within each 11-year solar cycle, the leading polarities (leading in the direction of solar rotation) of nearly all active regions on one hemisphere are the same and are opposite to those on the other hemisphere (see Figure 1), and the polarity order reverses on both hemispheres with the beginning of the next cycle. The magnetic fields at the solar north and south poles are also found to reverse sign every 11 years near sunspot maximum (i.e. near the middle of a solar cycle). Therefore, the complete magnetic cycle, which corresponds to the interval between successive appearances at mid-latitudes of active regions with the same polarity arrangement, is in fact 22 years.

Figure 1:
figure 1

A full disk magnetogram from the Kitt Peak Solar Observatory showing the line of sight magnetic flux density on the photosphere of the Sun on May 11, 2000. White (Black) color indicates a field of positive (negative) polarity.

Figure 2:
figure 2

A continuum intensity image of the Sun taken by the MDI instrument on board the SOHO satellite on the same day as Figure 1 . It shows the sunspots that are in some of the active regions in Figure 1 .

Figure 3:
figure 3

A full disk soft X-ray image of the solar coronal taken on the same day as Figure 1 from the soft X-ray telescope on board the Yohkoh satellite. Active regions appear as sites of bright X-ray emitting loops.

Besides their highly organized behavior during each solar cycle, active regions are found to possess some interesting asymmetries between their leading and following polarities. Observations show that the axis connecting the leading and the following polarities of each active region is nearly east-west oriented (or toroidal) but on average shows a small tilt relative to the east-west direction with the leading polarity of the region being slightly closer to the equator than the following (see Figure 1). This small mean tilt angle is found to increase approximately linearly with the latitude of the active region (Wang and Sheeley Jr, 1989, 1991; Howard, 1991a,b; Fisher et al., 1995). This observation of active region tilts is originally summarized in Hale et al. (1919) and is generally referred to as Joy’s law. Note that Joy’s law describes the statistical mean behavior of the active region tilts. The tilt angles of individual active regions also show a large scatter about the mean (Wang and Sheeley Jr, 1989, 1991; Howard, 1991a,b; Fisher et al., 1995). Another intriguing asymmetry is found in the morphology of the leading and the following polarities of an active region. The flux of the leading polarity tends to be concentrated in large well-formed sunspots, whereas the flux of the following polarity tends to be more dispersed and to have a fragmented appearance (see Bray and Loughhead, 1979). Observations also show that the magnetic inversion lines (the neutral lines separating the fluxes of the two opposite polarities) in bipolar active regions are statistically nearer to the main following polarity spot than to the main leading spot (van Driel-Gesztelyi and Petrovay, 1990; Petrovay et al., 1990). Furthermore for young growing active regions, there is an asymmetry in the east-west proper motions of the two polarities, with the leading polarity spots moving prograde more rapidly than the retrograde motion of the following polarity spots (see Chou and Wang, 1987; van Driel-Gesztelyi and Petrovay, 1990; Petrovay et al., 1990).

More recently, vector magnetic field observations of active regions on the photosphere have shown that active region magnetic fields have a small but statistically significant mean twist that is left-handed in the northern hemisphere and right-handed in the southern hemisphere (see Pevtsov et al., 1995, 2001). The twist is measured in terms of the quantity α ≡ 〈J z/B z〉, the ratio of the vertical electric current over the vertical magnetic field averaged over an active region. The measured α for individual solar active regions show considerable scatter, but there is clearly a statistically significant trend for negative α (left-handed field line twist) in the northern hemisphere and positive α (right-handed field line twist) in the southern hemisphere. In addition, soft X-ray observations of solar active regions sometimes show hot plasma of S or inverse-S shapes called “sigmoids” with the northern hemisphere preferentially showing inverse-S shapes and the southern hemisphere forward-S shapes (Rust and Kumar, 1996; Pevtsov et al., 2001, see Figure 4 for an example). This hemispheric preference of the sign of the active region field line twist and the direction of X-ray sigmoids do not change with the solar cycle (see Pevtsov et al., 2001).

Figure 4:
figure 4

A soft X-ray image of the solar coronal on May 27, 1999, taken by the Yohkoh soft X-ray telescope. The arrows point to two “sigmoids” at similar longitudes north and south of the equator showing an inverse-S and a forward-S shape respectively.

The cyclic large scale magnetic field of the Sun with a period of 22 years is believed to be sustained by a dynamo mechanism. The Hale polarity law of solar active regions indicates the presence of a large scale subsurface toroidal magnetic field generated by the solar dynamo mechanism. In the past decade, the picture of how and where the large scale solar dynamo operates has undergone substantial revision due in part to new knowledge from helioseismology regarding the solar internal rotation profile (see Deluca and Gilman, 1991; Gilman, 2000). Evidence now points to the tachocline, the thin shear layer at the base of the solar convection zone, where solar rotation changes from the latitudinal differential rotation of the solar convective envelope to the nearly solid-body rotation of the radiative interior, as the site for the generation and amplification of the large scale toroidal magnetic field from a weak poloidal magnetic field (see Charbonneau and MacGregor, 1997; Dikpati and Charbonneau, 1999; Dikpati and Gilman, 2001). Furthermore, with its stable (weakly) subadiabatic stratification, the thin overshoot region in the upper part of the tachocline layer (Gilman, 2000) allows storage of strong toroidal magnetic fields against their magnetic buoyancy for time scales comparable to the solar cycle period (Parker, 1975, 1979; van Ballegooijen, 1982; Moreno-Insertis et al., 1992; Fan and Fisher, 1996; Moreno-Insertis et al., 2002; Rempel, 2003). Thus with toroidal magnetic fields being generated and stored in the tachocline layer at the base of the solar convection zone, these fields need to traverse the entire convection zone before they can emerge at the photosphere to form the observed solar active regions.

High resolution observations have shown that magnetic fields on the solar photosphere are in a fibril state, i.e. in the form of discrete flux tubes of high field strength (B ≳ 103 G in equipartition with the thermal pressure) having a hierarchy of cross-sectional sizes that range from sunspots of active regions down to below the limit of observational resolution (see Zwaan, 1987; Stenflo, 1989; Domínguez Cerdeña et al., 2003; Khomenko et al., 2003; Socas-Navarro and Sánchez Almeida, 2003). It is thus likely that the subsurface magnetic fields in the solar convection zone are also concentrated into discrete flux tubes. One mechanism that can concentrate magnetic flux in a turbulent conducting fluid, such as the solar convection zone, into high field strength flux tubes is the process known as “flux expulsion”, i.e. magnetic fields are expelled from the interior of convecting cells into the boundaries. This process has been studied by MHD simulations of the interaction between convection and magnetic fields (see Galloway and Weiss, 1981; Nordlund et al., 1992). In particular, the 3D simulations of magnetic fields in convecting flows by Nordlund et al. (1992) show the formation of strong discrete flux tubes in the vicinity of strong downdrafts. In addition, Parker (1984) put forth an interesting argument that supports the fibril form of magnetic fields in the solar convection zone. He points out that although the magnetic energy is increased by the compression from a continuum field into the fibril state, the total energy of the convection zone (thermal + gravitational + magnetic) is reduced by the fibril state of the magnetic field by avoiding the magnetic inhibition of convective overturning. Assuming an idealized polytropic atmosphere, he was able to derive the filling factor of the magnetic fields that corresponds to the minimum total energy state of the atmosphere. By applying an appropriate polytropic index for the solar convection zone, he computed the filling factor which yielded fibril magnetic fields of about 1–5 kG, roughly in agreement with the observed fibril fields at the solar surface.

Since both observational evidence and theoretical arguments support the fibril picture of solar magnetic fields, the concept of isolated magnetic flux tubes surrounded by “field-free” plasmas has been developed and widely used in modeling magnetic fields in the solar convection zone (see Parker, 1979; Spruit, 1981; Vishniac, 1995a,b). The manner in which individual solar active regions emerge at the photosphere (see Zwaan, 1987) and the well-defined order of the active regions as described by the Hale polarity rule suggest that they correspond to coherent and discrete flux tubes rising through the solar convection zone and reaching the photosphere in a reasonably cohesive fashion, not severely distorted by convection. It is this process — the formation of buoyant flux tubes from the toroidal magnetic field stored in the overshoot region and their dynamic rise through the convection zone to form solar active regions — that is the central focus of this review.

The remainder of the review will be organized as follows.

  • Section 2 gives a brief overview of the simplifying models and computational approaches that have been applied to studying the dynamic evolution of magnetic flux tubes in the solar convection zone. In particular, the thin flux tube model is proven to be a very useful tool for understanding the global dynamics of emerging active region flux tubes in the solar convective envelope and (as discussed in the later sections) has produced results that explain the origin of several basic observed properties of solar active regions.

  • Section 3 discusses the storage and equilibrium properties of large scale toroidal magnetic fields in the stable overshoot region below the solar convection zone.

  • Section 4 focuses on the buoyancy instabilities associated with the equilibrium toroidal magnetic fields and the formation of buoyant flux tubes from the base of the solar convection zone.

  • Section 5 reviews results on the dynamic evolution of emerging flux tubes in the solar convection zone.

    • Section 5.1 discusses major findings from various thin flux tube simulations of emerging flux loops.

    • Section 5.2 discusses the observed hemispheric trend of the twist of the magnetic field in solar active regions and the models that explain its origin.

    • Section 5.3 reviews results from direct MHD simulations with regard to the minimum twist necessary for tube cohesion.

    • Section 5.4 discusses the kink evolution of highly twisted emerging tubes.

    • Section 5.5 reviews the influence of 3D stratified convection on the evolution of buoyant flux tubes.

  • Section 6 discusses results from 3D MHD simulations on the asymmetric transport of magnetic flux (or turbulent pumping of magnetic fields) by stratified convection penetrating into a stable overshoot layer.

  • Section 7 discusses an alternative mechanism of magnetic flux amplification by converting the potential energy associated with the stratification of the convection zone into magnetic energy.

  • Section 8 gives a brief overview of our current understanding of active region flux emergence at the surface and the post-emergence evolution of the subsurface fields.

  • Section 9 gives a summary of the basic conclusions.

2 Models and Computational Approaches

2.1 The thin flux tube model

The well-defined order of the solar active regions (see description of the observational properties in Section 1) suggests that their precursors at the base of the solar convection zone should have a field strength that is at least B eq, where B eq is the field strength that is in equipartition with the kinetic energy density of the convective motions: \(B_{{\rm eq}}^{2}/8\pi=\rho v_{c}^{2}/2\). If we use the results from the mixing length models of the solar convection zone for the convective flow speed v c, then we find that in the deep convection zone B eq is on the order of 104 G. In the past two decades, direct 3D numerical simulations have led to a new picture for solar convection that is non-local, driven by the concentrated downflow plumes formed by radiative cooling at the surface layer, and with extreme asymmetry between the upward and downward flows (see reviews by Spruit et al., 1990; Spruit, 1997). Hence it should be noted that the B eq derived based on the local mixing length description of solar convection may not really reflect the intensity of the convective flows in the deep solar convection zone. With this caution in mind, we nevertheless refer to B eq ∼ 104 G as the field strength in equipartition with convection in this review.

Assuming that in the deep solar convection zone the magnetic field strength for flux tubes responsible for active region formation is at least 104 G, and given that the amount of flux observed in solar active regions ranges from ∼ 1020 Mx to 1022 Mx (see Zwaan, 1987), then one finds that the cross-sectional sizes of the flux tubes are small in comparison to other spatial scales of variation, e.g. the pressure scale height. For an isolated magnetic flux tube that is thin in the sense that its cross-sectional radius a is negligible compared to both the scale height of the ambient unmagnetized fluid and any scales of variation along the tube, the dynamics of the flux tube may be simplified with the thin flux tube approximation (see Spruit, 1981; Longcope and Klapper, 1997) which corresponds to the lowest order in an expansion of MHD in powers of a/L, where L represents any of the large length scales of variation. Under the thin flux tube approximation, all physical quantities of the tube, such as position, velocity, field strength, pressure, density, etc. are assumed to be averages over the tube cross-section and they vary spatially only along the tube. Furthermore, because of the much shorter sound crossing time over the tube diameter compared to the other relevant dynamic time scales, an instantaneous pressure balance is assumed between the tube and the ambient unmagnetized fluid:

$$p + {{{B^2}} \over {8\pi }} = {p_{\rm{e}}}$$
(1)

where p is the tube internal gas pressure, B is the tube field strength, and p e is the pressure of the external fluid. Applying the above assumptions to the ideal MHD momentum equation, Spruit (1981) derived the equation of motion of a thin untwisted magnetic flux tube embedded in a field-free fluid. Taking into account the differential rotation of the Sun, \({{\bf{\Omega }}_{\rm{e}}}({\bf{r}}) = {\Omega _{\rm{e}}}({\bf{r}}){\bf{\hat z}}\), the equation of motion for the thin flux tube in a rotating reference frame of angular velocity \({\bf{\Omega }} = \Omega {\bf{\hat z}}\) is (Ferriz-Mas and Schüssler, 1993; Caligari et al., 1995)

$$\matrix{{\rho {{d{\bf{v}}} \over {dt}} = 2\rho ({\bf{v}} \times \Omega) + \rho ({\Omega ^2} - \Omega _{\rm{e}}^2)\varpi \hat \varpi + (\rho - {\rho _{\rm{e}}}){{\bf{g}}_{{\rm{eff}}}}} \hfill \cr {\quad \quad \quad + {\rm{\hat l}}{\partial \over {\partial s}}\left({{{{B^2}} \over {8\pi }}} \right) + {{{B^2}} \over {4\pi }}{\bf{k}} - {C_{\rm{D}}}{{{\rho _{\rm{e}}}\vert{{({{\bf{v}}_{{\rm{rel}}}})}_ \bot }\vert{{({{\bf{v}}_{{\rm{rel}}}})}_ \bot }} \over {\pi {{(\Phi/B)}^{1/2}}}},} \hfill \cr}$$
(2)

where

$${{\bf{g}}_{{\rm{eff}}}} = {\bf{g}} + \Omega _{\rm{e}}^2\varpi \hat \varpi,$$
(3)
$${({{\bf{v}}_{{\rm{rel}}}})_ \bot } = {[{\bf{v}} - ({{\bf{\Omega }}_{\rm{e}}} - {\bf{\Omega }}) \times {\bf{r}}]_ \bot }.$$
(4)

In the above, r, v, B, p, ρ, denote the position vector, velocity, magnetic field strength, plasma pressure and density of a Lagrangian tube element respectively, each of which is a function of time t and the arc-length s measured along the tube, ρ e(r) denotes the external density at the position r of the tube element, \({\bf{\hat z}}\) is the unit vector pointing in the direction of the solar rotation axis, \(\hat{\varpi}\) denotes the unit vector perpendicular to and pointing away from the rotation axis at the location of the tube element and ϖ denotes the distance to the rotation axis, \(\hat{\bf l}\equiv\partial{\bf r}\partial s\) is the unit vector tangential to the flux tube, k 2r/∂s 2 is the tubes curvature vector, the subscript ⊥ denotes the vector component perpendicular to the local tube axis, g is the gravitational acceleration, and C D is the drag coefficient. The drag term (the last term on the right hand side of the equation of motion (2)) is added to approximate the opposing force experienced by the flux tube as it moves relative to the ambient fluid. The term is derived based on the case of incompressible flows pass a rigid cylinder under high Reynolds number conditions, in which a turbulent wake develops behind the cylinder, creating a pressure difference between the up- and down-stream sides and hence a drag force on the cylinder (see Batchelor, 1967).

If one considers only the solid body rotation of the Sun, then the Equations (2), (3), and (4) can be simplified by letting Ω e = Ω. Calculations using the thin flux tube model (see Section 5.1) have shown that the effect of the Coriolis force 2ρ(v × Ω) acting on emerging flux loops can lead to east-west asymmetries in the loops that explain several well-known properties of solar active regions.

Note that in the equation of motion (2), the effect of the “enhanced inertia” caused by the back-reaction of the fluid to the relative motion of the flux tube is completely ignored. This effect has sometimes been incorporated by treating the inertia for the different components of Equation (2) differently, with the term ρ(dv/dt) on the left-hand-side of the perpendicular component of the equation being replaced by (ρ + ρ e)(dv/dt) (see Spruit, 1981). This simple treatment is problematic for curved tubes and the proper ways to treat the back-reaction of the fluid are controversial in the literature (Cheng, 1992; Fan et al., 1994; Moreno-Insertis et al., 1996). Since the enhanced inertial effect is only significant during the impulsive acceleration phases of the tube motion, which occur rarely in the thin flux tube calculations of emerging flux tubes, and the results obtained do not depend significantly on this effect, many later calculations have taken the approach of simply ignoring it (see Caligari et al., 1995, 1998; Fan and Fisher, 1996).

Equations (1) and (2) are to be complemented by the following equations to completely describe the dynamic evolution of a thin untwisted magnetic flux tube:

$${d \over {dt}}\left({{B \over \rho }} \right) = {B \over \rho }\left[ {{{\partial ({\bf{v}}\cdot{\rm{\hat l}})} \over {\partial s}} - {\bf{v}}\cdot{\bf{k}}} \right],$$
(5)
$${1 \over \rho }{{d\rho } \over {dt}} = {1 \over {\gamma p}}{{dp} \over {dt}} - {{{\nabla _{{\rm{ad}}}}} \over p}{{dQ} \over {dt}},$$
(6)
$$p = {{\rho RT} \over \mu },$$
(7)

where ∇ad ≡ ( ln T/∂ lnp)s. Equation (5) describes the evolution of the tube magnetic field and is derived from the ideal MHD induction equation (Spruit, 1981). Equation (6) is the energy equation for the thin flux tube (Fan and Fisher, 1996), in which dQ/dt corresponds to the volumetric heating rate of the flux tube by non-adiabatic effects, e.g. by radiative diffusion (Section 3.2). Equation (7) is simply the equation of state for an ideal gas. Thus the five Equations (1), (2), (5), (6), and (7) completely determine the evolution of the five dependent variables v (t, s), B (t, s), p (t, s), ρ (t, s), and T (t, s) for each Lagrangian tube element of the thin flux tube.

Spruit’s original formulation for the dynamics of a thin isolated magnetic flux tube as described above assumes that the tube consists of untwisted flux \({\bf B}=B\hat{\bf l}\). Longcope and Klapper (1997) extend the above model to include the description of a weak twist of the flux tube, assuming that the field lines twist about the axis at a rate q whose magnitude is 2π/L w, where L w is the distance along the tube axis over which the field lines wind by one full rotation and |qa| ≪ 1. Thus in addition to the axial component of the field B, there is also an azimuthal field component in each tube cross-section, which to lowest order in qa is given by B θ = qr B, where r denotes the distance to the tube axis. An extra degree of freedom for the motion of the tube element — the spin of the tube cross-section about the axis — is also introduced, whose rate is denoted by w (angle per unit time). By considering the kinematics of a twisted ribbon with one edge corresponding to the tube axis and the other edge corresponding to a twisted field line of the tube, Longcope and Klapper (1997) derived an equation that describes the evolution of the twist q in response to the motion of the tube:

$${{dq} \over {dt}} = - {{d\ln \delta s} \over {dt}}q + {{\partial \omega } \over {\partial s}} + ({\rm{\hat l}} \times {\bf{k}})\cdot{{d{\rm{\hat l}}} \over {dt}},$$
(8)

where δs denotes the length of a Lagrangian tube element. The first term on the right-hand-side describes the effect of stretching on q: Stretching the tube reduces the rate of twist q. The second term is simply the change of q resulting from the gradient of the spin along the tube. The last term is related to the conservation of total magnetic helicity which, for the thin flux tube structure, can be decomposed into a twist component corresponding to the twist of the field lines about the axis, and a writhe component corresponding to the “helicalness” of the axis (see discussion in Longcope and Klapper, 1997). It describes how the writhing motion of the tube axis can induce twist of the opposite sense in the tube.

Furthermore, by integrating the stresses over the surface of a tube segment, Longcope and Klapper (1997) evaluated the forces experienced by the tube segment. They found that for a weakly twisted (|qa| ≪ 1) thin tube (|a∂ s| ≪ 1), the equation of motion of the tube axis differs very little from that for an untwisted tube — the leading order term in the difference is O[qa 2 s] (see also Ferriz-Mas and Schüssler, 1990). Thus the equation of motion (2) applies also to a weakly twisted thin flux tube. By further evaluating the torques exerted on a tube segment, Longcope and Klapper (1997) also derived an equation for the evolution of the spin w:

$${{d\omega} \over {dt}} = - {2 \over a}{{da} \over {dt}}\omega + v_{\rm{a}}^2{{\partial q} \over {\partial s}},$$
(9)

where \(v_{\rm a}=B/\sqrt{4\pi\rho}\) is the Alfvén speed. The first term on the right hand side simply describes the decrease of spin due to the expansion of the tube cross-section as a result of the tendency to conserve angular momentum. The second term, in combination with the second term on the right hand side of Equation (8), describes the propagation of torsional Alfvén waves along the tube.

The two new Equations (8) and (9) — derived by Longcope and Klapper (1997) — together with the earlier Equations (1), (2), (5), (6), and (7) provide a description for the dynamics of a weakly twisted thin flux tube. Note that the two new equations are decoupled from and do not have any feedback on the solutions for the dependent variables described by the earlier equations. One can first solve for the motion of the tube axis using Equations (1), (2), (5), (6), and (7), and then apply the resulting motion of the tube axis to Equations (8) and (9) to determine the evolution of the twist of the tube. If the tube is initially twisted, then the twist q can propagate and re-distribute along the tube as a result of stretching (1st term on the right-hand-side of Equation (8)) and the torsional Alfvén waves (2nd term on the right-hand-side of Equation (8)). Twist can also be generated due to writhing motion of the tube axis (last term on the right-hand-side of Equation (8)), as required by the conservation of total helicity.

2.2 MHD simulations

The thin flux tube (TFT) model described above is physically intuitive and computationally tractable. It provides a description of the dynamic motion of the tube axis in a three-dimensional space, taking into account large scale effects such as the curvature of the convective envelope and the Coriolis force due to solar rotation. The Lagrangian treatment of each tube segment in the TFT model allows for preserving perfectly the frozen-in condition of the tube plasma. Thus there is no magnetic diffusion in the TFT model. However, the TFT model ignores variations within each tube cross-section. It is only applicable when the flux tube radius is thin (Section 2.1) and the tube remains a cohesive object (Section 5.3). Clearly, to complete the picture, direct MHD calculations that resolve the tube cross-section and its interaction with the surrounding fluid are needed. On the other hand, direct MHD simulations that discretize the spatial domain are subject to numerical diffusion. The need to adequately resolve the flux tube — so that numerical diffusion does not have a significant impact on the dynamical processes of interest (e.g. the variation of magnetic buoyancy) — severely limits the spatial extent of the domain that can be modeled. So far the MHD simulations cannot address the kinds of large scale dynamical effects that have been studied by the TFT model (Section 5.1). Thus the TFT model and the resolved MHD simulations complement each other.

For the bulk of the solar convection zone, the fluid stratification is very close to being adiabatic with δ ≪ 1, where δ ≡ ∇ − ∇ad is the non-dimensional superadiabaticity with ∇ = d ln T/d ln p and ∇ad = (d ln T/d ln p)ad denoting the actual and the adiabatic logarithmic temperature gradient of the fluid respectively, and the convective flow speed v c is expected to be much smaller than the sound speed c s: v c/c sδ 1/2 ≪ 1 (see Schwarzschild, 1958; Lantz, 1991). Furthermore, the plasma β defined as the ratio of the thermal pressure to the magnetic pressure (βp/(B 2/8π)) is expected to be very high (β ≫ 1) in the deep convection zone. For example for flux tubes with field strengths of order 105 G, which is significantly super-equipartition compared to the kinetic energy density of convection, the plasma β is of order 105. Under these conditions, a very useful computational approach for modeling subsonic magnetohydrodynamic processes in a pressure dominated plasma is the well-known anelastic approximation (see Gough, 1969; Gilman and Glatzmaier, 1981; Glatzmaier, 1984; Lantz and Fan, 1999). The main feature of the anelastic approximation is that it filters out the sound waves so that the time step of numerical integration is not limited by the stringent acoustic time scale which is much smaller than the relevant dynamic time scales of interest as determined by the flow velocity and the Alfvén speed.

Listed below is the set of anelastic MHD equations (see Gilman and Glatzmaier, 1981; Lantz and Fan, 1999, for details of the derivations):

$$\nabla \cdot ({\rho _0}{\bf{v}}) = 0,$$
(10)
$${\rho _0}\left[ {{{\partial {\bf{v}}} \over {\partial t}} + ({\bf{v}}\cdot\nabla){\bf{v}}} \right] = - \nabla {p_1} + {\rho _1}{\bf{g}} + {1 \over {4\pi }}(\nabla \times {\bf{B}}) \times {\bf{B}} + \nabla \cdot{\bf{\Pi }},$$
(11)
$${\rho _0}{T_0}\left[ {{{\partial {s_1}} \over {\partial t}} + ({\bf{v}}\cdot\nabla)({s_0} + {s_1})} \right] = \nabla \cdot(K{\rho _0}{T_0}\nabla {s_1}) + {1 \over {4\pi }}\eta \vert\nabla \times {\bf{B}}{\vert^2} + ({\bf{\Pi }}\cdot\nabla)\cdot{\bf{v}},$$
(12)
$$\nabla \cdot {\bf{B}} = 0,$$
(13)
$${{\partial {\bf{B}}} \over {\partial t}} = \nabla \times ({\bf{v}} \times {\bf{B}}) - \nabla \times (\eta \nabla \times {\bf{B}}),$$
(14)
$${{{\rho _1}} \over {{\rho _0}}} = {{{p_1}} \over {{p_0}}} - {{{T_1}} \over {{T_0}}},$$
(15)
$${{{s_1}} \over {{c_p}}} = {{{T_1}} \over {{T_0}}} - {{\gamma - 1} \over \gamma }{{{p_1}} \over {{p_0}}},$$
(16)

where s 0(z), p 0(z), ρ 0(z), and T 0(z) correspond to a time-independent, background reference state of hydrostatic equilibrium and nearly adiabatic stratification, and velocity v, magnetic field B, thermodynamic fluctuations s1, p 1, ρ 1, and T 1 are the dependent variables to be solved that describe the changes from the reference state. The quantity Π is the viscous stress tensor given by

$${\Pi _{ij}} \equiv \mu \left({{{\partial {v_i}} \over {\partial {x_j}}} + {{\partial {v_j}} \over {\partial {x_i}}} - {2 \over 3}(\nabla \cdot{\bf{v}}){\delta _{ij}}} \right),$$

and μ, K and η denote the dynamic viscosity, and thermal and magnetic diffusivity, respectively. The anelastic MHD equations (10)(16) are derived based on a scaled-variable expansion of the fully compressible MHD equations in powers of δ and β −1, which are both assumed to be quantities ≪ 1. To first order in δ, the continuity equation (10) reduces to the statement that the divergence of the mass flux equals to zero. As a result sound waves are filtered out, and pressure is assumed to adjust instantaneously in the fluid as if the sound speed was infinite. Although the time derivative of density no longer appears in the continuity equation, density ρ 1 does vary in space and time and the fluid is compressible but on the dynamic time scales (as determined by the flow speed and the Alfvén speed) not on the acoustic time scale, thus allowing convection and magnetic buoyancy to be modeled in the highly stratified solar convection zone. Fan (2001a) has shown that the anelastic formulation gives an accurate description of the magnetic buoyancy instabilities under the conditions of high plasma β and nearly adiabatic stratification.

Fully compressible MHD simulations have also been applied to study the dynamic evolution of a magnetic field in the deep solar convection zone using non-solar but reasonably large β values such as β ∼ 10 to 1000. In several cases comparisons have been made between fully compressible simulations using large plasma β and the corresponding anelastic MHD simulations, and good agreement was found between the results (see Fan et al., 1998a; Rempel, 2002). Near the top of the solar convection zone, neither the TFT model nor the anelastic approximation are applicable because the active region flux tubes are no longer thin (Moreno-Insertis, 1992) and the velocity field is no longer subsonic. Fully compressible MHD simulations are necessary for modeling flux emergence near the surface (Section 8).

3 Equilibrium Conditions of Toroidal Magnetic Fields Stored at the Base of the Solar Convection Zone

3.1 The mechanical equilibria for an isolated toroidal flux tube or an extended magnetic layer

The Hale’s polarity rule of solar active regions indicates a subsurface magnetic field that is highly organized, of predominantly toroidal direction, and with sufficiently strong field strength (super-equipartition compared to the kinetic energy density of convection) such that it is not subjected to strong deformation by convective motions. It is argued that the weakly subadiabatically stratified overshoot layer at the base of the solar convection zone is the most likely site for the storage of such a strong coherent toroidal magnetic field against buoyant loss for time scales comparable to the solar cycle period (see Parker, 1979; van Ballegooijen, 1982).

It is not clear if the toroidal magnetic field is in the state of isolated flux tubes or stored in the form of a more diffuse magnetic layer. Moreno-Insertis et al. (1992) have considered the mechanical equilibrium of isolated toroidal magnetic flux tubes (flux rings) in a subadiabatic layer using the thin flux tube approximation (Section 2.1). The forces experienced by an isolated toroidal flux ring at the base of the convection zone is illustrated in Figure 5(a). The condition of total pressure balance (1) and the presence of a magnetic pressure inside the flux tube require a lower gas pressure inside the flux tube compared to the outside. Thus either the density or the temperature inside the flux tube needs to be lower. If the flux tube is in thermal equilibrium with the surrounding, then the density inside needs to be lower and the flux tube is buoyant. The buoyancy force associated with a magnetic flux tube in thermal equilibrium with its surrounding is often called the magnetic buoyancy (Parker, 1975). It can be seen in Figure 5(a) that a radially directed buoyancy force has a component that is parallel to the rotation axis, which cannot be balanced by any other forces associated with the toroidal flux ring. Thus for the toroidal flux ring to be in mechanical equilibrium, the tube needs to be in a neutrally buoyant state with vanishing buoyancy force, and with the magnetic curvature force pointing towards the rotation axis being balanced by a Coriolis force produced by a faster rotational speed of the flux ring (see Figure 5(a)). Such a neutrally buoyant flux ring (with equal density between inside and outside) then requires a lower internal temperature than the surrounding plasma to satisfy the total pressure balance. If one starts with a toroidal flux ring that is initially in thermal equilibrium with the surrounding and rotates at the same ambient angular velocity, then the flux ring will move radially outward due to its buoyancy and latitudinally poleward due to the unbalanced poleward component of the tension force. As a result of its motion, the flux ring will lose buoyancy due to the subadiabatic stratification and attain a larger internal rotation rate with respect to the ambient field-free plasma due to the conservation of angular momentum, evolving towards a mechanical equilibrium configuration. The flux ring will undergo superposed buoyancy and inertial oscillations around this mechanical equilibrium state. It is found that the oscillations can be contained within the stably stratified overshoot layer and also within a latitudinal range of Δθ ≲ 20° to be consistent with the active region belt, if the field strength of the toroidal flux ring B ≲ 105 G and the subadiabaticity of the overshoot layer is sufficiently strong with δ ≡ ∇ − ∇ad ≲ −10−5, where ∇ ≡ d ln T/d ln P is the logarithmic temperature gradient and ∇ad is ∇ for an adiabatically stratified atmosphere. Flux rings with significantly larger field strength cannot be kept within the low latitude zones of the overshoot region.

Figure 5:
figure 5

Schematic illustrations based on Schüssler and Rempel ( 2002) of the various forces involved with the mechanical equilibria of an isolated toroidal flux ring (a) and a magnetic layer (b) at the base of the solar convection zone. In the case of an isolated toroidal ring (see the black dot in (a) indicating the location of the tube cross-section), the buoyancy force has a component parallel to the rotation axis, which cannot be balanced by any other forces. Thus mechanical equilibrium requires that the buoyancy force vanishes and the magnetic curvature force is balanced by the Coriolis force resulting from a prograde toroidal flow in the flux ring. For a magnetic layer (as indicated by the shaded region in (b)), on the other hand, a latitudinal pressure gradient can be built up, so that an equilibrium may also exist where a non-vanishing buoyancy force, the magnetic curvature force and the pressure gradient are in balance with vanishing Coriolis force (vanishing longitudinal flow).

Rempel et al. (2000) considered the mechanical equilibrium of a layer of an axisymmetric toroidal magnetic field of 105 G in a subadiabatically stratified region near the bottom of the solar convection zone in full spherical geometry. In this case, as illustrated in Figure 5(b), a latitudinal pressure gradient can be built up, allowing for force balance between a non-vanishing buoyancy force, the magnetic curvature force, and the pressure gradient without requiring a prograde toroidal flow. Thus a wider range of equilibria can exist. Rempel et al. (2000) found that under the condition of a strong subadiabatic stratification such as the radiative interior with δ ∼ − 0.1, the magnetic layer tends to establish a mechanical equilibrium where a latitudinal pressure gradient is built up to balance the poleward component of the magnetic tension, and where the net radial component of the buoyancy and magnetic tension forces is efficiently balanced by the strong subadiabaticity. The magnetic layer reaches this equilibrium solution in a time scale short compared to the time required for a prograde toroidal flow to set up for the Coriolis force to be significant. For this type of equilibrium where a latitudinal pressure gradient is playing a dominant role in balancing the poleward component of the magnetic curvature force, there is significant relative density perturbation (≫ 1/β) in the magnetic layer compared to the background stratification. On the other hand, under the condition of a very weak subadiabatic stratification such as that in the overshoot layer near the bottom of the convection zone with δ ∼ −10−5, the magnetic layer tends to evolve towards a mechanical equilibrium which resembles that of an isolated toroidal flux ring, where the relative density perturbation is small (≪C 1/β), and the magnetic curvature force is balanced by the Coriolis force induced by a prograde toroidal flow in the magnetic layer. Thus regardless of whether the field is in the state of an extended magnetic layer or isolated flux tubes, a 105 G toroidal magnetic field stored in the weakly subadiabatically stratified overshoot region is preferably in a mechanical equilibrium with small relative density perturbation and with a prograde toroidal flow whose Coriolis force balances the magnetic tension. The prograde toroidal flow necessary for the equilibrium of the 105 G toroidal field is about 200 ms−1, which is approximately 10% of the mean rotation rate of the Sun. Thus one may expect significant changes in the differential rotation in the overshoot region during the solar cycle as the toroidal field is being amplified (Rempel et al., 2000). Detecting these toroidal flows and their temporal variation in the overshoot layer via helioseismic techniques is a means by which we can probe and measure the toroidal magnetic field generated by the solar cycle dynamo.

3.2 Effect of radiative heating

Storage of a strong super-equipartition field of 105 G at the base of the solar convection zone requires a state of mechanical equilibrium since convective motion is not strong enough to counteract the magnetic stress (Section 5.5). For isolated flux tubes stored in the weakly subadiabatic overshoot layer, the mechanical equilibrium corresponds to a neutrally buoyant state with a lower internal temperature (Section 3.1). Therefore flux tubes will be heated by radiative diffusion due to the mean temperature difference between the tube and the surrounding field-free plasma (see Parker, 1979; van Ballegooijen, 1982). Moreover, it is not adequate to just consider this zeroth order contribution due to the mean temperature difference in evaluating the radiative heat exchange between the flux tube and its surroundings. Due to the convective heat transport, the temperature gradient in the overshoot region and the lower convection zone is very close to being adiabatic, deviating significantly from that of a radiative equilibrium, and hence there is a non-zero divergence of radiative heat flux (see Spruit, 1974; van Ballegooijen, 1982). Thus an isolated magnetic flux tube with internally suppressed convective transport should also experience a net heating due to this non-zero divergence of radiative heat flux, provided that the radiative diffusion is approximately unaffected within the flux tube (Fan and Fisher, 1996; Moreno-Insertis et al., 2002; Rempel, 2003). In the limit of a thin flux tube, the rate of radiative heating (per unit volume) experienced by the tube is estimated to be (Fan and Fisher, 1996)

$${{dQ} \over {dt}} = - \nabla \cdot ({{\bf{F}}_{{\rm{rad}}}}) - \kappa (x_1^2/{a^2})(\overline T - {\overline T _{\rm{e}}})$$

where F rad is the unperturbed radiative energy flux, κ is the unperturbed radiative conductivity, x 1 is the first zero of the Bessel function J 0(x), a is the tube radius, \(\overline{T}\) is the mean temperature of the flux tube, and \(\overline{T}_{\rm e}\) is the corresponding unperturbed temperature at the location of the tube. Under the conditions prevailing near the base of the solar convection zone and for flux tubes that are responsible for active region formation, the first term due to the non-vanishing divergence of the radiative heat flux is found in general to dominate the second term. In the overshoot region, it can be shown that for these flux tubes the time scale for the heating to significantly increase their buoyancy from an initial neutrally buoyant state is long compared to the dynamic time scale characterized by the Brunt-Väisälä frequency. Thus the radiative heating is found to cause a quasi-static rise of the toroidal flux tubes, during which the tubes remain close to being neutrally buoyant. The upward drift velocity is estimated to be ∼ 10−3|δ|−1 cms−1 which does not depend sensitively on the field strength of the flux tube (Fan and Fisher, 1996; Rempel, 2003). This implies that maintaining toroidal flux tubes in the overshoot region for a period comparable to the solar cycle time scale requires a strong subadiabaticity of δ < −10−4, which is significantly more subadiabatic than the values obtained by most of the overshoot models based on the non-local mixing length theory (see van Ballegooijen, 1982; Schmitt et al., 1984; Skaley and Stix, 1991).

On the other hand if the spatial filling factor of the toroidal flux tubes is large, or if the toroidal magnetic field is stored in the form of an extended magnetic layer, then the suppression of convective motion by the magnetic field is expected to alter the overall temperature stratification in the overshoot region. Rempel (2003) performed a 1D thermal diffusion calculation to model the change of the mean temperature stratification in the overshoot region when convective heat transport is being significantly suppressed. It is found that a reduction of the convective heat conductivity by a factor of 100 leads to the establishment of a new thermal equilibrium of significantly more stable temperature stratification with δ ∼ −10−4 in a time scale of a few months. Thus as the toroidal magnetic field is being amplified by the solar dynamo process, it may improve the conditions for its own storage by reducing the convective energy transport and increasing the subadiabaticity in the overshoot region.

4 Destabilization of a Toroidal Magnetic Field and Formation of Buoyant Flux Tubes

In the previous section, we have reviewed the equilibrium properties of a strong (∼ 105 G) toroidal magnetic field stored at the base of the solar convection zone. In this section we focus on the stability of the equilibria and the mechanisms by which the magnetic field can escape in the form of discrete buoyant flux tubes.

4.1 The buoyancy instability of isolated toroidal magnetic flux tubes

By linearizing the thin flux tube dynamic equations (1), (2), (5), (6), and (7), the stability of neutrally buoyant toroidal magnetic flux tubes to isentropic perturbations have been studied (see Spruit and van Ballegooijen, 1982a,b; Ferriz-Mas and Schüssler, 1993, 1995).

In the simplified case of a horizontal neutrally buoyant flux tube in a plane parallel atmosphere, ignoring the effects of curvature and solar rotation, the necessary and sufficient condition for instability is (Spruit and van Ballegooijen, 1982a,b)

$${k^2}H_p^2 < {{\beta/2} \over {1 + \beta }}(1/\gamma + \beta \delta),$$
(17)

where k is the wavenumber along the tube of the undulatory perturbation, H p is the local pressure scale height, βp/(B 2/8π) is the ratio of the plasma pressure divided by the magnetic pressure of the flux tube, δ = ∂ − ∂ad is the superadiabaticity, and γ is the ratio of the specific heats. If all values of k are allowed, then the condition for the presence of instability is

$$\beta \delta > - 1/\gamma.$$
(18)

Note that k → 0 is a singular limit. For perturbations with k = 0 which do not involve bending the field lines, the condition for instability becomes (Spruit and van Ballegooijen, 1982a)

$$\beta \delta > {2 \over \gamma }\left({{1 \over \gamma } - {1 \over 2}} \right)\sim 0.12$$
(19)

which is a significantly more stringent condition than (18), even more stringent than the convective instability for a field-free fluid (δ > 0). Thus the undulatory instability (with k ≠ 0) is of a very different nature and is easier to develop than the instability associated with uniform up-and-down motions of the entire flux tube. The undulatory instability can develop even in a convectively stable stratification with δ ∼ 0 as long as the field strength of the flux tube is sufficiently strong (i.e. β is of sufficiently small amplitude) such that |βδ| is smaller than 1/γ. In the regime of − 1/γ < βδ < (2/γ)(1/γ − 1/2) where only the undulatory modes with k ≠ 0 are unstable, a longitudinal flow from the crests to the troughs of the undulation is essential for driving the instability. Since the flux tube has a lower internal temperature and hence a smaller pressure scale height inside, upon bending the tube, matter will flow from the crests to the troughs to establish hydrostatic equilibrium along the field. This increases the buoyancy of the crests and destabilizes the tube (Spruit and van Ballegooijen, 1982a).

Including the curvature effect of spherical geometry, but still ignoring solar rotation, Spruit and van Ballegooijen (1982a,b) have also studied the special case of a toroidal flux ring in mechanical equilibrium within the equatorial plane. Since the Coriolis force due to solar rotation is ignored, the flux ring in the equatorial plane needs to be slightly buoyant to balance the inward tension force. For latitudinal motions out of the equatorial plane, the axisymmetric component is unstable, which corresponds to the poleward slip of the tube as a whole. But this instability can be suppressed when the Coriolis force is included (Ferriz-Mas and Schüssler, 1993). For motions within the equatorial plane, the conditions for instabilities are (Spruit and van Ballegooijen, 1982a,b)

$$\matrix{{{1 \over 2}\beta \delta > ({m^2} - 3 - s){f^2} + 2f/\gamma - 1(2\gamma)} \hfill \& {(m \geq 1),} \hfill \cr {{1 \over 2}\beta \delta > {f^2}(1 - s) - 2f/\gamma + {1 \over \gamma }\left({{1 \over \gamma } - {1 \over 2}} \right)} \hfill \& {(m = 0)} \hfill \cr}$$
(20)

where fH p/r 0 is the ratio of the pressure scale height over the radius of the bottom of the solar convection zone, m (having integer values 0, 1, …) denotes the azimuthal order of the undulatory mode of the closed toroidal flux ring, i.e. the wavenumber k = m/r 0, s is a parameter that describes the variation of the gravitational acceleration: gr s. Near the base of the solar convection zone, f ∼ 0.1, s ∼ −2. Thus conditions (20) show that it is possible for m = 0, 1, 2, 3, 4 modes to become unstable in the weakly subadiabatic overshoot region, and that the instabilities of m = 1, 2, 3 modes require less stringent conditions than the instability of m = 0 mode. Since Equation (20) is derived for the singular case of an equilibrium toroidal ring in the equatorial plane, its applicability is very limited.

The general problem of the linear stability of a thin toroidal flux ring in mechanical equilibrium in a differentially rotating spherical convection zone at arbitrary latitudes has been studied in detail by Ferriz-Mas and Schüssler (1993, 1995). For general non-axisymmetric perturbations, a sixth-order dispersion relation is obtained from the linearized thin flux tube equations. It is not possible to obtain analytical stability criteria. The dispersion relation is solved numerically to find instability and the growth rates of the unstable modes. The regions of instability in the (B 0, λ 0) plane (with B 0 being the magnetic field strength of the flux ring and λ 0 being the equilibrium latitude), under the conditions representative of the overshoot layer at the base of the solar convection zone are shown in Figure 6 (from Caligari et al., 1995). The basic parameters that determine the stability of an equilibrium toroidal flux ring are its field strength and the subadiabaticity of the external stratification. In the case δ ≡ ∇ − ∇ad = −2.6 × 10−6 (upper panel of Figure 6), unstable modes with reasonably short growth times (less than about a year) only begin to appear at sunspot latitudes for B 0 ≳ 1.2 × 105 G. These unstable modes are of m = 1 and 2. In case of a weaker subadiabaticity, δ ≡ ∇ − ∇ad = −1.9 × 10−7 (lower panel of Figure 6), reasonably fast growing modes (growth time less than a year) begin to appear at sunspot latitudes for B 0 ≳ 5 × 104 G, and the most unstable modes are of m = 1 and 2. These results suggest that toroidal magnetic fields stored in the overshoot layer at the base of the solar convection zone do not become unstable until their field strength becomes significantly greater than the equipartition value of 104 G.

Figure 6:
figure 6

From Caligari et al. (1995). Upper panel: Regions of unstable toroidal flux tubes in the (B 0, λ 0)-plane (with B 0being the magnetic field strength of the flux tubes and λ 0being the equilibrium latitude). The subadiabaticity at the location of the toroidal flux tubes is assumed to be δ ≡ ∇ − ∇ad = −2.6 × 10−6. The white area corresponds to a stable region while the shaded regions indicate instability. The degree of shading signifies the azimuthal wavenumber of the most unstable mode. The contours correspond to lines of constant growth time of the instability. Thicker lines are drawn for growth times of 100 days and 300 days. Lower panel: Same as the upper panel except that the subadiabaticity at the location of the toroidal tubes is δ ≡ ∇ − ∇ad = −1.9 × 10−7.

Thin flux tube simulations of the non-linear growth of the non-axisymmetric instabilities of initially toroidal flux tubes and the emergence of Ω-shaped flux loops through the solar convective envelope will be discussed in Section 5.1.

4.2 Breakup of an equilibrium magnetic layer and formation of buoyant flux tubes

It is possible that the toroidal magnetic field stored at the base of the convection zone is in the form of an extended magnetic layer, instead of individual magnetic flux tubes for which the thin flux tube approximation can be applied. The classic problem of the buoyancy instability of a horizontal magnetic field \({\bf B}=B(z)\hat{{\bf x}}\) in a plane-parallel, gravitationally stratified atmosphere with a constant gravity \(-g\hat{\bf z}\), pressure p (z), and density ρ (z), in hydrostatic equilibrium,

$${d \over {dz}}\left({p + {{{B^2}} \over {8\pi }}} \right) = - \rho g,$$
(21)

has been studied by many authors in a broad range of astrophysics contexts including

  • magnetic fields in stellar convection zones (see Newcomb, 1961; Parker, 1979; Hughes and Cattaneo, 1987),

  • magnetic flux emergence into the solar atmosphere (see Shibata et al., 1989),

  • stability of prominence support by a magnetic field (see Zweibel and Bruhwiler, 1992),

  • and the instability of the interstellar gas and magnetic field (see Parker, 1966).

The linear stability analysis of the above equilibrium horizontal magnetic layer (Newcomb, 1961) showed that the necessary and sufficient condition for the onset of the general 3D instability with non-zero wavenumbers (k x ≠ 0, k y ≠ 0) in both horizontal directions parallel and perpendicular to the magnetic field is that

$${{d\rho} \over {dz}} > - {{{\rho ^2}g} \over {\gamma p}},$$
(22)

is satisfied somewhere in the stratified fluid. On the other hand the necessary and sufficient condition for instability of the purely interchange modes (with k x = 0 and k y ≠ 0) is that

$${{d \rho} \over {dz}} > - {{{\rho ^2}g} \over {\gamma p + {B^2}/4\pi}}.$$
(23)

is satisfied somewhere in the fluid — a more stringent condition than (22). Note in Equations (22) and (23), p and ρ are the plasma pressure and density in the presence of the magnetic field. Hence the effect of the magnetic field on the instability criteria is implicitly included. As shown by Thomas and Nye (1975) and Acheson (1979), the instability conditions (22) and (23) can be alternatively written as

$${{v_{\rm{a}}^2} \over {c_{\rm{s}}^2}}{{d\ln B} \over {dz}} < - {1 \over {{c_p}}}{{ds} \over {dz}}$$
(24)

for instability of general 3D undulatory modes and

$${{v_{\rm{a}}^2} \over {c_{\rm{s}}^2}}{d \over {dz}}\left[ {\ln \left({{B \over \rho }} \right)} \right] < - {1 \over {{c_p}}}{{ds} \over {dz}}$$
(25)

for instability of purely 2D interchange modes, where v a is the Alfvén speed, c s is the sound speed, c p is the specific heat under constant pressure, and ds/dz is the actual entropy gradient in the presence of the magnetic field. The development of these buoyancy instabilities is driven by the gravitational potential energy that is made available by the magnetic pressure support. For example, the magnetic pressure gradient can “puff-up” the density stratification in the atmosphere, making it decrease less steeply with height (causing condition (22) to be met), or even making it top heavy. This raises the gravitational potential energy and makes the atmosphere unstable. In another situation, the presence of the magnetic pressure can support a layer of cooler plasma with locally reduced temperature embedded in an otherwise stably stratified fluid. This can also cause the instability condition (22) to be met locally in the magnetic layer. In this case the pressure scale height within the cooler magnetic layer is smaller, and upon bending the field lines, plasma will flow from the crests to the troughs to establish hydrostatic equilibrium, thereby releasing gravitational potential energy and driving the instability. This situation is very similar to the buoyancy instability associated with the neutrally buoyant magnetic flux tubes discussed in Section 4.1.

The above discussion on the buoyancy instabilities considers ideal adiabatic perturbations. It should be noted that the role of finite diffusion is not always stabilizing. In the solar interior, it is expected that ηK and νK, where η, ν, and K denote the magnetic diffusivity, the kinematic viscosity, and the thermal diffusivity respectively. Under these circumstances, it is shown that thermal diffusion can be destabilizing (see Gilman, 1970; Acheson, 1979; Schmitt and Rosner, 1983). The diffusive effects are shown to alter the stability criteria of Equations (24) and (25) by reducing the term ds/dz by a factor of η/K (see Acheson, 1979). In other words, efficient heat exchange can significantly “erode away” the stabilizing effect of a subadiabatic stratification. This process is called the doubly diffusive instabilities.

Direct multi-dimensional MHD simulations have been carried out to study the break-up of a horizontal magnetic layer by the non-linear evolution of the buoyancy instabilities and the formation of buoyant magnetic flux tubes (see Cattaneo and Hughes, 1988; Cattaneo et al., 1990; Matthews et al., 1995; Wissink et al., 2000; Fan, 2001a).

Cattaneo and Hughes (1988), Matthews et al. (1995), and Wissink et al. (2000) have carried out a series of 2D and 3D compressible MHD simulations where they considered an initial horizontal magnetic layer that supports a top-heavy density gradient, i.e. an equilibrium with a lower density magnetic layer supporting a denser plasma on top of it. It is found that for this equilibrium configuration, the most unstable modes are the Rayleigh-Taylor type 2D interchange modes. Two-dimensional simulations of the non-linear growth of the interchange modes (Cattaneo and Hughes, 1988) found that the formation of buoyant flux tubes is accompanied by the development of strong vortices whose interactions rapidly destroy the coherence of the flux tubes. In the non-linear regime, the evolution is dominated by vortex interactions which act to prevent the rise of the buoyant magnetic field. Matthews et al. (1995) and Wissink et al. (2000) extend the simulations of Cattaneo and Hughes (1988) to 3D allowing variations in the direction of the initial magnetic field. They discovered that the flux tubes formed by the initial growth of the 2D interchange modes subsequently become unstable to a 3D undulatory motion in the non-linear regime due to the interaction between neighboring counter-rotating vortex tubes, and consequently the flux tubes become arched. Matthews et al. (1995) and Wissink et al. (2000) pointed out that this secondary undulatory instability found in the simulations is of similar nature as the undulatory instability of a pair of counter-rotating (non-magnetic) line vortices investigated by Crow (1970). Wissink et al. (2000) further considered the effect of the Coriolis force due to solar rotation using a local f-plane approximation, and found that the principle effect of the Coriolis force is to suppress the instability. Further 2D simulations have also been carried out by Cattaneo et al. (1990) where they introduced a variation of the magnetic field direction with height into the previously unidirectional magnetic layer of Cattaneo and Hughes (1988). The growth of the interchange instability of such a sheared magnetic layer results in the formation of twisted, buoyant flux tubes which are able to inhibit the development of vortex tubes and rise cohesively.

On the other hand, Fan (2001a) has considered a different initial equilibrium state for a horizontal unidirectional magnetic layer, where the density stratification remains unchanged from that of an adiabatically stratified polytrope, but the temperature and the gas pressure are lowered in the magnetic layer to satisfy the hydrostatic condition. For such a neutrally buoyant state with no density change inside the magnetic layer, the 2D interchange instability is completely suppressed and only 3D undulatory modes (with non-zero wavenumbers in the field direction) are unstable. The strong toroidal magnetic field stored in the weakly subadiabatic overshoot region below the bottom of the convection zone is likely to be close to such a neutrally buoyant mechanical equilibrium state (see Section 3.1). Anelastic MHD simulations (Fan, 2001a) of the growth of the 3D undulatory instability of this horizontal magnetic layer show formation of significantly arched magnetic flux tubes (see Figure 7) whose apices become increasingly buoyant as a result of the diverging flow of plasma from the apices to the troughs. The decrease of the field strength B at the apex of the arched flux tube as a function of height is found to follow approximately the relation \(B/\sqrt \rho = {\rm{constant}}\), or, the Alfvén speed being constant, which is a significantly slower decrease of B with height compared to that for the rise of a horizontal flux tube without any field line stretching, for which case B/ρ should remain constant. The variation of the apex field strength with height following \(B/\sqrt \rho = {\rm{constant}}\) found in the 3D MHD simulations of the arched flux tubes is in good agreement with the results of the thin flux tube models of emerging Ω-loops (see Moreno-Insertis, 1992) during their rise through the lower half of the solar convective envelope where the stratification is very close to being adiabatic as is assumed in the 3D simulations.

Figure 7:
figure 7

mpg-Movie (255 KB) Still from a movie showing The formation of arched flux tubes as a result of the non-linear growth of the undulatory buoyancy instability of a neutrally buoyant equilibrium magnetic layer perturbed by a localized velocity field. From Fan ( 2001a). The images show the volume rendering of the absolute magnetic field strength |B|. Only one half of the wave length of the undulating flux tubes is shown, and the left and right columns of images show, respectively, the 3D evolution as viewed from two different angles. (For video see appendix)

5 Dynamic Evolution of Emerging Flux Tubes in the Solar Convection Zone

5.1 Results from thin flux tube simulations of emerging loops

Beginning with the seminal work of Moreno-Insertis (1986) and Choudhuri and Gilman (1987), a large body of numerical simulations solving the thin flux tube dynamic equations (1), (2), (5), (6), and (7) — or various simplified versions of them — have been carried out to model the evolution of emerging magnetic flux tubes in the solar convective envelope (see Choudhuri, 1989; D’Silva and Choudhuri, 1993; Fan et al., 1993, 1994; Schüssler et al., 1994; Caligari et al., 1995; Fan and Fisher, 1996; Caligari et al., 1998; Fan and Gong, 2000). The results of these numerical calculations have contributed greatly to our understanding of the basic properties of solar active regions and provided constraints on the field strengths of the toroidal magnetic fields at the base of the solar convection zone.

Most of the earlier calculations (see Choudhuri and Gilman, 1987; Choudhuri, 1989; D’Silva and Choudhuri, 1993; Fan et al., 1993, 1994) considered initially buoyant toroidal flux tubes by assuming that they are in temperature equilibrium with the external plasma. Various types of initial undulatory displacements are imposed on the buoyant tube so that portions of the tube will remain anchored within the stably stratified overshoot layer and other portions of the tube are displaced into the unstable convection zone which subsequently develop into emerging Ω-shaped loops.

Later calculations (see Schüssler et al., 1994; Caligari et al., 1995, 1998; Fan and Gong, 2000) considered more physically self-consistent initial conditions where the initial toroidal flux ring is in the state of mechanical equilibrium. In this state the buoyancy force is zero (neutrally buoyant) and the magnetic curvature force is balanced by the Coriolis force resulting from a prograde toroidal motion of the tube plasma. It is argued that this mechanical equilibrium state is the preferred state for the long-term storage of a toroidal magnetic field in the stably stratified overshoot region (Section 3.1). In these simulations, the development of the emerging Ω-loops is obtained naturally by the non-linear, adiabatic growth of the undulatory buoyancy instability associated with the initial equilibrium toroidal flux rings (Section 4.1). As a result there is far less degree of freedom in specifying the initial perturbations. The eruption pattern needs not be prescribed in an ad hoc fashion but is self-consistently determined by the growth of the instability once the initial field strength, latitude, and the subadiabaticity at the depth of the tube are given. For example Caligari et al. (1995) modeled emerging loops developed due to the undulatory buoyancy instability of initial toroidal flux tubes located at different depths near the base of their model solar convection zone which includes a consistently calculated overshoot layer according to the non-local mixing-length treatment. They choose values of initial field strengths and latitudes that lie along the contours of constant instability growth times of 100 days and 300 days in the instability diagrams (see Figure 6), given the subadiabaticity at the depth of the initial tubes. The tubes are then perturbed with a small undulatory displacement which consists of a random superposition of Fourier modes with azimuthal order ranging from m = 1 through m = 5, and the resulting eruption pattern is determined naturally by the growth of the instability.

On the other hand, non-adiabatic effects may also be important in the destabilization process. It has been discussed in Section 3.2 that isolated magnetic flux tubes with internally suppressed convective transport experience a net heating due to the non-zero divergence of radiative heat flux in the weakly subadiabatically stratified overshoot region and also in the lower solar convection zone. The radiative heating causes a quasi-static upward drift of the toroidal flux tube with a drift velocity ∼ 10−3|δ|−1 cms−1. Thus the time scale for a toroidal flux tube to drift out of the stable overshoot region may not be long compared to the growth time of its undulatory buoyancy instability. For example if the subadiabaticity δ is ∼ −10−6, the time scale for the flux tube to drift across the depth of the overshoot region is about 20 days, smaller than the growth times (∼ 100–300 days) of the most unstable modes for tubes of a ∼ 105 G field as shown in Figure 6. Therefore radiative heating may play an important role in destabilizing the toroidal flux tubes. The quasi-static upward drift due to radiative heating can speed-up the development of emerging Ω-loops (especially for weaker flux tubes) by bringing the tube out of the inner part of the overshoot region of stronger subadiabaticity, where the tube is stable or the instability growth is very slow, to the outer overshoot region of weaker subadiabaticity or even into the convection zone, where the growth of the undulatory buoyancy instability occurs at a much shorter time scale.

A possible scenario in which the effect of radiative heating helps to induce the formation of Ω-shaped emerging loops has been investigated by Fan and Fisher (1996). In this scenario, the initial neutrally buoyant toroidal flux tube is not exactly uniform, and lies at non-uniform depths with some portions of the tube lying at slightly shallower depths in the overshoot region. Radiative heating and quasi-static upward drift of this non-uniform flux tube bring the upward protruding portions of the tube first into the unstably stratified convection zone. These portions can become buoyantly unstable (if the growth of buoyancy overcomes the growth of tension) and rise dynamically as emerging loops. In this case the non-uniform flux tube remains close to a mechanical equilibrium state during the initial quasi-static rise through the overshoot region. The emerging loop develops gradually as a result of radiative heating and the subsequent buoyancy instability of the outer portion of the tube entering the convection zone.

In the following subsections we review the major findings and conclusions that have been drawn from the various thin flux tube simulations of emerging flux loops.

5.1.1 Latitude of flux emergence

Axisymmetric simulations of the buoyant rise of toroidal flux rings in a rotating solar convective envelope by Choudhuri and Gilman (1987) first demonstrate the significant influence of the Coriolis force on the rising trajectories. The basic effect is that the Coriolis force acting on the radial outward motion of the flux tube (or the tendency for the rising tube to conserve angular momentum) drives a retrograde motion of the tube plasma. This retrograde motion then induces a Coriolis force directed towards the Sun’s rotation axis which acts to deflect the trajectory of the rising tube poleward. The amount of poleward deflection by the Coriolis force depends on the initial field strength of the emerging tube, being larger for flux tubes with weaker initial field. For flux tubes with an equipartition field strength of B ∼ 104 G, the effect of the Coriolis force is so dominating that it deflects the rising tubes to emerge at latitudes poleward of the sunspot zones even though the flux tubes start out from low latitudes at the bottom of the convective envelope. In order for the rising trajectory of the flux ring to be close to radial so that the emerging latitudes are within the observed sunspot latitudes, the field strength of the toroidal flux ring at the bottom of the solar convection zone needs to be ∼ 105 G. This basic result is confirmed by later simulations of non-axisymmetric, Ω-shaped emerging loops rising through the solar convective envelope (see Choudhuri, 1989; D’Silva and Choudhuri, 1993; Fan et al., 1993; Schüssler et al., 1994; Caligari et al., 1995, 1998; Fan and Fisher, 1996).

Simulations by Caligari et al. (1995) of emerging loops developed self-consistently due to the undulatory buoyancy instability show that, for tubes with initial field strength ≳ 105 G, the trajectories of the emerging loops are primarily radial with poleward deflection no greater than 3°. For tubes with initial field strength exceeding 4 × 104 G, the poleward deflection of the emerging loops remain reasonably small (no greater than about 6°). However, for a tube with equipartition field strength of 104 G, the rising trajectory of the emerging loop is deflected poleward by about 20°. Such an amount of poleward deflection is too great to explain the observed low latitudes of active region emergence. Furthermore, it is found that with such a weak initial field the field strength of the emerging loop falls below equipartition with convection throughout most of the convection zone. Such emerging loops are expected to be subjected to strong deformation by turbulent convection and may not be consistent with the observed well defined order of solar active regions.

Fan and Fisher (1996) modeled emerging loops that develop as a result of radiative heating of non-uniform flux tubes in the overshoot region. The results on the poleward deflection of the emerging loops as a function of the initial field strength are very similar to that found in Caligari et al. (1995). Figure 8 shows the latitude of loop emergence as a function of the initial latitude at the base of the solar convection zone. It can be seen that tubes of 105 G emerge essentially radially with very small poleward deflection (< 3°), and for tubes with B ≳ 3 × 104 G, the poleward deflections remain reasonably small so that the emerging latitudes are within the observed sunspot zones.

Figure 8:
figure 8

From Fan and Fisher ( 1996). Latitude of loop emergence as a function of the initial latitude at the base of the solar convection zone, for tubes with initial field strengths B = 30, 60, and 100 kG and fluxes Φ = 1021and 1022 Mx.

5.1.2 Active region tilts

A well-known property of the solar active regions is the so called Joy’s law of active region tilts. The averaged orientation of bipolar active regions on the solar surface is not exactly toroidal but is slightly tilted away from the east-west direction, with the leading polarity (the polarity leading in the direction of rotation) being slightly closer to the equator than the following polarity. The mean tilt angle is a function of latitude, being approximately ∝ sin(latitude) (Wang and Sheeley Jr, 1989, 1991; Howard, 1991a,b; Fisher et al., 1995).

Using thin flux tube simulations of the non-axisymmetric eruption of buoyant Ω-loops in a rotating solar convective envelope, D’Silva and Choudhuri (1993) are the first to show that the active region tilts as described by Joy’s law can be explained by Coriolis forces acting on the flux loops. As the emerging loop rises, there is a relative expanding motion of the mass elements at the summit of the loop. The Coriolis force induced by this diverging, expanding motion at the summit is to tilt the summit clockwise (counter-clockwise) for loops in the northern (southern) hemisphere as viewed from the top, so that the leading side from the summit is tilted equatorward relative to the following side. Since the component of the Coriolis force that drives this tilting has a sin(latitude) dependence, the resulting tilt angle at the apex is approximately ∝ sin(latitude).

Caligari et al. (1995) studied tilt angles of emerging loops developed self-consistently due to the undulatory buoyancy instability of flux tubes located at the bottom as well as just above the top of their model overshoot region, with selected values of initial field strengths and latitudes lying along contours of constant instability growth times (100 days and 300 days). The resulting tilt angles at the apex of the emerging loops (see Figure 9) produced by these sets of unstable tubes (whose field strengths are within the range of 4 × 104 G to 1.5 × 105 G) show good agreement with the observed tilt angles for sunspot groups measured by Howard (1991b). They also found that loops formed from toroidal flux tubes with an equipartition field strength of 104 G develop a tilt angle (−9°) of the wrong sign at the loop apex.

Figure 9:
figure 9

From Caligari et al. ( 1995). Tilt angles at the apex of the emerging flux loops as a function of the emergence latitudes. The squares and the asterisks denote loops originating from initial toroidal tubes located at different depths with different local subadiabaticity (squares: δ ≡ ∇ − ∇ad = −2.6 × 10−6and field strength ranges between 105 G and 1.5 × 105 G; asterisks: δ ≡ ∇ − ∇ad = −1.9 × 10−7and field strength ranges between 4 × 104 G and 6 × 104 G). The shaded region indicates the range of the observed tilt angles of sunspot groups measured by Howard( 1991b).

Similar results of loop tilt angles (see Figure 10) are found by Fan and Fisher (1996) who considered formation of emerging loops by gradual radiative heating of non-uniform toroidal flux tubes initially in mechanical equilibrium in the overshoot region. It is found that emerging loops with initial field strength of the range 4 × 104 G to 105 G show tilt angles that are in good agreement with the observed tilt angles of sunspot groups. Tilts of the wrong sign or direction begin to appear for tubes with initial field strength ≲ 3 × 104 G.

Figure 10:
figure 10

Tilt angles at the apex of the emerging loops versus emerging latitudes resulting from thin flux tube simulations of Fan and Fisher ( 1996). The numbers used as data points indicate the corresponding initial field strength values in units of 10 kG (’X’s represent 100 kG). Also plotted (solid line) is the least squares fit: tilt angle = 15.7° × sin(latitude), obtained in Fisher et al. ( 1995) by fitting to the measured tilt angles of 24701 sunspot groups observed at Mt. Wilson, same data set studied in Howard ( 1991b).

Thin flux tube simulations also show that the tilt angles of the emerging loops tend to be smaller for tubes with smaller total flux or radius (see D’Silva and Choudhuri, 1993; Fan and Fisher, 1996). This is because smaller tubes are more influenced by the drag force, whose direct effect is to oppose the tilting motion of the loops. This prediction on the dependence of the tilt angle on active region flux is confirmed by Fisher et al. (1995), who studied the relation between the tilt angle and the mean separation distance between the leader and follower spots of a sunspot group. Wang and Sheeley Jr (1989) and Howard (1992) have shown that the mean polarity separation distance of an active region is a good proxy for the total magnetic flux in the region. Fisher et al. (1995) found that for a fixed latitude, the tilt angle of sunspot groups decreases with decreasing polarity separation, and hence decreasing total flux, consistent with the results from the thin flux tube calculations. This agreement adds support to the explanation that the Coriolis force acting on rising flux loops is the main cause of active region tilts and argues against the suggestion by Babcock (1961) that the active region tilt angles simply reflect the orientation of the underlying toroidal magnetic field stretched out by the latitudinal differential rotation.

Using Mount Wilson sunspot group data, Fisher et al. (1995) further studied the dispersion or scatter of spot-group tilts away from the mean tilt behavior described by Joy’s law. First they found that the magnitude of the tilt dispersion is significantly greater than the level expected from measurement errors, suggesting a solar origin of the tilt angle scatter. Furthermore, they found that the root-mean-square tilt scatter decreases with increasing polarity separation (or total flux) and does not vary with latitude (in contrast to the latitudinal dependence of the mean tilts). This result is consistent with the picture that scatter of active region tilts away from the mean Joy’s law behavior results from buffeting of emerging loops by convective motions during their rise through the solar convection zone (Longcope and Fisher, 1996).

Recent non-linear simulations of the two-dimensional MHD tachocline (Cally et al., 2003) show that bands of toroidal magnetic fields in the solar tachocline may become tipped relative to the azimuthal direction by an amount that is within +/−10° at sunspot latitudes due to the non-linear evolution of the 2D global joint instability of differential rotation and toroidal magnetic fields. This tipping may either enhance or reduce the observed tilt in bipolar active regions depending on from which part of the tipped band the emerging loops develop. Thus the basic consequence of the possible tipping of the toroidal magnetic fields in the tachocline is to contribute to the spread of the tilts of bipolar active regions.

5.1.3 Morphological asymmetries of active regions

An intriguing property of solar active regions is the asymmetry in morphology between the leading and following polarities. The leading polarity of an active region tends to be in the form of large sunspots, whereas the following polarity tends to appear more dispersed and fragmented; moreover, the leading spots tend to be longer lived than the following. Fan et al. (1993) offered an explanation for the origin of this asymmetry. In their thin flux tube simulations of the non-axisymmetric eruption of buoyant Ω-loops through a rotating model solar convective envelope, they found that an asymmetry in the magnetic field strength develops between the leading and following legs of an emerging loop, with the field strength of the leading leg being about 2 times that of the following leg. The field strength asymmetry develops because the Coriolis force, or the tendency for the tube plasma to conserve angular momentum, drives a counter-rotating flow of plasma along the emerging loop, which, in conjunction with the diverging flow of plasma from the apex to the troughs, gives rise to an effective asymmetric stretching of the two legs of the loop with a greater stretching and hence a stronger field strength in the leading leg. Fan et al. (1993) argued that the stronger field along the leading leg of the emerging loop makes it less subject to deformation by the turbulent convection and therefore explains the more coherent and less fragmented appearance of the leading polarity of an active region.

However, the calculations by Caligari et al. (1995) and Caligari et al. (1998) show that the field strength asymmetry found in Fan et al. (1993) depends on the choices of the initial conditions. They found that if one uses the more self-consistent mechanical equilibrium state (instead of the buoyant, temperature equilibrium state used in Fan et al. (1993)) for the initial toroidal tube, for which there exists an initial prograde toroidal motion, the subsequent differential stretching of the emerging loop is quite different. It is found that a consistently stronger field along the leading leg of the emerging loop only occurs for cases with relatively weak initial field strengths (∼ 104 G). For stronger fields (∼ 105 G), which seem to fit the observed properties of solar active regions better in all the other aspects, e.g. the latitude of flux emergence and the tilt angles, the field strength asymmetry becomes very small and even reverses its sense in the upper half of the solar convection zone.

Further work by Fan and Fisher (1996) examined the field strength asymmetry for emerging flux loops that form as a result of gradual radiative heating of non-uniform flux tubes initially in mechanical equilibrium. They find (see Figure 11) significantly stronger fields in the leading leg compared to the following in the lower solar convection zone for all initial field strengths, but nearly equal field strengths of the two legs in the upper convection zone for stronger initial fields (B ≳ 6 × 104 G). Clearly, further investigations, and perhaps fully resolved three-dimensional MHD simulations of emerging flux loops, are needed to understand the origin of the observed morphological asymmetries of the two polarities of solar active regions.

Figure 11:
figure 11

Plots of the magnetic field strength as a function of depth along the emerging loops calculated from the thin flux tube model of Fan and Fisher ( 1996) showing the asymmetry in field strength between the leading leg (solid curve) and the following leg (dash-dotted curve) of each loop. Panels (a), (b), and (c) correspond to the cases with initial toroidal field strengths of 3 × 104 G, 6 × 104 G, and 105 G respectively. The flux Φ = 1022 Mx and the initial latitude θ = 5° are the same for the three cases shown.

5.1.4 Geometrical asymmetry of emerging loops and the asymmetric proper motions of active regions

Another asymmetry in the emerging loop generated by the effect of the Coriolis force is the asymmetry in the east-west inclinations of the two sides of the loop. This asymmetry is first shown in the thin flux tube calculations of Moreno-Insertis et al. (1994) and Caligari et al. (1995) who modeled emerging loops that develop self-consistently as a result of the buoyancy instability of toroidal magnetic flux tubes initially in mechanical equilibrium. Moreno-Insertis et al. (1994) and Caligari et al. (1995) found that as the emerging loop rises, the Coriolis force, or the tendency for the tube to conserve angular momentum, drives a counter-rotating motion of the tube plasma, which causes the summit of the loop to move retrograde relative to the valleys, resulting in an asymmetry in the inclinations of the two legs of the loop with the leading leg being inclined more horizontally with respect to the surface than the following leg. This asymmetry in inclination can be clearly seen in Figure 12, which shows a view from the north pole of an asymmetric emerging loop obtained from a simulation by Caligari et al. (1995). Caligari et al. (1998) and Fan and Fisher (1996) show that the geometrical asymmetry is a robust result for models of emerging loops that originate from initial flux tubes in mechanical equilibrium. Loops with initial field strength ranging from 3 × 104 to 105 G consistently show this asymmetry.

Figure 12:
figure 12

From Caligari et al. ( 1995). A view from the north pole of the configuration of an emerging loop obtained from a thin flux tube simulation of a buoyantly unstable initial toroidal flux tube. The initial field strength is 1.2 × 105 G, and the initial latitude is 15°. Note the strong asymmetry in the east-west inclination of the two sides of the emerging loop.

The observational consequences of this geometric asymmetry are discussed in Moreno-Insertis et al. (1994) and Caligari et al. (1995). The emergence of such an eastward inclined loop is expected to produce apparent asymmetric east-west proper motions of the two polarities of the emerging region, with a more rapid motion of the leading polarity spots away from the emerging region compared to the motion of the following polarity spots. Such asymmetric proper motions are observed in young active regions and sunspot groups (see Chou and Wang, 1987; van Driel-Gesztelyi and Petrovay, 1990; Petrovay et al., 1990). Furthermore, the asymmetry in the inclination of the emerging loop may also explain the observation that the magnetic inversion line in bipolar regions is statistically nearer to the main following spot than to the main proceeding one (van Driel-Gesztelyi and Petrovay, 1990; Petrovay et al., 1990).

5.2 Hemispheric trend of the twist in solar active regions

Vector magnetic field observations of active regions on the photosphere have revealed that on average the solar active regions have a small but statistically significant mean twist that is left-handed in the northern hemisphere and right-handed in the southern hemisphere (see Pevtsov et al., 1995, 2001, 2003). What is being measured is the quantity α ≡ 〈J z/B z〉, the ratio of the vertical electric current over the vertical magnetic field averaged over the active region. When plotted as a function of latitude, the measured α for individual solar active regions show considerable scatter, but there is clearly a statistically significant trend for negative (positive) α in the northern (southern) hemisphere (see Figure 13). A linear least squares fit to the data of α as a function of latitude (Figure 13a) found that α = −2.7 × 10−1° θdeg m−1, where θdeg is latitude in degrees, and that the r.m.s. scatter of a from the linear fit is Δα = 1.28 × 10−8 m−1 (Longcope et al., 1998). The observed systematic a in solar active regions may reflect a systematic field line twist in the subsurface emerging flux tubes.

Figure 13:
figure 13

From Pevtsov et al. ( 2001). The figure shows the latitudinal profile of α best(see Pevtsov et al., 1995, for the exact way of determining α best) for (a) 203 active regions in cycle 22 (Longcope et al., 1998), and (b) 263 active regions in cycle 23. Error bars (when present) correspond to 1 standard deviation of the mean α bestfrom multiple magnetograms of the same active region. Points without error bars correspond to active regions represented by a single magnetogram. The solid line shows a least-squares best-fit linear function.

If the measured α values are a direct consequence of the emergence of twisted magnetic flux tubes from the interior, then it would imply subsurface emerging tubes with a field line twist of q = α/2 (Longcope and Klapper, 1997), where q denotes the angular rate of field line rotation about the axis over a unit axial distance along the tube. Several subsurface mechanisms for producing twist in emerging flux tubes have been proposed (see e.g. the review by Petrovay et al., 2006). The twist may be due to the current helicity in the dynamo generated toroidal magnetic field, from which buoyant flux tubes form at the base of the convection zone (Gilman and Charbonneau, 1999), or it may be acquired during the rise of the flux tubes through the solar convection zone (Longcope et al., 1998; Choudhuri, 2003; Choudhuri et al., 2004; Chatterjee et al., 2006).

Longcope et al. (1998) explain the origin of the observed twist in emerging active region flux tubes as a result of buffeting by the helical turbulence in the solar convection zone during the rise of the tubes. Applying the dynamic model of a weakly twisted thin flux tube (Section 2.1), Longcope et al. (1998) modeled the rise of a nearly straight, initially untwisted tube, buffeted by a random velocity field representative of the turbulent convection in the solar convection zone, which has a nonzero kinetic helicity due to the effect of solar rotation. The kinetic helicity causes helical distortion of the tube axis, which in turn leads to a net twist of the field lines about the axis in the opposite sense within the tube as a consequence of the conservation of magnetic helicity. This process is termed the Σ-effect by Longcope et al. (1998). Quantitative model calculations of Longcope et al. (1998) show that the Σ-effect can explain the hemispheric sign, magnitude, latitude variation, and the r.m.s. dispersion of the observed a of solar active regions. Furthermore, because the aerodynamic drag force acting on flux tubes by the convective flows has an a −2 dependence, where a is the tube radius, this model predicts that the resulting twist generated should also have an a −2 dependence, suggesting a systematic trend for greater twist in smaller flux tubes.

As discussed in Section 5.1.2, Ω-shaped emerging loops are themselves acted upon by the Coriolis force, developing a “tilt” of the loop. This helical deformation of the tube axis will then also induce a twist of the field lines of the opposite sense within the tube as a consequence of conservation of magnetic helicity. Calculations based on the weakly twisted thin flux tube model show that this twist generated by the large scale tilting of the emerging Ω-loop resulting from the Coriolis force has the right hemispheric sign and latitude dependence, but is of too small a magnitude to account for the observed twist in solar active regions (Longcope and Klapper, 1997; Fan and Gong, 2000).

Another interesting and natural explanation for the origin of twist in emerging flux tubes is the accretion of the background mean poloidal field onto the rising flux tube as it traverse through the solar convection zone (Choudhuri, 2003; Choudhuri et al., 2004; Chatterjee et al., 2006). In a Babcock-Leighton type dynamo, the dispersal of solar active regions with a slight mean tilt angle at the surface generates a mean poloidal magnetic field. The mean tilt angle of solar active regions is produced by the Coriolis force acting on the rising flux tubes (see Section 5.1.2). In the northern hemisphere, when a toroidal flux tube rises into a poloidal field that has been created due to the tilt of the same type of flux tubes emerged earlier, the poloidal field gets wrapped around the flux tube will produce a left-handed twist for the tube. This is illustrated in Figure 14. Using a circulation-dominated (or flux-transport) Babcock-Leighton type mean-field dynamo model, Choudhuri et al. (2004) did a rough estimate of the twist acquired by an emerging flux tube rising through the solar convection zone. Figure 15 shows the resulting butterfly diagram indicating the sign of α of the emerging regions as a function of latitude and time.

Figure 14:
figure 14

From Choudhuri ( 2003). This figure illustrates that in the northern hemisphere, when a toroidal flux tube (whose cross-section is the hashed area with a magnetic field going into the paper) rising into a region of poloidal magnetic field (in the clockwise direction) generated by the Babcock-Leighton type α-effect of earlier emerging flux tubes of the same type, the poloidal field gets wrapped around the cross-section of the toroidal tube and reconnects behind it, creating an emerging flux tube with left-handed twist. In this figure, the north-pole is to the left, equator to the right, and the dashed line indicating the solar surface. Note the α-effect for the Babcock-Leighton type solar dynamo model mentioned above is not to be confused with the α value measured in solar active region discussed in this section.

Figure 15:
figure 15

From Choudhuri et al. ( 2004). Simulated butterfly diagram of active region emergence based on a circulation-dominated mean-field dynamo model with Babcock-Leighton α-effect. The sign of the twist of the emerging active region flux tube is determined by considering poloidal flux accretion during its rise through the convection zone. Right handed twist (left handed twist) is indicated by plus signs (circles).

It is found that at the beginning of a solar cycle, there is a short duration where the sign of α is opposite to the preferred sign for the hemisphere. This is because of the phase relation between the toroidal and poloidal magnetic fields produced by this types of solar dynamo models. At the beginning of a cycle, the mean poloidal magnetic field in the convection zone is still dominated by that generated by the emerging flux tubes of the previous cycle, and toroidal flux tubes of the new cycle emerging into this poloidal field gives rise to a right-handed (left-handed) twist of the tube in the northern (southern) hemisphere. However, for the rest of the cycle starting from the solar maximum, the poloidal magnetic field changes sign and the twist for the emerging tubes becomes consistent with the hemispheric preference. For the whole cycle, it is found that about 67% of the emerging regions have a sign of α consistent with the hemispheric rule. The rough estimate also shows that the magnitude of the α values produced by poloidal flux accretion is consistent with the observed values, and that there is an a −2 dependence on the radius a of the emerging tube, i.e. smaller sunspots should have greater α values. This prediction is also made by the Σ-effect mechanism.

Given the frozen-in condition of the magnetic field, it is expected that the accreted poloidal flux be confined in a sheath at the outer periphery of the rising tube, and that in order to produce a twist within the tube, some form of turbulent diffusion needs to be invoked (Chatterjee et al., 2006). By solving the induction equation in a co-moving Lagrangian frame following the rising flux tube and using several simplifying assumptions, Chatterjee et al. (2006) modeled the evolution of the magnetic field in the rising tube cross-section as a result of poloidal flux accretion and penetration due to a field strength dependent turbulent diffusivity. They found that with plausible choices of assumptions and parameter values an α value comparable to the observations is obtained.

When a buoyant magnetic flux tube form at the base of the solar convection zone from the dynamo generated (pre-dominantly) toroidal magnetic field, it should already obtain an initial twist due to the weak poloidal mean field contained in the magnetic layer. This initial twist will then be further augmented or altered due to poloidal flux accretion and also due to the Σ-effect as the tube rises through the solar convection zone. MHD simulations of the formation and rise of buoyant magnetic flux tubes directly incorporating the mean field profiles from dynamo models (for both the fields at the base and in the bulk of the convection zone) as the initial state is necessary to quantify the initial twist and contribution from poloidal field accretion. Such simulations should be done for dynamo mean fields at different phases of the cycle to access the cycle variation of the twist in the emerging flux tubes.

Observationally, looking for any systematic variations of the the α value (or twist) of solar active regions with the solar cycle phase is helpful for identifying the main mechanisms for the origin of the twist. So far results in this area have been inconclusive (Bao et al., 2000; Pevtsov et al., 2001). On the other hand, observational studies of the correlation between active region twist (as measured by α) and tilt angles have revealed interesting results (Holder et al., 2004; Tian et al., 2005; Nandy, 2006). The Σ-effect predicts that the twist being generated in the tube is uncorrelated to the local tilt of the tube at the apex (Longcope et al., 1998). However, due to the Coriolis force, active region Ω-tubes acquire a mean tilt that has a well defined latitudinal dependence as described by the Joy’s law. The mean twist generated by the Σ-effect also has a latitudinal dependence that is consistent with the observed hemispheric rule. Thus, due to the mutual dependence of their mean values on latitude there should be a correlation between the tilt angle and twist of solar active regions and the correlation is expected to be positive if one assigns negative (positive) sign to a clockwise (counter-clockwise) tilt. However, Holder et al. (2004) found a statistically significant negative correlation between the twist and tilt for the 368 bipolar active regions studied, opposite to that expected from the mutual dependence on latitude of mean twist and mean tilt of active regions. Removing the effect of the mutual dependence of the mean tilt and mean twist on latitude (by either determining the correlation at a fixed latitude, or by subtracting off the fitted mean tilt and mean twist at the corresponding latitude), they found that the negative correlation is enhanced. Furthermore, it is found that the negative correlation is mainly contributed by those active regions (174 out of 368 regions) that deviate significantly from Joy’s law (by > 6σ), while regions that obey Joy’s law to within 6σ show no significant correlation between their twist and tilt. A separate study by Tian et al. (2005) found that a sample of 104 complex δ-configuration active regions, more than half of which have tilts that are opposite to the direction prescribed by Joy’s law, show a significant negative correlation between their twist and tilt (after correcting for their definition of the sign of tilt which is opposite to the definition used in Holder et al. (2004)). The results of Holder et al. (2004) and Tian et al. (2005) both indicate that there is a significant population (about one half in the case of Holder et al. (2004)) of solar active regions whose twist/tilt properties cannot be explained by the Σ-effect together with the effect of the Coriolis force alone. These active regions are consistent with the situation where the buoyant flux tube form at the base of the solar convection zone has acquired an initial twist, such that as it rises upward due to buoyancy into an Ω-tube, it develops a writhe that is of the same sense as the initial twist of the tube. In the extreme case, the twist can be so large that the flux tube becomes kink unstable (see section 5.4). The resulting tilt at the apex due to the writhe has the negative correlation with the twist as described in the above observations.

5.3 On the minimum twist needed for maintaining cohesion of rising flux tubes in the solar convection zone

As described in Section 5.1, simulations based on the thin flux tube approximation have revealed many interesting results with regard to the global dynamics of active region emerging flux loops in the solar convective envelope, which provide explanations for several basic observed properties of solar active regions. However, one major question ignored by the thin flux tube model is how a flux tube remains a discrete and cohesive object as it moves in the solar convection zone. The manner in which solar active regions emerge on the photosphere suggests that they are coherent flux bundles rising through the solar convection zone and reaching the photosphere in a reasonably cohesive fashion. To address this question, 3D MHD models that fully resolve the rising flux tubes are needed.

As a natural first step, 2D MHD simulations have been carried out to model buoyantly rising, infinitely long horizontal magnetic flux tubes in a stratified layer representing the solar convection zone, focusing on the dynamic evolution of the tube cross-section. The first of such calculations was done in fact much earlier by Schüssler (1979) and later, simulations of higher numerical resolutions have been performed (Moreno-Insertis and Emonet, 1996; Longcope et al., 1996; Fan et al., 1998a; Emonet and Moreno-Insertis, 1998). The basic result from these 2D models of buoyant horizontal flux tubes is that due to the vorticity generation by the buoyancy gradient across the flux tube cross-section, if the tube is untwisted, it quickly splits into a pair of vortex tubes of opposite circulations, which move apart horizontally and cease to rise. If on the other hand, the flux tube is sufficiently twisted such that the magnetic tension of the azimuthal field can effectively suppress the vorticity generation by the buoyancy force, then most of the flux in the initial tube is found to rise in the form of a rigid body whose rise velocity follows the prediction by the thin flux tube approximation. The result described above is illustrated in Figure 16 which shows a comparison of the evolution of the tube cross-section between the case where the buoyant horizontal tube is untwisted (upper panels) and a case where the twist of the tube is just above the minimum value needed for the tube to rise cohesively (lower panels). This minimum twist needed for tube cohesion can be estimated by considering a balance between the magnetic tension force from the azimuthal field and the magnetic buoyancy force. For a flux tube near thermal equilibrium whose buoyancy |Δρ/ρ| ∼ 1/Β, where ΔΡ = ΡΡ e denotes the density difference between the inside and the outside of the tube and β ≡ p/(B 2/8π) denotes the ratio of the gas pressure over the magnetic pressure, such an estimate (Moreno-Insertis and Emonet, 1996) yields the condition that the pitch angle ψ of the tube field lines on average needs to reach a value of order

$$\tan \Psi \equiv {{{B_\phi }} \over {{B_z}}} \mathbin{\lower.3ex\hbox{$\buildrel>\over{\smash{\scriptstyle\sim}\vphantom{_x}}$}} {\left({{a \over {{H_p}}}} \right)^{1/2}}.$$
(26)

In (26)B z and denote the axial and azimuthal field of the horizontal tube respectively, a is the characteristic radius of the tube, and H p is the local pressure scale height. The above result on the minimum twist can also be expressed in terms of the rate of field line rotation about the axis per unit length along the tube q. For a uniformly twisted flux tube, B φ = qr B z, where r is the radial distance to the tube axis. Then q needs to reach a value of order (Longcope et al., 1999)

$$q \mathbin{\lower.3ex\hbox{$\buildrel>\over{\smash{\scriptstyle\sim}\vphantom{_x}}$}} {\left({{1 \over {{H_p}a}}} \right)^{1/2}}$$
(27)

for the flux tube to maintain cohesion during its rise. Note that the conditions given by Equations (26) and (27) and also the 2D simulations described in this section all assume buoyant flux tubes with initial buoyancy |Δρ/ρ| ∼ 1/β. For tubes with lower level of buoyancy, the necessary twist is smaller with tan ψ and q both ∝ |Δρ/ρ|1/2 (see Emonet and Moreno-Insertis, 1998).

Figure 16:
figure 16

From Fan et al. ( 1998a). Upper panel: Evolution of a buoyant horizontal flux tube with purely longitudinal magnetic field. Lower panel: Buoyant rise of a twisted horizontal flux tube with twist that is just above the minimum value given by Equation (26) . The color indicates the longitudinal field strength and the arrows describe the velocity field. (For a corresponding movie showing the evolution of the tube for the untwisted and the twisted cases refer to Figure 17 .)

Figure 17:
figure 17

mpg-Movie (200 KB) Still from a movie showing The evolution of a flux tube. From Fan et al. ( 1998a) For a detailed description see Figure 16 (For video see appendix).

Longcope et al. (1999) pointed out that the amount of twist given by Equation (27) is about an order of magnitude too big compared to the twist deduced from vector magnetic field observations of solar active regions on the photosphere. They assumed that the averaged αJ z/B z (the ratio of the vertical electric current over the vertical magnetic field) measured in an active region on the photosphere directly reflects the twist in the subsurface emerging tube, i.e. q = α/2 (Longcope and Klapper, 1997; Longcope et al., 1998). If this is true then it seems that the measured twists in solar active regions directly contradict the condition for the cohesive rise of a horizontal flux tube with buoyancy as large as |Aρ/ρ| ∼ 1/β.

More recently, 3D simulations of Ω-shaped arched flux tubes have been carried out (Abbett et al., 2000, 2001; Fan, 2001a). Fan (2001a) performed 3D simulations of arched flux tubes which form from an initially neutrally buoyant horizontal magnetic layer as a result of its undulatory buoyancy instability (see Section 4.2 and Figure 7). It is found that without any initial twist the flux tubes that form rise through a distance of about one density scale height included in the simulation domain without breaking up. This significantly improved cohesion of the 3D arched flux tubes compared to the previous 2D models of buoyant horizontal tubes is not only due to the additional tension force made available by the 3D nature of the arched flux tubes, but also due largely to the absence of an initial buoyancy and a slower initial rise (Fan, 2001a). With a neutrally buoyant initial state, both the buoyancy force and the magnetic tension force grow self-consistently from zero as the flux tube arches. The vorticity source term produced by the growing magnetic tension as a result of bending and braiding the field lines is found to be able to effectively counteract the vorticity generation by the growing buoyancy force in the apex cross-section, preventing it from breaking up into two vortex rolls. The 2D models (Moreno-Insertis and Emonet, 1996; Longcope et al., 1996; Fan et al., 1998a; Emonet and Moreno-Insertis, 1998) on the other hand considered an initially buoyant flux tube for which there is an impulsive initial generation of vorticity by the buoyancy force. A significant initial twist is thus required to suppress this initial vorticity generation. Therefore the absence of an initial vorticity generation by buoyancy, and the subsequent magnetic tension force resulting from bending and braiding the field lines allow the arched tube with no net twist in Fan (2001a) to rise over a significantly greater distance without disruption.

Abbett et al. (2000) performed 3D simulations where an initial horizontal flux tube is prescribed with a non-uniform buoyancy distribution along the tube such that it rises into an Ω-shaped loop. As discussed above, due to the prescribed buoyancy in the initial horizontal tube, there is an impulsive initial generation of vorticity by the buoyancy force which breaks up the apex of the rising Ω-loop if there is no initial twist. However the separation of the two vortex fragments at the apex is reduced due to the three-dimensional effect (Abbett et al., 2000). By further including the effect of solar rotation using a local f-plane approximation, Abbett et al. (2001) found that the influence of the Coriolis force significantly suppresses the degree of fragmentation at the apex of the Ω-loop (Figure 18). They also found that the Coriolis force causes the emerging loop to become asymmetric about the apex, with the leading side (leading in the direction of rotation) having a shallower angle with respect to the horizontal direction compared to the following (Figure 18), consistent with the geometric asymmetry found in the thin flux tube calculations (Section 5.1.4).

Figure 18:
figure 18

The rise of a buoyant Ω-loop with an initial field strength B = 105 G in a rotating model solar convection zone at a local latitude of 15° (from Abbett et al. ( 2001)). The Ω-loop rises cohesively even though it is untwisted. The loop develops an asymmetric shape with the leading side (leading in the direction of rotation) having a shallower angle relative to the horizontal direction compared to the following side.

Another interesting possibility is suggested by the 3D simulations of Dorch and Nordlund (1998), who showed that a random or chaotic twist with an amplitude similar to that given by Equation (26) or (27) in the flux tube can ensure that the tube rises cohesively. Such a random twist may not be detected in the photosphere measurement of active region twists which is determined by taking some forms of average of the quantity α = J z/B z over the active region.

5.4 The rise of kink unstable magnetic flux tubes and the origin of delta-sunspots

Most of the solar active regions are observed to have very small twists, with an averaged value of α ≡ 〈J z/B z〉 (the averaged ratio of the vertical electric current over the vertical magnetic field) measured to be on the order of 0.01 Mm−1 (see Pevtsov et al., 1995, 2001). However there is a small but important subset of active regions, called the δ-sunspots, which are observed to be highly twisted with a reaching a few times 0.1 Mm−1 (Leka et al., 1996), and to have unusual polarity orientations that are sometimes reversed from Hale’s polarity rule (see Zirin and Tanaka, 1973; Zirin, 1988; Tanaka, 1991). These δ-sunspots are compact structures where umbrae of opposite polarity are contained within a common penumbra. They are found to be the most flare productive active regions (see Zirin, 1988). Through careful analysis of the evolution of flare-active δ-sunspot groups, Tanaka (1991) proposed a model of an emerging twisted flux rope with kinked or knotted geometry to explain the observed evolution of these regions.

Motivated by the observations of δ-sunspots, MHD calculations of the evolution of highly twisted, kink unstable magnetic flux tubes in the solar convection zone have been carried out (Linton et al., 1996, 1998, 1999; Fan et al., 1998b, 1999). For an infinitely long twisted cylindrical flux tube with axial field B z (r), azimuthal field B θ (r) = q (r) rB z (r), and plasma pressure p (r) in hydrostatic equilibrium where \(dp/dr=-(B_{\theta}^{2}/4\pi r)-d(B_{z}^{2}+B_{\theta}^{2})/dr\), a sufficient condition for the flux tube to be kink unstable is (see Freidberg, 1987)

$${r \over 4}{\left({{{q^{\prime}} \over q}} \right)^2} + {{8\pi p^{\prime}} \over {B_z^2}} < 0$$
(28)

to be true somewhere in the flux tube. In eq. (28) the superscript ′ denotes the derivative with respect to r. This is known as Suydam’s criterion. Note that condition (28) is sufficient but not necessary for the onset of the kink instability and hence there can be cases which are kink unstable but do not satisfy condition (28). One such example are the force-free twisted flux tubes which are shown to be always kink unstable without line-tying (i.e. infinitely long) (Anzer, 1968), but for which p′ = 0. Force-free fields are the preferred state for coronal magnetic fields under low plasma-β conditions and are not a likely state for magnetic fields in the high-β plasma of the solar interior. Linton et al. (1996) considered the linear kink instability of uniformly twisted cylindrical flux tubes with q = B θ/rB z being constant, confined in a high β plasma. They found that the equilibrium is kink unstable if q exceeds a critical value q cr = a −1, where a −2 is the coefficient for the r 2 term in the Taylor series expansion of the equilibrium axial magnetic field B z about the tube axis: B z(r) = B 0(1 − a −2r 2 + ⋯). This result is consistent with Suydam’s criterion. They further argued that an emerging, twisted magnetic flux loop will tend to have a nearly uniform q along its length since the rise speed through most of the solar convection zone is sub-Alfvénic and torsional forces propagating at the Alfvén speed will equilibrate quasi-statically. Meanwhile expansion of the tube radius at the apex as it rises will result in a decrease in the critical twist q cr = a −1 necessary for the instability. This implies that as a twisted flux tube rises through the solar interior, a tube that is initially stable to kinking may become unstable as it rises, and that the apex of the flux loop will become kink unstable first because of the expanded tube cross-section there (Parker, 1979; Linton et al., 1996).

The non-linear evolution of the kink instability of twisted magnetic flux tubes in a high-β plasma has been investigated by 3D MHD simulations (Linton et al., 1998, 1999; Fan et al., 1998b, 1999). Fan et al. (1998b, 1999) modeled the rise of a kink unstable flux tube through an adiabatically stratified model solar convection zone. In the case where the initial twist of the tube is significantly supercritical such that the e-folding growth time of the most unstable kink mode is smaller than the rise time scale, they found sharp bending of the flux tube as a result of the non-linear evolution of the kink instability. During the onset of the kink instability, the magnetic energy decreases while the magnetic helicity is approximately conserved. The writhing of the flux tube also significantly increases the axial field strength and hence enhances the buoyancy of the flux tube. The flux tube rises and arches upward at the portion where the kink concentrates, with a rotation of the tube orientation at the apex that exceeds 90° (see Figure 19). The emergence of this kinked flux tube can give rise to a compact magnetic bipole with polarity order inverted from the Hale polarity rule (Figure 20) as often seen in δ-sunspots. The conservation of magnetic helicity requires that the writhing of the tube due to the kink instability is of the same sense as the twist of the field lines. Hence for a kinked emerging tube, the rotation or tilt of the magnetic bipole from the east-west polarity orientation defined by the Hale’s polarity rule should be related to the twist of the tube. The rotation or tilt should be clockwise (counterclockwise) for right-hand-twisted (left-hand-twisted) flux tubes. This tilt-twist relation can be used as a means to test the model of kinked flux tubes as the origin of δ-sunspots (Tanaka, 1991; Leka et al., 1994, 1996; Lópezs Fuentes et al., 2003). Observations have found with both consistent and opposing cases (Leka et al., 1996; López Fuentes et al., 2003). A recent study (Tian et al., 2005) which includes a large sample (104) of complex δ-configuration active regions shows that 65–67% of these δ-regions have the same sign of twist and writhe, supporting the model of kinked flux tubes.

Figure 19:
figure 19

mpg-Movie (321 KB) (Still from a movie showing The rise of a kink unstable magnetic flux tube through an adiabatically stratified model solar convection zone (result from a simulation in Fan et al. ( 1999) with an initial right-handed twist that is 4 times the critical level for the onset of the kink instability). In this case, the initial twist of the tube is significantly supercritical so that the e-folding growth time of the most unstable kink mode is smaller than the rise time scale. The flux tube is perturbed with multiple unstable modes. The flux tube becomes kinked and arches upward at the center where the kink concentrates, with a rotation of the tube orientation at the apex that exceeds 90° (For video see appendix).

Figure 20:
figure 20

A horizontal cross-section near the top of the upward arching kinked loop shown in the last panel of Figure 19 . The contours denote the vertical magnetic field B zwith solid (dotted) contours representing positive (negative) B z. The arrows show the horizontal magnetic field. One finds a compact bipolar region with sheared transverse field at the polarity inversion line. The apparent polarity orientation (i.e. the direction of the line drawn from the peak of the positive pole to the peak of the negative pole) is rotated clockwise by about 145° from the +xdirection (the east-west direction) of the initial horizontal flux tube.

5.5 Buoyant flux tubes in a 3D stratified convective velocity field

To understand how active region flux tubes emerge through the solar convection zone, it is certainly important to understand how 3D convective flows in the solar convection zone affect the rise and the structure of buoyant flux tubes. The well-defined order of solar active regions as described by the Hale polarity rule suggests that the emerging flux tubes are not subject to strong deformation by the turbulent convection. One can thus consider the following two simplified order-of-magnitude estimates. First, the magnetic buoyancy of the flux tube should probably dominate the downward hydrodynamic force from the convective downflows:

$${{{B^2}} \over {8\pi {H_p}}} > {C_{\rm{D}}}{{\rho v_{\rm{c}}^2} \over {\pi a}} \Rightarrow B > \left({{{2{C_{\rm{D}}}} \over \pi }} \right){\left({{{{H_p}} \over a}} \right)^{1/2}}{B_{{\rm{eq}}}},$$
(29)

where B eq ≡ (4πρ)1/2v c is the field strength at which the magnetic energy density is in equipartition with the kinetic energy density of the convective downdrafts, v c is the flow speed of the downdrafts, H p is the local pressure scale height, a is the tube radius, and C D is the aerodynamic drag coefficient which is of order unity. In (29) we have used the aerodynamic drag force as an estimate for the magnitude of the hydrodynamic forces. The estimate (29) leads to the condition that the field strength of the flux tube needs to be significantly higher than the equipartition field strength by a factor of \(\sqrt{H_{p}/a}\). For flux tubes responsible for active region formation, H p/a ≥ 10 near the bottom of the solar convection zone. A second consideration is that the magnetic tension force resulting from bending the flux tube by the convective flows should also dominate the hydrodynamic force due to convection:

$${{{B^2}} \over {8\pi {l_{\rm{c}}}}} > {C_{\rm{D}}}{{\rho v_{\rm{c}}^2} \over {\pi a}} \Rightarrow B > \left({{{{C_{\rm{D}}}} \over \pi }} \right){\left({{{{l_{\rm{c}}}} \over a}} \right)^{1/2}}{B_{{\rm{eq}}}},$$
(30)

where l c denotes the size scales of the convective flows. This leads to a condition for the magnetic field strength very similar to Equation (29) if we consider the largest convective cell scale for l c to be comparable to the pressure scale height.

Fan et al. (2003) carried out direct 3D MHD simulations of the evolution of a buoyant magnetic flux tube in a stratified convective velocity field. The basic result is illustrated in Figure 21. They first computed a 3D convective velocity field in a superadiabatically stratified fluid, until the convection reaches a statistical steady state. The resulting velocity field (see top-left image in Figure 21) shows the typical features of overturning convection in a stratified fluid as found in many previous investigations. The surface layer displays a cellular pattern with patches of upflow region surrounded by narrow downflow lanes. In the bulk of the convecting domain, the downflows are concentrated into narrow filamentary plumes, some of which extend all the way across the domain, while the upflows are significantly broader and are of smaller velocity amplitude in comparison to the downdrafts. A uniformly buoyant, twisted horizontal magnetic flux tube having an entropy that is equal to the entropy at the base of the domain is inserted into the convecting box (see top-right image in Figure 21). In the case where the field strength of the tube is in equipartition to the kinetic energy density of the strongest downdraft (left column in the bottom panel of Figure 21), the magnetic buoyancy for this flux tube is weaker than the hydrodynamic force resulting from the convective downflows and the evolution of the tube depends sensitively on the local condition of the convective flows. Despite being buoyant, the portions of the tube in the paths of downdrafts are pushed downward and pinned down to the bottom, while the rise speed of sections within upflow regions is significantly boosted. The pinned-down flux is then further distorted and transported laterally by the horizontal diverging flow at the bottom. On the other hand in the case where the tube field strength is 10 times the equipartition value (right column in the bottom panel of Figure 21), the horizontal flux tube rises under its uniform buoyancy, nearly unaffected by the convection. In this case the horizontal flux tube is sufficiently twisted so that it does not break up into two vortex tubes.

Figure 21:
figure 21

The evolution of a uniformly buoyant magnetic flux tube in a stratified convective velocity field from the simulations of Fan et al. ( 2003). Top-left image: A snapshot of the vertical velocity of the 3D convective velocity field in a superadiabatically stratified fluid. The density ratio between the bottom and the top of the domain is 20. Top-right image: The velocity field (arrows) and the tube axial field strength (color image) in the vertical plane that contains the axis of the uniformly buoyant horizontal flux tube inserted into the convecting box. Lower panel: The evolution of the buoyant flux tube with B = B eq (left column) and with B = 10Beq (right column). The color indicates the absolute field strength of the flux tube scaled to the initial tube field strength at the axis. (For a corresponding movie see Figure 22 .)

Figure 22:
figure 22

mpg-Movie (336 KB) Still from a movie showing The evolution of a uniformly buoyant magnetic flux tube. From Fan et al. ( 2003). For a detailed description see Figure 21 (For video see appendix).

In case B = B eq, it is found that the random north-south tilting of the flux tube caused by convection is of the amplitude ∼ 30°, which is greater than the r.m.s. scatter of the active region tilts away from Joy’s law for large active regions (10°), but is not beyond the r.m.s. tilt scatter for small active regions (30°) (see Fisher et al., 1995). This indicates that the distortion of flux tubes of equipartition field strength by the convective flows during their buoyant rise through the solar convection zone is probably too large to be consistent with the observational constraint of tilt dispersion for large solar active regions. Furthermore it should be noted that the realistic convective flows in the solar convection zone is probably far more turbulent than that computed in Fan et al. (2003), containing flows of scales significantly smaller than the cross-section of the flux tube. Hence the flux tube distortion found in the simulations of Fan et al. (2003) is most likely a lower limit.

The distances between the major downflow plumes are also an important factor in determining the fate of the buoyant flux tubes of equipartition field strength. In Fan et al. (2003), the distance between neighboring downflow plumes can be as large as about 5H p, hence allowing the portion of the buoyant tube between the plumes to rise up to the top of the domain. If on the other hand, the distances between the downflow plumes are below ∼ 2H p, then the tension force for the tube between the pinned-down points will exceed the magnetic buoyancy force (∼ B 2/8πH p) and the entire flux tube will be prevented from emerging to the surface.

6 Turbulent Pumping of a Magnetic Field in the Solar Convection Zone

As discussed in Section 5.5, for buoyant flux tubes with significantly super-equipartition field strength B ≳ (H p/a)1/2B eq, where B eq is in equipartition with the kinetic energy density of the convective downflows, the magnetic buoyancy of the tubes dominates the hydrodynamic force from the convective downflows and the flux tubes can rise unaffected by convection. On the other hand if the field strength of the flux tubes is comparable to or smaller than the equipartition value B eq, the magnetic buoyancy is weaker than the hydrodynamic force from the convective downflows and the evolution of the tubes becomes largely controlled by the convective flows. In this regime of convection dominated evolution, due to the strong asymmetry between up- and downflows characteristic of stratified convection, it is found that magnetic flux is preferentially transported downward against its magnetic buoyancy out of the turbulent convection zone into the stably stratified overshoot region below. This process of “turbulent pumping” of a magnetic field has been demonstrated by several high resolution 3D compressible MHD simulations (see Tobias et al., 1998, 2001; Dorch and Nordlund, 2001).

Tobias et al. (2001) carried out a series of 3D MHD simulations to investigate the turbulent pumping of a magnetic field by stratified convection penetrating into a stably stratified layer. A thin slab of a unidirectional horizontal magnetic field is introduced into the middle of an unstably stratified convecting layer, which has a stable overshoot layer attached below. It is found that the fast, isolated downflow plumes efficiently pump flux from the convecting layer into the stable layer on a convective time scale. Tobias et al. (2001) quantify this flux transport by tracking the amount of flux in the unstable layer and that in the stable layer, normalized by the total flux. Within a convective turnover time, the flux in the stable layer is found to increase from the initial value of 0 to a steady value that is greater than 50% of the total flux (reaching 80% in many cases), i.e. more than half of the total flux is settled into the stable overshoot region in the final steady state. Moreover the stable overshoot layer is shown to be an effective site for the storage of a toroidal magnetic field. If the initial horizontal magnetic slab is put in the stable overshoot layer, the penetrative convection is found to be effective in pinning down the majority of the flux against its magnetic buoyancy, preventing it from escaping into the convecting layer. It should be noted however that the fully compressible simulations by Tobias et al. (2001) assumed a polytropic index of m = 1 for the unstable layer which corresponds to a value for the non-dimensional superadiabaticity δ of 0.1. Thus the convective flow speed v c for the strong downdrafts is on the order of \(\sqrt{\delta}\sim 0.3\) times the sound speed, i.e. not very subsonic. On the other hand, the range of initial magnetic field strength considered in the simulations corresponds to a plasma β ranging from 2 × 103 to 2 × 107. Thus comparing to the kinetic energy density of the convective flows, the strongest initial field considered is of the order 0.1 B eq, i.e. below equipartition, although a field with strength up to B eq may be generated as a result of amplification by the strong downflows during the evolution. Therefore the above results obtained with regard to the efficient turbulent pumping of a magnetic field out of the convection zone into the stable overshoot region against its magnetic buoyancy apply only to fields weaker than or at most comparable to B eq.

The turbulent pumping of magnetic flux with field strength BB eq out of the convection zone into the stable overshoot region demonstrated by the high resolution 3D MHD simulations has profound implications for the working of the interface mean field dynamo models (Parker, 1993; Charbonneau and MacGregor, 1997) as discussed by Tobias et al. (2001). The interface dynamo models require efficient transport of the large scale poloidal field generated by the α-effect of the cyclonic convection out of the convection zone into the stably stratified tachocline region where strong rotational shear generates and amplifies the large scale toroidal magnetic field. Turbulent pumping is shown to enhance both the transport of magnetic flux into the stable shear layer and the storage of the toroidal magnetic field there. It further implies that the transport of magnetic flux by turbulent pumping should not be simply treated as an enhanced isotropic turbulent diffusivity in the convection zone as is typically assumed in the mean field dynamo models. Tobias et al. (2001) noted that a better treatment would be to add an extra advective term to the mean field equation characterizing the effect of turbulent pumping, which would correspond to including the antisymmetric part of the α-tensor known as the γ-effect (see Moffatt, 1978).

7 Amplification of a Toroidal Magnetic Field by Conversion of Potential Energy

Thin flux tube models of emerging flux loops through the solar convective envelope (Section 5.1) have inferred a strong super-equipartition field strength of order 105 G for the toroidal magnetic field at the base of the solar convection zone. Generation of such a strong field is dynamically difficult since the magnetic energy density of a 105 G field is about 10–100 times the kinetic energy density of the differential rotation (Parker, 1994; Rempel and Schüssler, 2001). An alternative mechanism for amplifying the toroidal magnetic field has been proposed which converts the potential energy associated with the stratification of the convection zone into magnetic energy. It is shown that upflow of high entropy plasma towards the inflated and “exploding” top of a rising Ω-loop developed from an initial toroidal field of equipartition field strength (∼ 104 G) can significantly intensify the submerged part of the field by extracting plasma out of it (Parker, 1994; Moreno-Insertis et al., 1995; Rempel and Schüssler, 2001). This process is a barometric effect and is caused by the entropy gradient in the solar convection zone maintained by the energy transport.

In the thin flux tube simulations of rising Ω-loops in the solar convection zone, it is found that flux loops with a low initial field strength of ∼ 104 G do not reach the upper half of the solar convection zone before the apexes of the loops loose pressure confinement and effectively “explode” (Moreno-Insertis et al., 1995). This loss of pressure confinement at the top of the emerging loop occurs as plasma inside the tube establishes hydrostatic equilibrium along the tube which happens if the emerging loop rises sufficiently slowly (Moreno-Insertis et al., 1995). The loop rises adiabatically carrying the high entropy plasma from the base of the solar convection zone while the entropy outside the flux tube decreases with height in the superadiabatically stratified convection zone. Hydrostatic equilibrium therefore dictates that the plasma pressure inside the flux tube decreases with height slower than the outside and becomes equal to the external pressure at a certain height where the magnetic field can no longer be confined. This “explosion” height for the emerging loop is found to be a function of the initial tube field strength at the base of the solar convection zone (see Figure 23). For loops with an initial field strength ∼ 104 G, the explosion height is at about the middle of the solar convection zone. When the loop apex approaches the explosion height, it expands drastically and the buoyancy of the high entropy material in the tube is expected to drive an outflow which extracts plasma out of the lower part of the flux tube at the base of the solar convection zone. This process has been demonstrated by Rempel and Schüssler (2001), who performed MHD simulations of exploding magnetic flux sheets in two-dimensional Cartesian geometry.

Figure 23:
figure 23

Depth from the surface where the apex of an emerging flux loop with varying initial field strength rising from the bottom of the convection zone looses pressure confinement or “explodes” as a result of tube plasma establishing hydrostatic equilibrium along the tube. The explosion height is computed by considering an isentropic thin flux loop with hydrostatic equilibrium along the field lines (see Moreno-Insertis et al., 1995) in a model solar convection zone of Christensen-Dalsgaard (Christensen-Dalsgaard et al., 1993).

The simulations of Rempel and Schüssler (2001) start with a magnetic sheet with a higher value of entropy placed at the bottom of an adiabatically stratified layer (constant entropy layer). This setup avoids the complication of involving convective flows in the simulations while keeping the essential effect of the entropy decrease in the solar convection zone by assuming a constant entropy difference between the flux sheet and the isentropic layer. The central portion of the flux sheet is perturbed upward which subsequently forms a rising loop as a result of its magnetic buoyancy (see Figure 24). The apex of the rising loop explodes into a cloud of weak magnetic field as it crosses the predicted explosion height (middle panel) and high entropy plasma, driven by buoyancy, continues to flow out of the “stumps”, draining mass from the lower horizontal part of the flux sheet (bottom panel). The field strength of the horizontal part of the flux sheet is visibly intensified. It is found that the final field strength the horizontal part of the field can reach is roughly the value for which the explosion height is close to the top of the stratification. For the solar convection zone, this value corresponds to ∼ 105 G (see Figure 23). This implies that the large field strength of order 105 G may be achieved by the process of flux “explosion”, which draws energy from the potential energy associated with the stratification of the solar convection zone. A coherent picture maybe as follows. Differential rotation in the tachocline shear layer generates and amplifies the toroidal magnetic field to an equipartition value of about 104 G, at which point magnetic buoyancy becomes important dynamically, and magnetic buoyancy instability, radiative heating or other convective perturbations drive the formation of buoyant flux loops rising into the solar convection zone. These rising loops explode in the middle of the convection zone and fail to rise all the way to the surface. These “failed” eruptions “pump” out the plasma from the toroidal field at the base of the solar convection zone (Parker, 1994) and amplify the toroidal field until it reaches strongly super-equipartition field strength of order 105 G, whose eruptions then lead to the emergence of solar active regions at the surface.

Figure 24:
figure 24

From Rempel and Schüssler ( 2001). Evolution of magnetic field strength (gray scale: darker gray denotes stronger field) and velocity field (arrows) during the flux loop explosion. The horizontal part of the field is amplified by a factor of 3.

8 Flux Emergence at the Surface and Post-Emergence Evolution of Subsurface Fields

8.1 Evolution in the top layer of the solar convection zone

A complete 3D MHD model of emerging active region flux tubes, extending from the base of the solar convection zone up into to the visible solar atmosphere, is not yet possible. Thin flux tube calculations of rising flux tubes in the solar convection zone usually extend from the base of the solar convection zone up to roughly 20–30 Mm below the surface, at which depths the validity of the thin flux tube approximation breaks down since the diameter of the flux tube begins to exceed the local pressure scale height (Moreno-Insertis, 1992). The domain of fully 3D MHD simulations of rising magnetic flux tubes in a model convection zone typically contains a vertical stratification of up to 3 density scale heights (see Abbett et al., 2000, 2001; Fan et al., 2003), which corresponds to the density stratification over the range from the bottom of the solar convection zone to roughly 36 Mm below the photosphere. Because of the rapid decrease of the various scale heights in the top layer of the solar convection zone which demands increasing numerical resolution, it is not yet feasible to perform 3D MHD simulations that extend from the bottom of the convection zone all the way to the photosphere. Furthermore, there is an increased complexity in the physics of the top layer of the solar convection zone. The thermodynamics of the plasma is complicated by ionization effects and the radiative exchange is expected to play an important role in the heat transport (see Lantz, 1991; Stein and Nordlund, 2000). The anelastic approximation breaks down because the plasma flow speed is no longer slow compared to the sound speed. Fully compressible 3D MHD simulations of magneto-convection and emerging magnetic flux in the top few Mm layer of the solar convection zone and the overlying photospheric layer, incorporating realistic physics such as partial ionization of the dominant constituents and radiative transfer, have been carried out (see review by Stein and Nordlund, 2000, and the references therein). These simulations have produced results that can be directly compared with high resolution photospheric observations of the solar granulations and have found the formation of intense magnetic flux tubes of kG field strength with size scales up to that of a magnetic pore. However realistic simulations of the evolution of emerging active region scale flux tubes in the top ∼ 20 Mm layer of the solar convection zone are still beyond the reach of current numerical modeling capabilities. It remains an open question how the top of the emerging Ω-tube intensifies to form sunspots with kG field strength and β ∼ 1 at the photosphere. The intensification of small scale magnetic flux tubes to kG field strength at the photosphere can be explained by the process of “convective collapse” which results from the convective instability of plasma inside the vertical thin flux tubes in the top few hundred kilometers of the solar convection zone (see Parker, 1978; Spruit, 1979; Spruit and Zweibel, 1979). This process may also be responsible for the intensification of sunspot scale flux tubes, but there have been no realistic model calculations.

8.2 Flux emergence into the solar atmosphere

A fundamental and challenging question concerning the emergence of solar active regions is how buoyant magnetic fields from the convectively unstable interior emerge dynamically into the stably stratified, rarefied solar atmosphere and corona. Pioneering work in this area was carried out by Shibata and collaborators (see Shibata et al., 1989) who showed that the magnetic buoyancy instability is a mechanism through which magnetic flux reaching the photosphere can expand dynamically into the stably stratified solar atmosphere. Many 2D and 3D simulations of flux emergence into the atmosphere have since been carried out (see e.g. Magara, 2001; Murray et al., 2006, and the references therein).

Recently, a series of 3D MHD simulations (Fan, 2001b; Magara and Longcope, 2001, 2003; Magara, 2004; Manchester IV et al., 2004; Archontis et al., 2004, 2005, 2006; Galsgaard et al., 2005; Murray et al., 2006; Magara et al., 2005; Magara, 2006) have been carried out to model the emergence of a twisted magnetic flux tube through a multi-layered atmosphere that includes a polytropic layer representing the top of the solar convection zone, an isothermal layer representing the photosphere and chromosphere, and another isothermal layer representing the million-degree corona. In these simulations, a twisted flux tube is initially placed in the convection zone at a depth of about a few times the photospheric pressure scale height and the central segment of the twisted tube rises buoyantly towards the photosphere. When the top of the tube enters the photosphere the rise velocity at the upper boundary of the tube first slows down as it encounter the stable stratification of the photosphere. As a result, magnetic flux begins to pile up at the photosphere and a steep magnetic pressure gradient is established. It is found that subsequently, the magnetic flux undergoes a run-away expansion into the atmosphere due to the non-linear growth of the magnetic buoyancy instability. Archontis et al. (2004) and Murray et al. (2006) quantitatively examined their simulation data using the criterion for the onset of the undulatory magnetic buoyancy instability (Newcomb, 1961; Thomas and Nye, 1975; Acheson, 1979):

$$ - {H_p}{\partial \over {\partial z}}(\log B) > - {\gamma \over 2}\beta \delta + k_\parallel ^2\left({1 + {{k_ \bot ^2} \over {k_z^2}}} \right),$$
(31)

where z is height, H p is the local pressure scale height at the photosphere, B is the magnetic field strength, γ is the ratio of the specific heats, β is the ratio of the plasma pressure over the magnetic pressure, δ ≡ ∇ − ∇ad is the superadiabaticity of the atmosphere, which is −0.4 for the isothermal stratification of the photosphere, and k , k and k z are the three components of the local perturbation wave vector (normalized by 1/H p), with k and being the two horizontal components parallel and perpendicular to the local magnetic field direction, and k z being the z component. They found that the run-away expansion takes place at the time when the above critical condition is met. Furthermore, Murray et al. (2006) found that if either the field strength or the twist of the subsurface flux tube is too low, the magnetic pressure build-up at the photosphere may not achieve the critical condition given above and flux emergence into the atmosphere may fail to take place.

During the run-away expansion of the magnetic flux tube into the solar atmosphere, strong diverging downflows are found along the emerged field lines, forming shock fronts just above the photosphere (e.g. Fan, 2001b; Magara and Longcope, 2003; Archontis et al., 2004). Due to the magnetic tension force associated with the twisted field lines, a shear flow pattern develops on the photosphere with the plasma on the two sides of the polarity inversion line moving oppositely in the east-west direction (Manchester IV, 2001; Fan, 2001b; Magara and Longcope, 2003; Manchester IV et al., 2004; Archontis et al., 2004). The effect of the tension force driven shear flow is to transport the axial magnetic flux upward into the expanding portion of the flux tube in the solar atmosphere (Manchester IV et al., 2004). Newly developing active regions on the photosphere often exhibit such a shear flow pattern (Zirin, 1988; Strous et al., 1996). The two polarities of the emerging region on the photosphere initially emerge as north-south oriented, and over time they appear to shear along the polarity inversion line and separate in the east-west direction (e.g. Fan, 2001b; Magara and Longcope, 2003; Manchester IV et al., 2004; Archontis et al., 2004; Magara, 2006). A detailed study of both the apparent and plasma motions in the photosphere resulting from the emergence of a twisted flux tube (with different degree of twist for the initial subsurface flux tube) is carried out by Magara (2006). It is found that the apparent rotational motion exhibited by the two polarities are opposite to the direction of the actual rotational motion of the plasma.

The simulations also suggest that it is difficult for a twisted 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. Independent simulations by Fan (2001b) and Archontis et al. (2004), considering nearly identical initial conditions both find that that as the upper parts of the winding field lines of the twisted flux tube expand into the atmosphere and the corona, the axial field line of the tube settles to an equilibrium not much above the photosphere, and the lower U-shaped parts of the field lines largely remain trapped below the photosphere. On the other hand, Manchester IV et al. (2004) considered the case where the length of the emerging portion of the tube is shortened to a half of that in Fan (2001b), such that the field lines only perform about one turn about the axis over the emerging segment. In this case it is found that after the growth of the magnetic buoyancy instability which causes the expansion of the flux tube into the solar atmosphere, a secondary eruptive phase sets in where the upper emerged part of the flux tube begins to “pinch off” from the lower mass-laden part via magnetic reconnection, forming a new flux rope with a new central axial field line that rises well into the corona. Manchester IV et al. (2004) explains that the cause of the eruptive phase is the shear flow driven by the magnetic tension force, which continually transport axial flux and magnetic energy to the upper expanded portion of the magnetic field, driving its eruption. It is found that a sigmoid shaped current sheet develops in the chromospheric heights during the eruptive phase. Simulations by Magara and Longcope (2001), Magara (2004) Magara et al. (2005), and Magara (2006) also found formation of current concentration of similar sigmoid morphology in the chromospheric heights (see e.g. Figure 25), and the ascent of a flux rope structure into the atmosphere with an O-point of the magnetic field in the central cross-section reaching the base of the corona. It is found that the emerged field lines with footpoints rooted in regions of high electric current density in the chromosphere display an inverse-S shape for a left-hand-twisted emerging tube, in agreement with the X-ray sigmoid morphology and the sense of active region twist in the northern hemisphere (see Figure 25 and Figure 4).

Figure 25:
figure 25

Emerging magnetic field in the solar atmosphere resulting from the 3D simulation of the emergence of a left-hand-twisted magnetic flux tube by Magara ( 2004). The colors of the field lines represent the square value of the current density at their footpoints on a chromospheric plane located at z = 5. top left: Top view of the magnetic field lines. Note the inverse-S shape of the brighter field lines, which is consistent with the X-ray sigmoid morphology preferentially seen in the northern hemisphere. top right: The square of the current density (color image) and vertical magnetic flux (contours) at the chromospheric plane. bottom left: Side view of the magnetic field lines. bottom right: Another perspective view of the magnetic field.

The effects of the presence of a simple horizontal coronal magnetic field and the interaction and reconnection of the pre-existing coronal field with the emerging flux rope have been investigated in a series of work (Archontis et al., 2004, 2005, 2006; Galsgaard et al., 2005). It is found that the dynamic interaction and the rate of magnetic reconnection depends sensitively on the relative orientations between the upcoming emerging flux and the pre-existing coronal magnetic field. When the two flux systems coming into contact have a relative angle above 90°, the simulations show immediate and substantial magnetic reconnection, producing collimated high-speed and high-temperature jets from the reconnection site. As a result of the reconnection, most of the flux in the subsurface flux rope end up connecting to the coronal magnetic field. The interaction with the pre-existing coronal field tends to slow down the rise speed of the upper-boundary of the emerging flux tube into the corona, compared to the case where the flux tube emerges into a field-free atmosphere (Archontis et al., 2005). Isobe et al. (2005) performed high resolution 3D simulations of emerging flux and its reconnection with pre-existing magnetic field in the corona. In this calculation the emerging flux develops from an initial horizontal magnetic flux sheet (with uni-direction field) situated in the top of the convection zone. Due to the growth of the three-dimensional magnetic buoyancy instability of modes with high wave number in the direction perpendicular to the field (or magnetic Rayleigh-Taylor instability), undulating flux bundles with fine spacial scales rise into the corona and reconnect with the pre-existing coronal magnetic field in a spatially intermittent way. This results in the formation of filamentary structure resembling the observed arch filament systems in emerging flux regions (EFRs). The spatially intermittent reconnection and heating also explain the coexistence of many hot and cold loops and the jets being ejected from the loop footpoints in EFRS observed in EUV by the TRACE satellite.

Due to the need to resolve the photosphere pressure scale height (∼ 150 km), the 3D calculations of twisted emerging flux tubes described above are done for a domain size of up to a few tens of Mm, the size of a small active region. The implicit assumption is that the basic dynamical behavior of the smaller emerging flux tubes modeled in these simulations is representative of that of the larger active region scale emerging tubes. Furthermore, adiabatic evolution of an ideal gas is assumed in these simulations, which is not appropriate for the atmosphere layers considered. The adiabatic expansion of the tube plasma emerging into the atmosphere resulting in unphysically low temperatures. Radiative heat exchange, thermal conduction, and the still uncertain sources of coronal heating all play an important role in the thermal energy evolution of the plasma in emerging flux regions.

8.3 Post-emergence evolution of subsurface fields

This is so far a largely unexplored area with mostly speculations and very few quantitative calculations. One outstanding question is whether and to what extent sunspots and active region fields at the surface remain connected to the subsurface flux tubes in the deep solar convection zone. If the surface fields remain connected to the subsurface flux tubes, then based on the shape of the Ω-loop as shown in Figure 12, it seems difficult to understand why the separation of the two polarities of a solar active region stops at the scales of ≲ 100 Mm, without continuously increasing to at least the footpoint separation of the Ω-loop at the base of the convection zone, which is ∼ 1000 Mm! Due to the angles between the legs of the Ω-loop and the surface, the magnetic tension of the legs is expected to be pulling the two poles of the bipolar region apart until the legs become vertical. However there is also an attracting force between the two poles due to the emerged magnetic field above the photosphere. If we assume that the field above the photosphere is potential, then the attracting force of the two poles due to the field above can be estimated by considering the attraction between two opposite electric charges of ±q, where q = Φ/2π and Φ is the total magnetic flux of the active region. This yields (Φ/2π)2/d 2 for the attracting force, where d is the separation between the poles. On the other hand, the lateral force from the subsurface field that acts to pull the two poles apart can be estimated by integrating over the area of each pole the Maxwell stress exerted by the subsurface field from below. This leads to an estimated lateral force of \((B_{n}^{2}/4\pi){\rm ctg}\theta\Delta S\), where B n and ΔS are the photosphere normal field strength and area of each magnetic pole respectively (note that we have B nΔS = Φ), and θ is the angle between the slanted subsurface tube and the surface. Balancing the forces from above and below the photosphere (taking a representative value of θ = 45°), we find that it is possible for the two poles to reach a lateral force balance if πd 2 = ΔS. This suggests that the two poles of the bipolar region need to be very close, with the polarity separation no greater than the radius of each pole, in order for the attracting force from above to balance the force that drives the separation from below. Furthermore it appears that the force balance is an unstable equilibrium because of the 1/d 2 dependence of the attracting force which weakens in response to an increase in d.

One possible resolution of the above problem is that some other processes that can form Ω-loops of much shorter length scale compared to the undulatory buoyancy instability (Section 4.1) is responsible for the formation of active region emerging tubes, although it is not immediately clear what these processes might be. Another alternative is that the active region magnetic fields on the photosphere become dynamically disconnected from the interior flux tubes. Fan et al. (1994) speculated on the process of “dynamic disconnection” which has the same physical cause as that of the so called flux tube “explosion” described in Section 7. It can be seen from Figure 23, that even for an emerging flux loop with a field strength at the convection zone base that is as high as 105 G, the flux tube is expected to loose pressure confinement and hence “explode” at a height of about 10 Mm below the surface as a result of plasma establishing hydrostatic equilibrium (HE) along the tube. While a flux loop is rising, especially in the top few tens of Mm of the convection zone where the rise speed becomes comparable to the Alfvén speed, the variation with depth of the magnetic field strength can deviate significantly from the conditions of HE. However, after flux emergence at the surface, the plasma inside the submerged flux tube will try to establish HE. Tube plasma near the surface can cool through radiation and probably undergoes “convective collapse”, forming sunspots and pores (see Spruit and Zweibel, 1979; Stein and Nordlund, 2000) at the photosphere. However in deeper depths where plasma evolution is still nearly adiabatic, a catastrophic weakening of the magnetic field or flux tube “explosion” (Section 7) may occur, which is caused by an upflow of high entropy plasma as it tries to establish HE along the tube. This weakening of the field leads to the effective dynamic disconnection of the active region fields on the photosphere from the interior flux tubes. Fan et al. (1994) suggested that “dynamic disconnection” can explain

  1. 1.

    why the active region tilt persists once an active region has formed (see more recent work by Longcope and Choudhuri, 2002) and

  2. 2.

    why the 2D surface flux transport models based on diffusion by supergranular motion (see Wang et al., 1989) can explain well the observed dispersal of active regions and the migration of active region flux across latitudes.

A first quantitative calculation of the above process of “dynamic disconnection” has been carried out by Schüssler and Rempel (2005), using a 1-D self-similar vertical flux tube model whose top end has reached the photosphere. The model computes the quasi-static evolution of the flux tube under the effects of radiative cooling, convective energy transport, and a pressure buildup by a prescribed inflow at the bottom of the tube due to the high entropy tube plasma flowing upward to establish hydrostatic equilibrium along the field. The calculation shows that after emergence, the radiative losses near the surface drives an inward propagating cooling front accompanied by a downflow, which leads to a decrease of the gas pressure and an intensification of the magnetic field in the surface layers. In the mean time, the convergence of the radiative cooling driven downflow and the buoyancy driven upflow results in an increase of the gas pressure and a weakening of the magnetic field in the tube at a depth of a few Mms. It is found that for a reasonable range of the upflow speed, the magnetic field weakens to fall below the local equipartition value (i.e. dynamic disconnection takes place) at a depth between 2 and 6 Mm, in a time scale of up to 3 days. This is consistent with the time scale over which the evolution of a newly emerging active region changes from an “active phase” of growth with increasing polarity separation to a “passive phase” which can be well represented by transport by near surface flows and supergranular diffusion.

9 Summary

Summarizing the studies reviewed above, the following basic conclusions can be drawn:

  • Toroidal magnetic fields stored in the weakly subadiabatic overshoot region are preferably in an equilibrium state of neutral buoyancy, with the magnetic curvature force balanced by the Coriolis force due to a prograde toroidal flow (Section 3.1). This is true regardless of whether the field is in the state of an extended magnetic layer or isolated flux tubes. Detecting this prograde toroidal flow through helioseismology may be a way to probe and measure the toroidal magnetic fields stored in the tachocline region.

  • Isolated toroidal flux tubes stored in the weakly subadiabatic overshoot region experience a radiative heating due to the non-zero divergence of the radiative heat flux (Section 3.2). This radiative heating causes a quasi-static upward drift of the toroidal flux tubes. In order to maintain toroidal flux tubes in the overshoot region for a time scale comparable to the solar cycle period, a rather strong subadiabaticity of δ < −10−4 is needed, which is significantly more subadiabatic than the values obtained by most of the overshoot models based on the non-local mixing length theory. This strong subadiabaticity may be achieved as a result of some level of suppression of convective motions by the toroidal magnetic flux tubes themselves. Furthermore, a recent semi-analytical model of convective overshoot (Rempel, 2004) — based on the assumption that the convective energy transport is governed by coherent downflow plumes — shows that the overshoot region can have a subadiabaticity of δ ≲ −10−4 if the downflow filling factor at the base of the convection zone is ≲ 10−4.

  • Neutrally buoyant toroidal flux tubes stored in the weakly subadiabatic overshoot region are subject to the onset of the undulatory buoyancy instability depending on the field strength and the value of the subadiabaticity (Section 4.1). The toroidal flux tubes become buoyantly unstable and develop Ω-shaped emerging flux loops when the field strength becomes sufficiently large, or when the quasi-static upward drift due to radiative heating brings the tubes out to regions of sufficiently weak subadiabaticity or into the convection zone. A neutrally buoyant equilibrium magnetic layer is also subject to the same type of undulatory buoyancy instability and 3D MHD simulations show that arched buoyant flux tubes form as a result of the non-linear growth of the instability (Section 4.2).

  • Thin flux tube simulations of emerging flux loops in the solar convective envelope (Section 5.1) show that several asymmetries between the leading and the following sides of the emerging loops develop due to the effects of the Coriolis force. These asymmetries provide explanations for the observed active region tilts described by Joy’s law and the observed asymmetric proper motions of the two polarities of newly developing active regions. Results from these thin flux tube simulations also suggest that the field strength for the toroidal magnetic fields at the base of the solar convection zone is in the range of about 3 × 104 to 105 G, significantly higher than the equipartition field strength, in order for the emerging loops to be consistent with the observed properties of solar active regions.

  • Vector magnetic field observations on the photosphere show that on average solar active regions have a small but statistically significant mean twist that is left-handed sense in the northern hemisphere and right-handed sense in the southern hemisphere. The twist may be due to the current helicity in the dynamo generated toroidal magnetic field, from which buoyant flux tubes form at the base of the convection zone, and it may be acquired during the rise of the flux tubes through the convection zone as a result of buffeting by the helical convective motions (called the Σ-effect) and also through the accretion of the mean poloidal magnetic field in the convection zone onto the rising flux tube (Section 5.2). Observational studies of the correlation between active region twist and tilt angles indicate that the Σ-effect cannot be the only source for the twist (Section 5.2).

  • Two-dimensional MHD simulations of the rise of buoyant horizontal flux tubes show that a minimum twist is needed for the tube to rise cohesively (Section 5.3). This minimum twist depends on the initial buoyancy of the tube. For a flux tube initially in thermal equilibrium with its surrounding, the minimum twist needed is found to be an order of magnitude too large compared to the mean twist measured in the majority of solar active regions. Three-dimensional simulations of arched rising tubes that develop self-consistently from zero initial buoyancy indicate that the necessary twist needed for tube cohesion may be significantly reduced (Section 5.3). On the other hand, the rise of highly twisted, kink unstable magnetic flux tubes can produce kinked or knotted emerging tubes which may explain the origin of the unusual class of flare productive active regions called δ-sunspots which are observed to be highly twisted and often show polarity order inverted from the Hale polarity rule (Section 5.4).

  • In order for the magnetic buoyancy force of an emerging flux tube to dominate the hydrodynamic force from the strong downflows in a 3D stratified convective velocity field, the field strength of the flux tube needs to be greater than (H p/a)1/2B eq, where H p is the pressure scale height, a is the tube radius, and B eq is the field in equipartition with the kinetic energy density of the strong downdrafts. For emerging active region flux tubes in the deep convection zone (H p/a)1/2 ≲ 3. For a buoyant flux tube of equipartition field strength, the tube is pinned down at the locations of strong downdrafts although sections of the tube between downdrafts may still emerge (Section 5.5).

  • Magnetic fields of BB eq are found to be preferentially transported downward against their magnetic buoyancy out of the turbulent convection zone into the stably stratified overshoot region by compressible penetrative convection (Section 6). This turbulent pumping mechanism is complementary for storage of toroidal magnetic fields in the overshoot layer. It also provides an effective mechanism for transporting the weak poloidal mean field from the convection zone into the tachocline layer which is important for the working of the interface mean field dynamo models. However the magnetic buoyancy of significantly super-equipartition toroidal fields of e.g. ∼ 105 G cannot be balanced by convective “pummeling” (Section 5.5).

  • Toroidal fields of significantly super-equipartition field strengths of 3 × 104 to 105 G inferred from simulations of rising flux tubes are difficult to generate dynamically. However, because of the entropy gradient in the solar convection zone, the formation of emerging loops with equipartition initial field strengths are found to lead to a loss of pressure confinement or “explosion” in the middle of the convection zone. The “explosion” draws outflow of high entropy plasma from the tube, and hence intensifies the submerged flux tubes to significantly super-equipartition field strengths (Section 7).

  • The evolution of active region scale flux tubes in the top layer of the solar convection zone, the process of sunspot formation at the surface, and the post-emergence evolution of subsurface flux tubes are not well understood (Section 8). It is suggested that active region magnetic flux on the photosphere becomes “dynamically disconnected” from its roots in the deep convection zone soon after the rising flux tube reaches the photosphere (see Section 8.3. With regard to flux emergence into the solar atmosphere, it has been shown that the magnetic buoyancy instability is a mechanism through which magnetic flux reaching the photosphere can expand dynamically into the stably stratified solar atmosphere. Recent 3D MHD simulations of magnetic flux emergence into the solar atmosphere and reconnection of emerging flux with a pre-existing coronal field are able to reproduce many observed features in newly emerging active regions. The observed X-ray sigmoids in the corona may be caused by the emergence of twisted magnetic flux tubes into the corona. (Section 8.2).