Large-scale dynamos in rapidly rotating plane layer convection
P. J. Bushby, P. J. Käpylä, Y. Masada, A. Brandenburg, B. Favier, C. Guervilly, M. J. Käpylä
aa r X i v : . [ a s t r o - ph . S R ] J a n Astronomy&Astrophysicsmanuscript no. paper_arxiv c (cid:13)
ESO 2018January 16, 2018
Large-scale dynamos in rapidly rotating plane layer convection
P. J. Bushby , P. J. Käpylä , , , , Y. Masada , A. Brandenburg , , , , B. Favier , C. Guervilly , and M. J. Käpylä , School of Mathematics and Statistics, Newcastle University, Newcastle upon Tyne, NE1 7RU, UKe-mail: [email protected] Leibniz-Institut für Astrophysik Potsdam, An der Sternwarte 16, D-11482 Potsdam, Germany ReSoLVE Centre of Excellence, Department of Computer Science, Aalto University, PO Box 15400, FI-00076 Aalto, Finland Max-Planck-Institut für Sonnensystemforschung, Justus-von-Liebig-Weg 3, D-37077 Göttingen, Germany Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, SE-10691 Stockholm, Sweden Department of Physics and Astronomy, Aichi University of Education; Kariya, Aichi 446-8501, Japan JILA and Department of Astrophysical and Planetary Sciences, Box 440, University of Colorado, Boulder, CO 80303, USA Department of Astronomy, AlbaNova University Center, Stockholm University, SE-10691 Stockholm, Sweden Laboratory for Atmospheric and Space Physics, 3665 Discovery Drive, Boulder, CO 80303, USA Aix Marseille Univ, CNRS, Centrale Marseille, IRPHE UMR 7342, Marseille, FranceJanuary 16, 2018
ABSTRACT
Context.
Convectively-driven flows play a crucial role in the dynamo processes that are responsible for producing magnetic activityin stars and planets. It is still not fully understood why many astrophysical magnetic fields have a significant large-scale component.
Aims.
Our aim is to investigate the dynamo properties of compressible convection in a rapidly rotating Cartesian domain, focusingupon a parameter regime in which the underlying hydrodynamic flow is known to be unstable to a large-scale vortex instability.
Methods.
The governing equations of three-dimensional nonlinear magnetohydrodynamics (MHD) are solved numerically. Differentnumerical schemes are compared and we propose a possible benchmark case for other similar codes.
Results.
In keeping with previous related studies, we find that convection in this parameter regime can drive a large-scale dynamo.The components of the mean horizontal magnetic field oscillate, leading to a continuous overall rotation of the mean field. Whilst thelarge-scale vortex instability dominates the early evolution of the system, it is suppressed by the magnetic field and makes a negligiblecontribution to the mean electromotive force that is responsible for driving the large-scale dynamo. The cycle period of the dynamo iscomparable to the ohmic decay time, with longer cycles for dynamos in convective systems that are closer to onset. In these particularsimulations, large-scale dynamo action is found only when vertical magnetic field boundary conditions are adopted at the upper andlower boundaries. Strongly modulated large-scale dynamos are found at higher Rayleigh numbers, with periods of reduced activity(“grand minima”-like events) occurring during transient phases in which the large-scale vortex temporarily re-establishes itself, beforebeing suppressed again by the magnetic field.
Key words.
Convection – Dynamo – Instabilities – Magnetic fields – Magnetohydrodynamics (MHD) – Methods: numerical
1. Introduction
In a hydromagnetic dynamo the motions in an electrically con-ducting fluid continuously sustain a magnetic field against theaction of ohmic dissipation. Most astrophysical magnetic fields,including those in stars and planets, are dynamo-generated. Innumerical simulations, turbulent motions almost invariably pro-duce an intermittent magnetic field distribution that is correlatedon the scale of the flow. However, most astrophysical objects ex-hibit magnetism that is organised on much larger scales. Theselarge-scale fields might be steady or they may exhibit sometime dependence. In the case of the Sun, for example, obser-vations of surface magnetism (e.g. Stix 2002) indicate the pres-ence of a large-scale oscillatory magnetic field in the solar inte-rior that changes sign approximately every 11 years. Dependingupon their age and spectral type, similar magnetic activity cy-cles can also be observed in other stars (e.g. Brandenburg et al.1998). Whilst appearing to be comparatively steady on thesetime-scales, the Earth’s predominantly dipolar field does exhibitlong-term variations, occasionally even reversing its magneticpolarity (although rather irregular, a typical time span betweenreversals is of the order of × years, see, e.g., Jones 2011). Our understanding of the physical processes that are responsiblefor the production of large-scale magnetic fields in astrophysicalobjects relies heavily on mean-field dynamo theory (e.g. Moffatt1978). This approach, in which the small-scale physics is pa-rameterised in a plausible way, has had considerable success.However, despite much recent progress in this area, there areapparently contradictory findings, indicating that we still do notcompletely understand why large-scale magnetic fields appear tobe so ubiquitous in astrophysics.Convectively-driven flows are a feature of stellar and plan-etary interiors, where the effects of rotation can often play animportant dynamical role via the Coriolis force. In the rapidlyrotating limit, convective motions tend to be helical, leading tothe expectation of a strong α -effect (an important regenerativeterm for large-scale magnetic fields in mean-field dynamo the-ory, usually a parameterised effect in simplified mean-field mod-els). Many theoretical studies have therefore been motivated bythe question of whether rotationally influenced convective flowscan drive a large-scale dynamo in a fully self-consistent (i.e. non-parameterised) manner. Article number, page 1 of 16&Aproofs: manuscript no. paper_arxiv
Using the Boussinesq approximation, in which the elec-trically conducting fluid is assumed to be incompressible,Childress & Soward (1972) were the first to demonstrate thata rapidly rotating plane layer of weakly convecting fluidwas capable of sustaining a large-scale dynamo (see alsoSoward 1974). To be clear about terminology (here, andin what follows), we describe a plane layer dynamo as“large-scale” if it produces a magnetic field with a sig-nificant horizontally averaged component (such a magneticfield can also be described as “system-scale”, see, e.g.,Tobias et al. 2011). Building on the work of Childress & Soward(1972) and Soward (1974), many subsequent studies have ex-plored the dynamo properties of related Cartesian Boussinesqmodels (Fautrelle & Childress 1982; Meneguzzi & Pouquet1989; St. Pierre 1993; Jones & Roberts 2000; Rotvig & Jones2002; Stellmach & Hansen 2004; Cattaneo & Hughes 2006;Favier & Proctor 2013; Calkins et al. 2015) as well as the cor-responding weakly-stratified system (Mizerski & Tobias 2013).However, whilst near-onset rapidly rotating convection doesproduce a large-scale dynamo (see, e.g., Stellmach & Hansen2004), only small-scale dynamo action is observed in the rapidlyrotating turbulent regime (Cattaneo & Hughes 2006), contraryto the predictions of mean-field dynamo theory. Indeed, Tilgner(2014) was able to identify an approximate parametric threshold(based on the Ekman number and the magnetic Reynolds num-ber) above which small-scale dynamo action is preferred: lowlevels of turbulence and a rapid rotation rate are found to be es-sential for a large-scale convectively-driven dynamo.There have been many fewer studies of the correspondingfully compressible system, so parameter space has not yet beenexplored to the same extent in this case. Käpylä et al. (2009)were the first to demonstrate that it is possible to excite a large-scale dynamo in rapidly rotating compressible convection atmodestly supercritical values of the Rayleigh number (the keyparameter controlling the vigour of the convective motions).However, it again appears to be much more difficult to drive alarge-scale dynamo further from convective onset, with small-scale dynamos typically being reported in this comparativelyturbulent regime (Favier & Bushby 2013). This is in agreementwith the Boussinesq studies, such as Cattaneo & Hughes (2006).Moreover, as noted by Guervilly et al. (2015), the transition fromlarge-scale to small-scale dynamos in these compressible sys-tems appears to occur in a similar region of parameter space tothe Boussinesq transition that was identified by Tilgner (2014).In rapidly rotating Cartesian domains, hydrodynamic con-vective flows just above onset are characterised by a small hor-izontal spatial scale (Chandrasekhar 1961). However, at slightlyhigher levels of convective driving, these small-scale motionscan become unstable, leading to a large-scale vortical flow (Chan2007). The width of these large-scale vortices is limited by thesize of the computational domain; in such simulations, the cor-responding flow field has a negligible horizontal average (un-like the large-scale magnetic fields described above). With in-creasing rotation rate, Chan (2007) observed a transition fromcooler cyclonic vortices to warmer anticyclones, and subse-quent fully compressible studies have found similar behaviour(Mantere et al. 2011; Käpylä et al. 2011; Chan & Mayr 2013). Incorresponding Boussinesq calculations, Favier et al. (2014) andGuervilly et al. (2014) found a clear preference for cyclonic vor-tices, and dominant anticyclones are never observed, althoughStellmach et al. (2014) did find states consisting of cyclonesand anticyclones (of comparable magnitude) at higher rotationrates. As noted by Stellmach et al. (2014) and Kunnen et al.(2016), this large-scale vortex instability is inhibited when no- slip boundary conditions are adopted at the upper and lowerboundaries, so the formation of large-scale vortices depends tosome extent upon the use of stress-free boundary conditions.From the point of view of the convective dynamo problem,this large-scale vortex instability can play a very important role.In the rapidly rotating Boussinesq dynamo of Guervilly et al.(2015), the large-scale vortex leads directly to the production ofmagnetic fields at a horizontal wavenumber comparable to thatof the large-scale vortex. Although the total magnetic energy inthis case is less than of the total kinetic energy, the resultantmagnetic field is locally strong enough to inhibit the large-scaleflow. This temporarily suppresses the dynamo until the magneticfield becomes weak enough for the vortical instability to growagain. In some sense, the dynamo in this case switches on andoff as the energy in the large-scale flow fluctuates.In the fully compressible regime, the large-scale vortex in-stability has been shown to produce a different type of dynamo(Käpylä et al. 2013; Masada & Sano 2014a,b). As in the Boussi-nesq case that was considered by Guervilly et al. (2015), thelarge-scale flow produces a large-scale magnetic field which ex-hibits some time-dependence. However, whilst the large-scalevortex is again suppressed once the magnetic field becomes dy-namically significant, these dynamos are able to persist with-out the subsequent regeneration of these vortices (suggestingthat these dynamos may be a compressible analogue of the dy-namo considered by Childress & Soward 1972, albeit operatingin the strong field limit). Once established, the magnetic energyin these dynamos is comparable to the kinetic energy of the sys-tem, with the horizontal components of the large-scale magneticfield oscillating in a regular manner, with a phase shift of ap-proximately π/ between the two components, leading to a netrotation of the mean horizontal magnetic field. Although eachcomponent of the mean magnetic field certainly oscillates, be-cause the temporal variation in the mean field essentially takesthe form of a global rotation of the field orientation, it should bekept in mind that in a suitably corotating frame the mean fieldwould appear statistically stationary.The large-scale dynamo that was found by Käpylä et al.(2013) and Masada & Sano (2014a,b) is arguably the simplestknown example of a moderately supercritical convectively-driven dynamo in a rapidly-rotating Cartesian domain. However,to achieve this nonlinear magnetohydrodynamical state, any nu-merical code must successfully reproduce the large-scale vortexinstability of rapidly rotating hydrodynamic convection in orderto amplify a weak seed magnetic field. In the nonlinear regimeof the dynamo, the resultant large-scale magnetic field must thenbe sustained at a level that is (approximately) in equipartitionwith the local convective motions. As a result, this dynamo is anexcellent candidate for a benchmarking exercise. Correspond-ing benchmarks exist for convectively-driven dynamos in spher-ical geometry, both for Boussinesq (Christensen et al. 2001;Marti et al. 2014) and for anelastic fluids (Jones et al. 2011). Tothe best of our knowledge, there is no similar benchmark for afully compressible, turbulent, large-scale dynamo.The main aim of this paper is to further investigate the prop-erties of this large-scale dynamo, focusing particularly upon theeffects of varying the rotation rate and the convective driving,as well as the size of the computational domain. We will estab-lish the regions of parameter space in which this dynamo canbe sustained, looking at the ways in which the dynamo ampli-tude and cycle period depend upon the key parameters of thesystem. Most significantly, we will show that it is possible to in-duce strong temporal modulation in large-scale dynamos of thistype by increasing the level of convective driving at fixed rota- Article number, page 2 of 16. J. Bushby et al.: Large-scale dynamos in rapidly rotating convection tion rate. Finally, we carry out a preliminary code comparison(confirming the accuracy and validity of one particular solutionvia three independent codes) to assess whether or not this systemcould form the basis of a nonlinear Cartesian dynamo bench-mark, possibly involving broader participation from the dynamocommunity. In the next section, we set out our model and de-scribe the numerical codes. Our numerical results, are discussedin Sect. 3. In the final section, we present our conclusions. Thestrengths and weaknesses of the proposed benchmark solutionare described in Appendix A.
2. Governing equations and numerical methods
We consider a plane layer of electrically conducting, compress-ible fluid, which is assumed to occupy a Cartesian domain ofdimensions ≤ x ≤ λd , ≤ y ≤ λd and ≤ z ≤ d ,where λ is the aspect ratio. This layer of fluid is heated frombelow, and the whole domain rotates rigidly about the verticalaxis with constant angular velocity Ω = Ω e z . We define µ tobe the dynamical viscosity of the fluid, whilst K is the radiativeheat conductivity, η is the magnetic diffusivity, µ is the vacuumpermeability, whilst c P and c V are the specific heat capacities atconstant pressure and volume respectively (as usual, we define γ = c P /c V ) . All of these parameters are assumed to be constant,as is the gravitational acceleration g = − g e z ( z increases up-wards). The evolution of this system is then determined by theequations of compressible magnetohydrodynamics, which canbe expressed in the following form, ∂ A ∂t = U × B − ηµ J , (1) D ln ρDt = − ∇ · U , (2) D U Dt = g − Ω × U + 1 ρ (2 µ ∇ · S − ∇ p + J × B ) , (3) T DsDt = 1 ρ (cid:0) K ∇ T + 2 µ S + µ η J (cid:1) , (4)where A is the magnetic vector potential, U is the velocity, B = ∇ × A is the magnetic field, J = ∇ × B /µ is thecurrent density, ρ is the density, s is the specific entropy, T is thetemperature, p is the pressure and D/Dt = ∂/∂t + U · ∇ de-notes the advective time derivative. The fluid obeys the ideal gaslaw with p = ( γ − ρe , where e = c V T is the internal energy.The traceless rate of strain tensor S is given by S ij = ( U i,j + U j,i ) − δ ij ∇ · U , (5)whilst the magnetic field satisfies ∇ · B = 0 . (6)Stress-free impenetrable boundary conditions are used for thevelocity, U x,z = U y,z = U z = 0 on z = 0 , d, (7)and vertical field conditions for the magnetic field, i.e. B x = B y = 0 on z = 0 , d, (8)respectively ( ∇ · B = 0 then implies B z,z = 0 at z = 0 , d ).The temperature is fixed at the upper and lower boundaries. Weadopt periodic boundary conditions for all variables in each ofthe two horizontal directions. Dimensionless quantities are obtained by setting d = g = ρ m = c P = µ = 1 , (9)where ρ m is the initial density at z = z m = 0 . d . The units oflength, time, velocity, density, entropy, and magnetic field are [ x ] = d , [ t ] = p d/g , [ U ] = p dg , [ ρ ] = ρ m , (10) [ s ] = c P , [ B ] = p dgρ m µ . Having non-dimensionalised these equations, the behaviour ofthe system is determined by various dimensionless parameters.Quantifying the two key diffusivity ratios, the fluid and magneticPrandtl numbers are given by
Pr = ν m χ m , Pm = ν m η , (11)where ν m = µ/ρ m is the mean kinematic viscosity and χ m = K/ ( ρ m c P ) is the mean thermal diffusivity. Defining H P to bethe pressure scale height at z m , the midlayer entropy gradient inthe absence of motion is (cid:18) − c P d s d z (cid:19) m = ∇ − ∇ ad H P , (12)where ∇ − ∇ ad is the superadiabatic temperature gradient with ∇ ad = 1 − /γ and ∇ = ( ∂ ln T /∂ ln p ) z m . The strength ofthe convective driving can then be characterised by the Rayleighnumber, Ra = gd ν m χ m (cid:18) − c P d s d z (cid:19) m = gd ν m χ m (cid:18) ∇ − ∇ ad H P (cid:19) . (13)The amount of rotation is quantified by the Taylor number, Ta = 4Ω d ν (14)(which is related to the Ekman number, Ek , by Ek = Ta − / ).Since the critical Rayleigh number for the onset of hydrody-namic convection is proportional to Ta / in the rapidly rotatingregime (Chandrasekhar 1961), it is also useful to consider thequantity f Ra = RaTa / (cid:16) = Ra Ek / (cid:17) (15)(see, e.g., Julien et al. 2012). This rescaled Rayleigh number isa measure of the supercriticality of the convection that takes ac-count of the stabilising influence of rotation. We also quote theconvective Rossby number Ro c = (cid:18) RaPrTa (cid:19) / , (16)which is indicative of the strength of the thermal forcing com-pared to the effects of rotation.Whilst Ra , Ta and Pr are input parameters that must be spec-ified at the start of each simulation, it is possible to measure anumber of useful quantities based on system outputs. These areexpressed here in dimensional form for ease of reference. Wedefine the fluid and magnetic Reynolds numbers via Re = u rms ν m k f , Rm = u rms ηk f , (17) Article number, page 3 of 16&Aproofs: manuscript no. paper_arxiv where k f = 2 π/d is indicative of the vertical scale of varia-tion of the convective motions, and u rms is the time-average ofthe rms velocity during the saturated phase of the dynamo. Thetime-evolution of the rms velocity, U rms ( t ) , is also considered,but only its constant time-averaged value will be used to defineother diagnostic quantities. The quantity u rms k f is therefore anestimate of the inverse convective turnover time in the nonlinearphase of the dynamo. The Coriolis number, an alternative mea-sure of the importance of rotation (compared to inertial effects)is given by Co = 2 Ω u rms k f ≡ Ta / π Re . (18)All of the simulations described in this paper have Co & (and Ro c < . ) so are in a rotationally dominated regime . Finally,the equipartition magnetic field strength is defined by B eq ≡ h µ ρ U i / , (19)where angle brackets denote volume averaging. All of the simulations in this paper are initialised from a hy-drostatic state corresponding to a polytropic layer, for which p ∝ ρ /m , where m is the polytropic index. Assuming amonatomic gas with γ = 5 / we adopt a polytropic index of m = 1 throughout. This gives a superadiabatic temperature gra-dient of ∇ − ∇ ad = 1 / , so the layer is convectively unstable,as required. The degree of stratification is determined by speci-fying a density contrast of across the layer. To be as clear aspossible about our proposed benchmark case (see Appendix Afor details), it is useful to provide explicit functional forms forthe initial density, pressure and temperature profiles in our di-mensionless units. Recalling that the layer has a unit depth andthat ρ m = 1 , the initial density profile is given by ρ ( z ) = 25 (4 − z ) , whilst the initial pressure and temperature profiles are given by p ( z ) = 115 (4 − z ) and T ( z ) = 512 (4 − z ) , respectively. The fixed temperature boundary conditions implythat T (0) = 5 / and T (1) = 5 / (independent of x and y , forall time). These profiles are consistent with a dimensionless pres-sure scale height of H P = 1 / at the top of the domain, which isan alternative way of specifying the level of stratification. Sincethe governing equations are formulated in terms of the specificentropy, it is worth noting that these initial conditions are con-sistent with an initial entropy distribution of the form s = ln (cid:18) T ( z ) T m (cid:19) −
25 ln (cid:18) p ( z ) p m (cid:19) , Note that (2 π Co) − = u rms / d is equivalent to the standardRossby number, based on the layer depth and the rms velocity. We use Co here so as not to confuse this quantity with the convective Rossbynumber, Ro c . where T m = 25 / and p m = 5 / are the midlayer values ofthe temperature and pressure respectively (note that s m = 0 withthis normalisation).In all simulations, convection is initialised by weakly per-turbing this polytropic state in the presence of a low amplitudeseed magnetic field (which varies over short length scales, withzero net flux across the domain). The precise details of theseinitial perturbations do not strongly influence the nature of thefinal nonlinear dynamos, although it goes without saying thatthe early evolution of this system does depend upon the initialconditions that are employed (as illustrated in Appendix A). The P
ENCIL C ODE (Code 1) is a tool for solving partial dif-ferential equations on massively parallel architectures. We useit in its default configuration in which the MHD equations aresolved as stated in Eqs. (1)–(4). First and second spatial deriva-tives are computed using explicit centred sixth-order finite dif-ferences. Advective derivatives of the form U · ∇ are com-puted using a fifth-order upwinding scheme, which correspondsto adding a sixth-order hyperdiffusivity with the diffusion co-efficient | U | δx / (Dobler et al. 2006). For the time steppingwe use the low-storage Runge-Kutta scheme of Williamson(1980). Boundary conditions are applied by setting ghost zonesoutside the physical boundaries and computing all derivativeson and near the boundary in the same fashion. The non-dimensionalising scalings that have been described above cor-respond directly to those employed by Code 1. Whilst the othercodes that have been used (see below) employ different scalings,we have rescaled results from these to ensure direct comparabil-ity with the results from the P ENCIL C ODE .The second code that shall be used in this paper (Code 2) isan updated version of the code that was originally described byMatthews et al. (1995). The system of equations that is solvedis entirely equivalent to those presented in Sect. 2.1, but insteadof time-stepping evolution equations for the magnetic vector po-tential, the logarithm of the density, the velocity field and thespecific entropy, this code solves directly for the density, veloc-ity field and temperature (see, e.g., Matthews et al. 1995), whilsta poloidal-toroidal decomposition is used for the magnetic field.Due to the periodicity in both horizontal directions, horizontalderivatives are computed in Fourier space using standard fastFourier transform libraries. In the vertical direction, a fourth-order finite difference scheme is used, adopting an upwind sten-cil for the advective terms. The time-stepping is performed byan explicit third-order Adams-Bashforth technique, with a vari-able time-step. As noted above, this code adopts a different setof non-dimensionalising scalings to those described in Sect. 2.2(see Favier & Bushby 2012, for more details).Code 3 is based upon a second-order Godunov-type finite-difference scheme that employs an approximate MHD Riemannsolver with operator splitting (Sano et al. 1999; Masada & Sano2014a,b). The hydrodynamical part of the equations is solved bya Godunov method, using the exact solution of the simplifiedMHD Riemann problem. The Riemann problem is simplified byincluding only the tangential component of the magnetic field.The characteristic velocity is then that of the magneto-sonicwave alone, and the MHD Riemann problem can be solved ina way similar to the hydrodynamical one (Colella & Woodward1984). The piecewise linear distributions of flow quantities arecalculated with a monotonicity constraint following the method https://github.com/pencil-code Article number, page 4 of 16. J. Bushby et al.: Large-scale dynamos in rapidly rotating convection
Fig. 1.
Time evolution of U rms ( t ) and B rms ( t ) for the reference solu-tion. Time has been normalised by the inverse convective turnover timein the final nonlinear state. of van Leer (1979). The remaining terms, the magnetic tensioncomponent of the equation of motion and the induction equation,are solved by the Consistent MoC-CT method (Clarke 1996),guaranteeing ∇ · B = 0 to within round-off error throughout thecalculation (Evans & Hawley 1988; Stone & Norman 1992).
3. Results
In this section, we present a detailed description of a typical dy-namo run exhibiting large-scale dynamo action (carried out us-ing Code 2). This case, which will henceforth be referred to asthe reference solution, forms the basis for the proposed bench-mark calculation that is discussed in Appendix A (with the corre-sponding calculations denoted by cases A1, A2 and A3 in the up-per three rows of Table 1, where our parameter choices are alsosummarised). For numerical convenience, we fix
Pr = Pm = 1 for the reference solution, whilst our choice of
Ta = 5 · en-sures that rotation plays a significant role in the ensuing dynam-ics. Defining Ra crit to be the critical Rayleigh number at whichthe hydrostatic polytropic state becomes linearly unstable to con-vective perturbations, it can be shown that Ra crit = 6 . · for this set of parameters. This critical value (here quoted to threedecimal places) was determined using an independent Newton-Raphson-Kantorovich boundary value solver for the correspond-ing linearised system. At onset, the critical horizontal wavenum-ber for the convective motions is very large (in these dimension-less units, the preferred horizontal wavenumber at onset is givenby k ≈ ), indicating a preference for narrow convective cells.Our choice of Rayleigh number, Ra = 2 . · ≈ crit , ismoderately supercritical. However, we would still expect small-scale convective motions during the early phases of the simu-lation. Our choice of λ = 2 for the aspect ratio of the domainis small enough to enable us to properly resolve these motions,whilst also being large enough to allow the expected large-scalevortex instability to grow. Figure 1 shows the time evolution of U rms ( t ) and B rms ( t ) ,the rms velocity and magnetic field, for the reference solution.Whilst the rms quantities have been left in dimensionless form,time has been normalised by the inverse convective turnover time Fig. 2.
Early stages of the reference solution. Mid-layer temperature dis-tributions (at z m ) at times tu rms k f = 3 . ( top ) and tu rms k f = 362 . ( middle ), plotted as functions of x and y . In each plot, the mid-layertemperature is normalised by its mean (horizontally averaged) value,with the contours showing the deviations from the mean. Bottom: mid-layer kinetic energy spectra at the same times, normalised by the av-erage kinetic energy at z m at each time. Here k is the horizontalwavenumber (in the x direction), whilst k = 2 π/λ is the fundamen-tal mode. The inset highlights the low wavenumber behaviour, using anon-logarithmic abscissa to include k/k = 0 . during the final stages of the simulation, by which time the sys-tem is in a statistically steady state. Initially, there is a very briefperiod (spanning no more than a few convective turnover times)of exponential growth, followed by a similarly rapid reduction inthe magnitude of U rms ( t ) . At these very early times, the convec- Article number, page 5 of 16&Aproofs: manuscript no. paper_arxiv
Table 1.
Summary of the simulations in this paper.
Case Code Grid λ Pr Pm Ta Ra f Ra Ro c Re Rm Co
Large-scaleA1 1 · . · · . · · . · · . · · · × · · · · · · · · × · · × · · × · · × · . · · . · × · . · × · . · × . · . · · . · · . · · . · · . · Notes.
A1, A2 and A3 correspond to the reference solution (a detailed comparison of these cases is presented in Appendix A). Simulations B1–8were initialised from a polytropic state, whilst simulations C1–5 were started from case A2, gradually varying parameters along this solutionbranch. In Sets D and E the convective Rossby number is fixed at Ro c = 0 . . In E1 the density contrast is . in comparison to in the othersets. All input parameters are defined in the text. We recall that the Reynolds numbers, Re and Rm , and the Coriolis number, Co , are measuredduring the saturated phase of the dynamo (and the quoted values are time-averaged over this phase). The final column indicates whether or notthere is a large-scale dynamo. tive motions are (as expected from linear theory) characterisedby a small horizontal length scale. The solution then enters aprolonged phase of more gradual evolution, lasting until time tu rms k f ≈ , during which time U rms ( t ) tends to increase.This period of growth coincides with a gradual transfer of en-ergy from small to large scales, leading eventually to a convec-tive state that is dominated by larger scales of motion. This is anexample of the large-scale vortex instability that was discussedin Sect. 1.The transition from small to large-scale convection is clearlyillustrated by the time evolution of the corresponding tempera-ture distribution. The upper two plots in Fig. 2 show the mid-layer temperature distribution at time tu rms k f = 3 . , at whichpoint the motions have a small horizontal length scale, and time tu rms k f = 362 . , at which the temperature distribution variesover the largest available scale. The lower plot shows the corre-sponding kinetic energy spectra as a function of the horizontalwavenumber, k , at times tu rms k f = 3 . and tu rms k f = 362 . .As expected from linear theory, most of the power at time tu rms k f = 3 . is concentrated at relatively high wavenumbers.Defining k = 2 π/λ (in this domain, k = π ) to be the fun-damental mode, corresponding to one full oscillation across thewidth of the domain, the corresponding kinetic energy spectrumpeaks at around k/k = 9 at this early time. This implies afavoured horizontal wavenumber of k ≈ . in these dimen-sionless units. This is comparable to the critical wavenumber atconvective onset. At time tu rms k f = 362 . , there is still sig-nificant energy in the higher wavenumber, small-scale compo-nents of the flow (in fact, the spectrum is now broader, extend-ing to higher wavenumbers than before). However, most of the power in the kinetic energy spectrum has been transferred to the k/k = 1 component. This corresponds to the large-scale vor-tex. Even at tu rms k f = 3 . there are some indications of non-monotonicity in the power spectrum at low wavenumbers, so thisinstability sets in very rapidly as the system evolves.As U rms ( t ) (or, equivalently, the kinetic energy in the sys-tem) increases, it soon reaches a level above which the flow issufficiently vigorous to excite a dynamo. The seed magnetic fieldis initially very weak, with B rms ( t ) many orders of magnitudesmaller than the typical values of U rms ( t ) . However, as illus-trated in Fig. 1, once the large-scale vortex instability starts togrow, B rms ( t ) also starts to increase. With increasing U rms ( t ) (effectively, an ever-increasing magnetic Reynolds number dur-ing this growth phase), the growth of B rms ( t ) accelerates until itreaches a point at which the magnetic field is strong enough toexert a dynamical influence upon the flow. The sharp decrease in U rms ( t ) at time tu rms k f ≈ is an indicator of the effects ofthe Lorentz force upon the flow, with the large-scale mode be-ing strongly suppressed by the magnetic field. After this point,which marks the start of the nonlinear phase of the dynamo, amore gradual growth of B rms ( t ) is accompanied by a gradualdecrease in U rms ( t ) . By the time tu rms k f ≈ , the dynamoappears to reach a statistically-steady state in which B rms ( t ) iscomparable to U rms ( t ) , which itself is comparable to the rms ve-locity before the large-scale vortex started to develop. The ohmicdecay time, τ η (based on the depth of the layer), corresponds to ≈ convective turnover times. So having continued this cal-culation to tu rms k f ≈ , observing no significant changes in B rms or u rms over the latter half of the simulation, we can be Article number, page 6 of 16. J. Bushby et al.: Large-scale dynamos in rapidly rotating convection
Fig. 3.
Reference solution dynamo.
Top : normalised volume averagesof the squared horizontal magnetic field components, ˜ B x /B (blue)and ˜ B y /B (red) as functions of time; the black dotted line shows ( ˜ B x + ˜ B y ) /B . Middle : horizontally averaged horizontal magneticfields (normalised by the equipartition field strength) as functions of z and time. Bottom : mid-layer kinetic (black) and magnetic (red) energyspectra, normalised by the average kinetic energy at z m at the end of thedynamo run. The inset highlights the low wavenumber behaviour of thespectra, using a non-logarithmic abscissa to include k/k = 0 . very confident that this is a persistent dynamo solution that willnot decay over longer times.Having discussed the time evolution of the system, we nowturn our attention to the form of the magnetic field. Figure 3shows the time- and z -dependence of the horizontally aver-aged profiles for B x and B y , as well as the kinetic and mag- netic energy spectra at the end of the dynamo run. Here the B i (for i = x or y ) represent horizontally averaged magneticfield components, whilst the ˜ B i correspond to volume aver-ages. It is clear from these plots that the magnetic field distri-bution has no significant mean horizontal component (comparedto its equipartition value) when the dynamo enters the nonlin-ear phase (at tu rms k f ≈ ). However, as the dynamo evolvesfrom tu rms k f ≈ , a cyclically-varying mean horizontal fieldgradually emerges, eventually growing to a level at which thepeak mean horizontal field almost reaches the equipartition level(although it should be noted that the peak amplitude does varysomewhat, from one cycle to the next). This mean magnetic fieldis approximately symmetric about the mid-plane with no real in-dication of any propagation of activity either towards or awayfrom this mid-plane. The cycle period, τ cyc , is approximately convective turnover times, which is of a similar order ofmagnitude to the ohmic decay time quoted above ( τ cyc ≈ . τ η ).The dominance of the large-scale magnetic field is confirmed bythe presence of a pronounced peak at k/k = 0 in the mag-netic energy spectrum that is shown in the lower part of Fig. 3.It should be stressed again that the kinetic energy spectrum ispeaked at small scales at this stage, so the large-scale vortex(which plays a critical role in initialising the dynamo) no longerappears to have a significant role to play in this large-scale dy-namo. This is emphasised by the lack of a significant k/k = 1 peak in the kinetic energy spectrum. Defining E = U × B , and recalling that an overbar denotes ahorizontal average, it can be shown that ∂B x ∂t = − ∂ E y ∂z + η ∂ B x ∂z , (20)whilst ∂B y ∂t = ∂ E x ∂z + η ∂ B y ∂z . (21)With the given boundary and initial conditions, it is straightfor-ward to show that B z = 0 for all z and t . In the absence of asignificant mean flow, standard mean-field dynamo theory (e.g.Moffatt 1978) suggests that there should be a linear relationshipbetween the components of the mean electromotive force, E , andthe components and z -derivatives of B . To be more specific, wemight expect relations of the form E x ( z ) = α ( z ) B x ( z ) + β ( z ) ∂B y ∂z , E y ( z ) = α ( z ) B y ( z ) − β ( z ) ∂B x ∂z , (22)where α ( z ) and β ( z ) are scalar functions of height only. Whilstthe β ( z ) term would dominate near the boundaries, where themean magnetic field vanishes (and the field gradient is large),we might expect the α ( z ) term to dominate in a non-negligibleregion in the vicinity of the mid-plane, where the mean field gra-dient is relatively small. Having said that, it should be stressedthat the flow that is responsible for driving this large-scale dy-namo is the product of a fully nonlinear magnetohydrodynamicstate. Even before the large-scale magnetic field starts to grow,this flow is strongly influenced by magnetic effects. As noted byCourvoisier et al. (2009), this probably implies that the relation-ship between E and B is more complicated than that suggested Article number, page 7 of 16&Aproofs: manuscript no. paper_arxiv
Fig. 4.
Time and z dependence of B x ( top ) and E x ( middle ) for a dy-namo simulation with the reference solution parameters. Note that E x has been normalised by the product of u rms ( . ) and the correspond-ing rms magnetic field, b rms ( . ), both of which have been time-averaged over the given time period. Bottom: the results of temporallysmoothing E x , using a sliding (boxcar) average with a width of approx-imately convective turnover times. by Eq. (22). We defer detailed considerations of this question tofuture work, focusing here upon the form of E and the compo-nents of the flow that are responsible for producing it.Figure 4 shows contour plots of E x and (for comparison) B x as functions of z and t , taken from a simulation that duplicatesthe reference solution parameters, but is evolved from a statewith a weaker initial thermal perturbation. A quantitatively com-parable solution is obtained, indicating that this large-scale dy-namo is robust to such changes to the initial conditions. Clearly E x exhibits much stronger temporal fluctuations than B x , whilst E x is more asymmetric about the mid-plane than B x . Neverthe-less, in purely qualitative terms, there are some indications of asimple temporal correlation between E x and B x , with a similarlong time-scale variation apparent in both cases. This is not in-consistent with the relation suggested by Eq. (22). Although notshown here, B y and E y evolve in a similar way.As we have already noted, the large-scale vortex is stronglysuppressed by the magnetic field. To confirm that the lowwavenumber components of the flow make a negligible contribu-tion to this large-scale dynamo, we have investigated the effectsof removing the low wavenumber content from the flow beforecalculating E . It should be emphasised that we are not makingany changes to the dynamo calculation itself – all filtering is car-ried out in post-processing. The upper part of Fig. 5 shows theeffects of removing all Fourier modes with k/k ≤ beforecalculating E x (although not shown here, similar results are ob-tained for E y ). At least in qualitative terms, it is clear that this iscomparable to the corresponding unfiltered E x that is shown inFig. 4. A much more significant change is observed when modeswith k/k ≤ are removed: this filtered E x appears to be al- Fig. 5.
Contours of E x (again normalised by u rms b rms ) as a function of z and time, derived (in post-processing) from a filtered flow. Top: hori-zontal wavenumbers with k/k ≤ have been removed from the flowbefore calculating E x . Bottom: the filtering threshold is set at k/k ≤ . Fig. 6.
Time-averaged ( E x + E y ) / , normalised by u rms b rms , as afunction of z (solid black curve). The blue (dash-dotted) curve (veryclose to the black line) shows the same quantity but having removedlow horizontal wavenumber components ( k/k ≤ ) from the flow be-fore calculating E x and E y . The red (dashed) curve shows the effects ofremoving more Fourier modes ( k/k ≤ ). most antisymmetric about the mid-plane, with some propaga-tion of activity away from the mid-plane that is not apparent inthe reference solution. Intriguingly, there is a significant fractionof the domain over which this filtering process actually changesthe sign of the resulting E x (compared to the corresponding un-filtered case). This sign change occurs gradually as the filteringthreshold is increased from k/k = 2 to , so there is no abrupttransition. Further increments to this filtering threshold lead tono further qualitative changes in the distribution of E x , althoughits amplitude obviously decreases as more components are re-moved from the flow.Figure 6 shows the time-averaged values of ( E x + E y ) / as a function of z , for the filtered and unfiltered cases. Clearlythe first filtered case, with only the k/k ≤ modes removed,is quantitatively comparable to the unfiltered case, which rein-forces the notion that the low wavenumber content of the flowplays no significant role in the dynamo. On the other hand, theclear quantitative differences between the two filtered cases indi- Article number, page 8 of 16. J. Bushby et al.: Large-scale dynamos in rapidly rotating convection cate that components of the flow with horizontal wavenumbers inthe range ≤ k/k ≤ play a crucial role in driving this large-scale dynamo. It is worth highlighting the fact that the small-scale motions alone seem to be capable of producing a coherent,systematically-varying E with a peak value that is of compara-ble magnitude to the peak value of the unfiltered E . This suggeststhat it may be possible for the small-scale motions alone to sus-tain a dynamo. Having said that, it should also be emphasisedthat these small-scale motions are not independent of the large-scale flows and magnetic field in the system (we stress againthat all filtering is carried out in “post-processing”), and so thesmall-scale motions are almost certainly strongly influenced bythe fact that a large-scale dynamo is operating. Whilst it is tempt-ing to speculate that there may be some connection here with thenear-onset dynamos of (e.g.) Stellmach & Hansen (2004), whichrely purely on small-scale motions, we have not yet been able todemonstrate this in a conclusive manner.Having identified the components of the flow that are respon-sible for driving the large-scale dynamo, we conclude this sec-tion with a brief comment on the role of the magnetic boundaryconditions. Integrating Eqs. (20) and (21) over z , it is straightfor-ward to verify that our magnetic boundary conditions (in which B x = B y = 0 at z = 0 and z = 1 ) allow the net horizontalflux to vary in time (c.f. the near-onset study of Favier & Proctor2013). In particular, these boundary conditions allow for the dif-fusive transport of magnetic flux out of the domain. Recallingthat the dynamo oscillates on a timescale that is of the same or-der of magnitude as the ohmic decay time across the layer, itis probable that the cycle period of the dynamo reflects the rateat which horizontal magnetic flux can be ejected from the do-main. Assuming that the simple mean-field ansatz of Eq. (22)is a reasonable description of E , this process would be acceler-ated by a positive β ( z ) at the boundaries (which would enhancediffusion), so we should not be surprised to find examples ofthis dynamo with a considerably shorter cycle period. Havingsaid that, even a large β ( z ) at the boundaries can only help ifthe boundary conditions allow it to do so. If we were to modifythe magnetic boundary conditions to the Soward-Childress (per-fect conductor) conditions of B z = ∂B y /∂z = ∂B x /∂z = 0 at z = 0 and z = 1 , that would mean that the total horizon-tal flux would be invariant (initially set to zero). Under thesecircumstances, we have confirmed that a simulation with the ref-erence dynamo parameters produces only a small-scale dynamo.Of course, we cannot exclude the possibility that there are re-gions of parameter space in which a large-scale dynamo can beexcited via this large-scale vortex mechanism, with these per-fectly conducting magnetic boundary conditions. However, wecan say that these simulations suggest that such a configurationwould be less favourable for large-scale dynamo action than thevertical conditions that we have adopted here.As a final remark, it is worth noting that one way of cir-cumventing the dependence of this dynamo upon the magneticboundary conditions is to include one or more convectively sta-ble layers into the system (Käpylä et al. 2013), which could re-side either above or below the convective layer. Even if perfectlyconducting boundary conditions are applied, such a stable regioncan act as a “flux repository” for the system: if the convectivelayer can expel a sufficient quantity of magnetic flux into thisregion, the large-scale dynamo could continue to operate as de-scribed above (as is the case for run D1e in Käpylä et al. 2013).Whilst we emphasise that this is not a model for the solar dy-namo, the possible importance of an underlying stable layer act-ing as a flux repository has also been recognised in that context.Certain solar dynamo models (e.g. Parker 1993; Tobias 1996a; Fig. 7. Ek as a function of Ra Ek / (= f Ra) from various studies inthe literature where large-scale dynamos were present (blue symbols) orabsent (red). Data is shown from Käpylä et al. (2009) ( ⋄ ), Käpylä et al.(2013) ( ✳ ), Stellmach & Hansen (2004) ( + ), Favier & Bushby (2013)( (cid:3) ), Masada & Sano (2014b, 2016) ( △ ), Cattaneo & Hughes (2006)( × ), and Guervilly et al. (2015) ( (cid:13) ). The green and orange diamondscorrespond to the present study with (green) or without (orange) large-scale dynamos. The shaded area shows the parameter region wherelarge-scale vortices are present in the hydrodynamic regime. Thedarker area corresponds to the large-scale vortex region identified byGuervilly et al. (2015), whereas the small grey diamonds and the star re-fer to the simulations of Favier et al. (2014) and Stellmach et al. (2014),respectively. MacGregor & Charbonneau 1997) rely on the assumption thatthe bulk of the large-scale toroidal magnetic flux in the solarinterior resides in the stable layer just below the base of the con-vection zone. This can therefore be regarded another example ofa situation in which the addition of a stable layer can be benefi-cial for dynamo action.
As indicated in Table 1, we have carried out a range of simula-tions to assess the sensitivity of the reference solution dynamo tovariations in the parameters. Cases B1–8 were all initialised fromour standard polytropic state, whilst cases C1–5 were initialisedfrom the reference solution (A2). Cases D1–3 investigated theeffects of increasing Ra and Ta at fixed convective Rossby num-ber, Ro c . Finally, case E1 is a pseudo-Boussinesq calculationwith the initial density varying linearly between . and . ; allother parameters (including the polytropic index) were identicalto the reference solution.When comparing dynamos at different rotation rates, it isconvenient to consider the modified Rayleigh number, f Ra =Ra / Ta / . For ease of comparison with previous studies, it isworth recalling that the Ekman number, Ek , is related to the Tay-lor number by Ek = Ta − / . So f Ra is equivalent to Ra Ek / ,which tends to be the more usual form of this parameter in thegeodynamo literature. Figure 1 of Guervilly et al. (2015) sug-gests that f Ra must exceed a value of approximately for thelarge-scale vortex instability to operate (obviously rapid rotationis also required). For the reference solution, f Ra ≈ , whichis certainly consistent with that picture. The other point to notefrom Fig. 1 of Guervilly et al. (2015) is that (in the absence ofa large-scale vortex instability) large-scale dynamos tend to berestricted to a small region of parameter space in which the layer Article number, page 9 of 16&Aproofs: manuscript no. paper_arxiv is rapidly-rotating (i.e. Ta > ; equivalently, Ek < − )and the convection is only weakly supercritical (typically, f Ra is O (10) ).Following a similar approach to that of Guervilly et al.(2015), Fig. 7 shows how the simulations that are reported inthis paper compare with others in the literature. These dynamosimulations are classified according to f Ra = Ra Ek / and Ek (the Ekman number has been used here, rather than the Taylornumber, for ease of comparability with Guervilly et al. 2015).The shaded area indicates the approximate region of parame-ter space where the large-scale vortex instability has been ob-served, although the limits of this region do also depend to someextent on the aspect ratio of the corresponding simulation do-mains. The parameter regime in the lower right part of the plot isalso likely to support large-scale vortices but numerical simula-tions in this regime are currently beyond the available computa-tional resources. It should be stressed that there are many differ-ent types of convective dynamos on this plot. Some simulationsare Boussinesq rather than compressible, others have perfectly-conducting magnetic boundary conditions, others feature under-lying and/or overlying stable layers. These model differences be-come particularly important at the edges of the large-scale vor-tex region. Focusing on the upper left-hand corner of the shadedregion, the single-layer calculation of Favier & Bushby (2013)was just outside of the large-scale vortex region and so found asmall-scale dynamo. The large-scale vortex instability was, how-ever, present in the multiple-layer model of Käpylä et al. (2013),which explains the observation of large-scale dynamo action inthat case. Having said that, regardless of the details of the model,it is clear that the large-scale vortex instability seems to providea route by which large-scale dynamos can be found in moder-ately supercritical, rapidly rotating convection, outside of theirnormal operative region of parameter space.It is worth emphasising at this stage that not all simulationsin the large-scale vortex parameter region produce large-scaledynamos. In particular, we only find small-scale dynamos at lowaspect ratios (i.e. λ < ). For λ = 1 (see cases B6, B7 andB8, which correspond to the orange diamonds in Fig. 7), thisis true even at higher rotation rates, where there is a greaterseparation in scales between the domain size and the horizon-tal scale of the near-onset convective motions. The failure of thelarge-scale dynamo in this case can probably be attributed to thefact that the large-scale vortex instability, which plays a crucialrole in initialising the dynamo, is inhibited in smaller domains(see also Guervilly et al. 2015). Certainly, the large-scale vor-tex (measured, for example, by the rms velocity) stops growinglong before the magnetic field becomes dynamically significant,which suggests that geometrical effects are limiting its growth inthese low aspect ratio cases. A strong suppression of the large-scale vortex will also limit the production of motions at thoseintermediate scales that are responsible for sustaining the large-scale dynamo.Figure 8 shows the cycle frequency, ω cyc = 2 π/τ cyc (where τ cyc is the cycle period), for the successful large-scale dynamos,as a function of f Ra . To ensure some degree of comparabilityacross the simulations, the cycle frequency has been normalisedby the ohmic decay time, τ η , in each case. This choice of nor-malisation was motivated by the comparability of τ η and τ cyc for the reference solution (and reflects the important role playedby diffusive processes in the underlying dynamo mechanism).Case E1 has been excluded from this plot because it has not beenevolved for long enough to produce an accurate determination ofthe cycle period, for reasons that are discussed below. Across the × × × f Ra ω c y c τ η Set ASet BSet CSet D
Fig. 8.
Cycle frequency, ω cyc , of large-scale dynamo simulations, plot-ted against f Ra . In each case, the cycle frequency is normalised by theohmic decay time, τ η . other simulations, a clear trend is observed when moving closerto convective onset: cycle frequencies tend to decrease with de-creasing f Ra . So, as we move closer to onset, the dynamo cyclesbecome longer (c.f. Calkins et al. 2016, who found a similar re-sult in a related quasi-geostrophic dynamo model). In fact, thistrend explains the uncertainty regarding whether or not C5 is alarge-scale dynamo. If it were a large-scale dynamo, it wouldhave an extremely long cycle period, and we were unable to runthis calculation for long enough to establish whether or not thiswas the case. At higher convective driving, modulational effects(as discussed in the next section) make it more difficult to deter-mine cycle frequencies in some cases. However, it is clear thatthere is a greater degree of variability in the cycle frequenciesat higher values of f Ra . For example, despite having rather simi-lar values of f Ra , the normalised cycle frequencies for cases D2and B3a differ by more than a factor of two. This suggests somedependence of the cycle frequency upon some of the other pa-rameters in the system besides f Ra (a possible candidate is theconvective Rossby number, Ro c , which is fixed in cases D1–3but is free to vary in cases B1–8). These differences aside, al-though there is not a single best fit curve, there is still a clear tendof increasing cycle frequencies with increasing levels of convec-tive driving. The shortest period dynamos have a cycle periodthat is about % of τ η . Even in these cases, the cycle periodsare still long enough that diffusive effects are probably playing asignificant role.Whilst these large-scale dynamos do exhibit weak departuresaway from symmetry about the mid-plane, there are no obvioussuggestions that the stratification of the fluid is playing a ma-jor role in the operation of the dynamo. Indeed, as illustratedin Fig. 9, our quasi-Boussinesq case (E1) produces much thesame dynamo solution as the reference case. Only a relativelyshort time-series is shown in this figure (in this quasi-Boussinesqregime, the time-step constraints that are associated with acous-tic modes make long calculations prohibitively expensive), but itis clear that the behaviour of this solution is qualitatively sim-ilar to that of the corresponding stages of the reference solu-tion, so it is very unlikely that this large-scale dynamo will notpersist. In such cases, it is clearly more efficient to consider aproper Boussinesq model – whilst not shown here, we have con- Article number, page 10 of 16. J. Bushby et al.: Large-scale dynamos in rapidly rotating convection
Fig. 9.
As Fig. 3, but here for the weakly stratified (quasi-Boussinesq)case, E1. Note that the time-series here is relatively short, because thenumerical time-step constraints make longer runs very difficult to carryout. firmed that a Boussinesq calculation (using the code describedby Cattaneo et al. 2003) can produce a similar large-scale dy-namo with the reference solution parameter values. Due to theBoussinesq symmetries, E x and E y are anti-symmetric about themid-plane, but the mean horizontal magnetic fields are againsymmetric. The only significant change due to the reductionin the level of stratification is in the cycle period. The Boussi-nesq run suggests a cycle period (normalised by the convectiveturnover time) that is approximately twice that of the referencesolution. Whilst it is still evolving, the quasi-Boussinesq caseE1 appears to be consistent with this, with a longer cycle pe- riod compared to that seen in the corresponding plots in Fig. 3.The ohmic decay time in case E1 is approximately turnovertimes, which is similar to that of the reference solution, so thislonger cycle period is not simply a function of normalisation.Having said that, as indicated by Fig. 8, small changes in theconvective driving can lead to large changes in the cycle periodin the vicinity of the reference solution, so this discrepancy be-tween the cycle periods is probably not that significant. If any-thing, it is more remarkable that this seems to be the only signif-icant effect on the dynamo that can be attributed to a reductionin the stratification of the layer. As noted in Table 1, there are a number of cases in which thelarge-scale dynamo activity is described as “intermittent”. Re-calling that the critical Rayleigh number for the onset of con-vection is Ra crit = 6 . · for this set of parameters, thesecases are considerably more supercritical (in terms of their con-vective driving) than the reference solution. For cases B3a andB3b, Ra ≈ . crit , whilst Ra ≈ . crit for case B4 and Ra ≈ . crit for case B5. In each case, the large-scale vor-tex instability leads to vigorous convective flows.Even in the case of the large-scale dynamo in the referencesolution, there is some evidence of modulation in the mean mag-netic field, with the peak horizontally averaged magnetic fieldvarying noticeably from one cycle to the next. Small increasesin the Rayleigh number amplify this effect. This is illustratedin Fig. 10, which shows the large-scale dynamo for case B2( Ra = 3 × ≈ . crit ). In this simulation, the depth-averaged mean squared horizontal magnetic field (at ∼ convective turnover times) peaks at approximately % of thesquared equipartition field strength. The weakest cycles peak atapproximately % of this equipartition value. Some fluctuationsare visible in both the U rms ( t ) and B rms ( t ) time-series, althoughthese are fairly modest, particularly in the case of U rms ( t ) .Moving to higher Rayleigh numbers, we see much more dra-matic modulation. Figure 11 shows the evolution of the large-scale dynamo for case B3a. Whilst the large-scale dynamo isgenerally weaker (in terms of its equipartition value) than similardynamos at lower Rayleigh number, it is still possible to iden-tify well-defined oscillations in the mean horizontal magneticfield. However, the modulation is now extremely pronounced,with “active phases” of large-scale cyclic behaviour punctuatedby periods of negligible large-scale activity (during which thegenerated magnetic field is predominantly small scale). At evenhigher Rayleigh numbers, this modulational pattern reverses,with bursts of large-scale activity surrounded by long periodsof relative inactivity (see Fig. 12). For case B3a, large-scale ac-tivity cycles can be seen over approximately two thirds of thetime-series (mostly active phases, but with recurrent periods ofinactivity). This behaviour becomes increasingly intermittent atlarger values of the Rayleigh number: B4 is active for about onethird of the time series, whereas B5 is active for approximatelya quarter of the time.Certainly in the case of B3a, this modulational behaviour israther reminiscent of that observed in the solar cycle, where ex-tended phases of reduced magnetic activity (such as the MaunderMinimum, see, e.g. Eddy 1976), often referred to as “grand min-ima”, have been a recurrent feature of the activity pattern over atleast the last 10,000 years (see, e.g. McCracken et al. 2013). Westress again that we do not claim to be modelling the solar dy-namo in this paper. Nevertheless there may be some similaritiesbetween the modulational mechanism in these simulations with Article number, page 11 of 16&Aproofs: manuscript no. paper_arxiv
Fig. 10.
Dynamo evolution for Case B2.
Top: time evolution of the rmsvelocity and magnetic field (the inset highlights the modulation dur-ing the nonlinear phase).
Middle: volume-averaged squared horizontalfields, ˜ B x /B (blue) and ˜ B y /B (red), as functions of time; the blackdotted line shows ( ˜ B x + ˜ B y ) /B . Bottom: time and z -dependence ofthe mean horizontal magnetic fields. that of the Sun. As can clearly be seen in Figs. 11 and 12, the rmsvelocity in these simulations tends to increase during grand min-ima. This corresponds to the re-emergence of the large-scale vor-tex instability. During these inactive phases, the magnetic fieldis no longer strong enough to suppress this large-scale flow (thiscan be seen in the time-series of the rms magnetic field), so itcan again grow. It then becomes temporarily suppressed againonce the magnetic field reaches dynamically significant levels.A number of authors have proposed mean-field models of the Fig. 11.
As Fig. 10, but here for Case B3a. solar dynamo in which long-term modulation arises as the resultof the magnetic field inhibiting the large-scale differential rota-tion (see, e.g. Tobias 1996b; Brooke et al. 2002; Bushby 2006).In these models, the occurrence of grand minima depended uponthere being a separation in timescales between viscous and mag-netic diffusion, so that the perturbations to the flow velocity re-laxed over a much longer time-scale than the period of oscilla-tion of the dynamo. There is a similar separation in time-scales inthese modulated convective dynamo simulations, with the large-scale vortex growing on a much longer time-scale than the cycleperiod of the large-scale dynamos. The exchange of energy be-tween the magnetic field and the flow can therefore give rise tothis modulational behaviour.
Article number, page 12 of 16. J. Bushby et al.: Large-scale dynamos in rapidly rotating convection
Fig. 12.
As Fig. 10, but here for Case B4.
4. Conclusions and discussion
One of the great challenges for dynamo theorists is to explainthe origin of large-scale astrophysical magnetic fields. As antic-ipated from theory, helically-forced turbulence in an electricallyconducting fluid can drive a large-scale dynamo in a Cartesiandomain (Brandenburg 2001). Such idealised flows can never berealised in nature, although the effects of rotation are believedto give rise to helical convective motions in many astrophysi-cal bodies, so similar large-scale dynamos might be expected insuch cases. However, the dynamo properties of rapidly rotatingconvection appear to be rather subtle. Near-onset rapidly rotat-ing convection can drive a large-scale dynamo in a Cartesian do-main (Childress & Soward 1972) but, unless there is an imposed shear (which tends to promote large-scale dynamo action, seee.g. Käpylä et al. 2008; Hughes & Proctor 2009), dynamos tendto be small-scale in the more turbulent regime that is relevant forastrophysics.Previous work (Käpylä et al. 2013; Masada & Sano 2014a,b)has demonstrated that it is possible to find a large-scale dynamoin moderately turbulent, rapidly rotating convection in a Carte-sian domain (without shear). These simulations are in a param-eter regime in which the large-scale vortex instability can oper-ate. Building on this previous work, we have carried out a de-tailed analysis of the underlying dynamo mechanism and havedemonstrated that the large-scale dynamo is driven by the com-ponents of the flow with a horizontal wavenumber in the range ≤ k/k ≤ . In particular, this confirms that the large-scalevortex itself (i.e. the k = k mode, which becomes stronglydamped) does not play a direct role in sustaining the dynamo.Having said that, the structure of the flow at the driving scales isstill influenced (to some extent, at least) by the tendency for thelarge-scale vortex instability to transfer energy to larger scales.So, even though the large-scale vortex itself is damped, the ef-fects of the underlying instability should not be discounted. Inparticular the initialisation of the large-scale dynamo, in thisparameter regime, may well depend upon there being an effi-cient large-scale vortex instability in the first place. As a pre-lude to a possible benchmarking exercise, we have verified thatthis dynamo can be reproduced by three different codes, all ofwhom produce quantitatively comparable solutions in the non-linear regime.Provided that the magnetic boundary conditions (which ap-pear to be very important in this context) allow magnetic fluxto escape from the domain, we have shown that this large-scaledynamo is robust to moderate changes in the parameters. Thedynamo appears to be largely insensitive to the level of stratifi-cation within the domain, although the calculations that we havecarried out do suggest that large-scale dynamos in more weakly-stratified domains may tend to have longer cycle periods. Mov-ing towards convective onset, the form of the dynamo remainslargely the same, although the cycle period of the large-scale dy-namo increases (very dramatically at the lowest Rayleigh num-bers). As the level of convective driving is increased, the cycleperiod decreases, and the cycle becomes increasingly modulated.At moderate Rayleigh numbers, this modulation is almost solar-like, with active phases punctuated by inactive grand minima. Athigher driving, the large-scale dynamo becomes increasingly in-termittent. The modulation is driven by an exchange of energybetween the magnetic field and the flow.Future investigations in this field will focus on even moreturbulent regimes at higher Rayleigh numbers, exploring broadranges of magnetic and thermal Prandtl numbers and rotationrates. The dynamo mechanism could also be analysed from theperspective of mean-field dynamo theory. We have seen thatthere is a positive correlation between the components of themean horizontal magnetic field and the mean horizontal elec-tromotive force, just as expected for an α dynamo located inthe northern hemisphere of a rotating star, but more could bedone to clarify this. It would also be worthwhile to investigatepossible connections between these moderately-turbulent large-scale dynamos and the near-onset Boussinesq dynamos. We haveshown that vertical field boundary conditions seem to promotethis large-scale dynamo, but this does not rule out the possibilitythat a similar dynamo could operate in a turbulent regime withthe horizontal field boundary conditions adopted in most Boussi-nesq studies. Finally, there is the intriguing question of the extentto which these Cartesian dynamos are of relevance to particular Article number, page 13 of 16&Aproofs: manuscript no. paper_arxiv astrophysical bodies. Is it possible to find analogous dynamosin a sphere or spherical shell? Obviously the cycle periods ofthese large-scale dynamos are too long to be of direct relevanceto solar-type dynamos, but we can speculate that there may besome possible application to planetary dynamos (where, at leastin the case of the Earth, we know that there are polarity reversalsover long time-scales).
Acknowledgements.
We would welcome any interest from the dynamo commu-nity in the benchmark exercise that is proposed in the appendix. This work hasbeen supported in part by the National Science Foundation (grant AST1615100),the Research Council of Norway under the FRINATEK (grant 231444), theSwedish Research Council (grant 621-2011-5076), the University of Coloradothrough its support of the George Ellery Hale visiting faculty appointment (AB),the Academy of Finland (grant 272157) to the ReSoLVE Centre of Excel-lence (PJK, MJK) and the UK Natural Environment Research Council (grantNE/M017893/1, CG). Some of the simulations were performed using the super-computers hosted by CSC – IT Center for Science Ltd. in Espoo, Finland, whoare administered by the Finnish Ministry of Education. Much of this work madeuse of the facilities of N8 HPC Centre of Excellence, provided and funded bythe N8 consortium and EPSRC (Grant No.EP/K000225/1). The Centre is co-ordinated by the Universities of Leeds and Manchester. PB would also like toacknowledge a useful discussion with Michael Proctor that helped to clarify thenature of the underlying dynamo mechanism. We would also like to thank thereferee for their helpful comments and suggestions.
References
Brandenburg, A. 2001, ApJ, 550, 824Brandenburg, A., Saar, S. H., & Turpin, C. R. 1998, ApJ, 498, L51Brooke, J., Moss, D., & Phillips, A. 2002, A&A, 395, 1013Bushby, P. J. 2006, MNRAS, 371, 772Calkins, M. A., Julien, K., Tobias, S. M., & Aurnou, J. M. 2015, J. Fluid Mech.,780, 143Calkins, M. A., Julien, K., Tobias, S. M., Aurnou, J. M., & Marti, P. 2016, Phys.Rev. E, 93, 023115Cattaneo, F., Emonet, T., & Weiss, N. 2003, ApJ, 588, 1183Cattaneo, F. & Hughes, D. W. 2006, J. Fluid Mech., 553, 401Chan, K. L. 2007, Astron. Nachr., 328, 1059Chan, K. L. & Mayr, H. G. 2013, Earth Plan. Sci. Lett., 371, 212Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stability (OxfordUniversity Press)Childress, S. & Soward, A. M. 1972, Phys. Rev. Lett., 29, 837Christensen, U. R., Aubert, J., Cardin, P., et al. 2001, Phys. Earth Plan. Int., 128,25Clarke, D. A. 1996, ApJ, 457, 291Colella, P. & Woodward, P. R. 1984, J. Comp. Phys., 54, 174Courvoisier, A., Hughes, D. W., & Proctor, M. R. E. 2009, Proceedings of theRoyal Society of London Series A, 466, 583Dobler, W., Stix, M., & Brandenburg, A. 2006, ApJ, 638, 336Eddy, J. A. 1976, Science, 192, 1189Evans, C. R. & Hawley, J. F. 1988, ApJ, 332, 659Fautrelle, Y. & Childress, S. 1982, Geophys. Astrophys. Fluid Dynam., 22, 235Favier, B. & Bushby, P. J. 2012, Journal of Fluid Mechanics, 690, 262Favier, B. & Bushby, P. J. 2013, J. Fluid Mech., 723, 529Favier, B. & Proctor, M. R. E. 2013, Phys. Rev. E, 88, 053011Favier, B., Silvers, L. J., & Proctor, M. R. E. 2014, Phys. Fluids, 26, 096605Guervilly, C., Hughes, D. W., & Jones, C. A. 2014, J. Fluid Mech., 758, 407Guervilly, C., Hughes, D. W., & Jones, C. A. 2015, Phys. Rev. E, 91, 041001Hughes, D. W. & Proctor, M. R. E. 2009, Phys. Rev. Lett., 102, 044501Jones, C. A. 2011, Ann. Rev. Fluid Mech., 43, 583Jones, C. A., Boronski, P., Brun, A. S., et al. 2011, Icarus, 216, 120Jones, C. A. & Roberts, P. H. 2000, J. Fluid Mech., 404, 311Julien, K., Rubio, A. M., Grooms, I., & Knobloch, E. 2012, Geophysical andAstrophysical Fluid Dynamics, 106, 392Käpylä, P. J., Korpi, M. J., & Brandenburg, A. 2008, A&A, 491, 353Käpylä, P. J., Korpi, M. J., & Brandenburg, A. 2009, ApJ, 697, 1153Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2013, Geophys. Astrophys.Fluid Dynam., 107, 244Käpylä, P. J., Mantere, M. J., & Hackman, T. 2011, ApJ, 742, 34Kunnen, R. P. J., Ostilla-Mónico, R., van der Poel, E. P., Verzicco, R., & Lohse,D. 2016, Journal of Fluid Mechanics, 799, 413MacGregor, K. B. & Charbonneau, P. 1997, ApJ, 486, 484Mantere, M. J., Käpylä, P. J., & Hackman, T. 2011, Astron. Nachr., 332, 876Marti, P., Schaeffer, N., Hollerbach, R., et al. 2014, Geophysical Journal Inter-national, 197, 119 Masada, Y. & Sano, T. 2014a, PASJ, 66, 2Masada, Y. & Sano, T. 2014b, ApJ, 794, L6Masada, Y. & Sano, T. 2016, ApJ, 822, L22Matthews, P. C., Proctor, M. R. E., & Weiss, N. O. 1995, J. Fluid Mech., 305,281McCracken, K., Beer, J., Steinhilber, F., & Abreu, J. 2013, Space Sci. Rev., 176,59Meneguzzi, M. & Pouquet, A. 1989, J. Fluid Mech., 205, 297Mizerski, K. A. & Tobias, S. M. 2013, Geophys. Astrophys. Fluid Dynam., 107,218Moffatt, H. K. 1978, Magnetic field generation in electrically conducting fluids(Cambridge University Press)Parker, E. N. 1993, ApJ, 408, 707Rotvig, J. & Jones, C. A. 2002, Phys. Rev. E, 66, 056308Sano, T., Inutsuka, S., & Miyama, S. M. 1999, in Astrophysics and Space ScienceLibrary, Vol. 240, Numerical Astrophysics, ed. S. M. Miyama, K. Tomisaka,& T. Hanawa, 383Soward, A. M. 1974, Royal Society of London Philosophical Transactions SeriesA, 275, 611St. Pierre, M. G. 1993, in Solar and Planetary Dynamos, ed. M. R. E. Proctor,P. C. Matthews, & A. M. Rucklidge, 295–302Stellmach, S. & Hansen, U. 2004, Phys. Rev. E, 70, 056312Stellmach, S., Lischper, M., Julien, K., et al. 2014, Phys. Rev. Lett., 113, 254501Stix, M. 2002, The sun: an introduction (Springer)Stone, J. M. & Norman, M. L. 1992, ApJS, 80, 791Tilgner, A. 2014, Phys. Rev. E, 90, 013004Tobias, S. M. 1996a, ApJ, 467, 870Tobias, S. M. 1996b, A&A, 307, L21Tobias, S. M., Cattaneo, F., & Brummell, N. H. 2011, ApJ, 728, 153van Leer, B. 1979, J. Comp. Phys., 32, 101Williamson, J. H. 1980, J. Comp. Phys., 35, 48
Article number, page 14 of 16. J. Bushby et al.: Large-scale dynamos in rapidly rotating convection
Appendix A: A possible dynamo benchmark
This appendix contains further details on the quantitative com-parison that has been carried out between the three codes thatare described in Sect. 2.4. The aim here is to assess whether ornot this large-scale dynamo solution could form the basis for acommunity-wide nonlinear dynamo benchmark. We choose topresent the details here to encourage other researchers with sim-ilar codes to try this calculation. If there is sufficient enthusiasmfrom the community to carry out a detailed benchmark study, wewould agree on a common set of initial conditions and carry out adetailed analysis of the early evolution of the system (e.g. com-paring the mean growth rate of U rms ( t ) during the large-scalevortex phase) as well as the final nonlinear state. Whilst this so-lution is relatively complicated, we would (at the very least) ex-pect different codes to produce quantitatively comparable statis-tics. Here, we focus upon the simpler questions of whether ornot the final nonlinear solution is robust to small changes in theinitial conditions, confirming that three independent codes canproduce quantitatively comparable nonlinear dynamos.Fixing λ = 2 , Pr = Pm = 1 and
Ta = 5 · , it wasfirst confirmed that all three codes agree on the critical Rayleighnumber, Ra crit = 6 . · , with exponentially growing (de-caying) solutions being obtained for values of Ra that are frac-tionally above (below) this value. Once this agreement was con-firmed, nonlinear dynamo runs were carried out for the refer-ence solution value of Ra = 2 . · (see the upper three rowsof Table 1). As has already been described, this reference so-lution exhibits highly non-trivial behaviour, in which the initialconvective instability is subject to a secondary hydrodynamicalinstability (corresponding to the large-scale vortex). The result-ing flow then drives a nonlinear large-scale dynamo in whichthe total magnetic and kinetic energies are in a state of near-equipartition. This solution is therefore an excellent test of all as-pects of any compressible Cartesian MHD code. Although theywere all evolved from the same initial polytropic state, randominitial perturbations were applied. Some quantitative differencesbetween the codes are therefore to be expected during the earlystages of evolution. If, however, it is possible to confirm that allof the codes converge upon the same nonlinear dynamo, then thiscomparison can be deemed to be successful.Figures A.1–A.4 show the evolution of this system for eachof the three codes. Figure A.1 shows the time evolution of therms velocity and magnetic field. The time- and z -dependenceof the mean horizontal magnetic field components is illustratedin Fig. A.2, whilst Fig. A.3 shows the corresponding vertically-averaged values. Finally, Fig. A.4 shows the mid-layer kineticand magnetic energy spectra at the end of the three runs. It isimmediately apparent that all three codes are producing quanti-tatively comparable nonlinear dynamos. There are similar pro-nounced peaks at k = 0 in the magnetic energy spectra, whilstthe kinetic energy spectra all have broad peaks at around k/k =7 – . The peak amplitudes of the large-scale horizontal magneticfield components show similar levels of agreement (as indicatedby Fig. A.3). Table A.1 shows the (time-averaged) values of u rms and b rms , the rms velocity and magnetic field (respectively), Table A.1.
Details on the benchmark simulations.
Case Code Grid u rms b rms τ cyc A1 1 ∼ ∼ ∼ Fig. A.1.
Root mean square velocity ( top ) and magnetic field ( bottom )for Code 1 (solid line), Code 2 (dashed line) and Code 3 (dash-dottedline). The inset in the upper plot shows the early time behaviour, high-lighting the different initial conditions that have been used. from the nonlinear phase of the dynamo. It is clear that thereis quantitative comparability for both of these quantities acrossthe three codes. Whilst there is some variability in terms of themeasured cycle periods ( , and convective turnovertimes), it should be stressed that there is some intrinsic variationin the cycle duration as the dynamo progresses. This is one ofthe main contributors to the error bars in Fig. 8, and the variationacross the three codes could simply reflect this variability. Thisvariability in the cycle period may complicate a future bench-marking exercise, but this can probably be addressed satisfacto-rily by running longer calculations to average the cycle periodover more cycles.We should also note some other differences between thethree cases. As shown by Fig. A.1, the initial large-scale vortexgrowth phase is longer in Code 3 (with a lower growth rate) thanit is in the other two codes. This can be probably attributed todifferences in the initial conditions: the initial rms velocity startsfrom a much lower level in the case of Code 3, which apparentlydelays the onset of the large-scale vortex instability. It is alsoworth noting that the seed magnetic field is weaker in Code 3than it is in the others, so it takes longer to grow to a level whereit is able to influence the flow. This dependence upon the strengthof the seed field clearly indicates that it is the Lorentz force thatis eventually suppressing the large-scale vortex instability ratherthan a geometrical constraint due to the finite box size. Furtherdifferences can be seen in Figs. A.2 and A.3, where the large-scale dynamo emerges at different times for the different codes. Article number, page 15 of 16&Aproofs: manuscript no. paper_arxiv
Fig. A.2.
Horizontally averaged horizontal magnetic fields (normalised by the equipartition field strength) as functions of z and time for Code 1( left ), Code 2 ( middle ) and Code 3 ( right ). Note that the time-axes for Codes 2 and 3 have been shifted by turnover times and turnovertimes (respectively) for ease of comparability with Code 1. Fig. A.3.
Volume-averaged horizontal (squared) magnetic fields, ˜ B x /B (blue) and ˜ B y /B (red), as functions of time for Code 1 ( left ), Code 2( middle ) and Code 3 ( right ). For each code, the black dotted line shows ( ˜ B x + ˜ B y ) /B as a function of time. Note that the sampling time for themean fields is less frequent in the right-hand plot than in the other two. Fig. A.4.