Issue 
A&A
Volume 621, January 2019



Article Number  A61  
Number of page(s)  17  
Section  The Sun  
DOI  https://doi.org/10.1051/00046361/201833018  
Published online  09 January 2019 
Magnetic bipoles in rotating turbulence with coronal envelope
^{1}
Nordita, KTH Royal Institute of Technology and Stockholm University, 10691 Stockholm, Sweden
email: illa.rivero.losada@gmail.com
^{2}
Department of Astronomy, AlbaNova University Center, Stockholm University, 10691 Stockholm, Sweden
^{3}
Nordic Optical Telescope, La Palma, Canary Islands, Spain
^{4}
MaxPlanckInstitut für Sonnensystemforschung, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
^{5}
JILA and Department of Astrophysical and Planetary Sciences, University of Colorado, Box 440, Boulder, CO 80303, USA
^{6}
Laboratory for Atmospheric and Space Physics, University of Colorado, Box 590, Boulder, CO 80303, USA
^{7}
Department of Mechanical Engineering, BenGurion University of the Negev, POB 653, BeerSheva 84105, Israel
Received:
13
March
2018
Accepted:
16
September
2018
Context. The formation mechanism of sunspots and starspots is not yet fully understood. It is a major open problem in astrophysics.
Aims. Magnetic flux concentrations can be produced by the negative effective magnetic pressure instability (NEMPI). This instability is strongly suppressed by rotation. However, the presence of an outer coronal envelope was previously found to strengthen the flux concentrations and make them more prominent. It also allows for the formation of bipolar regions (BRs). We aim to understand the important issue of whether the presence of an outer coronal envelope also changes the excitation conditions and the rotational dependence of NEMPI.
Methods. We have used direct numerical simulations and meanfield simulations. We adopted a simple twolayer model of turbulence that mimics the jump between the convective turbulent and coronal layers below and above the surface of a star, respectively. The computational domain is Cartesian and located at a certain latitude of a rotating sphere. We investigated the effects of rotation on NEMPI by changing the Coriolis number, the latitude, the strengths of the imposed magnetic field, and the box resolution.
Results. Rotation has a strong impact on the process of BR formation. Even rather slow rotation is found to suppress BR formation. However, increasing the imposed magnetic field strength also makes the structures stronger and alleviates the rotational suppression somewhat. The presence of a coronal layer itself does not significantly reduce the effects of rotational suppression.
Key words: magnetohydrodynamics (MHD) / turbulence / dynamo / Sun: magnetic fields / Sun: rotation / Sun: activity
© ESO 2019
1. Introduction
The solar dynamo operates in hydromagnetic turbulence in the presence of strong stratification – especially near the surface. The stratification can lead to a secondary instability, in addition to the primary dynamo instability, and can concentrate the field further into spots. This instability was discovered by Kleeorin et al. (1989, 1990) and applied to explain sunspot formation and other hydromagnetic processes in the Sun (Kleeorin et al. 1993, 1996; Kleeorin & Rogachevskii 1994; Rogachevskii & Kleeorin 2007). In the last ten years, direct numerical simulations (DNS) have demonstrated that this instability – the negative effective magnetic pressure instability (NEMPI) – is able to concentrate magnetic fields in different physical environments.
Since its first detection in a stably stratified isothermal setup with a weak horizontal (Brandenburg et al. 2011) and a weak vertical (Brandenburg et al. 2013) magnetic field, NEMPI demonstrated its operation in polytropic stratification (Losada et al. 2014) and turbulent convection (Käpylä et al. 2012, 2016), in the presence of weak rotation (Losada et al. 2012, 2013), and with a coronal envelope (Warnecke et al. 2013a, 2016a). In this context, a “weak” magnetic field is one that is of subequipartition strength, but already dynamically important. Likewise, weak rotation means that the angular velocity is small compared with the inverse correlation time of the turbulence, although it is already dynamically important. The coronal envelope is crucial for producing a bipolar region (BR) in the domain. Here we define BRs as structures that are much larger than the size of individual turbulent eddies. The emergence of opposite polarities is a consequence of zero total vertical flux across a horizontal surface. This implies that regions with weak largescale magnetic fields are separated by regions with strong fields of opposite magnetic polarity.
In principle, spot structures could also be tripolar or quadrupolar (for example AR 11158), but those cases are rare. Interestingly, these tri and quadrupolar spot structures allow one to determine the magnetic helicity in the space above by only using the surface magnetic field (Bourdin & Brandenburg 2018). This property may help tying the nature of dynamogenerated subsurface magnetic fields to observations. A similarly useful tool is helioseismology, which may allow for a detection of subsurface magnetic fields in the days prior to active regions formation; see Singh et al. (2016) for such work providing evidence for a gradual buildup of active regions rather than a sudden buoyant emergence from deep down.
The coronal envelope allows the orientation of a weak imposed horizontal magnetic field to change locally due to the interface between the turbulence zone and the coronal envelope so that it can attain a vertical component. On the other hand, rotation plays against the instability (i.e., NEMPI), which cannot survive above a certain critical rotation rate. This critical rate is found to be surprisingly small (Losada et al. 2012, 2013).
Other types of surface magnetic flux concentrations develop in a number of different circumstances, all of which share the presence of a strong density stratification. There is first of all the phenomenon of magnetic flux segregation into weakly convecting magnetic islands within nearly fieldfree convecting regions (Tao et al. 1998). This has also been seen in several recent high resolution and high aspect ratio simulations (Käpylä et al. 2013, 2016) and perhaps also in those of Masada & Sano (2016). This process may explain the formation of flux concentrations seen in the simulation of Stein & Nordlund (2012), in which an unstructured magnetic field of 1 kG is allowed to enter the computational domain at the bottom. These simulations include realistic surface physics in a domain 96 Mm wide and 20 Mm deep, so, again, the aspect ratio is large and there is significant scale separation. Another approach is to let a flux tube rise from the bottom of the computational domain to simulate flux emergence at the surface (Fan et al. 1993; Fan 2001; Archontis et al. 2004, 2005; Fournier et al. 2017); see Fan (2009) for a review. Similar results, but with realistic surface physics, have been obtained by Cheung et al. (2010) and Rempel & Cheung (2014), who let a semitorus of magnetic field advect through the lower boundary. Their simulations show that the magnetic field is able to rise through the top 16 Mm of the convection zone to form spots; see the reviews by Cheung & Isobe (2014) and Schmieder et al. (2014). In their simulations, however, flux emergence is significantly faster than in the real Sun (Birch et al. 2016). Recently, Chen et al. (2017) have been able to reproduce a complex sunspot emergence using a modified magnetic flux bundle from the dynamo simulation of Fan & Fang (2014) in a spherical shell and inserting it into a setup similar to that of Rempel & Cheung (2014).
The production of magnetic flux concentrations in convection could also be related to the magnetic suppression of the convective heat flux, which, again, could lead to a largescale instability (Kitchatinov & Mazur 2000). On the other hand, Kitiashvili et al. (2010) explain the formation of magnetic flux concentrations in their radiationhydromagnetic simulations as being confined by the random vorticity associated with convective downdrafts. However, the process seen in their simulations seems to be similar to flux concentrations found in the downdraft of the axisymmetric simulations of Galloway & Weiss (1981). Likewise, there is a process known as convective collapse (Parker 1978), which can lead to a temporary concentration of field from a weaker, less concentrated state to a more concentrated collapsed state, for example, from 1270 G to 1650 G in the specific calculations of Spruit (1979). However, the collapsed state is not in thermal equilibrium, so the system will slowly return to an uncollapsed state. This effect may be related to the ionization physics, which can strongly enhance the resulting concentrations (Bhat & Brandenburg 2016).
Strong magnetic flux concentrations also appear in simulations where a largescale dynamo is responsible for generating magnetic field. The dynamo arises as a result of helically driven turbulence in the lower part of the domain, while in the upper part turbulence is nonhelically driven. In the upper part, the magnetic field forms strongly concentrated BRs. This was seen in DNS in both Cartesian domains (Mitra et al. 2014; Jabbari et al. 2016, 2017) and spherical shells (Jabbari et al. 2015).
The relation to NEMPI is unclear in some of these cases, because NEMPI can be excited when the magnetic field is below the equipartition value of the turbulence. A negative effective magnetic pressure is possible in somewhat deeper layers and at intermediate times, which may be important in initializing the formation of magnetic flux concentrations. In particular, it would lead to downward suction along vertical magnetic field lines that creates an underpressure in the upper parts and results in an inflow. The latter causes further concentrations in the upper parts. This was clearly seen in the axisymmetric meanfield simulations (MFS) of Brandenburg et al. (2014); see Brandenburg et al. (2016) and Losada et al. (2017) for recent reviews.
In the present work, we have considered a setup similar to that of Warnecke et al. (2013a, 2016a). There, turbulence of an isothermal gas is forced in the lower part of a horizontally periodic domain, while the upper part is left unforced and subject to the response from the dynamics of the lower part. This approach has been used to study the effect of a coronal envelope on the dynamo (e.g., Warnecke & Brandenburg 2010; Warnecke et al. 2011, 2013b, 2016b) and the formation of coronal ejections (Warnecke et al. 2012a,b). In these simulations, the simple treatment of the coronal envelope does not allow for a low plasma β_{c} as in the solar corona, where β_{c} is the ratio of magnetic pressure to gas pressure. However, in the solar corona too, the value of β_{c} is not extremely small (e.g., Peter et al. 2015) and plasma flows can play an important role for the formation of loop structures (Warnecke et al. 2017).
We have included the Coriolis force to examine the effects of rotation in the presence of our simplified corona. Then, we studied whether this facilitates the development and detectability of NEMPI, and whether it changes the critical growth rate above which NEMPI is suppressed. We have also studied the dependence on latitude, as well as the dependence on the numerical resolution. Finally, we compared our solutions with corresponding MFS, where a prescribed effective (meanfield) magnetic pressure operates only beneath the surface, but not in the coronal layer.
2. The model
2.1. DNS
We used the same twolayer model as Warnecke et al. (2013a, 2016a). We considered a Cartesian domain with forced turbulence in the lower part (referred to as turbulent layer), and a more quiescent upper part (referred to as coronal envelope). We further adopted an isothermal equation of state and solved the equations for the velocity U, the magnetic vector potential A, and the density ρ. We adopted units for the magnetic field such that the vacuum permeability is unity. Here we extend this model by including the presence of rotation with an angular velocity Ω,
where D/Dt = ∂/∂t + U ⋅ ∇ is the advective derivative, B = B_{0} + ∇ × A is the magnetic field, B_{0} = (0, B_{0}, 0) is a weak imposed uniform field in the y direction, J = ∇ × B is the current density, ρ^{−1} so that it reads F_{ν} = ρ^{−1} ∇ ⋅ (2νρ𝖲) is the viscous force, ν is the kinematic viscosity, η is the magnetic diffusivity, g is the gravitational acceleration, is the traceless rateofstrain tensor, and f is a forcing function that consists of random, whiteintime, plane, nonpolarized waves within a certain narrow interval around an average wavenumber k_{f}. It is modulated by a profile function Θ_{w}(z),
that ensures a smooth transition between unity in the lower layer and zero in the upper layer. Here w is the width of the transition. The angular velocity vector Ω is quantified by its modulus Ω and colatitude θ, such that
Following Warnecke et al. (2013a), the domain used in the DNS is L_{h} × L_{h} × L_{z}, where L_{h} = 2π and L_{z} = 3π with −π ≤ z ≤ 2π with −π ≤ z ≤ 0 for the turbulent layer and 0 ≤ z ≤ 2π for the coronal envelope. This defines the base horizontal wavenumber k_{1} = 2π/L_{h}, which is set to unity in our model. Our Cartesian coordinate system (x, y, z) corresponds to a local representation of a point on a sphere mapped to spherical coordinates (r, θ, ϕ)→(z, x, y), where r is radius, θ is colatitude, and ϕ is longitude. Similar to earlier work (Kemel et al. 2012a, 2013), in all cases we used k_{f} = 30 k_{1} and ν = 10^{−4} c_{s}/k_{1} with the sound speed c_{s}. The normalized gravity is given by , which is just slightly below the value of 1.2 that is found to maximize the amplification of magnetic field concentrations (Warnecke et al. 2016a). Here H_{ρ} is the density scale height. As in most of our earlier work, we use k_{1} H_{ρ} = 1, so the vertical density contrast is exp(L_{z}/H_{ρ}) = exp3π ≈ 12, 000 in the total box and exp(L_{z}/H_{ρ}) = expπ ≈ 23 in the turbulent layer. We set the width of the profile functions in the DNS and MFS to w = 0.05/k_{1}, as in most of the cases of Warnecke et al. (2016a).
For the rms velocity, we used the averaged value in the turbulent layer defined as: . We normalized the magnetic field by its equipartition value, , using either the z dependent value of B_{eq} or the value at the surface at z = 0, i.e., B_{eq0} ≡ B_{eq}(z = 0).
Our simulations are characterized by the magnetic and fluid Reynolds numbers,
respectively, and the Coriolis number
where τ = 1/u_{rms} k_{f} is the eddy turnover time, and the colatitude θ of our domain is positioned on the sphere. In the following, we use Re_{M} ≈ 14 and Re ≈ 29, so the magnetic Prandtl number is Pr_{M} = ν/η = 0.5. These values of the magnetic and fluid Reynolds numbers are based on the forcing wavenumber, which is rather high (k_{f}/k_{1} = 30). Thus, the values of Re_{M} and Re based on the wavenumber of the domain would be 420 and 870, respectively, and those based on the size of the domain, which are larger by another factor of 2π, would be 2640 and 5470, respectively. The definitions of these Reynolds numbers must therefore be kept in mind when comparing with other work. In our definition, the magnetic Reynolds number required for the effective magnetic pressure to be negative must be larger than a critical value of about three (Brandenburg et al. 2011, 2012; Käpylä et al. 2012; Warnecke et al. 2016a). For small Re_{M}, the effective magnetic pressure can only be positive (Rüdiger et al. 2012; Brandenburg et al. 2012). In our work, we chose Re_{M} to be about ten times supercritical. With the choice of Re ≈ 29 and Pr_{M} = 0.5, we exclude that the excitation of a smallscale dynamo influences our results. As shown in Warnecke et al. (2016a), Pr_{M} ≥ 1 is needed to excite a smallscale dynamo in this setup.
In this work, time is often expressed in units of the turbulent diffusive time, , where η_{t0} = u_{rms}/3k_{f} is an estimate of the turbulent magnetic diffusivity. As in our earlier work (Warnecke et al. 2013a, 2016a), we used the Fourier filtered magnetic and velocity fields as diagnostics for characterizing largescale properties of the solutions. Our Fourier filtered fields are denoted by an overbar and the superscript fil, i.e., for the vertical magnetic field. This includes contributions with horizontal wavenumbers below k_{f}/2. This filtering wavenumber is the same as that used in Warnecke et al. (2013a), but the cutoff wavenumber is three times larger than the filter value k_{f}/6 used by Brandenburg et al. (2013) and Warnecke et al. (2016a). This is appropriate here because, owing to the nature of our BRs, where the spots tend to be close together, they would not be well captured when the averaging scale is too large or the filtering wavenumber too small. The corresponding spectral magnetic energy contained in the vertical magnetic field, B_{z}, is and obeys . Of particular interest is the energy per logarithmic wavenumber interval, , which we usually evaluate at k_{*}/k_{1} = 2, where the energy reaches a maximum. (The factor 2 in front of compensates for the 1/2 factor in the definition of the energy.) For the velocity, however, we find that k_{f}/6 is the appropriate filtering wavenumber. Therefore, we filter the velocities on a larger scale than the magnetic field.
We computed growth rates and magnetic energies as averages over a certain time interval. We calculate error bars as the largest departure from any one third of the full time interval used for computing the average. In some cases, those error estimates were themselves unreliable. In such exceptional cases we have replaced it by the average error for other similar simulations.
We used resolutions between 192 × 192 × 384 and 1152 × 1152 × 2304 meshpoints in the x, y, and z directions, respectively. We adopted periodic boundary conditions in the xy plane, a stressfree perfect conductor condition on the bottom boundary, and a stressfree vertical field condition on the top boundary.
2.2. MFS
In this section, we state the relevant equations for the meanfield description of NEMPI in a system with coronal envelope in the presence of rotation. The pertinent equations have been obtained by Kleeorin et al. (1989, 1990, 1996) and Kleeorin & Rogachevskii (1994) through ensemble averaging. In practical applications, these averages should be replaced by spatial averages, but their precise nature depends on the physical circumstances and could be planar (e.g., horizontal) or azimuthal (e.g., around a flux tube), or some kind of spatial smoothing, although those obey the Reynolds rules only approximately. In the present case with inclined stratification, rotation vectors, and mean magnetic field vectors, we expect the formation of bipolar structures that cannot be described by simple planar or azimuthal averages. The only meaningful average is a smoothing operation that could preserve such structures. In the following, we denote the dependent variables of our MFS by overbars. For the analysis of our DNS, on the other hand, we sometimes used Fourier filtering and sometimes horizontal averaging and also denoted them with an overbar. In those cases, a corresponding comment is made, and in the particular case of Fourier filtering, the variable is additionally denoted by the superscript “fil”.
In the MFS, the equations for the mean velocity U̅, mean vector potential A̅, and mean density , are given by
where D̅/Dt = ∂/∂t + U̅ · ∇ is the advective derivative based on U̅, η_{t} is turbulent magnetic diffusivity,
is the total (turbulent plus microscopic) viscous force with ν_{t} being the turbulent viscosity, is the traceless rateofstrain tensor of the mean flow and, as in the DNS, we adopt units for the mean magnetic field such that the vacuum permeability is unity. The effective Lorentz force, , which takes into account the turbulence contributions, that is, the effective magnetic pressure (Kleeorin et al. 1996; Rogachevskii & Kleeorin 2007; Brandenburg et al. 2016) and an anisotropic contribution resulting from gravitational stratification, is given by
where q_{p} and q_{g} are functions that have previously been determined from DNS (Brandenburg et al. 2012). Warnecke et al. (2016a) found q_{g} to be negative for weak and moderate stratification (we note xthat the abscissa of their Fig. 6 shows and not , as was incorrectly written). Thus, there is the possibility of partial cancelation, which we model here by assuming q_{p} and q_{g} to have the same profile with
(13) Warnecke et al. (2016a) determined q_{p}β^{2} = −q_{g}β^{2} = −0.002, resulting in a_{g} = −1 for the same stratification as in this work (). We note here that the values of q_{p} and q_{g} are the result of averaging in time and space, so locally the values can be different, and therefore they do not need to cancel out locally. Furthermore, the error estimate by the spread is comparable to the averaged value, see Fig. 6 of Warnecke et al. (2016a). We modeled q_{p} and q_{g} as the product of a part that depends only on and a profile function Θ_{w}(z). The latter function varies only along the z direction and mimics the effects of the coronal layer, using the same error function as in Eq. (4), i.e.,
where
Here is a parameter that can be used alternatively to q_{p0} and has the advantage that the growth rate of NEMPI is predicted to be proportional to it (Kemel et al. 2013). In the following, we mainly use the parameters found by Losada et al. (2013), namely q_{p0} = 32 and β_{p} = 0.058, which corresponds to β_{∗} = 0.33. On one occasion, we also used another parameter combination that will be motivated by our results presented in Sect. 3.7 below. In hindsight, it might have been more physical to use in Eq. (14), but we know from experience that this would make a hardly noticeable difference.
We note that, while in both the DNS and the MFS, the coronal envelope is modeled with the same profile function Θ_{w}(z), in the MFS, it appears underneath the gradient (inside the effective Lorentz force), while in the DNS it does not. In the MFS, this leads to an additional term involving the derivative of Θ_{w}(z), which is a gaussian function. In one of the cases reported below, we checked that the presence of this term causes a very minor difference in the growth rates.
In addition to adopting the parameterization given by Eq. (14), the effects of turbulence are modeled in terms of turbulent viscosity ν_{t} and turbulent magnetic diffusivity η_{t}. Both can be expressed in terms of u_{rms} and k_{f}, whose values are known from the DNS and give ν_{t} = η_{t} ≈ 10^{−3} c_{s}H_{ρ}. Thus, 3η_{t}/η = Re_{M} and 3ν_{t}/ν = Re. Among the range of other possible meanfield effects, we have included here only the negative effective magnetic pressure functions q_{p} and q_{g}.
As already noted by Losada et al. (2012), the usual Coriolis number is not the relevant quantity characterizing the relative importance of rotation on NEMPI. A more meaningful quantity is the ratio of 2Ω/λ_{∗0}, where λ_{∗0} is the nominal value of the growth rate
which was found to be comparable to the growth rate of NEMPI (Losada et al. 2012). Therefore, the Coriolis number given by Eq. (7), can be rewritten as
With u_{rms} ≈ 0.1c_{s} and k_{1}H_{ρ} = 1, this means that λ_{∗0} = 0.033c_{s}k_{1}; see Eq. (16). This value has also been used to characterize the rotation rate of the DNS, which are well characterized by the parameter β_{∗} = 0.33; see Brandenburg et al. (2012). Therefore, 2Ω/λ_{∗0} is about 100 times larger than Co. We note also that in the MFS, (Jabbari et al. 2014), which results in the same estimate. The actual growth rates obtained from our DNS and MFS will be normalized either also by λ_{∗0} or by .
Both DNS and MFS simulations are done with the PENCIL CODE^{1}. It uses sixth order accurate finite differences and a thirdorder timestepping scheme. It comes with a special meanfield module that can be invoked for the MFS.
3. DNS results
3.1. Numerical resolution
Since NEMPI is a meanfield instability, which relies on smallscale turbulence for developing largescale structures, we begin by demonstrating the effects of changing the resolution on the formation of magnetic field concentrations. The growth rates are shown in Fig. 1 as a function of resolution. The corresponding simulations are listed in Table 1 where we have always used the same domain size of (2π)^{2} × 3π, but with different numbers of meshpoints. We also quote the mesh Reynolds number, Re_{grid} = u_{rms}/νk_{Ny}, where k_{Ny} = π/δx is the Nyquist wavenumber and δx is the mesh spacing. In all cases, this number is well below unity. Sometimes the quantity u_{rms}δx/ν is used in the literature; it is simply 2π times larger than our Re_{grid}.
Fig. 1.
Growth rate λ (black line), the peak spectral magnetic field (red), and the maximum of the Fourierfiltered vertical magnetic field (blue) for Runs A1–A5 at different resolutions. 
Summary of run with different resolutions.
We see formation of bipolar regions (denoted by BR in the table) in all the cases, but at the lowest resolution, the growth rate of the magnetic field is significantly smaller than at all higher resolutions. Based on this, we conclude that a resolution of 384^{2} × 768 meshpoints results in a good compromise between accuracy and computational cost. Therefore, we use simulations with this resolution for the following parameter study. Because of this, Warnecke et al. (2016a) double their resolution in their followup work of Warnecke et al. (2013a). Qualitatively, the coherence increases with increasing resolution; see Fig. 2. From Table 1, the structures become strongest at a resolution of 768^{2} × 1536 meshpoints, where the normalized Fourierfiltered vertical magnetic field at the surface, , reaches a value of 0.54 corresponding to only 0.4 at both lower and higher resolutions.
Fig. 2.
Vertical magnetic field B_{z} at the surface (z = 0) for runs with different resolution (Runs A1–A3 and A5) at the time when the structures are strongest. The resolution is indicated by the number of grid points at the top of each panel. 
3.2. Growth of BRs
We now discuss the main properties of BRs. Largescale BRs form during the first one or two turbulentdiffusive times. They are referred to as largescale structures, because their size is that of many turbulent eddies (about 2π/k_{f} in the horizontal plane). To average over these turbulent eddies and still resolve the largescale structure of the BRs, we Fourierfilter the magnetic field at the surface. To investigate the growth of these structures we then plot the maximum of the vertical Fourierfiltered magnetic field over time. In Fig. 3, which corresponds to Run A2, we clearly see an exponential growth with growth rate . A similar growth has been seen before in numerical experiments both with a horizontally imposed field (Brandenburg et al. 2011; Kemel et al. 2012a) and a vertical one (Brandenburg et al. 2013). A similar value has also been determined with the same setup, where BRs form (Warnecke et al. 2016a). Such exponential growth is suggestive of NEMPI.
Fig. 3.
Evolution of vs. t/τ_{td} for Run A2 showing exponential growth with growth rate (black solid line), compared with that of (red dashed line) and (blue solid line) at z = 0. We also show the unity line corresponding to B_{eq} (orange dotted line). 
In Fig. 3, we also show for comparison the evolution of the maximum of the surface vertical magnetic field , that is, not the filtered value. Its value is always close to the local equipartition field strength. No exponential growth phase can be seen in the temporal evolution of . Likewise, the rms value of surface vertical magnetic field is about 20 times smaller and shows no exponential growth.
The formation of BRs in our simulations are associated with largescale flows. As for the magnetic field, we performed Fourier filtering to averaged over the turbulent scales. In Fig. 4 we show these Fourierfiltered largescale flows that includes only wavenumbers below k_{f}/6. The two panels are for Runs A4 and A5 with two different resolutions, but otherwise the same as Run A2, and similar times. In both cases there are BRs, but in one case the two spots are more separated from each other. Nevertheless, in all cases the BRs are surrounded by a largescale inflow with relative rms value –0.18, where . Its maximum value is around 0.4. Given that the scale of the turbulent motions is much smaller than the scale of the spots, and that there is otherwise no mechanism producing largescale flow perturbations, the inflow can only be a consequence of the magnetic field itself, and not the other way around. This is indeed also what one expects for NEMPI.
Fig. 4.
Vectors of the largescale (Fourierfiltered) flow superimposed on a color scale representation of B_{z} at the surface at t/τ_{td} = 0.91 for Run A4 (upper panel) and t/τ_{td} = 1.12 for Run A5 (lower panel). Note that the coordinate system has been rotated by 90°, so y points to the right and x points downward. 
Another diagnostics for the formation of magnetic structures in our simulations are magnetic power spectra taken at the surface. In Fig. 5 we can see the time evolution of such power spectra. At early times when there are no structures yet, the spectrum peaks at the energy injection wavenumber, k_{f} = 30 k_{1}. As time evolves and structures start forming near the surface, magnetic energy is transported toward smaller wavenumbers. When we see the BRs forming at the surface, the magnetic power spectrum peaks at k = k_{∗} ≈ 2 k_{1}, and the amplitude at this wavenumber decreases until the structures disappear. The strength of the magnetic surface structures at this wavenumber is characterized by the value of at k = k_{∗} = 2 k_{1}; see Sect. 2.1. The resulting values are listed in Tables 1–3. In addition, we also give the spectral values averaged over the first four wavenumbers from k/k_{1} = 1 to 4. The averaging helps reducing the sensitivity to discretization noise that arises from looking at just one wavenumber k_{∗}. In the tables, we also judge BR formation as Yes, no, or weak. These attributes are based on the qualitative assessment of images of B_{z}.
Fig. 5.
Power spectra of B_{z} at the surface (z = 0) for Run A2 at several times around the maximum growth for zero rotation. The inset shows the time dependence (time in units of τ_{td}) for k = k_{∗} = 2 k_{1} and for the values averaged between k/k_{1} = 1 and 4. 
The energy transfer to larger scales is reminiscent of an inverse cascade. Brandenburg et al. (2014) have speculated that such a cascade might be a consequence of the conservation of cross helicity, , where the overbar denotes horizontal averaging, are the magnetic fluctuations and are the velocity fluctuations; see Rüdiger et al. (2011). (For horizontal averages, we usually have .) They studied the production of as a result of a mean magnetic field along the direction of gravity, so there exists a pseudoscalar g·B_{0} that has the same symmetry properties as and is also odd in the magnetic field. A similar result was also obtained by Kleeorin et al. (2003); see their Eq. (11) where they considered inhomogeneous turbulence. Furthermore, recent work of Zhang & Brandenburg (2018) has shown that the cross helicity spectrum shows a steep slope at large scales. This was interpreted as a potential signature of NEMPIlike effects.
3.3. Influence of rotation
The structure of the bipolar regions is strongly influenced by rotation; see Fig. 6. Faster rotation causes the magnetic flux concentrations to be weaker. In some cases, for example in Run B3 (middle panel of Fig. 6), the structure of BRs splits into three parts with one negative and two positive polarities. In our most rapidly rotating case (rightmost panel of Fig. 6), the structure appears rotated by 90° with respect to slower rotations.
Fig. 6.
Rotational dependency on the active region formation. We show the vertical magnetic field B_{z}/B_{eq} for Runs A2, B2 and F2 corresponding to Co = 0, 0.0012, 0.0076 at the time of the strongest BR. In the first panel for Co = 0, the value of θ is insignificant. 
Table 2 shows that, as the Coriolis number is increased from 0.0012 to 0.0076, BR formation gets almost entirely suppressed (see Runs B1 to F4). At θ = 0° (corresponding to the pole), BR formation tends to be slightly easier, but as the Coriolis number is increased the BRs become clearer at intermediate latitudes.
Summary of runs with different Coriolis numbers and colatitudes.
The rotational dependency of the growth rate λ for different colatitudes θ is shown in Fig. 7. For most of the colatitudes, the growth rate first increases for weak rotation and then gets reduced to around a third of the values for more rapid rotation. Even though we cannot visually detect any clear indication of BRs in the rapidly rotating runs, the growth rate of magnetic field is still positive. The magnetic field strength determined from the spectral energy (see middle panel of Fig. 7) shows actually an enhancement with rotation. Even for the rapidly rotating cases, the magnetic field strength in the flux concentrations is stronger than without any rotation. The strongest values are achieved with Co = 0.002. For most of the colatitudes, the maxima of the largescale magnetic field, as shown in the bottom panel of Fig. 7, decrease for increasing rotation, similarly as the λ (top panel). However, we see the maximum value for Co = 0.002 to 0.003.
Fig. 7.
Rotational dependency of the growth rate λ, the peak spectral magnetic field , and the maximum of the Fourierfiltered vertical magnetic field for colatitudes 0° (black line), 30° (blue line), 60° (purple line), and 90° (red line). The error bars are either the errors of the exponential fit (λ; see Table 2), or estimated as 10% of the actual value. 
The time when the Fourierfiltered magnetic field and the field of spectral energy become maximal depends on rotation, but there is no clear trend visible (see Tables 2 and 3). However there seems to be an indication that the time becomes longer, as expected for a smaller growth rate. As the Coriolis number increases, the structures become weaker and it gets more difficult to discern a clear rotation pattern.
Summary of runs with different values of imposed field and Coriolis number.
We note here that, contrary to convection, the energycarrying length scale is not influenced by rotation, because the driving scale is prescribed through the forcing function. Similarly, the kinetic energy is also only weakly influenced by rotation, so u_{rms} decreases only weakly for more rapid rotation.
3.4. Dependence on latitude
As we change the colatitude θ, the growth rates and the strengths of BRs change. This is demonstrated in Fig. 8, where we show λτ_{td}, and for Co = 0.0012 (corresponding to 2Ω/λ_{∗0} ≈ 0.15) and different values of θ. For θ = 0, which corresponds to the pole, the saturation magnetic field strength shows a maximum at t/τ_{td} ≈ 1. For larger values of θ, that is, closer to the equator, the maximum is slightly higher (about 0.8) and more longlived, for example, for 0.5 ≤ t/τ_{td} ≤ 2.5 at θ = 60°, when is above 0.6 and sometimes even 0.8.
Fig. 8.
Summary of the results for the different colatitudes showing the normalized growth rate of the Fourierfiltered zcomponent of the magnetic field (top panel) for B_{0}/B_{eq0} = 0.02, the averaged magnetic spectrum over k ≤ 4 (middle panel) and the maximum of the Fourierfiltered vertical magnetic field (bottom panel) for different values of the Coriolis number and colatitude. As elsewhere, error bars are either the errors of the exponential fit (λ; see Table 3), or estimated as 10% of the actual value. 
Bipolar structures are still fairly pronounced at θ = 60°, i.e., at 30° latitude; see Fig. 9. It is remarkable that for all values of θ, the inclination of BRs is approximately the same. The same feature is also seen in the MFS, except that for θ = 90° the structures are always aligned in the y direction. This may also be the case here, but the structures are so weak that this is hard to see. For somewhat faster rotation, when Co > 0.01, corresponding to 2Ω/λ_{∗0} ≈ 1, the structures disappear.
Fig. 9.
Vertical surface magnetic field B_{z}/B_{eq} for Runs C1–C3 and C5 with Co = 0.0012. 
3.5. Inclination of BRs
A systematic orientation of BRs was already seen in the original papers of Losada et al. (2012, 2013) both for MFS and DNS. For most of the runs, the BRs are either aligned with the imposed magnetic field or they are inclined by 45°. To compare with the Sun, we map our Cartesian coordinate system to spherical ones via (x, y, z) → (θ, ϕ, r), so the y coordinate points in the toroidal (eastward) direction and x corresponds to colatitude, which points southward. This explains the orientations of our surface visualizations where we plot −x versus y; see Fig. 9.
We usually find that rotation leads to a poleward tilt of the BR. We return to the question of the inclination angle in Sect. 4.4 where we study a similar phenomenon in MFS and show that the inclination angle is then not an artifact of the domain size. We also note that the orientation is the same in DNS and MFS. In fact, the orientation agrees with that found by Losada et al. (2012, 2013). It is interesting to note that the sense of inclination is the other way around than what is expected based on the buoyant rise of magnetic flux tubes, which gives rise to Joy’s law (Choudhuri & Gilman 1987).
The “antiJoy’s” law orientation of NEMPI structures is likely a consequence of the interaction of rotation with the concentration of flux as opposed to the expansion of flux, which is usually expected as a flux tube rises through a stratified layer. A similar phenomenon of a concentration of flux in stratified turbulence (as opposed to an expansion) has been found in turbulence driven by the magnetorotational instability, where this was argued to be the reason for an unconventional sign of the α effect (Brandenburg et al. 1995); see the more detailed discussion of Brandenburg & Campbell (1997).
3.6. Dependence on the imposed magnetic field strength
Rotation weakens the formation of structures for even smaller Coriolis numbers than those of previous studies (Losada et al. 2012, 2013). We can therefore increase the efficiency of the formation mechanism by increasing the imposed magnetic field, as was suggested by Warnecke et al. (2016a). The increase in the imposed field allows the BRs to be formed for even higher Coriolis numbers, up to the point where the magnetic field becomes so strong that the derivative becomes positive and NEMPI cannot be excited in the domain; for details, see Brandenburg et al. (2016) and the appendix of Kemel et al. (2013).
Table 3 shows simulations of domains located at the pole (θ = 0°) and gives their dependence on the magnetic field strength and angular velocity. As B_{0}/B_{eq0} is doubled from 0.07 to 0.14, BR formation becomes slightly easier: compare model G6 with model H4 (both are for Co = 0.018). Weak BR formation is only possible in model H4. For even stronger fields, however, this trend disappears. In model I4 with B_{0}/B_{eq0} = 0.3, BR formation is now very weak and for Co = 0.019, and in model J3 with B_{0}/B_{eq0} = 0.9, BR formation is impossible – even for Co = 0.014; see Table 4 for an overview.
Comparison of BR formation as a function of rotation in terms of Co and imposed magnetic field B_{0}/B_{eq0}.
The growth rate of the magnetic field shows a strong decrease for higher rotation rates. This can be compensated for to some extent by using a stronger imposed magnetic field; see the top panel of Fig. 10. For Co ≈ 0.006, the growth rates with imposed magnetic fields of B_{0}/B_{eq0} = 0.066 and B_{0}/B_{eq0} = 0.28 are indeed higher than for B_{0}/B_{eq0} = 0.026. A similar behavior can be found by looking at the magnetic field strength determined from the spectral energy (see middle panel of Fig. 10). For imposed magnetic field strengths between B_{0}/B_{eq0} = 0.066 and 0.85, the magnetic field is not much influenced by rotation and stays roughly constant above Co ≈ 0.006 at a level that is even higher than for smaller rotation. The largescale magnetic field for imposed magnetic fields larger than B_{0}/B_{eq0} = 0.06 does not show a strong rotational influence. The largescale magnetic field is higher for larger imposed magnetic field and keeps this level for a large range of rotation rates. This means that a higher imposed magnetic field can indeed prevent NEMPI from being quenched for larger rotation rates. Even with a large imposed magnetic field, rotational quenching takes place, but at larger rotation rates.
Fig. 10.
Rotational dependency of the growth rate λ, the peak spectral magnetic field , and the maximum of the Fourierfiltered vertical magnetic field for various averaged imposed magnetic field strengths B_{0}/B_{eq0} = 0.026 (blue line), 0.066 (red), 0.14 (green), 0.28 (cyan), and 0.85 (black). Error bars denote either the errors of the exponential fit (λ; see Table 2), or are estimated as 10% of the actual value. 
3.7. Effective magnetic pressure
The concentration of magnetic field in a NEMPI scenario is possible due to a turbulence effect on the effective magnetic pressure, which is the sum of nonturbulent and turbulent contributions to the largescale magnetic pressure. This effect results in a suppression of the total (hydrodynamic plus magnetic) turbulent pressure by the largescale (mean) magnetic field. This means that an increase of the mean magnetic field due to the instability will be accompanied by a decrease of the turbulent pressure and a reduction of the equipartition field strength, B_{eq}. This is because the hydrodynamic part of the total turbulent pressure, , as well as the turbulent kinetic energy density, decrease due to an increase of the turbulent magnetic energy density as well as the turbulent magnetic pressure through tangling magnetic fluctuations (Kleeorin et al. 1996; Rogachevskii & Kleeorin 2007).
Figure 11 shows the profiles of B_{eq} and B_{z} at the surface along y = 0 at t/τ_{td} = 1. This figure clearly shows that at the location of the maximum of the vertical magnetic field, the equipartition field B_{eq} is decreased. Most of this suppression comes from a local decrease in the turbulent velocity, while the local density through the structure varies only little; see the blue line of Fig. 11. We attribute this behavior to the operation of NEMPI in the simulation.
Fig. 11.
Profiles of B_{eq} and B_{z} at z = 0 (i.e., the surface of the turbulent region), for Run B1 with Co = 0.0015, θ = 0 at t/τ_{td} = 1 and along y = 0 (black line). The red solid line gives the Fourierfiltered profiles of B_{z}. The dotted orange line gives the values of B_{eq} through the magnetic structure based on the local turbulent pressure and the solid orange lines denotes its Fourierfiltered value. The blue line represents the value of B_{eq} based on the volumeaveraged velocity, but the local density, which does not show the local suppression of u_{rms}All values have been normalized by the volumeaveraged value of B_{eq}. 
Next, we determined the normalized effective magnetic pressure, , as
where q_{p} is obtained from (Kemel et al. 2013)
Overbars denote here xy averaging and the diagonal components of the total (Reynolds and Maxwell) stress tensor, , are obtained from the DNS as
where the subscript 0 refers to values for B_{0} = 0 and lower case letters denote fluctuations, i.e., and . No summation over the index i is assumed. In Fig. 12 we show the dependence of the minimum value of on the Coriolis number (see Table 2) and for different values of the imposed field (see Table 3). It turns out to be relatively insensitive to the value of Co, but it drops dramatically with increasing strength of the imposed field.
Fig. 12.
Dependence of the minimum value of the normalized effective magnetic pressure, , on the Coriolis number for different values of the imposed field strength, B_{0}/B_{eq0} (the curves have been labeled by an average value), for the runs listed in Tables 2 and 3. Error bars are estimated using the maximum difference of the total mean with the means of each third of the time series. 
Given that the values of shown in Fig. 12 for different values of Co are almost the same, we conclude that those coefficients do not depend on Ω. We can therefore obtain a new set of parameters that is representative of our model with coronal layer and for all values of Ω, which is different from our standard representation without a coronal envelope. This is shown in Fig. 13. The values of the parameters q_{p0} and β_{p} fit very well with the data of Warnecke et al. (2016a), shown as triangles.
Fig. 13.
Dependence of obtained from Fig. 12 by taking an average for different values of Co (filled circles). The blue line is our standard representation for q_{p0} = 32 and β_{p} = 0.058, which corresponds to β_{∗} = 0.33, while the red line is a better fit to the present data giving q_{p0} = 13 and β_{p} = 0.18, which corresponds to β_{∗} = 0.65. We note that the data point at β = 0.85 and is well outside any curve that could fit all the data and has therefore been discarded as an “outlier”. As a comparison, we plot the values of Warnecke et al. (2016a) as triangles. 
The new data can be fitted to Eq. (15) by determining the position and value of the minimum, β_{min} and , respectively. Looking at Fig. 13, we find β_{min} ≈ 0.3 and . We then obtain the fit parameters as (Kemel et al. 2012b)
This results in the following new set of parameters: q_{p0} = 13 and β_{p} = 0.18, which corresponds to β_{∗} = 0.65. In the next section we therefore also present MFS results based on our model with this parameter combination.
Along with q_{p} and , we also determine q_{s} and q_{g} in
resulting in
where g_{i} are components of g which, in our setup, has only a component in the negative z direction. In Fig. 14 we plot the dependencies of q_{p}, q_{s}, and q_{g} on Coriolis number and imposed magnetic field strength. The values for no and weak rotation are consistent with those obtained by Warnecke et al. (2016a) for a nonrotating setup and the same magnetic field strengths. Because of the large errors resulting from the strong variation in space and time, we cannot determine whether the sign is negative or positive. For larger imposed magnetic field strengths, the situation becomes more clear. There, the errors are significantly larger and q_{p}, q_{s}, and q_{g} do not strongly depend on rotation. However, we see a dependence on the imposed magnetic field strength. In particular, q_{p} is positive for B_{0}/B_{eq0} ≥ 0.066 and decreases from around 1.5 to 1.0 for increasing imposed magnetic field strength. Next, q_{s} is also positive for imposed magnetic fields B_{0}/B_{eq0} ≥ 0.066, but has an inconclusive dependence on the imposed magnetic field. Finally, q_{g} is also positive for B_{0}/B_{eq0} ≥ 0.066, it decreases with increasing imposed field until it is nearly zero for B_{0}/B_{eq0} = 0.85. For B_{0}/B_{eq0} = 0.066, q_{g} tends to decrease with increasing angular velocity, but not for the other magnetic field strengths. We should keep in mind that the q_{p} and q_{g} in Eq. (12) are multiplied by the horizontal averaged magnetic field, which will be larger for larger imposed magnetic fields. We also look at the latitudinal dependencies of q_{p}, q_{s}, and q_{g}; see Fig. 15. We find no latitudinal dependence of these coefficients, mostly because the errors are so large. In conclusion, we cannot explain the rotational dependence found in the DNS with just the rotational dependence of q_{p}, q_{s}, and q_{g}.
Fig. 14.
Dependence q_{p}, q_{s}, and q_{g} on the Coriolis number for different values of the imposed field strength, β_{0} = B_{0}/B_{eq0} (the curves have been labeled by an average value), for the runs listed in Tables 2 and 3. Error bars are estimated using the maximum difference of the total mean with the means of each third of the time series. 
Fig. 15.
Dependence of q_{p}, q_{s}, and q_{g} on the colatitude, θ, for Runs C1–C3 and C5, with Co=0.0023 and B_{0}/B_{eq0} = 0.026. Error bars are estimated using the maximum difference of the total mean with the means of each third of the time series. 
4. MFS results
4.1. General aspects
We now present meanfield calculations with a coronal envelope. In all cases, we use ν_{T} = η_{T} = 10^{−3}c_{s}H_{ρ} and u_{rms} = 0.1c_{s}, corresponding to k_{f}H_{ρ} = 33, which is similar to the value of 30 used in the DNS. We first considered the same domain size as in the DNS, i.e., (2π)^{2} × 3π. To alleviate finite domain size effects, we also considered a wider and deeper domain, but with a smaller coronal part that is reduced by half. We thus also considered L_{x} = L_{y} = 6π, L_{z} = 4π, and z_{top} = π.
The MFS lack smallscale turbulent motions, so fewer mesh points can be used. However, to resolve the vertical density contrast of around 12,000, we used 288 mesh points in the z direction. For our domains, are used 192^{2} × 288 meshpoints, but we found no differences in the results when using 96^{2} × 288 meshpoints. For the larger domains, we used 384^{3} meshpoints, but again, with fewer meshpoints the results would have been sufficiently accurate. In all cases, we used B_{0}/B_{eq0} = 0.1, because the DNS (Sect. 3.6) showed that this value maximizes the range of rotation rates where NEMPI is still excited.
Although we find that the parameters q_{p0} ≈ 13 and β_{p} ≈ 0.65 describe the DNS best, there are reasons to consider also other choices. First, there is no good reason why the parameters q_{p0} and β_{p} are so different in different circumstances. One would have expected them to reflect properties of the turbulence which is similar in all the different cases. Second, the response for vertical and horizontal magnetic fields turns out to be different; see Fig. 7 of Losada et al. (2014). There are also other differences between MFS and DNS that we address below. Thus, we cannot expect the two approaches to agree. One objective is therefore to find out just how well the MFS perform relative to the DNS.
4.2. Growth rates
In the DNS, the value of λ/λ_{∗0} is found to drop by about a factor of five as 2Ω/λ_{∗0} increases from 0.2 to 3; see Fig. 16. We compared the MFS for B_{0}/B_{eq0} = 0.1 with the DNS at B_{0}/B_{eq0} = 0.066 and 0.14; see Table 2. We also compared the MFS with the DNS of Losada et al. (2013) without coronal envelope using β_{∗} = 0.33 and, as in Losada et al. (2013), the value β_{∗} = 0.75, which was suitable for one of their sets of MFS. This now explains the ratio of λ/λ_{∗0} that is 2.3 times larger than in their Fig. 2. In addition, however, the growth rates of their DNS were incorrectly scaled by a factor (u_{rms}/c_{s})(k_{f}/k_{1}) ≈ 0.1 × 30 = 3, which is now corrected; see the green line of Fig. 16. Thus, even in the absence of a corona, the growth rates were by a factor of about seven larger in the MFS than in the DNS; see Appendix A with a corrected version of Fig. 2 of Losada et al. (2013).
Fig. 16.
Dependence of λ/λ_{∗0} on 2Ω/λ_{∗0} for B_{0}/B_{eq0} = 0.1 and θ = 0 for the MFS using the large (black line) and small domains (blue solid line), as well as the DNS for neighboring magnetic field strength of B_{0}/B_{eq0} = 0.066 and 0.14 (dashed and solid red lines, respectively). The fat blue line denotes the case of a small domain using a_{g} = −0.4 instead of 0. The green dashed line denotes the DNS of Losada et al. (2013) without coronal layer and B_{0}/B_{eq0} = 0.05. The orange line refers to MFS in a small domain with q_{p0} = 13 and β_{p} = 0.18, so β_{∗} = 0.65. The blue dashed line adjacent to the blue solid line denotes the results obtained when neglecting the derivative term of the profile function Θ_{w}(z). 
Again, the growth rates of the largescale instability in the MFS are significantly larger than those in the DNS. We speculate that this could be caused by a partial cancelation (i.e., a decrease of the effective magnetic pressure gradient) from the q_{g} term in Eq. (12). To illustrate this possibility, we have overplotted in Fig. 16 a case with a_{g} = −0.4; see Eq. (13). Given that in Eq. (12), only the q_{p} term comes with a 1/2 factor, the effective magnetic pressure gradient in the z direction is reduced. Thus, q_{p}/2 is replaced by q_{p}/2 + q_{g} = (1/2 + a_{g})q_{p}, so q_{p} is scaled by a factor 1 + 2a_{g} = 0.2 with respect to the z direction. We see that the growth rate is now suppressed, but also the maximum rotation rate for which NEMPI can operate is reduced. Obviously, a more accurate modeling requires a more detailed knowledge of the actual form of q_{g}, which is likely to be different from that of q_{p}. We also point out that, when calculating ∂_{z}q_{g}(β) in Eq. (12), one of the resulting terms involves the gradient of Θ_{w}(z) in Eq. (14). We have verified that neglecting this term causes only a very minor change in the resulting growth rates (cf. the solid and dashed blue lines in Fig. 16).
We now discuss the nonmonotonic behavior of the growth rate of the magnetic field as a function of the Coriolis number; see Fig. 16. In corresponding DNS without coronal layer (Jabbari et al. 2014), the increase in the growth rate at faster rotation rates (Co > 0.1 or 2Ω/λ_{∗0} > 10) has been explained as a result of largescale dynamo action; see Fig. 8 of Jabbari et al. (2014). In particular, at larger rotation rates, kinetic helicity is produced by a combined effect of uniform rotation and density stratified turbulence. It results in the excitation of an αΩ or α^{2}Ω dynamo instabilities and the generation of a largescale magnetic field. This causes an increase of the growth rate at larger Coriolis numbers, which is also observed in the DNS of Losada et al. (2013) and Jabbari et al. (2014). This implies that two different instabilities are excited in the system, meaning NEMPI at low Coriolis numbers and the meanfield dynamo instability at larger values of the Coriolis numbers. This causes a nonmonotonic behavior of the growth rate of magnetic field as the function of the Coriolis number, observed in Fig. 16.
A nontrivial evolution of the magnetic field in rotating turbulence with a coronal envelope is caused for the following reasons. In an earlier study by Losada et al. (2012) with rotation, no coronal envelope, and an imposed horizontal magnetic field, an analytical expression for the growth rate of NEMPI has been derived in the framework of a meanfield approach. During the magnetic field evolution in the presence of a coronal envelope, as in the simulation of Warnecke et al. (2013a, 2016a) and the present ones, there is a change in the direction of the largescale magnetic field from horizontal at t = 0 to nearly vertical after about one turbulent diffusive time. Therefore, in turbulence with a coronal envelope one can expect to find a mixture of effects caused by the horizontal and vertical magnetic fields.
The growth rates and saturation mechanisms of NEMPI for horizontal and vertical fields are very different (see review by Brandenburg et al. 2016). For horizontal fields, NEMPI saturates rapidly in the nonlinear stage of the magnetic field evolution due to the “potatosack” effect. This means that a local increase of the magnetic field causes a decrease of the negative effective magnetic pressure, which is compensated for by enhanced gas pressure. This leads to an enhanced gas density, so the gas is heavier than its surroundings and sinks. This effect removes horizontal magnetic flux structures from regions in which NEMPI is excited. For a vertical magnetic field, the heavier fluid moves downward along the field without affecting the flux tube, so that NEMPI is not stabilized by the potatosack effect. In this case the operation of NEMPI results in the formation of strong concentrations (see Brandenburg et al. 2013, 2016).
In the nonrotating cases, the growth rates of the largescale magnetic field measured in DNS with coronal envelope and horizontal initial field (Warnecke et al. 2016a), and with vertical magnetic field but without corona (Brandenburg et al. 2013), are the same. This is an indication that the growth rates of the largescale magnetic field in the simulations with coronal envelope are determined by the evolution of the vertical field. On the other hand, the fact that the BRs dissolve after a few turbulent diffusion times is more similar to the behavior of NEMPI with a horizontal imposed magnetic field. Therefore, the evolution of the magnetic structures in the system with coronal envelope is nontrivial and cannot easily be described with our current meanfield models. This also explains the rotational dependency of growth rates in this setup suffers from the mixture of effects.
4.3. Latitudinal dependence
At slower rotation, a decrease in λ/λ_{∗0} can be seen as θ increases from θ = 0 at the poles to θ = 90° at the equator; see Fig. 17. In both plots we also show the growth rates for the larger domain, which are found to be enhanced by a factor of two when 2Ω/λ_{∗0} ≈ 1. These values are about an order of magnitude larger than those for the DNS, but have otherwise a similar functional dependence on Ω and θ. The reason for the difference between DNS and MFS is not entirely clear. It is possible that the meanfield parameter β_{∗} is smaller than what was previously found for simulations with coronal envelope. There is also the possibility that β_{∗} decreases with increasing angular velocity, which is what Rüdiger (private communication) found, although our present simulations presented in Fig. 12 and earlier ones of Losada et al. (2013) did not give such indications.
Fig. 17.
Dependence of λ/λ_{∗0} on θ for B_{0}/B_{eq0} = 0.1 and Co = 0.006. 
4.4. Inclined surface structures
Next, we show crosssections of along the surface z = 0 at a chosen time t_{∗} during the linear growth phase of NEMPI for three values of Co using domain sizes of (2π)^{2} × 3π (Fig. 18) and (6π)^{2}×4π (Fig. 19). In the linear phase, when the magnetic field fluctuations are still growing exponentially in time, only relative values are of physical interest. We therefore present in the following the magnetic field normalized by its maximum value. As in the DNS, the imposed magnetic field points in the y direction. It turns out that rotation not only tends to make the structures inclined relative to the direction of the imposed magnetic field, but it also leads to higher wavenumbers of the structures. Figure 19 shows that the number of nodes in the x direction, which is perpendicular to the magnetic field, remains about constant (k_{x}L_{x}/2π = 4), while that in the y direction along the magnetic field increases from k_{y}L_{y}/2π = 1 (for Co = 0) to 2 (for Co = 0.006) and 4 (for Co = 0.018).
Fig. 18.
Crosssections of through the surface z = 0 at times t_{∗} during the linear growth phase of NEMPI for a MFS with B_{0}/B_{eq0} = 0.1 and three values of Co using the smaller domain size of (2π)^{2} × 3π. In each panel, the magnetic field is scaled to the maximum value. 
In both the DNS and the MFS, the orientation of the inclination is the same and it is opposite to what is seen in Joy’s law. The runs presented here apply only to the poles, but even at lower latitudes (e.g., at 30° latitude, corresponding to θ = 60°) do we find the same antiJoy’s law orientation of the tilt. This is shown in Fig. 20, where we compare two runs with θ = 80° (close to the equator) for a section of the large domain and the full section of the smaller domain, as well as a run with θ = 90° (at the equator). At the equator, the inclination angle with respect to the toroidal direction is 90°, which agrees with what was found in Losada et al. (2012). However, slightly away from the equator, at θ = 80°, the inclination angle is already 45°. This is not an artifact of having chosen a small domain, because even for a three times larger domain, the same inclination angle is found. We thus remain puzzled about this finding, and hope to be able to return to it as soon as further numerical and analytic results can be obtained and assessed.
Fig. 20.
Similar to Fig. 18, but for two runs with θ = 80° for a 1/3 section of the large domain and the full section of the smaller domain, as well as a run with θ = 90°. 
One might speculate that the reason for the difference to Joy’s law has to do with the expansion of rising structures whereas NEMPI structures are caused by contraction, which leads to the opposite tilt. Thinking again of possible applications to the Sun, one may therefore wonder whether the effect of flux concentrations in NEMPI, which must also be responsible for causing the tilt, operate on relatively small scales and might be responsible for causing sunspot rotation (Evershed 1909; Kempf 1910; Pevtsov 2012; Sturrock et al. 2015), for example. In a distributed solar dynamo scenario, the tilt of active regions is primarily caused either by differential rotation (Brandenburg 2005) or simply by the sign of the mean latitudinal field relative to that of the azimuthal field (Jabbari et al. 2015).
4.5. Emergence of solitary structures
The magnetic field patterns in Figs. 18–20 are more or less regular. This is basically because those pictures were taken from the linear phase of the run. In Fig. 21 we show a visualization of along with horizontal flow vectors through the surface for an arbitrarily chosen time t_{∗}/τ_{td} = 40 during the saturated state for the large domain using Co = 0.006. We can now clearly see solitary structures in the form of isolated spots. However, with a few exceptions, most of these structures lack the distinct bipolarity seen in the DNS. We also note that the mean flow is mostly circular around each spot rather than a convergent inflow, as seen in the DNS.
Fig. 21.
Crosssections of (color coded) together with velocity field vectors (white streak lines) through z_{∗} = 0 at t_{∗}/τ_{td} = 40 during the saturated state for a MFS with large domain and Co = 0.006. The white horizontal line through x = 0 marks the position of a BR near y = 0, which is discussed separately. The ellipse marks the position of the BR discussed in the text. 
In Fig. 21 we plot a white horizontal line in the toroidal direction through x = 0 the position of a BR near y = 0. Unlike the DNS, the separation of the two polarities is rather large (about 3H_{ρ}) and there are other spots in almost the same distance. Thus, it is not clear that these two polarities are connected to each other. There is also no clear indication of BRs. To examine this further, we present a side view of this BR in Sect. 4.6 below.
4.6. Side view of BRs
In Fig. 22 we show a longitudinal crosssection of the magnetic field through x_{∗} = 0 at the same time as Fig. 21 during the saturated state. Magnetic flux concentrations are seen to occur at a depth of z/H_{ρ} ≈ −6, which is well below the surface. Near the surface, on the other hand, there is only a relatively small number of vertical magnetic flux structures that seem to close upon themselves over relatively large horizontal distances. We recognize the positive and negative polarities at y/H_{ρ} ≈ ±1.5 and the negative one at y/H_{ρ} ≈ −8 in both Figs. 21 and 22. Conversely, over short distances, bipolar magnetic flux structures separate above the surface, which is consistent with them being the result of a localized subduction of a horizontal flux structure.
Fig. 22.
Crosssections of (color coded) together with magnetic field vectors (white streak lines) through x_{∗} = 0 at an arbitrarily chosen time during the saturated state for a MFS with a large domain. The surface at z = 0 is shown as a white horizontal line. 
To compare with DNS results, we show in Fig. 23 a similar plot of B_{z} together with Fourierfiltered magnetic field vectors in the same plane. A major difference to Fig. 22 is the absence of significant horizontal field in the deeper parts. However, since this horizontal field in the MFS is so deep down (z/H_{ρ} ≲ −8), it is unclear whether it plays any role in explaining the difference in, for example, the growth rates between DNS and MFS seen in Fig. 16. On the other hand, in the deeper parts of the MFS, there are magnetic structures of significant strength, which are not so prominent in the DNS. This is an important difference that would affect global comparisons of, for example, the growth rates of structures shown in Fig. 16.
Fig. 23.
Slice of B_{z}(x_{∗}, y, z, t) (color coded) together Crosssection of B_{z}(x_{∗}, y, z, t) (color coded) together with vectors of Fourierfiltered magnetic field vectors (k < k_{f}/2) superimposed for the DNS Run A4 at the time t/τ_{td} = 0.91 through x_{∗} = 0.5. 
5. Discussion
Our calculations have been performed using idealizing circumstances such as forced turbulence and an isothermal equation of state. This is in many ways different from turbulence in the Sun, which is driven by convection. Nevertheless, some tentative conclusions can be drawn regarding possible applications to sunspot formation. In the surface layers of the Sun, the turnover time τ_{to} = H_{p}/u_{rms} based on the pressure scale height H_{p} is around 5 min. Assuming k_{1}H_{p} = 1, the turbulent diffusive time scale is related to this via
In the simulations, we have k_{f}H_{p} ≈ k_{f}/k_{1} = 30, so τ_{td} would be about 90 times longer than τ_{to}, that is, about eight hours. This would appear suitable in view of applications to the Sun as well, where the formation time of sunspots is of a similar order.
Although the presence of our simplified corona has the effect of allowing BRs to form that reach a significant fraction of the equipartition value with respect to the turbulence, these structures can no longer form when the Coriolis number exceeds a critical value of about 0.02. This value is rather small. However, given that the growth rate of NEMPI does not scale with the inverse turnover time , but with the turbulentdiffusive time , a meaningful measure of the rotation rate in this context could also be the square root of the turbulent Taylor number, , where we took into account that the turbulent viscosity ν_{t} is equal to the turbulent magnetic diffusivity η_{t} (Yousef et al. 2003; Kleeorin & Rogachevskii 1994) and given by η_{t} ≈ u_{rms}/3k_{f} (Sur et al. 2008). The values of Ta^{1/2} are typically on the order of ten when NEMPI begins to be suppressed. Using our estimate of τ_{td} = 8hrs for the solar surface, we find Ta^{1/2} = 0.2, which is well below our critical value of ten.
In the standard mixing length theory of convection, the value of k_{f}H_{p} is estimated to be around six to seven (Kemel et al. 2013), but it is about 30 in the present simulations. Earlier work showed that NEMPI would not work for k_{f}H_{p} much below 15 (Brandenburg et al. 2012). On the other hand, the actual value in the Sun is unclear given that the findings of Hanasoge et al. (2012) did not confirm turbulent velocities in the Sun at the expected levels. A possible resolution to this problem might be the idea that the relevant scales of the energycarrying eddies is much smaller than the inverse pressure scale height (Brandenburg 2016). This would also help making NEMPI more powerful. However, this issue is controversial in view of results by Greer et al. (2015), which showed that the turbulent flows in the Sun might actually be just as large as assumed in standard mixing length theory. Thus, the possibility of NEMPI being responsible for the production of sunspots is being favored particularly in the scenario envisaged by Hanasoge et al. (2012).
An additional complication is that at the solar surface, radiation plays an important role. Meanfield models with radiation transport (Perri & Brandenburg 2018) have shown that the relevant length scale of NEMPI can drop significantly below the value found for isothermal and isentropic stratifications. It is possible, however, that this result is a consequence of not having included the convective flux in such a model. In the Sun, virtually 100% of the energy is transported by convection almost immediately beneath the surface, so radiation should be completely unimportant below the surface. In addition, ionization dynamics can strongly exaggerate the effects of cooling near the surface.
6. Conclusions
Our work has confirmed that NEMPI cannot be exited at Coriolis numbers above a critical value that can be as low as 0.02 or so. The presence of an upper coronal layer was previously found to make the appearance of structures more prominent. However, rotation seems to affect the growth rates more strongly with a coronal envelope than without. In the bulk of the solar convection zone, the Coriolis number is of the order of unity and above, but this is not the case in the surface layers, where the convective time scale is much shorter than the solar rotation period of 25 days. So this may not really be a problem for applications of NEMPI to sunspot formation.
A more severe problem for astrophysical applications of NEMPI are the moderate magnetic field strengths that can presently be achieved with NEMPI. This suggests that some essential physics is still missing. An important ingredient of sunspots physics is convection and its suppression in the presence of magnetic fields. A number of aspects such as radiation and ionization physics, taken in isolation, have as yet not produced more favorable conditions for NEMPI (Bhat & Brandenburg 2016; Perri & Brandenburg 2018).
The difficulty in explaining the spontaneous formation of sunspotlike magnetic flux concentrations is not really alleviated by invoking the rising flux tube scenario. The problem here is that any flux tube rising from some depth to the surface will expand and therefore weaken. This means that some mechanism for reamplification is needed. An important next step might be to invoke a suitable model for convection, or possibly the inclusion of a magnetic suppression of turbulent radiative diffusion as suggested by Kitchatinov & Mazur (2000).
Acknowledgments
J.W. acknowledges funding by the MaxPlanck/Princeton Center for Plasma Physics and funding from the People Program (Marie Curie Actions) of the European Union’s Seventh Framework Programmed (FP7/20072013) under REA grant agreement No. 623609. This work has been supported in part by the NSF Astronomy and Astrophysics Grants Program (grant 1615100), the Research Council of Norway under the FRINATEK (grant 231444), the Swedish Research Council (grant 20125797), and the University of Colorado through its support of the George Ellery Hale visiting faculty appointment. I.R. acknowledges the hospitality of NORDITA and Max Planck Institute for Solar System Research in Göttingen. 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 CNS0821794), 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. Additional simulations have been carried out on supercomputers at GWDG, on the Max Planck supercomputer at RZG in Garching, in the facilities hosted by the CSC – IT Center for Science in Espoo, Finland, which are financed by the Finnish ministry of education.
References
 Archontis, V., MorenoInsertis, F., Galsgaard, K., Hood, A., & O’Shea, E. 2004, A&A, 426, 1047 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Archontis, V., MorenoInsertis, F., Galsgaard, K., & Hood, A. W. 2005, ApJ, 635, 1299 [Google Scholar]
 Bhat, P., & Brandenburg, A. 2016, A&A, 587, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Birch, A. C., Schunker, H., Braun, D. C., et al. 2016, Sci. Adv., 2, e1600557 [NASA ADS] [CrossRef] [Google Scholar]
 Bourdin, P. A., & Brandenburg, A. 2018, ApJ, 869, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A. 2005, ApJ, 625, 539 [Google Scholar]
 Brandenburg, A. 2016, ApJ, 832, 6 [Google Scholar]
 Brandenburg, A., & Campbell, C. 1997, in Accretion Disks  New Aspects, eds. E. MeyerHofmeister, & H. Spruit (Berlin: Springer Verlag), Lecture Notes in Physics, 487, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A., Nordlund, A., Stein, R. F., & Torkelsson, U. 1995, ApJ, 446, 741 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A., Kemel, K., Kleeorin, N., Mitra, D., & Rogachevskii, I. 2011, ApJ, 740, L50 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A., Kemel, K., Kleeorin, N., & Rogachevskii, I. 2012, ApJ, 749, 179 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A., Kleeorin, N., & Rogachevskii, I. 2013, ApJ, 776, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A., Gressel, O., Jabbari, S., Kleeorin, N., & Rogachevskii, I. 2014, A&A, 562, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brandenburg, A., Rogachevskii, I., & Kleeorin, N. 2016, New J. Phys., 18, 125011 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, F., Rempel, M., & Fan, Y. 2017, ApJ, 846, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Cheung, M. C. M., & Isobe, H. 2014, Liv. Rev. Sol. Phys., 11, 3 [Google Scholar]
 Cheung, M. C. M., Rempel, M., Title, A. M., & Schüssler, M. 2010, ApJ, 720, 233 [NASA ADS] [CrossRef] [Google Scholar]
 Choudhuri, A. R., & Gilman, P. A. 1987, ApJ, 316, 788 [NASA ADS] [CrossRef] [Google Scholar]
 Evershed, J. 1909, MNRAS, 69, 454 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, Y. 2001, ApJ, 554, L111 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, Y. 2009, Liv. Rev. Sol. Phys., 6, 4 [Google Scholar]
 Fan, Y., & Fang, F. 2014, ApJ, 789, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, Y., Fisher, G. H., & Deluca, E. E. 1993, ApJ, 405, 390 [Google Scholar]
 Fournier, Y., Arlt, R., Ziegler, U., & Strassmeier, K. G. 2017, A&A, 607, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Galloway, D. J., & Weiss, N. O. 1981, ApJ, 243, 945 [NASA ADS] [CrossRef] [Google Scholar]
 Greer, B. J., Hindman, B. W., Featherstone, N. A., & Toomre, J. 2015, ApJ, 803, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasoge, S. M., Duvall, T. L., & Sreenivasan, K. R. 2012, Proc. Natl. Acad. Sci., 109, 11928 [Google Scholar]
 Jabbari, S., Brandenburg, A., Losada, I. R., Kleeorin, N., & Rogachevskii, I. 2014, A&A, 568, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jabbari, S., Brandenburg, A., Kleeorin, N., Mitra, D., & Rogachevskii, I. 2015, ApJ, 805, 166 [NASA ADS] [CrossRef] [Google Scholar]
 Jabbari, S., Brandenburg, A., Mitra, D., Kleeorin, N., & Rogachevskii, I. 2016, MNRAS, 459, 4046 [NASA ADS] [CrossRef] [Google Scholar]
 Jabbari, S., Brandenburg, A., Kleeorin, N., & Rogachevskii, I. 2017, MNRAS, 467, 2753 [NASA ADS] [Google Scholar]
 Käpylä, P. J., Brandenburg, A., Kleeorin, N., Mantere, M. J., & Rogachevskii, I. 2012, MNRAS, 422, 2465 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J., Brandenburg, A., Kleeorin, N., Mantere, M. J., & Rogachevskii, I. 2013, IAU Symp., 294, 283 [NASA ADS] [Google Scholar]
 Käpylä, P. J., Brandenburg, A., Kleeorin, N., Käpylä, M. J., & Rogachevskii, I. 2016, A&A, 588, A150 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kemel, K., Brandenburg, A., Kleeorin, N., Mitra, D., & Rogachevskii, I. 2012a, Sol. Phys., 280, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Kemel, K., Brandenburg, A., Kleeorin, N., & Rogachevskii, I. 2012b, Astron. Nachr., 333, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Kemel, K., Brandenburg, A., Kleeorin, N., Mitra, D., & Rogachevskii, I. 2013, Sol. Phys., 287, 293 [NASA ADS] [CrossRef] [Google Scholar]
 Kempf, P. 1910, Astron. Nachr., 185, 197 [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov, L. L., & Mazur, M. V. 2000, Sol. Phys., 191, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Kitiashvili, I. N., Kosovichev, A. G., Wray, A. A., & Mansour, N. N. 2010, ApJ, 719, 307 [Google Scholar]
 Kleeorin, N., & Rogachevskii, I. 1994, 50, 2716 [Google Scholar]
 Kleeorin, N. I., Rogachevskii, I. V., & Ruzmaikin, A. A. 1989, Sov. Astron. Lett., 15, 274 [NASA ADS] [Google Scholar]
 Kleeorin, N., Rogachevskii, I., & Ruzmaikin, A. 1990, JETP, 70, 878 [Google Scholar]
 Kleeorin, N., Mond, M., & Rogachevskii, I. 1993, Phys. Fluids B, 5, 4128 [NASA ADS] [CrossRef] [Google Scholar]
 Kleeorin, N., Mond, M., & Rogachevskii, I. 1996, A&A, 307, 293 [NASA ADS] [Google Scholar]
 Kleeorin, N., Kuzanyan, K., Moss, D., et al. 2003, A&A, 409, 1097 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Losada, I. R., Brandenburg, A., Kleeorin, N., Mitra, D., & Rogachevskii, I. 2012, A&A, 548, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Losada, I. R., Brandenburg, A., Kleeorin, N., & Rogachevskii, I. 2013, A&A, 556, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Losada, I. R., Brandenburg, A., Kleeorin, N., & Rogachevskii, I. 2014, A&A, 564, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Losada, I. R., Warnecke, J., Glogowski, K., et al. 2017, IAU Symp., 327, 46 [NASA ADS] [Google Scholar]
 Masada, Y., & Sano, T. 2016, ApJ, 822, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Mitra, D., Brandenburg, A., Kleeorin, N., & Rogachevskii, I. 2014, MNRAS, 445, 761 [NASA ADS] [CrossRef] [Google Scholar]
 Parker, E. N. 1978, ApJ, 221, 368 [NASA ADS] [CrossRef] [Google Scholar]
 Perri, B., & Brandenburg, A. 2018, A&A, 609, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Peter, H., Warnecke, J., Chitta, L. P., & Cameron, R. H. 2015, A&A, 584, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pevtsov, A. A. 2012, Astrophys. Space Sci. Proc., 30, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Rempel, M., & Cheung, M. C. M. 2014, ApJ, 785, 90 [Google Scholar]
 Rogachevskii, I., & Kleeorin, N. 2007, Phys. Rev. E, 76, 056307 [NASA ADS] [CrossRef] [Google Scholar]
 Rüdiger, G., Kitchatinov, L. L., & Brandenburg, A. 2011, Sol. Phys., 269, 3 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Rüdiger, G., Kitchatinov, L. L., & Schultz, M. 2012, Astron. Nachr., 333, 84 [NASA ADS] [CrossRef] [Google Scholar]
 Schmieder, B., Archontis, V., & Pariat, E. 2014, Space Sci. Rev., 186, 227 [NASA ADS] [CrossRef] [Google Scholar]
 Singh, N. K., Raichur, H., & Brandenburg, A. 2016, ApJ, 832, 120 [NASA ADS] [CrossRef] [Google Scholar]
 Spruit, H. C. 1979, Sol. Phys., 61, 363 [NASA ADS] [CrossRef] [Google Scholar]
 Stein, R. F., & Nordlund, Å. 2012, ApJ, 753, L13 [Google Scholar]
 Sturrock, Z., Hood, A. W., Archontis, V., & McNeill, C. M. 2015, A&A, 582, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sur, S., Brandenburg, A., & Subramanian, K. 2008, MNRAS, 385, L15 [NASA ADS] [CrossRef] [Google Scholar]
 Tao, L., Weiss, N. O., Brownjohn, D. P., & Proctor, M. R. E. 1998, ApJ, 496, L39 [NASA ADS] [CrossRef] [Google Scholar]
 Warnecke, J., & Brandenburg, A. 2010, A&A, 523, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Warnecke, J., Brandenburg, A., & Mitra, D. 2011, A&A, 534, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Warnecke, J., Brandenburg, A., & Mitra, D. 2012a, JSWSC, 2, A11 [Google Scholar]
 Warnecke, J., Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2012b, Sol. Phys., 280, 299 [NASA ADS] [CrossRef] [Google Scholar]
 Warnecke, J., Losada, I. R., Brandenburg, A., Kleeorin, N., & Rogachevskii, I. 2013a, ApJ, 777, L37 [NASA ADS] [CrossRef] [Google Scholar]
 Warnecke, J., Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2013b, ApJ, 778, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Warnecke, J., Losada, I. R., Brandenburg, A., Kleeorin, N., & Rogachevskii, I. 2016a, A&A, 589, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Warnecke, J., Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2016b, A&A, 596, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Warnecke, J., Chen, F., Bingert, S., & Peter, H. 2017, A&A, 607, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yousef, T. A., Brandenburg, A., & Rüdiger, G. 2003, A&A, 411, 321 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zhang, H., & Brandenburg, A. 2018, ApJ, 862, L17 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Comparison of growth rates in DNS and MFS of Losada et al. (2013)
In Losada et al. (2013), the growth rates are accidentally scaled by a factor (u_{rms}/c_{s})(k_{f}/k_{1}) ≈ 0.1 × 30 = 3. In addition, they use β_{∗} = 0.75, which was suitable for one of their sets of MFS, but not for the other. Therefore, the growth rates of their MFS exceeded those of their DNS by a factor of about seven. The corrected version of their Fig. 2 is shown in Fig. A.1.
Fig. A.1.
Corrected version of Fig. 2 of Losada et al. (2013) showing the dependence of λ/λ_{∗0} on 2Ω/λ_{∗0} for DNS (red dashed line), compared with MFS (i) where q_{p0} = 20 and β_{p} = 0.167 (black solid line), and MFS (ii) where q_{p0} = 32 and β_{p} = 0.058 (blue dashdotted line). 
All Tables
Comparison of BR formation as a function of rotation in terms of Co and imposed magnetic field B_{0}/B_{eq0}.
All Figures
Fig. 1.
Growth rate λ (black line), the peak spectral magnetic field (red), and the maximum of the Fourierfiltered vertical magnetic field (blue) for Runs A1–A5 at different resolutions. 

In the text 
Fig. 2.
Vertical magnetic field B_{z} at the surface (z = 0) for runs with different resolution (Runs A1–A3 and A5) at the time when the structures are strongest. The resolution is indicated by the number of grid points at the top of each panel. 

In the text 
Fig. 3.
Evolution of vs. t/τ_{td} for Run A2 showing exponential growth with growth rate (black solid line), compared with that of (red dashed line) and (blue solid line) at z = 0. We also show the unity line corresponding to B_{eq} (orange dotted line). 

In the text 
Fig. 4.
Vectors of the largescale (Fourierfiltered) flow superimposed on a color scale representation of B_{z} at the surface at t/τ_{td} = 0.91 for Run A4 (upper panel) and t/τ_{td} = 1.12 for Run A5 (lower panel). Note that the coordinate system has been rotated by 90°, so y points to the right and x points downward. 

In the text 
Fig. 5.
Power spectra of B_{z} at the surface (z = 0) for Run A2 at several times around the maximum growth for zero rotation. The inset shows the time dependence (time in units of τ_{td}) for k = k_{∗} = 2 k_{1} and for the values averaged between k/k_{1} = 1 and 4. 

In the text 
Fig. 6.
Rotational dependency on the active region formation. We show the vertical magnetic field B_{z}/B_{eq} for Runs A2, B2 and F2 corresponding to Co = 0, 0.0012, 0.0076 at the time of the strongest BR. In the first panel for Co = 0, the value of θ is insignificant. 

In the text 
Fig. 7.
Rotational dependency of the growth rate λ, the peak spectral magnetic field , and the maximum of the Fourierfiltered vertical magnetic field for colatitudes 0° (black line), 30° (blue line), 60° (purple line), and 90° (red line). The error bars are either the errors of the exponential fit (λ; see Table 2), or estimated as 10% of the actual value. 

In the text 
Fig. 8.
Summary of the results for the different colatitudes showing the normalized growth rate of the Fourierfiltered zcomponent of the magnetic field (top panel) for B_{0}/B_{eq0} = 0.02, the averaged magnetic spectrum over k ≤ 4 (middle panel) and the maximum of the Fourierfiltered vertical magnetic field (bottom panel) for different values of the Coriolis number and colatitude. As elsewhere, error bars are either the errors of the exponential fit (λ; see Table 3), or estimated as 10% of the actual value. 

In the text 
Fig. 9.
Vertical surface magnetic field B_{z}/B_{eq} for Runs C1–C3 and C5 with Co = 0.0012. 

In the text 
Fig. 10.
Rotational dependency of the growth rate λ, the peak spectral magnetic field , and the maximum of the Fourierfiltered vertical magnetic field for various averaged imposed magnetic field strengths B_{0}/B_{eq0} = 0.026 (blue line), 0.066 (red), 0.14 (green), 0.28 (cyan), and 0.85 (black). Error bars denote either the errors of the exponential fit (λ; see Table 2), or are estimated as 10% of the actual value. 

In the text 
Fig. 11.
Profiles of B_{eq} and B_{z} at z = 0 (i.e., the surface of the turbulent region), for Run B1 with Co = 0.0015, θ = 0 at t/τ_{td} = 1 and along y = 0 (black line). The red solid line gives the Fourierfiltered profiles of B_{z}. The dotted orange line gives the values of B_{eq} through the magnetic structure based on the local turbulent pressure and the solid orange lines denotes its Fourierfiltered value. The blue line represents the value of B_{eq} based on the volumeaveraged velocity, but the local density, which does not show the local suppression of u_{rms}All values have been normalized by the volumeaveraged value of B_{eq}. 

In the text 
Fig. 12.
Dependence of the minimum value of the normalized effective magnetic pressure, , on the Coriolis number for different values of the imposed field strength, B_{0}/B_{eq0} (the curves have been labeled by an average value), for the runs listed in Tables 2 and 3. Error bars are estimated using the maximum difference of the total mean with the means of each third of the time series. 

In the text 
Fig. 13.
Dependence of obtained from Fig. 12 by taking an average for different values of Co (filled circles). The blue line is our standard representation for q_{p0} = 32 and β_{p} = 0.058, which corresponds to β_{∗} = 0.33, while the red line is a better fit to the present data giving q_{p0} = 13 and β_{p} = 0.18, which corresponds to β_{∗} = 0.65. We note that the data point at β = 0.85 and is well outside any curve that could fit all the data and has therefore been discarded as an “outlier”. As a comparison, we plot the values of Warnecke et al. (2016a) as triangles. 

In the text 
Fig. 14.
Dependence q_{p}, q_{s}, and q_{g} on the Coriolis number for different values of the imposed field strength, β_{0} = B_{0}/B_{eq0} (the curves have been labeled by an average value), for the runs listed in Tables 2 and 3. Error bars are estimated using the maximum difference of the total mean with the means of each third of the time series. 

In the text 
Fig. 15.
Dependence of q_{p}, q_{s}, and q_{g} on the colatitude, θ, for Runs C1–C3 and C5, with Co=0.0023 and B_{0}/B_{eq0} = 0.026. Error bars are estimated using the maximum difference of the total mean with the means of each third of the time series. 

In the text 
Fig. 16.
Dependence of λ/λ_{∗0} on 2Ω/λ_{∗0} for B_{0}/B_{eq0} = 0.1 and θ = 0 for the MFS using the large (black line) and small domains (blue solid line), as well as the DNS for neighboring magnetic field strength of B_{0}/B_{eq0} = 0.066 and 0.14 (dashed and solid red lines, respectively). The fat blue line denotes the case of a small domain using a_{g} = −0.4 instead of 0. The green dashed line denotes the DNS of Losada et al. (2013) without coronal layer and B_{0}/B_{eq0} = 0.05. The orange line refers to MFS in a small domain with q_{p0} = 13 and β_{p} = 0.18, so β_{∗} = 0.65. The blue dashed line adjacent to the blue solid line denotes the results obtained when neglecting the derivative term of the profile function Θ_{w}(z). 

In the text 
Fig. 17.
Dependence of λ/λ_{∗0} on θ for B_{0}/B_{eq0} = 0.1 and Co = 0.006. 

In the text 
Fig. 18.
Crosssections of through the surface z = 0 at times t_{∗} during the linear growth phase of NEMPI for a MFS with B_{0}/B_{eq0} = 0.1 and three values of Co using the smaller domain size of (2π)^{2} × 3π. In each panel, the magnetic field is scaled to the maximum value. 

In the text 
Fig. 19.
Similar to Fig. 18, but for the large domain, (6π)^{2} × 4π. 

In the text 
Fig. 20.
Similar to Fig. 18, but for two runs with θ = 80° for a 1/3 section of the large domain and the full section of the smaller domain, as well as a run with θ = 90°. 

In the text 
Fig. 21.
Crosssections of (color coded) together with velocity field vectors (white streak lines) through z_{∗} = 0 at t_{∗}/τ_{td} = 40 during the saturated state for a MFS with large domain and Co = 0.006. The white horizontal line through x = 0 marks the position of a BR near y = 0, which is discussed separately. The ellipse marks the position of the BR discussed in the text. 

In the text 
Fig. 22.
Crosssections of (color coded) together with magnetic field vectors (white streak lines) through x_{∗} = 0 at an arbitrarily chosen time during the saturated state for a MFS with a large domain. The surface at z = 0 is shown as a white horizontal line. 

In the text 
Fig. 23.
Slice of B_{z}(x_{∗}, y, z, t) (color coded) together Crosssection of B_{z}(x_{∗}, y, z, t) (color coded) together with vectors of Fourierfiltered magnetic field vectors (k < k_{f}/2) superimposed for the DNS Run A4 at the time t/τ_{td} = 0.91 through x_{∗} = 0.5. 

In the text 
Fig. A.1.
Corrected version of Fig. 2 of Losada et al. (2013) showing the dependence of λ/λ_{∗0} on 2Ω/λ_{∗0} for DNS (red dashed line), compared with MFS (i) where q_{p0} = 20 and β_{p} = 0.167 (black solid line), and MFS (ii) where q_{p0} = 32 and β_{p} = 0.058 (blue dashdotted line). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.