DOUBLE-DIFFUSIVE INSTABILITIES OF A SHEAR-GENERATED MAGNETIC LAYER

The generation of magnetic field by velocity shear and the field's subsequent evolution are of great importance to an understanding of the operation of the solar dynamo. While the current dynamo paradigm contains many complex interacting components (see, e.g., Parker 1993; Rempel 2006), an integral part of all current large-scale solar dynamo models is the creation of strong toroidal magnetic structures in the tachocline, the Sun's thin region of strong radial differential rotation separating the latitudinally differentially rotating convection zone and the solid-body-rotating radiative interior. Strong toroidal magnetic field is thought to be induced by the stretching action of the differential rotation on any background poloidal field (the Ω-effect of mean-field dynamo theory; Steenbeck et al. 1966). Subsequently, magnetic buoyancy instabilities (Parker 1955) of the generated field are invoked as the mechanism for the creation of distinct magnetic structures and their subsequent rise toward eventual emergence at the solar surface as active regions.

In recent papers, Vasil & Brummell (2008, 2009) (hereinafter VB1, VB2) show that, when naively using the commonly accepted paradigm for the Ω-effect, it is surprisingly difficult to initiate magnetic buoyancy instabilities if thermodynamic adjustments are assumed to be adiabatic. Specifically, they concluded from analytic (VB2) and numeric (VB1) calculations that for magnetic buoyancy instabilities to operate in a thin tachocline, the velocity shear flow imposed to drive the system must be necessarily hydrodynamically unstable. Roughly speaking, VB2 found that Ri ≃ Δz/H_p is required for magnetic buoyancy, where Ri is the Richardson number, Δz is the vertical shear width, and H_p is the local pressure scale height. For the solar tachocline, the estimated Richardson number is very large (Ri ≃ 10³–10⁵) but the region is thin so that Δz/H_p < 1 (Gough 2007) and so the above condition for instability is therefore unlikely to be satisfied. The simulations of VB1 confirmed numerically that hydrodynamically unstable imposed shear flows could generate magnetic buoyancy instabilities in a thin shear region, but stable shear flows did not. Energetically, this can be thought of as the following conundrum; for magnetic buoyancy to occur, the shear flow must transfer enough energy into a toroidal magnetic field for it to overcome the stable background stratification. If the shear can only build (through an Alfvénic process) magnetic field to the level of equipartition with the flow, and the shear is constrained by the stratification (since it is hydrodynamically stable), then it is difficult for the shear-induced magnetic field to overcome the constraints of the stratification.

It is the case, however, that near the bottom of the convection zone the ratio between the magnetic and thermal diffusivities, ζ = η/κ, is very small. It has long been recognized (Gilman 1970; Acheson 1979) that, in such circumstances, instabilities can occur that rely on the much greater diffusion rate of the stabilizing thermal component (see Hughes (2007) for further discussions). Further, VB2 noted that instability might be enhanced by double-diffusive effects. Assuming isothermal (as opposed to adiabatic) adjustments VB2 obtained the much less stringent requirement that ζRi ≃ Δz/H_p. It is therefore possible that a more-solar-like (high Ri) shear-generated magnetic field can become buoyantly unstable in a double-diffusive manner when ζ ≪ 1.

In this Letter, we consider this possibility in a convectively stable atmosphere. We show that if ζ is sufficiently small, and the background Alfvénic timescales are sufficiently slow, then double-diffusive instabilities can exist for parameters that are relevant to the solar interior.

We consider a Cartesian domain (x, y, z), where z is depth, x is the toroidal/zonal direction, and (y, z) are the poloidal directions. Our system consists of an initially vertical uniform magnetic field, B_z,0, permeating a stratified layer of compressible fluid under the influence of a forcing designed to generate a target shear flow $\mbox{\boldmath {$u$}}=(U_{0}(z),0,0)$. This system is similar to that described in VB1 and in Silvers et al. (2009). We solve standard non-dimensionalized equations for forced compressible magnetohydrodynamics governing the evolution of velocity $\mbox{\boldmath {$u$}}=(u,v,w)$, magnetic field $\mbox{\boldmath {$B$}}=(B_{x},B_{y},B_{z})$, density ρ, temperature T, and pressure P based around a polytropic atmosphere, T₀(z) = 1 + θz,  ρ₀(z) = T₀(z)^m,  P₀(z) = T₀(z)^m+1. The non-dimensionalization uses T_*, the temperature at the top of the domain, ρ_* which is proportional to the total mass within the domain, P_* = (c_p − c_v)T_*ρ_*, the fiducial pressure (given the specific heat capacities, c_p, c_v), B_z,0, the imposed background magnetic field strength, the layer depth d, and time units of τ_* = d ρ^1/2_*/P^1/2_*. We impose stress-free velocity and vertical magnetic field boundary conditions at z = 0, 1 together with T(z = 0) = 1 and ∂_zT(z = 1) = θ. The horizontal directions are periodic.

The system is solved in a domain with aspect ratio 2:1:1 in a similar manner to that in VB1 and in Silvers et al. (2009) at a resolution of (roughly) 256 × 256 × 512. This aspect ratio is specified to allow for certain expected dynamics, e.g., we anticipate three-dimensional modes to have smaller scales in the y and z directions than in the x direction. The thinness of the forced shear region compared to the local scale height is the geometrical factor that identifies the model with the tachocline. The kinetic and magnetic Reynolds numbers for our calculations are roughly Re ∼ 2000, Re_M ∼ 1000, respectively (based on the forced shear magnitude and width). Estimates of the numerical degrees of freedom (Davidson 2004) imply that our resolution is reasonable. However, we further satisfied ourselves that the results were consistent at varying resolutions.

The important parameters in the problem are the dimensionless thermal diffusivity, C_K = Kτ_*/ρ_*c_pd², the Prandtl number σ = μc_p/K (μ representing dynamic viscosity), the inverse Roberts number ζ = ηc_pρ_*/K = η/κ, and α = B²_z,0τ²_*/μ₀ρ_*d² which gives a measure of the Alfvén speed along the background field in terms of the fundamental acoustic velocity scale. Note that ζ is the critical parameter governing double-diffusive instability (Hughes & Weiss 1995).

We add a forcing term, $\mbox{\boldmath {$F$}}=-\sigma {C_{K}}\partial _{z}^{2}U_{0}\mbox{\boldmath {$\hat{x}$}}$, in the x-momentum equation that would (in the absence of magnetic effects and instabilities) maintain a desired target velocity $\mbox{\boldmath {$u$}}=(U_{0}(z),0,0)$ against viscous decay. We choose

$$\begin{equation} U_{0}(z)=\frac{1}{20}\tanh \left[10\left(z-\frac{1}{2}\right)\right] \end{equation} \tag{ 1 }$$

to mimic the smooth radial shear transition believed to exist in the solar tachocline. The width is chosen sufficiently narrow that ∂_zU₀ ≈ 0 (to within numerical precision) at the boundaries.

An important derived parameter is the Richardson number,

$$\begin{eqnarray} Ri&=&\min _{0\le {z}\le 1}\left(\frac{N_{0}(z)^{2}}{[\partial _{z}U_{0}(z)]^{2}}\right), \end{eqnarray} \tag{ 2 }$$

where N₀(z) = [ − g ∂_zln (P₀(z)^1/γ/ρ₀(z))]^1/2 is the local Brunt–Väisälä frequency of the stable background atmosphere, ∂_zU₀(z) is the local turnover rate of the background shear, and γ = c_p/c_v = 5/3. Ri measures the relative tendency of a shear flow to overturn fluid vertically compared to gravity's tendency to restore it to its original position. A large Ri implies that gravity is strongly stabilizing, whereas a small value can mean that shear instabilities are possible (Drazin & Reid 2004). The goal in this work is to obtain a magnetic buoyancy instability at high Richardson number, a result not found in VB1, but here we investigate the low ζ regime where the effects of thermal stability are severely reduced.

While the Richardson number provides a useful rough measure of the stability properties of our system, we also define the respective thermal and magnetic Rayleigh numbers

$$\begin{equation} R_{T}(z)=-\frac{N_{0}(z)^{2}}{\nu (z)\kappa (z)} \end{equation} \tag{ 3 }$$

$$\begin{equation} R_{B}(z)=\frac{g\alpha \,B_{x}(z)^{2}}{\nu (z)\kappa (z)P_{0}(z)}\frac{d}{dz}\ln {B_{x}(z)}, \end{equation} \tag{ 4 }$$

where κ(z) = C_K/ρ₀(z) and ν(z) = σκ(z). These definitions are consistent with those given in Hughes & Weiss (1995) except that here the magnetic Rayleigh number, R_B, is given in a form that anticipates three-dimensional instabilities (Newcomb 1961). For the direct (non-oscillatory) type of instability, of interest here, both Rayleigh numbers are negative, and this corresponds to a thermal stratification that is stabilizing and magnetic gradient that is destabilizing. A Rayleigh number based on the total density stratification would correspond to R_total = R_T − R_B. The key result of VB2, derived under the adiabatic thermal adjustment assumption, is equivalent to saying that R_total is always negative (ostensibly stable) unless Ri ≪ 1. In the double-diffusive context, where thermal adjustments are closer to isothermal, the critical parameter is actually R_DD = R_T − ζ⁻¹R_B, where ζ ≪ 1. It is R_DD that is relevant to the stability of our system and it can potentially become positive even if Ri ≫ 1. This is equivalent to the isothermal result derived in VB2 that required ζRi ≃ Δz/H_p.

For the basic thermal parameters (see Table 1), we take the polytropic index m = 1.6 that enforces a convectively stable background polytropic stratification, and a non-dimensional lower-boundary heat flux of θ = 5. With the chosen shear flow, this produced a minimum Richardson number of Ri = 2.96. While this value is significantly lower than we expect for the solar tachocline, it is high enough to guarantee that U₀(z) hydrodynamically stable. This fact was checked numerically with a purely hydrodynamic run.

In order to see if double-diffusive instabilities can exist, we ran two simulations that differ only in the thermal diffusivity κ = K/c_pρ_* (varied through its non-dimensional counterpart C_K). Thus ζ and σ are varied commensurately to maintain fixed viscous and magnetic diffusivities μ = σC_K and η = ζC_K. We keep σC_K = 2.5 × 10⁻⁶ and ζC_K = 5.0 × 10⁻⁶ for all our cases. For our first case, C1, anticipating instability, we choose ζ = 5.0 × 10⁻⁴, σ = 2.5 × 10⁻⁴, and C_K = 0.01. In our second case, C2, anticipating stability, we choose ζ = 0.01, σ = 5.0 × 10⁻³, and C_K = 5.0 × 10⁻⁴. Thus, ζ in C1 is 20-fold smaller than in C2.

The parameter α = σζC²_KQ (where Q is the Chandrasekhar number measuring the ratio of magnetic and diffusive timescales) sets the background vertical magnetic field strength. For C1 and C2 we set α = 1.25 × 10⁻⁵. We desire α ≪ 1 so that Alfvén timescales are slow compared to the acoustic timescale, and Q ≫ 1 so that Alfvén timescales are fast compared to diffusion. C1 and C2 have Q = 1.0 × 10⁶. It turns out that the Alfvén timescale of the background magnetic field has important consequences for the dynamic stability of our evolving MHD configuration. To investigate this aspect, we therefore present a third case, C3, that is identical to C1 except that α = 5.0 × 10⁻⁵, (Q = 4.0 × 10⁶). C3 therefore has two-fold faster Alfvén timescales.

In all three simulations, we begin with a purely vertical magnetic field, and allow the induction of a toroidal field layer by the action of the imposed velocity forcing. In the absence of two- and three-dimensional effects, this induction process is governed by a one-dimensional (mean) set of MHD equations:

$$\begin{equation} \rho _{0}(z)\partial _{t}u=\alpha \partial _{z}B_{x}+\sigma {C_{K}}\partial _{z}^{2}u-\sigma {C_{K}}\partial _{z}^{2}U_{0} \end{equation} \tag{ 5 }$$

$$\begin{equation} \partial _{t}B_{x}=\partial _{z}u+\zeta {C_{K}}\partial _{z}^{2}B_{x}. \end{equation} \tag{ 6 }$$

Strictly speaking, the density profile in Equation (5) will evolve along with the shear flow and magnetic field. However, unlike the work in VB1-2 and Silvers et al. (2009), the target shear flow here, U₀(z), has a large Richardson number. To an excellent approximation, the background density maintains its initial polytropic profile ρ(t, z) = ρ₀(z) in this initial phase. The evolution, as shown in Figure 1, of the mean toroidal magnetic field is therefore essentially independent of the thermal properties of the system, unless two- or three-dimensional instabilities occur.

Zoom In Zoom Out Reset image size

As time progresses, the peak field grows in the region of maximum forced velocity shear, and the field gradients strengthen. The ultimate question is whether a weak shear (Ri > 1) can induce strong enough gradients for a magnetic buoyancy instability to occur. Below, we show that two initially identically evolving one-dimensional MHD configurations can have completely different stability properties, based solely on the magnitude of the thermal diffusivity.

In C1, the system does indeed initially evolve according to Equations (5) and (6). However, since ζ = 5.0 × 10⁻⁴ is small enough, the induction of mean toroidal field only proceeds for a finite time before an instability becomes apparent at t ≈ 88. Figure (2) shows volume-rendered images at this time of the vertical velocity w(x, y, z) and the fluctuating toroidal magnetic field $B_{x}(x,y,z)-\overline{B}_{x}(z)$. In the initial induction phase (governed by Equations (5) and (6)) both these quantities are zero. However, Figure (2) clearly shows that a wave-like perturbation in the vertical velocity and toroidal field has appeared in the region of strong toroidal magnetic field gradients near the top of the evolving toroidal magnetic layer. The instability appears in a quasi-two-dimensional manner with a wavevector primarily in the y–z plane, taking a form that is similar to a classical "interchange" instability (Cattaneo & Hughes 1988). Roll-like motions principally swap lines of toroidal magnetic field that pierce y–z planes without a large degree of bending along the toroidal x-direction. There is a high degree of correlation between vertical velocity and toroidal field perturbations, strongly implying a buoyancy-driven instability (also see Figure 3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We now compare C2 results with those of C1, which differs only the value of ζ. We estimate the value of ζ for C2 to be sufficiently large to render R_DD everywhere negative, and therefore anticipate that the instability found in C1, if double-diffusive, should not appear in this case. We run C2 for over twice the time that it takes for the instability in C1 to manifest (up to t ≈ 180) and the system simply continues to evolve according to the mean Equations (5) and (6). Eventually, these dynamics are benignly disrupted by Alfvénic processes without any buoyancy instabilities (as explained in VB2) and we therefore halt the computation. We conclude that the instability discovered in C1 is of a double-diffusive nature, since it depends critically on a sufficiently large thermal diffusivity (small ζ).

A subtle issue in this problem is the influence of the strength of the imposed background field. As argued in VB1-2, this essentially provides a timescale for disruptive Alfvénic processes. For instability to occur, the necessary conditions must be met before Alfvénic processes can disrupt the source. In VB1-2, this required a strong velocity shear (hydrodynamically unstable, Ri ≈ 0.03) in order to induce the required magnetic gradients sufficiently quickly for a regular (mean-density-deficit driven) magnetic buoyancy instability to occur. Things are more complex in the high Ri, double-diffusive regime considered here. To elucidate this issue, we examine case C3 where the imposed background field strength is twice that of C1 and C2. Figure (3) shows a vertical profile of the mean vertical transport of B_x for both C1 and C3, a good indicator of the existence and efficiency of any magnetic buoyancy instabilities. The measure plotted is

$$\begin{eqnarray} C(z)&=&\frac{\overline{w{B_{x}}}-\overline{w}\overline{B}_{x}}{\max _{z}|w_{\mathrm{rms}}B_{x,\mathrm{rms}}|}, \end{eqnarray} \tag{ 7 }$$

where the overline represents a horizontal average, and $w_{\mathrm{rms}}(z)=\overline{w^{2}}-\overline{w}^{2}$, $B_{x,\mathrm{rms}}(z)=\overline{B_{x}^{2}}-\overline{B}_{x}^{2}$. For ease of comparisons, the normalization factor in the denominator for all cases is that value calculated for C1 and B_x is rescaled to the same units. Figure (3) shows that the toroidal magnetic field perturbations correlate well with vertical velocity perturbations, confirming that the instabilities are most likely of a magnetic buoyancy type. However, the dotted curve in this figure shows that the transport is significantly weaker for C3 in the initial stages of the instability. This is a signature of the dynamic nature of the instability. Owing to the evolving background state, R_DD does not entirely determine stability. In the equilibrium stability problem for quasi-incompressible motion without an imposed shear considered in Hughes & Weiss (1995), the growth rate of the instability becomes positive when R_DD ⩾ 27π⁴/4. In the current problem, the portion of the domain that is unstable is perpetually moving as the magnetic layer evolves. Therefore, a locally unstable perturbation may not have enough time to amplify to a dynamically significant magnitude before it finds itself outside of the moving region of instability. Only if the growth rate of the instability is sufficiently large can an unstable mode amplify fast enough that the background appears effectively steady. Therefore, to exhibit instability in this evolving system unambiguously, the growth rate of the instability must be at least of the order of the Alfvénic evolution timescale for the background,

$$\begin{eqnarray} &&p_{\mathrm{crit}}(R_{T},R_{B},\sigma,\zeta)\approx \Delta {z}\sqrt{\sigma \zeta {Q}}, \end{eqnarray} \tag{ 8 }$$

where p_crit(R_T, R_B, σ, ζ) is the maximal growth rate as a function of the Rayleigh numbers as given in Hughes & Weiss (1995), and the Chandrasekhar number Q ≫ 1. Equation (8) implies that in the presence of a stronger background field, the instability requires a faster growth rate to proceed in the same manner. Even though C3 has the exact same background $\overline{B}_{x}(z)$ configuration as that for C1, the stronger background field strength, B_z,0, alters the nature of the instability. At the location where C1 grows strongly, C3 only reaches small amplitude (albeit in a shorter time: t ≈ 44 cf. t ≈ 88). However, the instability in C3 can still grow to significant amplitudes higher up in the domain as the system evolves further, as shown by the dashed curve in Figure (3). The amplitude is further modified since the stable background density stratification, and therefore R_DD changes slightly with height.

We have demonstrated that the presence of shear-generated double-diffusive magnetic buoyancy instabilities can be controlled primarily by the single parameter ζ. Given the dynamic nature of the instability here, the second most important parameter is the Chandresekhar number, Q, which dictates the background Alfvénic timescales. C3 provides good evidence that the efficiency of the instability is governed through a relation at least qualitatively similar to Equation (8). However, the question remains as to whether any other parameters, such as the magnetic Prandtl number, σ_M = σ/ζ, are critical to this process. Numerical constraints restrict our simulations to σ_M = 0.5, which is less than unity, thus preserving the correct ordering of timescales, but larger than the estimated solar value of σ_M ≈ 6 × 10⁻² (Gough 2007). To make some progress, we extend the analysis of Hughes & Weiss (1995), noting again that this discusses the double-diffusive instabilites of an equilibrium rather than a dynamically evolving system. It can be shown that Equation (8) embodies the solution to a 12th-order polynomial in p_crit. In the limit σ, ζ ≪ 1 and R_T ≪ 0 this can be simplified to a cubic polynomial depending solely on the reduced variables r = R_DD/|R_T|, q = Δz²Q/|R_T|, and σ_M:

$$\begin{eqnarray} &&4X_{r}^{3}+X_{q}X_{r}^{2}-18X_{q}X_{r}-X_{q}(27+4X_{q})=0, \end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray} &&X_{q}=\frac{q(1+\sigma _{\mathrm{M}})^{3}}{\sigma _{\mathrm{M}}^{2}},\quad {X}_{r}=\frac{(r-q)(1+\sigma _{\mathrm{M}})}{\sigma _{\mathrm{M}}}. \end{eqnarray} \tag{ 10 }$$

Figure (4) shows the resulting dependence of the critical value of r versus q for several values of σ_M. In particular, it can be seen that the instability is more easily obtained for smaller σ_M, which makes sense in terms of the lessening relative importance of viscosity, and, encouragingly, supports the viability of this process at more realistic solar parameters.

Zoom In Zoom Out Reset image size

The transport measure $\overline{w{B_{x}}}-\overline{w}\overline{B}_{x}$ has a more complicated dependence on σ_M. In the appropriate Rayleigh number regime, $\overline{w{B_{x}}}-\overline{w}\overline{B}_{x}\sim {p_{\mathrm{crit}}}+k_{\mathrm{crit}}^{2}$, where k_crit is the horizontal wavenumber of the fastest growing mode. As σ_M is decreased, k_crit decreases and p_crit increases, and hence the dependence on σ_M is not straightforward. However, in the limit as σ_M → 0, then $\overline{w{B_{x}}}-\overline{w}\overline{B}_{x}\sim {p_{\mathrm{crit}}}$, which, from Figure (4), is nonzero even at σ_M = 0. This argument implies that, if we were able to run a computation at a significantly lower value of σ_M than in C1 and C3 (an extreme numerical challenge given our other parameters), we might expect that the buoyant transport should decrease, but not disappear altogether.

After the pessimistic results of VB1-2, it is encouraging that when double-diffusive effects are taken into account, the dynamics intuitively expected in the tachocline appear to be possible, satisfying the constraints that we have anticipated. We know roughly that ζ ≈ 10⁻⁵ and Ri ≈ 10³–10⁵ in the tachocline, and it is therefore plausible that their product is smaller than Δz/H_p ≈ 1 (Gough 2007). However, we know very little about the configuration of any magnetic field in the solar tachocline or the associated Alfvénic timescales, and hence we do not know exactly how much smaller ζ must be for the growth rate of an instability to become large enough for it to proceed.

Even if instability occurs, it still remains unclear whether the nonlinear evolution of the system can produce magnetic fields that are sufficiently strong and coherent to be able to rise through the turbulent convection zone. While it now seems possible that high Ri shear flows could produce buoyant magnetic structures in the tachocline, one might expect that such buoyant structures would have magnetic energies that are (at most) in equipartition with the shear. If the turbulent convection zone acts locally in maintaining the tachocline's shear, then it is difficult to see how magnetic structures could obtain the significantly greater strengths that are required to survive the disrupting effects of the convection zone (Tobias et al. 2001; Cline 2003; Fan et al. 2003). Of course, this question depends critically on how the energy of the shear relates to that of the convection zone, which is not well understood.