THE RISE OF ACTIVE REGION FLUX TUBES IN THE TURBULENT SOLAR CONVECTIVE ENVELOPE

The Sun's cyclic large-scale magnetic field with a period of 22 years is believed to be sustained by a dynamo mechanism. The Hale polarity law (Hale et al. 1919) of solar active regions indicates the presence of a large-scale subsurface toroidal magnetic field generated by the solar dynamo. The current prevailing picture is that the large-scale toroidal magnetic field responsible for the formation of solar active regions is amplified and stored at the base of the solar convection zone (e.g., Gilman 2000; Charbonneau 2010). Thus answering the question of how active region flux tubes rise from the base of the convection zone to the solar surface is vitally important in the development of solar dynamo theory.

A large body of calculations based on a so-called thin flux tube approximation has provided important insights into the dynamic evolution of Ω-shaped rising loops in the solar convective envelope (see review by, e.g., Fan 2009). The thin flux tube model assumes that the radius of the tube is significantly smaller than all other relevant scales of variation such as the pressure scale height and the curvature of the tube. Given the observed scale of flux ∼10²² Mx in a large solar active region, and assuming that the field strength of the flux tubes at the base of the convection zone is at least 10 kG, which is the order of field strength in equipartition with convection, one finds that the above condition for the thin flux tube approximation is satisfied in the bulk of the solar convection zone. The thin flux tube model solves for the mean motion of each tube segment under the influence of various integrated forces acting on the tube segment (see, e.g., Spruit 1981; Longcope & Klapper 1997). Results from thin flux tube models without the influence of convective flows suggest that the field strength of the toroidal magnetic field at the base of the solar convection zone needs to be in the range of about 30 kG to about 100 kG in order for the latitude of emergence and the tilt angles of the emerging loops to be consistent with the observed properties of solar active regions. If the field strength is ⩽20 kG, the poleward deflections of the trajectories of the rising flux tubes by the Coriolis force are too large such that the emerging latitudes are inconsistent with the observed sunspot latitudes (e.g., Caligari et al. 1995).

The thin flux tube calculations show that the Coriolis force acting on the rising Ω-loops produces asymmetries between the leading and the following legs of the loop (the leading leg is the proceeding portion of the loop in the direction of solar rotation), which provide explanations for several observed asymmetries between the leading and following polarities of bipolar active regions. First, a slight tilt of the emerging loop is produced, with the leading leg being closer to the equator than the following, consistent with the observed Joy's law of active region tilt angles (e.g., D'silva & Choudhuri 1993). Second, an asymmetric geometric shape of the emerging loop is produced, where the leading side is inclined at a smaller angle to the surface (e.g., Caligari et al. 1995). This may give rise to an apparent faster proper motion of the leading polarity as the loop emerges through the surface, which is observed (van Driel-Gesztely & Petrovay 1990). A third asymmetry is in the field strength of the loop, with a stronger field strength in the leading side compared to the following, which may explain the observed more coherent morphology of the leading polarity of an active region (e.g., Fan et al. 1993). However, this field asymmetry is found to depend on the assumed initial condition of the flux tube (Caligari et al. 1998; Fan & Fisher 1996).

The studies discussed previously do not address the dynamical effects of turbulent convection on a rising flux tube. These are shown to be adequate simplifications when the flux tube's magnetic field strength B ≳ (H_p/a)^1/2B_eq ∼ 3–5B_eq, where B_eq ∼ 10 kG is the field strength for which the magnetic energy density is in equipartition with the kinetic energy density of convection, H_p is the local pressure scale height, and a is the tube radius (Fan et al. 2003). Solar cycle dynamo models which take into account the dynamic effects of the Lorentz force from the large-scale mean fields suggest that the toroidal magnetic field generated at the base of the convection zone is ∼15 kG, and most likely cannot exceed 30 kG (e.g., Rempel 2006a, 2006b). Therefore, it is important to understand how toroidal flux tubes in the range of 15 kG to 30 kG rise through a turbulent solar convective envelope. Thin flux tube simulations which neglect the effect of convection (e.g., see review by Fan 2009) show that tubes in this range of field strengths tend to be deflected significantly poleward by the Coriolis force during their rise, and so have difficulty reproducing the emergence of active regions at low latitudes, as well as Joy's law. At these field strengths, convective downdrafts are capable of pinning down the flux tube while helical upflows between the downdrafts can boost the rise of the flux tube, and therefore the resulting emerging loops would be significantly different from those produced in the absence of convection.

A few studies have been performed to investigate the buoyant rise of a flux tube in a turbulent convective velocity field (e.g., Fan et al. 2003; Jouve & Brun 2009). Fan et al. (2003) carried out three-dimensional magnetohydrodynamic (MHD) simulations of the evolution of a buoyant magnetic flux tube in a stratified convective velocity field in a local Cartesian geometry without the effect of solar rotation. It is found that buoyant flux tubes with B ⩽ (H_p/a)^1/2B_eq ∼ 30 kG are strongly influenced by convection, where portions of the tubes in convective downdrafts become pinned down to the bottom of the domain, while the rise speed of sections within upflow regions is significantly boosted. It is found that the evolution of the flux tube is no longer sensitive to the twist of the tube in this convection dominated regime. Jouve & Brun (2009) have carried out the first set of global anelastic MHD simulations of the buoyant rise of an initially toroidal flux ring in a rotating, fully convective spherical shell, possessing self-consistently generated mean flows such as meridional circulations and differential rotation, representative of the conditions of the solar convective envelope. Due to the limited numerical resolution and the relatively high magnetic diffusivity in these global three-dimensional simulations, flux tubes with a very large initial field strength (ranging from 150 kG to 600 kG) and a large radius, corresponding to a total flux on the order of a few times 10²³ Mx, significantly greater than the typical active region fluxes, are considered. These large values are needed such that the rise times of flux tubes are less than their diffusive timescales. Because most of the cases considered are essentially in the magnetic buoyancy dominated regime, the rising toroidal flux tube only develops rather moderate undulations by the influence of the convective flows, and Ω-shaped tubes with undulations extending the depth of the convection zone are not found.

In this study, we use the thin flux tube model, with the inclusion of a (separately computed) three-dimensional turbulent convective velocity field in a rotating model solar convective envelope, to study its effects via the aerodynamic drag force, together with the forces of magnetic buoyancy, tension, and the Coriolis force, on the dynamic evolution of the emerging Ω-loops. Since the thin flux tube model preserves the frozen-in magnetic field condition, it does not suffer the erosion of magnetic buoyancy and the tension force due to artificial diffusion of the magnetic field, as do the Eulerian-based multi-dimensional simulations. Although these thin flux tube models do not capture the internal structure and the possible fragmentation of the flux tube, they provide an initial step toward understanding the effects of convective flows in relation to other forces acting on the rising flux tubes.

This paper is organized as follows. The equations of thin flux tube dynamics and the separately computed convective flow are discussed in Section 2. Results of our simulations are presented in Section 3, focusing in particular on the overall influence of convection on the flux tubes, the latitude of emergence and tilt angles of the flux tube, the asymmetry in field strength and geometry of the emerging flux loops, and the rise time for the flux tubes to reach the surface. We perform these thin flux tube simulations for a range of magnetic fields (15–100 kG) from near-equipartition to super-equipartition initial field strengths at varying latitudes from 1° to 40° in both the northern and southern hemispheres. A conclusion and discussion of the results is given in Section 4.

A discussion of the thin flux tube model, and how it is derived from the MHD equations under the assumption that the tube radius is thin compared to all other relevant scales of variation can be found in many previous publications (see, e.g., Spruit 1981; Longcope & Klapper 1997; Fan 2009). For this work, the equations that describe the evolution of the thin flux tube are

$$\begin{eqnarray} \rho {d {\bf v} \over dt} & = & {-}2 \rho ({\bf \Omega _0} \times {\bf v}) -(\rho _e - \rho) [{\bf g} - {\bf \Omega _0} \times ({\bf \Omega _0} \times {\bf r})]\nonumber\\ &&+ {\bf l} {\partial \over \partial s} \left({ B^2 \over 8 \pi } \right) + {B^2 \over 4 \pi } {\bf k} \nonumber \\ && - C_d {\rho _e | ({\bf v }-{\bf v}_e)_{\perp } | ({\bf v}-{\bf v}_e)_{\perp } \over (\pi \Phi / B)^{1/2} }, \end{eqnarray} \tag{ 1 }$$

$$\begin{equation} {d \over dt} \left({B \over \rho } \right) = {B \over \rho } \left[ {\partial ({\bf v}\,{\bm\cdot}\,{\bf l}) \over \partial s} - {\bf v}\,{\bm\cdot}\,{\bf k} \right], \end{equation} \tag{ 2 }$$

$$\begin{equation} {1 \over \rho } {d \rho \over dt} = {1 \over \gamma p} {dp \over dt}, \end{equation} \tag{ 3 }$$

$$\begin{equation} p = {\rho R T \over \mu }, \end{equation} \tag{ 4 }$$

$$\begin{equation} p + {B^2 \over 8 \pi } = p_e, \end{equation} \tag{ 5 }$$

where, r, v, B, ρ, p, T, which are functions of the time t and arc length s measured along the tube, denote, respectively, the position, velocity, magnetic field strength, gas density, pressure, and temperature of a Lagrangian tube segment, l ≡ ∂r/∂s is the unit vector tangential to the flux tube and k ≡ ∂²r/∂s² is the tube's curvature vector, subscript "⊥" denotes the component perpendicular to the flux tube, Φ = 10²² Mx is the constant total flux of the tube, ρ_e, p_e, T_e, and μ, which are functions of depth only, are, respectively, the pressure, density, temperature, mean molecular weight of the surrounding external plasma, g is the gravitational acceleration and a function of depth, Ω₀ is the angular velocity of the reference frame co-rotating with the Sun, with Ω₀ set to 2.7 × 10⁻⁶ rad s⁻¹ in this calculation, C_d is the drag coefficient, set to 1 in this calculation, γ is the ratio of specific heats, and v_e(r, t) (discussed below) is a time-dependent velocity field (relative to the rotating frame of reference) that impacts the dynamics of the thin flux tube through the drag force term. In the above equations, we do not introduce an explicit magnetic diffusion or kinematic viscosity term. However, we do introduce a drag force, which is the last term in Equation (1). This describes the interaction of the external fluid with the flux tube in a high Reynolds number regime (e.g., Batchelor 1967). The thin flux tube is untwisted, and discretized with 800 uniformly spaced grid points along its arc length s. The numerical method used to solve for the flux tube evolution as determined by the above set of equations has been described in detail in Fan et al. (1993).

For the stratification of the external field free plasma, namely ρ_e, p_e, T_e, μ, g, γ, and the super-adiabaticity, we use the reference solar model by Christensen-Dalsgaard et al. (1996) for the solar convection zone with an extension of a simple polytropic, sub-adiabatically stratified thin overshoot layer, as described in Fan & Gong (2000). Profiles of T_e, ρ_e, p_e, and the adiabaticity δ = ∇ − ∇_ad, where ∇ = dln T_e/dln p_e and ∇_ad is the value of ∇ one obtains by considering local adiabatic perturbations, are shown in Figure 1. The bottom left panel of Figure 1 shows the sub-adiabaticity of the thin overshoot layer, which extends from 4.8 × 10¹⁰ cm to 5.026 × 10¹⁰ cm. The bottom right panel shows the logarithm of the super-adiabaticity in the convection zone, which extends from 5.026 × 10¹⁰ cm to 6.75 × 10¹⁰ cm.

Zoom In Zoom Out Reset image size

The main new ingredient of the current thin flux tube simulations is the inclusion into the drag force term of an external time-dependent convective velocity field v_e(r, t) relative to the rotating frame of reference, computed separately from a three-dimensional global convection simulation using the anelastic spherical harmonic (ASH) code, as described in Miesch et al. (2006, hereafter MBT06). In the anelastic approximation, convective flows and thermal variations are treated as a linear perturbation to a background state which is taken from a one-dimensional solar structure model. The computed v_e(r, t) captures giant-cell convection and the associated mean flows (including the differential rotation and meridional circulation) in the rotating solar convective envelope spanning r = 0.69 R_☉ to r = 0.97 R_☉ (4.8 × 10¹⁰ cm to 6.75 × 10¹⁰ cm), resolved by a grid of 129 points in r, 256 points in θ, and 512 points in ϕ. ASH is a pseudo-spectral code with horizontal and vertical basis functions given by spherical harmonics and Chebyshev polynomials, each de-aliased by keeping only the lowest 2/3 of modes (maximum spherical harmonic degree ℓ_max = 170 and Chebyshev degree n_max = 86).

Simulation parameters and boundary conditions are similar to Case AB3 in MBT06. In particular, the lower thermal boundary condition is the same as in Case AB3, with a latitudinal entropy gradient imposed in order to implicitly capture thermal coupling to the tachocline: S(θ, r₁) = C_P(a₂Y₂₀ + a₄Y₄₀), where S is the specific entropy per unit mass, r₁ is the inner boundary, C_P is the specific heat at constant pressure, Y_ℓm(θ, ϕ) is the spherical harmonic of degree ℓ and order m, a₂ = 1.7 × 10⁻⁶, and a₄ = −0.43 × 10⁻⁶. This helps promote a conical rotation profile but the strong rotational shear is attributed to the convective Reynolds stress rather than baroclinicity (see MBT06). However, the radial entropy gradient imposed at the outer boundary is steeper, more in line with solar structure models (e.g., Christensen-Dalsgaard et al. 1996); ∂S/∂r = −10⁻⁵ erg g⁻¹ K⁻¹ cm⁻¹ in this case compared to −10⁻⁷ erg g⁻¹ K⁻¹ cm⁻¹ in Case AB3. This implies a higher Rayleigh number, which is partially offset by higher values of the turbulent viscosity ν and thermal diffusivity κ. In this simulation the values of ν and κ at the outer boundary (r = 0.97R) are 2 × 10¹³ and 4 × 10¹³ cm² s⁻¹, respectively, and each decreases with depth in proportion to the inverse square root of the background density $\hat{\rho }^{-1/2}$. The density contrast across the domain is 69, which corresponds to 4.2 density scale heights (somewhat larger than AB3, due to a slightly deeper simulation domain). This yields a mid-convection zone Rayleigh number R_a of 5 × 10⁶ and Reynolds number R_e of order 50. The Rayleigh number is defined here as R_a = gr²dΔS/(νκC_P) where g is the gravitational acceleration and d = r₂ − r₁ is the depth of the layer. The Reynolds number is given by V_rmsd/ν where V_rms is the root mean square velocity relative to the rotating reference frame. This value of R_a is somewhat larger than that used in the flux tube simulations of Jouve & Brun (2009), but R_e is about a factor of two smaller.

This simulation is more laminar than some others done with the ASH code (e.g., Miesch et al. 2008; Jouve & Brun 2009), and we acknowledge that this may have an impact on our results. However, the focus of this paper is the fundamental mechanisms by which thin flux tubes interact with global convection and mean flows. The convection simulation possesses all the relevant features necessary to investigate this interaction, including columnar, asymmetric, rotationally aligned cells at low latitudes (density-stratified banana cells), a rapidly evolving downflow network at higher latitudes in the upper convection zone, dominated by helical plumes, and a strong, solar-like differential rotation. The convective turnover timescale is τ ∼ 2 × 10⁶ s, or about 23 days. Our work provides a baseline for future work which can further assess how the details of the convective flow affect the results. Since even the highest-resolution simulations exhibit similar basic features, we do not expect the essential results to change significantly. The principle effect we expect at higher Rayleigh and Reynolds numbers is an increase in the random scatter due to stochastic turbulent fluctuations. Decoherence of large-scale motions and turbulent drag from small-scale motions could also slightly reduce rise times but we expect this effect to be minor since we believe that large-scale, columnar banana-like cells must persist even in highly turbulent parameter regimes in order to provide the requisite Reynolds stresses to account for the solar differential rotation.

Figure 2 shows a snapshot of the radial velocity of the giant-cell convection at a depth of 25 Mm below the solar surface. An associated movie showing the evolution of the convective flow pattern at this depth for over a period of about 315 days is also available in the electronic journal. The convective flow pattern shows broad upflow cells surrounded by narrow and intense downflow lanes. The maximum downflow speed in the convective envelope reaches nearly 600 m s⁻¹ at a depth of about 86 Mm below the surface. Throughout most of the convection zone, the combined influence of the Coriolis force and the density stratification induces anti-cyclonic vorticity in expanding upflows and cyclonic vorticity in contracting downflows. This yields a mean kinetic helicity density $H_k = \left<\mbox{\boldmath $\omega \cdot v$}\right>$ which is negative in the northern hemisphere and positive in the southern hemisphere (Figure 3(a)), where ${\bm \omega}$ is the vorticity of the convective flow with velocity ${\bm v}$. There is a weak sign reversal of H_k in the lower convection zone where downflows expand and recirculate, inducing anti-cyclonic vorticity. Such a helicity pattern is typical for rotating, compressible convection (e.g., Miesch & Toomre 2009).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At low latitudes there is a preferential alignment of elongated downflow lanes with the rotation axis, reflecting the presence of the so-called banana cells (Figure 2). These features propagate in a prograde sense relative to polar regions, due in part to the differential rotation and in part to an intrinsic phase drift akin to traveling Rossby waves (e.g., Miesch & Toomre 2009). Such structures dominate the convective Reynolds stress, maintaining a strong differential rotation comparable to that inferred from helioseismic inversions. In particular, the total angular velocity Ω/2π (with respect to the inertia frame) decreases monotonically from about 470 nHz at the equator to about 330 nHz at the poles and exhibits nearly conical contours at mid latitudes (see Figure 3(b)), as in the solar convection zone (Thompson et al. 2003).

As is in Fan & Gong (2000), our simulations start with toroidal magnetic flux rings in mechanical equilibrium (neutral buoyancy), located at a radial distance to the center of the Sun r = r₀ = 5.05 × 10¹⁰ cm, slightly above the base of the solar convection zone at r = r_czb = 5.026 × 10¹⁰ cm. Figure 4 shows a snapshot of the radial velocity at an arbitrary azimuthal angle ϕ, with the green dot representing the radius (to scale) of the largest flux tube, which occurs for magnetic field strengths of 15 kG. This figure shows the flux tube at its initial starting position compared to the base of the convection zone. Note that the convective velocity field is allowed to penetrate into the overshoot region. To ensure initial neutral buoyancy, the internal temperature of the flux tube is reduced compared to the external temperature. We consider a range of initial field strengths, with B₀ = 15 kG, 30 kG, 40 kG, 50 kG, 60 kG, and 100 kG, and initial latitudes λ₀ ranging from 1° to 40° for the toroidal flux ring. Considering the rms of the convective downflows from the ASH simulation at the base of the convection zone which are on the order of 35 m s⁻¹, the equipartition magnetic field is on the order of B_eq ∼ 5 kG. In this case, we are investigating flux tubes on the order of 3–20 B_eq. In all cases, the flux of the tube is constant at 10²² Mx, on the order of large-scale solar active regions. For this study, we perform two sets of simulations sampling different time ranges of the ASH convective flow at each magnetic field strength for each initial latitude in both the northern and southern hemispheres. The toroidal ring in neutral buoyancy is perturbed with small undular motions which consist of a superposition of Fourier modes with azimuthal order ranging from m = 0 through m = 8 with random phase relations. In this work, the external time-dependent three-dimensional convective flow described above also impacts the flux tube through the drag force term, and we study its effect.

Zoom In Zoom Out Reset image size

3.1. Overview of Influence of Convection on Rising Flux Tubes

Figure 5 shows snapshots of the rising flux tube in the absence of convection, developed from an initial toroidal flux ring at a latitude of 6°, at a time when its apex is approaching the top boundary, with initial magnetic field strengths of 15 kG, 40 kG, and 100 kG, respectively. The top three images show the flux tube as if looking down on the Sun from the north pole, whereas the bottom three images show the flux tube as if looking directly at the equator. For these cases in the absence of convection, emerging loops develop as a result of the nonlinear growth of the magnetic buoyancy instability of the initial toroidal flux tube, and the time for the flux tube to rise to the surface decreases with increasing magnetic field strength. The rise time is 6.4, 2.1, and 0.25 yr for 15, 40, and 100 kG, respectively. The rise of the loop becomes more radial with increasing initial field strength, with emergence latitudes of 233, 188, and 78 for 15, 40, and 100 kG magnetic fields, respectively. In Figure 5, we note the appearance of predominately m = 2 modes for magnetic field strengths of 40 kG and 100 kG. The 15 kG flux tube shows a predominately m = 1 mode, superposed with a very small amplitude m = 3 mode. The results are similar to those in Fan & Gong (2000), where the toroidal flux ring develops the undular buoyancy instability (e.g., Caligari et al. 1995), with m = 1 and m = 2 being the dominant unstable modes, and Ω-shaped rising loops form.

Zoom In Zoom Out Reset image size

Figure 6 shows snapshots of the rising flux tubes with the same initial conditions as Figure 5, but with the influence of the external convective flow included. It can be seen that with a low initial field strength of 15 kG, the development of the rising loops is largely controlled by the convective flows. Convection drastically decreases the rise time from 6.4 yr to only 0.2 yr. On the other hand, at a large magnetic field strength of 100 kG, we find that the development of the rising loops is still mainly controlled by the development of the magnetic buoyancy instability, with the strongest convective downdrafts producing some moderate perturbations to the final emerging loops. In this case, the presence of convection also shortens the rise time by about one month. At 40 kG, the influence of convection is somewhere in between the above two cases. Again convection reduces the rise time, from 2.1 yr to 0.5 yr in this case. Interestingly, with convection, the flux tube with 40 kG initial field takes longer to rise than either cases with 15 kG or 100 kG. The predominately m = 2 mode for the 100 kG flux tube is still discernible. However, at lower magnetic field strengths the flux tube evolution is altered by convective flows such that a dominant mode can no longer be established.

Zoom In Zoom Out Reset image size

The differential rotation of the convective velocity field does have an influence on the evolution of the rising flux tubes, especially at low magnetic field strengths. Figure 7 plots the average velocity in the azimuthal direction ϕ of the mass elements in the legs of the emerging flux loop once it has reached the top of the simulation domain. We find that without convection, the average velocity of the mass elements of the emerging loop is always negative. This occurs because, as the tube rises, the conservation of angular momentum drives a retrograde motion of particles in the tube, producing a negative velocity in the ϕ direction. At a large magnetic field strength of 100 kG, the averaged azimuthal velocity of the emerging loop in the presence of convection is found to be centered around the azimuthal velocity for the tubes without convection. With decreasing magnetic fields, the azimuthal velocity in the presence of convection becomes preferentially faster than in the case without convection. This indicates that differential rotation is attempting to push the elements of the flux tube prograde compared to the case without convection, especially at low magnetic field strengths where tubes are more susceptible to convective flows.

Zoom In Zoom Out Reset image size

Figure 8 summarizes how long it takes for an Ω-shaped loop to emerge at the top of the simulation domain. Without convection, the time for the flux tube to emerge decreases with increasing field strength. This agrees with Moreno-Insertis (1983), where they find the same trend for flux tubes given sinusoidal perturbations. Also, the range of variation in rise time between tubes initiating from different latitudes (as represented by the size of the error bar centered on the mean in Figure 8) decreases with increasing initial magnetic field. When the influence of convection is included, a different trend emerges, and the rise times are no longer monotonic with increasing magnetic field strength. In all cases, the rise time is reduced with the addition of convection, especially for the cases with weaker fields. The rise time becomes the longest for mid-field strength cases (i.e., 40 kG). At these mid-field strengths, the effects due to the magnetic buoyancy and the average convective downflows influencing the rise of the flux tube are of similar magnitudes, and work against each other, keeping the flux tube in the mid-convection zone for a longer period of time than at other magnetic field strengths. These effects will be discussed more in detail in the following paragraphs. For all the field strengths investigated, convection is found to reduce the variation of rise times for flux tubes with different initial latitudes, as can be seen in the size of the error bars compared to the case without convection.

Zoom In Zoom Out Reset image size

To understand the importance of convection, we compare the magnitude of the magnetic buoyancy force with that of the drag force from the external convective flows, following Fan et al. (2003). For the drag force to dominate buoyancy

$$\begin{equation} \frac{C_{D} \rho _{e} v_{e}^2}{\pi a} > \frac{B^2}{8\pi H_{p}}, \end{equation} \tag{ 6 }$$

where H_p is the local pressure scale height, and a is the flux tube radius a = (Φ/πB)^1/2, with flux of the tube Φ = 10²² Mx. The pressure scale height ranges from 5.6 × 10⁹ cm at the base of the convection zone to 6.9 × 10⁸ cm at the top of the simulation domain. Also, the radius of the flux tube spans from 4.6 × 10⁸ cm for an initial magnetic field strength of 15 kG to 1.8 × 10⁸ cm for 100 kG. For an order of magnitude estimate of the magnetic buoyancy here, we have assumed thermal equilibrium between the flux tube and the external fluid. Assuming 2C_D/π is of order 1, then Equation (6) simplifies to

$$\begin{equation} v_{e} > v_{a} \bigg (\frac{a}{H_{p}} \bigg)^{1/2}, \end{equation} \tag{ 7 }$$

i.e., in order for convection to dominate, the convective flow speed v_e needs to be greater than the Alfvén speed v_a = B/(4πρ_e)^1/2 multiplied by (a/H_p)^1/2.

In Figure 9, we have plotted as a function of depth the peak downflow, peak upflow, the rms of the downflow, and the rms of the upflow of the convection velocity field at each constant r surface. In comparison, we have also plotted the right-hand side of Equation (7) evaluated at the apex (portion of flux tube with largest r value) of each of the flux tubes shown in Figure 6 with different initial field strengths, as it traverses the convection zone. It can be seen that for the 100 kG case, only the largest convective downflows are strong enough to impact the rise of the flux tube, whereas at 15 kG, even the rms of the convective flows, both upflows and downflows, can significantly affect the rise of the tube. In the case of 40 kG, all of the convective velocity fields except for the rms upflows contribute to the development of the rising loops. These explain the general behavior we see in Figure 6. Note in Figure 9 that there is a concentration of points for the 40 kG (green) and 15 kG (red) cases. This is due to the fact that the flux tube is continually buffeted by convection, with strong downflows pushing down some apices such that new apices will take over for portions of the flux tube evolution. Often, the apices are also aided by boosts from convective upflows. Eventually, one loop becomes buoyant enough such that it will finally emerge at the top of the simulation domain. There is a jump in the 15 kG curve that occurs near 6.3 × 10¹⁰ cm because one loop is boosted by a convective upflow and becomes significantly more buoyant than the previous apex. It is clear from Figure 9 that downflows in the convective velocity field dominate in amplitude, however their spatial extent is small compared to the upflows as indicated by the narrow downflow lanes shown in Figure 2. Strong downflows can pin the flux tubes to the base of the convection zone at the beginning of the simulation faster than the Fourier mode perturbations can in the case without convection. Especially at low magnetic field strengths where the tube is highly susceptible to convection deformation, many portions of the tube will become anchored as compared to the case without convection. Also, rising loops can be significantly boosted by broad upflows such that they emerge at the top of the simulation domain much faster than in the case without convection, where only buoyancy aids in driving the tube to the surface. Convection can also enhance buoyancy instabilities by introducing finite-amplitude perturbations and subsequent gravity induced draining of fluid from the flux tube apex. Figure 10 shows snapshots (one individual time instance) of the radial distance r of the flux tube from Sun center as a function of the azimuthal angle ϕ (black lines), as well as the external radial velocity experienced by the tube at the height r of each tube segment (red lines). These snapshots are for the same flux tubes as used in Figures 6 and 9, with the left column showing the snapshot at a time when the apex of the tube has reached 0.82 R_☉, and the right column showing the snapshot for the last time step when the flux tube apex has reached the top of the simulation domain at 0.97 R_☉. From this figure, it is evident that at large magnetic field strengths, only the strongest downflows can perturb the tube. However, at small magnetic field strengths, all flows are capable of deforming the flux tube.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2. Latitude of Emergence and Tilt Angles

In the absence of convection, we find that for weaker initial field strengths, there is a low latitude zone devoid of flux emergence due to the poleward deflection of the rising tube by the Coriolis force (e.g., Choudhuri & Gilman 1987; Caligari et al. 1995). This is shown in Figure 11, where we have plotted the deflection of the flux tube in latitude (emergence latitude minus initial latitude) as a function of the flux tube initial latitude. The red diamonds (blue plus signs) show the deflection of the flux tubes without (with) the influence of convection. This low latitude zone void of flux emergence in the case without convection shrinks as the initial magnetic field strength is increased. This is because at larger initial magnetic field strengths, the buoyancy force overpowers the Coriolis force such that tubes rise more radially. However, when the flux tubes are subjected to the convective flow, we find that they are able to emerge at lower latitudes near the equator, even for low field strength cases. In some cases, the convective velocity field is able to push the flux tube such that its apex will emerge closer to the equator than where the tube initially started. This is indicated in Figure 11 by negative plus sign points. We note that with convection, the poleward deflection of the flux tubes is reduced for the majority of cases. Flux tubes subjected to turbulent convection spend less time in the convection zone than the flux tubes without convection (see Figure 8). As such, the Coriolis force will have less time to deflect these flux tubes to higher latitudes. Also, increased anchoring of weaker flux tubes by convection reduces the fraction of the tube that moves away from the rotation axes, and therefore reduces the effect of the Coriolis force on the tube.

Zoom In Zoom Out Reset image size

In most cases without convection, the variation of the emergence latitude with initial latitude is smooth, except for the 15 kG flux tube, where there is an abrupt 14° jump in emergence latitude between the tubes originating at 5° and 10° (see diamond points in the bottom right panel of Figures 11 and 12). At this magnetic field strength in the absence of convection, flux tubes originating at 5° and 1° have troughs which become anchored to the overshoot region below the base of the convection zone. However, for the same tubes originating from 8° and above, the rising tubes do not become anchored. As such, the tube floats upward as a whole, and experiences a relatively greater deflection toward the pole and hence emerges at a higher latitude compared to the anchored cases. In these cases, the m = 1 to m = 3 modes grow slower than the m = 0 mode. This corresponds to a poleward slip of the tube as a whole. Thereafter, modes on the order of m = 1 to m = 3 begin to develop larger amplitudes, allowing loops to develop, though the troughs do not anchor at the base of the convection zone. With convection, at this same magnetic field strength, all of the flux tubes have troughs which become anchored because convective downflows pin them to the bottom of the convection zone from the beginning of the simulation. We note that the flux tubes tend to anchor (when they anchor at all) around 5.0 × 10¹⁰ cm in both the cases with and without convection. At this depth, the sub-adiabaticity is δ = −1.3 × 10⁻⁴. The flux tubes all originate at 5.05 × 10¹⁰ cm, 2.8 × 10⁸ cm above the base of the convection zone, and have a super-adiabaticity here of δ = 3.2 × 10⁻⁸, or very nearly zero. Reducing the sub-adiabaticity would reduce the chances for anchoring of the portions of the flux tube that are perturbed into the overshoot region. Increasing the sub-adiabaticity would help anchor these flux tubes if they were perturbed into the overshoot region and did not remain anchored at a reduced sub-adiabaticity. However, we find that the tubes which do not anchor in this study are never perturbed into the overshoot region. Instabilities of higher order with large enough amplitude to perturb the tube into the overshoot region do not develop fast enough to overcome the m = 0 mode.

Zoom In Zoom Out Reset image size

Figure 12 shows the tilt angle of the emerging flux loops as a function of the emergence latitude, for flux tubes with a range of initial latitudes (from 1° to 40° as indicated by the color of the points). The different panels in the figure show the results for initial magnetic field strengths of 15 kG, 30 kG, 40 kG, 50 kG, 60 kG, and 100 kG, respectively. In the plots, results with (without) the influence of convection are shown as plus signs (diamonds).

It is well known that bipolar active regions on the Sun tend to form with the leading polarity (in the direction of solar rotation) being slightly closer to the equator than the following polarity, and thus showing a tilt angle for the line connecting the following to the leading polarities with respect to the east–west direction, as described by Joy's law (e.g., Hale et al. 1919; Dasi-Espuig et al. 2010). On average, the tilt angle of the active region increases with its latitude. In Figure 12, the tilt angle is computed as the angle between the tangent vector at the apex of the emerging loop and the local east–west direction. We define a positive sign of tilt as a clockwise (anticlockwise) rotation of the tangent vector away from the east–west direction in the northern (southern) hemisphere, consistent with the direction of the observed mean tilt of active regions. If on the other hand the magnitude of the tilt angle exceeds 90°, then it means that there is an anti-Hale polarity arrangement for the active region.

It is important to note that we do not consider the effect of magnetic field line twist on the evolution of flux tubes. Fan (2008) considers the effect of twist on the tilt angle of emerging flux tubes for a three-dimensional MHD simulation in a local Cartesian geometry without an external velocity field. They find that in order for emerging flux tubes of highly super-equipartition field strength on the order of 10⁵ G to have tilts consistent with Joy's law, the initial twist rate of the flux tube needs to be smaller than about a half of that required for the tube to rise cohesively. However, in the case of flux tubes in a convecting box with equipartition magnetic field strengths on the order of 10⁴ G, it is found that convective downflows dominate over magnetic buoyancy effects, and the evolution of the tube is no longer sensitive to the initial twist of the tube (Fan et al. 2003). In this paper, we are primarily concerned with how convection alters the overall evolution of the rising flux tubes, and the thin flux tube model used here does not address the effect of the field line twist.

Dasi-Espuig et al. (2010) analyzed Mount Wilson and Kodaikanal sunspot group tilt angle and latitude data spanning solar cycles 15–21. In order to rederive Joy's law, they performed a fit for the equation: α = mλ to the data, where α, m, and λ represent the tilt angle, the slope, and the latitude, respectively. The fit is forced to go through zero because no tilt is expected for equatorial sunspot groups. They obtained a slope of 0.26 ± 0.05 and 0.28 ± 0.06 for Mount Wilson and Kodaikanal data, respectively.

In Figure 12, we show the same kind of linear least-squares fit for both the simulated tilt angles with and without the influence of the convective flow. The slope of the best-fit lines with uncertainties for 100 kG, 60 kG, 50 kG, 40 kG, 30 kG, and 15 kG without convection are, respectively: 0.25 ± 0.01, 0.30 ± 0.02, 0.22 ± 0.02, 0.12 ± 0.01, −0.07 ± 0.05, and −0.13 ± 0.13. Thus the slopes of the best-fit lines fall within the range found by Dasi-Espuig et al. (2010) only for flux tubes with initial field strengths of 100 kG, 60 kG, and 50 kG. For weaker initial field strengths (below about 40 kG) the emerging loops begin to show negative tilt angles, opposite to the sign of the active region mean tilts. This occurs because flow along the flux tube at the loop apex changes from diverging into converging as it enters the upper layers of the convection zone (e.g., Caligari et al. 1995; Fan & Fisher 1996), and the Coriolis force acting on the converging flow drives a tilt of the wrong sign. With the influence of the convective flow included, the slopes of the best linear fits to the tilt angles all increase due to the systematic effect of the kinetic helicity (see Figure 3(a)) associated with the rotating convective flows, which tends to drive rotation at the apex of the loops in the direction consistent with the sign of the observed active region mean tilts. The slope of the best-fit lines with convection for initial magnetic field strengths of 100 kG, 60 kG, 50 kG, 40 kG, 30 kG, and 15 kG are, respectively: 0.30 ± 0.08, 0.34 ± 0.05, 0.34 ± 0.12, 0.44 ± 0.14, 0.12 ± 0.18, and 0.15 ± 0.21. As a result, the tilt angles of emerging loops with initial field strengths of 100 kG, 60 kG, 50 kG, and 40 kG are all consistent with Joy's law. These field strengths all show a positive Joy's law trend, given the uncertainties in the fitted slopes. However, at lower field strengths, the tilt angles show large random scatters produced by convection, and do not show a significant systematic dependence on latitude as described by Joy's law. The uncertainty of the slope also tends to increase with decreasing magnetic field, reflecting the random scatter produced by convection. While the slopes of the best-fit lines for 30 kG and 15 kG are still positive, the uncertainties are too large to report a definitive positive Joy's law trend for these field strengths. There is a larger than expected uncertainty on the slope for 100 kG. This occurs because there is an outlier, which cannot be seen in Figure 12 as it falls outside the region graphed. The outlier corresponds to a flux tube which originates in the southern hemisphere, but emerges in the northern hemisphere such that the emergence of the tube would result in an anti-Hale polarity arrangement for the active region.

3.3. Asymmetry in Field Strength

Another well-known observed asymmetry is in the morphology of the leading and the following polarities of an active region, where the leading polarity flux tends to be concentrated into a well-formed sunspot, whereas the following polarity flux tends to appear more fragmented and dispersed (e.g., Bray & Loughhead 1979). The thin flux tube simulations of Fan et al. (1993) showed that the preceding leg of an emerging magnetic flux loop has a stronger magnetic field than the following leg as a result of the differential stretching of the rising loop due to the Coriolis force. This gives an explanation for the observed more coherent and less fragmented morphology for the leading polarity flux in an active region. However, subsequent simulations (e.g., Caligari et al. 1995; Fan & Fisher 1996) using the mechanical equilibrium initial state (which is more physical) as opposed to the temperature equilibrium as used by Fan et al. (1993) found that the leading leg of the emerging loop has a stronger magnetic field than the following only for flux tubes with an initial field strength that is below about 60 kG. For flux tubes with higher initial field strengths, the field strength asymmetry reverses at the top of the loop.

Here we investigate this magnetic field asymmetry by calculating dB/ds at the apex of the emerging flux loop as shown in Figure 13. If dB/ds is greater (less) than zero, then the leading (following) leg has a stronger magnetic field. As can be seen in Figure 13, at 50 kG and below, the majority of the emerging loops show stronger field in the leading leg than the following leg in the presence of convection. However, at 100 kG, with convection, only some of the emerging flux loops at higher latitudes have stronger field in the preceding leg. At 60 kG there are about equal number of positive and negative dB/ds cases. Thus with convection, loops with initial field ⩽50 kG tend to emerge with the appropriate magnetic field asymmetry (i.e., with a stronger field in the preceding leg and therefore can be expected to result in an emerging active region with a more coherent leading polarity).

Zoom In Zoom Out Reset image size

3.4. Asymmetry of Inclination

When looking down on the emerging loops from the pole, it is apparent that there is an asymmetry in the inclination of the leading and following legs with respect to the vertical direction, as is evident in the cases without convection shown in Figure 5. The following leg tends to have a steeper slope than the leading leg. This asymmetry is caused by the Coriolis force and the conservation of angular momentum as the tube rises through the convection zone (Moreno-Insertis et al. 1994; Caligari et al. 1995, 1998), and provides an explanation for the apparent asymmetric east–west proper motions of the two polarities of an emerging active region. As the tube rises above a constant r surface, inclination of the legs of the loop causes an apparent more rapid motion of the leading polarity spot away from the emerging region as compared to the motion of the following polarity spots.

Here we quantify the steepness of each leg of the emerging loop by calculating what we call the inclination angle. To do this, first we find the portion of the emerging loop which is concave downward. Then we find the best-fit line for each leg from the flux tube apex in the concave down portion. The angle between the best-fit line for each leg and $- {\hat{\bf r}}$ is the inclination angle. The smaller this angle, the steeper the slope of the leg. In Figure 14, we have plotted the average difference of the inclination angles (leading leg minus the following leg). Only the average inclination difference for each magnetic field strength is plotted because we find that the inclination difference does not vary systematically with latitude. The bars are the standard deviation, showing the spread in inclination angle dif- ference for each magnetic field strength. This spread is due primarily to convective effects. Positive inclination difference means the following leg is steeper, consistent with the results of Moreno-Insertis et al. (1994) and Caligari et al. (1995, 1998). We find that at all field strengths, the majority of the emerging loops develop a steeper slope for the following leg. The inclination angle differences are overall reduced with the inclusion of convection.

Zoom In Zoom Out Reset image size

It is worthwhile to note that convection aids in bringing the footpoints of the rising flux loops closer together, most notably in the 15 kG and 40 kG cases as can be seen in Figures 5 and 6. This can also be observed in Figure 10, where strong convective downflows pin portions of the tube down, creating many buoyant loops with very close footpoints at 15 kG.

In subjecting the thin flux tube model to a turbulent solar convective velocity field computed separately from a three-dimensional global convection simulation in a rotating spherical shell representing the solar convective envelope, we study how convection can influence the development and emergence of solar active region flux tubes. Although idealized, the thin flux tube approximation allows us to investigate flux emergence at moderate to strong field strengths (15–100 kG) in localized flux concentrations under perfect flux frozen-in conditions. We find that as the field strength of the initial toroidal flux ring placed near the bottom of the solar convective envelop is increased from 15 kG to 100 kG, the development and the evolution of the emerging flux loops change from being convection dominated to magnetic buoyancy dominated. At 15 kG initial field strength, the development of the emerging flux loops is largely controlled by the convective flows. Both the average upflows and downflows can significantly affect the rise of the flux tube. On the other hand, at 100 kG initial field strength, the development of the emerging loops is largely controlled by the growth of the magnetic buoyancy instability and only the strongest convective downdrafts can significantly impact the rising tubes.

With the inclusion of the convective flows, the previous issue of 15–30 kG flux tubes being significantly deflected poleward by the Coriolis force during their rise is resolved. Loops can now emerge at low latitudes all the way to the equator. The main reasons for this are the anchoring of the emerging flux loops by the convective downdrafts and also the faster rise of the loops propelled by the upflows. Convection is also found to produce random scatters in the tilts of the emerging loops, especially for the weaker field strength cases where the scatters are greater. In addition, because the convective flow in the rotating spherical envelope shows a mean kinetic helicity that is negative (positive) in the northern (southern) hemisphere, it on average tends to drive clockwise (anticlockwise) tilts for the rising loops in the northern (southern) hemisphere, consistent with the sign of the active region mean tilts. As a result, the inclusion of convection tends to increase the slope of the linear Joy's law fit of the tilt angles of the emerging loops as a function of emerging latitudes. Flux tubes with initial field ≳40 kG and ≲100 kG are found to produce a mean tilt angle dependence on latitude that is consistent with the observed Joy's law. However, the effect of the kinetic helicity in the convective flow is not sufficient, on average, to correct the tendency for loops to develop tilts of the wrong sign (opposite to that of the observed active region mean tilt) for tubes with initial field strength <40 kG.

Similar to Caligari et al. (1995) and Fan & Fisher (1996), we find that at or below 50 kG, the leading leg of the emerging loop tends to have a larger magnetic field than the following, which may provide an explanation for the observed better cohesion of the leading polarity of an emerging active region compared to the following polarity. This trend of asymmetry in field strength reverses for tubes with initial field >60 kG. Also, we find that the inclusion of convection does not change the general trend of the asymmetry in the geometry of the emerging loop found in previous studies (e.g., Caligari et al. 1995), where the leading leg tends to be more horizontally inclined than the following, which would provide an explanation for the apparent asymmetric east–west proper motions of the two polarities of an emerging active region. This trend of geometric asymmetry is not dependent on the magnetic field strength.

These results combined suggest that in order for the emerging loops to be consistent with the observed properties of solar active regions, the best range of field strengths for toroidal magnetic field at the bottom of the convection zone is ∼40–50 kG. Above this range, the magnetic field is strong enough to resist convection such that the flux tube properties are not changed significantly from the case without convection, but the magnetic field asymmetry of the emerging loops tends to be of the wrong sense. Below this range, convection dominates the development and evolution of the emerging loops. The loops are able to emerge at low latitudes, but their tilt angles show very large scatter. Increasing the number of simulations run to collect more data may result in a reduction of the uncertainty on the slope of the Joy's law best-fit lines, such that 15 kG and 30 kG flux tubes may also show a positive Joy's law trend.

The main limitation of the thin flux tube approach is that the model assumes the tube remains cohesive, and does not address the dynamic effect of the internal twist of the flux tube. Also, as the ASH convective flow is computed separately, the calculations do not address how the rising flux tube might affect the surrounding convective flow. Full three-dimensional spherical-shell MHD simulations (e.g., Jouve & Brun 2009) are necessary to address the above questions, but it is still beyond the current computational capabilities to adequately resolve an active region scale flux tube in a global convection simulation. It is important to note that these results are for flux tubes which have reached 21 Mm below the solar surface. Even though these flux tubes still must traverse a small portion of the convection zone and encounter the solar surface shear layer, we feel that the behavior of the flux tubes at this depth may give a good representation of the large-scale pattern of flux emergence at the solar surface. This paper serves as an initial study to understand how active region scale flux tubes of various strengths are affected by convective flows as they rise through the solar convective envelope.

This work is supported in part by NASA SHP grant NNX10AB81G to the National Center for Atmospheric Research (NCAR). NCAR is sponsored by the National Science Foundation. We also thank the anonymous referee for the very thorough report and helpful comments.