Astronomische Nachrichten, 6 January 2017

Small-scale dynamos in simulations of weakly stratified

turbulent convection

P. J. Käpylä^{1,2,3}^{,}^{⋆}, M. J. Käpylä^{3,2}, and A. Brandenburg^{4,5,6,7}

1 Leibniz-Institut für Astrophysik Potsdam, An der Sternwarte 16, D-11482 Potsdam, Germany

2 ReSoLVE Centre of Excellence, Department of Computer Science, Aalto University, PO Box 15400, FI-00076 Aalto,

Finland 3 Max-Planck-Institut für Sonnensystemforschung, Justus-von-Liebig-Weg 3, D-37077 Göttingen, Germany 4 NORDITA, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, SE-10691 Stock-

holm, Sweden 5 Department of Astronomy, AlbaNova University Center, Stockholm University, SE-10691 Stockholm, Sweden 6 JILA and Department of Astrophysical and Planetary Sciences, Box 440, University of Colorado, Boulder, CO 80303,

USA 7 Laboratory for Atmospheric and Space Physics, 3665 Discovery Drive, Boulder, CO 80303, USA

Received XXXX, accepted XXXX

Published online XXXX

Key words convection – turbulence – Sun: magnetic fields – stars: magnetic fields

Small-scale dynamo action is often held responsible for the generation of quiet-Sun magnetic fields. We aim to determine the excitation conditions and saturation level of small-scale dynamos in non-rotating turbulent convection at low magnetic Prandtl numbers. We use high resolution direct numerical simulations (DNS) of weakly stratified turbulent convection. We find that the critical magnetic Reynolds number for dynamo excitation increases as the magnetic Prandtl number is decreased. As a function of the magnetic Reynolds number (Rm), the growth rate λ of the dynamo is consistent with Re^{1/2} scaling. No evidence for a logarithmic increase of growth rate with Reynolds number is found.

Copyright line will be provided by the publisher

1 Introduction

Magnetic fields are ubiquitous in astrophysical systems. These fields are in most cases thought to be generated by some sort of dynamo process, involving either turbulent fluid motions or MHD-instabilities. In dynamo theory (e.g. Brandenburg & Subramanian 2005; Krause & Rädler 1980; Rüdiger & Hollerbach 2004) a distinction is made between large-scale and small-scale dynamos where the former produce fields whose length scale is greater than the scale of fluid motions whereas in the latter the two are comparable.

Small-scale dynamos have been found in direct numerical simulations of various types of flows provided that the magnetic Reynolds number (Rm) exceeds a critical value (Rm_{c}). However, in many astrophysical conditions molecular kinematic viscosity and magnetic diffusivity are vastly different implying that their ratio, which is the magnetic Prandtl number (Pm), is either very small or very large. For example, in the Sun Pm = 10^{−3} …10^{−6}. Numerical simulations of forced turbulence and other idealised flows indicate that Rm_{c} increases as Pm is decreased (Iskakov et al. 2007; Ponty et al. 2004; Schekochihin et al. 2004, 2005, 2007). Theoretical studies indicate a similar trend with an asymptotic value for Rm_{c} when Pm → 0 (Rogachevskii & Kleeorin 1997).

Simulations of turbulent convection have also been able to produce small-scale dynamos (e.g. Brandenburg et al. 1996; Cattaneo 1999; Favier & Bushby 2012; Nordlund et al. 1992; Pietarila Graham et al. 2010). However, the Pm dependence of the convection driven small-scale dynamo is so far unexplored.

2 Model

The numerical model is the same as in Käpylä et al. (2010) but without imposed shear or rotation. We use a Cartesian domain with dimensions L_{x} = L_{y} = 5d and L_{z} = d with 0 < z < d, where d is the depth of the convectively unstable layer. We solve the set of equations of magnetohydrodynamics

\frac{∂A}{∂t} = U × B − ημ_{0}J,  (1)
\frac{D lnρ}{Dt} = −∇ · U,  (2)
\frac{DU}{Dt} = −\frac{1}{ρ}∇p + g + \frac{1}{ρ}J × B + \frac{1}{ρ}∇ · 2νρS,  (3)
\frac{De}{Dt} = −\frac{p}{ρ}∇ · U + \frac{1}{ρ}∇ · K∇T + 2νS^{2} + \frac{μ_{0}η}{ρ}J^{2},  (4)

⋆ Corresponding author: pkapyla@aip.de

where D/Dt = ∂/∂t + U · ∇ is the advective time derivative, A is the magnetic vector potential, B = ∇ × A is the magnetic field, and J = μ_{0}^{−1}∇×B is the current density, μ_{0} is the vacuum permeability, η and ν are the magnetic diffusivity and kinematic viscosity, respectively, K is the heat conductivity, ρ is the density, U is the velocity, and g = −gẑ is the gravitational acceleration. The fluid obeys an ideal gas law p = ρe(γ − 1), where p and e are the pressure and internal energy, respectively, and γ = c_{P}/c_{V} = 5/3 is the ratio of specific heats at constant pressure and volume, respectively. The specific internal energy per unit mass is related to the temperature via e = c_{V}T. The rate of strain tensor S is given by S_{ij} = \frac{1}{2}(U_{i,j} + U_{j,i}) − \frac{1}{3}δ_{ij}∇ · U. (5)

In order to exclude complications due to overshoot and compressibility we employ a weak stratification: the density difference between the top and the bottom of the domain is twenty per cent and the average Mach number is always less than 0.1. The stratification in the associated hydrostatic initial state can be described by a polytrope with index m = 1.

The horizontal boundaries are periodic. We keep the temperature fixed at the top and bottom boundaries. For the velocity we apply impenetrable, stress-free conditions according to

∂_{z}U_{x} = ∂_{z}U_{y} = U_{z} = 0.  (6)

For the magnetic field we use vertical field conditions

B_{x} = B_{y} = 0.  (7)

2.1 Units, nondimensional quantities, and parameters

Dimensionless quantities are obtained by setting

d = g = ρ_{0} = c_{P} = μ_{0} = 1 ,  (8)

where ρ_{0} is the density at z_{m} = \frac{1}{2}d. The units of length, time, velocity, density, specific entropy, and magnetic field are then

[x] = d , [t] = \sqrt{d/g}, [U] = \sqrt{dg} ,

[ρ] = ρ_{0} , [s] = c_{P} , [B] = \sqrt{dgρ_{0}μ_{0}} .  (9)

The simulations are controlled by the following dimensionless parameters: thermal and magnetic diffusion in comparison to viscosity are measured by the Prandtl numbers

Pr = \frac{ν}{χ_{0}} , Pm = \frac{ν}{η},  (10)

where χ_{0} = K/(c_{P}ρ_{0}) is the reference value of the thermal diffusion coefficient, measured in the middle of the layer, z_{m}, in the non-convecting initial state. The efficiency of convection is measured by the Rayleigh number

Ra = \frac{gd^{4}}{νχ_{0}} (−\frac{1}{c_{P}} \frac{ds}{dz})_{z_{m}} ,  (11)

again determined from the initial non-convecting state at z_{m}. The entropy gradient can be presented in terms of logarithmic temperature gradients

(−\frac{1}{c_{P}} \frac{ds}{dz})_{z_{m}} = \frac{∇ − ∇_{ad}}{H_{P}} ,  (12)

Table 1 Summary of the runs, here Ma = u_{rms}/ \sqrt{dg}.

Run Pm Ma Re Rm λ̃[10^{−4}] δλ̃[10^{−4}] grid
A1 1 0.073 23 23  −52  47 64^{3}
A2 1 0.068 54 54  62  18 128^{3}
A3 1 0.064 101 101  162  20 128^{3}
A4 1 0.061 193 193  273  7 256^{3}
A5 1 0.056 354 354  453  2 512^{3}
B1 0.5 0.069 54 27 −128  24 128^{3}
B2 0.5 0.064 102 51  −27  12 128^{3}
B3 0.5 0.060 191 95  44  11 256^{3}
B4 0.5 0.056 360 180  155  18 512^{3}
B5 0.5 0.052 666 333  357  17 1024^{3}
C1 0.25 0.060 190 47 −144  22 256^{3}
C2 0.25 0.056 358 90 −35 28 512^{3}
C3 0.25 0.054 431 108 −46  10 512^{3}
C4 0.25 0.051 654 163 −27  11 1024^{3}
D1 0.1 0.057 360 36 −237  56 512^{3}
D2 0.1 0.052 664 65 −139  77 1024^{3}
D3 0.1 0.050 1057 106 −10  94 1024^{3}

with ∇ = (∂ lnT/∂ lnp)_{z_{m}}, ∇_{ad} = 1 − 1/γ, and H_{P} being the pressure scale height at z = z_{m}.

The effects of viscosity and magnetic diffusion are quantified respectively by the fluid and magnetic Reynolds numbers

Re = \frac{u_{rms}}{νk_{f}}, Rm = \frac{u_{rms}}{ηk_{f}} = PmRe,  (13)

where u_{rms} is the root-mean-square (rms) value of the velocity, and k_{f} = 2π/d is the wavenumber corresponding to the depth of the convectively unstable layer. Furthermore, we define the Péclet number as

Pe = \frac{u_{rms}}{χ_{0}k_{f}} = PrRe.  (14)

In most of our simulations we keep Pr = 1 and thus Pe = Re. The size of error bars is estimated by dividing the time series into three equally long parts. The largest deviation of the average for each of the three parts from that over the full time series is taken to represent the error.

The simulations were performed using the Pencil Code^{1}, which uses sixth-order explicit finite differences in space and a third-order accurate time stepping method. We use resolutions ranging from 64^{3} to 1024^{3}.

3 Results

We perform several sets of runs where we keep the magnetic Prandtl number fixed and vary Re and Rm, see Table 1. The lower resolution (64^{3}, 128^{3}, and 256^{3}) runs were started from a non-convecting state described in the previous section, whereas runs at 512^{3} and 1024^{3} were remeshed from saturated snapshots at lower resolutions. After the convection has reached a statistically saturated state we introduce

1 https://pencil-code.github.com/

Fig. 1 Upper row: specific entropy s/c_{P} near the surface z/d = 0.98 for Re = Pe = 23, 54, 101, and 193. Lower row: Re = Pe = 354, 666, and 1057.

Fig. 2 Growth rate γ of the rms-magnetic field normalized by the inverse convective turnover time τ^{−1} = u_{rms}k_{f} as a function of the magnetic Reynolds number Rm. The different symbols denote runs with Pm = 1 (triangles), Pm = 0.5 (stars), Pm = 0.25 (crosses), and Pm = 0.1 (diamonds). The red horizontal line denotes marginal stability, and the black dotted lines are curves proportional to Rm^{1/2}.

a weak random magnetic field of the order of 10^{−6}B_{eq}. After an initial transient the growth rate of the total rms-magnetic field is measured from

λ = ⟨d lnB_{rms}/dt⟩_{t},  (15)

where ⟨⋯ ⟩_{t} denotes time averaging. In the runs where the dynamo is clearly above or below critical, a short time series (few tens of turnover times) is sufficient to measure a statistically significant values of λ. The runs near the excitation threshold need to be run significantly longer (hundreds of turnover times). For the highest resolution runs at low Pm this is not feasible due to the computational cost and thus the error bars for these runs are typically somewhat greater than in the low resolution or Pm = 1 runs.

3.1 Growth rate in the kinematic regime

Figure 2 shows the growth rate of the magnetic field as a function of Rm for four magnetic Prandtl numbers. For reference we plot curves of the form

γ/τ^{−1} ≡ γ̃ = A + B Rm^{1/2},  (16)

where τ = (u_{rms}k_{f})^{−1} is the convective turnover time, the value of the constant A changes as the magnetic Prandtl number is changed, and where B = 3.5 · 10^{−3} for all values of Pm. The parameter A is negative, so the solutions will always decay for small values of Rm, but they increase with increasing values of Pm; see Fig. 2.

We find that the normalized growth rate for a given Rm decreases as Pm is decreased. Surprisingly, γ̃appears to follow a Rm^{1/2} trend for each value of Pm even in the cases when a dynamo is not excited. Such dependence is predicted by theory for high Rm, i.e. far away from excitation (Kleeorin & Rogachevskii 2012). However, given the relatively large error bars, the Rm^{1/2} scaling near the excitation threshold can at this point be only suggestive and far from definite. Indeed, analytic theory yields a different scaling in this regime (Kleeorin & Rogachevskii 2012).

3.2 Dependence on the box size

The dependence of the growth rate of the convection driven small scale dynamo on the horizontal size of the domain

Fig. 3 Growth rate γ of the rms-magnetic field normalized by the inverse convective turnover time τ^{−1} = u_{rms}k_{f} as a function of the horizontal box size for Pm = 1.

Fig. 4 Power spectra of velocity and magnetic field as functions of k̃ = k/k_{1} near the top of the domain from Run C4. E_{M} has been multiplied by 10^{14} for visualization purposes. The dotted line show a k^{−5/3} scaling for reference. The inset shows the velocity power spectrum compensated by k^{5/3}.

and the presence of mesogranulation has been raised in a recent paper by (Bushby et al. 2012). Figure 3 shows the growth rate of the magnetic field as a function of the horizontal box size for magnetic Prandtl number unity. Deviations from the constant trend are found for L_{H}/d = 1.25, and for L_{H}/d = 0.5 no dynamo action is found. Our standard box size L_{H}/d = 5 is thus adequate and does not seem to suffer from the issues raised by (Bushby et al. 2012).

3.3 Saturation level

Figure 5 shows the saturation field strength for Pm = 0.1 …1 and where the magnetic Reynolds number varies from 170 to 207 due to increasing u_{rms} when Pm decreases.

Fig. 5 Saturation field strength for runs with Pm = 1 …0.1 (Rm = 170 …207) as a function of the fluid Reynolds number.

4 Conclusions

Acknowledgements. The computations were performed on the facilities hosted by CSC – IT Center for Science Ltd. in Espoo, Finland, who are administered by the Finnish Ministry of Education. We also 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 National Supercomputer Centers in Linköping. This work was supported in part by Academy of Finland grants 121431, 136189 (PJK), and 112020 (MJK), the European Research Council under the AstroDyn Research Project 227952 and the Swedish Research Council grant 621-2007-4064.

References

Brandenburg, A., Jennings, R. L., Nordlund, �., et al. 1996, J.

Fluid Mech., 306, 325 Brandenburg, A. & Subramanian, K. 2005, Phys. Rep., 417, 1 Bushby, P. J., Favier, B., Proctor, M. R. E., & Weiss, N. O. 2012,

Geophysical and Astrophysical Fluid Dynamics, 106, 508 Cattaneo, F. 1999, ApJ, 515, L39 Favier, B. & Bushby, P. J. 2012, Journal of Fluid Mechanics, 690,

262 Iskakov, A. B., Schekochihin, A. A., Cowley, S. C., McWilliams,

J. C., & Proctor, M. R. E. 2007, Physical Review Letters, 98, 208501

Käpylä, P. J., Korpi, M. J., & Brandenburg, A. 2010, MNRAS,

402, 1458 Kleeorin, N. & Rogachevskii, I. 2012, Phys. Scr, 86, 018404 Krause, F. & Rädler, K.-H. 1980, Mean-field Magnetohydrody-

namics and Dynamo Theory (Oxford: Pergamon Press) Nordlund, A., Brandenburg, A., Jennings, R. L., et al. 1992, ApJ,

392, 647 Pietarila Graham, J., Cameron, R., & Schüssler, M. 2010, ApJ,

714, 1606 Ponty, Y., Politano, H., & Pinton, J.-F. 2004, Physical Review Let-

ters, 92, 144503 Rogachevskii, I. & Kleeorin, N. 1997, Phys. Rev. E, 56, 417 Rüdiger, G. & Hollerbach, R. 2004, The Magnetic Universe: Geo-

physical and Astrophysical Dynamo Theory (Weinheim: Wiley-VCH)

Schekochihin, A. A., Cowley, S. C., Taylor, S. F., Maron, J. L., &

McWilliams, J. C. 2004, ApJ, 612, 276 Schekochihin, A. A., Haugen, N. E. L., Brandenburg, A., et al.

2005, ApJ, 625, L115 Schekochihin, A. A., Iskakov, A. B., Cowley, S. C., et al. 2007,

New Journal of Physics, 9, 300