© ESO, 2011

1. Introduction

The magnetic fields in many astrophysical bodies have their origin in some kind of turbulent dynamo (Priest 1982). This means that a part of the kinetic energy of the turbulent motions is diverted to enhancing and maintaining a magnetic field. This magnetic field is generally also random, but under certain conditions a large-scale magnetic field can also emerge (Parker 1979). Here by large-scale we mean length scales larger than the energy containing scale of the fluid. In particular, large-scale fields may show up when the turbulence is helical, owing, e.g., to the simultaneous presence of rotation and stratification (Moffatt 1978).

The evolution of large-scale magnetic fields can be described by the mean-field equation that emerges upon averaging of the induction equation. In the process of averaging, new terms emerge (e.g., the α effect and turbulent diffusion) that result from correlations between small-scale velocity and magnetic fields (Krause & Rädler 1980). Here one usually considers the case where the magnetic fluctuations are caused by the fluctuating velocity acting on the mean field. However, under certain conditions it might well be the other way around. Imagine, for example, the case where initially no turbulent velocity is present, but there is instead a strong large-scale magnetic field the presence of which makes the initial state unstable. The magnetic field would then be responsible for driving velocity and magnetic fluctuations at the same time. This type of scenario was first simulated in the context of accretion discs where the magneto-rotational instability drives the turbulence (Brandenburg et al. 1995), and later in the context of the magnetic buoyancy instability with shear (Cline et al. 2003). The latter is relevant to the overshoot layer of the Sun. It had already been proposed by Moffatt (1978) that, once the dynamo-generated magnetic field in this layer reaches appreciable strengths, the magnetic buoyancy instability can set in and govern the dynamics thereafter. In this paper we consider a similar case in which an otherwise stable layer of rotating fluid is made unstable by the presence of a strong magnetic field. This is known as the magnetic buoyancy instability.

If the length scale of the initial magnetic field is larger than or of the same order as the height of the fluid layer, the magnetic Boussinesq approximation is often used. Linear stability analyses for such cases (Acheson 1979; Hughes 1985a,b) have shown that the initial state is unstable to both two-dimensional perturbations in the plane perpendicular to the initial magnetic field and three-dimensional perturbations. The two-dimensional perturbations do not bend the the initial magnetic field lines and are known as interchange modes. The three-dimensional perturbations are often called undular modes. Here we concentrate on the other limit where the length scale of initial magnetic field is significantly smaller than the height of the fluid layer. Such a magnetic field is also unstable to both the interchange and undular modes. The linear phase of this instability with stratification and rotation has been studied in detail by Schmitt (1984, 1985) and Ferriz-Mas et al. (1994) for ideal MHD. A necessary but not sufficient condition for instability is (Eq. (3.3) of Acheson 1979) $\frac{\mathit{\partial }}{\mathit{\partial z}}\log \left(\frac{\mathit{B}}{\mathit{\rho }}\right)\mathit{<}\mathrm{0}\mathit{,}$(1)which essentially means that the magnetic field modulus B decreases faster with height z than the density ρ. Fan (2001) has studied the linear stability of the same problem in absence of rotation, yet in a parameter range where only the undular modes are unstable. Gilman (1970) considered the influence of large-scale shear and the case where the thermal diffusivity is large compared to viscosity and magnetic diffusivity. Silvers et al. (2009) have studied the same but in the absence of rotation. While the focus of the Gilman (1970) study has been on the formation of flux tubes from a pre-existing toroidal magnetic layer in a stably stratified atmosphere, in Silvers et al. (2009) a magnetic layer was generated from an initially vertical magnetic field in presence of strong shear. It was further shown that, when the ratio of magnetic to thermal diffusivities is sufficiently low, magnetic buoyancy can still operate.

In this paper we are not only interested in the linear growth of this instability but also in the turbulent transport coefficients characterizing the turbulence driven by such an instability. There has been several earlier attempts to calculate the coefficient α for the magnetic buoyancy instability. Brandenburg & Schmitt (1998) performed numerical calculations in presence of rotation and determined α for the resulting turbulence. Thelen (2000) has used the linearly unstable eigenmodes to calculate α, but this calculation is not applicable to the fully-developed turbulent regime which is our main interest here. Wissink et al. (2000) calculated α_yy in the nonlinear stage of the instability, but both works did not take into account the turbulent diffusivity. More recently Davies & Hughes (2011) have used the same technique as Thelen (2000) to compute the turbulent electromotive force (EMF). They, however, claim that the complexity of their results renders a discussion of the EMF in terms of α effect and turbulent diffusivity impossible. All these works rely upon fully compressible numerical calculations and hence are closest to the present one. In the concluding section of this paper we shall compare their results with ours.

The focus of this work is twofold. Firstly, we study the nature of the magnetic buoyancy instability in its initial linear stage. In particular, we study its dependence on various parameters such as magnetic and thermal Prandtl numbers, angular velocity, strength of the initial field, etc., and compare against the linear theory and previous numerical work. Secondly, to gain some insights on whether the instability constitutes a viable dynamo process, we measure the mean-field transport coefficients, namely the tensors α and η using the quasi-kinematic test-field (QKTF) method (Schrinner et al. 2005, 2007). Note that, with one exception (Vermersch & Brandenburg 2009), the QKTF method has not been applied previously to the calculation of transport coefficients for the case of an inhomogeneous turbulence induced by the mean magnetic field. Therefore we validate the QKTF method for this problem by reconstruction of the turbulent EMF from the turbulent transport coefficients. For a review on the transport coefficients and their determination using test fields; see Brandenburg et al. (2010).

Fig. 1 The Cartesian simulation domain with respect to spherical coordinates.

2. The model

We consider a setup similar to that described in Brandenburg & Schmitt (1998). The computational domain is a cuboid with constant gravity, g_z, pointing in the negative z direction, and rotation Ω making an angle θ with the vertical. The box may be thought to be placed at a colatitude θ on the surface of a sphere with its unit vectors pointing along the local θ,   φ,   r directions, respectively, as shown in Fig. 1.

We solve the following set of MHD equations. The continuity equation is given by $\frac{\mathrm{D}\ln \mathit{\rho }}{\mathrm{D}\mathit{t}}\mathrm{=}\mathrm{-}\mathbf{\nabla }\mathrm{·}{U}\mathit{,}$(2)where D/Dt ≡ ∂/∂t + U·∇ denotes the Lagrangian derivative with respect to the local velocity of the gas U. Assuming an ideal gas, we express the pressure in terms of density, specific entropy s, and sound speed c_s, which, in turn, is a function of ρ and s. Thus the momentum equation in a frame of reference rotating with angular velocity Ω reads $\hspace{0.17em}\begin{array}{}\begin{array}{}\end{array}\frac{\mathrm{D}{U}}{\mathrm{D}\mathit{t}}\mathrm{=}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathbf{\nabla }(\frac{\mathit{s}}{{\mathit{c}}_{\mathrm{p}}}\mathrm{+}\ln \mathit{\rho })\mathrm{-}\mathrm{2}\mathbf{\Omega }\mathrm{\times }{U}\mathrm{+}{\mathit{g}}_{\mathit{z}}{ẑ}\mathrm{+}\frac{{J}\mathrm{\times }{B}}{\mathit{\rho }}\begin{array}{}\end{array}\\ \begin{array}{}\end{array}\mathrm{+}\mathit{\nu }\left({\mathrm{\nabla }}^{\mathrm{2}}{U}\mathrm{+}\frac{\mathrm{1}}{\mathrm{3}}\mathbf{\nabla }\mathbf{\nabla }\mathrm{·}{U}\mathrm{+}\mathrm{2}\mathbf{S}\mathrm{·}\mathbf{\nabla }\ln \mathit{\rho }\right)\mathit{,}\end{array}$(3)where J is the current density, B is the magnetic field, ν is the constant kinematic viscosity, and S is the traceless rate-of-strain tensor. The sound speed is related to temperature by ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\mathrm{(}{\mathit{c}}_{\mathrm{p}}\mathrm{-}{\mathit{c}}_{\mathrm{v}}\mathrm{)}\mathit{\gamma T}$ with c_p and c_v the specific heat at constant pressure and constant volume, respectively, and γ = c_p/c_v is here fixed to 5/3. The induction equation is solved in terms of the magnetic vector potential A, such that ∇ × A = B, hence $\frac{\mathit{\partial }{A}}{\mathit{\partial t}}\mathrm{=}{U}\mathrm{\times }{B}\mathrm{+}\mathit{\eta }{\mathrm{\nabla }}^{\mathrm{2}}{A}\mathit{,}$(4)where η denotes constant molecular magnetic diffusivity.

Finally, we have for the entropy equation with temperature T and constant radiative (thermal) conductivity K$\mathit{\rho T}\frac{\mathrm{D}\mathit{s}}{\mathrm{D}\mathit{t}}\mathrm{=}\mathbf{\nabla }\mathrm{·}\mathrm{\left(}\mathit{K}\mathbf{\nabla }\mathit{T}\mathrm{\right)}\mathrm{+}\mathit{\eta }{\mathit{\mu }}_{\mathrm{0}}{J}{\mathrm{2}}^{}\mathrm{+}\mathrm{2}\mathit{\rho \nu }\mathbf{S}{\mathrm{2}}^{}\mathit{,}$(5)where the temperature is related to the specific entropy by $\mathit{s}\mathrm{=}{\mathit{s}}_{\mathrm{0}}\mathrm{+}{\mathit{c}}_{\mathit{v}}\ln \frac{\mathit{T}\mathit{/}{\mathit{T}}_{\mathrm{0}}}{\mathrm{(}\mathit{\rho }\mathit{/}{\mathit{\rho }}_{\mathrm{0}}{\mathrm{)}}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}}\mathrm{·}$(6)We use the fully compressible Pencil Code¹ for all our calculations.

For all quantities, periodic boundary conditions in the x and y directions are adopted. In the z direction we use the no-slip boundary condition for the velocity and the vertical field condition B_x = B_y = 0, that is, J_z = 0, for the magnetic field as a proxy for vacuum boundaries. We keep the temperature at the top and the (radiative) heat flux at the bottom fixed. Their values were chosen to conform with the initial temperature profile of the (not magnetically modified) polytrope described below.

2.1. Initial state

The base state is a polytrope that is, p = Cρ^Γ, with index m = 1/(Γ − 1) = 3. The initial z profiles of density, pressure, temperature and entropy are given by $\hspace{0.17em}\begin{array}{}\begin{array}{}\end{array}{\mathit{\rho }}_{\mathrm{i}}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}{\mathrm{\Phi }}^{\mathrm{3}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{,} {\mathit{p}}_{\mathrm{i}}\mathrm{=}{\mathit{p}}_{\mathrm{0}}{\mathrm{\Phi }}^{\mathrm{4}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{,} {\mathit{T}}_{\mathrm{i}}\mathrm{=}{\mathit{T}}_{\mathrm{0}}\mathrm{\Phi }\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{,}\begin{array}{}\end{array}\\ \begin{array}{}\end{array}{\mathit{s}}_{\mathrm{i}}\mathrm{=}{\mathit{s}}_{\mathrm{0}}\mathrm{-}{\mathit{c}}_{\mathit{v}}\ln \mathrm{\Phi }\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{,}\end{array}$(7)where Φ is a non-dimensional gravitational potential given by $\mathrm{\Phi }\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{4}}\frac{{\mathit{g}}_{\mathit{z}}}{\mathrm{(}{\mathit{c}}_{\mathit{p}}\mathrm{-}{\mathit{c}}_{\mathit{v}}\mathrm{)}{\mathit{T}}_{\mathrm{0}}}\mathrm{(}\mathit{z}\mathrm{-}{\mathit{z}}_{\mathrm{0}}\mathrm{)}\mathit{,}$with the reference point z₀ chosen to be at the bottom of the domain and the values at this point given by ρ₀, ${\mathit{p}}_{\mathrm{0}}\mathrm{=}{\mathit{c}}_{\mathrm{s}\mathrm{0}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}\mathit{/}\mathit{\gamma }$, ${\mathit{T}}_{\mathrm{0}}\mathrm{=}{\mathit{c}}_{\mathrm{s}\mathrm{0}}^{\mathrm{2}}\mathit{/}\mathrm{(}{\mathit{c}}_{\mathit{p}}\mathrm{-}{\mathit{c}}_{\mathit{v}}\mathrm{)}\mathit{\gamma }$ and s₀. Here, c_s0 is the reference sound speed to which we also refer to when calculating Mach numbers.

As the adiabatic index here is m_ad = 1/(γ − 1) = 3/2, the subadiabaticity in the domain is very large, namely ∂lnT/∂lnP − (∂lnT/∂lnP)_ad = −0.15. Thus, the initial stratification is highly stable to convection in the absence of any magnetic field, guaranteeing that turbulence is generated solely by the buoyancy instability.

The initial magnetic field is a horizontal layer of thickness 2H_B, where B_y has the profile ${\mathit{B}}_{\mathit{y}\mathrm{0}}\mathrm{=}{\mathit{B}}_{\mathrm{0}}{\mathit{H}}_{\mathit{B}}\frac{\mathit{\partial }}{\mathit{\partial }\mathit{z}}\tanh \left(\frac{\mathit{z}\mathrm{-}{\mathit{z}}_{\mathit{B}}}{{\mathit{H}}_{\mathit{B}}}\right)\mathit{,}$(8)and the reference Alfvén speed is defined by ${\mathit{v}}_{\mathrm{A}\mathrm{0}}\mathrm{=}{\mathit{B}}_{\mathrm{0}}\mathit{/}\sqrt{{\mathit{\rho }}_{\mathrm{0}}{\mathit{\mu }}_{\mathrm{0}}}$ with μ₀ being the vacuum permeability. If not indicated otherwise, the initial magnetic field strength is fixed to v_A0/c_s0 = 0.5. In order to satisfy the condition (1)initially, we have to ensure 2H_B < H_ρ(z_B), where H_ρ(z) = |∇lnρ(z)|^-1 is the local density scale height. When choosing z_B − z₀ = 0.3L_z this is satisfied for 2H_B < 0.1L_z + 4T₀(c_p − c_v)/3|g_z| which is surely true for the choice H_B = 0.05L_z.

Upon addition of a magnetic field, we modify the base state such that the density profile remains unchanged. In order to obey magnetostatic equilibrium, pressure and temperature are adjusted in the following way: ${\mathit{p}}_{\mathrm{i}}\mathrm{\Rightarrow }{\mathit{p}}_{\mathrm{i}}\mathrm{-}\frac{{\mathit{B}}_{\mathit{y}\mathrm{0}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}}\mathit{,} {\mathit{T}}_{\mathrm{i}}\mathrm{\Rightarrow }{\mathit{T}}_{\mathrm{i}}\mathrm{-}\frac{{\mathit{B}}_{\mathit{y}\mathrm{0}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}}\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{i}}\mathrm{(}{\mathit{c}}_{\mathit{p}}\mathrm{-}{\mathit{c}}_{\mathit{v}}\mathrm{)}}\mathrm{·}$(9)The entropy is then re-calculated from Eq. (6). This leads to a local minimum in the temperature profile which would smooth out over a diffusion time scale ${\mathit{H}}_{\mathit{B}}^{\mathrm{2}}\mathit{/}{\mathit{\chi }}_{\mathrm{0}}$, being in any case much larger than the e-folding time of the instability. Another common choice is to keep the temperature profile unchanged, but this implies a decrease of the density in the magnetic layer, making it immediately buoyant in presence of small perturbations.

The initial velocity components U_x and U_y are specified such that they contain about 20 localized eddies in the plane z = z_B with Mach numbers of about 10^-5. Also the initial vertical velocity, U_z is Gaussian random noise with the same Mach number and the rms value of the initial kinetic helicity, scaled with the product of initial rms velocity and vorticity, $\mathrm{(}\sqrt{\mathrm{⟨}\mathrm{(}{W}\mathrm{·}{U}{\mathrm{)}}^{\mathrm{2}}\mathrm{⟩}}\mathit{/}{\mathit{U}}_{\mathrm{rms}}{\mathit{W}}_{\mathrm{rms}}\mathrm{)}\mathrm{(}\mathit{t}\mathrm{=}\mathrm{0}\mathrm{)}$, is 4 × 10^-6. Here and in the following angular brackets mean volume averaging.

2.2. Control parameters, non-dimensional quantities, and computational grid

The problem posed by (2)through (5)is governed by five independent dimensionless parameters; (i) the Prandtl number Pr = ν/χ₀, with the temperature conductivity χ₀ = K/ρ₀c_p; (ii) the magnetic Prandtl number Pr_M = ν/η; (iii) the “magnetic Taylor number” ${\mathrm{Ta}}_{\mathrm{M}}\mathrm{=}\mathrm{2}{\mathrm{\Omega }}^{\mathrm{2}}{\mathit{L}}_{\mathit{y}}^{\mathrm{4}}\mathit{/}{\mathit{\eta }}^{\mathrm{2}}$; (iv) the rotational inclination (colatitude), θ, and (v) the normalized gravitational acceleration ${\mathit{g}}_{\mathit{z}}{\mathit{L}}_{\mathit{y}}^{\mathrm{3}}\mathit{/}{\mathit{\eta }}^{\mathrm{2}}$. In addition there are two independent parameters of the initial equilibrium (vi) the normalized pressure scale height at the bottom, ${\mathit{H}}_{\mathit{P}\mathrm{0}}\mathit{/}{\mathit{L}}_{\mathit{z}}\mathrm{=}{\mathit{c}}_{\mathrm{s}\mathrm{0}}^{\mathrm{2}}\mathit{/}\mathit{\gamma }{\mathit{g}}_{\mathit{z}}{\mathit{L}}_{\mathit{z}}$ and (vii) the initial Lundquist number, Lu₀ = v_A0H_B/η, based upon the thickness of the magnetic layer. In addition to this, we also have the non-dimensional sound speed, c_s0L_y/η. In this paper we shall keep the normalized pressure scale height and the sound speed fixed, while varying both Prandtl numbers, Ta_M, θ and Lu₀. The definitions as well as the values or ranges of the control parameters are summarized in Table 1. We have also included in the same table two dependent parameters namely the Roberts number Rb = Pr_M/Pr = χ₀/η and the modified initial plasma beta $\hspace{0.17em}{\begin{array}{}{}^{\mathrm{􏽥}}\\ \mathit{\beta }\end{array}}_{\mathrm{0}}$ in the midplane of the magnetic layer. It is defined as the ratio of the total pressure p_tot = p + p_M to the magnetic pressure ${\mathit{p}}_{\mathrm{M}}\mathrm{=}{\mathit{B}}_{\mathit{y}\mathrm{0}}^{\mathrm{2}}\mathit{/}\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}$, because $\hspace{0.17em}{\begin{array}{}{}^{\mathrm{􏽥}}\\ \mathit{\beta }\end{array}}_{\mathrm{0}}$ then adopts a simple 1/B² dependence on the magnetic field; cf. Eq. (9).

The computational domain is defined by |x| ≤ L_x/2, |y| ≤ L_y/2, −L_z/4 ≤ z ≤ 3L_z/4, L_x = L_z = L_y/3, thus its aspect ratio is 1:3:1. The results will be presented in non-dimensional form, velocity in units of the reference Alfvén speed, v_A0, time in units of the corresponding Alfvén travel time in the y direction, t_A0 = L_y/v_A0, and magnetic field in units of B₀ or the rms value ${\mathrm{(}}^{\mathrm{\int }}{}_{\mathit{z}}\mathit{B}_{\mathit{y}\mathrm{0}}^{\mathrm{2}}\mathrm{d}\mathit{z}\mathit{/}{\mathit{L}}_{\mathit{z}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$.

It is instructive to look upon the relevant definitions of the fluid Reynolds number, Re, and the magnetic Reynolds number, R_m, for this problem where the turbulence is driven solely by the instability of the magnetic layer. From first principles, Re characterizes the ratio of the advective term ⟨ (U·∇U)² ⟩ ^1/2 and the approximate viscous term ⟨ (ν∇²U)² ⟩ ^1/2 in the Navier-Stokes equation, while R_m characterizes the ratio of ⟨ (∇ × (U × B))² ⟩ ^1/2 and ⟨ (η∇²B)² ⟩ ^1/2 in the induction equation with the angular brackets representing volume averaging. Let us denote these ab initio definitions as “term based” and refer to them by Re ^∗ and ${\mathit{R}}_{\mathrm{m}}^{\mathrm{\ast }}$. Note that with the term-based definitions R_m/Re may well deviate from Pr_M. Alternatively, we can define a length scale L_U = U_rms/2πW_rms from the rms values of velocity and vorticity and define the more conventional length-based Reynolds numbers Re = U_rmsL_U/ν and R_m = U_rmsL_U/η.

The calculations were carried out on equidistant grids with resolutions of either 64³ or 128³. For numerical testing we have also performed a few runs with 256³ or 128² × 256 resolutions.

2.3. The test-field method

Here, for the sake of completeness, we give a concise description of the test-field method we use to calculate the turbulent transport coefficients. For a more detailed discussion we refer the reader to Brandenburg et al. (2008c) and Rheinhardt & Brandenburg (2010).

Let us define mean magnetic and velocity fields, $\overline{{B}}$ and $\overline{{U}}$, where overbars denote Reynolds averaging. Fluctuations are defined correspondingly as ${b}\mathrm{=}{B}\mathrm{-}\overline{{B}}$ and ${u}\mathrm{=}{U}\mathrm{-}\overline{{U}}$. Following these conventions, the induction equation may be correspondingly averaged, resulting in $\frac{\mathit{\partial }\overline{{B}}}{\mathit{\partial }\mathit{t}}\mathrm{=}\mathbf{\nabla }\mathrm{\times }\mathrm{(}\overline{{U}}\mathrm{\times }\overline{{B}}{\mathrm{)}}^{}\mathrm{+}\mathbf{\nabla }\mathrm{\times }\overline{\mathbf{ℰ}}\mathrm{+}\mathit{\eta }{\mathrm{\nabla }}^{\mathrm{2}}\overline{{B}}\mathit{,}$(10)where $\overline{\mathbf{ℰ}}\mathrm{\equiv }\overline{{u}\mathrm{\times }{b}}$ is the mean EMF resulting from the fluctuating fields u and b. The essence of mean-field magneto-hydrodynamics is to provide an expression for $\overline{\mathbf{ℰ}}$ as a functional of the large scale magnetic field and its derivatives. A simple ansatz reads $\overline{\mathbf{ℰ}}\mathrm{=}{\alpha }\overline{{B}}\mathrm{-}{\eta }\mathbf{\nabla }\overline{{B}}\mathit{,}$(11)where α and η are called transport coefficients. Note that a much more general representation of $\overline{\mathbf{ℰ}}$ is given by the convolution integral $\overline{\mathbf{ℰ}}\mathrm{\left(}{x}\mathit{,t}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{{\mathit{t}}_{\mathrm{0}}}^{\mathit{t}}\mathrm{\int }{G}\mathrm{(}{x}\mathit{,}{x}\mathrm{-}{x}\mathrm{\prime }\mathit{,t,t}\mathrm{-}{\mathit{t}}^{\mathrm{\prime }}\mathrm{)}\hspace{0.17em}\overline{{B}}\mathrm{\left(}{x}\mathrm{\prime }\mathit{,}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathrm{3}}^{}{\mathit{x}}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}$(12)with an appropriate tensorial kernel G or Green’s function. When considering Fourier transforms with respect to x′ (and now dropping the time dependence) we have $\overline{\mathbf{ℰ}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}\mathrm{\int }\hspace{0.17em}\hat{{G}}\mathrm{\left(}{x}\mathit{,}{k}\mathrm{\right)}\hspace{0.17em}\hspace{0.17em}\widehat{\overline{{B}}}\mathrm{\left(}{k}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}{x}}\hspace{0.17em}\mathrm{d}{\mathrm{3}}^{}\mathit{k}$(13)with $\hspace{0.17em}\hat{{G}}\mathrm{\left(}{x}\mathit{,}{k}\mathrm{\right)}$ being the Fourier transform of G(x,x′) with respect to x′. As shown in Brandenburg et al. (2008c) a decomposition in the form $\hspace{0.17em}\mathit{Ĝ}\mathit{ij}\mathrm{\left(}{x}\mathit{,}{k}\mathrm{\right)}\mathrm{=}\mathit{\alpha ̂}\mathit{ij}\mathrm{\left(}{x}\mathit{,}{k}\mathrm{\right)}\mathrm{-}\mathrm{i}\mathit{\eta ̂}\mathit{ijk}\mathrm{\left(}{x}\mathit{,}{k}\mathrm{\right)}{\mathit{k}}_{\mathit{k}}$(14)with real $\hat{{\alpha }}$ and $\hat{{\eta }}$ is always achievable. Further, the traditional (local) coefficients are given by $\hat{{\alpha }}\mathrm{(}{k}\mathrm{=}\mathbf{0}\mathrm{)}$ and ${\lim }_{{k}\mathrm{\to }\mathbf{0}}\hat{{\eta }}\mathrm{\left(}{k}\mathrm{\right)}$.

The aim of the test-field method is to provide an expression for G as a functional of properties of the mean and fluctuating fluid velocity, $\overline{{U}}$ and u by utilizing direct numerical simulations. Subtracting the averaged induction Eq. (10)from the original one (4), written in terms of B rather than A, we obtain the following equation for the fluctuating magnetic field b. $\frac{\mathit{\partial }\hspace{0.17em}{b}{\mathit{pq}}^{}}{\mathit{\partial }\mathit{t}}\mathrm{=}\mathbf{\nabla }\mathrm{\times }\mathrm{(}\overline{{U}}\mathrm{\times }{b}{\mathit{pq}}^{}\mathrm{+}{u}\mathrm{\times }\overline{{B}}{\mathit{pq}}^{}\mathrm{+}{e}{\mathit{pq}}^{\mathrm{)}}\mathrm{+}\mathit{\eta }{\mathrm{\nabla }}^{\mathrm{2}}{b}{\mathit{pq}}^{}\mathit{,}$(15)with ${e}{}^{\mathit{pq}}\mathrm{=}{u}\mathrm{\times }{b}{\mathit{pq}}^{}\mathrm{-}\overline{\mathbf{ℰ}}{\mathit{pq}}^{}$. The superscripts pq now indicate that this equation is solved for several suitably chosen test fields $\overline{{B}}{\mathit{pq}}^{}$ with p,q = 1,2. This is sufficient as we define the mean as horizontal (xy) average. Thus ${\overline{\mathit{B}}}_{\mathit{z}}\mathrm{=}\mathrm{0}$, and only the components α_ij and η_ij3, i,j = 1,2, are of interest. In the following we will hence refer to the rank-2 tensor η_ij = −η_ik3ϵ_jk3 only which allows to rewrite the second term in (11)into $\mathrm{-}{\eta }\hspace{0.17em}\mathrm{curl}\hspace{0.17em}\overline{{B}}$. Equation (15) is the central relation invoked by the test-field method. The test-field suite of the Pencil Code has the provision for using either harmonic test fields, i.e., $\begin{array}{ccc}\hspace{0.17em}\begin{array}{}\begin{array}{}\end{array}{\overline{\mathit{B}}}^{\mathrm{11}}\mathrm{=}\mathrm{\left(}\cos \mathit{k,}\mathrm{0}\mathit{,}\mathrm{0}\mathrm{\right)}\mathit{,} {\overline{\mathit{B}}}^{\mathrm{12}}\mathrm{=}\mathrm{\left(}\mathrm{0}\mathit{,}\cos \mathit{k,}\mathrm{0}\mathrm{\right)}\mathit{,}\begin{array}{}\end{array}\\ \begin{array}{}\end{array}{\overline{\mathit{B}}}^{\mathrm{21}}\mathrm{=}\mathrm{\left(}\sin \mathit{k,}\mathrm{0}\mathit{,}\mathrm{0}\mathrm{\right)}\mathit{,} {\overline{\mathit{B}}}^{\mathrm{22}}\mathrm{=}\mathrm{\left(}\mathrm{0}\mathit{,}\sin \mathit{k,}\mathrm{0}\mathrm{\right)}\mathit{,}\end{array}& & \end{array}$(16)or linear test fields i.e., $\begin{array}{ccc}\hspace{0.17em}\begin{array}{}\begin{array}{}\end{array}{\overline{\mathit{B}}}^{\mathrm{11}}\mathrm{=}\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathrm{\right)}\mathit{,} {\overline{\mathit{B}}}^{\mathrm{12}}\mathrm{=}\mathrm{\left(}\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{0}\mathrm{\right)}\mathit{,}\begin{array}{}\end{array}\\ \begin{array}{}\end{array}{\overline{\mathit{B}}}^{\mathrm{21}}\mathrm{=}\mathrm{\left(}\mathit{z,}\mathrm{0}\mathit{,}\mathrm{0}\mathrm{\right)}\mathit{,} {\overline{\mathit{B}}}^{\mathrm{22}}\mathrm{=}\mathrm{\left(}\mathrm{0}\mathit{,z,}\mathrm{0}\mathrm{\right)}\mathit{.}\end{array}& & \end{array}$(17)Harmonic test fields with varying wavenumber k allow us to calculate for each k the 2 × 2 tensors $\hat{{\alpha }}$ and $\hat{{\eta }}$ by solving an algebraic system of 8 equations given by $\overline{\mathrm{ℰ}}{\mathit{i}}_{\mathit{pq}}^{}\mathrm{=}{\mathit{\alpha ̂}}_{\mathit{ij}}\overline{\mathit{B}}\begin{array}{}\mathit{pq}\\ \mathit{j}\end{array}\mathrm{-}{\mathit{\eta ̂}}_{\mathit{ij}}\overline{\mathit{J}}\begin{array}{}\mathit{pq}\\ \mathit{j}\end{array}\mathrm{,with}\mathit{p,q,i,j}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{2.}$(18)As shown above, expressing the mean EMF by $\hat{{\alpha }}$ and $\hat{{\eta }}$ avoids the shortcomings of a truncation like in and provides the desired non-local relationship between $\overline{\mathbf{ℰ}}$, $\overline{{B}}$ and $\overline{{J}}$.

In the present context, the buoyancy-driven turbulence and hence the transport coefficients necessarily show an intrinsic inhomogeneity both due to stratification and the mean magnetic field itself. No specific complication is induced by it as $\hat{{\alpha }}$ and $\hat{{\eta }}$ emerge straightforwardly from the testfield method as functions of z and k.

In the nonlinear situation, the Green’s function approach remains valid if $\overline{\mathbf{ℰ}}$ is considered as a functional of U and $\overline{{B}}$ which is then linear and homogeneous in the latter. However, we have to label G by the $\overline{{B}}$ actually acting upon U, that is, ${G}\mathrm{\left(}{x}\mathit{,}{x}{}^{\mathrm{\prime }}\mathrm{;}\overline{{B}}\mathrm{\right)}$, and can thus only make statements about the transport tensors for just the particular $\overline{{B}}$ at hand. Hence, the tensors have to be labelled likewise: $\hat{{\alpha }}\mathrm{\left(}\mathit{z,k}\mathrm{;}\overline{{B}}\mathrm{\right)}$, $\hat{{\eta }}\mathrm{\left(}\mathit{z,k}\mathrm{;}\overline{{B}}\mathrm{\right)}$.

3. Results

3.1. Nature of the instability

To start with we have performed a number of runs with different values of Pr and Pr_M, but all other dimensionless parameters held fixed, see Table 2. In particular, we have used a value of Ta_M = 3.24 × 10¹⁰ for the magnetic Taylor number and Lu₀ = 500 for the initial Lundquist number. Table 2 shows the Reynolds numbers according to the two alternative definitions provided in Sect. 2.2. Note that with the exception of the run B128f, Re from the length-based and the term-based definitions are in agreement. Also the ratio ${\mathit{R}}_{\mathrm{m}}^{\mathrm{\ast }}\mathit{/}\mathrm{Re}{}^{\mathrm{\ast }}$ from the term-based definitions approaches Pr_M reasonably.

We first show the temporal evolution of the magnetic field for a few representative cases in Fig. 2. In all of them, we can clearly distinguish a first stage of exponential growth from a subsequent saturation phase. The x and z components of the magnetic field are generated at the expense of its y component. Although there exists a persistent energy source in the form of a constant heat flux into the domain, the final saturated stage always undergoes a slow decay. This decay is most clearly visible in B_y. Thus the instability is not able to maintain a dynamo on its own.

Fig. 2 Time evolution of the runs in Table 2. Upper panel: mean squared values of velocity and generated magnetic field components Bx, Bz of run B128c (Pr = PrM = 1) scaled by $\mathrm{⟨}{\mathit{B}}_{\mathit{y}\mathrm{0}}^{\mathrm{2}}\mathit{/}{\mathit{\mu }}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{0}}\mathrm{⟩}$ and $\mathrm{⟨}{\mathit{B}}_{\mathit{y}\mathrm{0}}^{\mathrm{2}}\mathrm{⟩}$ respectively. Note the clear exponential growth until t ≈ 1.4tA0. Fast oscillations in ⟨ U2 ⟩ until t ≈ tA0 indicate g-modes originating from the initial velocity perturbation. Middle panel: mean squared values of generated magnetic field components for different runs. For legend see lower panel. Lower panel: rms value of By. Prandtl numbers indicated as (Pr, PrM).

We suppose that the magnetic layer formed by B_x, though having a vertical scale suited to maintain the instability, is eventually not strong enough to take over the role of the initial magnetic layer. Let us first discuss the initial linear stage of the instability.

Fig. 3 Top: velocity components Uy (in color), Ux and Uz (vectors) in the plane y = 0. Bottom: Ux (in color), Uy and Uz (vectors) in the plane x = 0. Both during the linear evolution phase of the run B128a. The initial condition for U was set to be harmonic in x and y with kx = 8 × 2π/Lx and ky = 2π/Ly.

3.2. Linear stage

At first we verify that the instability is indeed driven by magnetic buoyancy. As the coefficients in Eqs. (2)–(5)are constant, the initial state (7)depends only on z, and the boundary conditions in the x and y directions are periodic, all eigensolutions ψ = (ρ,u,b,s) of the linearized problem must have the form ${\psi }\mathrm{\left(}{x}\mathit{,t}\mathrm{\right)}\mathrm{=}\hspace{0.17em}\hat{{\psi }}\mathrm{\left(}\mathit{z}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{2}\mathit{\pi }\mathrm{i}\mathrm{(}\mathit{mx}\mathit{/}{\mathit{L}}_{\mathit{x}}\mathrm{+}\mathit{ny}\mathit{/}{\mathit{L}}_{\mathit{y}}\mathrm{)}\mathrm{-}\mathrm{i}\mathit{\omega t}}\mathit{,}$(19)where m and n are integers and ω = ω_R + iω_I. For n ≠ 0 the undular modes are obtained. Corresponding dispersion relations ω(m,n) have been found by Fan (2001), using variational minimization of energy, and by Schmitt (1985), using eigenmode decomposition of linearized MHD equations. The former considers the non-rotating case, and the latter the rotating one. Both works use the anelastic approximation and their results are strictly applicable only for ideal MHD. Fan (2001) reports only non-oscillatory modes. Schmitt (1985), who further uses the magnetostrophic approximation, finds that the growing modes also have an oscillatory part, with the ratio ω_R/ω_I decreasing with latitude. Note however that the results of Schmitt (1985) do not directly apply to the present case because our boundary conditions are different².

For the runs in Table 2 we find that the early exponential growth phase is dominated by a mixture of several eigenmodes having almost the same growth rates, but different x wavenumbers, 5 ≤ m ≤ 8, whereas for all of them n = 1. This is further substantiated by the occasional occurrence of a beating in the x dependence with a modulation period given by m = 1 which indicates the presence of two almost equally strong eigenmodes with m and m ± 1. By setting the initial conditions to be harmonic in x and y it is also possible to select single eigenmodes cleanly. Figure 3 shows the velocity pattern of the fastest growing mode with m = 8, n = 1. According to the terminology of Hughes (1985a) we may qualify our eigenmodes as undular as they change periodically in the direction of the initial magnetic field. Our results are consistent with the findings of Thelen (2000) where for moderate rotation the fastest growing mode had always the smallest possible (non-vanishing) wavenumber in the direction of the field while the wavenumber perpendicular to the field was high.

The growth rates ω_I = ℑ(ω) presented in Table 2, have to be considered as average values for a (not precisely known) number of eigenmodes. While the ω_I could be easily identified from the averaged quantities shown in Fig. 2, it is difficult to access the oscillation frequencies ℜ(ω). This is because they are small compared to the growth rates and saturation sets in too early to allow for the observation of a complete oscillation period. Nevertheless, some indications for temporal variations in the eigenmode geometries have been found. Generally, we observe an increase of the growth rate with increasing magnetic Prandtl number, but a decrease with increasing Prandtl number. We find that the growth rate increases with the Roberts number as shown in Fig. 4. This means that increasing the efficiency of heat conduction in comparison to magnetic diffusion destabilizes the sub-adiabatic stratification in the system in agreement with the destabilizing effect of thermal diffusion studied by Acheson (1979).

Fig. 4 Dependence of growth rate ωI on the inverse Roberts number derived from the runs in Table 2. Solid line: best linear fit. Size of circles codes for the value of Rm (length based, see Table 2).

3.3. Dependence on initial magnetic field and rotation

Another piece of evidence for the magnetic character of the instability is its dependence on the initial magnetic field strength. From Fig. 5 we see a clear increase of the growth rate and saturation level with decreasing $\hspace{0.17em}{\begin{array}{}{}^{\mathrm{􏽥}}\\ \mathit{\beta }\end{array}}_{\mathrm{0}}$, that is, increasing Lu₀, while keeping the rotation rate fixed at Ta_M = 3.24 × 10¹⁰ (see also the last three runs in Table 3 below). Schmitt (2003) predicted in the magnetostrophic approximation a growth rate ∝ ${\mathit{v}}_{\mathrm{A}\mathrm{0}}^{\mathrm{2}}\mathit{/}\mathrm{\Omega }$ inversely proportional to $\mathit{\zeta }\mathrm{=}\hspace{0.17em}{\begin{array}{}{}^{\mathrm{􏽥}}\\ \mathit{\beta }\end{array}}_{\mathrm{0}}\sqrt{\mathrm{T}{\mathrm{a}}_{\mathrm{M}}}$.

Next we keep $\hspace{0.17em}{\begin{array}{}{}^{\mathrm{􏽥}}\\ \mathit{\beta }\end{array}}_{\mathrm{0}}$ constant at 1.51 and decrease ζ gradually from 4.09 × 10⁵ to 0. Inspecting Fig. 5, we find that growth rate and saturation level of $\mathrm{⟨}{\mathit{B}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{B}}_{\mathit{z}}^{\mathrm{2}}\mathrm{⟩}$ increase monotonically and reach their maxima at ζ = 0 (Ω = 0) while the saturation time is decreasing. The impeding effect of rotation onto the instability at large Ω is plausible in view of the Taylor-Proudman theorem because the unstable eigenmodes do show pronounced z gradients in U, see Fig. 3.

Fig. 5 Dependence of the instability on initial magnetic field strength, expressed by $\hspace{0.17em}{\begin{array}{}{}^{\mathrm{􏽥}}\\ \mathit{\beta }\end{array}}_{\mathrm{0}}$, and rotation, expressed by $\mathit{\zeta }\mathrm{=}\hspace{0.17em}{\begin{array}{}{}^{\mathrm{􏽥}}\\ \mathit{\beta }\end{array}}_{\mathrm{0}}{{\mathrm{Ta}}_{\mathrm{M}}}^{\mathrm{1}\mathit{/}\mathrm{2}}$. Upper panel: $\mathrm{⟨}{\mathit{B}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{B}}_{\mathit{z}}^{\mathrm{2}}\mathrm{⟩}$. Lower panel: $\mathrm{⟨}{\mathit{B}}_{\mathit{y}}^{\mathrm{2}}\mathrm{⟩}$. Legend shows $\mathrm{(}\hspace{0.17em}{\begin{array}{}{}^{\mathrm{􏽥}}\\ \mathit{\beta }\end{array}}_{\mathrm{0}}\mathit{,}{\mathrm{10}}^{-5}\mathit{\zeta }\mathrm{)}$. Pr = PrM = 1, colatitude θ = 30°, resolution 643 throughout. Note that the normalization time tA0 is not the same for all curves, but proportional to $\hspace{0.17em}\begin{array}{}{}^{\mathrm{􏽥}}\\ \mathit{\beta }\end{array}{\mathrm{0}}_{\mathrm{1}\mathit{/}\mathrm{2}}^{}$.

3.4. Saturated stage

At later time the instability reaches saturation, characterized by turbulent magnetic, velocity, density and temperature fields, that decay slowly thereafter. However, in most of the analysis below, this decay will be ignored and the turbulence considered as approximately statistically stationary. The turbulence is necessarily both inhomogeneous and anisotropic and we shall further show that it is also helical. Under such conditions we expect the emergence of a mean EMF. Indeed magnetic fields perpendicular to the initial magnetic layer are produced having non-vanishing horizontal averages.

In order to give a better idea of the 3D geometry of the magnetic field we provide in Fig. 6 a volume rendering of B_y at a time after t^sat for the runs B128e and B128f (see Table 2) which differ only in their magnetic Prandtl numbers. Notice how the magnetic layer breaks into flux tubes – similar to what is seen in Fig. 3 of Fan (2001) and also in Matthews et al. (1995). More relevantly, it is also similiar to the findings of Wissink et al. (2000) who include rotation. The difference between the two runs is most striking in the nature of corrugation in the surface shown. We attribute it to larger twist in the rising tubular structures in the run B128e compared to B128f which becomes clearly visible in the field line pictures also depicted in Fig. 6 (right).

Figure 7 demonstrates the breakup of the magnetic layer into tubular structures of concentrated magnetic field which are also regions of low density, hence rising. Notice also the high density regions just above and below these tubular structures. They show a significantly lower temperature than their surroundings (Fig. 7, bottom).

Studies of rising flux tubes in simulations of the solar convection zone have shown that their ability to rise depends crucially on how twisted they are (see, e.g., Emonet & Moreno-Insertis 1998). Hence it behooves us to measure the twist of the large-scale magnetic field in our simulations. For a quantitative measurement we utilize the dimensionless parameter ϵ_J = ⟨ J·B ⟩ /J_rmsB_rms, the relative current helicity, essentially measuring the overall degree of alignment between B and J. Here, angular brackets denote volume averages. The corresponding horizontally averaged quantity is ${\mathit{\epsilon }}_{\mathit{J}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\overline{{J}\mathrm{·}{B}}\mathit{/}{\mathit{J}}_{\mathrm{rms}}{\mathit{B}}_{\mathrm{rms}}$. Figure 8 shows the normalized current helicity density J·B/J_rmsB_rms (filled contours) as well as B_y in the plane y = 0 for run B128e. Notice that the contours are bend leftward because of the Coriolis force with Ω·ẑ > 0. The contour plots of B_y in this figure also show the formation of rising tubular structures from the magnetic layer.

In Fig. 9 we show the dependence of ϵ_J on Rb and the profile of ε_J(z) for some selected runs. The latter quantity can be as strong as 30% near to the initial location of the magnetic sheet, but the total helicity reaches only values of a few per cent. The clear dependence of ϵ_J on Rb is in contrast to the only weak dependences on Pr_M and Pr individually. This is an important result from this section. Our conjecture is that, at large Rb, the magnetic buoyancy instability may play an important role in the formation of twisted flux tubes in the Sun, where Rb ≫ 1 is expected.

Fig. 6 Left: volume rendering of the By = 0.1B0 isosurface. Right: field lines, colored according to the value of vorticity Wy for runs B128f (top) and B128e (bottom) at t = 2tA0 (saturated stage). The contours of By are shown on side walls.

To demonstrate the emergence of a mean magnetic field we present in Fig. 10 time-height plots of ${\overline{\mathit{B}}}_{\mathit{y}}$ and ${\overline{\mathit{B}}}_{\mathit{x}}$ for the run B128g (note that ${\overline{\mathit{B}}}_{\mathit{z}}\mathrm{=}\mathrm{0}$). There, t ≈ 1.6   t_A0 marks the end of the exponential growth phase after which a strong growth of ${\overline{\mathit{B}}}_{\mathit{x}}$, obviously at the expense of ${\overline{\mathit{B}}}_{\mathit{y}}$, sets in. ${\overline{\mathit{B}}}_{\mathit{x}}$ reaches its maximum around t ≈ 3t_A0 and is then subject to the overall decay. Note the strong vertical concentration of ${\overline{\mathit{B}}}_{\mathit{x}}$, approximately antisymmetric about the midplane of the magnetic sheet.

3.5. Calculation of turbulent transport coefficients

The turbulence resulting from the buoyancy instability generates a mean magnetic field component ${\overline{\mathit{B}}}_{\mathit{x}}$ from an initial ${\overline{\mathit{B}}}_{\mathit{y}}$ which is also modified compared to its initial shape (see Fig. 10). We employ the quasi-kinematic test-field method to calculate transport coefficients like the α and η tensors which describe this process. So far, the method has mostly been applied in situations where a hydrodynamic background was already present in absence of the mean magnetic field (see, e.g. Brandenburg et al. 2008a,b,c). Here, in contrast, the (magnetohydrodynamic) turbulence results entirely from the instability of a pre-existing mean magnetic field, B_y0(z). In other words, our simulations do not posses a kinematic stage in which the influence of $\overline{{B}}$ would be negligible. The applicability of QKTF to such a situation has been questioned by Courvoisier et al. (2010). We emphasize, however, that Eq. (15)continues to be valid also in that case and hence all conclusions drawn from it, because the decisive criterion for its validity is that there exists no magnetic turbulence in the absence of the mean magnetic field, that is, b → 0 if $\overline{{B}}\mathrm{\to }\mathbf{0}$ (see Rheinhardt & Brandenburg 2010, for a more thorough explanation). As this is the case here the QKTF method is applicable. Furthermore, we explicitly check the applicability of QKTF to this problem by reconstructing the turbulent EMF from the turbulent transport coefficients we calculate. The only peculiarity occurring is the fact that all components of α and η vanish for $\mathrm{0}\mathrm{\le }{\overline{\mathit{B}}}_{\mathrm{rms}}\mathrm{\le }{\overline{\mathit{B}}}_{\mathrm{threshold}}$, because fluctuating velocity and magnetic fields develop only after the instability has set in. Another aspect not considered in most previous test-field studies is the strong intrinsic inhomogeneity of the turbulence not only as a consequence of the strong z dependence of $\overline{{B}}$, but also due to the stratified density background. Thus the transport coefficients need to be determined as z dependent quantities.

In the next section we demonstrate that the test-field method works reasonably well for our problem. Note that to calculate the transport coefficients in addition to the usual MHD equations four additional evolution equations of the form (15)for four independent test-fields have to be solved. Hence the test-field runs are computationally almost thrice as expensive. We have thus reduced resolution to 64³ grid points for all these runs.

Fig. 7 Relative density perturbation, δρ/ρi (top) and relative temperature perturbation δT/Ti (bottom), with δρ = ρ − ρi, δT = T − Ti and ρi(z), Ti(z) taken from Eq. (7), in the plane y = 0 at t = tsat ~ 2tA0 for the run B128h. Both plots overlaid with contours of By (solid lines).

Fig. 8 Top: relative current helicity J·B/JrmsBrms for run B128d at t = tsat ~ 2tA0; arrows: vx, vz. Bottom: By/B0 for the same run; arrows: Bx, Bz. Both panels show the plane y = 0.

Fig. 9 Top: dependence of the total relative current helicity ϵJ on inverse Roberts number. Size of circles codes for value of Rm (length based). Bottom: dependence of ${\mathit{\epsilon }}_{\mathit{J}}\mathrm{=}\overline{{J}\mathrm{·}{B}}\mathit{/}{\mathit{J}}_{\mathrm{rms}}{\mathit{B}}_{\mathrm{rms}}$ on z after saturation for runs B128e (solid, Rb-1 = 0.25), B128h (dashed, Rb-1 = 0.5), and B128c (dotted, Rb-1 = 1).

Fig. 10 Time-height diagram for ${\overline{\mathit{B}}}_{\mathit{x,y}}$ normalized on B0 for run B128c in Table 2 (z extent of the box clipped). Note the difference of two orders of magnitude between ${\overline{\mathit{B}}}_{\mathit{y}}$ and ${\overline{\mathit{B}}}_{\mathit{x}}$.

3.5.1. Reconstruction of the mean EMF

To validate the test-field method we first confirm that the quantity $\overline{\mathbf{ℰ}}$, taken directly from the direct numerical simulation (DNS), can be reproduced by employing the relation (11)between $\overline{\mathbf{ℰ}}$ and $\overline{{B}}$ with the tensors $\hat{{\alpha }}$ and $\hat{{\eta }}$ determined using the quasi-kinematic test-field method. In mathematical terms, $\overline{\mathrm{ℰ}}{\mathit{i}}_{\mathrm{R}}^{}\mathrm{\left(}\mathit{z}\mathrm{;}\overline{{B}}\mathrm{\right)}\mathrm{=}\sum _{{\mathit{k}}^{\mathrm{c}}}\mathrm{𝒦̂}\begin{array}{}\mathrm{c}\\ \mathrm{⟩}\mathrm{}\end{array}\mathrm{(}\mathrm{‡}\mathit{,}{\mathrm{\parallel }}^{\mathrm{c}}\mathrm{;}\overline{{B}}\mathrm{)}\hspace{0.17em}\mathrm{ℬ̂}\begin{array}{}\mathrm{c}\\ \mathrm{}\end{array}\mathrm{\left(}{\mathrm{\parallel }}^{\mathrm{c}}\mathrm{\right)}\mathrm{+}\sum _{{\mathrm{\parallel }}^{\mathrm{s}}}\mathrm{𝒦̂}\begin{array}{}\mathrm{s}\\ \mathrm{⟩}\mathrm{}\end{array}\mathrm{(}\mathrm{‡}\mathit{,}{\mathrm{\parallel }}^{\mathrm{s}}\mathrm{;}\overline{{B}}\mathrm{)}\hspace{0.17em}\mathrm{ℬ̂}\begin{array}{}\mathrm{s}\\ \mathrm{}\end{array}\mathrm{\left(}{\mathrm{\parallel }}^{\mathrm{s}}\mathrm{\right)}$(20)with $\begin{array}{ccc}& & \hspace{0.17em}\begin{array}{}\begin{array}{}\end{array}\mathrm{𝒦̂}\begin{array}{}\mathrm{c}\\ \mathrm{⟩}\mathrm{}\end{array}\mathrm{=}\mathit{\alpha ̂}\mathit{ij}\mathrm{\left(}\mathit{z,}{\mathit{k}}^{\mathrm{c}}\mathrm{;}\overline{{B}}\mathrm{\right)}\cos \mathrm{\left(}{\mathit{k}}^{\mathrm{c}}\hspace{0.17em}\begin{array}{}\mathrm{􏽥}\\ \mathit{z}\end{array}\mathrm{\right)}\begin{array}{}\end{array}\mathrm{-}\mathit{\eta ̂}\mathit{il}\mathrm{\left(}\mathit{z,}{\mathit{k}}^{\mathrm{c}}\mathrm{;}\overline{{B}}\mathrm{\right)}\hspace{0.17em}{\mathit{ϵ}}_{\mathit{lj}\mathrm{3}}\hspace{0.17em}\hspace{0.17em}{\mathit{k}}^{\mathrm{c}}\sin \mathrm{\left(}{\mathit{k}}^{\mathrm{c}}\hspace{0.17em}\begin{array}{}\mathrm{􏽥}\\ \mathit{z}\end{array}\mathrm{\right)}\mathit{,}\\ \begin{array}{}\end{array}\mathrm{𝒦̂}\begin{array}{}\mathrm{s}\\ \mathrm{⟩}\mathrm{}\end{array}\mathrm{=}\mathit{\alpha ̂}\mathit{ij}\mathrm{\left(}\mathit{z,}{\mathit{k}}^{\mathrm{s}}\mathrm{;}\overline{{B}}\mathrm{\right)}\sin \mathrm{\left(}{\mathit{k}}^{\mathrm{s}}\hspace{0.17em}\begin{array}{}\mathrm{􏽥}\\ \mathit{z}\end{array}\mathrm{\right)}\begin{array}{}\end{array}\mathrm{+}\mathit{\eta ̂}\mathit{il}\mathrm{\left(}\mathit{z,}{\mathit{k}}^{\mathrm{s}}\mathrm{;}\overline{{B}}\mathrm{\right)}\hspace{0.17em}{\mathit{ϵ}}_{\mathit{lj}\mathrm{3}}\hspace{0.17em}\hspace{0.17em}{\mathit{k}}^{\mathrm{s}}\cos \mathrm{\left(}{\mathit{k}}^{\mathrm{s}}\hspace{0.17em}\begin{array}{}\mathrm{􏽥}\\ \mathit{z}\end{array}\mathrm{\right)}\mathit{,}\end{array}\\ & & \\ & & {\mathit{k}}^{\mathrm{c}}\mathrm{=}\frac{\mathrm{(}\mathrm{2}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{\pi }}{{\mathit{L}}_{\mathit{z}}}\mathit{,} {\mathit{k}}^{\mathrm{s}}\mathrm{=}\frac{\mathrm{2}\mathit{k\pi }}{{\mathit{L}}_{\mathit{z}}}\mathit{,} \mathit{k}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathit{...}\mathit{,}\\ & & \hspace{0.17em}\begin{array}{}\mathrm{􏽥}\\ \mathit{z}\end{array}\mathrm{=}\mathit{z}\mathrm{-}{\mathit{z}}_{\mathrm{0}}\mathrm{-}\frac{{\mathit{L}}_{\mathit{z}}}{\mathrm{2}}\mathit{,}\end{array}$(21)where the superscript R indicates reconstruction. Here, the boundary condition for B gives rise to the selection of discrete cosine and sine modes with wavenumbers k^c and k^s, respectively. The additional argument $\overline{{B}}$ is to indicate that the quantities , related to $\hspace{0.17em}\hat{{G}}\mathrm{\left(}{x}\mathit{,}{k}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}{x}}$ in Eq. (13), as well as the tensors $\hat{{\alpha }}$ and $\hat{{\eta }}$, are valid just for that mean field $\overline{{B}}$ which is present in the main run. As a consequence, the reconstruction of the mean EMF can be successful only when employing exactly this $\overline{{B}}$ in (21). That is, the mean field representation of the turbulence in terms of $\hat{{\alpha }}$ and $\hat{{\eta }}$ has, at this level, merely descriptive rather than predictive potential. In order to possess the latter, α and η have to be formulated as functionals of $\overline{{B}}$ and z. Consequently, the mean EMF will appear as a nonlinear functional of $\overline{{B}}$. Since it is additionally also non-local, we are confronted with a task of considerable mathematical complication which is not addressed here.

Let us denote $\overline{\mathbf{ℰ}}{\mathrm{R}}^{}$ as the reconstructed EMF according to Eq. (20) truncated by k′ ≤ K′, with k^c,s = 2πk′/L_z. Here k′ can now take both integer and half-integer values where the integer (half-integer) values of k′ correspond to the family of sine (cosine) modes in Eq. (21). An initial estimate of K′ required for a reasonable reconstruction of $\overline{\mathbf{ℰ}}$ was obtained from the power spectra of both ${\overline{\mathit{B}}}_{\mathit{x}}$ and ${\overline{\mathit{B}}}_{\mathit{y}}$. It turned out that ${\overline{\mathit{B}}}_{\mathit{x}}$ has significant spectral power up until k′ = 16, whereas for ${\overline{\mathit{B}}}_{\mathit{y}}$ the power spectrum has levelled off already at k′ = 8.

The components of the tensors $\hat{{\alpha }}$ and $\hat{{\eta }}$ also show rather different spectral behavior, both in the midplane of the magnetic layer and near the midplane of the box as seen in Fig. 11. Given that the spectrum of the mean field clearly converges to zero, this guarantees that the mean EMF will also spectrally converge. The tensor components either appear to converge to a constant value or to zero. For coefficients like and , converging to a constant, this asymptotic, say $\mathit{\alpha ̂}\begin{array}{}\mathrm{0}\\ \mathrm{12}\end{array}$, can in principle be separated and would appear as a part $\mathit{\alpha ̂}\begin{array}{}\mathrm{0}\\ \mathrm{12}\end{array}{\mathit{\delta }}^{\mathrm{3}}\mathrm{(}{x}\mathrm{-}{x}{}^{\mathrm{\prime }}\mathrm{)}$ of the convolution kernel of (12), representing a local part of the relation between $\overline{\mathbf{ℰ}}$ and $\overline{{B}}$. This is, however, not done here. Instead we simply employ the k dependent coefficients up to the maximum required by the spectra of $\overline{{B}}$, that is, k′ = 16. Note that the values for k′ = 0 are not relevant here as, due to the boundary conditions, $\overline{{B}}$ does not possess a k′ = 0 contribution.

The α tensor can be decomposed in symmetric and antisymmetric parts, the latter being representable by a turbulent pumping velocity γ, which gives rise to the term ${\gamma }\mathrm{\times }\overline{{B}}$ in the mean EMF. Since ${\overline{\mathit{B}}}_{\mathit{z}}\mathrm{=}\mathrm{0}$, the only relevant component of γ is here . Analytic results indicate that in a wide range of situations, the turbulent pumping is directed away from the region of strong turbulence (“turbulent diamagnetism”, see Krause & Rädler 1980). is presented in Fig. 11 together with in one panel. As decreases in modulus beyond k′ = 1 and even changes sign at k′ = 8 while increases in modulus, but preserves its sign, we observe a sign change in at k′ = 4 and k′ = 3 in the midplane of the magnetic layer and above it, respectively (see also Fig. 18 below). A similar dependence of turbulent pumping on wavenumber has been found by Käpylä et al. (2009) in DNS of convection. Consequently, magnetic fields with scales up to ≈ (5...6)H_B, that is, k′ > 4 are pumped into the layer, but those with larger scales are pumped away from it the former being contrary to the standard concept of “turbulent diamagnetism”. It is thus difficult to comment on the transport of the total $\overline{{B}}$ by $\hat{{\gamma }}$. Only if the pumping were oriented away from the magnetic layer for all the wavenumbers of the dominating constituents in $\overline{{B}}$ it would lead to a broadening of the initial layer i.e., a reduction of $\mathit{\partial }{\overline{\mathit{B}}}_{\mathit{y}}\mathit{/}\mathit{\partial z}$ and would hence inhibit the instability. With respect to its saturation, however, the strong turbulent magnetic diffusion given by is likely to be more important, as for low k′ it reaches ≈ 20 times the molecular value, cf. Fig. 11.

The result of the assembly of $\overline{\mathbf{ℰ}}{\mathrm{R}}^{}$ from (20)with (21), is presented in Fig. 12, middle column. From simple visual inspection we find it to be a faithful reproduction of $\overline{\mathbf{ℰ}}$ from the DNS shown in the left column. Clearly, a naive application of the test-field procedure with harmonic test fields with only the lowest k′ = 0.5 results in an inadequate description as shown in the right column.

Fig. 11 Dependence of $\hat{{\alpha }}$ and $\hat{{\eta }}$ time averaged between t = tsat to t = tsat + tA0 on the test-field wavenumber k′ in the midplane of the magnetic layer (z = zB) and near the midplane of the box (z = 0.21Lz) for run TF30+ of Table 3. Thick lines in panel for : turbulent pumping coefficient . Rescaling of to the molecular value η by the factor UrmsHB/η = 21.5. Integer and half-integer values of k′ belong to sine and cosine modes in $\overline{{B}}$, respectively. Note that k′ = 0 refers to constant and linear test fields and that the corresponding coefficients do not enter the $\overline{\mathbf{ℰ}}$ – $\overline{{B}}$ relation for the given setup.

We define two measures for the quality of the mean EMF reconstruction at a given K′ namely ${\mathit{\chi }}_{{\mathit{K}}^{\mathrm{\prime }}}^{\mathrm{2}}$ and the correlation coefficient r_K′ defined as ${\mathit{\chi }}_{{\mathit{K}}^{\mathrm{\prime }}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{⟨}\mathrm{(}\overline{\mathrm{ℰ}}{\mathit{x,y}}_{}\mathrm{-}\overline{\mathrm{ℰ}}{\mathit{x,y}}_{\mathrm{R}}^{}{\mathrm{)}}^{\mathrm{2}}{\mathrm{⟩}}_{\mathit{z,t}}}{\mathrm{⟨}\overline{\mathrm{ℰ}}{\mathit{x,y}}_{\mathrm{2}}^{}{\mathrm{⟩}}_{\mathit{z,t}}}\mathit{,}{\mathit{r}}_{{\mathit{K}}^{\mathrm{\prime }}}\mathrm{=}\frac{\mathrm{⟨}\overline{\mathrm{ℰ}}{\mathit{x,y}}_{}\mathrm{\times }\overline{\mathrm{ℰ}}{\mathit{x,y}}_{\mathrm{R}}^{}{\mathrm{⟩}}_{\mathit{z,t}}}{\sqrt{\mathrm{⟨}\overline{\mathrm{ℰ}}{}_{\mathit{x,y}}^{\mathrm{2}}{\mathrm{⟩}}_{\mathit{z,t}}\mathrm{⟨}\overline{\mathrm{ℰ}}{}_{\mathit{x,y}}^{{\mathrm{R}}^{\mathrm{2}}}{\mathrm{⟩}}_{\mathit{z,t}}}}\mathit{,}$(22)where the subscript “z,t” denotes that the averaging has been carried out over the vertical coordinate z as well as over the temporal range 1.2t_A0 ≤ t ≤ 3.4t_A0. The relative error of the reconstruction, ${\mathit{\chi }}_{{\mathit{K}}^{\mathrm{\prime }}}^{\mathrm{2}}$, and the correlation coefficient, r_K′, are plotted in Fig. 13 as a function of the truncation wavenumber K′. The ${\mathit{\chi }}_{{\mathit{K}}^{\mathrm{\prime }}}^{\mathrm{2}}$ reach a minimum value of about 35% and 30% for $\overline{\mathrm{ℰ}}{\mathit{x}}_{}$ and $\overline{\mathrm{ℰ}}{\mathit{y}}_{}$, respectively, and level off around K′ = 8. This implies that including higher harmonic test fields beyond k′ = 7 does not improve the reconstructed EMF. Similarly, r_K′ for $\overline{\mathrm{ℰ}}{\mathit{x}}_{\mathrm{R}}^{}$$\mathrm{\left(}\overline{\mathrm{ℰ}}{\mathit{y}}_{\mathrm{R}}^{}\mathrm{\right)}$ converges to a value of 0.98 (0.93) at K′ = 4 (8). It is important to note that even though the tensor components and do not converge with increasing k′, the reconstructed EMFs do as a consequence of the sufficiently fast convergence of $\overline{{B}}$ with k′. Also, including transport coefficients for k′ ≥ 8 does not improve the reconstruction any further. Clearly, one reason behind the discrepancies is that we have neglected memory effects in the turbulent transport coefficients (see Hubbard & Brandenburg 2009). It can be particularly important in the present situation as we are not in a statistically stationary regime. Apart from this, enhancement of resolution might further improve the results.

Fig. 12 Reconstruction of the mean EMF for the run TF30+ using $\hat{{\alpha }}$ and $\hat{{\eta }}$ from the test-field method. Top: $\overline{\mathrm{ℰ}}{\mathit{x}}_{}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}$, bottom: $\overline{\mathrm{ℰ}}{\mathit{y}}_{}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}$, both scaled by 10-4vA0B0. Left: directly from $\overline{{u}\mathrm{\times }{b}}$. Middle: reconstruction using all contributions k′ = 0.5,1,1.5,...16 in (20). Right: same as before, but using only the k′ = 0.5 contribution.

Fig. 13 Quality of the mean EMF reconstruction as a function of the truncation wavenumber K′: ${\mathit{\chi }}_{{\mathit{K}}^{\mathrm{\prime }}}^{\mathrm{2}}$ (left) and correlation rK′ (right) calculated for $\overline{\mathrm{ℰ}}{\mathit{x}}_{\mathrm{R}}^{}$ (solid) and $\overline{\mathrm{ℰ}}{\mathit{y}}_{\mathrm{R}}^{}$ (dashed) using Eq. (22).

3.5.2. Dependence of the transport tensors on inclination

In view of the solar dynamo problem it is important to look at $\hat{{\alpha }}$ and $\hat{{\eta }}$ as functions of the rotational inclination θ or latitude λ = 90° − θ. We expect the growth rate of the instability to increase from the equator to the pole (Schmitt 2003). This can be explained by the buoyant nature of the turbulence, for which vertical motions are essential: At the poles, the effect of the Coriolis force on them is weakest whereas they are strongly deflected at the equator. This is indeed confirmed by Fig. 14, where the growth rate is seen to decrease continuously when changing θ from 0° (pole) towards 90° (equator). The corresponding runs are summarized in Table 3.

Let us now consider the symmetry properties of our solution with respect to λ = 0, which is the solar equator. Moving from the northern hemisphere at λ to the southern at − λ, that is changing θ to 180° − θ, but keeping all other problem parameters invariant, is equivalent to inverting the sign of Ω_z. The same can be accomplished by reflecting the corresponding rigid rotation about the plane x = 0. Hence we can construct the solution (ρ,U,B,s) of (2)–(5)for − λ simply by reflecting the solution for λ properly about the same plane. Under this reflection, polar vectors like velocity transform as, $\mathrm{\left\{}{\mathit{U}}_{\mathit{x}}\mathit{,}{\mathit{U}}_{\mathit{y}}\mathit{,}{\mathit{U}}_{\mathit{z}}\mathrm{\right\}}\mathrm{\left(}\mathit{x,y,z}\mathrm{\right)}\mathrm{\to }\mathrm{\{}\mathrm{-}{\mathit{U}}_{\mathit{x}}\mathit{,}{\mathit{U}}_{\mathit{y}}\mathit{,}{\mathit{U}}_{\mathit{z}}\mathrm{\}}\mathrm{(}\mathrm{-}\mathit{x,y,z}\mathrm{)}\mathit{,}$(23)and axial vectors like the magnetic field as $\mathrm{\left\{}{\mathit{B}}_{\mathit{x}}\mathit{,}{\mathit{B}}_{\mathit{y}}\mathit{,}{\mathit{B}}_{\mathit{z}}\mathrm{\right\}}\mathrm{\left(}\mathit{x,y,z}\mathrm{\right)}\mathrm{\to }\mathrm{\{}{\mathit{B}}_{\mathit{x}}\mathit{,}\mathrm{-}{\mathit{B}}_{\mathit{y}}\mathit{,}\mathrm{-}{\mathit{B}}_{\mathit{z}}\mathrm{\}}\mathrm{(}\mathrm{-}\mathit{x,y,z}\mathrm{)}\mathit{.}$(24)(Note that the gravitational acceleration is invariant under this reflection.) Hence, for the initial magnetic field, B_y0(z), the transition to − λ requires only a sign inversion. However, since the induction equation is linear in B, and Lorentz force as well as Ohmic dissipation are quadratic, inverting the sign of B_y0(z) would just transform the solution { ρ,U,B,s } to { ρ,U, − B,s } , that is, would leave the turbulence essentially unchanged and can be omitted. Moreover, as the transport coefficients which express the correlation properties of the turbulent velocity u are functions of z only, the reflection operation cannot change their magnitudes. With respect to possible sign inversions in the coefficients upon reflection, we note that $\overline{\mathbf{ℰ}}$ and $\overline{{J}}$, being polar vectors, invert the sign of their x components under reflection, but keep their y components unchanged. The axial vector $\overline{{B}}$ behaves just the opposite way. Thus, we have for i = 1,2 (no summation) and for i ≠ j, whereas for i = 1,2 and for i ≠ j when moving from λ to − λ. Consequently, it appears that the results for the southern hemisphere can be derived from those for the northern by simple operations. Strictly speaking however, this is only true when the initial condition for U is also reflected upon the transition from λ to − λ. From a naive point of view we might suppose that omitting this reflection can hardly be of any importance, because the initial conditions are anyway random. But we have found this not to be true.

In Fig. 15 we show z–λ plots of $\overline{{B}}$ and $\overline{\mathbf{ℰ}}$, constructed from the runs with Ω ≠ 0. We expect ${\overline{\mathit{B}}}_{\mathit{x}}$ and $\overline{\mathrm{ℰ}}{\mathit{y}}_{}$ to be antisymmetric about the equator. This is indeed true for higher latitudes but not in the neighbourhood of the equator where they are rather symmetric. Violation of the expected symmetry properties can also be inferred from the extrema of $\overline{\mathrm{ℰ}}{\mathit{y}}_{}$ given in Table 3.

Symmetry can formally be restored as follows. From Fig. 15 we construct another set of plots by “symmetrizing” ${\overline{\mathit{B}}}_{\mathit{y}}$ and $\overline{\mathrm{ℰ}}{\mathit{x}}_{}$, but “anti-symmetrizing” ${\overline{\mathit{B}}}_{\mathit{x}}$ and $\overline{\mathrm{ℰ}}{\mathit{y}}_{}$, that is, setting ${\overline{\mathit{B}}}_{\mathit{y}}\mathrm{\left(}\mathit{z,\lambda }\mathrm{\right)}$ to $({\overline{\mathit{B}}}_{\mathit{y}}\mathrm{\left(}\mathit{z,\lambda }\mathrm{\right)}\mathrm{+}{\overline{\mathit{B}}}_{\mathit{y}}\mathrm{(}\mathit{z,}\mathrm{-}\mathit{\lambda }\mathrm{)})\mathit{/}\mathrm{2}$ and ${\overline{\mathit{B}}}_{\mathit{x}}\mathrm{\left(}\mathit{z,\lambda }\mathrm{\right)}$ to $({\overline{\mathit{B}}}_{\mathit{x}}\mathrm{\left(}\mathit{z,\lambda }\mathrm{\right)}\mathrm{-}{\overline{\mathit{B}}}_{\mathit{x}}\mathrm{(}\mathit{z,}\mathrm{-}\mathit{\lambda }\mathrm{)})\mathit{/}\mathrm{2}$, analogously for $\overline{\mathbf{ℰ}}$. The result is shown in Fig. 16. In the same way, we have obtained the transport coefficients shown in Figs. 17 and 18 for the test-field wavenumber k′ = 0.5, that is, for the largest possible wavelength (cosine mode). We find that the moduli of all turbulent transport coefficients increase monotonically for all z as we go from the equator to the poles, except and which exhibit local extrema at |λ| ≈ 30°.

Fig. 14 Dependence of the instability on rotational inclination θ in terms of rms value of generated field components $\mathrm{⟨}{\mathit{B}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{B}}_{\mathit{z}}^{\mathrm{2}}\mathrm{⟩}$ for the runs TF0+, TF30+, TF60+ and TF90+ in Table 3.

A careful look at the “symmetrized” and “anti-symmetrized” plots reveals a discontinuity at the equator best seen in the plots of (Figs. 17). This phenomenon is counterintuitive at first: if we regard the computational box as a local area within a spherical body, changing λ from positive to negative values implies travelling along a meridian across the equator. If quantities like ${\overline{\mathit{B}}}_{\mathit{x}}$ and $\overline{\mathrm{ℰ}}{\mathit{y}}_{}$ change their sign there, they should vanish for continuity reasons and the same should hold true for and . Instead of looking at the “symmetrized” quantities we can reflect the initial velocity ${U}{\mathrm{0}}_{}$ according to (23)upon the transition λ   →   − λ. This restores exact symmetry/antisymmetry about the equator and we again observe a discontinuity at the equator.

In conclusion, near the equator the initial conditions play a crucial role in determining the saturated values of ${\overline{\mathit{B}}}_{\mathit{x}}$ and $\overline{\mathrm{ℰ}}{\mathit{y}}_{}$ and consequently the values of the turbulent transport coefficients. As the coefficients in are typically connected with the kinetic helicity density h_K, we have calculated it for the saturated states. Indeed, at the equator the solutions resulting from the initial condition ${U}{\mathrm{0}}_{}$ and its reflected counterpart show opposite signs of h_K. This is an example of spontaneous symmetry breaking with bifurcation into helical solutions of opposite handedness which enables finite values of ${\overline{\mathit{B}}}_{\mathit{x}}$ and $\overline{\mathrm{ℰ}}{\mathit{y}}_{}$ as well as of at the equator. (Note, that are not related to the non-vanishing helicity at the equator and consequently show no discontinuity in their anti-symmetrized profiles in Fig. 17.) As a prerequisite for such a bifurcation, the unstable eigenmodes of the buoyancy instability must be degenerate. That is, two linearly independent modes with opposite signs of helicity, but identical growth rate must exist. This has been noticed for Ω ≠ 0 at the equator (Thelen 2000) and also holds true for Ω = 0³. A detailed discussion of this phenomenon for Ω = 0 is the subject of a separate paper Chatterjee et al. (2011).

Fig. 15 ${\overline{\mathit{B}}}_{\mathit{x,y}}\mathrm{\left(}\mathit{z,\lambda }\mathrm{\right)}$, $\overline{\mathrm{ℰ}}{\mathit{x,y}}_{}\mathrm{\left(}\mathit{z,\lambda }\mathrm{\right)}$ averaged between t = tsat and tsat + tA0 from the runs in Table 3 with Ω ≠ 0. (Initial velocity identical for all latitudes.) Dashed line: initial position of the magnetic layer.

Fig. 16 Same as Fig. 15, but with symmetrized ${\overline{\mathit{B}}}_{\mathit{y}}\mathrm{\left(}\mathit{z,\lambda }\mathrm{\right)}$, $\overline{\mathrm{ℰ}}{\mathit{x}}_{}\mathrm{\left(}\mathit{z,\lambda }\mathrm{\right)}$ and anti-symmetrized ${\overline{\mathit{B}}}_{\mathit{x}}\mathrm{\left(}\mathit{z,\lambda }\mathrm{\right)}$, $\overline{\mathrm{ℰ}}{\mathit{y}}_{}\mathrm{\left(}\mathit{z,\lambda }\mathrm{\right)}$. For explanation see text.

Fig. 17 $\hat{{\alpha }}\mathrm{\left(}\mathit{z,\lambda }\mathrm{\right)}$ and $\hat{{\eta }}\mathrm{\left(}\mathit{z,\lambda }\mathrm{\right)}$ from the runs in Table 3 with Ω ≠ 0, calculated using test fields with k′ = 0.5 and averaged between t = tsat and tsat + tA0. symmetrized, anti-symmetrized analogously to $\overline{{B}}$, $\overline{\mathbf{ℰ}}$ in Fig. 16. $\hat{{\alpha }}$ scaled by Urms for the run TF30+, $\hat{{\eta }}$ scaled by the molecular diffusivity η. Dashed line: initial position of the magnetic layer.

Fig. 18 Pumping velocity scaled with Urms for the run TF30+, as a function of height and latitude for k′ = 0.5 (left) and k′ = 4.5 (right).

3.5.3. Dependence on initial magnetic field strength

In standard mean-field theory for a prescribed hydrodynamic background the turbulent transport coefficients usually decrease as the mean magnetic field increases (“quenching”). The present problem is different, however, because the instability and hence the turbulence is just caused by the initial (mean) magnetic field. In Fig. 19 we show the z averages of and at t^sat as functions of the initial magnetic field ${\mathit{B}}_{\mathrm{0}}\mathrm{\propto }\hspace{0.17em}\begin{array}{}{}^{\mathrm{􏽥}}\\ \mathit{\beta }\end{array}{\mathrm{0}}_{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}^{}$. These components have been selected because they essentially do not invert sign in the domain during the entire evolution. Because of its importance for the generation of ${\overline{\mathit{B}}}_{\mathit{x}}$ from ${\overline{\mathit{B}}}_{\mathit{y}}$ we have also given in this figure. However, as this coefficient shows a significant sign change with respect to z, its maximum value with respect to z and the time span 0 < t < 2t^sat was plotted.

Clearly, , , and increase with the initial magnetic field strength supporting earlier ideas of a possible “anti-quenching” in the case of the buoyancy instability (see, e.g., Brandenburg et al. 1998). Note, however, that the dependence on the initial field strength in Fig. 19 might differ from the dependence on the local field strength to which the term “quenching” usually refers.

Fig. 19 Vertical averages of (solid) and (dashed) for k′ = 0.5 at t = tsat and maximum of with respect to all z and t (dotted) as functions of $\hspace{0.17em}\begin{array}{}{}^{\mathrm{􏽥}}\\ \mathit{\beta }\end{array}\begin{array}{}\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}\\ \mathrm{0}\end{array}\mathrm{\propto }{\mathit{B}}_{\mathrm{0}}$.

4. Conclusions

We have studied in detail the generation of the α effect due to the buoyancy instability of a horizontal magnetic layer in a stratified atmosphere by using direct numerical simulations. We find that both the magnetic energy and the current helicity (twist) in the system increase monotonically with the ratio of thermal conductivity to magnetic diffusivity, i.e. the Roberts number Rb (Fig. 4). This agrees with earlier analytical work of Gilman (1970) and Acheson (1979), as well as numerical work of Silvers et al. (2009), who find that efficient thermal diffusion or heat exchange can destabilize a stable stratification. The increase of the resulting twist with Rb is an important result since the buoyancy instability would thus produce twisted flux tubes from a magnetic layer, if it existed in the overshoot layer of the Sun where Rb ≫ 1. Vasil & Brummel (2008) also report the formation of twisted flux tubes from a horizontal magnetic layer, but in their case it is due to the action of shear on a weak vertical magnetic field. We further find that the growth rate of the buoyancy instability is reduced in presence of rotation compared to the case with Ω = 0.

Most of the earlier work on the magnetic buoyancy instability resorted to either the magnetic Boussinesq approximation or the anelastic approximation. Here we have performed fully compressible numerical calculations (for differences between anelastic and compressible calculations see also Berkoff et al. 2010). Let us now relate our results to corresponding earlier calculations. The work of Wissink et al. (2000) comes closest to our own as they not only included rotation, but also employed a thin magnetic layer. However, both thermal and magnetic boundary conditions were different from ours and the initial magnetostatic equilibrium was established by lowering the density in the layer, thus making it immediately Rayleigh-Taylor unstable. Moreover, their simulations were done for different parameters (our values in brackets): Rb = 17 [0.25...1] (bottom), density contrast 5.8 [223] , Ta_M ≈ 10⁷ [10¹⁰] and Lu₀ ≈ 50 [500...600]. Hence, a quantitative comparison is not possible, although the occurrence of similar tubular structures in the saturated stage (see our Fig. 6 and their Fig. 4) indicates that those might be a robust feature of magnetic-buoyancy driven turbulence.

Thelen (2000) has also performed fully compressible calculations and provides the α coefficient, but only as a result from individual unstable eigenmodes in the linear stage. In contrast, our calculations of α and magnetic diffusivity were done in the saturated stage and are hence not comparable to the results of Thelen (2000). Furthermore, in his analysis, the relation α_yy ⟨ B_0y ⟩ = ⟨ u × b ⟩ _y was employed with ⟨ · ⟩ being either a volume or horizontal average. In the first case it yields the correct α_yy, but this is valid only for uniform mean fields. In the second case the result is not a measure of the actual α_yy, because the contribution of $\mathit{\partial }{\overline{\mathit{B}}}_{\mathrm{0}\mathit{y}}\mathit{/}\mathit{\partial z}$ to $\overline{\mathbf{ℰ}}$ is ignored. Hence, in both cases no description of the mean EMF that would be useful for a mean-field dynamo model is provided. Further, he admits only a local relationship between $\overline{\mathbf{ℰ}}$ and $\overline{{B}}$.

Kersalé et al. (2007) considered a non-rotating box with a linear vertical profile of the horizontal magnetic field, which is supported by the boundary condition. Perturbations to this state were subjected to perfect conductor boundary conditions, and the background state was a polytrope with index 1.6 along with Pr = 1/Rb = 0.02. Hence, both their model and the employed parameters are too different from ours to allow a meaningful comparison of the results. Nevertheless, their magnetic field structure for random initial conditions (their Fig. 4) resembles our finding in Fig. 6.

We have run our simulations only until the time taken by the initial magnetic layer to break up due to the back-reaction of the unstable modes and ohmic diffusion. In the absence of any other forcing such as a strong shear, the buoyancy instability is found to be incapable of sustaining itself past the break-up phase since the scale height of the magnetic layer becomes comparable to that of the density. We may hence say that strong shear is not imperative to the production of tubular structures from the magnetic layer, but is likely to play a key role in keeping the layer from breaking up. It may also be possible that turbulent downward pumping (Nordlund et al. 1992) arrests the decay of a magnetic layer in the overshoot region. However, it is not yet clear if such a layer exists at all in the real Sun and, moreover, if it is subject to the buoyancy instability there.

For mean fields defined by horizontal averaging we have “measured” the turbulent transport coefficients using the technique of the quasi-kinematic test-field method. In order to prove that the $\hat{{\alpha }}$ and $\hat{{\eta }}$ tensors obtained from this method are reasonably accurate, we have verified the agreement between $\overline{{u}\mathrm{\times }{b}}$ from the DNS and the representation of $\overline{\mathbf{ℰ}}$ by the convolution ${\alpha }\mathrm{◦}\overline{{B}}\mathrm{-}{\eta }\mathrm{◦}\overline{{J}}$, with corresponding tensors $\hat{{\alpha }}\mathrm{\left(}\mathit{k}\mathrm{\right)}$ and $\hat{{\eta }}\mathrm{\left(}\mathit{k}\mathrm{\right)}$ obtained for harmonic test fields with wavenumbers 0 ≤ (L_z/2π)k ≤ 16. A technique for the quantitative assessment of this agreement is presented. We find that, even in presence of magnetically driven turbulence, the obtained $\hat{{\alpha }}$ and $\hat{{\eta }}$ provide a reasonably accurate description of the turbulent EMF; see Figs. 12 and 13 (although the relative errors are still about 30   ...     35%). This is an important outcome of our study.

We find that $\overline{\mathrm{ℰ}}{\mathit{y}}_{}$, determined using a harmonic test field with the lowest wavenumber that fits in the vertical extent of the box, already comprises a considerable part of the total EMF. Hence it can be enlightening to look at the turbulent transport coefficients obtained by the QKTF method for this wavenumber only. At finite Ω their dependence on (solar or stellar) latitude λ is of particular interest. The component contributes to the generation of ${\overline{\mathit{B}}}_{\mathit{x}}$ from the strong initial field ${\overline{\mathit{B}}}_{\mathit{y}}$ in the layer. The off-diagonal components of $\hat{{\alpha }}$ contribute to a vertical turbulent pumping velocity directed away from the region of turbulence surrounding the magnetic layer. The effect of this pumping systematically broadens magnetic structures with increasing latitude. We find that all transport coefficients except and increase with latitude and are significantly reduced near the equator due to the suppressing effect of the Coriolis force on the instability. For the first time the turbulent magnetic diffusivity given by the diagonal components of $\hat{{\eta }}$ has been computed for the magnetic buoyancy instability; see Fig. 17. In particular, near the magnetic layer, the diagonal component is 25 times larger than the molecular value η.

The buoyancy instability has the property that the EMF, as indicated by the extrema of $\overline{\mathrm{ℰ}}{\mathit{y}}_{}\mathit{/}{\mathit{v}}_{\mathit{A}\mathrm{0}}{\mathit{B}}_{\mathrm{0}}$, increases progressively with the magnitude of the magnetic field in the horizontal layer (compare solid and dashed lines in Fig. 5). Indeed, as a remarkable result, we demonstrate that the most significant coefficients , , and increase with the initial field strength; see Fig. 19. This feature makes the buoyancy instability an attractive candidate for generating an α-effect inside the Sun. Unlike turbulent convection which yields an α that is quenched for strong mean magnetic fields, this one increases with increasing field strength. Our findings support suggestions by Brandenburg et al. (1998) that, if turbulent transport coefficients are caused by flows that are magnetically driven like in the present case or, e.g., in the magneto-rotational instability, then both α and η may increase with magnetic field strength. This trend is sometimes referred to as “anti-quenching” and could be useful to explain the steepness of the observational relation between the ratio of dynamo cycle to rotation frequencies, ω_cyc/Ω and inverse Rossby number for stellar data (Brandenburg et al. 1998; Saar & Brandenburg 1999).

Although highly desirable, modelling of α and η as functionals of position and the mean magnetic field for use in a predictive mean field dynamo model is a difficult proposition that needs to be postponed to future work.

Acknowledgments

We thank Alexander Hubbard for reading the manuscript carefully. The computations have been carried out on the National Supercomputer Centre in Linköping and the Center for Parallel Computers at the Royal Institute of Technology in Sweden. This work was supported in part by the European Research Council under the AstroDyn Research Project No. 227952 and the Swedish Research Council Grant No. 621-2007-4064.

References

All Tables

All Figures

	Fig. 1 The Cartesian simulation domain with respect to spherical coordinates.
In the text

Fig. 2 Time evolution of the runs in Table 2. Upper panel: mean squared values of velocity and generated magnetic field components Bx, Bz of run B128c (Pr = PrM = 1) scaled by $\mathrm{⟨}{\mathit{B}}_{\mathit{y}\mathrm{0}}^{\mathrm{2}}\mathit{/}{\mathit{\mu }}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{0}}\mathrm{⟩}$ and $\mathrm{⟨}{\mathit{B}}_{\mathit{y}\mathrm{0}}^{\mathrm{2}}\mathrm{⟩}$ respectively. Note the clear exponential growth until t ≈ 1.4tA0. Fast oscillations in ⟨ U2 ⟩ until t ≈ tA0 indicate g-modes originating from the initial velocity perturbation. Middle panel: mean squared values of generated magnetic field components for different runs. For legend see lower panel. Lower panel: rms value of By. Prandtl numbers indicated as (Pr, PrM).

In the text

	Fig. 3 Top: velocity components Uy (in color), Ux and Uz (vectors) in the plane y = 0. Bottom: Ux (in color), Uy and Uz (vectors) in the plane x = 0. Both during the linear evolution phase of the run B128a. The initial condition for U was set to be harmonic in x and y with kx = 8 × 2π/Lx and ky = 2π/Ly.
In the text

	Fig. 4 Dependence of growth rate ωI on the inverse Roberts number derived from the runs in Table 2. Solid line: best linear fit. Size of circles codes for the value of Rm (length based, see Table 2).
In the text

Fig. 5 Dependence of the instability on initial magnetic field strength, expressed by $\hspace{0.17em}{\begin{array}{}{}^{\mathrm{􏽥}}\\ \mathit{\beta }\end{array}}_{\mathrm{0}}$, and rotation, expressed by $\mathit{\zeta }\mathrm{=}\hspace{0.17em}{\begin{array}{}{}^{\mathrm{􏽥}}\\ \mathit{\beta }\end{array}}_{\mathrm{0}}{{\mathrm{Ta}}_{\mathrm{M}}}^{\mathrm{1}\mathit{/}\mathrm{2}}$. Upper panel: $\mathrm{⟨}{\mathit{B}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{B}}_{\mathit{z}}^{\mathrm{2}}\mathrm{⟩}$. Lower panel: $\mathrm{⟨}{\mathit{B}}_{\mathit{y}}^{\mathrm{2}}\mathrm{⟩}$. Legend shows $\mathrm{(}\hspace{0.17em}{\begin{array}{}{}^{\mathrm{􏽥}}\\ \mathit{\beta }\end{array}}_{\mathrm{0}}\mathit{,}{\mathrm{10}}^{-5}\mathit{\zeta }\mathrm{)}$. Pr = PrM = 1, colatitude θ = 30°, resolution 643 throughout. Note that the normalization time tA0 is not the same for all curves, but proportional to $\hspace{0.17em}\begin{array}{}{}^{\mathrm{􏽥}}\\ \mathit{\beta }\end{array}{\mathrm{0}}_{\mathrm{1}\mathit{/}\mathrm{2}}^{}$.

In the text

	Fig. 6 Left: volume rendering of the By = 0.1B0 isosurface. Right: field lines, colored according to the value of vorticity Wy for runs B128f (top) and B128e (bottom) at t = 2tA0 (saturated stage). The contours of By are shown on side walls.
In the text

	Fig. 7 Relative density perturbation, δρ/ρi (top) and relative temperature perturbation δT/Ti (bottom), with δρ = ρ − ρi, δT = T − Ti and ρi(z), Ti(z) taken from Eq. (7), in the plane y = 0 at t = tsat ~ 2tA0 for the run B128h. Both plots overlaid with contours of By (solid lines).
In the text

	Fig. 8 Top: relative current helicity J·B/JrmsBrms for run B128d at t = tsat ~ 2tA0; arrows: vx, vz. Bottom: By/B0 for the same run; arrows: Bx, Bz. Both panels show the plane y = 0.
In the text

Fig. 9 Top: dependence of the total relative current helicity ϵJ on inverse Roberts number. Size of circles codes for value of Rm (length based). Bottom: dependence of ${\mathit{\epsilon }}_{\mathit{J}}\mathrm{=}\overline{{J}\mathrm{·}{B}}\mathit{/}{\mathit{J}}_{\mathrm{rms}}{\mathit{B}}_{\mathrm{rms}}$ on z after saturation for runs B128e (solid, Rb-1 = 0.25), B128h (dashed, Rb-1 = 0.5), and B128c (dotted, Rb-1 = 1).

In the text

	Fig. 10 Time-height diagram for ${\overline{\mathit{B}}}_{\mathit{x,y}}$ normalized on B0 for run B128c in Table 2 (z extent of the box clipped). Note the difference of two orders of magnitude between ${\overline{\mathit{B}}}_{\mathit{y}}$ and ${\overline{\mathit{B}}}_{\mathit{x}}$.
In the text

Fig. 11 Dependence of $\hat{{\alpha }}$ and $\hat{{\eta }}$ time averaged between t = tsat to t = tsat + tA0 on the test-field wavenumber k′ in the midplane of the magnetic layer (z = zB) and near the midplane of the box (z = 0.21Lz) for run TF30+ of Table 3. Thick lines in panel for : turbulent pumping coefficient . Rescaling of to the molecular value η by the factor UrmsHB/η = 21.5. Integer and half-integer values of k′ belong to sine and cosine modes in $\overline{{B}}$, respectively. Note that k′ = 0 refers to constant and linear test fields and that the corresponding coefficients do not enter the $\overline{\mathbf{ℰ}}$ – $\overline{{B}}$ relation for the given setup.

In the text

Fig. 12 Reconstruction of the mean EMF for the run TF30+ using $\hat{{\alpha }}$ and $\hat{{\eta }}$ from the test-field method. Top: $\overline{\mathrm{ℰ}}{\mathit{x}}_{}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}$, bottom: $\overline{\mathrm{ℰ}}{\mathit{y}}_{}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}$, both scaled by 10-4vA0B0. Left: directly from $\overline{{u}\mathrm{\times }{b}}$. Middle: reconstruction using all contributions k′ = 0.5,1,1.5,...16 in (20). Right: same as before, but using only the k′ = 0.5 contribution.

In the text

	Fig. 13 Quality of the mean EMF reconstruction as a function of the truncation wavenumber K′: ${\mathit{\chi }}_{{\mathit{K}}^{\mathrm{\prime }}}^{\mathrm{2}}$ (left) and correlation rK′ (right) calculated for $\overline{\mathrm{ℰ}}{\mathit{x}}_{\mathrm{R}}^{}$ (solid) and $\overline{\mathrm{ℰ}}{\mathit{y}}_{\mathrm{R}}^{}$ (dashed) using Eq. (22).
In the text

	Fig. 14 Dependence of the instability on rotational inclination θ in terms of rms value of generated field components $\mathrm{⟨}{\mathit{B}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{B}}_{\mathit{z}}^{\mathrm{2}}\mathrm{⟩}$ for the runs TF0+, TF30+, TF60+ and TF90+ in Table 3.
In the text

	Fig. 15 ${\overline{\mathit{B}}}_{\mathit{x,y}}\mathrm{\left(}\mathit{z,\lambda }\mathrm{\right)}$, $\overline{\mathrm{ℰ}}{\mathit{x,y}}_{}\mathrm{\left(}\mathit{z,\lambda }\mathrm{\right)}$ averaged between t = tsat and tsat + tA0 from the runs in Table 3 with Ω ≠ 0. (Initial velocity identical for all latitudes.) Dashed line: initial position of the magnetic layer.
In the text

Fig. 16 Same as Fig. 15, but with symmetrized ${\overline{\mathit{B}}}_{\mathit{y}}\mathrm{\left(}\mathit{z,\lambda }\mathrm{\right)}$, $\overline{\mathrm{ℰ}}{\mathit{x}}_{}\mathrm{\left(}\mathit{z,\lambda }\mathrm{\right)}$ and anti-symmetrized ${\overline{\mathit{B}}}_{\mathit{x}}\mathrm{\left(}\mathit{z,\lambda }\mathrm{\right)}$, $\overline{\mathrm{ℰ}}{\mathit{y}}_{}\mathrm{\left(}\mathit{z,\lambda }\mathrm{\right)}$. For explanation see text.

In the text

Fig. 17 $\hat{{\alpha }}\mathrm{\left(}\mathit{z,\lambda }\mathrm{\right)}$ and $\hat{{\eta }}\mathrm{\left(}\mathit{z,\lambda }\mathrm{\right)}$ from the runs in Table 3 with Ω ≠ 0, calculated using test fields with k′ = 0.5 and averaged between t = tsat and tsat + tA0. symmetrized, anti-symmetrized analogously to $\overline{{B}}$, $\overline{\mathbf{ℰ}}$ in Fig. 16. $\hat{{\alpha }}$ scaled by Urms for the run TF30+, $\hat{{\eta }}$ scaled by the molecular diffusivity η. Dashed line: initial position of the magnetic layer.

In the text

	Fig. 18 Pumping velocity scaled with Urms for the run TF30+, as a function of height and latitude for k′ = 0.5 (left) and k′ = 4.5 (right).
In the text

	Fig. 19 Vertical averages of (solid) and (dashed) for k′ = 0.5 at t = tsat and maximum of with respect to all z and t (dotted) as functions of $\hspace{0.17em}\begin{array}{}{}^{\mathrm{􏽥}}\\ \mathit{\beta }\end{array}\begin{array}{}\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}\\ \mathrm{0}\end{array}\mathrm{\propto }{\mathit{B}}_{\mathrm{0}}$.
In the text