Free Access
Issue
A&A
Volume 568, August 2014
Article Number A112
Number of page(s) 11
Section The Sun
DOI https://doi.org/10.1051/0004-6361/201423499
Published online 02 September 2014

© ESO, 2014

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 100 kG 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 40 kG 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 10 kG field strength is injected at the bottom of their local domain (7.5 Mm below the surface). On the other hand, the deep solar simulations of Stein & Nordlund (2012) develop a bipolar active region with just 1 kG 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, 1996; Kleeorin & Rogachevskii 1994; 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 Co = 2Ωτ 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 urms and the wavenumber kf of the energy-carrying eddies via τ = (urmskf)-1. For the solar convection zone, the Coriolis number, (1)varies from 2 × 10-3 (at the surface using τ = 5min) to 5 (at the bottom of the convection zone using τ = 10 days). The value Co = 0.03 corresponds to a turnover time as short as two hours, which is the case at a depth of 10 Mm.

The strength of stratification, on the other hand, is quantified by the nondimensional parameter (2)where is the density scale height, cs is the sound speed, and g is the gravitational acceleration. In the cases considered by Losada et al. (2012, 2013), the stratification parameter was Gr = 0.03, which is rather small compared with the estimated solar value of Gr = 0.16 (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 U, the magnetic vector potential A, and the density ρ in the presence of rotation Ω, where D/Dt = /∂t + U· is the advective derivative, is the angular velocity, (6)determines the hydrostatic force balance, ν is the kinematic viscosity, η is the magnetic diffusivity due to Spitzer conductivity of the plasma, are the modified vorticity and the current density, respectively, where the vacuum permeability μ0 has been set to unity, (9)is the total magnetic field, B0 = (0,B0,0) is the imposed uniform field, and (10)is the traceless rate-of-strain tensor. The forcing function f consists of random, white-in-time, plane, nonpolarized waves with a certain average wavenumber kf. The turbulent rms velocity is approximately independent of z with . The gravitational acceleration g = (0,0, − g) is chosen such that k1Hρ = 1, so the density contrast between bottom and top is exp(2π) ≈ 535 in a domain πk1zπ. Here, is the density scale height and k1 = 2π/L is the smallest wavenumber that fits into the cubic domain of size L3. We adopt Cartesian coordinates (x,y,z), with periodic boundary conditions in the x and y directions and stress-free, perfectly conducting boundaries at top and bottom (z = ± Lz/ 2). In most of the calculations, we use a scale separation ratio kf/k1 of 30, so Gr = 0.03 is still the same as in earlier calculations. We use a fluid Reynolds number Re ≡ urms/νkf of 36, and a magnetic Prandtl number PrM = ν/η of 0.5. The magnetic Reynolds number is therefore ReM = PrMRe = 18. These values are a compromise between having both kf and Re large enough for NEMPI to develop at an affordable numerical resolution. The value of B0 is specified in units of , where ρ0 = ⟨ ρ is the volume-averaged density, which is constant in time. The local equipartition field strength is . In our units, k1 = cs = μ0 = ρ0 = 1. However, time is specified as the turbulent-diffusive time , where ηt0 = urms/ 3kf 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 code1, which uses sixth-order explicit finite differences in space and a third-order accurate time-stepping method. We use a numerical resolution of 2563 mesh points, which was found to be sufficient for the parameter regime specified above.

thumbnail Fig. 1

Visualization of and together with effective magnetic pressure for different times. Here Ω = 0.15, Co = 0.09, Gr = 0.033, and kf/k1 = 30.

Open with DEXTER

thumbnail Fig. 2

Like Fig. 1, but for a wider domain.

Open with DEXTER

thumbnail Fig. 3

Like Fig. 1, but for Ω = 0.35, so Co = 0.22.

Open with DEXTER

2.2. At the verge between NEMPI and dynamo

The work of Losada et al. (2013) suggested that for Gr = 0.03 and Co ≥ 0.03, 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 Co = 0.09. 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 y direction with a strength of 0.05Beq0, 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 y direction and over a certain time interval Δt. 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 π/ 2 (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 2B0. Furthermore, varies in the range ± 2B0. In Fig. 2, the x extent of the domain is twice as big: −2π<k1x< 2π. In Fig. 3 we show the result for Co = 0.22, where a Beltrami-type field with a π/ 2 phase shift between and is well developed. For smaller values of Co, there are structures (e.g., for t/τ = 1.8 at x/Hρ ≈ 1.5 and for t/τ = 2.4 at x/Hρ ≈ 1.5 and −2) 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 α2 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 xx component of the total stress from the fluctuating velocity and magnetic fields as (11)where the subscript 0 refers to the case with B0 = 0. We then calculate (Brandenburg et al. 2012a) (12)Here, qp(β) is a function of . We then calculate , which is the effective magnetic pressure divided by . We note that shows a systematic z dependence and is negative in the upper part. Variations in the x 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 Bk of a suitably low wavenumber Fourier mode (k/k1 = 1 or 2), which is based on the magnetic field in the upper part of the domain, 2 ≤ z/Hρπ. In Fig. 4 we compare the evolution of Bk/Beq0 for runs with different values of Co. For comparison, we also reproduce the first few runs of Losada et al. (2013) for Co = 0.0060.13, where we used k/k1 = 1 in all cases. It turns out that for the new cases with Co = 0.09 and 0.22, the growth of Bk/Beq0 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 k/k1 = 2, while for Co = 0.22 they are still characterized by k/k1 = 1. The growth for all three cases (Co = 0.09, both for normal and wider domains, as well as Co = 0.22) is similar. However, given that the typical NEMPI structures are not clearly seen for Co = 0.22, 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 Bk/Beq0, 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 α2-dynamo, which is characterized by the dynamo number (13)where α is the typical value of the α effect (here assumed spatially constant), ηT = ηt + η is the sum of turbulent and microphysical magnetic diffusivities, and k1 is the lowest wavenumber of the magnetic field that can be fitted into the domain. For isotropic turbulence, α and ηt 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) (14)where τ = (urmskf)-1, so that (Blackman & Brandenburg 2002; Candelaresi & Brandenburg 2013) (15)Here, ϵk is a free parameter characterizing possible dependencies on the forcing wavenumber, and ϵf is a measure for the relative kinetic helicity. Simulations of Brandenburg et al. (2012b) and Losada et al. (2013) showed that (16)where ϵf0 is yet another nondimensional parameter on the order of unity that may depend weakly on the scale separation ratio, kf/k1, and is slightly different with and without imposed field. In the absence of an imposed field, Brandenburg et al. (2012b) found ϵf0 ≈ 2 using kf/k1 = 5. However, both an imposed field and a larger value of kf/k1 lead to a slightly increased value of ϵf0. 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 ϵf0 ≈ 2.8. For Gr Co ≳ 0.5, the results of Brandenburg et al. (2012b) show a maximum with a subsequent decline of ϵf 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 Cα 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 kf/k1.

2.4. Test-field results

Our estimate for Cα is based on the reference values α0 and ηt0 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 ηt, we now use the test-field method (TFM). It consists of solving auxiliary equations describing the evolution of magnetic fluctuations, bpq, resulting from a set of several prescribed mean or test fields, . We solve for the corresponding vector potential with , (17)where is the fluctuating part of the electromotive force and (18)are the four test fields, which can show a cosine or sine variation with z, while and are unit vectors in the two horizontal coordinate directions. The resulting bpq are used to compute the electromotive force, , which is then expressed in terms of and as (19)By doing this for all four test field vectors, the x and y components of each of them gives eight equations for the eight unknowns, α11, α12, ..., η22 (for details see Brandenburg 2005b).

thumbnail Fig. 4

Comparison of the evolution of Bk/Beq0 for runs with different values of Co. In the first panel k/k1 = 1, while in the second panel k/k1 = 2 for the two runs with Co = 0.09 (label W refers to the wider box in the x direction), and k/k1 = 1 for the run with Co = 0.22.

Open with DEXTER

thumbnail Fig. 5

Dependence of ϵf on Gr Co obtained in DNS with imposed field (open symbols, red) and without (closed symbols, blue), for kf/k1 = 30. The black symbols connected by a dotted line correspond to the values of Brandenburg et al. (2012b) for kf/k1 = 5. The horizontal lines correspond to the dynamo threshold for the two values of kf/k1.

Open with DEXTER

With the TFM, we obtain the kernels αij and ηij, from which we compute We normalize α and ηt 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 urms and present . In our normalization, the molecular value is given by η/η0 = 3/ReM.

thumbnail Fig. 6

TFM coefficients versus scale separation ratio, k/kf, for Co = 0.59, ReM = 18, B0y/Beq0 = 0.05, , and ηk1/cs = 2 × 10-4.

Open with DEXTER

thumbnail Fig. 7

TFM coefficients versus Coriolis number, Co, for k/k1 = 1, ReM = 18, B0y/Beq0 = 0.05, , and ηk1/cs = 2 × 10-4.

Open with DEXTER

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

In Fig. 6 we show the dependence of the coefficients on the normalized wavenumber of the test field, k/kf. The three coefficients , , and show the same behavior of the form of (22)for , , or , while for we use (23)where , , and γ = 2.5. These results have been obtained for Co = 0.59 and B0y/Beq0 = 0.05. 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 (24)where and . Here and in the following, we keep k/kf = 1/30. For the same value of k/kf, 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 by 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 ηt that have been obtained now with the TFM are somewhat different from previous TFM studies that originally estimated ( and ). The TFM results now suggest ϵk = 0.6 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 kf/k1 that is now much larger than before (30 compared to 5 in most previous studies), which explains the reason for our choice of the subscript in ϵk.

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 nonuniform 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 ReM 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 Cα, 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 where is the advective derivative, (27)is the mean-field hydrostatic force balance, ηT = ηt + η and νT = νt + ν 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, (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, qp depends on the local field strength and is approximated by (Kemel et al. 2012a) (29)where qp0, βp, and are constants, is the normalized mean magnetic field, and is the equipartition field strength. For ReM ≲ 60, Brandenburg et al. (2012a) found β ≈ 0.33 and βp ≈ 1.05/ReM. We use as our reference model the parameters for ReM = 18, also used by Losada et al. (2013), which yields (30)In some cases we also compare with β = 0.44, 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) (31)to replace kf = urms/ 3ηt, so (32)and (Losada et al. 2013) (33)We now consider separately cases where we vary either Co or Gr. In addition, we also vary the scale separation ratio kf/k1, which is essentially a measure of the inverse turbulent diffusivity, i.e., (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 (35)As discussed above, the normalized growth rate λ/λ∗ 0 shows first a decline with increasing values of Co, but then an increase for Co > 0.13, which was argued to be a result of the dynamo effect (Losada et al. 2013). This curve has a minimum at Co ≈ 0.13. 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 α2 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 λ∗ 0. Kemel et al. (2013) proposed an empirical formula replacing λ by λ + ηtk2, where k is the wavenumber. This would lead to (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 (37)For a cubic domain, large-scale dynamo action occurs for Cα> 1, which was confirmed by Losada et al. (2013), who found the typical Beltrami fields for two supercritical cases. They used the parameters Gr = 0.033 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 β = 0.44, referred to as MFS(ii).

thumbnail Fig. 8

Nondimensional growth rate of NEMPI versus Co for MFS(i) with β = 0.33 and MFS(ii) with β = 0.44, as well as DNS for Gr = 0.033 and β0 = 0.05.

Open with DEXTER

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 λ/λ∗ 0 ≈ 1. To check this, we performed two-dimensional MFS in a squared domain of size (2π)2. The result is shown in Fig. 9 for the model parameters given in Eq. (30). When Gr is small, we find that λ/λ∗ 0 ≈ 0.3, which is below the expected value. As we increase Gr, λ/λ∗ 0 decreases until NEMPI can no longer be detected for Gr ≳ 1.2.

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 B0 (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 B0. 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 g or ηt (=νt). This suggests that our results for large values of g shown in Fig. 9 may in fact be accurate. To illustrate this more clearly, we rewrite (38)where we have defined (39)Figure 9 shows that λ/λ∗ 0 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 λ/λ∗ 0 ≈ 0.44 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 λ/urmskf versus itself. This shows that the growth rate (in units of the inverse turnover time) increases with increasing when is small. However, 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 λ/λ∗ 0 with increasing values of Gr can be approximated by the formula (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 ηtk2. In that case, contrary to our MFS, the normalized growth rate would actually increase with increasing values of Gr (see Eq. (36)).

thumbnail Fig. 9

Normalized growth rate of NEMPI versus stratification parameter Gr that varies with changing gravity, g, for Co = 0 with constant ( black filled symbols and blue open symbols), or with changing ηt = νt for constant (red open symbols). The dash-dotted line shows the approximate fit given by Eq. (40). The inset shows the growth rate normalized by the turnover time as a function of .

Open with DEXTER

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 g and varying ηt 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, 0.12, 0.21, and 1 (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 Gr = 0.12 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.

thumbnail Fig. 10

Normalized growth rate of the combined NEMPI and dynamo instability (solid lines) together with cases with pure dynamo instability (no imposed field, dashed lines) versus Co for three different values of Gr; Gr = 0.12 (blue), Gr = 0.21 (red), and Gr = 1.0 (black). In these simulations and (blue line), , (red line), and , (black line).

Open with DEXTER

In Fig. 10 we see that the dependence of λ/λ∗ 0 on Gr is opposite for small and large values of Co. When Co ≲ 0.05, an increase in Gr leads to a decrease in λ/λ∗ 0 (compare the Gr = 1 line with that for 0.21 along a cut through Co = 0.05 in Fig. 10), while for Co ≳ 0.2, an increase in Gr leads to an increase in λ/λ∗ 0 (compare all three lines in Fig. 10 along a cut through Co = 0.3). The second case is caused by the increase of the dynamo number Cα, 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 kf/k1 or Cα, 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 40 kG field strength can be produced in the solar convection zone. This is almost as strong as the 100 kG 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 10 kG fields at a depth of 10 Mm 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 Co ≳ 0.03? 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 Co ≈ 0.1. 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 Ω ≈ 3 × 10-6 s-1, a value of Co ≡ 2Ωτ = 0.1 corresponds to τ = 5 h, which is longer than the earlier MFS values of 2 h for a horizontal field (Losada et al. 2013) and 30min for a vertical field (Brandenburg et al. 2014).

Using the TFM, we have confirmed earlier findings regarding α and ηt, although for our new simulations both coefficients are somewhat larger, which is presumably due to the larger scale separation. The ratio between α and ηt 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 nonhelically 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.


Acknowledgments

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

All Figures

thumbnail Fig. 1

Visualization of and together with effective magnetic pressure for different times. Here Ω = 0.15, Co = 0.09, Gr = 0.033, and kf/k1 = 30.

Open with DEXTER
In the text
thumbnail Fig. 2

Like Fig. 1, but for a wider domain.

Open with DEXTER
In the text
thumbnail Fig. 3

Like Fig. 1, but for Ω = 0.35, so Co = 0.22.

Open with DEXTER
In the text
thumbnail Fig. 4

Comparison of the evolution of Bk/Beq0 for runs with different values of Co. In the first panel k/k1 = 1, while in the second panel k/k1 = 2 for the two runs with Co = 0.09 (label W refers to the wider box in the x direction), and k/k1 = 1 for the run with Co = 0.22.

Open with DEXTER
In the text
thumbnail Fig. 5

Dependence of ϵf on Gr Co obtained in DNS with imposed field (open symbols, red) and without (closed symbols, blue), for kf/k1 = 30. The black symbols connected by a dotted line correspond to the values of Brandenburg et al. (2012b) for kf/k1 = 5. The horizontal lines correspond to the dynamo threshold for the two values of kf/k1.

Open with DEXTER
In the text
thumbnail Fig. 6

TFM coefficients versus scale separation ratio, k/kf, for Co = 0.59, ReM = 18, B0y/Beq0 = 0.05, , and ηk1/cs = 2 × 10-4.

Open with DEXTER
In the text
thumbnail Fig. 7

TFM coefficients versus Coriolis number, Co, for k/k1 = 1, ReM = 18, B0y/Beq0 = 0.05, , and ηk1/cs = 2 × 10-4.

Open with DEXTER
In the text
thumbnail Fig. 8

Nondimensional growth rate of NEMPI versus Co for MFS(i) with β = 0.33 and MFS(ii) with β = 0.44, as well as DNS for Gr = 0.033 and β0 = 0.05.

Open with DEXTER
In the text
thumbnail Fig. 9

Normalized growth rate of NEMPI versus stratification parameter Gr that varies with changing gravity, g, for Co = 0 with constant ( black filled symbols and blue open symbols), or with changing ηt = νt for constant (red open symbols). The dash-dotted line shows the approximate fit given by Eq. (40). The inset shows the growth rate normalized by the turnover time as a function of .

Open with DEXTER
In the text
thumbnail Fig. 10

Normalized growth rate of the combined NEMPI and dynamo instability (solid lines) together with cases with pure dynamo instability (no imposed field, dashed lines) versus Co for three different values of Gr; Gr = 0.12 (blue), Gr = 0.21 (red), and Gr = 1.0 (black). In these simulations and (blue line), , (red line), and , (black line).

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.