Magnetic flux concentrations from dynamo-generated fields

# Magnetic flux concentrations from dynamo-generated fields

S. Jabbari Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden Department of Astronomy, AlbaNova University Center, Stockholm University, 10691 Stockholm, Sweden    A. Brandenburg Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden Department of Astronomy, AlbaNova University Center, Stockholm University, 10691 Stockholm, Sweden    I. R. Losada Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden Department of Astronomy, AlbaNova University Center, Stockholm University, 10691 Stockholm, Sweden    N. Kleeorin Department of Mechanical Engineering, Ben-Gurion University of the Negev, POB 653, Beer-Sheva 84105, Israel Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden Department of Radio Physics, N. I. Lobachevsky State University of Nizhny Novgorod, Russia    I. Rogachevskii Department of Mechanical Engineering, Ben-Gurion University of the Negev, POB 653, Beer-Sheva 84105, Israel Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden Department of Radio Physics, N. I. Lobachevsky State University of Nizhny Novgorod, Russia
July 8, 2019,  Revision: 1.148
###### Key Words.:
Sun: sunspots – Sun: dynamo – turbulence – magnetohydrodynamics (MHD) – hydrodynamics
###### Abstract

Context: The mean-field theory of magnetized stellar convection gives rise to two distinct instabilities: the large-scale dynamo instability, operating in the bulk of the convection zone and a negative effective magnetic pressure instability (NEMPI) operating in the strongly stratified surface layers. The latter might be important in connection with magnetic spot formation. However, the growth rate of NEMPI is suppressed with increasing rotation rates. On the other hand, recent direct numerical simulations (DNS) have shown a subsequent increase in the growth rate.

Aims: We examine quantitatively whether this increase in the growth rate of NEMPI can be explained by an mean-field dynamo, and whether both NEMPI and the dynamo instability can operate at the same time.

Methods: We use both DNS and mean-field simulations (MFS) to solve the underlying equations numerically either with or without an imposed horizontal field. We use the test-field method to compute relevant dynamo coefficients.

Results: DNS show that magnetic flux concentrations are still possible up to rotation rates above which the large-scale dynamo effect produces mean magnetic fields. The resulting DNS growth rates are quantitatively reproduced with MFS. As expected for weak or vanishing rotation, the growth rate of NEMPI increases with increasing gravity, but there is a correction term for strong gravity and large turbulent magnetic diffusivity.

Conclusions: Magnetic flux concentrations are still possible for rotation rates above which dynamo action takes over. For the solar rotation rate, the corresponding turbulent turnover time is about 5 hours, with dynamo action commencing in the layers beneath.

## 1 Introduction

The appearance of surface magnetic fields in the Sun presents some peculiar characteristics, such as being strongly concentrated into discrete spots. The origin and depth of such magnetic flux concentrations has long been the subject of considerable speculation. A leading theory by Parker (1955) interprets the emergence of such spots as the result of magnetically buoyant flux tubes at a depth of some 20 Mm. This magnetic field must be the result of a dynamo, but magnetic buoyancy also leads to the buoyant rise and subsequent loss of those magnetic structures. It was therefore thought that the dynamo should operate mainly at or even below the bottom of the convection zone where magnetic buoyancy could be stabilized by a subadiabatic temperature gradient (Parker, 1975). This led eventually to the idea that sunspots might be a direct consequence of dynamo-generated flux tubes that rise all the way from the bottom of the convection zone to the surface (e.g., Caligari et al., 1995). However, Schüssler (1980, 1983) emphasized early on that such fields would easily lose their systematic east–west orientation while ascending through the turbulent convection zone. D’Silva & Choudhuri (1993) estimated that a magnetic field strength of about would be needed to preserve the overall east–west orientation (Hale et al., 1919) and also to produce the observed tilt angle of active regions known as Joy’s law.

A great deal of effort has gone into determining the conditions under which magnetic flux ropes may or may not be able to rise buoyantly across the convection zone. Emonet et al. (1998) determined for the first time the basic minimum twist thresholds for the survival of twisted magnetic flux ropes during the rise. Subsequent studies were based on different types of numerical simulations, which tested the underlying hypotheses and looked for other effects, such as the robustness against background convective motions (Jouve et al., 2013) and magnetic flux erosion by reconnection with the background dynamo field (Pinto et al., 2013). These studies, as well as many others (see, e.g., Fan, 2008, 2009, and references therein) specifically looked at which flux-rope configurations are able to reproduce the observed emergent polarity tilt angles (Joy’s law).

The observed variation in the number of sunspots in time and latitude is expected to be linked to some kind of large-scale dynamo, as was modeled by Leighton et al. (1969) and Steenbeck & Krause (1969) long ago. This led Schüssler (1980) to propose a so-called flux-tube dynamo approach that would couple the buoyant rise of thin flux tubes to their regeneration. However, even today the connection between dynamos and flux tubes is done by hand (see, e.g., Choudhuri et al., 2007; Miesch & Dikpati, 2014), which means that an ad hoc procedure is invoked to link flux tube emergence to a mean-field dynamo. Of course, such tubes, or at least bipolar regions, should ultimately emerge from a sufficiently well-resolved and realistic simulation of solar convection. While global convective dynamo simulations of Nelson et al. (2011, 2013, 2014) show magnetically buoyant magnetic flux tubes of field strength, they do not yet address bipolar region formation. Indeed, solar surface simulations of Cheung et al. (2010) and Rempel & Cheung (2014) demonstrate that bipolar spots do form once a magnetic flux tube of field strength is injected at the bottom of their local domain ( below the surface). On the other hand, the deep solar simulations of Stein & Nordlund (2012) develop a bipolar active region with just magnetic field injected at the bottom of their domain. While these simulations taken together outline what might occur in the Sun, they do not necessarily support the description of spots as a direct result of thin flux tubes piercing the surface (e.g. Caligari et al., 1995).

A completely different suggestion is that sunspots develop locally at the solar surface, and that their east–west orientation would reflect the local orientation of the mean magnetic field close to the surface. The tilt angle would then be determined by latitudinal shear producing the observed orientation of the meridional component of the magnetic field (Brandenburg, 2005a). One of the possible mechanisms of local spot formation is the negative effective magnetic pressure instability (NEMPI; see Kleeorin et al., 1989, 1990; Kleeorin & Rogachevskii, 1994; Kleeorin et al., 1996; Rogachevskii & Kleeorin, 2007). Another potential mechanism of flux concentration is related to a thermo-magnetic instability in turbulence with radiative boundaries caused by the suppression of turbulent heat flux through the large-scale magnetic field (Kitchatinov & Mazur, 2000). The second instability has so far only been found in mean-field simulations (MFS), but not in direct numerical simulations (DNS) nor in large-eddy simulations (LES). By contrast, NEMPI has recently been found in DNS (Brandenburg et al., 2011) and LES (Brandenburg et al., 2014) of strongly stratified fully developed turbulence.

As demonstrated in earlier work (Brandenburg et al., 2013, 2014), NEMPI can lead to the formation of equipartition-strength magnetic spots, which are reminiscent of sunspots. Even bipolar spots can form in the presence of a horizontal magnetic field near the surface (see Warnecke et al., 2013). For this idea to be viable, NEMPI and the dynamo instability would need to operate in reasonable proximity to each other, so that the dynamo can supply the magnetic field that would be concentrated into spots, as was recently demonstrated by Mitra et al. (2014). In studying this process in detail, we have a chance of detecting new joint effects resulting from the two instabilities, which is one of the goals of the present paper. However, these two instabilities may also compete against each other, as was already noted by Losada et al. (2013). The large-scale dynamo effect relies on the combined presence of rotation and stratification, while NEMPI requires stratification and large enough scale separation. On the other hand, even a moderate amount of rotation suppresses NEMPI. In fact, Losada et al. (2012) found significant suppression of NEMPI when the Coriolis number is larger than about 0.03. Here, is the angular velocity and the turnover time of the turbulence, which is related to the rms velocity and the wavenumber of the energy-carrying eddies via . For the solar convection zone, the Coriolis number,

 \rm Co=2Ω/urmskf, (1)

varies from (at the surface using ) to (at the bottom of the convection zone using ). The value corresponds to a turnover time as short as two hours, which is the case at a depth of .

The strength of stratification, on the other hand, is quantified by the nondimensional parameter

 \rm Gr=g/c2skf≡(kfHρ)−1, (2)

where is the density scale height, is the sound speed, and is the gravitational acceleration. In the cases considered by Losada et al. (2012, 2013), the stratification parameter was , which is rather small compared with the estimated solar value of (see the conclusions of Losada et al., 2013). One can expect that larger values of Gr would result in correspondingly larger values of the maximum permissible value of Co, for which NEMPI is still excited, but this has not yet been investigated in detail.

The goal of the present paper is to study rotating stratified hydromagnetic turbulence in a parameter regime that we expect to be at the verge between NEMPI and dynamo instabilities. We do this by performing DNS and MFS. In MFS, the study of combined NEMPI and dynamo instability requires suitable parameterizations of the negative effective magnetic pressure and effects using suitable turbulent transport coefficients.

## 2 DNS study

We begin by reproducing some of the DNS results of Losada et al. (2013), who found the suppression of the growth rate of NEMPI with increasing values of Co and a subsequent enhancement at larger values, which they interpreted as being the result of dynamo action in the presence of an externally applied magnetic field. We also use DNS to determine independently the expected efficiency of the dynamo by estimating the effect from kinetic helicity measurements and by computing both effect and turbulent diffusivity directly using the test-field method (TFM).

### 2.1 Basic equations

In DNS of an isothermally stratified layer (Losada et al., 2013), we solve the equations for the velocity , the magnetic vector potential , and the density in the presence of rotation ,

 D\boldmathUDt = 1ρ\boldmathJ×B−2Ω×\boldmathU−ν\boldmathQ+\boldmathF+% \boldmathf, (3) ∂\boldmathA∂t = \boldmathU×B−η\boldmathJ, (4) ∂ρ∂t = −\boldmath∇⋅ρ\boldmathU, (5)

where is the advective derivative, is the angular velocity,

 \boldmathF=\boldmathg−c2s\boldmath∇lnρ (6)

determines the hydrostatic force balance, is the kinematic viscosity, is the magnetic diffusivity due to Spitzer conductivity of the plasma,

 −\boldmathQ = ∇2\boldmathU+\boldmath∇% \boldmath∇⋅\boldmathU/3+2\boldmathS\boldmath∇lnρ, (7) −\boldmathJ = ∇2\boldmathA−\boldmath∇% \boldmath∇⋅\boldmathA, (8)

are the modified vorticity and the current density, respectively, where the vacuum permeability has been set to unity,

 B=B0+\boldmath∇×\boldmathA (9)

is the total magnetic field, is the imposed uniform field, and

 Sij=12(∂jUi+∂iUj)−13δij\boldmath∇⋅\boldmath% U (10)

is the traceless rate-of-strain tensor. The forcing function consists of random, white-in-time, plane, nonpolarized waves with a certain average wavenumber . The turbulent rms velocity is approximately independent of with . The gravitational acceleration is chosen such that , so the density contrast between bottom and top is in a domain . Here, is the density scale height and is the smallest wavenumber that fits into the cubic domain of size . We adopt Cartesian coordinates , with periodic boundary conditions in the and directions and stress-free, perfectly conducting boundaries at top and bottom (). In most of the calculations, we use a scale separation ratio of 30, so is still the same as in earlier calculations. We use a fluid Reynolds number of 36, and a magnetic Prandtl number of 0.5. The magnetic Reynolds number is therefore . These values are a compromise between having both and Re large enough for NEMPI to develop at an affordable numerical resolution. The value of is specified in units of , where is the volume-averaged density, which is constant in time. The local equipartition field strength is . In our units, . However, time is specified as the turbulent-diffusive time , where is the estimated turbulent diffusivity. We also use DNS to compute these values more accurately with the TFM. The simulations are performed with the Pencil Code (http://pencil-code.googlecode.com), which uses sixth-order explicit finite differences in space and a third-order accurate time-stepping method. We use a numerical resolution of mesh points, which was found to be sufficient for the parameter regime specified above.

### 2.2 At the verge between NEMPI and dynamo

The work of Losada et al. (2013) suggested that for and , NEMPI becomes strongly suppressed, and that for still larger values, the growth rate increases again. This was tentatively associated with dynamo action, but it was not investigated in further detail. We now consider such a case with . This is a value that resulted in a rather low growth rate for NEMPI, while the estimated growth rate would be still subcritical for dynamo action. Following the work of Losada et al. (2013), we impose here a horizontal magnetic field in the direction with a strength of , which was previously found to be in the optimal range for NEMPI to develop (Kemel et al., 2012a).

To bring out the structures more clearly, it was found to be advantageous to present mean magnetic fields by averaging over the direction and over a certain time interval . We denote such averages by an overbar, e.g., . Once a dynamo develops, we expect a Beltrami-type magnetic field with phase shifted relative to by (Brandenburg, 2001). These are force-free fields with such as , for example.

Figure 1 shows visualizations of and together with the effective magnetic pressure, (defined below), at different times for a value of Co that is around the point where we expect onset of dynamo action. As in earlier work without rotation (Kemel et al., 2013), varies between 0 to . Furthermore, varies in the range . In Fig. 2, the extent of the domain is twice as big: . In Fig. 3 we show the result for , where a Beltrami-type field with a phase shift between and is well developed. For smaller values of Co, there are structures (e.g., for at and for at and ) that are reminiscent of those associated with NEMPI. This can be seen by comparing our Fig. 1 with Fig. 4 of Kemel et al. (2013) or Fig. 3 of Losada et al. (2013). When the domain is twice as wide, the number of structures simply doubles. A similar phenomenon was also seen in the simulations of Kemel et al. (2012b). For larger values of Co, NEMPI is suppressed and the dynamo, which generates mean magnetic field of a Beltrami-type structure, becomes more strongly excited.

The effective magnetic pressure shown in Figs. 13 is estimated by computing the component of the total stress from the fluctuating velocity and magnetic fields as

 Δ¯¯¯¯Πfxx=¯¯¯ρ(¯¯¯¯¯¯u2x−¯¯¯¯¯¯¯¯u20x)+12(¯¯¯¯¯b2−¯¯¯¯¯b20)−(¯¯¯¯¯b2x−¯¯¯¯¯¯¯b20x), (11)

where the subscript refers to the case with . We then calculate (Brandenburg et al., 2012a)

 qp=−2Δ¯¯¯¯Πfxx/¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathB2. (12)

Here, is a function of . We then calculate , which is the effective magnetic pressure divided by . We note that shows a systematic dependence and is negative in the upper part. Variations in the direction are comparatively weak and therefore do not show a clear correspondence with the horizontal variations of .

As in earlier work (Brandenburg et al., 2011), we characterize the strength of resulting structures by an amplitude of a suitably low wavenumber Fourier mode ( or 2), which is based on the magnetic field in the upper part of the domain, . In Fig. 4 we compare the evolution of for runs with different values of Co. For comparison, we also reproduce the first few runs of Losada et al. (2013) for , where we used in all cases. It turns out that for the new cases with and , the growth of is not as strong as for the cases with smaller Co. Furthermore, as is also evident from Figs. 1 and 2, the structures are now characterized by , while for they are still characterized by . The growth for all three cases (, both for normal and wider domains, as well as ) is similar. However, given that the typical NEMPI structures are not clearly seen for , it is possible that the growth of structures is simply overwhelmed by the much stronger growth due to the dynamo, which is not reflected in the growth of , whose growth is still mainly indicative of NEMPI. In this sense, there is some evidence of the occurrence of NEMPI in both cases.

### 2.3 Kinetic helicity

To estimate the effect and study its relation to kinetic helicity we begin by considering a fixed value of Gr equal to that used by Losada et al. (2013) and vary Co. For small values of Co, their data agreed with the MFS of Losada et al. (2012). For faster rotation, one eventually crosses the dynamo threshold. This is also the point at which the growth rate begins to increase again, although it now belongs to a different instability than for small values of Co. The underlying mechanism is the -dynamo, which is characterized by the dynamo number

 Cα=α/ηTk1, (13)

where is the typical value of the effect (here assumed spatially constant), is the sum of turbulent and microphysical magnetic diffusivities, and is the lowest wavenumber of the magnetic field that can be fitted into the domain. For isotropic turbulence, and are respectively proportional to the negative kinetic helicity and the mean squared velocity (Moffatt, 1978; Krause & Rädler, 1980; Rädler et al., 2003; Kleeorin & Rogachevskii, 2003)

 α≈α0≡−13τ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯% \boldmathω⋅\boldmathu,ηt≈ηt0≡13τ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathu2, (14)

where , so that (Blackman & Brandenburg, 2002; Candelaresi & Brandenburg, 2013)

 Cα=−ϵkϵfkf/k1. (15)

Here, is a free parameter characterizing possible dependencies on the forcing wavenumber, and is a measure for the relative kinetic helicity. Simulations of Brandenburg et al. (2012b) and Losada et al. (2013) showed that

 (16)

where is yet another non-dimensional parameter on the order of unity that may depend weakly on the scale separation ratio, , and is slightly different with and without imposed field. In the absence of an imposed field, Brandenburg et al. (2012b) found using . However, both an imposed field and a larger value of lead to a slightly increased value of . Our results are summarized in Fig. 5 for cases with and without imposed magnetic fields. Error bars are estimated as the largest departure of any one third of the full time series. The relevant points of Losada et al. (2013) give . For , the results of Brandenburg et al. (2012b) show a maximum with a subsequent decline of with increasing values of Co. However, although it is possible that the position of this maximum may be different for other values of Gr, it is unlikely to be relevant to our present study where we focus on smaller values of near dynamo onset. Thus, in conclusion, Eq. (16) seems to be a useful approximation that has now been verified over a range of different values of .

### 2.4 Test-field results

Our estimate for is based on the reference values and that are defined in Eq. (14) and represent approximations obtained from earlier simulations of helically forced turbulence (Sur et al., 2008). In the present study, helicity is self-consistently generated from the interaction between rotation and stratification. As an independent way of computing and , we now use the test-field method (TFM). It consists of solving auxiliary equations describing the evolution of magnetic fluctuations, , resulting from a set of several prescribed mean or test fields, . We solve for the corresponding vector potential with ,

 ∂\boldmathapq∂t=\boldmathu×¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathBpq+¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathU×bpq+(\boldmathu×bpq)′+η∇2\boldmathapq, (17)

where is the fluctuating part of the electromotive force and

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathBic=^\boldmathxicoskz,¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathBis=^\boldmathxisinkz,i=1,2, (18)

are the four test fields, which can show a cosine or sine variation with , while and are unit vectors in the two horizontal coordinate directions. The resulting are used to compute the electromotive force, , which is then expressed in terms of and as

 ¯¯¯Epqi=αij¯¯¯¯Bpqj−ηij¯¯¯¯Jpqj. (19)

By doing this for all four test field vectors, the and components of each of them gives eight equations for the eight unknowns, , , …, (for details see Brandenburg, 2005b).

With the TFM, we obtain the kernels and , from which we compute

 α=12(α11+α22),ηt=12(η11+η22), (20)
 γ=12(α21−α12),δ=12(η21−η12). (21)

We normalize and by their respective values obtained for large magnetic Reynolds numbers defined in Eq. (14), and denote them by a tilde, i.e., and . We use the latter normalization also for , i.e., , but expect its value to vanish in the limit of zero angular velocity. No standard turbulent pumping velocity is expected (Krause & Rädler, 1980; Moffatt, 1978), because the rms turbulent velocity is independent of height. However, this is not quite true. To show this, we normalize by and present . In our normalization, the molecular value is given by

.

We consider test fields that are constant in time and vary sinusoidally in the direction. We choose certain values of between and (= and also vary the value of Co between 0 and about 1.06 while keeping fixed. In all cases where the scale separation ratio is held fixed, we used , which is larger than what has been used in earlier studies (Brandenburg et al., 2008b), where was typically 5.

In Figure 6 we show the dependence of the coefficients on the normalized wavenumber of the test field, . The three coefficients , , and show the same behavior of the form of

 ~σ=~σ0/(1+ℓ2σk2) (22)

for , , or , while for we use

 ~γ=~γ0+~γ2ℓ2γk2/(1+ℓ2γk2), (23)

where , , and . These results have been obtained for and . Again, error bars are estimated as the largest departure of any one third of the full time series.

Most of the coefficients are only weakly dependent on the value of Co, except and . The former varies approximately as

 ~γ=~γ0+~γ2Ω\rm Co2, (24)

where and . Here and in the following, we keep . For the same value of , the functional form for shows a linear increase with Co, i.e., where . Figure 7 shows that is nearly independent of the Coriolis number. This result is in agreement with that obtained Kleeorin & Rogachevskii (2003), where a theory of versus Coriolis number was developed for large fluid and magnetic Reynolds numbers. It turns out that the new values of and that have been obtained now with the TFM are somewhat different from previous TFM studies that originally estimated ( and ). The TFM results now suggest in Eq. (15). The reason for the departure from unity cannot just be the fact that helicity is now self-consistently generated, because this was also the case in the earlier work of Brandenburg et al. (2012b). The only plausible reason is the large value of that is now much larger than before ( compared to in most previous studies), which explains the reason for our choice of the subscript in .

The origin of weak pumping found in Figs. 6 and 7 is unclear. For a weak mean magnetic field, pumping of the magnetic field can cause not only inhomogeneous distributions of the velocity fluctuations (Krause & Rädler, 1980; Moffatt, 1978) or magnetic fluctuations (Rädler et al., 2003), but also non-uniform distribution of the fluid density in the presence either of small-scale dynamo or turbulent convection (Rogachevskii & Kleeorin, 2006). In our simulations there is no small-scale dynamo effect, because is too low. There is also no turbulent convection possible in our setup. The pumping effect is also not connected with nonlinear effects; see Fig. 2 in Rogachevskii & Kleeorin (2004).

## 3 MFS study

We now want to see whether the suppression of NEMPI and the subsequent increase in the resulting growth rate can be reproduced in MFS. In addition to a parameterization for the negative effective magnetic pressure in the momentum equation, we add one for the electromotive force. The important terms here are the effect and the turbulent magnetic diffusivity, whose combined effect is captured by the quantity , which is defined in Eq. (13) and related to DNS parameters in Eq. (15). In contrast to DNS, the advantage of MFS is that they can more easily be extended to astrophysically interesting conditions of large Reynolds numbers and more complex geometries.

### 3.1 The model

Our MFS model is in many ways the same as that of Jabbari et al. (2013), where parameterizations for negative effective magnetic pressure and electromotive force where, for the first time, considered in combination with each other. Their calculations were performed in spherical shells without Coriolis force, while here we apply instead Cartesian geometry and do include the Coriolis force. The evolution equations for mean velocity , mean vector potential , and mean density , are thus

 D¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathUDt = 1¯¯¯ρ⎛⎝¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathJ×¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathB+\boldmath∇qp¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathB22⎞⎠−2Ω×¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathU−νT¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathQ+¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathF, (25) ∂¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathA∂t = ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathU×¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmath% B+α¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathB−ηT¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯% \boldmathJ, (26) D¯¯¯ρDt = −¯¯¯ρ\boldmath∇⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯% \boldmathU,

where is the advective derivative,

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathF=\boldmathg−c2s% \boldmath∇ln¯¯¯ρ (27)

is the mean-field hydrostatic force balance, and are the sums of turbulent and microphysical values of magnetic diffusivity and kinematic viscosities, respectively, is the aforementioned coefficient in the effect, is the mean current density,

 −¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathQ=∇2¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathU+13\boldmath∇\boldmath∇⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathU+2¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathS\boldmath∇ln¯¯¯ρ (28)

is a term appearing in the viscous force, where is the traceless rate of strain tensor of the mean flow with components , and finally determines the turbulent contribution to the mean Lorentz force. Here, depends on the local field strength and is approximated by (Kemel et al., 2012a)

 qp(β)=qp01+β2/β2p=β2⋆β2p+β2, (29)

where , , and are constants, is the normalized mean magnetic field, and is the equipartition field strength. For , Brandenburg et al. (2012a) found and . We use as our reference model the parameters for , also used by Losada et al. (2013), which yields

 βp=0.058,β⋆=0.33(reference model). (30)

In some cases we also compare with , which was found to match more closely the measured dependence of the effective magnetic pressure on by Losada et al. (2013). For vertical magnetic fields, MFS for a range of model parameters have been given by Brandenburg et al. (2014). In the MFS, we use (Sur et al., 2008)

 ηt≈ηt0≡urms/3kf (31)

to replace , so

 \rm Gr=3ηt/urmsHρ (32)

and (Losada et al., 2013)

 \rm Co=2Ω/urmskf=6Ωηt/u2rms. (33)

We now consider separately cases where we vary either Co or Gr. In addition, we also vary the scale separation ratio , which is essentially a measure of the inverse turbulent diffusivity, i.e.,

 kf/k1=urms/3ηtk1 (34)

(see Eq. (31)).

### 3.2 Fixed value of Gr

The work of Losada et al. (2012) has shown that the growth rate of NEMPI, , decreases with increasing values of the rotation rate. They found it advantageous to express in terms of the quantity

 λ∗0=β⋆urms/Hρ. (35)

As discussed above, the normalized growth rate shows first a decline with increasing values of Co, but then an increase for , which was argued to be a result of the dynamo effect (Losada et al., 2013). This curve has a minimum at . As rotation is increased further, the combined action of stratification and rotation leads to increased kinetic helicity and thus eventually to the onset of mean-field dynamo action.

Owing to the effects of turbulent diffusion, the actual value of the growth rate of NEMPI is always expected to be less than . Kemel et al. (2013) proposed an empirical formula replacing by , where is the wavenumber. This would lead to

 λ/λ∗0∝1−\rm Gr∗/\rm Gr, (36)

with a coefficient . However, as we will see below, this expression is not found to be consistent with our numerical data.

The onset of the dynamo instability is governed by the dynamo number

 Cα=ϵf0\rm Gr\rm Cokf/k1. (37)

For a cubic domain, large-scale dynamo action occurs for , which was confirmed by Losada et al. (2013), who found the typical Beltrami fields for two supercritical cases. They used the parameters and values of Co up to 0.6. Here we present MFS in two- and three-dimensional domains for the same values of Gr and a similar range of Co values. In Fig. 8, we compare the DNS of Losada et al. (2013) with our reference model defined through Eq. (30) and referred to as MFS(i) as well as with the case , referred to as MFS(ii).

### 3.3 Larger stratification, smaller scale separation

The expected theoretical maximum growth rate of NEMPI is given by Eq. (35). At zero rotation, we thus expect . To check this, we performed two-dimensional MFS in a squared domain of size . The result is shown in Fig. 9 for the model parameters given in Eq. (30). When Gr is small, we find that , which is below the expected value. As we increase Gr, decreases until NEMPI can no longer be detected for .

It is conceivable that this decrease may have been caused by the following two factors. First, the growth rate is expected to increase with Gr, but for fixed scale separation, the resulting density contrast becomes huge. Finite resolution might therefore have caused inaccuracies. Second, although the growth rate should not depend on (Kemel et al., 2012a), we need to make sure that the mode is fully contained within the domain. In other words, we are interested in the largest growth rate as we vary the value of . Again, to limit computational expense, we tried only a small number of runs, keeping the size of the domain the same. This may have caused additional uncertainties. However, it turns out that our results are independent of whether Gr is changed by changing or (). This suggests that our results for large values of shown in Fig. 9 may in fact be accurate. To illustrate this more clearly, we rewrite

 \rm Gr=3ηturmsHρ=3~ηtk1csurmsgc2s=3~ηt~g/\rm Ma, (38)

where we have defined

 ~ηt=ηtk1/cs,~g=g/c2sk1≡(k1Hρ)−1. (39)

Figure 9 shows that is indeed independent of the individual values of and as long as Gr is the same. For small values of and large diffusivity (), the velocity evolves in an oscillatory fashion with a rapid growth and a gradual subsequent decline. In Fig. 9, the isolated data point at reflects the speed of growth during the periodic rise phase, but it is unclear whether or not it is related to NEMPI.

In the inset, we plot versus itself. This shows that the growth rate (in units of the inverse turnover time) increases with increasing when is small. However, as follows from theoretical analysis, the growth rate decreases with increasing . When is larger (corresponding to smaller scale separation), the growth rate of NEMPI is reduced for the same value of and it decreases with when .

The decrease of with increasing values of Gr can be approximated by the formula

 λ/λ∗0≈0.3/[1+2\rm Gr+(4\rm Gr)2], (40)

which is shown in Fig. 9 as a dash-dotted line. This expression is qualitatively different from the earlier, more heuristic expression proposed by Kemel et al. (2013) where the dimensional growth rate was simply modified by an ad hoc diffusion term of the form . In that case, contrary to our MFS, the normalized growth rate would actually increase with increasing values of Gr (see Eq. (36)).

### 3.4 Co dependence at larger stratification

We consider the normalized growth rate of the combined NEMPI and dynamo instabilities as a function of Co for different values of Gr. As is clear from Fig. 9, using a fixed value of and varying gives us the possibility to increase Gr to larger values of up to 1. In the following we use this procedure to compare the behavior of the growth rate versus Co for three values of Gr, , , and (see Fig. 10). It can be seen that the behavior of the curves is independent of the values of Gr, but the points where the minima of the curves occur moves toward bigger values of Co as Gr increases. This also happens in the case when there is only dynamo action without imposed magnetic field (dashed lines in Fig. 10). One also sees that the increase of the growth rate with increasing Co is much stronger in the case of larger Gr (compare the lines for with those for 0.21 and 1). Finally, comparing runs with and without imposed magnetic field, but the same value of Gr, the growth rate of NEMPI is in most cases below that of the coupled system with NEMPI and dynamo instability.

In Fig. 10 we see that the dependence of on Gr is opposite for small and large values of Co. When , an increase in Gr leads to a decrease in (compare the line with that for 0.21 along a cut through in Fig. 10), while for , an increase in Gr leads to an increase in (compare all three lines in Fig. 10 along a cut through ). The second case is caused by the increase of the dynamo number , which is directly proportional to Gr (see Eq. (37)). On the other hand, for small values of Co, only NEMPI operates, but if Gr in Eq. (38) is increased by increasing rather than , the dynamo is suppressed by enhanced turbulent diffusion (see also Fig. 9). This is related to the fact that the properties of the system depend not just on Gr and Co, but also on or , which is proportional to all three parameters (see Eq. (37)).

## 4 Discussion and conclusions

The present work has brought us one step closer to being able to determine whether the observable solar activity such as sunspots and active regions could be the result of surface effects associated with strong stratification. A particularly important aspect has been the interaction with a dynamo process that must ultimately be responsible for generating the overall magnetic field. Recent global convective dynamo simulations of Nelson et al. (2011, 2013, 2014) have demonstrated that flux tubes with field strength can be produced in the solar convection zone. This is almost as strong as the magnetic flux tubes anticipated from earlier investigations of rising flux tubes requiring them to not break up and to preserve their east–west orientation (D’Silva & Choudhuri, 1993). Would we then still need surface effects such as NEMPI to produce sunspots? The answer might well be yes, because the flux ropes that have been isolated in the visualizations of Nelson et al. (2011, 2013, 2014) appear to have cross sections that are much larger than sunspots at the solar surface. Further concentration into thinner tubes would be required if they were to explain sunspots by just letting them pierce the surface.

Realistic hydromagnetic simulations of the solar surface are now beginning to demonstrate that fields at a depth of can produce sunspot-like appearances at the surface (Rempel & Cheung, 2014). However, we have to ask about the physical process contributing to this phenomenon. A purely descriptive analysis of simulation data cannot replace the need for a more prognostic approach that tries to reproduce the essential physics using simpler models. Although Rempel & Cheung (2014) propose a mechanism involving mean-field terms in the induction equation, they do not show that their model equations can actually describe the process of magnetic flux concentration. In fact, their description is somewhat reminiscent of flux expulsion, which was invoked earlier by Tao et al. (1998) to explain the segregation of magneto-convection into magnetized and unmagnetized regions. In this context, NEMPI provides such an approach that can be used prognostically rather than diagnostically. However, this approach has problems of its own, some of which are addressed in the present work. Does NEMPI stop working when ? How does it interact with the underlying dynamo? Such a dynamo is believed to control the overall sunspot number and the concentration of sunspots to low latitudes.

Our new DNS suggest that, although rotation tends to suppress NEMPI, magnetic flux concentrations can still form at Coriolis numbers of . This is slightly larger than what was previously found from MFS both with horizontal and vertical magnetic fields and the same value of Gr. For the solar rotation rate of , a value of corresponds to , which is longer than the earlier MFS values of for a horizontal field (Losada et al., 2013) and for a vertical field (Brandenburg et al., 2014).

Using the TFM, we have confirmed earlier findings regarding and , although for our new simulations both coefficients are somewhat larger, which is presumably due to the larger scale separation. The ratio between and determines the dynamo number and is now about 40% below previous estimates. There is no evidence of other important mean-field effects that could change our conclusion about a cross-over from suppressed NEMPI to increased dynamo activity. We now confirm quantitatively that the enhanced growth past the initial suppression of NEMPI is indeed caused by mean-field dynamo action in the presence of a weak magnetic field. The position of the minimum in the growth rate coincides with the onset of mean-field dynamo action that takes the effect into account.

For weak or no rotation, we find that the normalized NEMPI growth rate is described by a single parameter Gr, which is proportional to the product of gravity and turbulent diffusivity, where the latter is a measure of the inverse scale separation ratio. This normalization takes into account that the growth rate increases with increasing gravity. The growth rate compensated in this way shows a decrease with increasing gravity and turbulent diffusivity that is different from an earlier, more heuristic, expression proposed by Kemel et al. (2013). The reason for this departure is not quite clear. One possibility is some kind of gravitational quenching, because the suppression is well described by a quenching factor that becomes important when Gr exceeds a value of around 0.5. This quenching is probably not important for stellar convection where the estimated value of Gr is 0.17 (Losada et al., 2013). It might, however, help explain mismatches with the expected theoretical growth rate that was found to be proportional to Gr (Kemel et al., 2013) and that was determined from recent DNS (Brandenburg et al., 2014).

An important question is whether NEMPI will really be strong enough to produce sunspots with super-equipartition strength. It has always been clear that NEMPI can only work for a magnetic field strength that is a small fraction of the local equipartition field value. However, super-equipartition fields are produced if the magnetic field is vertical (Brandenburg et al., 2013). Subsequent work showed quantitatively that NEMPI does indeed work at subequipartition field strengths, but since mass flows mainly along magnetic field lines, the reduced pressure leads to suction which tends to evacuate the upper parts of the tube (Brandenburg et al., 2014). This is similar to the “hydraulic effect” envisaged by Parker (1976), who predicted such downflows along flux tubes. In a later paper Parker (1978), gives more realistic estimates, but the source of downward flows remained unclear. Meanwhile, the flux emergence simulations of Rempel & Cheung (2014) show at first upflows in their magnetic spots (see their Fig. 5), but as the spots mature, a downflow develops (see their Fig. 7). In their case, because they have convection, those downflows can also be ascribed to supergranular downflows, as was done by Stein & Nordlund (2012). Nevertheless, in the isothermal simulations of Brandenburg et al. (2013, 2014), this explanation would not apply. Thus, we now know that the required downflows can be caused by NEMPI, but we do not know whether this is also what happens in the Sun.

Coming back to our paper, where NEMPI is coupled to a dynamo, the recent work of Mitra et al. (2014) is relevant because it shows that intense bipolar spots can be generated in an isothermal simulation with strongly stratified non-helically driven turbulence in the upper part and a helical dynamo in the lower part. The resulting surface structure resembles so-called spots that have previously only been found in the presence of strongly twisted and kink-unstable flux tubes (Linton et al., 1998). While the detailed mechanism of this work is not yet understood, it reminds us that it is too early to draw strong conclusions about NEMPI as long as not all its aspects have been explored in sufficient detail.

###### Acknowledgements.
This work was supported in part by the European Research Council under the AstroDyn Research Project No. 227952, the Swedish Research Council under the grants 621-2011-5076 and 2012-5797, the Research Council of Norway under the FRINATEK grant 231444 (AB, IR), and a grant from the Government of the Russian Federation under contract No. 11.G34.31.0048 (NK, IR). We acknowledge the allocation of computing resources provided by the Swedish National Allocations Committee at the Center for Parallel Computers at the Royal Institute of Technology in Stockholm and the Nordic Supercomputer Center in Reykjavik.

## References

• Blackman & Brandenburg (2002) Blackman, E. G., & Brandenburg, A. 2002, ApJ, 579, 359
• Brandenburg (2001) Brandenburg, A. 2001, ApJ, 550, 824
• Brandenburg (2005a) Brandenburg, A. 2005a, ApJ, 625, 539
• Brandenburg (2005b) Brandenburg, A. 2005b, Astron. Nachr., 326, 787
• Brandenburg et al. (2014) Brandenburg, A., Gressel, O., Jabbari, S., Kleeorin, N., & Rogachevskii, I. 2014, A&A, 562, A53
• Brandenburg et al. (2011) Brandenburg, A., Kemel, K., Kleeorin, N., Mitra, D., & Rogachevskii, I. 2011, ApJ, 740, L50
• Brandenburg et al. (2012a) Brandenburg, A., Kemel, K., Kleeorin, N., & Rogachevskii, I. 2012a, ApJ, 749, 179
• Brandenburg et al. (2013) Brandenburg, A., Kleeorin, N., & Rogachevskii, I. 2013, ApJ, 776, L23
• Brandenburg et al. (2012b) Brandenburg, A., Rädler, K.-H., & Kemel, K. 2012b, A&A, 539, A35
• Brandenburg et al. (2008a) Brandenburg, A., Rädler, K.-H., Rheinhardt, M., & Käpylä, P. J. 2008a, ApJ, 676, 740
• Brandenburg et al. (2008b) Brandenburg, A., Rädler, K.-H., & Schrinner, M. 2008b, A&A, 482, 739
• Caligari et al. (1995) Caligari, P., Moreno-Insertis, F., & Schüssler, M. 1995, ApJ, 441, 886
• Candelaresi & Brandenburg (2013) Candelaresi, S., & Brandenburg, A. 2013, Phys. Rev. E, 87, 043104
• Cheung et al. (2010) Cheung, M. C. M., Rempel, M., Title, A. M., & Schüssler, M. 2010, ApJ, 720, 233
• Choudhuri et al. (2007) Choudhuri, A. R., Chatterjee, P., & Jiang, J. 2007, Phys. Rev. Lett., 98, 131103
• D’Silva & Choudhuri (1993) D’Silva, S., & Choudhuri, A. R. 1993, A&A, 272, 621
• Emonet et al. (1998) Emonet, T., & Moreno-Insertis, F. 1998, ApJ, 492, 804
• Fan (2008) Fan, Y. 2008, ApJ, 676, 680
• Fan (2009) Fan, Y. 2009, Living Rev. Solar Phys., 6, 4
• Hale et al. (1919) Hale, G. E., Ellerman, F., Nicholson, S. B., & Joy, A. H. 1919, ApJ, 49, 153
• Jabbari et al. (2013) Jabbari, S., Brandenburg, A., Kleeorin, N., Mitra, D., & Rogachevskii, I. 2013, A&A, 556, A106
• Jouve et al. (2013) Jouve, L., Brun, A. S., & Aulanier, G. 2013, ApJ, 762, 23
• Kemel et al. (2012a) Kemel, K., Brandenburg, A., Kleeorin, N., & Rogachevskii, I. 2012a, Astron. Nachr., 333, 95
• Kemel et al. (2012b) Kemel, K., Brandenburg, A., Kleeorin, N., Mitra, D., & Rogachevskii, I. 2012b, Solar Phys., 280, 321
• Kemel et al. (2013) Kemel, K., Brandenburg, A., Kleeorin, N., Mitra, D., & Rogachevskii, I. 2013, Solar Phys., 287, 293
• Kitchatinov & Mazur (2000) Kitchatinov, L. L. & Mazur, M. V. 2000, Solar Phys., 191, 325
• Kleeorin & Rogachevskii (1994) Kleeorin, N., & Rogachevskii, I. 1994, Phys. Rev. E, 50, 2716
• Kleeorin & Rogachevskii (2003) Kleeorin, N., & Rogachevskii, I. 2003, Phys. Rev. E, 67, 026321
• Kleeorin et al. (1996) Kleeorin, N., Mond, M., & Rogachevskii, I. 1996, A&A, 307, 293
• Kleeorin et al. (1989) Kleeorin, N.I., Rogachevskii, I.V., & Ruzmaikin, A.A. 1989, Sov. Astron. Lett., 15, 274
• Kleeorin et al. (1990) Kleeorin, N. I., Rogachevskii, I. V., Ruzmaikin, A. A. 1990, Sov. Phys. JETP, 70, 878
• Krause & Rädler (1980) Krause, F., & Rädler, K.-H. 1980, Mean-field Magnetohydrodynamics and Dynamo Theory (Oxford: Pergamon Press)
• Leighton et al. (1969) Leighton, R. B. 1969, ApJ, 156, 1
• Linton et al. (1998) Linton, M. G., Dahlburg, R. B., Fisher, G. H., & Longcope, D. W. 1998, ApJ, 507, 404
• Losada et al. (2012) Losada, I. R., Brandenburg, A., Kleeorin, N., Mitra, D., & Rogachevskii, I. 2012, A&A, 548, A49
• Losada et al. (2013) Losada, I. R., Brandenburg, A., Kleeorin, N., & Rogachevskii, I. 2013, A&A, 556, A83
• Miesch & Dikpati (2014) Miesch, M. S., & Dikpati, M. 2014, ApJ, 785, L8
• Mitra et al. (2014) Mitra, D., Brandenburg, A., Kleeorin, N., & Rogachevskii, I. 2014, MNRAS, submitted, arXiv:1404.3194
• Moffatt (1978) Moffatt, H.K. 1978, Magnetic Field Generation in Electrically Conducting Fluids (Cambridge: Cambridge Univ. Press)
• Nelson et al. (2011) Nelson, N. J., Brown, B. P., Brun, A. S., Miesch, M. S., & Toomre, J. 2011, ApJ, 739, L38
• Nelson et al. (2013) Nelson, N. J., Brown, B. P., Brun, A. S., Miesch, M. S., & Toomre, J. 2013, ApJ, 762, 73
• Nelson et al. (2014) Nelson, N. J., Brown, B. P., Brun, A. S., Miesch, M. S., & Toomre, J. 2014, Solar Phys., 289, 441
• Parker (1955) Parker, E. N. 1955, ApJ, 121, 491
• Parker (1975) Parker, E. N. 1975, ApJ, 198, 205
• Parker (1976) Parker, E. N. 1976, ApJ, 210, 816
• Parker (1978) Parker, E. N. 1978, ApJ, 221, 368
• Pinto et al. (2013) Pinto, R. F., & Brun, A. S. 2013, ApJ, 772, 55
• Rädler et al. (2003) Rädler, K.-H., Kleeorin, N., & Rogachevskii, I. 2003, Geophys. Astrophys. Fluid Dyn., 97, 249
• Rempel & Cheung (2014) Rempel, M., & Cheung, M. C. M. 2014, ApJ, 785, 90
• Roberts & Soward (1975) Roberts, P. H. & Soward, A. M. 1975, Astron. Nachr., 296, 49
• Rogachevskii & Kleeorin (2004) Rogachevskii, I., & Kleeorin, N. 2004, Phys. Rev. E, 70, 046310
• Rogachevskii & Kleeorin (2006) Rogachevskii, I. & Kleeorin, N., 2006, Geophys. Astrophys. Fluid Dyn., 100, 243
• Rogachevskii & Kleeorin (2007) Rogachevskii, I., & Kleeorin, N. 2007, Phys. Rev. E, 76, 056307
• Schüssler (1980) Schüssler, M. 1980, Nature, 288, 150
• Schüssler (1983) Schüssler, M. 1983, in Solar and stellar magnetic fields: origins and coronal effects, ed. J. O. Stenflo (Reidel), 213
• Steenbeck & Krause (1969) Steenbeck, M., & Krause, F. 1969, Astron. Nachr., 291, 49
• Stein & Nordlund (2012) Stein, R. F., & Nordlund, Å. 2012, ApJ, 753, L13
• Sur et al. (2008) Sur, S., Brandenburg, A., & Subramanian, K. 2008, MNRAS, 385, L15
• Tao et al. (1998) Tao, L., Weiss, N. O., Brownjohn, D. P., & Proctor, M. R. E. 1998, ApJ, 496, L39
• Warnecke et al. (2013) Warnecke, J., Losada, I. R., Brandenburg, A., Kleeorin, N., & Rogachevskii, I. 2013, ApJ, 777, L37
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters