Spontaneous flux concentrations from the negative effective magnetic pressure instability beneath a radiative stellar surface

TL;DR

This study investigates the negative magnetic pressure instability (NEMPI) in a radiative, stratified model to understand magnetic flux concentration formation beneath stellar surfaces, revealing oscillatory flux structures that are weaker and narrower than sunspots.

Contribution

It introduces a radiative, stratified model with a free surface to study NEMPI, highlighting effects of radiation on flux concentration formation in a more realistic stellar context.

Findings

NEMPI produces oscillatory flux concentrations beneath the surface.

Flux structures are narrower and weaker than sunspots.

Radiation influences the development and characteristics of NEMPI.

Abstract

The formation of sunspots requires the concentration of magnetic flux near the surface. The negative magnetic pressure instability (NEMPI) might be a possible mechanism for accomplishing this, but it has mainly been studied in simple systems using an isothermal equation of state without a natural free surface. We study NEMPI in a stratified Cartesian mean-field model where turbulence effects are parameterized. We use an ideal equation of state and include radiation transport, which establishes selfconsistently a free surface. We use a Kramers-type opacity with adjustable exponents chosen such that the deeper layers are approximately isentropic. No convection is therefore possible in this model, allowing us to study NEMPI with radiation in isolation. We restrict ourselves to two-dimensional models. We use artificially enhanced mean-field coefficients to allow NEMPI to develop, making it therefore possible to study the reason why it is much harder to excite in the presence of radiation. NEMPI yields moderately strong magnetic flux concentrations a certain distance beneath the surface where optical depth is unity. The instability is oscillatory and in the form of upward travelling waves. This seems to be a new effect that has not been found in earlier models without radiative transport. The horizontal wavelength is about ten times smaller than what has been found previously in more idealized isothermal models. In our models, NEMPI saturates at field strengths too low to explain sunspots. Furthermore, the structures appear too narrow and too far beneath the surface to cause significant brightness variations at the radiative surface. We speculate that the failure to reproduce effects resembling sunspots may be related to the neglect of convection.

Tables

Table 1. Table 2: Summary of Runs A–H.

	$L_x$	$B_0$	${\kappa }_0$	${\eta }_t$	$N$	$\lambda$	$\omega$	$z_B$
A	$1.25$	100	${10}^7$	${\mathrm{5\hspace{0.33em}10}}^{-3}$	${288}^2$	0.011	0.79	3.0
A′	$2.5$	100	${10}^7$	${\mathrm{5\hspace{0.33em}10}}^{-3}$	${288}^2$	0.009	0.89	2.5
B	$1.25$	200	${10}^7$	${\mathrm{5\hspace{0.33em}10}}^{-3}$	${288}^2$	0.030	1.44	2.7
B′	$2.5$	200	${10}^7$	${\mathrm{5\hspace{0.33em}10}}^{-3}$	${288}^2$	0.017	1.42	2.5
B′′	$5.0$	200	${10}^7$	${\mathrm{5\hspace{0.33em}10}}^{-3}$	${288}^2$	0.049	1.21	2.5
C′	$2.5$	500	${10}^7$	${\mathrm{5\hspace{0.33em}10}}^{-3}$	${288}^2$	0.021	0.90	1.7
D	$1.25$	200	${\mathrm{2\hspace{0.33em}10}}^7$	${\mathrm{5\hspace{0.33em}10}}^{-3}$	${288}^2$	0.022	1.51	2.8
E	$1.25$	200	${\mathrm{5\hspace{0.33em}10}}^7$	${\mathrm{5\hspace{0.33em}10}}^{-3}$	${288}^2$	0.043	1.55	3.2
F	$1.25$	200	${\mathrm{2\hspace{0.33em}10}}^7$	${\mathrm{2\hspace{0.33em}10}}^{-3}$	${576}^2$	0.094	0.70	2.7
G	$1.25$	200	${\mathrm{5\hspace{0.33em}10}}^7$	${\mathrm{2\hspace{0.33em}10}}^{-3}$	${576}^2$	0.101	0.83	2.8
H	$1.25$	200	${10}^8$	${\mathrm{2\hspace{0.33em}10}}^{-3}$	${576}^2$	0.152	0.70	3.1

Equations

$$\frac{\partial \overline{\mathbf{A}}}{\partial t}=\overline{\mathbf{U}}\times \overline{\mathbf{B}}+{\eta }_{\mathrm{T}}{\nabla }^2\overline{\mathbf{A}},$$

$$\overline{\rho}{{\rm D}{}\overline{\bm{U}}\over{\rm D}{}t}=-{\bm{\nabla}}\left(\overline{P}-{q_{\rm p}{\overline{\bm{B}}}^{2}\over 2\mu_{0}}\right)+{\overline{\bm{J}}}\times{\overline{\bm{B}}}+\overline{\rho}\bm{g}+{\bm{\nabla}}\cdot(2\nu_{\rm T}\overline{\rho}\,\overline{\mbox{\boldmath${\mathsf{S}}$}}{}),\;$$

$$\overline{\rho}\overline{T}{{\rm D}{}\overline{s}\over{\rm D}{}t}=-{\bm{\nabla}}\cdot(\bm{F}_{\rm rad}+\bm{F}_{\rm conv})+2\nu_{\rm T}\overline{\rho}\,\overline{\mbox{\boldmath${\mathsf{S}}$}}{}^{2},$$

$$\frac{\mathrm{D}\ln \overline{\rho }}{\mathrm{D}t}=-\nabla \cdot \overline{\mathbf{U}},$$

$$\hat{\mathbf{n}}\cdot \nabla I=-\kappa \overline{\rho }(I-S)$$

$$\nabla \cdot {\mathbf{F}}_{\mathrm{rad}}=-\kappa \overline{\rho }{\oint }_{4\pi }(I-S)\mathrm{d}\Omega ,$$

$$q_{\rm p}(\beta)={q_{\rm p0}\over 1+\beta^{2}/\beta_{\rm p}^{2}}={\beta_{\star}^{2}\over\beta_{\rm p}^{2}+\beta^{2}},\quad\mbox{where $\beta_{\star}=\beta_{\rm p}q_{\rm p0}^{1/2}$}.$$

$$\kappa ={\kappa }_0(\overline{\rho }/{\rho }_0{)}^a(\overline{T}/T_0{)}^b$$

$$K=K_0\frac{(\overline{T}/T_0{)}^{3-b}}{(\overline{\rho }/{\rho }_0{)}^{1+a}}=K_0{\left[\frac{(\overline{T}/T_0{)}^n}{\overline{\rho }/{\rho }_0}\right]}^{1+a},$$

$$K_0=16{\sigma }_{\mathrm{SB}}{T}_0^3/3{\kappa }_0{\rho }_0$$

$$n=(3-b)/(1+a)$$

$$\partial\overline{U}_{x}/\partial z=\partial\overline{U}_{y}/\partial z=\overline{U}_{z}=0\quad\mbox{on $\;z=0$, $L_{z}$},$$

$$\partial\overline{A}_{x}/\partial z=\partial\overline{A}_{y}/\partial z=\overline{A}_{z}=0\quad\mbox{on $\;z=0$, $L_{z}$}.$$

$$\overline{T}=T_{0}\quad\mbox{on $z=0$},$$

$$\partial\overline{T}/\partial z=-\sigma_{\rm SB}\overline{T}^{4}\quad\mbox{(on $z=d$)}.$$

$$d=c_{\mathrm{p}}(T_0-T_1){\nabla }_{\mathrm{ad}}/(g\nabla ),$$

$$T_1={\mathcal{K}}^{1/4}{T}_0\text{mbox}with\mathcal{K}=\frac{g{K}_0}{c_{\mathrm{p}}{\sigma }_{\mathrm{SB}}{T}_{\mathrm{bot}}^4}\frac{\nabla }{{\nabla }_{\mathrm{ad}}},$$

$$k_{\mathrm{f}}/k_1=u_{\mathrm{rms}}/3{\eta }_{\mathrm{t}}{k}_1.$$

$$\ell =(\kappa \rho {)}^{-1},$$

$${\mathcal{P}}_{\mathrm{eff}}(\beta )=\frac{1}{2}[1-q_{\mathrm{p}}(\beta )]{\beta }^2,$$

Full text

*Aims. *We study NEMPI in a stratified Cartesian mean-field model where turbulence effects are parameterized. We use an ideal equation of state and include radiation transport, which establishes selfconsistently a free surface.

*Methods. *We use a Kramers-type opacity with adjustable exponents chosen such that the deeper layers are approximately isentropic. No convection is therefore possible in this model, allowing us to study NEMPI with radiation in isolation. We restrict ourselves to two-dimensional models. We use artificially enhanced mean-field coefficients to allow NEMPI to develop, making it therefore possible to study the reason why it is much harder to excite in the presence of radiation.

*Results. *NEMPI yields moderately strong magnetic flux concentrations a certain distance beneath the surface where the optical depth is unity. The instability is oscillatory and in the form of upward travelling waves. This seems to be a new effect that has not been found in earlier models without radiative transport. The horizontal wavelength is about ten times smaller than what has previously been found in more idealized isothermal models.

*Conclusions. *In our models, NEMPI saturates at field strengths too low to explain sunspots. Furthermore, the structures appear too narrow and too far beneath the surface to cause significant brightness variations at the radiative surface. We speculate that the failure to reproduce effects resembling sunspots may be related to the neglect of convection.

1 Introduction

Sunspots are a highly intermittent manifestation of strong magnetic flux concentrations at the solar surface. The underlying magnetic fields are produced by hydromagnetic turbulence in the convection zone beneath the surface (Brun et al., 2004; Brown et al., 2011; Käpylä et al., 2012a; Augustson et al., 2015). Numerous simulations of turbulence and turbulent convection have displayed such magnetic field production by a dynamo process (Brandenburg and Subramanian, 2005). This alone, however, does not explain the occasional concentration into spots. On the other hand, more realistic simulations that include the effects of strong density stratification near the surface, as well as radiation and ionization, have been able to demonstrate the appearance of magnetic spots (Stein & Nordlund, 2012). Furthermore, it is known that strong density stratification can lead to a large-scale instability of an initially unstructured random magnetic field (Kleeorin et al., 1989, 1990). This instability leads to magnetic flux concentrations and even magnetic spots (Brandenburg et al., 2013; Warnecke et al., 2013) through the negative magnetic pressure instability (NEMPI). NEMPI has been associated with sunspot formation by Kleeorin et al. (1995, 1996), following a series of earlier work on its theoretical foundations (Kleeorin et al., 1993; Kleeorin and Rogachevskii, 1994; Rogachevskii & Kleeorin, 2007).

Numerical simulations have also displayed types of magnetic flux concentrations that are not straightforwardly associated with NEMPI. This tends to be the case when the magnetic field is produced by a large-scale dynamo some distance beneath the surface (Mitra et al., 2014; Jabbari et al., 2015, 2016). Nevertheless, also in those cases strong stratification was shown to be essential, as has been demonstrated by comparing with weakly stratified cases.

In the case of NEMPI, the underlying instability can well be modeled using mean-field magnetohydrodynamics, where the negative effective magnetic pressure is parameterized in terms of the mean magnetic field; see Brandenburg et al. (2016) for a review. In some of those cases there is good quantitative agreement between direct numerical simulations (DNS) and mean field simulations (MFS), as has been demonstrated in several papers (Kemel et al., 2013; Losada et al., 2013).

The main difference between MFS and DNS is the inclusion of the parameterization of the small-scale unresolved motions $\bm{u}=\bm{U}-\overline{\bm{U}}$ and magnetic fields $\bm{b}=\bm{B}-{\overline{\bm{B}}}$ in the MFS. Here, the overbar denotes a suitably defined average, which, in practice, could be a spatial average. The evolution equations for $\overline{\bm{U}}$ involve correlations of the form $\overline{u_{i}u_{j}}$ and $\overline{b_{i}b_{j}}$ that need to be expressed in terms of $\overline{\bm{U}}$ and ${\overline{\bm{B}}}$. They are similar to the parameterization in terms of the rate-of-strain tensor of the mean flow involving turbulent viscosity, but there are also contributions that are quadratic in ${\overline{\bm{B}}}$. Similar parameterizations also exist for the Maxwell stresses in the momentum equation and the electromotive force in the induction equation.

Before we can think of applying NEMPI to real sunspot formation, we must begin to address the effects of radiation, ionization, and other potentially important surface effects. Here, we focus on radiative transfer. Radiation has two important effects. On the one hand, it leads to the establishment of a natural surface from which most of the observed radiation is emitted and above which the density drops off sharply. On the other hand, radiation also leads to the equilibration of temperature differences between neighboring fluid elements. Earlier investigations have suggested that this may indeed be the case and that NEMPI may be difficult to excite in the presence of radiation (Barekat, 2013; Bhat & Brandenburg, 2016). This is also the reason why we focus here on mean-field simulations, because they allow us to artificially exaggerate the effects of NEMPI by choosing unrealistically large mean-field parameters, which allows us to study the properties of NEMPI in that case and helps us determining the conditions under which NEMPI may still operate.

Most of the earlier investigations of NEMPI have been carried out in an isothermally stratified layer using an isothermal equation of state. This means that no energy equation was solved. This was also true in simulations with an outer coronal envelope (Warnecke et al., 2013), where the interface was characterized by a layer above which the driving of turbulence was turned off. The aim of the present paper is therefore to study NEMPI in a simple model with radiative heating and cooling included.

2 The model

2.1 Mean-field equations and radiative transfer

We consider the mean-field equations in Cartesian coordinates, but restrict ourselves to including only the effects of turbulent magnetic diffusion, turbulent viscosity, and the negative effective magnetic pressure effect, which means that the ordinary magnetic pressure from the mean field, ${\overline{\bm{B}}}^{2}/2\mu_{0}$, is modified and becomes $(1-q_{\rm p}){\overline{\bm{B}}}^{2}/2\mu_{0}$, where $q_{\rm p}=q_{\rm p}({\overline{\bm{B}}})$ depends on the local magnetic field strength. We write the mean magnetic field as ${\overline{\bm{B}}}=\bm{B}_{0}+{\bm{\nabla}}\times{\overline{\bm{A}}}$, where $\bm{B}_{0}=(0,0,B_{0})$ is an imposed vertical field, and ${\overline{\bm{A}}}$ is the mean magnetic vector potential. We thus solve the equations for ${\overline{\bm{A}}}$, the mean velocity $\overline{\bm{U}}$, the mean specific entropy $\overline{s}$, and the mean density $\overline{\rho}$ in the form

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\eta_{\rm T}=\eta+\eta_{\rm t}$ is the total magnetic diffusivity consisting of a microphysical and a turbulent value, $\nu_{\rm T}=\nu+\nu_{\rm t}$ is the total viscosity consisting again of a microphysical and a turbulent value, $\overline{\sf S}_{ij}={\textstyle{1\over 2}}(\overline{U}_{i,j}+\overline{U}_{j,i})-{\textstyle{1\over 3}}\delta_{ij}{\bm{\nabla}}\cdot\overline{\bm{U}}$ is the traceless rate-of-strain tensor, ${\overline{\bm{J}}}={\bm{\nabla}}\times{\overline{\bm{B}}}/\mu_{0}$ is the Lorentz force from the mean fields (without the effects of turbulence that are being parameterized through $q_{\rm p}$), $\bm{g}=(0,0,-g)$ is the gravitational acceleration, $\overline{P}$ is the mean gas pressure, $\overline{T}$ is the mean temperature, ${\rm D}{}/{\rm D}{}t=\partial/\partial t+\overline{\bm{U}}\cdot{\bm{\nabla}}$ is the advective derivative, $\bm{F}_{\rm rad}$ is the radiative flux, and $\bm{F}_{\rm conv}$ is the convective flux, but it will be neglected in our present exploratory work.

The radiative flux divergence is obtained by solving the radiative transfer equations for the intensity $I(\bm{x},t,\hat{\bm{n}})$ in the gray approximation in the form (Nordlund, 1982)

[TABLE]

along a set of rays in different directions $\hat{\bm{n}}$, where $\kappa$ is the opacity and $S=(\sigma_{\rm SB}/\pi)\,\overline{T}^{4}$ is the source function with $\sigma_{\rm SB}$ being the Stefan–Boltzmann constant. The radiative flux divergence is found by integrating Eq. (5) over all directions, i.e.,

[TABLE]

We adopt the equation of state for a perfect gas, i.e., $\overline{P}=({\cal R}/\mu)\overline{T}\,\overline{\rho}$, where ${\cal R}$ is the universal gas constant and $\mu$ the mean specific weight. The mean specific entropy is, up to an irrelevant additive constant, given by $\overline{s}/c_{p}=(\ln\overline{P})/\gamma-\ln\overline{\rho}$, where $\gamma=c_{\rm p}/c_{\rm v}$ is the ratio of specific heats at constant pressure and constant density, respectively, and ${\cal R}/\mu=c_{\rm p}-c_{\rm v}$. In the following, we take $\gamma=5/3$ which is appropriate for a monatomic gas and in the absence of ionization. The pressure scale height, $H_{\rm p}=-{\rm d}{}\ln\overline{P}/{\rm d}{}z$, is then given by $H_{\rm p}={\cal R}\overline{T}/\mu g$. In the isothermal part near the top, pressure and density scale heights are equal, i.e., $H_{\rho}=H_{\rm p}$, where $H_{\rho}=-{\rm d}{}\ln\overline{\rho}/{\rm d}{}z$. However, in the deeper isentropic parts, we have $H_{\rho}=\gamma H_{\rm p}$.

2.2 Parameterizations

Turbulence effects such as NEMPI depend on the relative importance of the magnetic field to the equipartition field strength with respect to the turbulent energy, that is, on $\beta\equiv|{\overline{\bm{B}}}|/B_{\rm eq}$. Here, the equipartition field strength $B_{\rm eq}$ is given by $B_{\rm eq}^{2}(z)=\mu_{0}\overline{\rho}u_{\rm rms}^{2}$. The effective magnetic pressure is characterized by the functional form of $q_{\rm p}=q_{\rm p}(\beta)$, for which we assume (Kemel et al., 2012)

[TABLE]

In addition, we have to specify $\eta_{\rm T}$ and $\nu_{\rm T}$, which we assume to be constant and equal to each other, i.e., we assume the turbulent magnetic Prandtl number $\mbox{\rm Pr}_{M}=\nu_{\rm T}/\eta_{\rm T}$ to be unity (Yousef et al., 2003). We define a fiducial model where we take $q_{\rm p0}=300$ and $\beta_{\rm p}=0.05$. Earlier work of Kemel et al. (2013) showed that the growth rate is mainly dependent on the parameter $\beta_{\star}$, whose value is then 0.87. For comparison, Kemel et al. (2012) and Käpylä et al. (2012b) used the parameter combination $q_{\rm p0}=40$ and $\beta_{\rm p}=0.05$, which then yields about a third for $\beta_{\star}=0.32$. Our value of $\beta_{\star}$ is thus much higher than what has been assumed before, which should help us to study the effects of radiation in the development of NEMPI. Following earlier work of Barekat (2013) and Barekat & Brandenburg (2014), we assume a Kramers-like opacity law for $\kappa$ of the form

[TABLE]

with constant coefficients $\kappa_{0}$, $\rho_{0}$, and $T_{0}$, and given exponents $a$ and $b$. The resulting radiative conductivity is then given by (Barekat & Brandenburg, 2014)

[TABLE]

where

[TABLE]

is a constant and

[TABLE]

is, for $n>-1$, related to the polytropic index of the resulting stratification. The radiative diffusivity is $\chi=K/\overline{\rho}c_{\rm p}$. The optical depth is $\tau(z)=\int_{z}^{\infty}\kappa\overline{\rho}\,{\rm d}{}z^{\prime}$. The region where $\tau\ll 1$ is optically thin, while the region where $\tau\gg 1$ is optically thick, which corresponds to the convection zone in the Sun; $\tau=1$ represents thus the solar surface.

2.3 Boundary conditions and numerical aspects

We adopt impenetrable stress-free vertical field boundary conditions in the $z$ direction, so the velocity obeys

[TABLE]

where $L_{z}$ is the vertical extent of the computational domain and the bottom boundary is at $z=0$. For the magnetic field we adopt the vertical field condition,

[TABLE]

We assume zero incoming intensity at the top, and compute the incoming intensity at the bottom from a quadratic Taylor expansion of the source function, which implies that the diffusion approximation is obeyed; see Appendix A of Heinemann et al. (2006) for details. As in Barekat & Brandenburg (2014), we fix the temperature at the bottom,

[TABLE]

while the temperature at the top is allowed to evolve freely. There is no boundary condition on the density, but since no mass is flowing in or out, the volume-averaged density is automatically constant; see Appendix C of Barekat & Brandenburg (2014).

To reduce the computational expense, we solve Eqs. (1)–(6) in two spatial dimensions. In an earlier investigation of NEMPI, Losada et al. (2012) found that this simplification can lead to about two times smaller growth rates, but the qualitative dependencies on various input parameters were still reproduced correctly. In the present model, we use either $288$ or $576$ meshpoints in the $z$ direction. The number of mesh points in the $x$ direction depends on the domain size and is constrained such that the mesh spacings $\delta x$ and $\delta z$ are equal in the two directions. We employ the Pencil Code111 https://github.com/pencil-code, where all relevant terms are readily implemented. The code uses a high-order finite-difference scheme. The radiation module was implemented by Heinemann et al. (2006).

2.4 Comparison with the optically thick approximation

It will be instructive to compare with the more familiar case in which $F_{\rm rad}$ is computed in the optically thick approximation as $\bm{F}_{\rm rad}=-K{\bm{\nabla}}\overline{T}$ in a domain $0\leq z\leq d$, where $d$ is less than the $L_{z}$ used in the general case with full radiative transfer. At $z=d$, we apply a radiative boundary condition

[TABLE]

The value of $d$ is computed from (see Sect. 3.12 of Barekat & Brandenburg, 2014) as

[TABLE]

where $\nabla_{\rm ad}=1-1/\gamma$ and $\nabla=1+1/n$ with $n$ given by Eq. (11), and

[TABLE]

where $K_{0}$ is given by Eq. (10). The quantities in Eqs. (16)–(17) are fully determined by the parameters of the radiative model. We emphasize that the temperature at the top is close to $T_{1}$, but it allowed to evolve freely subject to Eq. (15). Computationally, the optically thick approximation is by about a factor of two cheaper, but it is more restrictive, because the values of $d$ and $T_{1}$ are intimately tied to the choice of $\kappa_{0}$ and cannot be varied independently.

2.5 Scale separation ratio

In our mean-field model, turbulence is parameterized in terms of a magnetic turbulent diffusivity, which is estimated to be $\eta_{\rm t}=u_{\rm rms}/3k_{\rm f}$ (Sur et al., 2008), where $k_{\rm f}$ is the wavenumber of the energy-carrying motions. We compare this with the reference wavenumber $k_{1}=2\pi/L_{z}$ based on our domain of height $L_{z}$. We refer to $k_{\rm f}/k_{1}$ as the scale separation ratio. Thus, we have (Jabbari et al., 2014)

[TABLE]

This ratio must be large enough for NEMPI to be excited (Brandenburg et al., 2012). Early DNS of Brandenburg et al. (2011), where NEMPI was excited, used $k_{\rm f}/k_{1}=15$, but with $k_{\rm f}/k_{1}=30$, NEMPI became much more pronounced (Käpylä et al., 2012b).

Following the work of Barekat & Brandenburg (2014), we measure length in $\,{\rm Mm}$, velocity in $\,{\rm km}\,{\rm s}^{-1}$, and density in $\,{\rm g}\,{\rm cm}^{-3}$. We choose $L_{z}=5\,{\rm Mm}$. We adopt a squared domain, $L_{x}=L_{z}$, and assume for the turbulent small-scale velocity $u_{\rm rms}=1\,{\rm km/s}$. Thus, we have $\eta_{\rm t}=5\times 10^{-3}\,{\rm Mm}\,{\rm km/s}$, so we have $k_{\rm f}/k_{1}=53$, which should be large enough for NEMPI to be excited (Brandenburg et al., 2012). In some models with larger resolution ($576^{2}$ meshpoints), we used $\eta_{\rm t}=2\times 10^{-3}\,{\rm Mm}\,{\rm km/s}$, corresponding to $k_{\rm f}/k_{1}=133$; see Table LABEL:units for the conversion of several quantities from code units to cgs units. Following earlier work (Brandenburg et al., 2011), we also define the general turbulent–diffusive time $\tau_{\rm td}=(\eta_{\rm t}k_{1}^{2})^{-1}$.

3 Results

We design the model such that it has an isentropic deeper part. The stratification in our model is similar to Run B7 of Barekat & Brandenburg (2014) with $a=1$ and $b=0$, which, as discussed above, yields $n=1.5$. In particular, we use $\kappa_{0}=10^{7}\,{\rm Mm}^{-1}\,{\rm cm}^{3}\,{\rm g}^{-1}$, which results in a surface temperature of around $5000\,{\rm K}$. As in Barekat & Brandenburg (2014), we compute a hydrostatic equilibrium solution ($\bm{u}={\bf 0}$) by solving Eqs. (1)–(6) only in the $z$ direction in one dimension. The result is shown in Fig. 1, where we plot the $z$ dependence of $\overline{\rho}$, $\overline{s}$, $\overline{T}$, and $\chi$. In the deeper parts, where $\tau\gg 1$, $\overline{T}$ increases linearly with depth and, because $\overline{s}$ is nearly constant in that part, $\overline{\rho}(z)\propto\overline{T}^{3/2}$, which is in agreement with the expected polytropic stratification. Above the surface, $\overline{T}(z)$ is approximately constant, so $\overline{\rho}(z)$ falls off exponentially with height, as expected for an isothermal stratification. We begin by discussing in some detail a run with $200\,{\rm G}$, which will later also be referred to as Run B*′′*; see Table 2. The presence of an imposed field changes the stratification, but this change is small: $T$ decreases by $\approx 4\,{\rm K}$ for $B_{0}=200\,{\rm G}$.

3.1 Early evolution into saturation

In the early phase of the evolution, structures form where $\max(|{\overline{\bm{B}}}|)$ is at $z=z_{B}\approx 2.5\,{\rm Mm}$ and a horizontal wavenumber $k=4\,k_{1}$; see Fig. 2. These structures gradually move downward, disappear, and new ones form at $z\approx 3\,{\rm Mm}$. Those structures then also move downward, and so on. The structures occur well below the $\tau=1$ line and are close to the $\tau=100$ line. Here, the photon mean-free path,

[TABLE]

is about $0.05\,{\rm Mm}$, while at $\tau=1$, it is about $0.14\,{\rm Mm}$. The downward motions are associated with a local field enhancements, as can clearly be seen from field lines getting more concentrated in some locations. At later times, the field becomes more irregular, but retains a typical horizontal wavenumber of $4\,k_{1}$. In some cases, however, we found that, in the late nonlinear stage, $k$ can decrease from four to three.

The growth of structures can be characterized both by the typical velocities $\overline{\bm{U}}$ in the domain and the departures from the imposed field $\Delta{\overline{\bm{B}}}={\overline{\bm{B}}}-\bm{B}_{0}$. In Fig. 3, we show for two independent realizations of Run B*′′* (with $L_{x}=L_{z}$) the evolution of the rms values, $\overline{U}_{\rm rms}$ and $\Delta\overline{B}_{\rm rms}$, with different seeds for the random initial velocity perturbations. We clearly see an oscillatory growth of both quantities, as can also be seen by showing a plot compensated by $\exp(-\lambda t)$, where $\lambda\approx 0.048\,{\rm ks}^{-1}$ is the growth rate, as determined during the exponential growth phase of the instability. In the following, we measure the period $P_{\rm osc}$ as the volume-integrated rms velocity of the mean field and record the frequency $\omega=2\pi/P_{\rm osc}$. The frequency of the actual (signed) magnetic field is half that value.

The growth rate is independent of the initial seed for the random number generator, but the detailed nonlinear evolution does depend on it (compare the lines in each of the panels of Fig. 3). This suggests that the evolution of NEMPI is chaotic in the nonlinear regime. Animations show that the field lines are constantly swinging back and forth. This type of time-dependence of NEMPI is new and has not previously been seen – neither in isothermal nor in polytropic calculations. It may therefore be an effect related to the presence of radiation. The slight apparent difference in oscillation amplitudes of the compensated plots in the insets is caused by the fact that both have been compensated by the same factor, but the amplitudes were slightly different by the time the eigenfunction begins to be established.

The spatio-temporal evolution of NEMPI is seen more clearly in Fig. 4, where we show $\overline{U}_{z}(x_{\ast},z,t)$ for $x_{\ast}=-1.7\,{\rm Mm}$. This position $x_{\ast}$ is where $\overline{U}_{z}$ has an anti-node during the linear growth phase. At $t=1740\,{\rm ks}$, the slope in the $tz$ diagram corresponds to a pattern speed of about $0.2\,{\rm km}\,{\rm s}^{-1}$, which is small compared with the sound speed $c_{\rm s}\approx 20\,{\rm km}\,{\rm s}^{-1}$ at $z\approx 2.5\,{\rm Mm}$, and it is also small compared with the turbulent rms velocity of $1\,{\rm km}\,{\rm s}^{-1}$, but it agrees with the typical NEMPI-produced downflow speeds found earlier for isothermal NEMPI experiments (Brandenburg et al., 2014).

3.2 Dependence on control parameters

We now consider the dependence of NEMPI on $B_{0}$, $\kappa_{0}$, and $\eta_{\rm t}$. We revisit some of these dependencies later in more detail. Several input and output parameters of our runs are summarized in Table 2. Although most of the runs discussed in this paper are performed for a domain with $L_{x}=L_{z}=5\,{\rm Mm}$, several aspects can also be reproduced in narrower domains with $L_{x}/L_{z}=0.5$ and $0.25$. The growth rate $\lambda$ is rather sensitive to this, while the oscillation frequency $\omega=2\pi/P_{\rm osc}$ and the position $z_{B}$ of the magnetic field maximum are less sensitive. For Runs A–C with $\kappa_{0}=10^{7}\,{\rm Mm}^{-1}\,{\rm cm}^{3}\,{\rm g}^{-1}$, the growth rates are roughly in the range between $\lambda=0.01\,{\rm ks}^{-1}$ and $0.03\,{\rm ks}^{-1}$ and do not seem to be systematically dependent on the value of $B_{0}$. This is mainly related to the fact that NEMPI can develop deeper down as $B_{0}$ is increased; see Kemel et al. (2012). This is characterized by the value of $z_{B}$ given in Table 2; compare especially with the value for Run C. Our results thus confirm that the structures develop at larger depths when the field becomes stronger. This is in agreement with earlier work (Kemel et al., 2012; Losada et al., 2014).

3.3 Comparison with earlier work

In units of $\tau_{\rm td}$ (defined in Sect. 2.5), the growth rate is $\tilde{\lambda}\equiv\lambda\tau_{\rm td}=\lambda/\eta_{\rm t}k_{1}^{2}$, which is about 6 for Run B*′′*. However, if we normalize instead by actual horizontal wavenumber $k$ of the structures, which is 4 times larger than $k_{1}$ (see Fig. 2), we have $\lambda/\eta_{\rm t}k^{2}\approx 0.4$. This value is rather low and comparable to the value in the first DNS of Brandenburg et al. (2011), where the scale separation ratio was much lower ($k_{\rm f}/k_{1}=15$ compared to $53$ in the present case).

Earlier work using isothermal layers has shown that the horizontal wavenumber of the instability is comparable to the inverse density scale height; Kemel et al. (2013) found $kH_{\rho}=1.1$–$1.5$. Subsequent work showed that during the nonlinear evolution of the instability, $kH_{\rho}$ can decrease from about 0.8 to 0.2. This has been associated with an inverse cascade-type behavior (Brandenburg et al., 2014). The polytropic simulations of Losada et al. (2014) gave larger values: $kH_{\rho}=1$ in the upper layers and $kH_{\rho}=2$ in deeper ones; see their Fig. 12. In the present case, at the height where the instability based on the absolute field strength is strongest ($z=3\,{\rm Mm}$), and for $k/k_{1}=4$, we find $kH_{\rm p}=5$; see the dashed line in Fig. 5a. This is a striking difference between the present models and the earlier ones using an isothermal equation of state.

Figure 5b shows that, in the region where NEMPI develops, the equipartition field strength $B_{\rm eq}(z)$ is around $3000$–$4000\,{\rm G}$. This is about 20 times larger than the strength of the imposed field, which is typical of NEMPI and in agreement with earlier results (Losada et al., 2014; Brandenburg et al., 2014).

3.4 Magnetic field dependence

As alluded to above, there are several other aspects of NEMPI that can be compared with what has been found earlier. We now compare our results with Fig. 6 of Brandenburg et al. (2014), where the vertical dependence of the maximum field in the structures, $B_{\max}(z)$, was plotted, normalized either by $B_{0}$ or by $B_{\rm eq}(z)$. The corresponding result for our present simulations is shown in Fig. 6. The local maxima in $B_{\max}(z)$ are caused by the spatial wave-like structures seen in Figs. 2 and 4.

Unlike the earlier work for isothermal layers, where the slope of $B_{\rm eq}/B_{0}$ was constant, it varies in the present case. More importantly. the magnetic field drops significantly near the surface and does not cross the $B_{\rm eq}/B_{0}$ line. This means that, unlike the earlier work with imposed vertical fields (Brandenburg et al., 2014), the field in the vertical flux tubes never exceeds $B_{\rm eq}$.

The magnetic field strengths of the flux concentrations are obviously much weaker than what is expected for the Sun. More surprising is perhaps the fact that they are also much weaker than in the earlier isothermal models. For $200\,{\rm G}$, the ratio $B_{\max}/B_{0}$ reaches 1.1, while for $500\,{\rm G}$, it reaches 0.94. In the isothermal case, this value could easily reach 50. We can also observe that, when we increase the external field, $B_{\max}$ becomes smaller.

3.5 Effective magnetic pressure

In Fig. 7 we plot the normalized effective magnetic pressure,

[TABLE]

versus $z$. We compare this with Fig. 9 of Losada et al. (2014), which was a polytropic run with $\gamma=5/3$. In the present work, the values of ${\cal P}_{\rm eff}$ are about ten times larger than in the earlier polytropic models. This is probably related to the rather large values of $q_{\rm p0}$ and $\beta_{\star}$. However, the shapes of the curves are similar in those two models. Our values of the relative strength of the imposed field are similar: for $B_{0}=200\,{\rm G}$ we have $B/B_{\rm eq}=0.05$, which is comparable to the value of Losada et al. (2014). Our value of $100\,{\rm G}$ corresponds to their ratio 0.01, while $500\,{\rm G}$ corresponds to 0.07. The results change slightly when replacing $\beta_{0}\equiv B_{0}/B_{\rm eq}$ by the value for the actual magnetic field $\beta=|\bm{B}|/B_{\rm eq}$. Furthermore, the changes in the effective magnetic pressure caused by the induced magnetic field are rather strong; see the dashed lines in Fig. 7. We also see regular variations in the vertical direction, which are associated with corresponding (time-dependent) extrema in the actual magnetic field.

3.6 Dependence on $\kappa_{0}$

Increasing $\kappa_{0}$ means decreasing the radiative diffusive in the deeper parts, which tends to let NEMPI appear sooner and grow faster. It also reduces the temperature near the top of the surface and therefore also the density scale height $H_{\rho 0}$.

To see whether radiation has a noticeable effect on NEMPI, we compare in Fig. 8 vertical profiles of the relative magnetic and temperature fluctuations $\langle\delta\ln\bm{B}^{2}\rangle^{1/2}$ and $\langle\delta\ln T\rangle$ for runs with different values of $\kappa_{0}=2\times 10^{7}$, $5\times 10^{7}$, and $10^{8}\,{\rm Mm}^{-1}\,{\rm cm}^{3}\,{\rm g}^{-1}$. The effect is surprisingly small. The magnetic fluctuations are of the order of unity (and somewhat larger for $\kappa_{0}=2\times 10^{7}\,{\rm Mm}^{-1}\,{\rm cm}^{3}\,{\rm g}^{-1}$), while the relative temperature fluctuations are at most $5\times 10^{-4}$.

The vigor of the temporal variation of the field increases considerably as we increase $\kappa_{0}$, even though the relative strength of the variations and the effect on the temperature remain comparable.

3.7 Dependence on $q_{\rm p0}$ and $\beta_{\star}$

When the value of $q_{\rm p0}$ is below 250, keeping $\beta_{\rm p}=0.05$ fixed, so $\beta_{\star}=079$, NEMPI is found to be no longer excited and thus no magnetic structures are created. This remains true even when we increase $\kappa_{0}$ to $5\times 10^{7}\,{\rm Mm}^{-1}\,{\rm cm}^{3}\,{\rm g}^{-1}$, which is generally more favorable to the onset of NEMPI. This may indicate that there is a threshold for $\beta_{\star}$ for the excitation of NEMPI in the presence of radiation, which would be somewhere between 0.7 and 0.8.

4 Comparison with simpler models

To trace the origin of the difference to earlier results, we compare with models without radiative transfer. The next closest to those fully radiative models is that described in Sect. 2.4, in which the dynamics is optically thick, but a radiative boundary condition (15) is adopted at the top. The height where this condition is applied is $z=d$, which corresponds to the position where $\tau=1$ in the fully radiative model; see Fig. 1. This is at $z=d=4.3\,{\rm Mm}$, where the mean-free path is $\ell=0.14\,{\rm Mm}$, so structures that are smaller than that experience reduced radiative heat exchange with the surroundings in the fully radiative model, but not in the optically thick treatment.

Another type of simplified model is one where $\overline{s}={\rm const}{}$ in space and time. This is a strictly isentropic case, where Eq. (3) is ignored. Other than that, it has the same height and density stratification as both the optically thick model and the fully radiative one.

4.1 Optically thick case

To shed some light on the occurrence of small horizontal length scales of NEMPI in our radiative transfer models, we now compare with the optically thick approximation discussed in Sect. 2.4. The result is shown in Fig. 9 for a model that is comparable to Run B with $\kappa_{0}=10^{7}$. In that case, Eqs. (16)–(17) yield $d=4.3\,{\rm Mm}$, $T_{1}=4998\,{\rm K}$, and ${\cal K}=2.7\times 10^{-4}$. It turns out that structures now develop at $z\approx 4\,{\rm Mm}$, which is close to the top of the domain; see Fig. 9. With radiative transfer, by comparison, structures typically develop deeper down at $z\approx 3\,{\rm Mm}$. However, the structures still have very small length scales comparable to those in the models with radiative transfer. By comparing with Fig. 2 is is evident that in the models with optically thin radiative transfer, the formation of structures at $z\approx 4\,{\rm Mm}$ appears to be suppressed. The mean free path is only about $\ell=0.14\,{\rm Mm}$ for our structures with $k/k_{1}=4$; see the red dashed line in Fig. 5a. This is rather small and can therefore not be an explanation for the suppression of structures in the models with optically thin radiative transfer. There is, however, another difference between the models with optically thin radiative transfer and the optically thick approximation that does not have to do with NEMPI. All models with optically thick radiative transfer have a stably stratified layer at the top, where the entropy increases with height. Therefore, a downdraft pulls with it high entropy material, contrary to the case with a radiative boundary condition at $z=d$, where downdrafts always have low entropy. This difference was already noted by Barekat & Brandenburg (2014). It explains why NEMPI does not develop near the $\tau=1$ surface at $z=4.3\,{\rm Mm}$ in the optically thin radiative transfer model. However, it does not explain the small size of NEMPI structures. We should also point out here that, in the optically thick model, NEMPI is no longer oscillatory.

4.2 Isentropic case

In Fig. 10, we show the same model as in Fig. 9, but now with fixed mean specific entropy, so $\overline{s}={\rm const}{}$, i.e., Eq. (3) is not solved. This means that the negative buoyancy is just the result of the negative effective magnetic pressure, without any influence from changes in specific entropy and temperature. By contrast, when temperature and entropy are allowed to change, this can either enhance or diminish the effect of NEMPI. The answer discussed below is not completely straightforward.

In a stratified layer, a downdraft, even if it is initiated by NEMPI (instead of thermal buoyancy, for example), will always be compressed, so its density increases. This leads to adiabatic heating, and the corresponding radiation causes a loss of entropy, so those structures become even more negatively buoyant. This happens most efficiently at the scale of the photon mean free path or at the radiative diffusion scale. Both scale are rather small and this might explain the observed tendency for developing small structures in our model. At the same time, however, those small length scales also make NEMPI less efficient. In this sense, radiation both promotes NEMPI by enhancing buoyancy effects (both negative and positive ones), but it also counteracts NEMPI, because it operates on progressively smaller length scales.

5 Conclusions

We have presented here the first calculations of NEMPI with radiation. Within the limitations of our simplified model, NEMPI would not have been excited had we chosen the previously determined control parameters for the negative effective magnetic pressure effect, i.e., $\beta_{\star}$ and $\beta_{\rm p}$. By using a nearly three times larger value of $\beta_{\star}$, we were able to study the reason behind this. It turned out that in our model with radiation, the horizontal wavelength of the instability is dramatically decreased. As a consequence, turbulent and radiative diffusion have much stronger effects, suppressing therefore the instability. Nevertheless, even with a strongly enhanced value of $\beta_{\star}$, the resulting magnetic structures are still far too weak to form sunspots.

We found for the first time that NEMPI can display oscillatory behavior during the linear phase of the instability. These oscillations are associated with travelling waves moving upward with a speed of $0.2\,{\rm km}\,{\rm s}^{-1}$. The oscillations have a period of about $4$–$9\,{\rm ks}$ in the volume-integrated velocity, but since the period of the actual (signed) magnetic field is twice as long, so the recurrence time of pronounced downward flows is $8$–$18\,{\rm ks}$.

We do not yet know enough about the nature of the oscillations and whether they could also exist in reality. To address this question further, we have to focus on the limitations associated with the small horizontal length scales of NEMPI in the presence of radiation. Given that the oscillations occur only in the presence of a stably stratified layer above, it is possible that they are related to buoyancy oscillations in a thin upper radiative layer, where the stratification is sufficiently stable, while still being coupled to NEMPI in the deeper layers through suction along magnetic field lines.

The treatment of turbulent magnetic diffusion as a multiplicative factor in front of a Laplacian diffusion operator becomes invalid on small length scales, so the actual diffusion will be smaller; see Brandenburg et al. (2008). It is also possible that the opacity is still not large enough, and therefore the radiative diffusivity is too large. This is another unrealistic limitation of our present model. On the other hand, in the deeper layers, the radiative diffusivity is already now smaller than the turbulent magnetic diffusivity. One would therefore not have expected this to be the limiting factor. Most important is perhaps the limitation associated with the neglect of turbulent convection in the deeper parts. Convection would imply the presence of a strongly negative entropy gradient just below the surface. Therefore, the stabilizing effect from the top layers encountered in the present model would be absent. However, NEMPI would still lead to small length scales, except that now turbulent convection leads to an effective thermal diffusivity that is much larger than the radiative one. Moreover, the transition between a radiative surface above and strong turbulence with small-scale convection beneath the surface would be very abrupt. Given that NEMPI is most effective for large scale separation (small-scale turbulence) and the stratification is strongest near the surface, it might still be a viable alternative for the formation of sunspots. Extending our model by including convection in parameterized form would therefore be a first task to be addressed in a follow-up investigation.

Ultimately, the aim is to model the formation of sunspots, where convective heat transport is either suppressed by the magnetic field (Biermann, 1941) or the cooling enhanced (Parker, 1974a). The former effect may lead to its own instability, which was modelled by Kitchatinov & Mazur (2000) using a mean-field approach. This instability could be strengthened further by the effects of ionization and would therefore be another urgent target for subsequent investigations.

Acknowledgements.

We thank the referee for useful comments and Sacha Brun for support and encouragement. Support through the NSF Astrophysics and Astronomy Grant Program (grant 1615100) and the Research Council of Norway (FRINATEK grant 231444) are gratefully acknowledged. 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. This work utilized the Janus supercomputer, which is supported by the National Science Foundation (award number CNS-0821794), the University of Colorado Boulder, the University of Colorado Denver, and the National Center for Atmospheric Research. The Janus supercomputer is operated by the University of Colorado Boulder.