###### Abstract

The negative effective magnetic pressure instability discovered recently in direct numerical simulations (DNS) may play a crucial role in the formation of sunspots and active regions in the Sun and stars. This instability is caused by a negative contribution of turbulence to the effective mean Lorentz force (the sum of turbulent and non-turbulent contributions) and results in formation of large-scale inhomogeneous magnetic structures from initial uniform magnetic field. Earlier investigations of this instability in DNS of stably stratified, externally forced, isothermal hydromagnetic turbulence in the regime of large plasma beta are now extended into the regime of larger scale separation ratios where the number of turbulent eddies in the computational domain is about 30. Strong spontaneous formation of large-scale magnetic structures is seen even without performing any spatial averaging. These structures encompass many turbulent eddies. The characteristic time of the instability is comparable to the turbulent diffusion time, , where is the turbulent diffusivity and is the scale of the domain. DNS are used to confirm that the effective magnetic pressure does indeed become negative for magnetic field strengths below the equipartition field. The dependence of the effective magnetic pressure on the field strength is characterized by fit parameters that seem to show convergence for larger values of the magnetic Reynolds number.

## Section 1 Introduction

The 11-year solar activity cycle manifests itself through the periodic variation of the sunspot number. Sunspots consist of vertical magnetic fields with a strength of up to (see, e.g. Parker, 1979; Ossendrijver, 2003). It is generally believed that these fields continue in a similarly concentrated fashion also beneath the surface in the form of magnetic flux tubes or fibers (Priest, 1982). It was thought that fibral magnetic fields constitute a lower energy state and might therefore be preferred. Tube-like magnetic fields were frequently seen in hydromagnetic turbulence simulations (Nordlund et al., 1992; Brandenburg et al., 1995, 1996), but those tubes were similar to the vortex tubes in hydrodynamic turbulence. Those vortex tubes have typical diameters comparable to the viscous length. Similarly, the aforementioned magnetic tubes have a thickness comparable to the resistive length. For the Sun, however, the resulting tube thickness would be too small to be relevant for sunspots.

At the surface, strong magnetic flux concentrations also form in regions of strong flow convergence, but the size of these regions is too small for sunspots, because sunspots are usually much bigger than a single granular convection cell. In fact, a typical sunspot can have a diameter of some 30 pressure scale heights. The tremendous size of sunspots has therefore been used as an argument that they are not made near the surface, but at much greater depth far away from the surface. At the bottom of the convection zone the convection cells are big enough and could in principle be responsible for producing much bigger flux concentrations. Magnetic flux tubes can also form through the action of shear as shown by various simulations (Cline et al., 2003; Guerrero & Käpylä, 2011). While shear is likely an important ingredient of the solar dynamo, it remains unclear whether the resulting magnetic tubes are really able to produce sunspots as a result of them piercing through the solar surface and, more importantly, whether one should consider them being tied to the deep shear layers of the Sun.

In this paper we discuss an alternative scenario in which sunspots are shallow phenomena that are not anchored at the bottom of the convection zone. Various mechanisms have been discussed. Of particular interest here are mechanisms that are based on the suppression of turbulence by magnetic fields. In the mechanism of Kitchatinov & Mazur (2000) it is the suppression of the turbulent heat flux while in the mechanism of Rogachevskii & Kleeorin (2007) it is the suppression of the turbulent hydromagnetic pressure. Both mechanisms lead to a linear large-scale instability in a stratified medium. However, these mechanisms may be of different importance in different layers.

In this study we focus mainly on the second mechanism, which has recently been studied in direct numerical simulations (DNS) as well as in mean-field calculations. This mechanism is called the negative effective magnetic pressure instability (or NEMPI for short). NEMPI is a convective type instability that is similar to the interchange instability in plasmas (Tserkovnikov, 1960; Newcomb, 1961; Priest, 1982) and the magnetic buoyancy instability (Parker, 1966). However, the free energy in interchange and magnetic buoyancy instabilities is due to the gravitational field, while in NEMPI it is provided by the small-scale turbulence. NEMPI is caused by the suppression of turbulent hydromagnetic pressure by the mean magnetic field. When fluid Reynolds numbers is larger than 1 and the mean magnetic field is less than the equipartition magnetic fields, the negative turbulent contribution to the mean Lorentz force is enough large so that the effective mean magnetic pressure (the sum of turbulent and non-turbulent contributions) become negative (Kleeorin et al., 1990; Rogachevskii & Kleeorin, 2007). This is the main reason for the excitation of the large-scale instability that results in formation of large-scale inhomogeneous magnetic structures.

The effect of the suppression of the turbulent heat flux has not yet been studied as extensively as NEMPI since the original paper of Kitchatinov & Mazur (2000), although Kitchatinov & Olemskoy (2006) have used that model to study the decay of sunspots. Moreover, there is now some evidence that the effects anticipated by Kitchatinov & Mazur (2000) may have already been in operation in various simulations of solar convection where the spontaneous formation of pores has been seen (Stein et al., 2011). Another example is the large-eddy simulation of Kitiashvili et al. (2010) where an imposed vertical large-scale magnetic field gets concentrated into giant vortices. This result is reminiscent of that of Tao et al. (1998), who found a segregation into magnetized and nearly unmagnetized regions in stratified convection simulations. However, it has not yet been possible to obtain the large-scale magnetic structures resembling the sunspots, in models with a strong imposed magnetic flux tube structure at the bottom of the domain (see, e.g. Rempel et al., 2009a, b; Cheung et al., 2010; Rempel, 2011a, b).

In the rest of this paper, we focus on NEMPI, which was first found in mean-field calculations of a stratified layer (Kleeorin et al., 1996; Rogachevskii & Kleeorin, 2007; Brandenburg et al., 2010, 2011b). However, those results remained unconvincing until NEMPI was also discovered in DNS (Brandenburg et al., 2011a, hereafter BKKMR). It it therefore most appropriate to begin our discussion with the latter.

## Section 2 The model

Following the earlier work of BKKMR, we solve the equations for the velocity , the magnetic vector potential , and the density ,

(\theequation) |

(\theequation) |

(\theequation) |

where is the kinematic viscosity, is the magnetic diffusivity due to Spitzer conductivity of the plasma, is the magnetic field, is the imposed uniform field, is the current density, is the vacuum permeability, is the traceless rate of strain tensor, and commas denote partial differentiation. The forcing function consists of random, white-in-time, plane non-polarized waves with a certain average wavenumber. The forcing strength is such that the turbulent rms velocity is approximately independent of with . The gravitational acceleration is chosen such that , which leads to a density contrast between bottom and top of . Here, is the density scale height. We consider a domain of size in Cartesian coordinates , with periodic boundary conditions in the and directions and stress-free perfectly conducting boundaries at top and bottom (). The volume-averaged density is therefore constant in time and equal to its initial value, .

Our simulations are characterized by the scale separation ratio, , the fluid Reynolds number , the magnetic Prandtl number and the magnetic Reynolds number . Following earlier work (Brandenburg et al., 2011b), it is clear that NEMPI is more effective for small values of , so here we choose and in the range –. The magnetic field is expressed in units of the local equipartition field strength near the top, , while is specified in units of the volume averaged value, . In addition to visualizations of the actual magnetic field, we also monitor , where is an average over and a certain time interval . Time is expressed in eddy turnover times, . For comparison, we also consider the turbulent-diffusive timescale, , where is the estimated turbulent magnetic diffusivity. Another diagnostic quantity is the rms magnetic field in the Fourier mode of , referred to as , which is here taken as an average over , and is close to the top at . (Note that does not include the imposed field at .)

The simulations are performed with the Pencil Code,^{1}^{1}1http://pencil-code.googlecode.com
which uses sixth-order explicit finite differences in space and a
third-order accurate time stepping method.
We use numerical resolutions of and mesh points
when , and when .
To capture mean-field effects on the slower turbulent-diffusive
timescale, which is times smaller than
the dynamical timescale, we perform simulations for several thousand
turnover times.

## Section 3 Results

In simulations, the clearest indication of a spontaneous flux concentration is seen when the scale separation ratio is large. In BKKMR, we only used . Here we consider calculations where this ratio is twice as large. A useful diagnostic is the magnetic field averaged along the direction of the imposed field, i.e., along the direction. In particular, we shall be looking at the component of the field, i.e., we look at . To see the effect even more clearly, we perform an additional time average over about one hundred turnover times. This average is then referred to as . In Figure 1 we show for . The other parameters are and . An inhomogeneous magnetic structure forms first near the surface (at ), but then the structure propagates downward. This is consistent with our interpretation that this is caused by negative effective magnetic pressure operating on the scale of many turbulent eddies. Indeed, a local decrease of the effective magnetic pressure must be compensated by an increase in gas pressure, which implies higher density, so the structure becomes heavier and sinks in the nonlinear stage of NEMPI. This is also seen in three-dimensional visualizations without averaging; see Figure 2 with the same parameters as in Figure 1.

To confirm that NEMPI really is a large-scale instability, we would expect to see an exponential growth phase. This is shown in the right-hand panel of Figure 3, where we show the growth of versus time. We give time both in turbulent-diffusive times (lower abscissa) as well as in eddy turnover times (upper abscissa). We do indeed see that there is an exponential growth phase which lasts for about one turbulent-diffusive time, i.e., the growth rate is comparable to , where is the expected turbulent magnetic diffusivity. However, compared with the eddy turnover time, , the turbulent-diffusive time scale is times slower. This illustrates that NEMPI is indeed a very slow process compared with, for example, the saturation of the overall rms magnetic field (left-hand panel of Figure 3).

While the rate of NEMPI may be too slow to explain the relatively rapid appearance of sunspots on a time scale of days, it would explain how one can produce magnetic structures on length and time scales much bigger than the naturally occurring scales in the upper layers of the Sun. This discrepancy is exactly one of the reasons why one normally places the formation of active regions at the bottom of the convection zone (Golub et al., 1981). Here we see a clear example that this conclusion may be not correct.

Next, we turn to a simulation where the scale separation ratio, , is only half as big, i.e., . In that case we also see an exponential growth phase, but the growth is slower (even in terms of turbulent-diffusive times), lasts longer, and amplifies only by about one order of magnitude; see Figure 4.

Again, in agreement with related work involving the suppression of turbulent heat flux, the most unstable mode has a horizontal wavelength comparable to the vertical scale height of the layer; see Fig. 2 of Kitchatinov & Mazur (2000); hereafter referred to as KM. This is also the case for NEMPI, as is shown in Figure 5, where we compare an instantaneous plot of with a time averaged one, , making the appearance of large-scale structures even more pronounced. Next, in Figure 6 we compare cross-sections of (for fixed values of , , and ; top panel) with corresponding averages (middle panel) and averages (bottom panel). Here, , , and . Without averaging, no clear magnetic structure is seen yet, but the structures become clearly more pronounced with and averaging. Runs with similar parameters have been shown in BKKMR, for a computational domain whose extent is instead of .

The large-scale flux concentrations have an amplitude of only and are therefore not easily seen in single snapshots, where the field reaches peak strengths comparable to . Furthermore, as for any linear instability, the flux concentrations form a repetitive pattern, and are in that sense similar to flux concentrations seen in the calculations of KM that were based on the magnetic suppression of the turbulent heat flux. However, there are indications that at larger values of , flux concentrations occur more rarely, which might be more realistic in view of astrophysical applications.

## Section 4 Quantifying the negative effective magnetic pressure effect

An important condition for the formation of structures by the mechanisms of KM and Rogachevskii & Kleeorin (2007) is sufficient scale scale separation. In other words, both the suppression of convective heat flux and the suppression of total effective pressure have in common that they only work if a substantial number of eddies is involved in a turbulent structure under consideration. If this is not the case, turbulent transport coefficients become increasingly inefficient. This is a natural feature of large-scale effects such as this one. Here, we mean by turbulent transport coefficients any of the mean-field coefficients that relate correlation functions of small-scale quantities in terms of large-scale quantities.

Small- and large-scale quantities refer here simply to a suitably defined averaging procedure so that the velocity , for example, can be split into a mean (or large-scale) quantity and a fluctuation . An important example of a turbulent transport coefficient is the turbulent viscosity that emerges when relating the Reynolds stress to the spatial derivatives of the mean flow in the averaged momentum equation,

(\theequation) |

Here, is the average density, and correlations with density fluctuations are neglected. The simplest parameterization for is

(\theequation) |

where is the turbulent viscosity and is the turbulent bulk viscosity. This relation is also known as the Boussinesq ansatz, especially when contrasted with representations where the effect is included, which is responsible for producing differential rotation in the Sun (Rüdiger, 1989). Here, however, we shall focus on magnetic effects.

In general, when there are magnetic fields, the right-hand side of equation of motion (4) has to be replaced by the sum of Reynolds and Maxwell stresses

(\theequation) |

where the superscript f indicates contributions from the fluctuating field. Furthermore, in the presence of a mean magnetic field, symmetry arguments allow one to write down additional components, in particular those proportional to and .

Expressing in terms of the mean field, the leading terms are

(\theequation) |

where the dots indicate the presence of additional terms such as the turbulent viscosity term mentioned earlier and terms that enter when the stratification affect the anisotropy of the turbulence further. Note in particular the definition of the signs of the terms involving the functions and . This becomes clear when writing down the mean Maxwell stress resulting from both mean and fluctuating fields, i.e.,

(\theequation) |

Thus, the signs are defined such that for positive and the effects of magnetic stress and magnetic pressure are reduced and the signs of the net effects may even change. Equations (4) and (4) have been derived using the spectral relaxation approach (Kleeorin et al., 1990, 1996; Rogachevskii & Kleeorin, 2007) and the renormalization procedure (Kleeorin & Rogachevskii, 1994).

A broad range of different DNS in stratified turbulence (Brandenburg et al., 2011a, b) or turbulent convection (Käpylä et al., 2011) have now confirmed that is positive for , but is small and negative. A positive value of (but with large error bars) was originally reported for unstratified turbulence (Brandenburg et al., 2010). Later, stratified simulations with isothermal stable stratification (Brandenburg et al., 2011b) and convectively unstable stratification (Käpylä et al., 2011) show that it is small and negative. Nevertheless, is consistently larger than unity provided and above unity while is below a certain critical value that is around 0.5. This is shown in Figure 7, where we plot the effective magnetic pressure,

(\theequation) |

for different values of using and . Note that the minimum of is deeper for the case with and becomes then shallower.

It should be noted that the coefficients for the fit formulae are chosen such that the plot is most accurate for small values of , because this is what is most important for the onset of NEMPI. Especially for larger values of the fit formula does no longer reproduce correctly the critical value for which changes sign.

Note that is well below unity. This implies that it is probably not possible to produce flux concentrations stronger than half the equipartition field strength. So, making sunspots with this mechanism alone is maybe unlikely, and other effects such as that of KM may be needed. It is possible that this mechanisms works preferentially in the uppermost layers, provided enough flux has already been accumulated. This may then be achieved with NEMPI which also works in somewhat lower layers.

To compare the resulting functions in a systematic fashion for different parameters, we use the fit formula (Kemel et al., 2012)

(\theequation) |

The resulting dependencies , , and are shown in Figure 8. We see that varies relatively little between 0.1 and 0.2 and is typically around 0.15. For small values of , drops from 1 to 0.1 and stays then approximately constant, while rises proportional to for and then levels off at a value around 40.

The significance of a positive value comes from mean-field simulations with indicating the formation of three-dimensional (non-axisymmetric) flux concentrations (Brandenburg et al., 2010). This result was later identified to be a direct consequence of having (Kemel et al., 2012). Before making any further conclusions, it is important to assess the effect of other terms that have been neglected. Two of them are related to the vertical stratification, i.e. additional terms in Eq. (4) that are proportional to and with being gravity. The coefficient of the former term seems to be small (Käpylä et al., 2011) and the second only has an effect when there is a vertical or inclined imposed magnetic field. However, there could be other terms such as as well as and when the scale separation is not enough large.

## Section 5 Conclusions

In the present paper we have performed detailed investigations of NEMPI detected recently by BKKMR. Most notably, we have extended the values of the scale separation ratio, , from 15 to 30. In this case, the spontaneous formation of magnetic structures becomes particularly evident and can clearly be noticed even without any averaging. Whether or not the particular structures seen in DNS really have a correspondence to phenomena in the Sun, cannot be answered at the moment, because our model is still quite unrealistic in many respects. For example in the Sun, and change with depth, which is not currently taken into account in DNS. Also, of course, the stratification is not isothermal, but convectively unstable. However, DNS in turbulent convection by Käpylä et al. (2011) have shown that still has a negative minimum in that case.

With regards to the production of sunspots, it is likely that NEMPI will shut off before the magnetic energy density has reached values comparable with the internal energy of the gas, as is the case in sunspots. Thus, some other mechanism is still needed to push the field of flux concentrations into that regime. One likely candidate is the mechanism of KM where the suppression of convective heat flux by the magnetic field is crucial. This impression is further justified by resent calculations of Stein et al. (2011) where pores are seen to form spontaneously in a simulation where horizontal magnetic fields are injected at the bottom of the domain.

Pores are small sunspots, so something else is needed to make these structures bigger and to amplify this mechanism further. Again, the answer could be related to scale separation. Thus, it will now be important to undertake detailed investigations of instabilities in strongly stratified layers with finite heat flux and finite magnetic field.

### Acknowledgements

We acknowledge the NORDITA dynamo programs of 2009 and 2011 for providing a stimulating scientific atmosphere. Computing resources provided by the Swedish National Allocations Committee at the Center for Parallel Computers at the Royal Institute of Technology in Stockholm and the High Performance Computing Center North in Umeå. This work was supported in part by the European Research Council under the AstroDyn Research Project No. 227952.

## References

- Brandenburg (2005) Brandenburg, A. 2005, ApJ, 625, 539
- Brandenburg et al. (1996) Brandenburg, A., Jennings, R. L., Nordlund, Å., Rieutord, M., Stein, R. F., & Tuominen, I. 1996, J. Fluid Mech., 306, 325
- Brandenburg et al. (2011a) Brandenburg, A., Kemel, K., Kleeorin, N., Mitra, D., & Rogachevskii, I. 2011a, ApJ, 740, L50 (BKKMR)
- Brandenburg et al. (2011b) Brandenburg, A., Kemel, K., Kleeorin, N., & Rogachevskii, I. 2011b, ApJ, submitted (arXiv:1005.5700)
- Brandenburg et al. (2010) Brandenburg, A., Kleeorin, N., & Rogachevskii, I. 2010, Astron. Nachr., 331, 5
- Brandenburg et al. (1995) Brandenburg, A., Procaccia, I., & Segel, D. 1995, Phys. Plasmas, 2, 1148
- Cheung et al. (2010) Cheung, M. C. M., Rempel, M., Title, A. M., & Schüssler, M. 2010, ApJ, 720, 233
- Cline et al. (2003) Cline, K. S., Brummell, N. H., & Cattaneo, F. 2003, ApJ, 599, 1449
- Golub et al. (1981) Golub, L., Rosner, R., Vaiana, G. S., & Weiss, N. O. 1981, ApJ, 243, 309
- Guerrero & Käpylä (2011) Guerrero, G., & Käpylä, P. J. 2011, A&A, 533, A40
- Käpylä et al. (2011) Käpylä, P. J., Brandenburg, A., Kleeorin, N., Mantere, M. J., & Rogachevskii, I. 2011, MNRAS, submitted (arXiv:1105.5785)
- Kemel et al. (2012) Kemel, K., Brandenburg, A., Kleeorin, N., & Rogachevskii, I. 2012, Astron. Nachr., in press (arXiv:1107.2752)
- Kitchatinov & Mazur (2000) Kitchatinov, L.L., & Mazur, M.V. 2000, Solar Phys., 191, 325 (KM)
- Kitchatinov & Olemskoy (2006) Kitchatinov, L. L., & Olemskoy, S. V. 2006, Astron. Lett., 32, 320
- Kitiashvili et al. (2010) Kitiashvili, I. N., Kosovichev, A. G., Wray, A. A., & Mansour, N. N. 2010, ApJ, 719, 307
- Kleeorin et al. (1996) Kleeorin, N., Mond, M., & Rogachevskii, I. 1996, A&A, 307, 293
- Kleeorin & Rogachevskii (1994) Kleeorin, N., & Rogachevskii, I. 1994, Phys. Rev. E, 50, 2716
- Kleeorin et al. (1990) Kleeorin, N.I., Rogachevskii, I.V., & Ruzmaikin, A.A. 1990, Sov. Phys. JETP, 70, 878
- Newcomb (1961) Newcomb, W.A. 1961, Phys. Fluids, 4, 391
- Nordlund et al. (1992) Nordlund, Å., Brandenburg, A., Jennings, R. L., Rieutord, M., Ruokolainen, J., Stein, R. F., & Tuominen, I. 1992, ApJ, 392, 647
- Ossendrijver (2003) Ossendrijver, M. 2003 Astron. Astrophys. Rev. 11, 287
- Parker (1966) Parker, E.N. 1966, ApJ, 145, 811
- Parker (1979) Parker, E. N. 1979, Cosmical magnetic fields (Oxford University Press, New York)
- Priest (1982) Priest, E.R. 1982, Solar Magnetohydrodynamics (D. Reidel Publ. Co., Dordrecht)
- Rempel (2011a) Rempel, M. 2011a, ApJ, 729, 5
- Rempel (2011b) Rempel, M. 2011b, ApJ, 740, 15
- Rempel et al. (2009a) Rempel, M., SchÃ¼ssler, M., & Knölker, M. 2009a, ApJ, 691, 640
- Rempel et al. (2009b) Rempel, M., Schüssler, M., Cameron, R. H., & Knölker, M. 2009b, Science, 325, 171
- Rogachevskii & Kleeorin (2007) Rogachevskii, I., & Kleeorin, N. 2007, Phys. Rev. E, 76, 056307
- Rüdiger (1989) Rüdiger, G. 1989, Differential rotation and stellar convection: Sun and solar-type stars (Gordon & Breach, New York)
- Rüdiger et al. (2011) Rüdiger, G., Kitchatinov, L. L., & Brandenburg, A. 2011, Solar Phys., 269, 3
- Stein et al. (2011) Stein, R. F., Lagerfjärd, A., Nordlund, Å., & Georgobiani, D. 2011, Solar Phys., 268, 271
- Tao et al. (1998) Tao, L., Weiss, N. O., Brownjohn, D. P., & Proctor, M. R. E. 1998, ApJ, 496, L39
- Tserkovnikov (1960) Tserkovnikov, Y. A. 1960, Sov. Phys. Dokl., 5, 87