Issue 
A&A
Volume 564, April 2014



Article Number  A78  
Number of page(s)  13  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201322801  
Published online  10 April 2014 
Topology and field strength in spherical, anelastic dynamo simulations^{⋆}
MAG (ENS/IPGP), LRA/LERMA, Département de Physique, École Normale
Supérieure, 24 rue Lhomond,
75252
Paris Cedex 5,
France
email:
martin@schrinner.eu
Received:
4
October
2013
Accepted:
30
December
2013
Context. Numerical modelling of convection driven dynamos in the Boussinesq approximation revealed fundamental characteristics of the dynamogenerated magnetic fields and the fluid flow. Because these results were obtained for an incompressible fluid of constant density, their validity for gas planets and stars remains to be assessed. A common approach is to take some density stratification into account with the socalled anelastic approximation.
Aims. The validity of previous results obtained in the Boussinesq approximation is tested for anelastic models. We point out and explain specific differences between both types of models, in particular, with respect to the field geometry and the field strength, but we also compare scaling laws for the velocity amplitude, the magnetic dissipation time, and the convective heat flux.
Methods. Our investigation is based on a systematic parameter study of spherical dynamo models in the anelastic approximation. We make use of a recently developed numerical solver and provide results for the test cases of the anelastic dynamo benchmark.
Results. The dichotomy of dipolar and multipolar dynamos identified in Boussinesq simulations is also present in our sample of anelastic models. Dipolar models require that the typical length scale of convection is an order of magnitude larger than the Rossby radius. However, the distinction between both classes of models is somewhat less explicit than in previous studies. This is mainly due to two reasons: we found a number of models with a considerable equatorial dipole contribution and an intermediate overall dipole field strength. Furthermore, a large density stratification may hamper the generation of dipole dominated magnetic fields. Previously proposed scaling laws, such as those for the field strength, are similarly applicable to anelastic models. It is not clear, however, if this consistency necessarily implies similar dynamo processes in both settings.
Key words: convection / magnetohydrodynamics (MHD) / dynamo / methods: numerical / stars: magnetic field
Appendices B and C are available in electronic form at http://www.aanda.org
© ESO, 2014
1. Introduction
Magnetic fields of lowmass stars and planets are maintained by currents resulting from the motion of a conducting fluid (or gas) in their interiors. Because the magnetic field acts back on the flow via the Lorentz force, the hydrodynamic dynamo problem is intrinsically nonlinear. Moreover, a sufficiently complex flow and magnetic field geometry has to be assumed in order to enable dynamo action. Further complications result from tiny diffusivities, such as small kinematic viscosities, which introduce small dynamic length scales compared to stellar or planetary radii. Thus, selfconsistent simulations of natural dynamos are not only threedimensional, but a vast range of spatial and temporal scales has to be resolved. These difficulties prevented a direct numerical treatment of the dynamo problem for a significant time. Only for the past 20 years, increasing computer power made global, direct numerical simulations feasible, in particular for the geodynamo problem (e.g., Glatzmaier & Roberts 1995; Kageyama & Sato 1997; Kuang & Bloxham 1997; Christensen et al. 1998; Sarson et al. 1998; Katayama et al. 1999; Buffett 2000; Dormy et al. 2000); for an early cylindrical annulus model of the geodynamo see Busse (1975). In these simulations, an incompressible conducting fluid was considered, and the Boussinesq approximation was applied. Intensive and systematic parameter studies revealed fundamental properties of these models related to their field and flow topologies, their field strength, their velocity amplitudes, their advective heat transport, or their time dependence (e.g. Christensen et al. 1999; Grote et al. 2000; Kutzner & Christensen 2002; Busse & Simitev 2006; Christensen & Aubert 2006; Sreenivasan & Jones 2006; Busse & Simitev 2010; Hori et al. 2010; Landeau & Aubert 2011; Schrinner et al. 2012; Yadav et al. 2013a). The simplifying assumption of constant density in Boussinesq models, however, is probably not justified for gas planets or stars, in which the density typically varies over many scale heights. An alternative approach, which takes compressibility into account, is the socalled anelastic approximation (Ogura & Phillips 1962; Gough 1969; Gilman & Glatzmaier 1981). In anelastic models, the density varies with radius, but its time derivative is neglected in the continuity equation and the mass flux is solenoidal (e.g., Glatzmaier 1984; Braginsky & Roberts 1995; Lantz & Fan 1999; Miesch et al. 2000; Brun et al. 2004; Browning et al. 2004). Consequently, fast travelling sound waves are filtered out and, compared to fully compressible models, larger time steps in the discretisation scheme may be reached. In this article, we carry out a systematic parameter study of global dynamo simulations in the anelastic approximation guided by well known results of Boussinesq models. In this way, we intend to point out specific differences between anelastic and Boussinesq models and assess the validity of previous findings obtained in the Boussinesq approximation. A similar approach was followed by Gastine et al. (2012) and Yadav et al. (2013b). We compare our results with their findings and discuss some differences we obtained for our varied sample of models. The conditions for the generation of largescale, dipolar fields (Sect. 3) and the test of the fluxbased scalinglaw for the magnetic field strength (Sect. 4), originally proposed by Christensen & Aubert (2006) for Boussinesq models, are revisited . We argue that the typical length scale of convection relative to the Rossby radius is of crucial importance for the resulting field topology (see also Schrinner et al. 2012) and show that larger magnetic Prandtl numbers are required to obtain dipolar solutions with increasing density contrast. Furthermore, the fluxbased scaling laws derived for Boussinesq models seem to hold in the anelastic approximation as well. However, because of their general validity, the fluxbased scaling laws might not be appropriate to distinguish between different conditions for magnetic field generation.
The paper is organised as follows: Sect. 2 introduces the anelastic models considered here and the recently developed numerical solver; results for a numerical benchmark are given in Appendix A. In Sect. 3, we present new evidence of the existence of a class of dynamos dominated by an axial dipole and a class of models with a more variable magnetic field geometry. Various scaling laws originally derived for Boussinesq models are tested and discussed in Sect. 4, and we give some conclusions in Sect. 5.
2. Dynamo calculations
2.1. The anelastic approximation
Convection of a gas or a compressible fluid in the interior of planets and stars takes place on a vast range of spatial and temporal scales. Sound waves excited in convection zones, for example, have very short oscillation periods compared to the turnover time of convection or the magnetic diffusion time relevant for the generation of magnetic fields. Thus, extremely small timesteps would be required to resolve these waves in numerical dynamo models. To avoid this problem, simplifications of the governing equations are often applied. The anelastic approximation used in this study advantageously filters out sound waves (Ogura & Phillips 1962; Gough 1969; Gilman & Glatzmaier 1981); it is motivated by the idea that the superadiabatic temperature gradient driving convection in planetary or stellar convection zones is tiny. The thermodynamic variables are then decomposed into the sum of (close to adiabatic) reference values, denoted here by an overbar, and perturbations, denoted by a prime, $\mathit{\varrho}\mathrm{=}\overline{)\mathit{\varrho}}\mathrm{+}{\mathit{\varrho}}^{\mathrm{\prime}}\mathit{,}\u2001\mathit{T}\mathrm{=}\overline{)\mathit{T}}\mathrm{+}{\mathit{T}}^{\mathrm{\prime}}\mathit{,}\u2001\mathit{P}\mathrm{=}\overline{)\mathit{P}}\mathrm{+}{\mathit{P}}^{\mathrm{\prime}}\mathit{.}$(1)Subsequently, the anelastic equations result from the “thermodynamic linearization” around the reference state. It should be stressed that a number of different formulations of the anelastic problem can be found in the literature (e.g., Glatzmaier 1984; Braginsky & Roberts 1995; Lantz & Fan 1999; Miesch et al. 2000; Brun et al. 2004; Rogers & Glatzmaier 2005; Jones & Kuzanyan 2009; Alboussière & Ricard 2013). We follow here the approach introduced by Lantz & Fan (1999) and Braginsky & Roberts (1995), also known as LBRapproximation. They noticed that the only relevant thermodynamic variable in the equation of motion is the entropy, if the reference state is assumed to be close to adiabatic. Further advantages of the LBRequations over others are that they give a mathematically consistent, asymptotic limit of the full, general equations (Jones et al. 2011) and guarantee the conservation of energy (Brown et al. 2012). Moreover, the LBRequations were used to formulate anelastic dynamo benchmarks (Jones et al. 2011). The presentation of the equations given here follows the benchmark paper and Jones et al. (2009).
2.2. Basic assumptions
We consider a perfect, electrically conducting gas in a rotating spherical shell with an inner boundary at r = r_{i} and an outer boundary at r = r_{o}. The aspect ratio of the shell is then defined by χ = r_{i}/r_{o}. Convection in our simulations is driven by an imposed entropy difference, Δs, between the inner and the outer boundary. As discussed above, Δs is assumed to be small. This implies small convective velocities compared to the speed of sound. For consistency, we also require that the Alfvén velocity of the magnetic field is small. Moreover, the kinematic viscosity ν, the thermal diffusivity κ, and the magnetic diffusivity η are constants throughout in this paper. Following Jones et al. (2011), we represent the heat flux in our models in terms of the entropy gradient instead of the temperature gradient. This assumption relies on widespread ideas about turbulent mixing (Braginsky & Roberts 1995) but does not follow from first principles. Applying this simplification allows us to consider the entropy as the only relevant thermodynamic variable in the formulation of the anelastic problem.
2.3. The reference state
The reference state of our models is a solution of the hydrostatic equations for an adiabatic atmosphere. Moreover, the centrifugal acceleration is neglected and we assume that gravity varies radially, , with G being the gravitational constant and M the central mass of the star or the planet. This admits a polytropic solution for the reference atmosphere, $\overline{)\mathit{P}}\mathrm{=}{\mathit{P}}_{\mathrm{c}}\hspace{0.17em}{\mathit{w}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathit{,}\u2001\overline{)\mathit{\varrho}}\mathrm{=}{\mathit{\varrho}}_{\mathrm{c}}\hspace{0.17em}{\mathit{w}}^{\mathit{n}}\mathit{,}\u2001\overline{)\mathit{T}}\mathrm{=}{\mathit{T}}_{\mathrm{c}}\hspace{0.17em}\mathit{w,}\u2001\mathit{w}\mathrm{=}{\mathit{c}}_{\mathrm{0}}\mathrm{+}\frac{{\mathit{c}}_{\mathrm{1}}\mathit{d}}{\mathit{r}}\mathit{,}$(2)with the polytropic index n and d = r_{o} − r_{i}. We note that n defines the value of the adiabatic exponent γ, or the ratio of specific heats c_{p}/c_{v} via γ = (n + 1)/n. The values P_{c}, ϱ_{c}, and T_{c} are taken midway between the inner and the outer boundary and serve as units for the referencestate variables. Moreover, the constants c_{0} and c_{1} in (2) are defined as ${\mathit{c}}_{\mathrm{0}}\mathrm{=}\frac{\mathrm{2}{\mathit{w}}_{\mathrm{0}}\mathrm{}\mathit{\chi}\mathrm{}\mathrm{1}}{\mathrm{1}\mathrm{}\mathit{\chi}}\mathit{,}\u2001{\mathit{c}}_{\mathrm{1}}\mathrm{=}\frac{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\chi}\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{}{\mathit{w}}_{\mathrm{o}}\mathrm{)}}{\mathrm{(}\mathrm{1}\mathrm{}\mathit{\chi}{\mathrm{)}}^{\mathrm{2}}}\mathit{,}$(3)with ${\mathit{w}}_{\mathrm{0}}\mathrm{=}\frac{\mathit{\chi}\mathrm{+}\mathrm{1}}{\mathit{\chi}\mathrm{exp}\mathrm{(}{\mathit{N}}_{\mathit{\varrho}}\mathit{/}\mathit{n}\mathrm{)}\mathrm{+}\mathrm{1}}\mathit{,}\u2001{\mathit{w}}_{\mathrm{i}}\mathrm{=}\frac{\mathrm{1}\mathrm{+}\mathit{\chi}\mathrm{}{\mathit{w}}_{\mathrm{o}}}{\mathit{\chi}}\mathit{,}$(4)and N_{ϱ} = ln(ϱ_{i}/ϱ_{o}), where ϱ_{i} and ϱ_{o} denote the reference state density at the inner and outer boundary, respectively. We emphasise again that convection in our models is not driven by the reference state, or by the choice of a particular polytropic index n, but by an imposed entropy difference Δs between the boundaries.
2.4. The nondimensional equations
The use of nondimensional equations minimizes the number of free parameters and is a prerequisite for a systematic parameter study. We choose the shell width d as the fundamental length scale of our models, time is measured in units of d^{2}/η, and Δs is the unit of entropy. The magnetic field is then measured in units of $\sqrt{\mathrm{\Omega}{\mathit{\varrho}}_{\mathrm{c}}\mathit{\mu \eta}}$, where Ω is the rotation rate and μ the magnetic permeability. Finally, our dynamo models are solutions for the velocity v, the magnetic field B, and the entropy s of the following, nondimensional equations, $\begin{array}{ccc}\frac{\mathit{\partial}{v}}{\mathit{\partial t}}\mathrm{+}{v}\mathrm{\xb7}\mathrm{\nabla}{v}& \mathrm{=}& \mathit{P}\mathrm{m}\hspace{0.17em}[\mathrm{}\frac{\mathrm{1}}{\mathit{E}}\mathrm{\nabla}\frac{{\mathit{P}}^{\mathrm{\prime}}}{{\mathit{w}}^{\mathit{n}}}\mathrm{+}\frac{\mathit{P}\mathrm{m}}{\mathit{Pr}}\mathit{Ra}\frac{\mathit{s}}{{\mathit{r}}^{\mathrm{2}}}{r\u0302}\mathrm{}\frac{\mathrm{2}}{\mathit{E}}\hspace{0.17em}{z\u0302}\mathrm{\times}{v}\\ & & \mathrm{+}{{F}}_{\mathit{\nu}}\mathrm{+}\frac{\mathrm{1}}{\mathit{E}\hspace{0.17em}{\mathit{w}}^{\mathit{n}}}\mathrm{(}\mathrm{\nabla}\mathrm{\times}{B}\mathrm{)}\mathrm{\times}{B}]\mathit{,}\\ \frac{\mathit{\partial}{B}}{\mathit{\partial t}}& \mathrm{=}& \mathrm{\nabla}\mathrm{\times}\mathrm{(}{v}\mathrm{\times}{B}\mathrm{)}\mathrm{+}{\mathrm{\nabla}}^{\mathrm{2}}{B}\mathit{,}\\ \frac{\mathit{\partial s}}{\mathit{\partial t}}\mathrm{+}{v}\mathrm{\xb7}\mathrm{\nabla}\mathit{s}& \mathrm{=}& {\mathit{w}}^{\mathrm{}\mathit{n}\mathrm{}\mathrm{1}}\frac{\mathit{P}\mathrm{m}}{\mathit{Pr}}\mathrm{\nabla}\mathrm{\xb7}\mathrm{(}{\mathit{w}}^{\mathit{n}\mathrm{+}\mathrm{1}}\hspace{0.17em}\mathrm{\nabla}{\mathit{s}}^{\mathrm{)}}\\ & & \mathrm{+}\frac{\mathit{Di}}{\mathit{w}}\mathrm{[}{\mathit{E}}^{1}{\mathit{w}}^{\mathrm{}\mathit{n}}\mathrm{(}\mathrm{\nabla}\mathrm{\times}{B}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{Q}}_{\mathit{\nu}}\mathrm{]}\mathit{,}\\ \mathrm{\nabla}\mathrm{\xb7}\left({\mathit{w}}^{\mathit{n}}{v}\right)& \mathrm{=}& \mathrm{0}\mathit{,}\\ \mathrm{\nabla}\mathrm{\xb7}{B}& \mathrm{=}& \mathrm{0.}\end{array}$In (5), we used the viscous force F_{ν} = w^{−n}∇S with the rate of strain tensor ${\mathit{S}}_{\mathit{ij}}\mathrm{=}\mathrm{2}{\mathit{w}}^{\mathit{n}}\left({\mathit{e}}_{\mathit{ij}}\mathrm{}\frac{\mathrm{1}}{\mathrm{3}}{\mathit{\delta}}_{\mathit{ij}}\mathrm{\nabla}\mathrm{\xb7}{v}\right)\mathit{,}\u2001{\mathit{e}}_{\mathit{ij}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\mathit{\partial}{\mathit{v}}_{\mathrm{i}}}{\mathit{\partial}{\mathit{x}}_{\mathit{j}}}\mathrm{+}\frac{\mathit{\partial}{\mathit{v}}_{\mathit{j}}}{\mathit{\partial}{\mathit{x}}_{\mathrm{i}}}\right)\mathrm{\xb7}$(10)Moreover, the dissipation parameter Di and the viscous heating Q_{ν} in (7) are given by $\mathit{Di}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{1}}\mathit{Pr}}{\mathit{P}\mathrm{m}\mathit{Ra}}\mathit{,}$(11)and ${\mathit{Q}}_{\mathit{\nu}}\mathrm{=}\mathrm{2}\left[{\mathit{e}}_{\mathit{ij}}{\mathit{e}}_{\mathit{ij}}\mathrm{}\frac{\mathrm{1}}{\mathrm{3}}\mathrm{(}\mathrm{\nabla}\mathrm{\xb7}{v}{\mathrm{)}}^{\mathrm{2}}\right]\mathit{.}$(12)The system of Eqs. (5)−(9) is governed by a number of dimensionless parameters. These are the Rayleigh number Ra, the Ekman number E, the Prandtl number Pr, and the magnetic Prandtl number Pm. With the aspect ratio χ, the polytropic index n, and the number of density scale heights N_{ϱ} defining the reference state, our models are therefore fully determined by seven dimensionless parameters: $\begin{array}{ccc}\mathit{Ra}& \mathrm{=}& \frac{\mathit{GMd}\mathrm{\Delta}\mathit{s}}{\mathit{\nu \kappa}{\mathit{c}}_{\mathit{p}}}\mathit{,}\u2001\mathit{Pr}\mathrm{=}\frac{\mathit{\nu}}{\mathit{\kappa}}\mathit{,}\u2001\mathit{P}\mathrm{m}\mathrm{=}\frac{\mathit{\nu}}{\mathit{\eta}}\mathit{,}\u2001\mathit{E}\mathrm{=}\frac{\mathit{\nu}}{\mathrm{\Omega}{\mathit{d}}^{\mathrm{2}}}\hspace{0.17em}\mathit{,}\\ {\mathit{N}}_{\mathit{\varrho}}& \mathrm{=}& \mathrm{ln}\left(\frac{{\mathit{\varrho}}_{\mathrm{i}}}{{\mathit{\varrho}}_{\mathrm{o}}}\right)\mathit{,}\u2001\mathit{n,}\u2001\mathit{\chi}\mathrm{=}\frac{{\mathit{r}}_{\mathrm{i}}}{{\mathit{r}}_{\mathrm{o}}}\mathrm{\xb7}\end{array}$(13)Following Jones et al. (2009), we also considered a linearized form of Eqs. (5) and (7) to calculate some critical Rayleigh numbers for the onset of convection. These are listed in Table A.1.
2.5. Boundary conditions
The mechanical boundary conditions are impenetrable and stress free on both boundaries, ${\mathit{v}}_{\mathit{r}}\mathrm{=}\frac{\mathit{\partial}}{\mathit{\partial r}}\left(\frac{{\mathit{v}}_{\mathit{\theta}}}{\mathit{r}}\right)\mathrm{=}\frac{\mathit{\partial}}{\mathit{\partial r}}\left(\frac{{\mathit{v}}_{\mathit{\phi}}}{\mathit{r}}\right)\mathrm{=}\mathrm{0}\u2001\mathrm{on}\u2001\mathit{r}\mathrm{=}{\mathit{r}}_{\mathrm{i}}\u2001\mathrm{and}\u2001\mathit{r}\mathrm{=}{\mathit{r}}_{\mathrm{o}}\mathit{.}$(14)Furthermore, the magnetic field matches a potential field outside the fluid shell. The choice of these boundary conditions requires that the total angular momentum is conserved (Jones et al. 2011). Finally, the entropy is fixed on the inner and the outer boundary with $\mathit{s}\mathrm{=}\mathrm{\Delta}\mathit{s}\u2001\mathrm{on}\u2001\mathit{r}\mathrm{=}{\mathit{r}}_{\mathrm{i}}\u2001\mathrm{and}\u2001\mathit{s}\mathrm{=}\mathrm{0}\u2001\mathrm{on}\u2001\mathit{r}\mathrm{=}{\mathit{r}}_{\mathrm{o}}\mathit{.}$(15)
2.6. Output parameters
We use a number of nondimensional output parameters to characterize our numerical dynamo models. These are mostly based on the kinetic and magnetic energy densities, ${\mathit{E}}_{\mathrm{k}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\hspace{0.17em}\mathit{V}}{\mathrm{\int}}_{\mathit{V}}\hspace{0.17em}{\mathit{w}}^{\mathit{n}}{{v}}^{\mathrm{2}}\hspace{0.17em}\mathrm{d}\mathit{v}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{and}\u2001{\mathit{E}}_{\mathrm{m}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\hspace{0.17em}\mathit{V}}\frac{\mathit{P}\mathrm{m}}{\mathit{E}}{\mathrm{\int}}_{\mathit{V}}\hspace{0.17em}{{B}}^{\mathrm{2}}\hspace{0.17em}\mathrm{d}\mathit{v,}$(16)where the integrals are taken over the volume of the fluid shell V. A nondimensional measure for the velocity amplitude is then the magnetic Reynolds number, $\mathit{Rm}\mathrm{=}\sqrt{\mathrm{2}{\mathit{E}}_{\mathrm{k}}}$, or the Rossby number, $\mathit{Ro}\mathrm{=}\sqrt{\mathrm{2}{\mathit{E}}_{\mathrm{k}}}\mathit{E}\mathit{/}\mathit{P}\mathrm{m}$. To distinguish models with different field geometries, it turned out to be useful to introduce also a local Rossby number, Ro_{ℓ} = Ro_{c}ℓ_{c}/π. Here, ℓ_{c} stands for the mean harmonic degree of the velocity component v_{c} from which the mean zonal flow has been subtracted (Schrinner et al. 2012), ${\mathit{\ell}}_{\mathrm{c}}\mathrm{=}\sum _{\mathit{\ell}}\mathit{\ell}\frac{\u27e8{\mathit{w}}^{\mathit{n}}\hspace{0.17em}\mathrm{(}{{v}}_{\mathrm{c}}{\mathrm{)}}_{\mathit{\ell}}\mathrm{\xb7}\mathrm{(}{{v}}_{\mathrm{c}}{\mathrm{)}}_{\mathit{\ell}}\u27e9}{\mathrm{\u27e8}{\mathit{w}}^{\mathit{n}}\hspace{0.17em}{{v}}_{\mathrm{c}}\mathrm{\xb7}{{v}}_{\mathrm{c}}\mathrm{\u27e9}}\mathrm{\xb7}$(17)The brackets in (17) denote an average over time and radii. Also, Ro_{c} is adapted consistently and stands for the Rossby number based on the kinetic energy density without the contribution from the mean zonal flow. The definition of Ro_{ℓ} given here is different from Christensen & Aubert (2006), as it is not based on the total velocity and tries to avoid any dependence on the mean zonal flow.
The amplitude of the average magnetic field in our simulations is measured in terms of the Lorentz number, $\mathit{Lo}\mathrm{=}\sqrt{\mathrm{2}{\mathit{E}}_{\mathrm{m}}}\mathit{E}\mathit{/}\mathit{P}\mathrm{m}$, which was previously used to derive a power law for the field strength in Boussinesq simulations (Christensen & Aubert 2006). The topology of the field is characterized by the relative dipole field strength, f_{dip}, defined as the timeaverage ratio on the outer shell boundary of the dipole field strength to the total field strength.
The total amount of heat transported in and out of the fluid shell relative to the conductive heat flux is quantified by the Nusselt number, $\begin{array}{ccc}\mathit{N}{\mathit{u}}_{\mathrm{bot}}& \mathrm{=}& \mathrm{}\frac{\mathrm{(}\mathrm{exp}\mathrm{\left(}{\mathit{N}}_{\mathit{\varrho}}\mathrm{\right)}\mathrm{}\mathrm{1}\mathrm{)}{\mathit{w}}_{\mathrm{i}}{\mathit{r}}_{\mathrm{i}}^{\mathrm{2}}}{\mathrm{4}\mathit{\pi n}{\mathit{c}}_{\mathrm{1}}}{\mathrm{\int}}_{{\mathit{S}}_{\mathrm{i}}}\frac{\mathit{\partial s}}{\mathit{\partial r}}\mathrm{sin}\mathit{\theta}\hspace{0.17em}\mathrm{d}\mathit{\theta}\mathrm{d}\mathit{\phi ,}\\ \mathit{N}{\mathit{u}}_{\mathrm{top}}& \mathrm{=}& \mathrm{}\frac{\mathrm{(}\mathrm{1}\mathrm{}\mathrm{exp}\mathrm{(}\mathrm{}{\mathit{N}}_{\mathit{\varrho}}\mathrm{)}\mathrm{)}{\mathit{w}}_{\mathrm{o}}{\mathit{r}}_{\mathrm{o}}^{\mathrm{2}}}{\mathrm{4}\mathit{\pi n}{\mathit{c}}_{\mathrm{1}}}{\mathrm{\int}}_{{\mathit{S}}_{\mathrm{o}}}\frac{\mathit{\partial s}}{\mathit{\partial r}}\mathrm{sin}\mathit{\theta}\hspace{0.17em}\mathrm{d}\mathit{\theta}\mathrm{d}\mathit{\phi}\mathit{.}\end{array}$The integrals are taken here over the spherical surface at radius r_{i} and radius r_{o}, respectively. For a steady equilibrium state, Nu_{bot} and Nu_{top} are identical if time averaged. For later use, we also define a Nusselt number based on the advective heat flux alone, $\mathit{N}{\mathit{u}}^{\mathit{\star}}\mathrm{=}\mathrm{(}\mathit{N}{\mathit{u}}_{\mathrm{bot}}\mathrm{}\mathrm{1}\mathrm{)}\frac{\mathit{E}}{\mathit{Pr}}\mathit{,}$(20)and accordingly a quantity usually referred as the flux based Rayleigh number, $\mathit{R}{\mathit{a}}_{\mathit{Q}}\mathrm{=}\mathrm{(}\mathit{N}{\mathit{u}}_{\mathrm{bot}}\mathrm{}\mathrm{1}\mathrm{)}\hspace{0.17em}\frac{\mathit{Ra}\hspace{0.17em}{\mathit{E}}^{\mathrm{3}}}{{\mathit{r}}_{\mathrm{o}}^{\mathrm{2}}\hspace{0.17em}\mathit{P}{\mathit{r}}^{\mathrm{2}}}\mathrm{\xb7}$(21)The energy balance plays a crucial role in the classical derivation of scaling laws for the saturation level of the magnetic field. In particular, the fraction of ohmic to total dissipation, f_{ohm} = D/P, is introduced because it determines the available power used for the magnetic field generation. In an equilibrium state, the total dissipation equals the power released by buoyancy, $\mathit{P}\mathrm{=}\frac{\mathit{Ra}\hspace{0.17em}{\mathit{E}}^{\mathrm{3}}}{\mathit{Pr}\hspace{0.17em}\mathit{P}\mathrm{m}}{\mathrm{\int}}_{\mathit{V}}\hspace{0.17em}{\mathit{w}}^{\mathit{n}}\hspace{0.17em}{\mathit{v}}_{\mathit{r}}\hspace{0.17em}\mathit{s}\hspace{0.17em}\mathrm{d}\mathit{v,}$(22)and the ohmic dissipation is given by $\mathit{D}\mathrm{=}{\left(\frac{\mathit{E}}{\mathit{P}\mathrm{m}}\right)}^{\mathrm{2}}{\mathrm{\int}}_{\mathit{V}}\hspace{0.17em}\mathrm{(}\mathrm{\nabla}\mathrm{\times}{B}{\mathrm{)}}^{\mathrm{2}}\hspace{0.17em}\mathrm{d}\mathit{v}\mathit{.}$(23)In (22) and (23), we scaled P and D by ϱ_{c}Ω^{3}d^{5}.
2.7. Equation for a tracer field
For some of our models, we calculated the evolution of a magnetic tracer field simultaneously with (5)−(9), $\frac{\mathit{\partial}{{B}}_{\mathrm{Tr}}}{\mathit{\partial t}}\mathrm{=}\mathrm{\nabla}\mathrm{\times}\mathrm{(}{v}\mathrm{\times}{{B}}_{\mathrm{Tr}}\mathrm{)}\mathrm{+}{\mathrm{\nabla}}^{\mathrm{2}}{{B}}_{\mathrm{Tr}}\mathit{.}$(24)In the above equation, B_{Tr} is a passive vector field, which was advanced at each time step kinematically with the quenched velocity field, but B_{Tr} did not contribute to the Lorentz force (Cattaneo & Tobias 2009; Tilgner & Brandenburg 2008; Schrinner et al. 2010). Schrinner et al. (2010) found that B_{Tr} grows exponentially for multipolar dynamos but stays stable for models dominated by a dipole field. The stability of B_{Tr} also serves in this study to distinguish between different classes of dynamos (Schrinner et al. 2012).
2.8. Numerical implementation
The numerical solver used to compute solutions of Eqs. (5)−(9) is a recently developed anelastic version of PaRoDy (Dormy et al. 1998 and further developments). The code uses a poloidaltoroidal expansion and a pseudospectral spherical harmonic expansion. The numerical method is similar in these aspects to the one originally introduced in Glatzmaier (1984). The radial discretisation, however, is based on finite differences on a stretched grid (allowing for a parallelization by a radial domain decomposition). Moreover, the pressure term has been eliminated by taking twice the curl of the momentum equation. The anelastic benchmark results obtained with PaRoDy are presented in Appendix A.
3. Field topology
3.1. Dipolar and multipolar dynamos
Parameter studies for Boussinesq simulations revealed two distinct classes of dynamo models. They can be distinguished by their field geometry and are therefore referred to as “dipolar” and “multipolar” models (Kutzner & Christensen 2002; Christensen & Aubert 2006). The spatial variability of multipolar dynamos is a direct consequence of dynamo action in a turbulent environment and has to be expected. The class of dipolar dynamos, however, is more peculiar. Schrinner et al. (2011b) showed that these models are singlemode dynamos, that is, except for the fundamental mode, all more structured magnetic eigenmodes are highly damped. The singlemode property leads to further characteristic differences between both classes of dynamos apart from their different field geometries. Whereas the dipole axis is stable for models with a dominant axial dipole field, multipolar models show frequent polarity reversals (Kutzner & Christensen 2002) or oscillations (Goudard & Dormy 2008; Schrinner et al. 2012). A third fundamental difference between dipolar and multipolar models is related to their saturation mechanism. If a magnetic tracer field is advanced kinematically with the selfconsistent, quenched velocity field stemming from the full dynamo simulation, the tracer field grows exponentially for multipolar but not for dipolar models. Dipolar dynamos are “kinematically stable” and in this numerical experiment, the tracer field becomes aligned with the actual, selfconsistent magnetic field after some initial transition period (Schrinner et al. 2010). Finally, dipolar and multipolar dynamos follow slightly different scaling laws for the magnetic field (Christensen 2010; Schrinner et al. 2012; Yadav et al. 2013a). This aspect is further discussed in Sect. 4.1.
Christensen & Aubert (2006) proposed a criterion based on a local Rossby number to separate dipolar from multipolar dynamos. We adopt this criterion in a slightly altered form (Schrinner et al. 2012). It says that dipolar dynamos may be found if the typical length scale of convection, ℓ, is at least an order of magnitude larger than the Rossby radius, or, Ro_{ℓ} = v/(Ωℓ) < 0.12 (in which v is a typical rms velocity). Our criterion is different from Christensen & Aubert (2006), and entirely based on convection and not influenced by the mean zonal flow. This helped to generalize the Rossby number rule to models with different aspect ratios and mechanical boundary conditions (Schrinner et al. 2012). Moreover, our reinterpretation assumes that the magnetic field is generated only by convection and therefore explains why the Rossby number criterion is not applicable to models for which differential rotation plays an essential role.
Fig. 1 Relative dipole fieldstrength versus the local Rossby number for our sample of models. Filled symbols stand for dipolar, open symbols for multipolar dynamos. The symbol shape indicates the number of density scale heights: N_{ϱ} = 0.5: circle; N_{ϱ} = 1: upward triangle; N_{ϱ} = 1.5: downward triangle; N_{ϱ} = 2: diamond; N_{ϱ} = 2.5: square; N_{ϱ} = 3,3.5,4: star. A cross inscribed in some open symbols means that the field of these models exhibits a strong equatorial dipole component. 
Figure 1 shows the relative dipole field strength versus the local Rossby number for all anelastic models considered here. Gastine et al. (2012) presented a similar plot but with f_{dip} based on the magnetic energy density instead of the field strength. This leads to considerably lower values of f_{dip} for multipolar dynamos. As for Boussinesq simulations, only multipolar models are found for Ro_{ℓ}> 0.12 (Christensen & Aubert 2006), and the multipolar branch extends into the dipolar regime in the form of a bistable region, where both solutions are possible, depending on the initial conditions (Schrinner et al. 2012). However, in contrast to comparable parameter studies of Boussinesq models (Christensen & Aubert 2006; Schrinner et al. 2012), dipolar and multipolar dynamos are hardly distinguishable from each other in terms of their relative dipole field strength. Contrary to previous results, models with an intermediate dipolarity (f_{dip} ≈ 0.5) lead to a fairly smooth transition of f_{dip} in Fig. 1. These are in particular those models with a high equatorial dipole contribution denoted by a cross that is inscribed in the plotting symbol. Because the dipole field strength alone is not conclusive to classify our models in Fig. 1, their timedependence, their kinematic stability, and their scaling behaviour (see Sect. 4.1) were additionally considered to assign them to one of both classes.
As in the case of Boussinesq simulations, only multipolar models were found to exhibit polarity reversals or oscillatory dynamo solutions. An example of a coherent dynamo wave for model3m (N_{ϱ} = 3) is given in Fig. 2. The period of these oscillatory dynamo modes and the poleward propagation direction of the resulting wave can be surprisingly well explained by Parker’s plane layer formalism (Parker 1955; Busse & Simitev 2006; Goudard & Dormy 2008; Schrinner et al. 2011a; Gastine et al. 2012). However, the recent claim that dynamo waves could migrate towards the equator if there is a considerable density stratification (Käpylä et al. 2013) was not confirmed by our simulations.
Fig. 2 Contour plot of the azimuthally averaged radial magnetic field of model3m versus time and colatitude. The contour plot was normalised by the maximum absolute value at each time step. The greyscale coding ranges from −1, white, to +1, black. 
Moreover, we tested 13 arbitrarily chosen models (see the caption of Table B.1) for kinematic stability and found the dipolar models to be kinematically stable, whereas all multipolar models considered exhibited at least periods of instability. Figure 3 shows as an example the evolution of the kinematically advanced tracer field for model2m and model54d. For the first, the tracer field grows exponentially but it stays stable for the latter although it has been permanently perturbed during the simulation.
Fig. 3 Evolution of the energy of the tracer field normalised by the energy of the actual magnetic field for model2m (dashed line) and model54d (solid line). 
Fig. 4 Relative dipole field strength versus Ro_{ℓ} for a sequence of models with E = 3 × 10^{4}, Ra = 4Ra_{c}, Pm = 3, and Pr = 1. The meaning of the symbols is defined in the caption of Fig. 1. 
A transition from the dipolar to the multipolar regime can be triggered by a decrease in the rotation rate or the dynamical length scale (possibly associated with a change in the aspect ratio), or an increase in the velocity amplitude. These three quantities influence the local Rossby number directly. In Fig. 4, we show that a transition towards the multipolar regime may also be forced by increasing N_{ϱ}. A higher density stratification with all the other parameters fixed causes smaller length scales and larger velocity amplitudes. This leads to an increase of Ro_{ℓ} and to a decrease of f_{dip} at Ro_{ℓ} ≈ 0.12 in Fig. 4.
Fig. 5 Contour plot (equatorial cut) of the radial magnetic field of model2m at a given time. 
3.2. Equatorial dipole
An example of a model strongly influenced by an equatorial dipole mode is presented in Fig. 5. A strong mean zonal flow often present in these models seems to be in conflict with the generation of nonaxisymmetric fields. Figure 6 demonstrates that the strong equatorial dipole mode of model5m is indeed maintained and rebuilt by the columnar convection and damped by the differential rotation. In Fig. 6 the mean zonal kinetic energy normalised by an arbitrary value (dotted line) and the ratio of the axisymmetric magnetic energy to the total magnetic energy (solid line) are displayed. The action of the mean zonal flow, or more precisely the differential rotation, tends to damp nonaxisymmetric components of the magnetic field. Thus, a burst of the mean zonal kinetic energy is followed by a maximum of the axisymmetric and a dip in the nonaxisymmetric magnetic energy. Subsequently, the mean zonal flow is quenched by the axisymmetric field, the axisymmetric field decays and the nonaxisymmetric field is rebuilt. The interaction between the mean zonal flow and the magnetic field observed in this model is still fairly weak, although the mean zonal flow contributes already 58% to the total kinetic energy. Therefore, the magnetic field of model5m stays on average highly nonaxisymmetric. We note that this is very different from the Sun, for instance, where probably an even more efficient differential rotation causes a predominantly axisymmetric largescale magnetic field (Charbonneau 2010), but also nonaxisymmetric stellar magnetic fields were reported (Donati & Landstreet 2009).
Fig. 6 Dotted line: axisymmetric kinetic energy of model5m normalised by an arbitrary value. Solid line: ratio of axisymmetric to total magnetic energy. 
3.3. Discussion
The fundamental cause of the high dipolarity of dynamo models in the low Rossby number regime is an outstanding question. Schrinner et al. (2012) argued that cylindrical convection in a spherical fluid domain leads to a characteristic pattern of the axisymmetric toroidal field for Boussinesq models, which eventually results in the clear preference of only one, dipolar eigenmode. The argument relies on the idea that a line of fluid elements moving towards the outer spherical boundary has to shorten and causes a converging flow towards the equatorial plane. The toroidal field is then advected and markedly shaped by this flow component (see also Olson et al. 1999). This advection process could be rigorously identified and quantified as a strong γeffect in a corresponding meanfield description (Schrinner et al. 2007, 2012). In addition, the recent finding that the dichotomy of dipolar and multipolar dynamos seems to be absent in convective dynamo simulations in Cartesian geometry (Tilgner 2012) is consistent with this argument and points again to the significance of the underlying symmetry constraints.
What has been said above about Boussinesq models largely applies to anelastic models, too. However, geometrical constraints are somewhat relaxed for a compressible fluid. Therefore, compressibility might damp the advection of the mean toroidal field towards the equatorial plane (γeffect), and we hypothesize that this results in at least two specific differences.
First, depending on the density contrast applied, it is more difficult to obtain dipolar solutions for anelastic models than for Boussinesq ones, even if Ro_{ℓ} < 0.12. However, unlike Gastine et al. (2012), we did not find that dipolar solutions become impossible if N_{ϱ} exceeds a certain threshold. Instead, we observe that for a given N_{ϱ}, Ekman and Prandtl number, there seems to exist a critical magnetic Prandtl number for dipolar dynamos. For E = 10^{4}, Pr = 1, and N_{ϱ} ≥ 1.5, we found Pm_{crit} = 2N_{ϱ} − 2, as apparent from Fig. 7. We emphasize again that the results of Fig. 7 depend of course on E and Pr; the data of our numerical study indicate that decreasing E and increasing Pr is favorable to dipolar dynamo models.
Second, magnetic field configurations dominated by an equatorial dipole seem to be more easily realized in anelastic than in Boussinesq simulations. For the latter, only a few examples under very specific conditions were reported (Aubert & Wicht 2004; Gissinger et al. 2012). The preference of nonaxisymmetric modes is well known from dynamo models based on columnar convection (e.g. Ruediger 1980; Tilgner 1997), it is also the case of the Karlsruhe dynamo experiment (Müller & Stieglitz 2002). This agrees with our reasoning on the importance of the γeffect in the axial dipole generation mechanisms (see also Schrinner et al. 2012). Indeed, the γeffect vanishes in the above examples, as the geometrical constraints are relaxed.
Fig. 7 Magnetic Prandtl number versus N_{ϱ} for models with E = 10^{4} and Pr = 1 and variable Rayleigh numbers. Filled circles stand for parameters for which dipolar solutions were obtained. 
4. Scaling laws
Because of computational limitations, very small length scales and time scales associated with extreme parameter values that are relevant for planets and stars cannot be resolved in global direct numerical dynamo simulations. Therefore, numerical models are in general not directly comparable to planetary or stellar dynamos. Instead, scaling laws, in particular for the field strength, have been derived from theory and simulations and then extrapolated to realistic parameter regimes (see Christensen 2010, and references therein).
Subsequently, their predictions may be compared with planetary or stellar magneticfield data obtained from observations (Christensen et al. 2009; Christensen 2010; Davidson 2013). By this consistency test, scaling laws may provide some evidence about the reliability of numerical dynamo models.
Moreover, different scaling laws typically represent different force balances or dynamo mechanisms, and their investigation might enable us to better distinguish between different types of dynamo models. It is this second aspect in particular, which is of interest in the following. We adopt here the approach by Christensen & Aubert (2006) and derive scaling laws for the field strength, the velocity, the magnetic dissipation time, and the convective heat transport and compare them with previous results from Boussinesq simulations. A similar study was recently published by Yadav et al. (2013b) based on a somewhat different sample of models. Similarities and differences with their findings are discussed.
Most of the proposed scaling laws are independent of diffusivities, which are thought to be negligible under astrophysical conditions (Christensen 2010). However, present global dynamo simulations run in parameter regimes where diffusivities still influence the overall dynamics and weak dependencies on the magnetic Prandtl number seem to persist in purely empirically derived scalings (Christensen & Tilgner 2004; Christensen & Aubert 2006; Christensen 2010; Yadav et al. 2013a; Stelzer & Jackson 2013). In this study we do not attempt to resolve this secondary dependence on Pm because the magnetic Prandtl number varies only between 1 and 5 in our sample of models.
4.1. Magnetic field scaling
The magnetic field strength measured in terms of the Lorentz number scales with the available energy flux to the power of approximately 1/3. For the dipolar dynamos of our sample we find $\frac{\mathit{Lo}}{{\mathit{f}}_{\mathrm{ohm}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathrm{\simeq}\mathrm{1.58}\hspace{0.17em}\mathit{R}{\mathit{a}}_{\mathit{Q}}^{\mathrm{0.35}}\mathit{,}$(25)and for the multipolar models, $\frac{\mathit{Lo}}{{\mathit{f}}_{\mathrm{ohm}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathrm{\simeq}\mathrm{1.19}\hspace{0.17em}\mathit{R}{\mathit{a}}_{\mathit{Q}}^{\mathrm{0.34}}\mathit{.}$(26)Except for somewhat larger exponential prefactors, this is in good agreement with previous results from Boussinesq simulations (Christensen 2010; Schrinner et al. 2012; Yadav et al. 2013a) and very similar to the magnetic field scaling given by Yadav et al. (2013b). Unlike Yadav et al. (2013b), we note, however, that we scale the Lorentz number with the fluxbased Rayleigh number Ra_{Q} and not directly with the power released by buoyancy forces. Of course, both should be closely related to each other. The same remark applies for the velocity scaling discussed below.
Models on the multipolar branch exhibit lower field strengths compared to their dipolar counterparts. This is not only apparent by the smaller prefactor in the multipolar scaling, but also the dynamo efficiency f_{ohm} for multipolar models is systematically lower than for the corresponding dipolar ones. The latter indicates that the bistable behaviour for models at Ro_{ℓ} ≤ 0.12 is caused by different dynamo mechanisms. This was already seen in Boussinesq simulations (Schrinner et al. 2012) and later confirmed by Gastine et al. (2012) for anelastic models.
Apart from a few exceptions, the shift between the two scalings in Fig. 8 may serve to separate dipolar from multipolar dynamos. In agreement with Yadav et al. (2013b), we obtained several models with dipole field strengths up to f_{dip} ≈ 0.5, which, nevertheless, clearly follow the multipolar scaling and belong to the multipolar class of dynamos.
Fig. 8 Lorentz number compensated by f_{ohm} versus the fluxbased Rayleigh number for our sample of models. Filled symbols correspond to dipolar models, open symbols are multipolar models and the symbol shape indicates N_{ρ} as explained in the caption of Fig. 1. 
4.2. Velocity scaling
There is an ongoing discussion about the velocity scaling in dynamo models (Christensen 2010; Davidson 2013; Yadav et al. 2013b). It is probably not surprising that the velocity measured in terms of the Rossby number scales with the flux based Rayleigh number, but the correct exponent and its theoretical justification is debated. The lower bound is set by the assumption of a balance between inertia and buoyancy forces (mixing length balance), which leads to an exponent of 1/3 (Christensen 2010). If, however, the predominant force balance is assumed to be between the Lorentz force, the buoyancy, and the Coriolis force (MACbalance) the exponent is closer to 1 /2 (Christensen 2010; Davidson 2013). As most previous studies (Christensen & Aubert 2006; Christensen 2010; Yadav et al. 2013a; Stelzer & Jackson 2013; Yadav et al. 2013b), we obtained an exponent in between these two values for our sample of models, $\mathit{Ro}\mathrm{=}\mathrm{1.66}\hspace{0.17em}\mathit{R}{\mathit{a}}_{\mathit{Q}}^{\mathrm{0.42}}\mathit{.}$(27)The scatter in Fig. 9 is considerable, but the standard error is of the same order as for Boussinesq models with stressfree mechanical boundary conditions (Yadav et al. 2013a). Compressible effects do not seem to deteriorate the scaling.
Fig. 9 Rossby number versus the fluxbased Rayleigh number for our sample of models. 
However, as in Yadav et al. (2013b), we are not able to distinguish between dipolar and multipolar models in our velocity scaling, which is contrary to what has been previously reported by Yadav et al. (2013a) for Boussinesq models.
Scaling laws for anelastic and Boussinesq models.
4.3. Scaling of Ohmic dissipation time
The scaling of magnetic dissipation time, ${\mathit{\tau}}_{\mathrm{diss}}\mathrm{=}{\mathit{E}}_{\mathrm{M}}\mathit{/}\mathit{D}\mathrm{=}{\mathit{\ell}}_{\mathit{B}}^{\mathrm{2}}\mathit{/}\mathit{\eta ,}$(28)is used to evaluate the characteristic length scale ℓ_{B} of the magnetic field. Christensen & Tilgner (2004) originally identified a linear dependence of τ_{diss} on the inverse Rossby number provided that time is measured in units of Ω^{1}. Their finding was supported by dipoledominated Boussinesq models with noslip mechanical boundary conditions and the evaluation of the Ohmic dissipation time in the Karlsruhe dynamo experiment. The best fit for our data points in Fig. 10, however, gives an exponent with a significantly lower absolute value, ${\mathit{\tau}}_{\mathrm{diss}}\mathrm{=}\mathrm{0.75}\hspace{0.17em}\mathit{R}{\mathit{o}}^{0.76}\mathit{.}$(29)An almost identical result was found by Yadav et al. (2013b) from their somewhat more diverse and scattered data set. Apparently, the application of stressfree boundary conditions and possibly the compressible effects flatten the slope of τ_{diss} as a function of the Rossby number. Moreover, it would seem plausible that τ_{diss} followed different scaling relations for dipolar and multipolar models. Indeed, the dissipation time for bistable pairs is systematically larger for dipolar than for multipolar models. However, separate least square fits for all dipolar and all multipolar models of our sample lead to very similar results.
Fig. 10 Ohmic dissipation time versus Rossby number for all models considered in this study. 
4.4. Nusselt number scaling
The convective heat transport in dynamo models is very sensitive to rotation and depends on the magnetic field, boundary conditions, or the geometry of the fluid domain to a much lower degree (Christensen 2002; Christensen & Aubert 2006; Aurnou 2007; Schmitz & Tilgner 2009; Busse & Simitev 2011; Gastine & Wicht 2012; Yadav et al. 2013a; Stelzer & Jackson 2013). The power law for the Nusselt number inferred from Fig. 11, $\mathit{N}{\mathit{u}}^{\mathit{\star}}\mathrm{=}\mathrm{0.25}\hspace{0.17em}\mathit{R}{\mathit{a}}_{\mathit{Q}}^{\mathrm{0.59}}\mathit{,}$(30)is consistent with previous results and also confirms this finding for anelastic dynamo models; the exponent of 0.59 is very close to the value of 5/9 established by the above mentioned references. However, the scaling is somewhat more scattered than for Boussinesq models (Yadav et al. 2013a). We excluded in a test all models for which convection is only marginaly above the onset (Nu^{⋆} < 2), but this reselection of models did not improve the quality of the fit.
Fig. 11 Nusselt number versus the fluxbased Rayleigh number for our sample of models. 
4.5. Discussion
In an overall view, the scaling relations for Boussinesq and anelastic models are very similar (see Table 1). Beyond that, there is no obvious effect of compressibility on the scaling results and they might be even considered as consistent irrespective of the density stratification of the underlying models (Yadav et al. 2013b). However, the reason for the good agreement could be that the fluxbased scaling laws are insensitive to different physical conditions. Using the example of the magnetic field scaling, we argue in the following that differences in the dynamo processes might not be visible in the scaling relation, and some caution is needed in generalizing results from Boussinesq simulations.
If the magnetic energy density follows a simple power law in terms of the convective energy flux, an exponent of 2 /3 is already required for dimensional reasons (e.g., Christensen 2010). Moreover, the fluxbased scaling law for the magnetic field is composed of the scalings for the velocity and the magnetic dissipation time. By definition, we have E_{M} ~ f_{ohm}τ_{diss}P, and with Ro ~ P^{α} and τ_{diss} ~ Ro^{β}, we find E_{M} ~ f_{ohm}PP^{αβ}. Dimensional arguments require αβ = −1/3, which establishes relations (25) and (26). Whereas the exponent in the fluxbased scaling law for the magnetic field is fix, α and β are to some extend variable and may change according to the specific physical conditions. This reflects the outcome of more and more extended parameter studies: The exponent of 1/3 in the magneticfield scaling is reliably reproduced, but the values for α and β seem to be less certain and are under debate.
In addition, scaling relations (25) and (26) require that the field strength, measured by Lo, is compensated by the square root of f_{ohm} (interpreted as dynamo efficiency in Schrinner 2013). However, f_{ohm} is probably a complicated function of several control parameters and might depend strongly on the specific physical conditions. The often made assumption that f_{ohm} → 1 for Pm ≪ 1 (e.g. Davidson 2013) is probably too simple. For example, Schrinner (2013) recently demonstrated that f_{ohm} in dynamo models might depend strongly on the rotation rate. The dynamo efficiency dropped by two orders of magnitude as the rotation rate of these models was decreased. A further counterexample could be the solar dynamo. Independent estimates result in f_{ohm} ~ O(10^{3}), (Schrinner 2013; Rempel 2006) although the magnetic Prandtl number is thought to be much smaller than one in the solar interior.^{1} In other words, the flux based scaling laws probably do not distinguish between different types of dynamos because differences in the field strength are absorbed by changes in f_{ohm}.
5. Conclusions
Our study revealed a number of similarities between Boussinesq and anelastic dynamo models. The dichotomy between dipolar and multipolar models seems to extend to anelastic models, and the fluxbased scaling laws originally proposed for Boussineq models appear to hold similarly for models in the anelastic approximation. Thus, large scale, dipolar magnetic fields for both types of models can only be produced if rotation is important (as measured by the local Rossby number), and the magnetic field strength is directly related to the energy flux via (25) and (26) (see Fig. 8).
However, we also pointed out some significant differences between Boussinesq and anelastic dynamo simulations. Magnetic field configurations with a significant equatorial dipole contribution are less typical for Boussinesq than for anelastic models. Moreover, a large density stratification in anelastic models may inhibit the generation of magnetic fields dominated by an axial dipole. The above claimed consistency of the scalings for Boussinesq and anelastic simulations partly relies on the very general formulation of the fluxbased scaling laws and does not necessarily imply similar dynamo processes. We also stress that the assumption of a radially varying conductivity may introduce additional effects, which were not examined here. Whereas Yadav et al. (2013b) obtained very similar scaling laws for models with variable conductivities, Duarte et al. (2013) reported that the field topology of some models depends on the radial conductivity profile. A meanfield analysis (Schrinner et al. 2007; Schrinner 2011) of numerical dynamo models in the anelastic approximation might give more detailed insight in relevant dynamo processes and is envisaged for a future study.
Online material
Appendix A: Critical Rayleigh numbers for the onset of convection
Overview of the critical Rayleigh numbers and the corresponding critical azimuthal wavenumbers for the onset of convection determined, as explained in Jones et al. (2009).
Appendix B: Numerical models
Overview of the simulations carried out, ordered with respect to their local Rossby number.
In Schrinner (2013), a wrong mean solar density has been used to estimate f_{ohm} which lead to f_{ohm} ≈ 10^{2}. We correct this error here.
Acknowledgments
M.S. is grateful for financial support from the DFG fellowship SCHR 1299/11. Computations were performed at CINES, CEMAG, and GWDG computing centres. Moreover, this work was granted access to the HPC resources of MesoPSL financed by the Région Ile de France and the project Equip@Meso (reference ANR10EQPX2901) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche.
References
 Alboussière, T., & Ricard, Y. 2013, J. Fluid Mech., 725, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Aubert, J., & Wicht, J. 2004, Earth Planetary Sci. Lett., 221, 409 [Google Scholar]
 Aurnou, J. M. 2007, Geophys. Astrophys. Fluid Dyn., 101, 327 [NASA ADS] [CrossRef] [Google Scholar]
 Braginsky, S. I., & Roberts, P. H. 1995, Geophysical and Astrophysical Fluid Dynamics, 79, 1 [Google Scholar]
 Brown, B. P., Vasil, G. M., & Zweibel, E. G. 2012, ApJ, 756, 109 [Google Scholar]
 Browning, M. K., Brun, A. S., & Toomre, J. 2004, ApJ, 601, 512 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., Miesch, M. S., & Toomre, J. 2004, ApJ, 614, 1073 [NASA ADS] [CrossRef] [Google Scholar]
 Buffett, B. A. 2000, Science, 288, 2007 [NASA ADS] [CrossRef] [Google Scholar]
 Busse, F. H. 1975, Geophys. J. Int., 42, 437 [Google Scholar]
 Busse, F. H., & Simitev, R. D. 2006, Geophys. Astrophys. Fluid Dyn., 100, 341 [NASA ADS] [CrossRef] [Google Scholar]
 Busse, F. H., & Simitev, R. D. 2010, in Turbulence in the Atmosphere and Oceans, IUTAM Bookseries (Springer), 28, 181 [Google Scholar]
 Busse, F., & Simitev, R. 2011, Geophys. Astrophys. Fluid Dyn., 105, 234 [NASA ADS] [CrossRef] [Google Scholar]
 Cattaneo, F., & Tobias, S. M. 2009, J. Fluid Mech., 621, 205 [NASA ADS] [CrossRef] [Google Scholar]
 Charbonneau, P. 2010, Liv. Rev. Sol. Phys., 7 [Google Scholar]
 Christensen, U. R. 2002, J. Fluid Mech., 470, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Christensen, U. R. 2010, Space Sci. Rev., 152, 565 [NASA ADS] [CrossRef] [Google Scholar]
 Christensen, U. R., & Aubert, J. 2006, Geophy. J. Int., 166, 97 [Google Scholar]
 Christensen, U. R., & Tilgner, A. 2004, Nature, 429, 169 [NASA ADS] [CrossRef] [Google Scholar]
 Christensen, U., Olson, P., & Glatzmaier, G. A. 1998, Geophys. Res. Lett., 25, 1565 [NASA ADS] [CrossRef] [Google Scholar]
 Christensen, U., Olson, P., & Glatzmaier, G. A. 1999, Geophys. J. Int., 138, 393 [NASA ADS] [CrossRef] [Google Scholar]
 Christensen, U. R., Holzwarth, V., & Reiners, A. 2009, Nature, 457, 167 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Davidson, P. A. 2013, Geophys. J. Int., 195, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Donati, J.F., & Landstreet, J. D. 2009, ARA&A, 47, 333 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Dormy, E., Cardin, P., & Jault, D. 1998, Earth Planet. Sci. Lett., 160, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Dormy, E., Valet, J.P., & Courtillot, V. 2000, Geochem. Geophys. Geosyst., 1, 1037 [NASA ADS] [CrossRef] [Google Scholar]
 Duarte, L. D. V., Gastine, T., & Wicht, J. 2013, Phys. Earth Planet. Inter., 222, 22 [Google Scholar]
 Gastine, T., & Wicht, J. 2012, Icarus, 219, 428 [NASA ADS] [CrossRef] [Google Scholar]
 Gastine, T., Duarte, L., & Wicht, J. 2012, A&A, 546, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gilman, P. A., & Glatzmaier, G. A. 1981, ApJS, 45, 335 [NASA ADS] [CrossRef] [Google Scholar]
 Gissinger, C., Petitdemange, L., Schrinner, M., & Dormy, E. 2012, Phys. Rev. Lett., 108, 234501 [NASA ADS] [CrossRef] [Google Scholar]
 Glatzmaier, G. A. 1984, J. Comput. Phys., 55, 461 [Google Scholar]
 Glatzmaier, G. A., & Roberts, P. H. 1995, Nature, 377, 203 [NASA ADS] [CrossRef] [Google Scholar]
 Goudard, L., & Dormy, E. 2008, Europhys. Lett., 83, 59001 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gough, D. O. 1969, J. Atmos. Sci., 26, 448 [NASA ADS] [CrossRef] [Google Scholar]
 Grote, E., Busse, F. H., & Tilgner, A. 2000, Phys. Earth Planet. Inter., 117, 259 [NASA ADS] [CrossRef] [Google Scholar]
 Hori, K., Wicht, J., & Christensen, U. R. 2010, Phys. Earth Planet. Inter., 182, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, C. A., & Kuzanyan, K. M. 2009, Icarus, 204, 227 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, C. A., Kuzanyan, K. M., & Mitchell, R. H. 2009, J. Fluid Mech., 634, 291 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, C. A., Boronski, P., Brun, A. S., et al. 2011, Icarus, 216, 120 [NASA ADS] [CrossRef] [Google Scholar]
 Kageyama, A., & Sato, T. 1997, Phys. Rev. E, 55, 4617 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J., Mantere, M. J., Cole, E., Warnecke, J., & Brandenburg, A. 2013, ApJ, 778, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Katayama, J. S., Matsushima, M., & Honkura, Y. 1999, Phys. Earth Planet. Inter., 111, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Kuang, W., & Bloxham, J. 1997, Nature, 389, 371 [NASA ADS] [CrossRef] [Google Scholar]
 Kutzner, C., & Christensen, U. R. 2002, Phys. Earth Planet. Inter., 131, 29 [Google Scholar]
 Landeau, M., & Aubert, J. 2011, Phys. Earth Planet. Inter., 185, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Lantz, S. R., & Fan, Y. 1999, ApJS, 121, 247 [Google Scholar]
 Miesch, M. S., Elliott, J. R., Toomre, J., et al. 2000, ApJ, 532, 593 [NASA ADS] [CrossRef] [Google Scholar]
 Müller, U., & Stieglitz, R. 2002, Nonlinear Processes in Geophysics, 9, 165 [NASA ADS] [CrossRef] [Google Scholar]
 Ogura, Y., & Phillips, N. A. 1962, J. Atmos. Sci., 19, 173 [NASA ADS] [CrossRef] [Google Scholar]
 Olson, P., Christensen, U. R., & Glatzmaier, G. A. 1999, J. Geophys. Res., 104, 10383 [NASA ADS] [CrossRef] [Google Scholar]
 Parker, E. N. 1955, ApJ, 122, 293 [NASA ADS] [CrossRef] [Google Scholar]
 Rempel, M. 2006, ApJ, 647, 662 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, T. M., & Glatzmaier, G. A. 2005, ApJ, 620, 432 [NASA ADS] [CrossRef] [Google Scholar]
 Ruediger, G. 1980, Astron. Nachr., 301, 181 [NASA ADS] [CrossRef] [Google Scholar]
 Sarson, G. R., Jones, C. A., & Longbottom, A. W. 1998, Geophys. Astrophys. Fluid Dyn., 88, 225 [NASA ADS] [CrossRef] [Google Scholar]
 Schmitz, S., & Tilgner, A. 2009, Phys. Rev. E, 80, 015305 [NASA ADS] [CrossRef] [Google Scholar]
 Schrinner, M. 2011, A&A, 533, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schrinner, M. 2013, MNRAS, 431, L78 [NASA ADS] [CrossRef] [Google Scholar]
 Schrinner, M., Rädler, K.H., Schmitt, D., Rheinhardt, M., & Christensen, U. R. 2007, Geophys. Astrophys. Fluid Dyn., 101, 81 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Schrinner, M., Schmitt, D., Cameron, R., & Hoyng, P. 2010, Geophys. J. Int., 182, 675 [NASA ADS] [CrossRef] [Google Scholar]
 Schrinner, M., Petitdemange, L., & Dormy, E. 2011a, A&A, 530, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schrinner, M., Schmitt, D., & Hoyng, P. 2011b, Phys. Earth Planet. Inter., 188, 185 [NASA ADS] [CrossRef] [Google Scholar]
 Schrinner, M., Petitdemange, L., & Dormy, E. 2012, ApJ, 752, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Sreenivasan, B., & Jones, C. A. 2006, Geophys. J. Int., 164, 467 [NASA ADS] [CrossRef] [Google Scholar]
 Stelzer, Z., & Jackson, A. 2013, Geophys. J. Int., 193, 1265 [NASA ADS] [CrossRef] [Google Scholar]
 Tilgner, A. 1997, Phys. Lett. A, 226, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Tilgner, A. 2012, Phys. Rev. Lett., 109, 248501 [NASA ADS] [CrossRef] [Google Scholar]
 Tilgner, A., & Brandenburg, A. 2008, MNRAS, 391, 1477 [NASA ADS] [CrossRef] [Google Scholar]
 Yadav, R. K., Gastine, T., & Christensen, U. R. 2013a, Icarus, 225, 185 [NASA ADS] [CrossRef] [Google Scholar]
 Yadav, R. K., Gastine, T., Christensen, U. R., & Duarte, L. D. V. 2013b, ApJ, 774, 6 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Benchmark results
Hydrodynamic benchmark (E = 10^{3},N_{ϱ} = 5,χ = 0.35,Ra = 3.52 × 10^{5},Pr = 1,n = 2).
Steady dynamo benchmark ( E = 2 × 10^{3},N_{ϱ} = 3,χ = 0.35,Ra = 8.00 × 10^{4},Pr = 1,Pm = 50,n = 2).
Unsteady dynamo benchmark (E = 5 × 10^{5},N_{ϱ} = 3,χ = 0.35,Ra = 2.50 × 10^{7},Pr = 2,Pm = 2,n = 2).
All Tables
Overview of the critical Rayleigh numbers and the corresponding critical azimuthal wavenumbers for the onset of convection determined, as explained in Jones et al. (2009).
Overview of the simulations carried out, ordered with respect to their local Rossby number.
Hydrodynamic benchmark (E = 10^{3},N_{ϱ} = 5,χ = 0.35,Ra = 3.52 × 10^{5},Pr = 1,n = 2).
Steady dynamo benchmark ( E = 2 × 10^{3},N_{ϱ} = 3,χ = 0.35,Ra = 8.00 × 10^{4},Pr = 1,Pm = 50,n = 2).
Unsteady dynamo benchmark (E = 5 × 10^{5},N_{ϱ} = 3,χ = 0.35,Ra = 2.50 × 10^{7},Pr = 2,Pm = 2,n = 2).
All Figures
Fig. 1 Relative dipole fieldstrength versus the local Rossby number for our sample of models. Filled symbols stand for dipolar, open symbols for multipolar dynamos. The symbol shape indicates the number of density scale heights: N_{ϱ} = 0.5: circle; N_{ϱ} = 1: upward triangle; N_{ϱ} = 1.5: downward triangle; N_{ϱ} = 2: diamond; N_{ϱ} = 2.5: square; N_{ϱ} = 3,3.5,4: star. A cross inscribed in some open symbols means that the field of these models exhibits a strong equatorial dipole component. 

In the text 
Fig. 2 Contour plot of the azimuthally averaged radial magnetic field of model3m versus time and colatitude. The contour plot was normalised by the maximum absolute value at each time step. The greyscale coding ranges from −1, white, to +1, black. 

In the text 
Fig. 3 Evolution of the energy of the tracer field normalised by the energy of the actual magnetic field for model2m (dashed line) and model54d (solid line). 

In the text 
Fig. 4 Relative dipole field strength versus Ro_{ℓ} for a sequence of models with E = 3 × 10^{4}, Ra = 4Ra_{c}, Pm = 3, and Pr = 1. The meaning of the symbols is defined in the caption of Fig. 1. 

In the text 
Fig. 5 Contour plot (equatorial cut) of the radial magnetic field of model2m at a given time. 

In the text 
Fig. 6 Dotted line: axisymmetric kinetic energy of model5m normalised by an arbitrary value. Solid line: ratio of axisymmetric to total magnetic energy. 

In the text 
Fig. 7 Magnetic Prandtl number versus N_{ϱ} for models with E = 10^{4} and Pr = 1 and variable Rayleigh numbers. Filled circles stand for parameters for which dipolar solutions were obtained. 

In the text 
Fig. 8 Lorentz number compensated by f_{ohm} versus the fluxbased Rayleigh number for our sample of models. Filled symbols correspond to dipolar models, open symbols are multipolar models and the symbol shape indicates N_{ρ} as explained in the caption of Fig. 1. 

In the text 
Fig. 9 Rossby number versus the fluxbased Rayleigh number for our sample of models. 

In the text 
Fig. 10 Ohmic dissipation time versus Rossby number for all models considered in this study. 

In the text 
Fig. 11 Nusselt number versus the fluxbased Rayleigh number for our sample of models. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.