Magnetic flux concentrations in a polytropic atmosphere
aa r X i v : . [ a s t r o - ph . S R ] F e b Astronomy&Astrophysicsmanuscript no. paper c (cid:13)
ESO 2018March 5, 2018
Magnetic flux concentrations in a polytropic atmosphere
I. R. Losada , , A. Brandenburg , , N. Kleeorin , , , and I. Rogachevskii , , 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 Department of Mechanical Engineering, Ben-Gurion University of the Negev, POB 653, Beer-Sheva 84105, Israel Department of Radio Physics, N. I. Lobachevsky State University of Nizhny Novgorod, RussiaMarch 5, 2018, Revision: 1.267
ABSTRACT
Context.
Strongly stratified hydromagnetic turbulence has recently been identified as a candidate for explaining the spontaneousformation of magnetic flux concentrations by the negative effective magnetic pressure instability (NEMPI). Much of this work hasbeen done for isothermal layers, in which the density scale height is constant throughout.
Aims.
We now want to know whether earlier conclusions regarding the size of magnetic structures and their growth rates carry overto the case of polytropic layers, in which the scale height decreases sharply as one approaches the surface.
Methods.
To allow for a continuous transition from isothermal to polytropic layers, we employ a generalization of the exponentialfunction known as the q -exponential. This implies that the top of the polytropic layer shifts with changing polytropic index suchthat the scale height is always the same at some reference height. We used both mean-field simulations (MFS) and direct numericalsimulations (DNS) of forced stratified turbulence to determine the resulting flux concentrations in polytropic layers. Cases of bothhorizontal and vertical applied magnetic fields were considered. Results.
Magnetic structures begin to form at a depth where the magnetic field strength is a small fraction of the local equipartitionfield strength with respect to the turbulent kinetic energy. Unlike the isothermal case where stronger fields can give rise to magneticflux concentrations at larger depths, in the polytropic case the growth rate of NEMPI decreases for structures deeper down. Moreover,the structures that form higher up have a smaller horizontal scale of about four times their local depth. For vertical fields, magneticstructures of super-equipartition strengths are formed, because such fields survive downward advection that causes NEMPI with hori-zontal magnetic fields to reach premature nonlinear saturation by what is called the “potato-sack” effect. The horizontal cross-sectionof such structures found in DNS is approximately circular, which is reproduced with MFS of NEMPI using a vertical magnetic field.
Conclusions.
Results based on isothermal models can be applied locally to polytropic layers. For vertical fields, magnetic flux con-centrations of super-equipartition strengths form, which supports suggestions that sunspot formation might be a shallow phenomenon.
Key words. magnetohydrodynamics (MHD) – hydrodynamics – turbulence – Sun: dynamo
1. Introduction
In a turbulent medium, the turbulent pressure can lead to dy-namically important effects. In particular, a stratified layer canattain a density distribution that is significantly altered com-pared to the nonturbulent case. In addition, magnetic fields canchange the situation further, because it can locally suppressthe turbulence and thus reduce the total turbulent pressure (thesum of hydrodynamic and magnetic turbulent contributions). Onlength scales encompassing many turbulent eddies, this totalturbulent pressure reduction must be compensated for by addi-tional gas pressure, which can lead to a density enhancementand thus to horizontal magnetic structures that become heavierthan the surroundings and sink (Brandenburg et al., 2011). Thisis quite the contrary of magnetic buoyancy, which is still ex-pected to operate on the smaller scale of magnetic flux tubesand in the absence of turbulence. Both effects can lead to in-stability: the latter is the magnetic buoyancy or interchangeinstability (Newcomb, 1961; Parker, 1966), and the formeris now often referred to as negative effective magnetic pres-sure instability (NEMPI), which has been studied at the levelof mean-field theory for the past two decades (Kleeorin et al.,1989, 1990, 1993, 1996; Kleeorin & Rogachevskii, 1994;Rogachevskii & Kleeorin, 2007). These are instabilities of astratified continuous magnetic field, while the usual magnetic buoyancy instability requires nonuniform and initially sepa-rated horizontal magnetic flux tubes (Parker, 1955; Spruit, 1981;Sch¨ussler et al., 1994).Unlike the magnetic buoyancy instability, NEMPI occurs atthe expense of turbulent energy rather than the energy of thegravitational field. The latter is the energy source of the mag-netic buoyancy or interchange instability. NEMPI is caused bya negative turbulent contribution to the effective mean magneticpressure (the sum of nonturbulent and turbulent contributions).For large Reynolds numbers, this turbulent contribution to theeffective magnetic pressure is larger than the nonturbulent one.This results in the excitation of NEMPI and the formation oflarge-scale magnetic structures – even from an originally uni-form mean magnetic field.Direct numerical simulations (DNS) have recently demon-strated the operation of NEMPI in isothermally stratified layers(Brandenburg et al., 2011; Kemel et al., 2012b). This is a partic-ularly simple case in that the density scale height is constant; i.e.,the computational burden of covering large density variation isdistributed over the depth of the entire layer. In spite of this sim-plification, it has been argued that NEMPI is important for ex-plaining prominent features in the manifestation of solar surfaceactivity. In particular, it has been associated with the formationof active regions (Kemel et al., 2013; Warnecke et al., 2013) andsunspots (Brandenburg et al., 2013, 2014). However, it is now
1. R. Losada et al.: Magnetic flux concentrations in a polytropic atmosphere important to examine the validity of conclusions based on suchsimplifications using more realistic models. In this paper, wetherefore now consider a polytropic stratification, for which thedensity scale height is smallest in the upper layers, and the den-sity variation therefore greatest.NEMPI is a large-scale instability that can be excited in strat-ified small-scale turbulence. This requires (i) sufficient scale sep-aration in the sense that the maximum scale of turbulent motions, ℓ , must be much smaller than the scale of the system, L , and (ii)strong density stratification such that the density scale height H ρ is much smaller than L ; i.e., ℓ ≪ H ρ ≪ L. (1)However, both the size of turbulent motions and the typical sizeof perturbations due to NEMPI can be related to the densityscale height. Furthermore, earlier work of Kemel et al. (2013)using isothermal layers shows that the scale of perturbations dueto NEMPI exceeds the typical density scale height. Unlike theisothermal case, in which the scale height is constant, it de-creases rapidly with height in a polytropic layer. It is then un-clear how such structures could fit into the narrow space left bythe stratification and whether the scalings derived for the isother-mal case can still be applied locally to polytropic layers.NEMPI has already been studied previously for polytropiclayers in mean-field simulations (MFS) (Brandenburg et al.,2010; K¨apyl¨a et al., 2012; Jabbari et al., 2013), but no system-atic comparison has been made with NEMPI in isothermal or inpolytropic layers with different values of the polytropic index.This will be done in the present paper, both in MFS and DNS.Those two complementary types of simulations have provedto be a good tool for understanding the underlying physics ofNEMPI. An example are the studies of the effects of rotationon NEMPI (Losada et al., 2012, 2013), where MFS have beenable to give quantitatively useful predictions before correspond-ing DNS were able to confirm the resulting dependence.
2. Polytropic stratification
We discuss here the equation for the vertical profile of the fluiddensity in a polytropic layer. In a Cartesian plane-parallel layerwith polytropic stratification, the temperature gradient is con-stant, so the temperature goes linearly to zero at z ∞ . The tem-perature, T , is proportional to the square of the sound speed, c ,and thus also to the density scale height H ρ ( z ) , which is givenby H ρ = c /g for an isentropic stratification, where g is the ac-celeration due to the gravity. For a perfect gas, the density ρ isproportional to T n , and the pressure p is proportional to T n +1 ,such that p/ρ is proportional to T , where n is the polytropic in-dex. Furthermore, we have p ( z ) ∝ ρ ( z ) Γ , where Γ = ( n + 1) /n is another useful coefficient.For a perfect gas, the specific entropy can be defined (up toan additive constant) as s = c v ln( p/ρ γ ) , where γ = c p /c v is theratio of specific heats at constant pressure and constant density,respectively. For a polytropic stratification, we have exp( s/c v ) = p/ρ γ ∝ ρ Γ − γ , (2)so s is constant when Γ = γ , which is the case for an isen-tropic stratification. In the following, we make this assumptionand specify from now only the value of γ . For a monatomic gas,we have γ = 5 / , which is relevant for the Sun, while for a di-atomic molecular gas, we have γ = 7 / , which is relevant forair. In those cases, a stratification with Γ = γ can be motivatedby assuming perfect mixing across the layer. The isothermal case Fig. 1.
Isothermal and polytropic relations for different values of γ when calculated using the conventional formula ρ ∝ ( z ∞ − z ) n with n = 1 / ( γ − . Fig. 2.
Polytropes [Eq. (4)] with γ = 1 (solid line), 1.2 (dash-dotted), 1.4 (dotted), and 5/3 (dashed) and density scale height[Eq. (5)] for − ≤ − z/H ρ ≤ . . The total density contrast issimilar for γ = 1 and 5/3.with γ = 1 can be motivated by assuming rapid heating/coolingto a constant temperature.Our aim is to study the change in the properties of NEMPIin a continuous fashion as we go from an isothermally strati-fied layer to a polytropic one. In the latter case, the fluid densityvaries in a power law fashion, ρ ∝ ( z ∞ − z ) n , while in the for-mer, it varies exponentially, ρ ∝ exp( z ∞ − z ) . This is shown inFig. 1 where we compare the exponential isothermal atmospherewith a family of polytropic atmospheres with γ = 1 . , 1.4,and 5/3. Clearly, there is no continuous connection between the
2. R. Losada et al.: Magnetic flux concentrations in a polytropic atmosphere
Fig. 3.
Sketch showing the expected size distribution of thenearly circular NEMPI eigenfunction structures at differentheights.isothermal case and the polytropic one in the limit γ → . Thiscannot be fixed by rescaling the isothermal density stratification,because in Fig. 1 its values would still lie closer to 5/3 than to1.4 or 1.2. Another problem with this description is that for poly-tropic solutions the density is always zero at z = z ∞ , but finite inthe isothermal case. These different behaviors between isother-mal and polytropic atmospheres can be unified by using the gen-eralized exponential function known as the “ q -exponential” (see,e.g., Yamano, 2002), which is defined as e q ( x ) = [1 + (1 − q ) x ] / (1 − q ) , (3)where the parameter q is related to γ via q = 2 − γ . Thisgeneralization of the usual exponential function was originallyintroduced by Tsallis (1988) in connection with a possiblegeneralization of the Boltzmann-Gibbs statistics. Its usefulnessin connection with stellar polytropes has been employed byPlastino & Plastino (1993). Thus, the density stratification isgiven by ρρ = (cid:20) γ − (cid:18) − zH ρ (cid:19)(cid:21) / ( γ − = (cid:18) − znH ρ (cid:19) n , (4)which reduces to ρ/ρ = exp( − z/H ρ ) for isothermal stratifi-cation with γ → and n → ∞ . The density scale height is thengiven by H ρ ( z ) = H ρ − ( γ − z = H ρ − z/n. (5)In the following, we measure lengths in units of H ρ = H ρ (0) .In Fig. 2 we show the dependencies of ρ ( z ) and H ρ given byEqs. (4) and (5) for different values of γ . Compared with Fig. 1,where z ∞ is held fixed, in Fig. 2 it is equal to z ∞ = nH ρ = H ρ / ( γ − . The total density contrast is roughly the same inall four cases for different γ , but for increasing values of γ , thevertical density gradient becomes progressively stronger in theupper layers.Assuming that the radius R of the resulting structures is pro-portional to H ρ , we sketch in Fig. 3 a situation in which R ishalf the depth. With ρ ∝ ( z ∞ − z ) n , the density scale height is H ρ ( z ) = ( z ∞ − z ) /n . Thus, Fig. 3 applies to a case in which R = ( z ∞ − z ) / n/ H ρ ( z ) . The solar convection zoneis nearly isentropic and well described by n = 3 / . This meansthat the structures of Fig. 3 have R = (3 / H ρ ( z ) . The re-sults of Kemel et al. (2013) and Brandenburg et al. (2014) sug-gest that the horizontal wavenumber of structures, k ⊥ , formed by NEMPI, is less than or about H − ρ , so their horizontal wave-length is < ∼ πH ρ . One wavelength corresponds to the dis-tance between two nodes, i.e., the distance between two spheres,which is R . Thus, in the isothermal case we have R/H ρ =2 π/ ≈ . , which implies that such a structure would not fitinto the isentropic atmosphere described above. This providesan additional motivation for our present work.
3. DNS study
In this section we study NEMPI in DNS for the polytropic layer.Corresponding MFS are presented in Sect. 4.
We solve the equations for the magnetic vector potential A , thevelocity U , and the density ρ , in the form ∂ A ∂t = U × B − ηµ J , (6) D U D t = 1 ρ J × B + f − ν Q − ∇ H, (7) D ρ D t = − ρ ∇ · U , (8)where D / D t = ∂/∂t + U · ∇ is the advective derivative withrespect to the actual (turbulent) flow, B = B + ∇ × A is themagnetic field, B the imposed uniform field, J = ∇ × B /µ the current density, µ the vacuum permeability, H = h + Φ thereduced enthalpy, h = c p T the enthalpy, Φ is the gravitationalpotential, − Q = ∇ U + ∇∇ · U + 2 S ∇ ln ρ (9)is a term appearing in the viscous force − ν Q , S is the tracelessrate-of-strain tensor with components S ij = ( ∇ j U i + ∇ i U j ) − δ ij ∇ · U , (10) ν is the kinematic viscosity, and η is the magnetic diffusioncoefficient caused by electrical conductivity of the fluid. As inLosada et al. (2012), z corresponds to radius, x to colatitude,and y to azimuth. The forcing function f consists of random,white-in-time, plane, nonpolarized waves with a certain averagewavenumber k f . In the DNS we use stress-free boundary conditions for the ve-locity at the top and bottom; i.e., ∇ z U x = ∇ z U y = U z = 0 .For the magnetic field we use either perfect conductor boundaryconditions, A x = A y = ∇ z A z = 0 , or vertical field conditions, ∇ z A x = ∇ z A y = A z = 0 , again at both the top and bottom.All variables are assumed periodic in the x and y directions.The turbulent rms velocity is approximately independent of z with u rms = h u i / ≈ . c s . The gravitational accelera-tion g = (0 , , − g ) is chosen such that k H ρ = 1 , where k = 2 π/L and L is the size of the domain. With one ex-ception (Sect. 3.5), we always use the value k f /k = 30 forthe scale separation ratio. For B we choose either a horizon-tal field pointing in the y direction or a vertical one pointing inthe z direction. The latter case, B = (0 , , B ) , is usually com-bined with the use of the vertical field boundary condition, whilethe former one, B = (0 , B , , is combined with using per-fect conductor boundary conditions. The strength of the imposed
3. R. Losada et al.: Magnetic flux concentrations in a polytropic atmosphere
Fig. 4.
Snapshots of B y from DNS for γ = 5 / and β = 0 . (upper row), . (middle row), and . (lower row) at differenttimes (indicated in turbulent-diffusive times, increasing from left to right) in the presence of a horizontal field using the perfectconductor boundary condition.field is expressed in terms of B eq0 = B eq ( z = 0) , which is theequipartition field strength at z = 0 . Here, the equipartition field B eq ( z ) = ( µ ρ ( z )) / u rms . The imposed field is normalizedby B eq0 and denoted as β = B /B eq0 , while β = | B | /B eq isthe modulus of the normalized mean magnetic field. Time is ex-pressed in terms of the turbulent-diffusive time, τ td = H ρ /η t0 ,where η t0 = u rms / k f (Sur et al., 2008) is an estimate for theturbulent magnetic diffusivity used in the DNS.Our values of ν and η are characterized by specifying thekinetic and magnetic Reynolds numbers,Re = u rms /νk f , R m = u rms /ηk f . (11)In most of this paper (except in Sect. 3.5) we use Re = 36 and R m = 18 , which are also the values used by Kemel et al. (2013).The DNS are performed with the P ENCIL C ODE , http://pencil-code.googlecode.com , which usessixth-order explicit finite differences in space and a third-orderaccurate time stepping method. We use a numerical resolutionof mesh points. We focus on the case γ = 5 / and show in Fig. 4 visualizationsof B y at different instants for three values of the imposed hori-zontal magnetic field strength. For β = 0 . , a magnetic struc-ture is clearly visible at t/τ td = 1 . , while for β = 0 . struc-tures are already fully developed at t/τ td = 0 . . In that case( β = 0 . ), at early times ( t/τ td = 0 . ), there are two struc-tures, which then begin to merge at t/τ td = 1 . . The growthrate of the magnetic structure is found to be λ ≈ η t0 /H ρ forthe runs shown in Fig. 4. This is less than the value of λ ≈ η t0 /H ρ found earlier for the isothermal case (Kemel et al.,2012b). For γ = 5 / and β ≤ . , the magnetic structures be-come smaller ( k ⊥ H ρ = 2 ) near the surface. In the nonlin-ear regime, i.e., at late times, the structures move downwarddue to the so-called “potato-sack” effect, which was first seenin MFS (Brandenburg et al., 2010) and later confirmed in DNS(Brandenburg et al., 2011). The magnetic structures sink in thenonlinear stage of NEMPI, because an increase in the mean mag-netic field inside the horizontal magnetic flux tube increasesthe absolute value of the effective magnetic pressure. On theother hand, a decrease in the negative effective magnetic pres-sure is balanced out by increased gas pressure, which in turnleads to higher density, so the magnetic structures become heav-ier than the surroundings and sink. This potato-sack effect hasbeen clearly observed in the present DNS with the polytropiclayer (see the right column in Fig. 4). Recent DNS using isothermal layers have shown that strong cir-cular flux concentrations can be produced in the case of a verti-cal magnetic field (Brandenburg et al., 2013, 2014). This is alsoobserved in the present study of a polytropic layer; see Fig. 5,where we show the evolution of B z on the periphery of the com-putational domain for γ = 5 / and β = 0 . at different times.A difference to the DNS for γ = 1 (Brandenburg et al., 2013)seems to be that for γ = 5 / the magnetic structures are shal-lower than for γ = 1 ; see Fig. 6, where we show xy and xz slices of B z through the spot. Owing to periodicity in the xy plane, we have shifted the location of the spot to x = y = 0 .We note also that the field lines of the averaged magnetic fieldshow a structure rather similar to the one found in MFS ofBrandenburg et al. (2014). The origin of circular structures is as-sociated with a cylindrical symmetry for the vertical magneticfield. The growth rate of the magnetic field in the spot is found
4. R. Losada et al.: Magnetic flux concentrations in a polytropic atmosphere
Fig. 5.
Snapshots from DNS showing B z on the periphery of the computational domain for γ = 5 / and β = 0 . at differenttimes for the case of a vertical field using the vertical field boundary condition. Fig. 6.
Cuts of B z /B eq ( z ) in (a) the xy plane at the top boundary ( z/H ρ = 1 . ) and (b) the xz plane through the middle of thespot at y = 0 for γ = 5 / and β = 0 . . In the xz cut, we also show magnetic field lines and flow vectors obtained by numericallyaveraging in azimuth around the spot axis.to be λ ≈ . η t0 /H ρ , which is similar to the value of 1.3 foundearlier for the isothermal case (Brandenburg et al., 2013). As pointed out in Sect. 1, the main reason for the formationof strongly inhomogeneous large-scale magnetic structures isthe negative contribution of turbulence to the large-scale mag-netic pressure, so that the effective magnetic pressure (the sumof turbulent and nonturbulent contributions) can be negative athigh magnetic Reynolds numbers. The effective magnetic pres-sure has been determined from DNS for isothermally stratified forced turbulence (Brandenburg et al., 2010, 2012) and for tur-bulent convection (K¨apyl¨a et al., 2012). To see whether the na-ture of polytropic stratification has any influence on the effectivemagnetic pressure, we use DNS.We first explain the essence of the effect of turbulence on theeffective magnetic pressure. We consider the momentum equa-tion in the form ∂∂t ρ U i = − ∂∂x j Π ij + ρ g i , (12)where Π ij = ρ U i U j + δ ij (cid:0) p + B / µ (cid:1) − B i B j /µ − νρ S ij (13)
5. R. Losada et al.: Magnetic flux concentrations in a polytropic atmosphere is the momentum stress tensor and δ ij the Kronecker tensor.Ignoring correlations between velocity and density fluctua-tions for low-Mach number turbulence, the averaged momentumequation is ∂∂t ρ U i = − ∂∂x j Π ij + ρ g i , (14)where an overbar means xy averaging, Π ij = Π m ij + Π f ij is themean momentum stress tensor, split into contributions resultingfrom the mean field (indicated by superscript m) and the fluc-tuating field (indicated by superscript f). The tensor Π m ij has thesame form as Eq. (13), but all quantities have now attained anoverbar; i.e., Π m ij = ρ U i U j + δ ij (cid:0) p + B / µ (cid:1) − B i B j /µ − νρ S ij . (15)The contributions, Π f ij , which result from the fluctuations u = U − U and b = B − B of velocity and magnetic fields, respec-tively, are determined by the sum of the Reynolds and Maxwellstress tensors: Π f ij = ρ u i u j + δ ij b / µ − b i b j /µ . (16)This contribution, together with the contribution from the meanfield, Π m ij , comprises the total mean momentum tensor. The con-tribution from the fluctuating fields is split into a contributionthat is independent of the mean magnetic field Π f , ij (which de-termines the turbulent viscosity and background turbulent pres-sure) and a contribution that depends on the mean magnetic field Π f ,Bij . The difference between the two, ∆Π f ij = Π f ,Bij − Π f , ij , iscaused by the mean magnetic field and is parameterized in theform ∆Π f ij = µ − (cid:0) q s B i B j − q p δ ij B / − q g ˆ g i ˆ g j B (cid:1) , (17)where the coefficient q p represents the isotropic turbulence con-tribution to the mean magnetic pressure, the coefficient q s rep-resents the turbulence contribution to the mean magnetic ten-sion, while the coefficient q g is the anisotropic turbulence con-tribution to the mean magnetic pressure, and it characterizesthe effect of vertical variations of the magnetic field causedby the vertical magnetic pressure gradient. Here, ˆ g i is the unitvector in the direction of the gravity field (in the vertical di-rection). We consider cases with horizontal and vertical meanmagnetic fields separately. Analytically, the coefficients q p , q g , and q s have been obtained using both the spectral τ approach(Rogachevskii & Kleeorin, 2007) and the renormalization ap-proach (Kleeorin & Rogachevskii, 1994). The form of Eq. (17)is also obtained using simple symmetry arguments; e.g., for ahorizontal field, the linear combination of three independent truetensors, δ ij , ˆ g i ˆ g j and B i B j , yields Eq. (17), while for the verti-cal field, the linear combination of two independent true tensors, δ ij and B i B j , yields this ansatz.Previous DNS studies (Brandenburg et al., 2012) haveshown that q s and q g are negligible for forced turbulence. Toavoid the formation of magnetic structures in the nonlinear stageof NEMPI, which would modify our results, we use here a lowerscale separation ratio, k f /k = 5 , keeping k H ρ = 1 , and us-ing Re = 140 and R m = 70 , as in Brandenburg et al. (2012).To determine q p ( β ) , it is sufficient to measure the three diagonalcomponents of Π f ij both with and without an imposed magneticfield. Fig. 7.
Effective magnetic pressure obtained from DNS in apolytropic layer with different γ for horizontal (H, red curves)and vertical (V, blue curves) mean magnetic fields.In Fig. 7 we present the results for forced turbulence in thepolytropic layer with different γ for horizontal and vertical meanmagnetic fields. It turns out that the normalized effective mag-netic pressure, P eff = (1 − q p ) β , (18)has a minimum value P min at β min . Following Kemel et al.(2012a), the function q p ( β ) is approximated by: q p ( β ) = q p0 β /β = β ⋆ β + β , (19)where q p0 , β p , and β ⋆ = β p q / are constants. This equation canbe understood as a quenching formula for the effective magneticpressure; see Jabbari et al. (2013). The coefficients β p and β ⋆ arerelated to P min at β min via (Kemel et al., 2012a) β p = β .p − P min , β ⋆ = β p + p − P min . (20)In Fig. 8 we show these fitting parameters for the function q p ( β ) for polytropic layer with different γ for horizontal and verticalmean magnetic fields. The effects of negative effective magneticpressure are generally stronger for vertical magnetic fields ( q p0 is larger and β p smaller) than for horizontal ones ( q p0 is smallerand β p larger), but the values β ⋆ = β p q / are similar in bothcases, and increasing from 0.35 (for γ = 1 ) to 0.6 (for γ = 5 / );see Fig. 8.
4. Mean-field study
We now consider two sets of parameters that we refer to asModel I (with q p0 = 32 and β p = 0 . corresponding to β ⋆ = 0 . ) and Model II ( q p0 = 9 and β p = 0 . correspondingto β ⋆ = 0 . ). These cases are representative of the strong (large β ⋆ ) and weak (small β ⋆ ) effects of NEMPI. Following earlierstudies (Brandenburg et al., 2012), we find q s to be compatiblewith zero. We thus neglect this coefficient in the following. The purpose of this section is to summarize the findings for theisothermal case in MFS. One of the key results is the prediction
6. R. Losada et al.: Magnetic flux concentrations in a polytropic atmosphere
Fig. 8.
Parameters q p0 , β p , and β ⋆ for the function q p ( β ) [seeEq. (19)] versus γ for horizontal (red line) and vertical (blueline) mean magnetic fields.of the growth rate of NEMPI. The work of Kemel et al. (2013)showed that in the ideal case (no turbulent diffusion), the growthrate λ is approximated well by λ ≈ β ⋆ u rms /H ρ (no turbulent diffusion) . (21)However, turbulent magnetic diffusion, η t , can clearly not beneglected and is chiefly responsible for shutting off NEMPIif the turbulent eddies are too big and η t too large. This wasdemonstrated in Fig. 17 of Brandenburg et al. (2012). A heuris-tic ansatz, which is motivated by similar circumstances in mean-field dynamo theory (Krause & R¨adler, 1980), is to add a term − η t k to the righthand side of Eq. (21), where k is the wavenum-ber of NEMPI.To specify the expression for λ , we normalize the wavenum-ber of the perturbations by the inverse density scale height anddenote this by κ ≡ kH ρ . The wavenumber of the energy-carrying turbulent eddies k f is in nondimensional form κ f ≡ k f H ρ , and the normalized horizontal wavenumber of the re-sulting magnetic structures is referred to as κ ⊥ = k ⊥ H ρ .For NEMPI, these values have been estimated to be κ ⊥ =0 . –1.0, and can even be smaller for vertical magnetic fields(Brandenburg et al., 2014). Using an approximate aspect ratioof unity for magnetic structures, we have κ = √ κ ⊥ ≈ . –1.4. In stellar mixing length theory (Vitense, 1953), the mixinglength is ℓ f = α mix H ρ , where α mix ≈ . is a nondimensionalmixing length parameter. Since k f = 2 π/ℓ f , we arrive at thefollowing estimate: κ f = 2 πγ/α mix ≈ . . [Owing to a confu-sion between pressure and density scale heights, this value wasunderestimated by Kemel et al. (2013) to be . , although an in- dependent calculation of this value from turbulent convectionsimulations would still be useful.] Using η t ≈ u rms / k f , theturbulent magnetic diffusive rate for an isothermal atmosphereis given by η t k = u rms H ρ κ κ f , (22)and the growth rate of NEMPI in that normalization is λη t k = 3 β ⋆ κ f κ − . (23)Using β ⋆ = 0 . , which is the relevant value for high mag-netic Reynolds numbers (Brandenburg et al., 2012), we find λ/η t k ≈ . – . . However, since Fig. 8 shows an increase of β ⋆ with increasing polytropic index, one might expect a corre-sponding increase in the growth rate of NEMPI for a polytropiclayer, in which H ρ varies strongly with height z . Indeed, in apolytropic atmosphere, H ρ is proportional to depth. Thus, at anygiven depth there is a layer beneath, where the stratification isless strong and the growth rate of NEMPI is lower. In addition,there is a thinner, more strongly stratified layer above, whereNEMPI might grow faster if only the structures generated byNEMPI have enough room to develop before they touch the topof the atmosphere at z ∞ , where the temperature vanishes. In the following, we consider MFS and compare it withDNS. We also compare our MFS results with those of DNS(Brandenburg et al., 2011; Kemel et al., 2013) using a similarpolytropic setup.The governing equations for the mean quantities (denoted byan overbar) are fairly similar to those for the original equations,except that in the MFS viscosity and magnetic diffusivity arereplaced by their turbulent counterparts, and the mean Lorentzforce is supplemented by a parameterization of the turbulent con-tribution to the effective magnetic pressure.The evolution equations for mean vector potential A , meanvelocity U , and mean density ρ , are ∂ A ∂t = U × B − η T µ J , (24) D U D t = 1 ρ (cid:2) J × B + ∇ ( q p B / µ ) (cid:3) − ν T Q − ∇ H, (25) D ρ D t = − ρ ∇ · U , (26)where D / D t = ∂/∂t + U · ∇ is the advective derivative with re-spect to the mean flow, ρ the mean density, H = h + Φ the meanreduced enthalpy, h = c p T the mean enthalpy, T ∝ ρ γ − themean temperature, Φ the gravitational potential, η T = η t + η ,and ν T = ν t + ν are the sums of turbulent and microphysicalvalues of magnetic diffusivity and kinematic viscosities, respec-tively. Also, J = ∇ × B /µ is the mean current density, µ isthe vacuum permeability, − Q = ∇ U + ∇∇ · U + 2 S ∇ ln ρ (27)is a term appearing in the mean viscous force − ν T Q , where S isthe traceless rate-of-strain tensor of the mean flow with compo-nents S ij = ( ∇ j U i + ∇ i U j ) − δ ij ∇ · U , and finally the term ∇ ( q p B / µ ) on the righthand side of Eq. (25) determines the
7. R. Losada et al.: Magnetic flux concentrations in a polytropic atmosphere
Fig. 9.
Comparison of P eff ( z ) and λ ( z ) in MFS for polytropiclayers with four values of γ (top to bottom) and three values of β (different line types).turbulent contribution to the effective magnetic pressure. Here, q p depends on the local field strength; see Eq. (19). This termenters with a plus sign, so positive values of q p correspond toa suppression of the total turbulent pressure. The net effect ofthe mean magnetic field leads to an effective mean magneticpressure that becomes negative for q p > . This can indeedbe the case for magnetic Reynolds numbers well above unity(Brandenburg et al., 2012); see also Fig. 7 for a polytropic atmo-sphere.The boundary conditions for MFS are the same as for DNS,i.e., stress-free for the mean velocity at the top and bottom. Forthe mean magnetic field, we use either perfect conductor bound-ary conditions (for horizontal, imposed magnetic fields) or ver-tical field conditions (for vertical, imposed fields) at the top andbottom. All mean-field variables are assumed to be periodic inthe x and y directions. The MFS are performed again with theP ENCIL C ODE , which is equipped with a mean-field module forsolving the corresponding equations.
To get an idea about the vertical dependence of NEMPI, wenow consider the resulting dependencies of P eff ( z ) ; see the left-hand panels of Fig. 9. We note that P eff is just a function of β (Eqs. (18) and (19)), which allows us to approximate the local growth rates as (Rogachevskii & Kleeorin, 2007; Kemel et al.,2013) λ = v A H ρ (cid:18) − P eff d β (cid:19) , (28)which are plotted in the righthand panels of Fig. 9 for Model I. To analyze the kinematic stage of MFS, we measure the value ofthe maximum downflow speed, | U | downmax at each height. We thendetermine the time interval during which the maximum down-flow speed increases exponentially and when the height of thepeak is constant and equal to z B . This yields the growth rate ofthe instability as λ = d ln | U | downmax / d t .In Figs. 10 and 11 we plot, respectively for Models I and II, λ (in units of η t /H ρ ) and z B versus horizontal imposed magneticfield strength, B /B eq . The maximum growth rates for γ = 5 / and γ = 1 are similar for both models (4–5 η t /H ρ for Model Iand 16–20 η t /H ρ for Model II). It turns out that for γ = 5 / ,the growth rate λ attains a maximum at some value B = B max ,and then it decreases with increasing B , while in an isothermalrun λ is nearly constant for greater field strength, except nearthe surface where the proximity to the boundary is too small.This close proximity reduces the growth rate. For Model I with γ = 1 , the decline of λ (toward weaker fields on the lefthandside of the plot) begins when the distance to the top boundary( z top − z B ≈ . H ρ for β = 0 . ) is less than the radius ofmagnetic structures ( R ≈ π/k ⊥ ≈ . H ρ using k ⊥ H ρ = 1 ).In Model II with γ = 1 , the decline of λ occurs for strongerfields, but the distance to the top boundary ( ≈ . H ρ ) is stillnearly the same as for Model I.In an isothermal layer, the height where the eigenfunctionpeaks is known to decrease with increasing field strength; seeFigure 6 of Kemel et al. (2012a). One might have expected thisdecrease to be less steep in the polytropic case, because the op-timal depth where NEMPI occurs cannot easily be decreasedwithout suffering a dramatic decrease of the growth rate. Thisis however not the case, and we find that the optimal depth ofNEMPI is now falling off more quickly in Model I, but is moresimilar for Model II; see the second panels of Figs. 10 and 11.This means that in a polytropic layer, NEMPI works more ef-fectively, and its growth rate is fastest when the magnetic field isnot too strong. At the same time, the optimal depth of NEMPI in-creases, i.e., the resulting value of z B increases as B decreases.The resulting growth rates are somewhat less for Model I andsomewhat higher for Model II than those of earlier mean-fieldcalculations of Kemel et al. (2012a) (their Fig. 6), who found λH /η t ≈ . in a model with β ⋆ = 0 . and β p = 0 . .These differences in the growth rates are plausibly explained bydifferences in the mean-field parameters.Visualizations of the resulting horizontal field structures areshown in Fig. 12 for two values of γ and B . Increase in theparameter γ results in a stronger localization of the instability atthe surface layer, where the density scale height is minimum andthe growth of NEMPI is strongest. In the presence of a vertical field, the early evolution of the in-stability is similar to that for a horizontal field. In both cases,the maximum field strength occurs at a somewhat larger depth
8. R. Losada et al.: Magnetic flux concentrations in a polytropic atmosphere
Fig. 10.
Dependence of growth rate and height where theeigenfunction attains its maximum value (the optimal depth ofNEMPI) on field strength from MFS for different values of γ inthe presence of a horizontal field for Model I. Fig. 11.
Same as Fig. 10 (horizontal field), but for Model II.when saturation is reached, except that shortly before satura-tion there is a brief interval during which the location of max-imum field strength rises slightly upwards in the vertical fieldcase. In the saturated case, however, the flux concentrations fromNEMPI are much stronger compared to the case of a horizon-tal field and it leads to the formation of magnetic flux concen-
Fig. 12.
Snapshots of B y from MFS during the kinematicgrowth phase for different values of γ and B /B eq0 for γ = 1 (top) and / (bottom) with β = 0 . (left) and 0.1 (right) inthe presence of a horizontal field for Model I.trations of equipartition field strength (Brandenburg et al., 2013,2014). This is possible because the resulting vertical flux tube isnot advected downward with the flow that develops as a con-sequence of NEMPI. The latter effect is the aforementioned“potato-sack” effect, which acts as a nonlinear saturation mech-anism of NEMPI with a horizontal field.In Figs. 13 and 14 we plot the growth rates λ and the heightswhere the eigenfunction attains its maximum values for different β = B /B eq0 for Models I and II, respectively. For γ = 5 / ,the maximum growth rate is higher larger than for γ = 1 . Thisis true for Models I and II, where they are 8–10 η t /H ρ for γ = 5 / and 5–7 η t /H ρ for γ = 1 . The nonmonotonous be-havior seen in the dependence of λ on B is suggestive of thepresence of different mode structures, although a direct inspec-tion of the resulting magnetic field did not show any obvious dif-ferences. However, this irregular behavior may be related to ar-tifacts resulting from a finite domain size and were not regardedimportant enough to justify further investigation.Next we focus on a comparison of the growth rates ob-tained from MFS for horizontal and vertical fields. The valuesof β , z B , and β ( z B ) = B /B eq ( z B ) for horizontal and ver-tical fields are compared in Table 1 for both models. We seethat NEMPI is most effective in regions where the mean mag-netic field is a small fraction of the local equipartition field andtypically slightly less for γ = 5 / than for γ = 1 . Indeed,for Model I, β ( z B ) is 3–4% for horizontal fields and 3–6% forvertical fields, while for Model II, β ( z B ) is 12–17% for hori-zontal fields and 8–22% for vertical fields. Here, we have used β ( z B ) = β e − γ ( z B / H ρ ) , where e q ( x ) is the q -exponentialfunction defined by Eq. (3).We expect that higher values of β ⋆ will lead to greater growthrates. To verify this, we compare in Fig. 15 the time evolutionsof U rms for Models I (with β ⋆ = 0 . ) and Model II (with β ⋆ = 0 . ). The growth rate has now increased by a factor of2.4 (from λH ρ /η t = 3 . to 9.5), which is slightly more thanwhat is expected from β ⋆ , which has increased by a factor of1.9. This change in the growth rate can also be seen in Fig. 11
9. R. Losada et al.: Magnetic flux concentrations in a polytropic atmosphere
Fig. 13.
Dependence of growth rate and optimal depth ofNEMPI on field strength from MFS for different values of γ inthe presence of a vertical field for Model I. Fig. 14.
Same as Fig. 13 (vertical field), but for Model II.and Fig. 14 for Model II (in comparison with Figs. 10 and 13 forModel I). The dependence of the growth rate on the magneticfield strength is qualitatively similar for Models I and II. In par-ticular, it becomes constant for γ = 1 , but declines for γ = 5 / as the field increases. The increase of z B with B is, however,less strong for Model II. Fig. 15.
Comparison of U rms from MFS for Models I (solidline) and II (dashed line), showing exponential growth followedby nonlinear saturation. In both cases we have γ = 5 / and animposed vertical magnetic field with β = 0 . .Snapshots of B z from MFS for γ = 5 / and β = 0 . atdifferent times for Model I are shown in Fig. 16. Comparisonwith the results of MFS for Model II (see Fig. 17) shows thatModel II fits the DNS better. This is also seen by comparingFig. 17 with Fig. 5. However, our basic conclusions formulatedin this paper are not affected.
5. Conclusions
The present work has demonstrated that in a polytropic layer,both in MFS and DNS, NEMPI develops primarily in the upper-most layers, provided the mean magnetic field is not too strong.If the field gets stronger, NEMPI can still develop, but the mag-netic structures now occur at greater depths and the growth rateof NEMPI is lower. However, at some point when the mag-netic field gets too strong, NEMPI is suppressed in the case ofa polytropic layer, while it would still operate in the isothermalcase, provided the domain is deep enough. The slow down ofNEMPI is not directly a consequence of a longer turnover timeat greater depths, but it is related to stratification being too weakfor NEMPI to be excited.By and large, the scaling relations determined previously forisothermal layers with constant scale height still seem to applylocally to polytropic layers with variable scale heights. In partic-ular, the horizontal scale of structures was previously determinedto be about – H ρ (Kemel et al., 2013; Brandenburg et al.,2014). Looking now at Fig. 4, we see that for β = 0 . and γ = 5 / , the structures have a wavelength of ≈ H ρ , but thisis at a depth where H ρ ≈ . H ρ . Thus, locally we have a wave- Table 1.
Comparison of the optimal depth z B and the corre-sponding normalized magnetic field strength β ( z B ) for threevalues of γ for imposed horizontal and vertical magnetic fieldsof normalized strengths β = B /B eq0 , for Model I. horizontal field vertical fieldMod γ β z B /H ρ β ( z B ) β z B /H ρ β ( z B ) I 1 0.05 − . . . . . − . . Fig. 16.
Snapshots from MFS showing B z on the periphery of the computational domain for γ = 5 / and β = 0 . at differenttimes for Model I for the case of a vertical field. Fig. 17.
Similar to Fig. 16 of MFS, but for Model II at times similar to those in the DNS of Fig. 5. There are now more structuresthan in the earlier MFS of Fig. 16, and they develop more rapidly.length of ≈ H ρ . The situation is similar in the next panel ofFig. 4 where the wavelength is ≈ H ρ , and the structures are ata depth where H ρ ≈ . H ρ , so locally we have a wavelengthof ≈ H ρ . We can therefore conclude that our earlier results forisothermal layers can still be applied locally to polytropic layers.A new aspect, however, that was not yet anticipated at thetime, concerns the importance of NEMPI for vertical fields.While NEMPI with horizontal magnetic field still leads to down-flows in the nonlinear regime (the “potato-sack” effect), ourpresent work now confirms that structures consisting of verticalfields do not sink, but reach a strength comparable to or in ex-cess of the equipartition value (Brandenburg et al., 2013, 2014).This makes NEMPI a viable mechanism for spontaneouslyproducing magnetic spots in the surface layers. Our presentstudy therefore supports ideas about a shallow origin for activeregions and sunspots (Brandenburg, 2005; Brandenburg et al.,2010; Kitiashvili et al., 2010; Stein & Nordlund, 2012), con-trary to common thinking that sunspots form near the bottomof the convective zone (Parker, 1975; Spiegel & Weiss, 1980;D’Silva & Choudhuri, 1993). More specifically, the studies ofLosada et al. (2013) point toward the possibility that magneticflux concentrations form in the top , i.e., in the upper partof the supergranulation layer.There are obviously many other issues of NEMPI that needto be understood before it can be applied in a meaningful wayto the formation of active regions and sunspots. One questionis whether the hydrogen ionization layer and the resulting H − opacity in the upper layers of the Sun are important in provid-ing a sharp temperature drop and whether this would enhancethe growth rate of NEMPI, just like strong density stratifica-tion does. Another important question concerns the relevance of a radiating surface, which also enhances the density con-trast. Finally, of course, one needs to verify that the assump-tion of forced turbulence is useful in representing stellar con-vection. Many groups have considered magnetic flux concen-trations using realistic turbulent convection (Kitiashvili et al.,2010; Rempel, 2011; Stein & Nordlund, 2012). However, only atsufficiently large resolution can one expect strong enough scaleseparation between the scale of the smallest eddies and the sizeof magnetic structures. That is why forced turbulence has an ad-vantage over convection. Ultimately, however, such assumptionsshould no longer be necessary. On the other hand, if scale sepa-ration is poor, our present parameterization might no longer beaccurate enough, and one would need to replace the multipli-cation between q p and B in Eq. (25) by a convolution. Thispossibility is fairly speculative and requires a separate investiga-tion. Nevertheless, in spite of these issues, it is important to em-phasize that the qualitative agreement between DNS and MFS isalready surprisingly good. Acknowledgements.
This work was supported in part by the European ResearchCouncil under the AstroDyn Research Project No. 227952, by the SwedishResearch Council under the project grants 621-2011-5076 and 2012-5797 (IRL,AB), by EU COST Action MP0806, by the European Research Council under theAtmospheric Research Project No. 227915, and by a grant from the Governmentof the Russian Federation under contract No. 11.G34.31.0048 (NK, IR). Weacknowledge the allocation of computing resources provided by the SwedishNational Allocations Committee at the Center for Parallel Computers at theRoyal Institute of Technology in Stockholm and the National SupercomputerCenters in Link¨oping, the High Performance Computing Center North in Ume˚a,and the Nordic High Performance Computing Center in Reykjavik.
References
Brandenburg, A. 2005, ApJ, 625, 539
11. R. Losada et al.: Magnetic flux concentrations in a polytropic atmosphere