Free Access
Issue
A&A
Volume 567, July 2014
Article Number A107
Number of page(s) 9
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/201423902
Published online 22 July 2014

© ESO, 2014

1. Introduction

Dynamo action, i.e. the self-amplification of a magnetic field by the flow of an electrically conducting fluid, is considered to be the main mechanism for generating magnetic fields in the universe for a variety of systems, including planets, stars, and galaxies (e.g. Dormy & Soward 2007). Dynamo action is an instability by which a conducting fluid transfers part of its kinetic energy to magnetic energy. Because of the difficulty simulating turbulent fluid motions, one must resort to some approximations to model the fluid flow, whose convective motions are assumed to be driven by the temperature difference between a hot inner core and a cooler outer surface. A strong simplification can be achieved when applying the Boussinesq approximation (Boussinesq 1903), which performs well in so far as variations in pressure scarcely affect the density of the fluid. However, in essence, this approximation will not be adequate for describing convection in highly stratified systems, such as stars or gas giants. A common approach to overcoming this difficulty is then to use the anelastic approximation, which allows for a reference density profile while filtering out sound waves for faster numerical integration. This approximation was first developed to study atmospheric convection (Ogura & Phillips 1962; Gough 1969). It has then been used to model convection in the Earth core or in stars and is found in the literature under slightly different formulations (Gilman & Glatzmaier 1981; Braginsky & Roberts 1995; Lantz & Fan 1999; Anufriev et al. 2005; Berkoff et al. 2010; Jones et al. 2011; Alboussière & Ricard 2013). Nevertheless, the starting point in the anelastic approximation is always to consider convection as a perturbation of a stratified reference state that is assumed to be close to adiabatic.

Observations of low mass stars have revealed very different magnetic field topologies from small scale fields to large scale dipolar fields (Donati & Landstreet 2009; Morin et al. 2010), and highlight possible correlations between differential rotation and magnetic field topologies (Reinhold et al. 2013). Boussinesq models partly reproduce this diversity (Busse & Simitev 2006; Morin et al. 2011; Sasaki et al. 2011; Schrinner et al. 2012). Moreover, the dichotomy between dipolar and “non-dipolar” (or multipolar) dynamos seems to hold for anelastic models (Gastine et al. 2012; Yadav et al. 2013b). However, we show in Schrinner et al. (2014) that this distinction may somewhat be less clear than it was with Boussinesq models. Indeed, in a systematic parameter study we found a large number of models with both a high equatorial dipole contribution and an intermediate dipole field strength. Only a few examples of equatorial dipoles have been reported from Boussinesq spherical dynamo simulations (Aubert & Wicht 2004; Gissinger et al. 2012). At the same time, observations have shown that for some planets, such as Uranus or Neptune, the dipole axis can make an angle up to π/2 with respect to the rotation axis, owing to a significant contribution from the equatorial dipole (Jones 2011).

In this paper, we aim to clarify the reasons likely for the emergence of an equatorial dipole contribution when measuring the dipole field strength at the surface of numerical models. Since our approach closely follows previous methodology for studying the link with Boussinesq results, we decided to focus in more detail on one important change that comes with the anelastic approximation as formulated in Jones et al. (2011), assuming that all mass is concentrated inside the inner sphere to determine the gravity profile. In contrast, as proposed by the Boussinesq dynamo benchmark Christensen et al. (2001), it was common for geodynamo studies to assume that the density is homogeneously distributed. This leads to different gravity profiles, the first being proportional to 1/r2, whereas the second is proportional to r. According to Duarte et al. (2013), Gastine et al. (2012) show that both gravity profiles lead to very similar results. Contrary to this statement, we show that the choice of the gravity profile may have strong consequences on the dynamo-generated field topology. We briefly recall the anelastic equations in Sect. 2 before presenting our results in Sect. 3. In Appendix A, we give the fit coefficients obtained for the scalings of the magnetic and velocity fields and the convective heat flux. A summary of the numerical simulations carried out is given in Appendix B.

2. Governing equations

2.1. The non-dimensional anelastic equations

We rely on the LBR-formulation, named after Lantz & Fan (1999) and Braginsky & Roberts (1995), as it is used in the dynamo benchmarks proposed by Jones et al. (2011). It guarantees the energy conservation, unlike other formulations (see Brown et al. 2012). A more detailed presentation of the equations can be found in our preceding paper Schrinner et al. (2014).

Let us consider a spherical shell of width d and aspect ratio χ, rotating about the z axis at angular velocity Ω and filled with a perfect, electrically conducting gas with kinematic viscosity ν, thermal diffusivity κ, specific heat cp, and magnetic diffusivity η (all assumed to be constant). In contrast to the usual Boussinesq framework, convection is driven by an imposed entropy difference Δs between the inner and the outer boundaries, and the gravity is given by g=GMrˆ/r2\hbox{${\vec g}=-GM{\hat{\vec r}}/r^2$} where G is the gravitational constant and M the central mass, assuming that the bulk of the mass is concentrated inside the inner sphere.

The reference state is given by the polytropic equilibrium solution of the anelastic system P=Pcwn+1,ϱ=ϱcwn,T=Tcw,w=c0+c1dr,c0=2w0χ11χ,c1=(1+χ)(1wo)(1χ)2,\begin{eqnarray} \label{ref_state} \mP&=&P_{\rm c}\,w^{n+1},\quad\mrho=\varrho_{\rm c}\,w^n,\quad \mT=T_{\rm c}\,w,\quad w=c_0+\frac{c_1 d}{r}, \\ \label{c0_c1} c_0&=&\frac{2w_0-\aspectratio-1}{1-\aspectratio},\quad c_1=\frac{(1+\aspectratio)(1-w_{\rm o})}{(1-\aspectratio)^2}, \end{eqnarray}with w0=χ+1χexp(Nϱ/n)+1,wi=1+χwoχ·\begin{eqnarray} w_0=\frac{\aspectratio+1}{\aspectratio\exp(N_\varrho/n)+1},\quad w_{\rm i}=\frac{1+\aspectratio-w_{\rm o}}{\aspectratio}\cdot \label{w1_w0} \end{eqnarray}(3)In the above expressions, n is the polytropic index and Nϱ the number of density scale heights, defined by Nϱ = ln(ϱi/ϱo), where ϱi and ϱo denote the reference state density at the inner and outer boundaries, respectively. The values Pc, ϱc, and Tc are the reference-state density, pressure, and temperature midway between the inner and outer boundaries, and serve as units for these variables.

We adopt the same non-dimensional form as Jones et al. (2011): length is scaled by the shell width d, time by the magnetic diffusion time d2/η, and entropy by the imposed entropy difference Δs. The magnetic field is measured in units of Ωϱcμη\hbox{$\sqrt{\Omega\varrho_{\rm c}\mu\eta}$} where μ is the magnetic permeability. Then, the equations governing the system are v∂t+v·v=Pm[1EPwn+PmPrRasr2rˆ2Ezˆ×v+Fν+1Ewn(×B)×B],B∂t=×(v×B)+2B,∂s∂t+v·s=wn1PmPr·(wn+1s)+Diw[E-1wn(×B)2+Qν],·(wnv)=0,·B=0.\begin{eqnarray} \label{mhd1} \frac{\partial{\vec v}}{\partial t} + {\vec v}\cdot\nabla{\vec v} &=& Pm\, \bigg[-\frac{1}{E}\nabla\frac{P'}{w^n} +\frac{Pm}{Pr}Ra\frac{s}{r^2} {\hat{\vec r}} -\frac{2}{E}\, {\hat{\vec z}} \times{\vec v} \nonumber\\ &&\quad + {\vec F}_\nu+\frac{1}{E\,w^n}(\nabla\times{\vec B})\times{\vec B} \bigg]\,, \\ \label{mhd2} \frac{\partial{\vec B}}{\partial t} &= & \nabla\times({\vec v}\times{\vec B})+\nabla^2{\vec B} ,\\ \label{mhd3} \frac{\partial s}{\partial t}+{\vec v}\cdot\nabla s &= & w^{-n-1}\frac{Pm}{Pr}\nabla\cdot\left(w^{n+1}\,\nabla s\right)\nonumber \\ \label{mhd4} &&\quad + \frac{Di}{w}\left[E^{-1}w^{-n}(\nabla\times {\vec B})^2+Q_\nu\right], \\ \label{mhd5} \nabla\cdot \left(w^n{\vec v} \right) &=& 0\,,\\ \nabla\cdot{\vec B} &=& 0. \end{eqnarray}The viscous force Fν in (4)is given by Fν = wnS, where S is the rate of strain tensor Sij=2wn(eij13δij·v),eij=12(vixj+vjxi)·\begin{eqnarray} S_{ij}=2w^n\left(e_{ij}-\frac{1}{3}\delta_{ij}\nabla\cdot {\vec v}\right),\quad e_{ij}=\frac{1}{2}\left(\frac{\partial v_{\rm i}}{\partial x_j}+ \frac{\partial v_j}{\partial x_{\rm i}}\right) \cdot \end{eqnarray}(9)Moreover, the expressions of the dissipation parameter Di and the viscous heating Qν in (6) are Di=c1PrPmRa,\begin{eqnarray} Di=\frac{c_1Pr}{Pm Ra}, \end{eqnarray}(10)and Qν=2[eijeij13(·v)2]·\begin{eqnarray} Q_\nu=2\left[e_{ij}e_{ij}-\frac{1}{3}(\nabla\cdot{\vec v})^2\right]\cdot \end{eqnarray}(11)The boundary conditions are the following. We impose stress free boundary conditions for the velocity field at both the inner and the outer sphere, the magnetic field matches a potential field inside and outside the fluid shell, and the entropy is fixed at the inner and outer boundaries.

The system of Eqs. (4)–(7) involves seven control parameters, namely the Rayleigh number Ra = GMdΔs/(νκcp), the Ekman number E = ν/(Ωd2), the Prandtl number Pr = ν/κ, and the magnetic Prandtl number Pm = ν/η, together with the aspect ratio χ, the polytropic index n, and the number of density scale heights Nϱ that define the reference state. In this study, we restrict our investigation of the parameter space keeping E = 10-4, Pr = 1, χ = 0.35, and n = 2 for all simulations. Furthermore, to differentiate the effects related to the change in gravity profile from those related to the anelastic approximation, we decided to perform low Nϱ simulations so that we can assume that stratification no longer influences the dynamo process. In practice, we chose Nϱ = 0.1, which means that the density contrast between the inner and outer spheres is only 1.1, and the simulations are thus very close to the Boussinesq limit. To further ensure the lack of stratification effects, we also checked in a few cases that the results do not differ from purely Boussinesq simulations.

We have integrated our system at least on one magnetic diffusion time with the anelastic version of parody (Dormy et al. (1998) and further developments), which reproduces the anelastic dynamo benchmarks (see Schrinner et al. 2014). The vector fields are transformed into scalars using the poloidal-toroidal decomposition. The equations are then discretized in the radial direction with a finite-difference scheme; on each concentric sphere, variables are expanded using a spherical harmonic basis. The coefficients of the expansion are identified with their degree and order m. Typical resolutions use from 256 to 288 points in the radial direction and a spectral decomposition truncated at 80 ≤ max ~ mmax ≤ 116.

2.2. Diagnostic parameters

The quantities used to analyse our simulations first rely on the kinetic and magnetic energy densities Ek and Em, Ek=12VVwnv2dvandEm=12VPmEVB2dv.\begin{eqnarray} \Ek = \frac{1}{2\, V}\int_V\,w^n{\vec v}^2 \mathrm{d}v \quad \text{and} \quad \Em = \frac{1}{2\, V}\frac{Pm}{E}\int_V\,{\vec B}^2 \mathrm{d}v. \label{emag} \end{eqnarray}(12)We define the corresponding Rossby number Ro=2EkE/Pm\hbox{$Ro=\sqrt{2 \Ek} E/Pm$} and Lorentz number Lo=2EmE/Pm\hbox{$\Lo=\sqrt{2 \Em}E/Pm$}. The latter is a non-dimensional measure of the magnetic field amplitude, while the former is a non-dimensional measure of the velocity field amplitude. A measure of the mean zonal flow is the zonal Rossby number Roz, whose definition is based on the averaged toroidal axisymmetric kinetic energy density. To distinguish between dipolar and multipolar dynamo regimes, we know from Boussinesq results that it is useful to measure the balance between inertia and Coriolis force, which can be approximated in terms of a local Rossby number Ro = Rocc/π, which depends on the characteristic length scale of the flow rather than on the shell thickness (Christensen & Aubert 2006; Olson & Christensen 2006; Schrinner et al. 2012). Again, we emphasize that our definition of the local Rossby number tries to avoid any dependence on the mean zonal flow and thus differs from the original definition introduced by Christensen & Aubert (2006)(see Schrinner et al. 2012, App. A for a discussion). Our typical convective length scale is based on the mean harmonic degree c of the velocity component vc from which the mean zonal flow has been subtracted, c=wn(vc)·(vc)wnvc·vc·\begin{eqnarray} \lc = \sum_\ell \ell\frac{\langle w^n\, ({\vec v}_{\rm c})_\ell\cdot({\vec v}_{\rm c})_\ell\rangle} {\langle w^n \, {\vec v}_{\rm c}\cdot{\vec v}_{\rm c} \rangle}\cdot \end{eqnarray}(13)Where the brackets denote an average of time and radii. Consistently, the contribution of the mean zonal flow is removed for calculating Roc.

The dipolarity of the magnetic field is characterized by the relative dipole field strength, fdip, originally defined as the time-average ratio on the outer shell boundary So of the dipole field strength to the total field strength, fdip=SoB2=1m={0,1}sinθdθdφSoB2sinθdθdφt·\begin{eqnarray} \fdip = \left\langle \sqrt{\frac {\int_{S_{\!o}} {{\vec B}^2}_{\,\ell\,=\,1}^{\,m=\left \{0,1\right \}} \, \sin \theta \mathrm{d}\theta \mathrm{d}\phi }{\int_{S_{\!o}} {\vec B}^2 \, \sin \theta \mathrm{d}\theta \mathrm{d}\phi } } \right\rangle_t\cdot \end{eqnarray}(14)We also define a relative axial dipole field strength filtering out non-axisymmetric contributions fdipax=SoB2=1m=0sinθdθdφSoB2sinθdθdφt·\begin{eqnarray} \fdipAX = \left\langle \sqrt{\frac{\int_{S_{\!o}} {{\vec B}^2}_{\,\ell\,=\,1}^{\,m\,=\,0} \, \sin \theta \mathrm{d}\theta \mathrm{d}\phi }{\int_{S_{\!o}} {\vec B}^2 \, \sin \theta \mathrm{d}\theta \mathrm{d}\phi } } \right\rangle_t\cdot\label{e_fdipAX} \end{eqnarray}(15)This definition of fdipax is similar to the relative dipole field strength used by Gastine et al. (2012), except for the square root, which explains the lower values for the dipolarity found in Gastine et al. (2012).

To further characterize the topology of the magnetic field, we introduce a time-averaged measure of the departure from a pure equatorial dipole solution θ=2π(Θ(t)π2)2t,\begin{eqnarray} \theta = \frac{2}{\pi} \left\langle \sqrt{\left(\Theta(t) - \frac{\pi}{2} \right)^2}\right\rangle_t \,,\label{e_theta} \end{eqnarray}(16)where Θ is the tilt angle of the dipole axis. A low value of θ indicates that the tilt angle of the dipole fluctuates close to Θ = π/2.

3. Results

3.1. Bifurcations between dynamo branches

thumbnail Fig. 1

Dipolar (black circles) and multipolar (white squares) dynamos as a function of Ra/Rac and Pm, for a central mass a) and a uniform mass distribution b). Crosses indicate the absence of a self-sustained dynamo.

thumbnail Fig. 2

a) Elsasser number Λ as a function of Ra/Rac, for Pm = 1 (green) and Pm = 3 (red). The meaning of the symbol shapes is defined in the caption of Fig. 1. A grey marker indicates that the solution loses its stability. b) Detail of the bifurcation close to the dynamo threshold for Pm = 1. Error bars represent the standard deviations.

In our simulations we recover the two distinct dynamo regimes observed with both Boussinesq (Kutzner & Christensen 2002; Christensen & Aubert 2006; Schrinner et al. 2012; Yadav et al. 2013a) and anelastic models (Gastine et al. 2012; Schrinner et al. 2014). These are characterized by different magnetic field configurations: dipolar dynamos are dominated by a strong axial dipole component, whereas “non-dipolar” dynamos usually present a more complex geometry with higher spatial and temporal variability. The branches are easily identified by continuing simulations performed with other parameters, for which the dipolar/multipolar characteristic was previously established. Figure 1a shows the regime diagram we obtained, as a function of the Rayleigh and magnetic Prandtl numbers. For comparison, we use the data from Schrinner et al. (2012) and show in Fig. 1b the same regime diagram obtained for Boussinesq models with a uniform mass distribution. For Pm = 1, the transition from the dipolar to the multipolar branch can be triggered by an increase in Ra. In that case, the transition is due to the increasing role of inertia as revealed by Ro. Alternatively, the transition from multipolar to dipolar dynamo can be triggered by increasing Pm. Then, the multipolar branch is lost when the saturated amplitude of the mean zonal flow becomes too small to prevent the growth of the dipolar solution (see Schrinner et al. 2012). It is worth noting that the two branches overlap for a restricted parameter range for which dipolar and multipolar dynamos may coexist. In that case, the observed solution strongly depends of the initial magnetic field, so we tested both weak and strong field initial conditions for all our models to delimit the extent of the bi-stable zone with greater accuracy. Actually, multipolar dynamos are favoured by the stronger zonal wind that may develop with stress-free boundary conditions, allowing for this hysteretic transition (Schrinner et al. 2012). Finally, we see that the dynamo threshold is lower for multipolar models, which allows the multipolar branch to extend below the dipolar branch at low Rayleigh and magnetic Prandtl numbers. We see in Fig. 1b that this is different from Boussinesq models with a uniform mass distribution.

To investigate the different transitions between the different dynamo branches, we plot the Elsasser number Λ=Brms2/(Ωϱcμη)\hbox{$\Elsasser = B^2_\text{rms}/(\Omega \varrho_{\rm c} \mu \eta) $} in Fig. 2a (related to the Lorentz number by Λ = Lo2Pm/E) as a function of the distance to the threshold for models at Pm = 1 and Pm = 3. We see in Fig. 2b that the bifurcation for multipolar branch is supercritical. When decreasing the Rayleigh number, the dipolar branch loses its stability for Ra/Rac ~ 20, when the magnetic field strength becomes too weak.

thumbnail Fig. 3

Λdip/Rozmul\hbox{$\Lambda^\mathrm{dip}/\Roz^\mathrm{mul}$} as a function of Ra/Rac, for Pm = 1 (green diamonds), Pm = 2 (blue pentagons) and Pm = 3 (red stars). The point marked with the grey star has been computed with the model corresponding to the grey square in Fig. 2.

For higher magnetic Prandtl numbers, the bifurcation of the multipolar branch still seems to be supercritical. Interestingly, one notes in Fig. 2a for Pm = 3 that the multipolar branch loses its stability when increasing the Rayleigh number. A physical explanation for this behaviour is that the mean zonal flow does not grow fast enough as the field strength increases, and the dynamo switches to the dipolar solution. This simple physical scenario can be illustrated by comparing the variation in the field strength of the dipolar branch, as measured by Λdip, and the zonal shear of the multipolar branch, as measured by Rozmul\hbox{$\Roz^\mathrm{mul}$}. Indeed, we see in Fig. 3 that the higher the magnetic Prandtl number, the faster the growth of the ratio between Λdip and Rozmul\hbox{$\Roz^\mathrm{mul}$}. This explains why the multipolar branch destabilizes at large forcing for larger Pm (Pm = 3, red dashed line in Fig. 2a), while it remains stable at smaller Pm (Pm = 1, dashed green line in Fig. 2a). Because of computational limitations, we were not able to find for Pm> 1 the Rayleigh numbers for which the dipolar branch should disappear.

thumbnail Fig. 4

a) Relative dipole field strength fdip versus local Rossby number. b) Relative axial dipole field strength fdipax versus local Rossby number. The meaning of the symbol shapes is defined in the caption of Fig. 1.

3.2. Equatorial dipole

Schrinner et al. (2014) show that dipolar and multipolar dynamos in anelastic simulations were no longer distinguishable from each other in terms of fdip, contrary to Boussinesq models. This smoother transition has been attributed to the presence of dynamos with a high equatorial dipole contribution, which leads to intermediate values for fdip. However, Fig. 4a shows that this tendency already exits at low Nϱ, and thus cannot be accounted for only in terms of anelastic effects. Furthermore, when the equatorial dipole component is removed to compute the relative dipole field strength, we recover a more abrupt transition, as we can see in Fig. 4b which shows the relative axial dipole field strength fdipax. Dipolar dynamos are left unchanged by this new definition, whereas multipolar dynamos of intermediate dipolarity are no longer observed, which confirms that the increase in fdip is due to a significant equatorial dipole component. The quantity fdipax therefore provides a robust criterion to distinguish the dipolar and the multipolar branches.

thumbnail Fig. 5

a) Evolution of the modified tilt angle θ defined by Eq. (16)for multipolar dynamos as a function of Ra/Rac and Pm. Colour scale ranges from white (θ = 0) to black (θ = 0.7). b) θ as a function of Ra/Rac for Pm = 1. Upper x axis corresponds to the values of Ro for the multipolar branch. The meaning of the symbol shapes is defined in the caption of Fig. 1. Error bars represent the standard deviations.

To further characterize the emergence of multipolar dynamos with a significant equatorial dipole contribution, we plot the values of the modified tilt angle θ in the parameter space in Fig. 5a. Low values of θ are characteristic of an equatorial dipole on the surface of the outer sphere and they appear to be preferably localized close to the dynamo threshold of the multipolar branch, at low Rayleigh and magnetic Prandtl numbers. In our case, dynamos with a stronger equatorial dipole component belong to the class of multipolar dynamos, but since they are always close to the threshold, fewer modes are likely to be excited. As the Rayleigh number or the magnetic Prandtl number is increased, the dipole axis is not stable anymore but fluctuates in the interval [ 0 ], which is typical of polarity reversals for multipolar dynamos (Kutzner & Christensen 2002). This evolution is illustrated in Fig. 5b for dynamos at Pm = 1. For this subset of models, we computed the percentage of the non-axisymmetric magnetic energy density with respect to the total magnetic energy density Em and saw that it tends to increase from 85% on average for multipolar dynamos up to 93% as the Rayleigh number is decreased.

thumbnail Fig. 6

Equatorial cross-section of the radial component of the velocity, with E = 10-4, Pr = 1. a)gr and Ra/Rac = 9.0 , Pm = 1. b)g ∝ 1/r2 and Ra/Rac = 9.9, Pm = 1.2.

thumbnail Fig. 7

Equatorial cross-section of the radial component of the magnetic field,with E = 10-4, Pr = 1. a) gr and Ra/Rac = 9.0 , Pm = 1. b), c)g ∝ 1/r2 and Ra/Rac = 9.9, Pm = 1.2. Colour in panel c) has been rescaled to highlight the emergence of a m = 1 mode at the outer sphere.

Part of the changes reported in Schrinner et al. (2014) about anelastic dynamos simulations do not seem to come from the stratified reference density profile, but from the choice of a gravity profile proportional to 1/r2. This profile differs from the gravity profile proportional to r that was used for Boussinesq simulations and is actually the only significant difference between previous studies and our low Nϱ simulations. As a consequence, convection cells form and stay closer to the inner sphere, as we can see in Fig. 6. We compare here equatorial cuts of the radial component of the velocity for both choices of gravity profile. This strong difference in the flow reflects on the localization of the active dynamo regions, as we can see in the corresponding cuts of the radial component of the magnetic field in Fig. 7. With a gravity profile proportional to 1/r2, the magnetic field is mainly generated close to the inner sphere where the convection cells form. Consequently, our measure of the dipole field strength fdip at the surface of the outer sphere appears to be biased, since it will essentially be sensitive to the less diffusive large scale modes. This filter effect is likely to be responsible for the increase in fdip we reported in some anelastic dynamo models. However, for higher density stratification Nϱ = 3 and Prandtl numbers Pr = 2 and Pm = 4, Schrinner et al. (2014) identify equatorial dipole dynamos with a m = 1 component that is not localized on the outer sphere and for which the present mechanism will not be relevant.

4. Conclusion

In this paper, we focussed on very weakly stratified anelastic dynamo models with a central mass distribution. We investigated the bifurcations between the dipolar and multipolar dynamo branches and recovered in parts the behaviour that has been observed for Boussinesq models with a uniform mass distribution. In addition, we show that the dipolar branch can now lose its stability and switch to the multipolar branch at low Rayleigh and magnetic Prandtl numbers. The multipolar dynamos that are observed in this restricted parameter regime usually present a stronger equatorial dipole component at the surface of the outer sphere. When increasing the Rayleigh number at higher Pm, it seems that the mean zonal flow does not grow fast enough to maintain the multipolar solution, and we identified several cases where the multipolar branch switches to the dipolar solution.

This study has shed interesting light on the recent systematic parameter study of spherical anelastic dynamo models started by Schrinner et al. (2014). Focussing on weakly stratified dynamo models, we showed that magnetic field configurations with a significant equatorial dipole contribution can already be observed in the Boussinesq limit. In the parameter regime we investigate, our study reveals that the choice of gravity profile has a strong influence on the fluid flow and thus on the dynamo generated magnetic field, depending whether one considers a uniform or a central mass distribution. In the parameter space, we showed that multipolar dynamos with a significant equatorial dipole contribution are preferably observed close to the dynamo threshold. This is reminiscent of the results of Aubert & Wicht (2004), who studied the competition between axial and equatorial dipolar dynamos when varying the aspect ratio of the spherical shell in geodynamo models.

However, the filter effect we highlight here focusses on the topology of the magnetic field at the outer surface of the models. As a consequence, it will not be able to explain the stronger equatorial dipole component for all models in Schrinner et al. (2014). Independently, Cole et al. (2014) also report the discovery of an azimuthal dynamo wave of a m = 1 mode in numerical simulations corresponding to higher density stratification. Their Coriolis number plays a similar role to the inverse of our local Rossby number. Upper x axis of Fig. 5b shows that equatorial dipole configurations are favoured with the decrease in Ro.

Observational results from photometry (Hackman et al. 2013) and spectropolarimetry (Kochukhov et al. 2013) of rapidly rotating cool active stars reveal that the surface magnetic field of these objects can be highly non-axisymmetric. Further investigation of direct numerical simulations is therefore required to better understand the influence of the Prandtl number and the density stratification on the magnetic field topology.

Online material

Appendix A: Scaling laws

Table A.1

Summary of the coefficients obtained for the different scaling laws, with their standard error from the linear regression.

For our samples of models we compute the usual scaling laws that have been derived for the magnetic and velocity fields and the convective heat flux (Christensen & Aubert 2006; Yadav et al. 2013a; Stelzer & Jackson 2013; Yadav et al. 2013b). As in Schrinner et al. (2014), we do not attempt to solve any secondary dependence on Pm because we do not vary this parameter on a wide enough range. We transform our problem to a linear one by taking the logarithm and look for a law of the form lnŷ = β + αlnx. To quantify the misfit between data and fitted values, we follow Christensen & Aubert (2006) and compute the relative misfit χrel=1ni=1n(yiyiˆyi)2,\appendix \setcounter{section}{1} \begin{eqnarray} \chir = \sqrt{\frac{1}{n} \sum_{i=1}^n \left(\frac{y_{\rm i}-\hat{y_{\rm i}}}{y_{\rm i}} \right)^2}, \end{eqnarray}(A.1)where yiˆ\hbox{$\hat{y_{\rm i}}$} stands for predicted values, yi for measured values, and n is the number of data. Our results are summarized in Table A.1 and compared with those found in Schrinner et al. (2014) and Yadav et al. (2013a) in Table A.2. We did not find significant differences between the anelastic scalings and the scaling we obtained in the Boussinesq limit with the same mass distribution. The coefficients we obtained seem closer on average to the coefficients obtained by Yadav et al. (2013a) with Boussinesq models with a uniform mass distribution. However, it is not possible to deduce from our data set any influence of Nϱ on the different coefficients of the scaling laws. In our models, the Nusselt number evaluated at the surface of the inner sphere Si is defined by Nubot=(exp(Nϱ)1)wiri24πnc1Si∂s∂rsinθdθdφ,\appendix \setcounter{section}{1} \begin{eqnarray} Nu_\mathrm{bot} = -\frac{(\exp{(N_\varrho)}-1)w_{\rm i}r_{\rm i}^2}{4\pi nc_1}\int_{S_{\!i}}\frac{\partial s}{\partial r}\sin\theta \mathrm{d}\theta \mathrm{d}\phi, \label{nu_bot} \end{eqnarray}(A.2)which enables a flux based Rayleigh number to be defined, RaQ=NuRaE2ro2PrwithNu=(Nubot1)EPr,\appendix \setcounter{section}{1} \begin{eqnarray} Ra_Q=Nu^\star\,\frac{Ra\, E^2}{\ro^2\, Pr} \quad \text{with} \quad Nu^\star = (Nu_\mathrm{bot}-1)\,\frac{E}{Pr}, \label{def_raq} \end{eqnarray}(A.3)which is used in derivating scaling laws, together with the fraction of ohmic to total dissipation fohm.

Table A.2

Comparison between the different scaling laws obtained with different dynamo models.

Appendix B: Numerical models

Table B.1

Overview of the simulations carried out with E = 10-4, Pr = 1, χ = 0.35, n = 2, and Nϱ = 0.1.

Acknowledgments

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 ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the French Agence Nationale pour la Recherche. Numerical simulations were also carried out at CEMAG and TGCC computing centres (GENCI project x2013046698). L. P. acknowledges financial support from “Programme National de Physique Stellaire” (PNPS) of CNRS/INSU, France.

References

  1. Alboussière, T., & Ricard, Y. 2013, J. Fluid Mech., 725, 1 [NASA ADS] [CrossRef] [Google Scholar]
  2. Anufriev, A. P., Jones, C. A., & Soward, A. M. 2005, Phys. Earth Planet. Inter., 152, 163 [NASA ADS] [CrossRef] [Google Scholar]
  3. Aubert, J., & Wicht, J. 2004, Earth Planet. Sci. Lett., 221, 409 [NASA ADS] [CrossRef] [Google Scholar]
  4. Berkoff, N. A., Kersale, E., & Tobias, S. M. 2010, Geophys. Astrophys. Fluid Dyn., 104, 545 [NASA ADS] [CrossRef] [Google Scholar]
  5. Boussinesq, J. 1903, Théorie analytique de la chaleur, Vol. 2 (Gauthier-Villars) [Google Scholar]
  6. Braginsky, S. I., & Roberts, P. H. 1995, Geophys. Astrophys. Fluid Dyn., 79, 1 [NASA ADS] [CrossRef] [Google Scholar]
  7. Brown, B. P., Vasil, G. M., & Zweibel, E. G. 2012, ApJ, 756, 109 [Google Scholar]
  8. Busse, F. H., & Simitev, R. D. 2006, Geophys. Astrophys. Fluid Dyn., 100, 341 [NASA ADS] [CrossRef] [Google Scholar]
  9. Christensen, U. R., & Aubert, J. 2006, Geophys. J. Int., 166, 97 [NASA ADS] [CrossRef] [Google Scholar]
  10. Christensen, U. R., Aubert, J., Cardin, P., et al. 2001, Phys. Earth Planet. Inter., 128, 25 [Google Scholar]
  11. Cole, E., Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2014, ApJ, 780, L22 [NASA ADS] [CrossRef] [Google Scholar]
  12. Donati, J.-F., & Landstreet, J. D. 2009, ARA&A, 47, 333 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  13. Dormy, E., Cardin, P., & Jault, D. 1998, Earth Planet. Sci. Lett., 160, 15 [NASA ADS] [CrossRef] [Google Scholar]
  14. Dormy, E., & Soward, A. M. 2007, Mathematical aspects of natural dynamos (CRC Press/Taylor & Francis) [Google Scholar]
  15. Duarte, L. D. V., Gastine, T., & Wicht, J. 2013, Phys. Earth Planet. Inter., 222, 22 [Google Scholar]
  16. Gastine, T., Duarte, L., & Wicht, J. 2012, A&A, 546, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Gilman, P. A., & Glatzmaier, G. A. 1981, ApJS, 45, 335 [NASA ADS] [CrossRef] [Google Scholar]
  18. Gissinger, C., Petitdemange, L., Schrinner, M., & Dormy, E. 2012, Phys. Rev. Lett., 108, 234501 [NASA ADS] [CrossRef] [Google Scholar]
  19. Gough, D. O. 1969, J. Atmos. Sci., 26, 448 [NASA ADS] [CrossRef] [Google Scholar]
  20. Hackman, T., Pelt, J., Mantere, M. J., et al. 2013, A&A, 553, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Jones, C. A. 2011, Ann. Rev. Fluid Mech., 43, 583 [NASA ADS] [CrossRef] [Google Scholar]
  22. Jones, C. A., Boronski, P., Brun, A. S., et al. 2011, Icarus, 216, 120 [NASA ADS] [CrossRef] [Google Scholar]
  23. Kochukhov, O., Mantere, M. J., Hackman, T., & Ilyin, I. 2013, A&A, 550, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Kutzner, C., & Christensen, U. R. 2002, Phys. Earth Planet. Inter., 131, 29 [Google Scholar]
  25. Lantz, S. R., & Fan, Y. 1999, ApJS, 121, 247 [Google Scholar]
  26. Morin, J., Donati, J.-F., Petit, P., et al. 2010, MNRAS, 407, 2269 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  27. Morin, J., Dormy, E., Schrinner, M., & Donati, J.-F. 2011, MNRAS, 418, L133 [CrossRef] [Google Scholar]
  28. Ogura, Y., & Phillips, N. A. 1962, J. Atmos. Sci., 19, 173 [NASA ADS] [CrossRef] [Google Scholar]
  29. Olson, P., & Christensen, U. R. 2006, Earth Planet. Sci. Lett., 250, 561 [NASA ADS] [CrossRef] [Google Scholar]
  30. Reinhold, T., Reiners, A., & Basri, G. 2013, A&A, 560, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Sasaki, Y., Takehiro, S.-i.,Kuramoto, K., & Hayashi, Y.-Y. 2011, Phys. Earth Planet. Inter., 188, 203 [NASA ADS] [CrossRef] [Google Scholar]
  32. Schrinner, M., Petitdemange, L., & Dormy, E. 2012, ApJ, 752, 121 [NASA ADS] [CrossRef] [Google Scholar]
  33. Schrinner, M., Petitdemange, L., Raynaud, R., & Dormy, E. 2014, A&A, 564, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Stelzer, Z., & Jackson, A. 2013, Geophys. J. Int., 193, 1265 [NASA ADS] [CrossRef] [Google Scholar]
  35. Yadav, R. K., Gastine, T., & Christensen, U. R. 2013a, Icarus, 225, 185 [NASA ADS] [CrossRef] [Google Scholar]
  36. Yadav, R. K., Gastine, T., Christensen, U. R., & Duarte, L. D. V. 2013b, ApJ, 774, 6 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table A.1

Summary of the coefficients obtained for the different scaling laws, with their standard error from the linear regression.

Table A.2

Comparison between the different scaling laws obtained with different dynamo models.

Table B.1

Overview of the simulations carried out with E = 10-4, Pr = 1, χ = 0.35, n = 2, and Nϱ = 0.1.

All Figures

thumbnail Fig. 1

Dipolar (black circles) and multipolar (white squares) dynamos as a function of Ra/Rac and Pm, for a central mass a) and a uniform mass distribution b). Crosses indicate the absence of a self-sustained dynamo.

In the text
thumbnail Fig. 2

a) Elsasser number Λ as a function of Ra/Rac, for Pm = 1 (green) and Pm = 3 (red). The meaning of the symbol shapes is defined in the caption of Fig. 1. A grey marker indicates that the solution loses its stability. b) Detail of the bifurcation close to the dynamo threshold for Pm = 1. Error bars represent the standard deviations.

In the text
thumbnail Fig. 3

Λdip/Rozmul\hbox{$\Lambda^\mathrm{dip}/\Roz^\mathrm{mul}$} as a function of Ra/Rac, for Pm = 1 (green diamonds), Pm = 2 (blue pentagons) and Pm = 3 (red stars). The point marked with the grey star has been computed with the model corresponding to the grey square in Fig. 2.

In the text
thumbnail Fig. 4

a) Relative dipole field strength fdip versus local Rossby number. b) Relative axial dipole field strength fdipax versus local Rossby number. The meaning of the symbol shapes is defined in the caption of Fig. 1.

In the text
thumbnail Fig. 5

a) Evolution of the modified tilt angle θ defined by Eq. (16)for multipolar dynamos as a function of Ra/Rac and Pm. Colour scale ranges from white (θ = 0) to black (θ = 0.7). b) θ as a function of Ra/Rac for Pm = 1. Upper x axis corresponds to the values of Ro for the multipolar branch. The meaning of the symbol shapes is defined in the caption of Fig. 1. Error bars represent the standard deviations.

In the text
thumbnail Fig. 6

Equatorial cross-section of the radial component of the velocity, with E = 10-4, Pr = 1. a)gr and Ra/Rac = 9.0 , Pm = 1. b)g ∝ 1/r2 and Ra/Rac = 9.9, Pm = 1.2.

In the text
thumbnail Fig. 7

Equatorial cross-section of the radial component of the magnetic field,with E = 10-4, Pr = 1. a) gr and Ra/Rac = 9.0 , Pm = 1. b), c)g ∝ 1/r2 and Ra/Rac = 9.9, Pm = 1.2. Colour in panel c) has been rescaled to highlight the emergence of a m = 1 mode at the outer sphere.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.