© ESO, 2014

1. Introduction

In a turbulent medium, the turbulent pressure can lead to dynamically important effects. In particular, a stratified layer can attain a density distribution that is significantly altered compared to the nonturbulent case. In addition, magnetic fields can change the situation further, because it can locally suppress the turbulence and thus reduce the total turbulent pressure (the sum of hydrodynamic and magnetic turbulent contributions). On length scales encompassing many turbulent eddies, this total turbulent pressure reduction must be compensated for by additional gas pressure, which can lead to a density enhancement and thus to horizontal magnetic structures that become heavier than the surroundings and sink (Brandenburg et al. 2011). This is quite the contrary of magnetic buoyancy, which is still expected to operate on the smaller scale of magnetic flux tubes and in the absence of turbulence. Both effects can lead to instability: the latter is the magnetic buoyancy or interchange instability (Newcomb 1961; Parker 1966), and the former is now often referred to as negative effective magnetic pressure instability (NEMPI), which has been studied at the level of mean-field theory for the past two decades (Kleeorin et al. 1989, 1990, 1993, 1996; Kleeorin & Rogachevskii 1994; Rogachevskii & Kleeorin 2007). These are instabilities of a stratified continuous magnetic field, while the usual magnetic buoyancy instability requires nonuniform and initially separated horizontal magnetic flux tubes (Parker 1955; Spruit 1981; Schüssler et al. 1994).

Unlike the magnetic buoyancy instability, NEMPI occurs at the expense of turbulent energy rather than the energy of the gravitational field. The latter is the energy source of the magnetic buoyancy or interchange instability. NEMPI is caused by a negative turbulent contribution to the effective mean magnetic pressure (the sum of nonturbulent and turbulent contributions). For large Reynolds numbers, this turbulent contribution to the effective magnetic pressure is larger than the nonturbulent one. This results in the excitation of NEMPI and the formation of large-scale magnetic structures – even from an originally uniform mean magnetic field.

Direct numerical simulations (DNS) have recently demonstrated the operation of NEMPI in isothermally stratified layers (Brandenburg et al. 2011; Kemel et al. 2012b). This is a particularly simple case in that the density scale height is constant; i.e., the computational burden of covering large density variation is distributed over the depth of the entire layer. In spite of this simplification, it has been argued that NEMPI is important for explaining prominent features in the manifestation of solar surface activity. In particular, it has been associated with the formation of active regions (Kemel et al. 2013; Warnecke et al. 2013) and sunspots (Brandenburg et al. 2013, 2014). However, it is now important to examine the validity of conclusions based on such simplifications using more realistic models. In this paper, we therefore now consider a polytropic stratification, for which the density scale height is smallest in the upper layers, and the density variation therefore greatest.

NEMPI is a large-scale instability that can be excited in stratified small-scale turbulence. This requires (i) sufficient scale separation in the sense that the maximum scale of turbulent motions, ℓ, must be much smaller than the scale of the system, L; and (ii) strong density stratification such that the density scale height H_ρ is much smaller than L; i.e., $\begin{array}{ccc}\mathit{\ell }\mathrm{\ll }{\mathit{H}}_{\mathit{\rho }}\mathrm{\ll }\mathit{L}\mathit{.}& & \end{array}$(1)However, both the size of turbulent motions and the typical size of perturbations due to NEMPI can be related to the density scale height. Furthermore, earlier work of Kemel et al. (2013) using isothermal layers shows that the scale of perturbations due to NEMPI exceeds the typical density scale height. Unlike the isothermal case, in which the scale height is constant, it decreases rapidly with height in a polytropic layer. It is then unclear how such structures could fit into the narrow space left by the stratification and whether the scalings derived for the isothermal case can still be applied locally to polytropic layers.

NEMPI has already been studied previously for polytropic layers in mean-field simulations (MFS; Brandenburg et al. 2010; Käpylä et al. 2012; Jabbari et al. 2013), but no systematic comparison has been made with NEMPI in isothermal or in polytropic layers with different values of the polytropic index. This will be done in the present paper, both in MFS and DNS. Those two complementary types of simulations have proved to be a good tool for understanding the underlying physics of NEMPI. An example are the studies of the effects of rotation on NEMPI (Losada et al. 2012, 2013), where MFS have been able to give quantitatively useful predictions before corresponding DNS were able to confirm the resulting dependence.

2. Polytropic stratification

We discuss here the equation for the vertical profile of the fluid density in a polytropic layer. In a Cartesian plane-parallel layer with polytropic stratification, the temperature gradient is constant, so the temperature goes linearly to zero at z_∞. The temperature, T, is proportional to the square of the sound speed, ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$, and thus also to the density scale height H_ρ(z), which is given by ${\mathit{H}}_{\mathit{\rho }}\mathrm{=}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{/}\mathit{g}$ for an isentropic stratification, where g is the acceleration due to the gravity. For a perfect gas, the density ρ is proportional to Tⁿ, and the pressure p is proportional to Tⁿ⁺¹, such that p/ρ is proportional to T, where n is the polytropic index. Furthermore, we have p(z) ∝ ρ(z)^Γ, where Γ = (n + 1)/n is another useful coefficient.

For a perfect gas, the specific entropy can be defined (up to an additive constant) as s = c_vln(p/ρ^γ), where γ = c_p/c_v is the ratio of specific heats at constant pressure and constant density, respectively. For a polytropic stratification, we have $\exp \mathrm{(}\mathit{s}\mathit{/}{\mathit{c}}_{\mathit{v}}\mathrm{)}\mathrm{=}\mathit{p}\mathit{/}{\mathit{\rho }}^{\mathit{\gamma }}\mathrm{\propto }{\mathit{\rho }}^{\mathrm{\Gamma }\mathrm{-}\mathit{\gamma }}\mathit{,}$(2)so s is constant when Γ = γ, which is the case for an isentropic stratification. In the following, we make this assumption and specify from now only the value of γ. For a monatomic gas, we have γ = 5/3, which is relevant for the Sun, while for a diatomic molecular gas, we have γ = 7/5, which is relevant for air. In those cases, a stratification with Γ = γ can be motivated by assuming perfect mixing across the layer. The isothermal case with γ = 1 can be motivated by assuming rapid heating/cooling to a constant temperature.

Fig. 1 Isothermal and polytropic relations for different values of γ when calculated using the conventional formula ρ ∝ (z∞ − z)n with n = 1/(γ − 1).

Our aim is to study the change in the properties of NEMPI in a continuous fashion as we go from an isothermally stratified layer to a polytropic one. In the latter case, the fluid density varies in a power law fashion, ρ ∝ (z_∞ − z)ⁿ, while in the former, it varies exponentially, ρ ∝ exp(z_∞ − z). This is shown in Fig. 1 where we compare the exponential isothermal atmosphere with a family of polytropic atmospheres with γ = 1.2, 1.4, and 5/3. Clearly, there is no continuous connection between the isothermal case and the polytropic one in the limit γ → 1. This cannot be fixed by rescaling the isothermal density stratification, because in Fig. 1 its values would still lie closer to 5/3 than to 1.4 or 1.2. Another problem with this description is that for polytropic solutions the density is always zero at z = z_∞, but finite in the isothermal case. These different behaviors between isothermal and polytropic atmospheres can be unified by using the generalized exponential function known as the “q-exponential” (see, e.g., Yamano 2002), which is defined as ${\mathit{e}}_{\mathit{q}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}{\left[\mathrm{1}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{q}\mathrm{)}\mathit{x}\right]}^{\mathrm{1}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{q}\mathrm{)}}\mathit{,}$(3)where the parameter q is related to γ via q = 2 − γ. This generalization of the usual exponential function was originally introduced by Tsallis (1988) in connection with a possible generalization of the Boltzmann-Gibbs statistics. Its usefulness in connection with stellar polytropes has been employed by Plastino & Plastino (1993). Thus, the density stratification is given by $\frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{0}}}\mathrm{=}{\left[\mathrm{1}\mathrm{+}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}\left(\mathrm{-}\frac{\mathit{z}}{{\mathit{H}}_{\mathit{\rho }\mathrm{0}}}\right)\right]}^{\mathrm{1}\mathit{/}\mathrm{(}\mathit{\gamma }\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{1}\mathrm{)}}\mathrm{=}{\left(\mathrm{1}\mathrm{-}\frac{\mathit{z}}{\mathit{n}{\mathit{H}}_{\mathit{\rho }\mathrm{0}}}\right)}^{\mathit{n}}\mathit{,}$(4)which reduces to ρ/ρ₀ = exp(−z/H_ρ0) for isothermal stratification with γ → 1 and n → ∞. The density scale height is then given by ${\mathit{H}}_{\mathit{\rho }}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\mathit{H}}_{\mathit{\rho }\mathrm{0}}\mathrm{-}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}\mathit{z}\mathrm{=}{\mathit{H}}_{\mathit{\rho }\mathrm{0}}\mathrm{-}\mathit{z}\mathit{/}\mathit{n}\mathit{.}$(5)In the following, we measure lengths in units of H_ρ0 = H_ρ(0).

In Fig. 2 we show the dependencies of ρ(z) and H_ρ given by Eqs. (4) and (5) for different values of γ. Compared with Fig. 1, where z_∞ is held fixed, in Fig. 2 it is equal to z_∞ = nH_ρ0 = H_ρ0/(γ − 1). The total density contrast is roughly the same in all four cases for different γ, but for increasing values of γ, the vertical density gradient becomes progressively stronger in the upper layers.

Fig. 2 Polytropes (Eq. (4)) with γ = 1 (solid line), 1.2 (dash-dotted), 1.4 (dotted), and 5/3 (dashed) and density scale height (Eq. (5)) for −4 ≤ −z/Hρ0 ≤ 1.2. The total density contrast is similar for γ = 1 and 5/3.

Fig. 3 Sketch showing the expected size distribution of the nearly circular NEMPI eigenfunction structures at different heights.

Assuming that the radius R of the resulting structures is proportional to H_ρ, we sketch in Fig. 3 a situation in which R is half the depth. With ρ ∝ (z_∞ − z)ⁿ, the density scale height is H_ρ(z) = (z_∞ − z)/n. Thus, Fig. 3 applies to a case in which R = (z_∞ − z)/2 = (n/2) H_ρ(z). The solar convection zone is nearly isentropic and well described by n = 3/2. This means that the structures of Fig. 3 have R = (3/4) H_ρ(z). The results of Kemel et al. (2013) and Brandenburg et al. (2014) suggest that the horizontal wavenumber of structures, k_⊥, formed by NEMPI, is less than or about ${\mathit{H}}_{\mathit{\rho }}^{-1}$, so their horizontal wavelength is ≲2πH_ρ. One wavelength corresponds to the distance between two nodes, i.e., the distance between two spheres, which is 4 R. Thus, in the isothermal case we have R/H_ρ = 2π/4 ≈ 1.5, which implies that such a structure would not fit into the isentropic atmosphere described above. This provides an additional motivation for our present work.

3. DNS study

In this section we study NEMPI in DNS for the polytropic layer. Corresponding MFS are presented in Sect. 4.

3.1. The model

We solve the equations for the magnetic vector potential A, the velocity U, and the density ρ, in the form

$\begin{array}{ccc}& & \frac{\mathit{\partial }{A}}{\mathit{\partial t}}\mathrm{=}{U}\mathrm{\times }{B}\mathrm{-}\mathit{\eta }{\mathit{\mu }}_{\mathrm{0}}{J}\mathit{,}\\ & & \frac{\mathrm{D}{U}}{\mathrm{D}\mathit{t}}\mathrm{=}\frac{\mathrm{1}}{\mathit{\rho }}{J}\mathrm{\times }{B}\mathrm{+}{f}\mathrm{-}\mathit{\nu }{Q}\mathrm{-}{\nabla }\mathit{H,}\\ & & \frac{\mathrm{D}\mathit{\rho }}{\mathrm{D}\mathit{t}}\mathrm{=}\mathrm{-}\mathit{\rho }{\nabla }\mathrm{·}{U}\mathit{,}\end{array}$where D/Dt = ∂/∂t + U·∇ is the advective derivative with respect to the actual (turbulent) flow, B = B₀ + ∇ × A is the magnetic field, B₀ the imposed uniform field, J = ∇ × B/μ₀ the current density, μ₀ the vacuum permeability, H = h + Φ the reduced enthalpy, h = c_pT the enthalpy, Φ is the gravitational potential, $\begin{array}{ccc}\mathrm{-}{Q}\mathrm{=}{\mathrm{\nabla }}^{\mathrm{2}}{U}\mathrm{+}\frac{\mathrm{1}}{\mathrm{3}}{\nabla }{\nabla }\mathrm{·}{U}\mathrm{+}\mathrm{2}{S}{\nabla }\ln \mathit{\rho }& & \end{array}$(9)is a term appearing in the viscous force −νQ, S is the traceless rate-of-strain tensor with components $\begin{array}{ccc}{\mathrm{S}}_{\mathit{ij}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathrm{\nabla }}_{\mathit{j}}{\mathit{U}}_{\mathit{i}}\mathrm{+}{\mathrm{\nabla }}_{\mathit{i}}{\mathit{U}}_{\mathit{j}}\mathrm{)}\mathrm{-}\frac{\mathrm{1}}{\mathrm{3}}{\mathit{\delta }}_{\mathit{ij}}{\nabla }\mathrm{·}{U}\mathit{,}& & \end{array}$(10)ν is the kinematic viscosity, and η is the magnetic diffusion coefficient caused by electrical conductivity of the fluid. As in Losada et al. (2012), z corresponds to radius, x to colatitude, and y to azimuth. The forcing function f consists of random, white-in-time, plane, nonpolarized waves with a certain average wavenumber k_f.

3.2. Boundary conditions and parameters

In the DNS we use stress-free boundary conditions for the velocity at the top and bottom; i.e., ∇_zU_x = ∇_zU_y = U_z = 0. For the magnetic field we use either perfect conductor boundary conditions, A_x = A_y = ∇_zA_z = 0, or vertical field conditions, ∇_zA_x = ∇_zA_y = A_z = 0, again at both the top and bottom. All variables are assumed periodic in the x and y directions.

Fig. 4 Snapshots of ${\overline{\mathit{B}}}_{\mathit{y}}$ from DNS for γ = 5/3 and β0 = 0.02 (upper row), 0.05 (middle row), and 0.1 (lower row) at different times (indicated in turbulent-diffusive times, increasing from left to right) in the presence of a horizontal field using the perfect conductor boundary condition.

The turbulent rms velocity is approximately independent of z with ${\mathit{u}}_{\mathrm{rms}}\mathrm{=}\mathrm{⟨}{u}{\mathrm{2}}^{}{\mathrm{⟩}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\approx }\mathrm{0.1}\hspace{0.17em}{\mathit{c}}_{\mathrm{s}}$. The gravitational acceleration g = (0,0, −g) is chosen such that k₁H_ρ0 = 1, where k₁ = 2π/L and L is the size of the domain. With one exception (Sect. 3.5), we always use the value k_f/k₁ = 30 for the scale separation ratio. For B₀ we choose either a horizontal field pointing in the y direction or a vertical one pointing in the z direction. The latter case, B₀ = (0,0,B₀), is usually combined with the use of the vertical field boundary condition, while the former one, B₀ = (0,B₀,0), is combined with using perfect conductor boundary conditions. The strength of the imposed field is expressed in terms of B_eq0 = B_eq(z = 0), which is the equipartition field strength at z = 0. Here, the equipartition field ${\mathit{B}}_{\mathrm{eq}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\left({\mathit{\mu }}_{\mathrm{0}}\overline{\mathit{\rho }}\mathrm{\left(}\mathit{z}\mathrm{\right)}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{u}}_{\mathrm{rms}}$. The imposed field is normalized by B_eq0 and denoted as β₀ = B₀/B_eq0, while $\mathit{\beta }\mathrm{=}\mathrm{\left|}\overline{{B}}\mathrm{\right|}\mathit{/}{\mathit{B}}_{\mathrm{eq}}$ is the modulus of the normalized mean magnetic field. Time is expressed in terms of the turbulent-diffusive time, ${\mathit{\tau }}_{\mathrm{td}}\mathrm{=}{\mathit{H}}_{\mathit{\rho }\mathrm{0}}^{\mathrm{2}}\mathit{/}{\mathit{\eta }}_{\mathrm{t}\mathrm{0}}$, where η_t0 = u_rms/3k_f (Sur et al. 2008) is an estimate for the turbulent magnetic diffusivity used in the DNS.

Our values of ν and η are characterized by specifying the kinetic and magnetic Reynolds numbers, $\begin{array}{ccc}\mathrm{Re}\mathrm{=}{\mathit{u}}_{\mathrm{rms}}\mathit{/}\mathit{\nu }{\mathit{k}}_{\mathrm{f}}\mathit{,} {\mathit{R}}_{\mathrm{m}}\mathrm{=}{\mathit{u}}_{\mathrm{rms}}\mathit{/}\mathit{\eta }{\mathit{k}}_{\mathrm{f}}\mathit{.}& & \end{array}$(11)In most of this paper (except in Sect. 3.5) we use Re = 36 and R_m = 18, which are also the values used by Kemel et al. (2013).

The DNS are performed with the Pencil Code, http://pencil-code.googlecode.com, which uses sixth-order explicit finite differences in space and a third-order accurate time stepping method. We use a numerical resolution of 256³ mesh points.

3.3. Horizontal fields

We focus on the case γ = 5/3 and show in Fig. 4 visualizations of ${\overline{\mathit{B}}}_{\mathit{y}}$ at different instants for three values of the imposed horizontal magnetic field strength. For β₀ = 0.02, a magnetic structure is clearly visible at t/τ_td = 1.43, while for β₀ = 0.05 structures are already fully developed at t/τ_td = 0.42. In that case (β₀ = 0.05), at early times (t/τ_td = 0.42), there are two structures, which then begin to merge at t/τ_td = 1.38. The growth rate of the magnetic structure is found to be $\mathit{\lambda }\mathrm{\approx }\mathrm{2}{\mathit{\eta }}_{\mathrm{t}\mathrm{0}}\mathit{/}{\mathit{H}}_{\mathit{\rho }\mathrm{0}}^{\mathrm{2}}$ for the runs shown in Fig. 4. This is less than the value of $\mathit{\lambda }\mathrm{\approx }\mathrm{5}{\mathit{\eta }}_{\mathrm{t}\mathrm{0}}\mathit{/}{\mathit{H}}_{\mathit{\rho }\mathrm{0}}^{\mathrm{2}}$ found earlier for the isothermal case (Kemel et al. 2012b).

Fig. 5 Snapshots from DNS showing ${\overline{\mathit{B}}}_{\mathit{z}}$ on the periphery of the computational domain for γ = 5/3 and β0 = 0.05 at different times for the case of a vertical field using the vertical field boundary condition.

Fig. 6 Cuts of Bz/Beq(z) in (a) the xy plane at the top boundary (z/Hρ0 = 1.2) and (b) the xz plane through the middle of the spot at y = 0 for γ = 5/3 and β0 = 0.05. In the xz cut, we also show magnetic field lines and flow vectors obtained by numerically averaging in azimuth around the spot axis.

For γ = 5/3 and β₀ ≤ 0.02, the magnetic structures become smaller (k_⊥H_ρ0 = 2) near the surface. In the nonlinear regime, i.e., at late times, the structures move downward due to the so-called “potato-sack” effect, which was first seen in MFS (Brandenburg et al. 2010) and later confirmed in DNS (Brandenburg et al. 2011). The magnetic structures sink in the nonlinear stage of NEMPI, because an increase in the mean magnetic field inside the horizontal magnetic flux tube increases the absolute value of the effective magnetic pressure. On the other hand, a decrease in the negative effective magnetic pressure is balanced out by increased gas pressure, which in turn leads to higher density, so the magnetic structures become heavier than the surroundings and sink. This potato-sack effect has been clearly observed in the present DNS with the polytropic layer (see the right column in Fig. 4).

3.4. Vertical fields

Recent DNS using isothermal layers have shown that strong circular flux concentrations can be produced in the case of a vertical magnetic field (Brandenburg et al. 2013, 2014). This is also observed in the present study of a polytropic layer; see Fig. 5, where we show the evolution of ${\overline{\mathit{B}}}_{\mathit{z}}$ on the periphery of the computational domain for γ = 5/3 and β₀ = 0.05 at different times. A difference to the DNS for γ = 1 (Brandenburg et al. 2013) seems to be that for γ = 5/3 the magnetic structures are shallower than for γ = 1; see Fig. 6, where we show xy and xz slices of B_z through the spot. Owing to periodicity in the xy plane, we have shifted the location of the spot to x = y = 0. We note also that the field lines of the averaged magnetic field show a structure rather similar to the one found in MFS of Brandenburg et al. (2014). The origin of circular structures is associated with a cylindrical symmetry for the vertical magnetic field. The growth rate of the magnetic field in the spot is found to be $\mathit{\lambda }\mathrm{\approx }\mathrm{0.9}{\mathit{\eta }}_{\mathrm{t}\mathrm{0}}\mathit{/}{\mathit{H}}_{\mathit{\rho }\mathrm{0}}^{\mathrm{2}}$, which is similar to the value of 1.3 found earlier for the isothermal case (Brandenburg et al. 2013).

3.5. Effective magnetic pressure

As pointed out in Sect. 1, the main reason for the formation of strongly inhomogeneous large-scale magnetic structures is the negative contribution of turbulence to the large-scale magnetic pressure, so that the effective magnetic pressure (the sum of turbulent and nonturbulent contributions) can be negative at high magnetic Reynolds numbers. The effective magnetic pressure has been determined from DNS for isothermally stratified forced turbulence (Brandenburg et al. 2010, 2012) and for turbulent convection (Käpylä et al. 2012). To see whether the nature of polytropic stratification has any influence on the effective magnetic pressure, we use DNS.

We first explain the essence of the effect of turbulence on the effective magnetic pressure. We consider the momentum equation in the form $\begin{array}{ccc}\frac{\mathit{\partial }}{\mathit{\partial t}}\mathit{\rho }\hspace{0.17em}{\mathit{U}}_{\mathit{i}}\mathrm{=}\mathrm{-}\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{x}}_{\mathit{j}}}{\mathrm{\Pi }}_{\mathit{ij}}\mathrm{+}\mathit{\rho }\hspace{0.17em}{\mathit{g}}_{\mathit{i}}\mathit{,}& & \end{array}$(12)where $\begin{array}{ccc}{\mathrm{\Pi }}_{\mathit{ij}}\mathrm{=}\mathit{\rho }\hspace{0.17em}{\mathit{U}}_{\mathit{i}}{\mathit{U}}_{\mathit{j}}\mathrm{+}{\mathit{\delta }}_{\mathit{ij}}\mathrm{(}\mathit{p}\mathrm{+}{{B}}^{\mathrm{2}}\mathit{/}\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}\mathrm{)}\mathrm{-}{\mathit{B}}_{\mathit{i}}{\mathit{B}}_{\mathit{j}}\mathit{/}{\mathit{\mu }}_{\mathrm{0}}\mathrm{-}\mathrm{2}\mathit{\nu \rho }\hspace{0.17em}{\mathrm{S}}_{\mathit{ij}}& & \end{array}$(13)is the momentum stress tensor and δ_ij the Kronecker tensor.

Ignoring correlations between velocity and density fluctuations for low-Mach number turbulence, the averaged momentum equation is $\begin{array}{ccc}\frac{\mathit{\partial }}{\mathit{\partial t}}\overline{\mathit{\rho }}\hspace{0.17em}\overline{{U}}{}_{\mathit{i}}\mathrm{=}\mathrm{-}\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{x}}_{\mathit{j}}}{\overline{\mathrm{\Pi }}}_{\mathit{ij}}\mathrm{+}\overline{\mathit{\rho }}\hspace{0.17em}{\mathit{g}}_{\mathit{i}}\mathit{,}& & \end{array}$(14)where an overbar means xy averaging, ${\overline{\mathrm{\Pi }}}_{\mathit{ij}}\mathrm{=}\overline{\mathrm{\Pi }}\begin{array}{c}\mathrm{m}\\ \mathit{ij}\end{array}\mathrm{+}\overline{\mathrm{\Pi }}\begin{array}{c}\mathrm{f}\\ \mathit{ij}\end{array}$ is the mean momentum stress tensor, split into contributions resulting from the mean field (indicated by superscript m) and the fluctuating field (indicated by superscript f). The tensor $\overline{\mathrm{\Pi }}\begin{array}{c}\mathrm{m}\\ \mathit{ij}\end{array}$ has the same form as Eq. (13), but all quantities have now attained an overbar; i.e., $\begin{array}{ccc}\overline{\mathrm{\Pi }}\begin{array}{c}\mathrm{m}\\ \mathit{ij}\end{array}\mathrm{=}\overline{\mathit{\rho }}\hspace{0.17em}{\overline{\mathit{U}}}_{\mathit{i}}{\overline{\mathit{U}}}_{\mathit{j}}\mathrm{+}{\mathit{\delta }}_{\mathit{ij}}\mathrm{(}\overline{\mathit{p}}\mathrm{+}\overline{{B}}{}^{\mathrm{2}}\mathit{/}\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}\mathrm{)}\mathrm{-}{\overline{\mathit{B}}}_{\mathit{i}}{\overline{\mathit{B}}}_{\mathit{j}}\mathit{/}{\mathit{\mu }}_{\mathrm{0}}\mathrm{-}\mathrm{2}\mathit{\nu }\overline{\mathit{\rho }}\hspace{0.17em}{\overline{\mathrm{S}}}_{\mathit{ij}}\mathit{.}& & \end{array}$(15)The contributions, $\overline{\mathrm{\Pi }}\begin{array}{c}\mathrm{f}\\ \mathit{ij}\end{array}$, which result from the fluctuations ${u}\mathrm{=}{U}\mathrm{-}\overline{{U}}$ and ${b}\mathrm{=}{B}\mathrm{-}\overline{{B}}$ of velocity and magnetic fields, respectively, are determined by the sum of the Reynolds and Maxwell stress tensors: $\begin{array}{ccc}\overline{\mathrm{\Pi }}\begin{array}{c}\mathrm{f}\\ \mathit{ij}\end{array}\mathrm{=}\overline{\mathit{\rho }}\hspace{0.17em}\overline{{\mathit{u}}_{\mathit{i}}{\mathit{u}}_{\mathit{j}}}\mathrm{+}{\mathit{\delta }}_{\mathit{ij}}\overline{{{b}}^{\mathrm{2}}}\mathit{/}\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}\mathrm{-}\overline{{\mathit{b}}_{\mathit{i}}{\mathit{b}}_{\mathit{j}}}\mathit{/}{\mathit{\mu }}_{\mathrm{0}}\mathit{.}& & \end{array}$(16)This contribution, together with the contribution from the mean field, $\overline{\mathrm{\Pi }}\begin{array}{c}\mathrm{m}\\ \mathit{ij}\end{array}$, comprises the total mean momentum tensor. The contribution from the fluctuating fields is split into a contribution that is independent of the mean magnetic field $\overline{\mathrm{\Pi }}\begin{array}{c}\mathrm{f}\mathit{,}\mathrm{0}\\ \mathit{ij}\end{array}$ (which determines the turbulent viscosity and background turbulent pressure) and a contribution that depends on the mean magnetic field $\overline{\mathrm{\Pi }}\begin{array}{c}\mathrm{f}\mathit{,}\overline{\mathit{B}}\\ \mathit{ij}\end{array}$. The difference between the two, $\mathrm{\Delta }\overline{\mathrm{\Pi }}\begin{array}{c}\mathrm{f}\\ \mathit{ij}\end{array}\mathrm{=}\overline{\mathrm{\Pi }}\begin{array}{c}\mathrm{f}\mathit{,}\overline{\mathit{B}}\\ \mathit{ij}\end{array}\mathrm{-}\overline{\mathrm{\Pi }}\begin{array}{c}\mathrm{f}\mathit{,}\mathrm{0}\\ \mathit{ij}\end{array}$, is caused by the mean magnetic field and is parameterized in the form $\begin{array}{ccc}\mathrm{\Delta }\overline{\mathrm{\Pi }}\begin{array}{c}\mathrm{f}\\ \mathit{ij}\end{array}\mathrm{=}{\mathit{\mu }}_{\mathrm{0}}^{-1}\mathrm{(}{\mathit{q}}_{\mathrm{s}}{\overline{\mathit{B}}}_{\mathit{i}}{\overline{\mathit{B}}}_{\mathit{j}}\mathrm{-}{\mathit{q}}_{\mathrm{p}}\hspace{0.17em}{\mathit{\delta }}_{\mathit{ij}}\overline{{B}}{}^{\mathrm{2}}\mathit{/}\mathrm{2}\mathrm{-}{\mathit{q}}_{\mathrm{g}}\hspace{0.17em}\mathit{ĝ}\mathit{i}\hspace{0.17em}\mathit{ĝ}\mathit{j}\overline{{B}}{}^{\mathrm{2}}\mathrm{)}\mathit{,}& & \end{array}$(17)where the coefficient q_p represents the isotropic turbulence contribution to the mean magnetic pressure, the coefficient q_s represents the turbulence contribution to the mean magnetic tension, while the coefficient q_g is the anisotropic turbulence contribution to the mean magnetic pressure, and it characterizes the effect of vertical variations of the magnetic field caused by the vertical magnetic pressure gradient. Here, ĝ_i is the unit vector in the direction of the gravity field (in the vertical direction). We consider cases with horizontal and vertical mean magnetic fields separately. Analytically, the coefficients q_p, q_g, and q_s have been obtained using both the spectral τ approach (Rogachevskii & Kleeorin 2007) and the renormalization approach (Kleeorin & Rogachevskii 1994). The form of Eq. (17) is also obtained using simple symmetry arguments; e.g., for a horizontal field, the linear combination of three independent true tensors, δ_ij,ĝ_iĝ_j and ${\overline{\mathit{B}}}_{\mathit{i}}{\overline{\mathit{B}}}_{\mathit{j}}$, yields Eq. (17), while for the vertical field, the linear combination of only two independent true tensors, δ_ij and ${\overline{\mathit{B}}}_{\mathit{i}}{\overline{\mathit{B}}}_{\mathit{j}}$.

Previous DNS studies (Brandenburg et al. 2012) have shown that q_s and q_g are negligible for forced turbulence. To avoid the formation of magnetic structures in the nonlinear stage of NEMPI, which would modify our results, we use here a lower scale separation ratio, k_f/k₁ = 5, keeping k₁H_ρ = 1, and using Re = 140 and R_m = 70, as in Brandenburg et al. (2012). To determine q_p(β), it is sufficient to measure the three diagonal components of $\overline{\mathrm{\Pi }}\begin{array}{c}\mathrm{f}\\ \mathit{ij}\end{array}$ both with and without an imposed magnetic field using ${\mathit{q}}_{\mathrm{p}}\mathrm{=}\mathrm{-}\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}\mathrm{\Delta }\overline{\mathrm{\Pi }}\begin{array}{c}\mathrm{f}\\ \mathit{xx}\end{array}\mathit{/}{\overline{{B}}}^{\mathrm{2}}$.

Fig. 7 Effective magnetic pressure obtained from DNS in a polytropic layer with different γ for horizontal (H, red curves) and vertical (V, blue curves) mean magnetic fields.

In Fig. 7 we present the results for forced turbulence in the polytropic layer with different γ for horizontal and vertical mean magnetic fields. It turns out that the normalized effective magnetic pressure, $\begin{array}{ccc}{\mathrm{𝒫}}_{\mathrm{eff}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{q}}_{\mathrm{p}}\mathrm{)}{\mathit{\beta }}^{\mathrm{2}}\mathit{,}& & \end{array}$(18)has a minimum value at β_min. Following Kemel et al. (2012a), the function q_p(β) is approximated by: $\begin{array}{ccc}{\mathit{q}}_{\mathrm{p}}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathrm{=}\frac{{\mathit{q}}_{\mathrm{p}\mathrm{0}}}{\mathrm{1}\mathrm{+}{\mathit{\beta }}^{\mathrm{2}}\mathit{/}{\mathit{\beta }}_{\mathrm{p}}^{\mathrm{2}}}\mathrm{=}\frac{{\mathit{\beta }}_{\mathit{\star }}^{\mathrm{2}}}{{\mathit{\beta }}_{\mathrm{p}}^{\mathrm{2}}\mathrm{+}{\mathit{\beta }}^{\mathrm{2}}}\mathit{,}& & \end{array}$(19)where q_p0, β_p, and ${\mathit{\beta }}_{\mathit{\star }}\mathrm{=}{\mathit{\beta }}_{\mathrm{p}}{\mathit{q}}_{\mathrm{p}\mathrm{0}}^{\mathrm{1}\mathit{/}\mathrm{2}}$ are constants. This equation can be understood as a quenching formula for the effective magnetic pressure; see Jabbari et al. (2013). The coefficients β_p and β_⋆ are related to at β_min via (Kemel et al. 2012a) ${\mathit{\beta }}_{\mathrm{p}}\mathrm{=}{\mathit{\beta }}_{\min }^{\mathrm{2}}􏼮\sqrt{\mathrm{-}\mathrm{2}{\mathrm{𝒫}}_{\min }}\mathit{,}{\mathit{\beta }}_{\mathit{\star }}\mathrm{=}{\mathit{\beta }}_{\mathrm{p}}\mathrm{+}\sqrt{\mathrm{-}\mathrm{2}{\mathrm{𝒫}}_{\min }}\mathit{.}$(20)In Fig. 8 we show these fitting parameters for the function q_p(β) for polytropic layer with different γ for horizontal and vertical mean magnetic fields. The effects of negative effective magnetic pressure are generally stronger for vertical magnetic fields (q_p0 is larger and β_p smaller) than for horizontal ones (q_p0 is smaller and β_p larger), but the values ${\mathit{\beta }}_{\mathit{\star }}\mathrm{=}{\mathit{\beta }}_{\mathrm{p}}{\mathit{q}}_{\mathrm{p}\mathrm{0}}^{\mathrm{1}\mathit{/}\mathrm{2}}$ are similar in both cases, and increasing from 0.35 (for γ = 1) to 0.6 (for γ = 5/3); see Fig. 8.

Fig. 8 Parameters qp0, βp, and β⋆ for the function qp(β) (see Eq. (19)) versus γ for horizontal (red line) and vertical (blue line) mean magnetic fields.

4. Mean-field study

We now consider two sets of parameters that we refer to as Model I (with q_p0 = 32 and β_p = 0.058 corresponding to β_⋆ = 0.33) and Model II (q_p0 = 9 and β_p = 0.21 corresponding to β_⋆ = 0.63). These cases are representative of the strong (large β_⋆) and weak (small β_⋆) effects of NEMPI. Following earlier studies (Brandenburg et al. 2012), we find q_s to be compatible with zero. We thus neglect this coefficient in the following.

4.1. Governing parameters and estimates

The purpose of this section is to summarize the findings for the isothermal case in MFS. One of the key results is the prediction of the growth rate of NEMPI. The work of Kemel et al. (2013) showed that in the ideal case (no turbulent diffusion), the growth rate λ is approximated well by $\mathit{\lambda }\mathrm{\approx }{\mathit{\beta }}_{\mathit{\star }}{\mathit{u}}_{\mathrm{rms}}\mathit{/}{\mathit{H}}_{\mathit{\rho }} \mathrm{\left(}\mathrm{noturbulentdiffusion}\mathrm{\right)}\mathit{.}$(21)However, turbulent magnetic diffusion, η_t, can clearly not be neglected and is chiefly responsible for shutting off NEMPI if the turbulent eddies are too big and η_t too large. This was demonstrated in Fig. 17 of Brandenburg et al. (2012). A heuristic ansatz, which is motivated by similar circumstances in mean-field dynamo theory (Krause & Rädler 1980), is to add a term −η_tk² to the righthand side of Eq. (21), where k is the wavenumber of NEMPI.

To specify the expression for λ, we normalize the wavenumber of the perturbations by the inverse density scale height and denote this by κ ≡ kH_ρ. The wavenumber of the energy-carrying turbulent eddies k_f is in nondimensional form κ_f ≡ k_fH_ρ, and the normalized horizontal wavenumber of the resulting magnetic structures is referred to as κ_⊥ = k_⊥H_ρ. For NEMPI, these values have been estimated to be κ_⊥ = 0.8–1.0, and can even be smaller for vertical magnetic fields (Brandenburg et al. 2014). Using an approximate aspect ratio of unity for magnetic structures, we have $\mathit{\kappa }\mathrm{=}\sqrt{\mathrm{2}}{\mathit{\kappa }}_{\mathrm{\perp }}\mathrm{\approx }\mathrm{1.1}$–1.4. In stellar mixing length theory (Vitense 1953), the mixing length is ℓ_f = α_mixH_ρ, where α_mix ≈ 1.6 is a nondimensional mixing length parameter. Since k_f = 2π/ℓ_f, we arrive at the following estimate: κ_f = 2πγ/α_mix ≈ 6.5. [Owing to a confusion between pressure and density scale heights, this value was underestimated by Kemel et al. (2013) to be 2.4, although an independent calculation of this value from turbulent convection simulations would still be useful.] Using η_t ≈ u_rms/3k_f, the turbulent magnetic diffusive rate for an isothermal atmosphere is given by ${\mathit{\eta }}_{\mathrm{t}}{\mathit{k}}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{u}}_{\mathrm{rms}}}{\mathrm{3}{\mathit{H}}_{\mathit{\rho }}}\frac{{\mathit{\kappa }}^{\mathrm{2}}}{{\mathit{\kappa }}_{\mathrm{f}}}\mathit{,}$(22)and the growth rate of NEMPI in that normalization is $\frac{\mathit{\lambda }}{{\mathit{\eta }}_{\mathrm{t}}{\mathit{k}}^{\mathrm{2}}}\mathrm{=}\mathrm{3}{\mathit{\beta }}_{\mathit{\star }}\frac{{\mathit{\kappa }}_{\mathrm{f}}}{{\mathit{\kappa }}^{\mathrm{2}}}\mathrm{-}\mathrm{1.}$(23)Using β_⋆ = 0.23, which is the relevant value for high magnetic Reynolds numbers (Brandenburg et al. 2012), we find λ/η_tk² ≈ 2.7−1.3 for κ ≈ 1.1−1.4. However, since Fig. 8 shows an increase of β_⋆ with increasing polytropic index, one might expect a corresponding increase in the growth rate of NEMPI for a polytropic layer, in which H_ρ varies strongly with height z. Indeed, in a polytropic atmosphere, H_ρ is proportional to depth. Thus, at any given depth there is a layer beneath, where the stratification is less strong and the growth rate of NEMPI is lower. In addition, there is a thinner, more strongly stratified layer above, where NEMPI might grow faster if only the structures generated by NEMPI have enough room to develop before they touch the top of the atmosphere at z_∞, where the temperature vanishes.

4.2. Mean-field equations

In the following, we consider MFS and compare it with DNS. We also compare our MFS results with those of DNS (Brandenburg et al. 2011; Kemel et al. 2013) using a similar polytropic setup.

The governing equations for the mean quantities (denoted by an overbar) are fairly similar to those for the original equations, except that in the MFS viscosity and magnetic diffusivity are replaced by their turbulent counterparts, and the mean Lorentz force is supplemented by a parameterization of the turbulent contribution to the effective magnetic pressure.

The evolution equations for mean vector potential $\overline{{A}}$, mean velocity $\overline{{U}}$, and mean density $\overline{\mathit{\rho }}$, are

$\begin{array}{ccc}& & \frac{\mathit{\partial }\overline{{A}}}{\mathit{\partial t}}\mathrm{=}\overline{{U}}\mathrm{\times }\overline{{B}}\mathrm{-}{\mathit{\eta }}_{\mathrm{T}}{\mathit{\mu }}_{\mathrm{0}}\overline{{J}}\mathit{,}\\ & & \frac{\overline{\mathrm{D}}\hspace{0.17em}\overline{{U}}}{\overline{\mathrm{D}}\mathit{t}}\mathrm{=}\frac{\mathrm{1}}{\overline{\mathit{\rho }}}\mathrm{[}\overline{{J}}\mathrm{\times }\overline{{B}}\mathrm{+}{\nabla }\mathrm{(}{\mathit{q}}_{\mathrm{p}}\overline{{B}}{}^{\mathrm{2}}\mathit{/}\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{]}}\mathrm{-}{\mathit{\nu }}_{\mathrm{T}}\overline{{Q}}\mathrm{-}{\nabla }\overline{\mathit{H}}\mathit{,}\\ & & \frac{\overline{\mathrm{D}}\hspace{0.17em}\overline{\mathit{\rho }}}{\overline{\mathrm{D}}\mathit{t}}\mathrm{=}\mathrm{-}\overline{\mathit{\rho }}{\nabla }\mathrm{·}\overline{{U}}\mathit{,}\end{array}$Awhere $\overline{\mathrm{D}}\mathit{/}\overline{\mathrm{D}}\mathit{t}\mathrm{=}\mathit{\partial }\mathit{/}\mathit{\partial t}\mathrm{+}\overline{{U}}\mathrm{·}{\nabla }$ is the advective derivative with respect to the mean flow, $\overline{\mathit{\rho }}$ the mean density, $\overline{\mathit{H}}\mathrm{=}\overline{\mathit{h}}\mathrm{+}\mathrm{\Phi }$ the mean reduced enthalpy, $\overline{\mathit{h}}\mathrm{=}{\mathit{c}}_{\mathit{p}}\overline{\mathit{T}}$ the mean enthalpy, $\overline{\mathit{T}}\mathrm{\propto }{\overline{\mathit{\rho }}}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}$ the mean temperature, Φ the gravitational potential, η_T = η_t + η, and ν_T = ν_t + ν are the sums of turbulent and microphysical values of magnetic diffusivity and kinematic viscosities, respectively. Also, $\overline{{J}}\mathrm{=}{\nabla }\mathrm{\times }\overline{{B}}\mathit{/}{\mathit{\mu }}_{\mathrm{0}}$ is the mean current density, μ₀ is the vacuum permeability, $\begin{array}{ccc}\mathrm{-}\overline{{Q}}\mathrm{=}{\mathrm{\nabla }}^{\mathrm{2}}\overline{{U}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{3}}{\nabla }{\nabla }\mathrm{·}\overline{{U}}\mathrm{+}\mathrm{2}\overline{{S}}{\nabla }\ln \overline{\mathit{\rho }}& & \end{array}$(27)is a term appearing in the mean viscous force $\mathrm{-}{\mathit{\nu }}_{\mathrm{T}}\overline{{Q}}$, where $\overline{{S}}$ is the traceless rate-of-strain tensor of the mean flow with components ${\overline{\mathrm{S}}}_{\mathit{ij}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathrm{\nabla }}_{\mathit{j}}{\overline{\mathit{U}}}_{\mathit{i}}\mathrm{+}{\mathrm{\nabla }}_{\mathit{i}}{\overline{\mathit{U}}}_{\mathit{j}}\mathrm{)}\mathrm{-}\frac{\mathrm{1}}{\mathrm{3}}{\mathit{\delta }}_{\mathit{ij}}{\nabla }\mathrm{·}\overline{{U}}$, and finally the term ${\nabla }\mathrm{(}{\mathit{q}}_{\mathrm{p}}\overline{{B}}{\mathrm{2}}^{}\mathit{/}\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}\mathrm{)}$ on the righthand side of Eq. (25) determines the turbulent contribution to the effective magnetic pressure. Here, q_p depends on the local field strength; see Eq. (19). This term enters with a plus sign, so positive values of q_p correspond to a suppression of the total turbulent pressure. The net effect of the mean magnetic field leads to an effective mean magnetic pressure that becomes negative for q_p > 1. This can indeed be the case for magnetic Reynolds numbers well above unity (Brandenburg et al. 2012); see also Fig. 7 for a polytropic atmosphere.

The boundary conditions for MFS are the same as for DNS, i.e., stress-free for the mean velocity at the top and bottom. For the mean magnetic field, we use either perfect conductor boundary conditions (for horizontal, imposed magnetic fields) or vertical field conditions (for vertical, imposed fields) at the top and bottom. All mean-field variables are assumed to be periodic in the x and y directions. The MFS are performed again with the Pencil Code, which is equipped with a mean-field module for solving the corresponding equations.

4.3. Expected vertical dependence of NEMPI

Fig. 9 Comparison of and λ0(z) in MFS for polytropic layers with four values of γ (top to bottom) and three values of β (different line types).

Fig. 10 Dependence of growth rate and height where the eigenfunction attains its maximum value (the optimal depth of NEMPI) on field strength from MFS for different values of γ in the presence of a horizontal field for Model I.

To get an idea about the vertical dependence of NEMPI, we now consider the resulting dependencies of ; see the lefthand panels of Fig. 9. We note that is just a function of β (Eqs. (18) and (19)), which allows us to approximate the local growth rates as (Rogachevskii & Kleeorin 2007; Kemel et al. 2013) ${\mathit{\lambda }}_{\mathrm{0}}\mathrm{=}\frac{{\mathit{v}}_{\mathrm{A}}}{{\mathit{H}}_{\mathit{\rho }\mathrm{0}}}{\left(\mathrm{-}\mathrm{2}\frac{\mathrm{d}{\mathrm{𝒫}}_{\mathrm{eff}}}{\mathrm{d}{\mathit{\beta }}^{\mathrm{2}}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathit{,}$(28)which are plotted in the righthand panels of Fig. 9 for Model I.

4.4. Horizontal fields

To analyze the kinematic stage of MFS, we measure the value of the maximum downflow speed, $\mathrm{|}\overline{\mathit{U}}{\mathrm{|}}_{\max }^{\mathrm{down}}$ at each height. We then determine the time interval during which the maximum downflow speed increases exponentially and when the height of the peak is constant and equal to z_B. This yields the growth rate of the instability as $\mathit{\lambda }\mathrm{=}\mathrm{d}\ln \mathrm{|}\overline{\mathit{U}}{\mathrm{|}}_{\max }^{\mathrm{down}}\mathit{/}\mathrm{d}\mathit{t}$.

In Figs. 10 and 11 we plot, respectively for Models I and II, λ (in units of ${\mathit{\eta }}_{\mathrm{t}}\mathit{/}{\mathit{H}}_{\mathit{\rho }\mathrm{0}}^{\mathrm{2}}$) and z_B versus horizontal imposed magnetic field strength, B₀/B_eq. The maximum growth rates for γ = 5/3 and γ = 1 are similar for both models (4–5 ${\mathit{\eta }}_{\mathrm{t}}\mathit{/}{\mathit{H}}_{\mathit{\rho }\mathrm{0}}^{\mathrm{2}}$ for Model I and 16–20 ${\mathit{\eta }}_{\mathrm{t}}\mathit{/}{\mathit{H}}_{\mathit{\rho }\mathrm{0}}^{\mathrm{2}}$ for Model II). It turns out that for γ = 5/3, the growth rate λ attains a maximum at some value B₀ = B_max, and then it decreases with increasing B₀, while in an isothermal run λ is nearly constant for greater field strength, except near the surface where the proximity to the boundary is too small. This close proximity reduces the growth rate. For Model I with γ = 1, the decline of λ (toward weaker fields on the lefthand side of the plot) begins when the distance to the top boundary (z_top − z_B ≈ 1.2 H_ρ0 for β₀ = 0.02) is less than the radius of magnetic structures (R ≈ 2π/k_⊥ ≈ 1.5 H_ρ0 using k_⊥H_ρ0 = 1). In Model II with γ = 1, the decline of λ occurs for stronger fields, but the distance to the top boundary (≈1.0 H_ρ0) is still nearly the same as for Model I.

Fig. 11 Same as Fig. 10 (horizontal field), but for Model II.

Fig. 12 Snapshots of ${\overline{\mathit{B}}}_{\mathit{y}}$ from MFS during the kinematic growth phase for different values of γ and B0/Beq0 for γ = 1 (top) and 5/3 (bottom) with β0 = 0.02 (left) and 0.1 (right) in the presence of a horizontal field for Model I.

In an isothermal layer, the height where the eigenfunction peaks is known to decrease with increasing field strength; see Fig. 6 of Kemel et al. (2012a). One might have expected this decrease to be less steep in the polytropic case, because the optimal depth where NEMPI occurs cannot easily be decreased without suffering a dramatic decrease of the growth rate. This is however not the case, and we find that the optimal depth of NEMPI is now falling off more quickly in Model I, but is more similar for Model II; see the second panels of Figs. 10 and 11. This means that in a polytropic layer, NEMPI works more effectively, and its growth rate is fastest when the magnetic field is not too strong. At the same time, the optimal depth of NEMPI increases, i.e., the resulting value of z_B increases as B₀ decreases.

The resulting growth rates are somewhat less for Model I and somewhat higher for Model II than those of earlier mean-field calculations of Kemel et al. (2012a, their Fig. 6), who found $\mathit{\lambda }{\mathit{H}}_{\mathrm{p}\mathrm{0}}^{\mathrm{2}}\mathit{/}{\mathit{\eta }}_{\mathrm{t}}\mathrm{\approx }\mathrm{9.7}$ in a model with β_⋆ = 0.32 and β_p = 0.05. These differences in the growth rates are plausibly explained by differences in the mean-field parameters.

Visualizations of the resulting horizontal field structures are shown in Fig. 12 for two values of γ and B₀. Increase in the parameter γ results in a stronger localization of the instability at the surface layer, where the density scale height is minimum and the growth of NEMPI is strongest.

Fig. 13 Dependence of growth rate and optimal depth of NEMPI on field strength from MFS for different values of γ in the presence of a vertical field for Model I.

4.5. Vertical fields

In the presence of a vertical field, the early evolution of the instability is similar to that for a horizontal field. In both cases, the maximum field strength occurs at a somewhat larger depth when saturation is reached, except that shortly before saturation there is a brief interval during which the location of maximum field strength rises slightly upwards in the vertical field case. In the saturated case, however, the flux concentrations from NEMPI are much stronger compared to the case of a horizontal field and it leads to the formation of magnetic flux concentrations of equipartition field strength (Brandenburg et al. 2013, 2014). This is possible because the resulting vertical flux tube is not advected downward with the flow that develops as a consequence of NEMPI. The latter effect is the aforementioned “potato-sack” effect, which acts as a nonlinear saturation mechanism of NEMPI with a horizontal field.

In Figs. 13 and 14 we plot the growth rates λ and the heights where the eigenfunction attains its maximum values for different β₀ = B₀/B_eq0 for Models I and II, respectively. For γ = 5/3, the maximum growth rate is higher larger than for γ = 1. This is true for Models I and II, where they are 8–10 ${\mathit{\eta }}_{\mathrm{t}}\mathit{/}{\mathit{H}}_{\mathit{\rho }\mathrm{0}}^{\mathrm{2}}$ for γ = 5/3 and 5–7 ${\mathit{\eta }}_{\mathrm{t}}\mathit{/}{\mathit{H}}_{\mathit{\rho }\mathrm{0}}^{\mathrm{2}}$ for γ = 1. The nonmonotonous behavior seen in the dependence of λ on B₀ is suggestive of the presence of different mode structures, although a direct inspection of the resulting magnetic field did not show any obvious differences. However, this irregular behavior may be related to artifacts resulting from a finite domain size and were not regarded important enough to justify further investigation.

Fig. 14 Same as Fig. 13 (vertical field), but for Model II.

Fig. 15 Comparison of ${\overline{\mathit{U}}}_{\mathrm{rms}}$ from MFS for Models I (solid line) and II (dashed line), showing exponential growth followed by nonlinear saturation. In both cases we have γ = 5/3 and an imposed vertical magnetic field with β0 = 0.05.

Fig. 16 Snapshots from MFS showing ${\overline{\mathit{B}}}_{\mathit{z}}$ on the periphery of the computational domain for γ = 5/3 and β0 = 0.05 at different times for Model I for the case of a vertical field.

Fig. 17 Similar to Fig. 16 of MFS, but for Model II at times similar to those in the DNS of Fig. 5. There are now more structures than in the earlier MFS of Fig. 16, and they develop more rapidly.

Next we focus on a comparison of the growth rates obtained from MFS for horizontal and vertical fields. The values of β₀, z_B, and β(z_B) = B₀/B_eq(z_B) for horizontal and vertical fields are compared in Table 1 for both models. We see that NEMPI is most effective in regions where the mean magnetic field is a small fraction of the local equipartition field and typically slightly less for γ = 5/3 than for γ = 1. Indeed, for Model I, β(z_B) is 3–4% for horizontal fields and 3–6% for vertical fields, while for Model II, β(z_B) is 12–17% for horizontal fields and 8–22% for vertical fields. Here, we have used β(z_B) = β₀e_2 − γ(z_B/2H_ρ0), where e_q(x) is the q-exponential function defined by Eq. (3).

We expect that higher values of β_⋆ will lead to greater growth rates. To verify this, we compare in Fig. 15 the time evolutions of ${\overline{\mathit{U}}}_{\mathrm{rms}}$ for Models I (with β_⋆ = 0.33) and Model II (with β_⋆ = 0.63). The growth rate has now increased by a factor of 2.4 (from $\mathit{\lambda }{\mathit{H}}_{\mathit{\rho }\mathrm{0}}^{\mathrm{2}}\mathit{/}{\mathit{\eta }}_{\mathrm{t}}\mathrm{=}\mathrm{3.9}$ to 9.5), which is slightly more than what is expected from β_⋆, which has increased by a factor of 1.9. This change in the growth rate can also be seen in Fig. 11 and Fig. 14 for Model II (in comparison with Figs. 10 and 13 for Model I). The dependence of the growth rate on the magnetic field strength is qualitatively similar for Models I and II. In particular, it becomes constant for γ = 1, but declines for γ = 5/3 as the field increases. The increase of z_B with B₀ is, however, less strong for Model II.

Snapshots of ${\overline{\mathit{B}}}_{\mathit{z}}$ from MFS for γ = 5/3 and β₀ = 0.05 at different times for Model I are shown in Fig. 16. Comparison with the results of MFS for Model II (see Fig. 17) shows that Model II fits the DNS better. This is also seen by comparing Fig. 17 with Fig. 5. However, our basic conclusions formulated in this paper are not affected.

5. Conclusions

The present work has demonstrated that in a polytropic layer, both in MFS and DNS, NEMPI develops primarily in the uppermost layers, provided the mean magnetic field is not too strong. If the field gets stronger, NEMPI can still develop, but the magnetic structures now occur at greater depths and the growth rate of NEMPI is lower. However, at some point when the magnetic field gets too strong, NEMPI is suppressed in the case of a polytropic layer, while it would still operate in the isothermal case, provided the domain is deep enough. The slow down of NEMPI is not directly a consequence of a longer turnover time at greater depths, but it is related to stratification being too weak for NEMPI to be excited.

By and large, the scaling relations determined previously for isothermal layers with constant scale height still seem to apply locally to polytropic layers with variable scale heights. In particular, the horizontal scale of structures was previously determined to be about 6–8 H_ρ (Kemel et al. 2013; Brandenburg et al. 2014). Looking now at Fig. 4, we see that for β₀ = 0.02 and γ = 5/3, the structures have a wavelength of ≈3 H_ρ0, but this is at a depth where H_ρ ≈ 0.3 H_ρ0. Thus, locally we have a wavelength of ≈10 H_ρ. The situation is similar in the next panel of Fig. 4 where the wavelength is ≈6 H_ρ0, and the structures are at a depth where H_ρ ≈ 1.5 H_ρ0, so locally we have a wavelength of ≈9 H_ρ. We can therefore conclude that our earlier results for isothermal layers can still be applied locally to polytropic layers.

A new aspect, however, that was not yet anticipated at the time, concerns the importance of NEMPI for vertical fields. While NEMPI with horizontal magnetic field still leads to downflows in the nonlinear regime (the “potato-sack” effect), our present work now confirms that structures consisting of vertical fields do not sink, but reach a strength comparable to or in excess of the equipartition value (Brandenburg et al. 2013, 2014). This makes NEMPI a viable mechanism for spontaneously producing magnetic spots in the surface layers. Our present study therefore supports ideas about a shallow origin for active regions and sunspots (Brandenburg 2005; Brandenburg et al. 2010; Kitiashvili et al. 2010; Stein & Nordlund 2012), contrary to common thinking that sunspots form near the bottom of the convective zone (Parker 1975; Spiegel & Weiss 1980; D’Silva & Choudhuri 1993). More specifically, the studies of Losada et al. (2013) point toward the possibility that magnetic flux concentrations form in the top 6 Mm, i.e., in the upper part of the supergranulation layer.

There are obviously many other issues of NEMPI that need to be understood before it can be applied in a meaningful way to the formation of active regions and sunspots. One question is whether the hydrogen ionization layer and the resulting H⁻ opacity in the upper layers of the Sun are important in providing a sharp temperature drop and whether this would enhance the growth rate of NEMPI, just like strong density stratification does. Another important question concerns the relevance of a radiating surface, which also enhances the density contrast. Finally, of course, one needs to verify that the assumption of forced turbulence is useful in representing stellar convection. Many groups have considered magnetic flux concentrations using realistic turbulent convection (Kitiashvili et al. 2010; Rempel 2011; Stein & Nordlund 2012). However, only at sufficiently large resolution can one expect strong enough scale separation between the scale of the smallest eddies and the size of magnetic structures. That is why forced turbulence has an advantage over convection. Ultimately, however, such assumptions should no longer be necessary. On the other hand, if scale separation is poor, our present parameterization might no longer be accurate enough, and one would need to replace the multiplication between q_p and $\overline{{B}}{\mathrm{2}}^{}$ in Eq. (25) by a convolution. This possibility is fairly speculative and requires a separate investigation. Nevertheless, in spite of these issues, it is important to emphasize that the qualitative agreement between DNS and MFS is already surprisingly good.

Acknowledgments

This work was supported in part by the European Research Council under the AstroDyn Research Project No. 227952, by the Swedish Research Council under the project grants 621-2011-5076 and 2012-5797 (I.R.L., A.B.), by EU COST Action MP0806, by the European Research Council under the Atmospheric Research Project No. 227915, and by a grant from the Government of the Russian Federation under contract No. 11.G34.31.0048 (N.K., I.R.). We acknowledge the allocation of computing resources provided by the Swedish National Allocations Committee at the Center for Parallel Computers at the Royal Institute of Technology in Stockholm and the National Supercomputer Centers in Linköping, the High Performance Computing Center North in Umeå, and the Nordic High Performance Computing Center in Reykjavik.

References

All Tables

All Figures

	Fig. 1 Isothermal and polytropic relations for different values of γ when calculated using the conventional formula ρ ∝ (z∞ − z)n with n = 1/(γ − 1).
In the text

	Fig. 2 Polytropes (Eq. (4)) with γ = 1 (solid line), 1.2 (dash-dotted), 1.4 (dotted), and 5/3 (dashed) and density scale height (Eq. (5)) for −4 ≤ −z/Hρ0 ≤ 1.2. The total density contrast is similar for γ = 1 and 5/3.
In the text

	Fig. 3 Sketch showing the expected size distribution of the nearly circular NEMPI eigenfunction structures at different heights.
In the text

	Fig. 4 Snapshots of ${\overline{\mathit{B}}}_{\mathit{y}}$ from DNS for γ = 5/3 and β0 = 0.02 (upper row), 0.05 (middle row), and 0.1 (lower row) at different times (indicated in turbulent-diffusive times, increasing from left to right) in the presence of a horizontal field using the perfect conductor boundary condition.
In the text

	Fig. 5 Snapshots from DNS showing ${\overline{\mathit{B}}}_{\mathit{z}}$ on the periphery of the computational domain for γ = 5/3 and β0 = 0.05 at different times for the case of a vertical field using the vertical field boundary condition.
In the text

	Fig. 6 Cuts of Bz/Beq(z) in (a) the xy plane at the top boundary (z/Hρ0 = 1.2) and (b) the xz plane through the middle of the spot at y = 0 for γ = 5/3 and β0 = 0.05. In the xz cut, we also show magnetic field lines and flow vectors obtained by numerically averaging in azimuth around the spot axis.
In the text

	Fig. 7 Effective magnetic pressure obtained from DNS in a polytropic layer with different γ for horizontal (H, red curves) and vertical (V, blue curves) mean magnetic fields.
In the text

	Fig. 8 Parameters qp0, βp, and β⋆ for the function qp(β) (see Eq. (19)) versus γ for horizontal (red line) and vertical (blue line) mean magnetic fields.
In the text

	Fig. 9 Comparison of and λ0(z) in MFS for polytropic layers with four values of γ (top to bottom) and three values of β (different line types).
In the text

	Fig. 10 Dependence of growth rate and height where the eigenfunction attains its maximum value (the optimal depth of NEMPI) on field strength from MFS for different values of γ in the presence of a horizontal field for Model I.
In the text

	Fig. 11 Same as Fig. 10 (horizontal field), but for Model II.
In the text

	Fig. 12 Snapshots of ${\overline{\mathit{B}}}_{\mathit{y}}$ from MFS during the kinematic growth phase for different values of γ and B0/Beq0 for γ = 1 (top) and 5/3 (bottom) with β0 = 0.02 (left) and 0.1 (right) in the presence of a horizontal field for Model I.
In the text

	Fig. 13 Dependence of growth rate and optimal depth of NEMPI on field strength from MFS for different values of γ in the presence of a vertical field for Model I.
In the text

	Fig. 14 Same as Fig. 13 (vertical field), but for Model II.
In the text

	Fig. 15 Comparison of ${\overline{\mathit{U}}}_{\mathrm{rms}}$ from MFS for Models I (solid line) and II (dashed line), showing exponential growth followed by nonlinear saturation. In both cases we have γ = 5/3 and an imposed vertical magnetic field with β0 = 0.05.
In the text

	Fig. 16 Snapshots from MFS showing ${\overline{\mathit{B}}}_{\mathit{z}}$ on the periphery of the computational domain for γ = 5/3 and β0 = 0.05 at different times for Model I for the case of a vertical field.
In the text

	Fig. 17 Similar to Fig. 16 of MFS, but for Model II at times similar to those in the DNS of Fig. 5. There are now more structures than in the earlier MFS of Fig. 16, and they develop more rapidly.
In the text