Quasi-geostrophic convection-driven dynamos in a thick spherical shell

We present dynamos computed using a hybrid QG-3D numerical scheme in a thick spherical shell geometry. Our model is based on a quasi-geostrophic convection code extended with a 3D treatment of heat transport and magnetic induction. We find a collection of self-sustained, multipolar, weak field dynamos with magnetic energy one or two orders of magnitude lower than the kinetic energy. The poloidal magnetic energy is weak and, by construction, there is a lack of equatorially anti-symmetric components in the Buoyancy and Lorentz forces. This leads to configurations where the velocity field is only weakly impacted by the magnetic field, similar to dynamos found in 3D simulations where zonal flows and the $\Omega$-effect dominate. The time-dependence of these dynamos is characterised by quasi-periodic oscillations that we attribute to dynamo waves. The QG-3D dynamos found so far are not Earth-like. The inability of our setup to produce strong, dipole-dominated, magnetic fields likely points to a missing ingredient in our QG flows, and a related lack of helicity and $\alpha$-effect. The models presented here may be more relevant for studying stellar dynamos where zonal flows are known to dominate. This study was carried out at modest control parameters, however moving to lower Ekman numbers, when smaller values of both the magnetic and hydrodynamic Prandtl numbers can be of interest, our approach will be able to gain in efficiency by using relatively coarse grids for the 3D magnetic and temperature fields and a finer grid for the QG velocity field.


INTRODUCTION
The magnetic field of the Earth is produced and sustained in Earth's outer core by turbulent motions of the liquid metal.These motions are widely believed to be driven by thermal and chemical convection and can be described by the equations of magnetohydrodynamics (hereafter MHD).The rotation of the Earth strongly influences the dynamics of the outer core as is evident from the extremely small values of its Ekman number, that characterises the ratio between the viscous and the Coriolis forces, Ek = ν/Ωd 2 ∼ 10 −15 -where ν is the kinematic viscosity, Ω is the Earth's rotation rate and d is the thickness of the outer core -and its Rossby number, which characterises the ratio between the inertia and the Coriolis forces, Ro = U/Ωd ∼ 10 −6 -where U is a typical velocity of the fluid.In addition to the inertia and the viscosity, which are dwarfed by the Coriolis force, the remaining forces in the system are the Buoyancy and the Lorentz force.The fact that Ek and Ro are so small in the outer core has lead to suggestions that it may be in a Quasi-Geostrophic (hereafter QG) dynamical balance at the 0 th order, i.e. a balance between the pressure and the Coriolis force (Davidson 2013;Calkins 2018).Other authors have argued that a Magnetostrophic balance could hold at leading order, i.e. an equi-librium between pressure, Coriolis and Lorentz forces (Roberts & Scott 1965;Dormy 2016); the relevant balance may depend on the length-scale (Aurnou & King 2017;Schwaiger et al. 2019) and part of the fluid volume (Schaeffer et al. 2017) considered, the QGbalance being dominant in most cases at large length-scales and outside of the tangent cylinder and of the boundary layers.
Numerical simulations of the primitive equations governing core dynamics have proven to be powerful tools for investigating dynamo mechanisms in a spherical shell geometry and their parameter dependencies (e.g., Glatzmaier & Roberts 1995;Christensen & Aubert 2006).Despite being restricted to a region of parameter space still remote from that of the Earth's core -Ek ≲ 10 −7 and Re ≲ 5 × 10 3 to be compared with Ek ∼ 10 −15 and Re ∼ 10 9 , where Re = U d/ν is the Reynolds number characterising the ratio between the inertia and the viscous force -progress has gradually been made in moving towards the Earth-like regime and in understanding the generic mechanisms at work (Schaeffer et al. 2017;Sheyko et al. 2018;Aubert et al. 2017;Schwaiger et al. 2019).In particular, following a specific path in parameter space, holding some parameters fixed at Earth-like values and gradually moving others toward the desired values, Aubert (2019Aubert ( , 2023) ) has been able to approach the conditions of Earth's core, by focusing on large arXiv:2312.09946v1[physics.geo-ph]15 Dec 2023 length-scales and employing hyperdiffusion.Such 3D simulations are nonetheless a major computational undertaking, especially for the most extreme parameters.
The possibility that the Coriolis force dominates and that core motions are columnar has motivated the development of reduced QG models to study core dynamics.These involve a 2D projection of the 3D MHD equations with the flow dynamics of the system being constrained to the equatorial plane (Busse 1970;Cardin & Olson 1994;Gillet & Jones 2006).Despite some limitations, especially if the temperature is also treated as 2D and contributions from the thermal wind are not included (Gillet & Jones 2006), this approach has proven effective at mimicking the behaviour of the full 3D system, at least in the purely hydrodynamic case (Aubert et al. 2003;Guervilly et al. 2019; Barrois et al. 2022).In the past few years, some interesting extensions of the QG method, mainly aiming at better accounting for the dynamics near the outer boundary, have been proposed (Labbé et al. 2015;Maffei et al. 2017;Jackson & Maffei 2020;Gerick et al. 2020).
A small number of studies have tried to include the effect of the magnetic field within a QG framework.QG-MHD models have been implemented in the context of mechanical forcing and only considering the largest scales of the Lorentz force (Schaeffer & Cardin 2006) or within a kinematic dynamo framework with convection including the effects of a 3D thermal wind (Guervilly 2010).Schaeffer & Cardin (2006) argued that a combination of flow time-dependence and the β-effect -due to the Coriolis force acting on fluid columns in spherical geometry -is sufficient to produce a QG dynamo independent of Ekman pumping.Both studies managed to obtain dynamos but found that the onset of the dynamo action as a function of the magnetic Reynolds number Rm = U d/λ, which characterises the ratio of magnetic diffusion over the convection overturn timescales -where λ is the magnetic diffusivity -was about five to ten times higher than that expected for 3D dynamos (e.g., Petitdemange 2018), with critical values of Rmc ∼ 500 in the QG compared with Rmc ∼ 50 in the full 3D case.
Several studies have investigated eigenmodes in QG-MHD systems, including the effect of the Lorentz force, considering small perturbations about an imposed background magnetic field.Canet et al. (2014) considered only the dynamics of axially invariant magnetic fields within a purely QG model, while more recently hybrid models have been developed, considering QG flow but a 3D magnetic field and projecting the 3D quantities on a QG basis (Gerick et al. 2021).Lately Jackson & Maffei (2020) have described a more complete 2D model QG-MHD based on quadratic magnetic quantities.
Using mean-field electrodynamics theory (Steenbeck et al. 1966;Krause & Radler 1980), it is possible to characterise dynamo action by considering azimuthally-averaged effects (see e.g.Schrinner et al. 2007).In this context, the terminology α-effect refers to the mean electrodynamics effect of helical flow generating poloidal magnetic energy from toroidal magnetic energy (or toroidal magnetic energy from poloidal energy) while Ω-effect refers to the production of toroidal from poloidal magnetic energy through an axisymmetric shear flow (Parker 1955;Moffatt 1978;Hollerbach 1996).Dynamo action in 3D convection-driven models of the geodynamo that produce strong dipolar fields is usually classified as being of α 2 type (Olson et al. 1999), at least when considering field generation by convective motions outside the inner core tangent cylinder.Inside the tangent cylinder the Ω-effect can also play a role, especially in strongly-driven cases (Schaeffer et al. 2017).In contrast, when strong zonal winds dominate con-vection outside the tangent cylinder the dynamo mechanism is typically found to be of αΩ type with the resulting poloidal magnetic fields being weak and multipolar (Schrinner et al. 2012).It seems there is a trade-off between strong zonal winds and strong dipolar magnetic fields.
QG models can efficiently simulate the dynamics of strong zonal flows (Schaeffer & Cardin 2005;Gastine 2019), so they might be expected to be relevant for studying dynamos where the Ω-effect is important, for instance in the context of stellar magnetic fields (e.g., Grote & Busse 2000;Goudard & Dormy 2008), or gas giants (e.g., Gastine et al. 2012).It is however less obvious whether or not QG dynamo models are relevant to terrestrial planets such as the Earth (e.g., Aubert et al. 2013).Schaeffer et al. (2016) have shown, within a kinematic dynamo framework, that adding magnetic pumping (an additional source of helicity related to the action of the Lorentz force, see Sreenivasan & Jones 2011) enables simple, observation-based, QG flows to generate dipole-dominated dynamos.The question of whether dynamically-consistent QG flow models, driven by convection and including feedback from the Lorentz force, can result in Earth-like dynamos is central to our study.
Our main objective here is to develop a hybrid QG-3D model based on QG convection in a thick spherical shell geometry (Gastine 2019), incorporating a 3D temperature field and thermal wind effects (Guervilly & Cardin 2017;Barrois et al. 2022), treating the magnetic field and its time evolution through the magnetic induction equation in 3D, and exploring the type of dynamos that are possible in this configuration.We compute the Lorentz force in 3D then z-average to obtain the impact on the QG flows.Building on previous QG-dynamos studies (Guervilly 2010;Schaeffer & Cardin 2006;Schaeffer et al. 2016), we present here an attempt to produce fully-resolved self-consistent convection-driven dynamos.
We describe our method and then the equations used in Section 2, present our main results in Section 3 and we conclude with a brief discussion and summary in Section 4. Tables with diagnostics and benchmarks of our method can be found in the Appendix section.

Hybrid QG-3D model formulation
Our hybrid QG-3D model builds on earlier work by Schaeffer & Cardin (2005); Gillet & Jones (2006); Guervilly & Cardin (2016) and Gastine (2019).We adopt the QG model formulation and notations of Barrois et al. (2022), and use the cylindrical coordinates system (s, ϕ, z) -with unit vectors (es, e ϕ , ez) -in a spherical shell between the inner and outer radii, si and so respectively, that rotates about the z-axis at a constant angular velocity Ω.We take η = si/so = 0.35 suitable for a thick shell such as the Earth's outer core.We solve the dimensionless equations of our problem under the Boussinesq approximation for the velocity field u, the magnetic field B and the temperature field T3D ≡ T cond 3D + ϑ3D.The last two fields are fully treated in 3D using the spherical coordinates system system (r, θ, ϕ3D), with unit vectors (er, e θ , e ϕ3D ).Both boundaries of the spherical shells are considered as electrically insulating, mechanically rigid, and we impose a fixed temperature contrast ∆T = Ti − To = T3D(ri) − T3D(ro) which drives convection.
In order to non-dimensionalise our variables, we use the shell thickness d = so − si as the reference length-scale, the viscous dif-fusion time d 2 /ν as the reference time-scale, the temperature contrast between the boundaries ∆T as the reference for temperature, and √ ρµ0λΩ -where ρ and λ are respectively the density and the magnetic diffusivity of the fluid and µ0 is the magnetic permeability of the vacuum -as the reference for the magnetic field.In such context, our system is controlled by four dimensionless parameters: the Ekman number, the Rayleigh number, the Prandtl number and the magnetic Prandtl number which are respectively defined by where αT is the thermal expansion coefficient, go = g(ro) is the gravity at the outer boundary, and κ is the thermal diffusivity.Note that the magnetic Prandtl number can equally be thought of as the ratio between the magnetic diffusion time τ λ and the viscous diffusion time τν , i.e.P m = τ λ /τν .We further assume that the dynamics is well described by the evolution of the axial vorticity averaged in the z direction ωz, such that the dynamics is restricted to that in the equatorial plane of the spherical shell (e.g., Maffei et al. 2017).Thus, the horizontal components of the velocity field u ⊥ , perpendicular to the rotation axis, are assumed to be mostly invariant along the rotation axis, i.e. u ⊥ ∼ (us, u ϕ , 0), where us and u ϕ are respectively the radial and azimuthal velocities.The axial velocity uz is considered as varying linearly with z in the direction of the rotation axis, including mass conservation at the spherical outer boundary, and the Ekman pumping contribution P (Schaeffer & Cardin 2005;Gastine 2019), yields where the Ekman pumping term P(Ek, u ⊥ , ωz) is deduced from Greenspan's formula (Greenspan et al. 1968) in a rigid sphere, i.e.
with β = 1 h dh ds = − s h 2 , and h ≡ √ s 2 o − s 2 , the half-height of a cylinder aligned with the rotation axis at a radius s.Note that the singularity of β at s = so is not an issue since mechanical boundary conditions enforce u = 0 there.
In this framework, the continuity equation from which it follows that there is a streamfunction ψ, which satisfies which accounts for the non-axisymmetric QG-velocity.u ϕ is the remaining axisymmetric zonal flow component, with the overbar x denoting the azimuthal average of any quantity x, i.e.
The dynamics of the non-axisymmetric motions are then described by the time-evolution of the z-averaged axial vorticity ωz ≡ ⟨(∇ × u) • ez⟩, where the angular brackets ⟨x⟩ refer to the axial average of any quantity x, such that ⟨x⟩ ≡ 1 2h The axial vorticity can be expressed in our framework as and its time evolution reads where the subscript ⊥ corresponds to the horizontal part of the operators and j ≡ ∇ × B.
Compared to the classical QG axial vorticity model, we have here followed the hybrid approach of Guervilly & Cardin (2016) and Barrois et al. (2022) and used the full 3D temperature and magnetic fields.The above equation ( 9) is thus coupled with the 3D temperature equation and with the 3D magnetic induction equation where u3D = (ur, u θ , u ϕ3D ) is the 3D-velocity in spherical coordinates.
In the above equations T cond 3D is the conducting temperature profile, a solution of ∇ 2 T cond 3D = 0.For a fixed temperature contrast between ri and ro without internal heating this takes the form We reconstruct the 3D-velocity field u3D from the QG velocity field u ⊥ using the conversion between cylindrical and spherical coordinate systems, where the cylindrical quantities are cast onto the 3D-grid using a bi-linear extrapolation (see Barrois et al. 2022, Appendix D for more details), such that where uz is obtained from Eq. (2), and with an additional contribution to the axisymmetric azimuthal motions u ϕ3D of the thermal wind, Tw(r, θ), which satisfies the relation and is integrated between the position z and the height of the column above the equator h, i.e.
where g(r) = r/ro is the dimensionless 3D gravity field.Note that inside the tangent cylinder, apart from the thermal wind contribution, the 3D velocitiy components are set to zero and thus mainly temperature or magnetic diffusion occurs in that region.
Finally, the z-averaged axial vorticity equation ( 9) has to be supplemented by an equation to account for the axisymmetric motions, that is where the last term of the right-hand-side is the Ekman-pumping contribution to the axisymmetric motions.

Computation of the QG-Lorentz Force
As seen in equations ( 9) and ( 16), our method requires that we compute the following two quantities related to the 3D Lorentz force averaged over axial direction We can expand Eq. ( 17) as and using the identity es = sin θ er + cos θ e θ this becomes Recalling that the s, ϕ and z components are orthogonal, that the z-averaging operator is linear so ⟨u + v⟩ = ⟨u⟩ + ⟨v⟩, and adopting the more compact notation f • ex ≡ fx for any field f and coordinate x, the Lorentz force terms become Regarding practical implementation, we find it useful to make use of the Leibniz's rule to switch the order of the s-derivative and z-integration steps in the first term of the curl in Eq. ( 21), where the surface term (j × B) ϕ (±h) cancels when B matches to a potential field at the outer boundary.This avoids the need to compute the s-derivative on the 3D physical grid using a finite-difference scheme, before the z-averaging step.Appendix B1 presents validation and benchmark tests that have been carried out to verify our computations of these QG-Lorentz force terms.

Numerics
The calculations presented in this study have been carried out using an extension of the open-source pseudo-spectral spherical QG code pizza (Gastine 2019; Barrois et al. 2022), written in Fortran and freely available at https://github.org/magic-sph/pizza/tree/hybrid_QG-3D under the GNU GPL v3 license.The 2D quantities are expanded in Fourier series up to the degree Nm in the azimuthal direction and in Chebyshev polynomials up to degree Ns in the radial direction.The 3D fields are expanded in Spherical Harmonics up to the degree and order ℓmax in the angular (θ, ϕ3D) directions and in Chebyshev polynomials with Nr collocation grid points in the radial direction.The open-source SHTns ⋆ library is employed to handle the Spherical Harmonic Transforms (Schaeffer 2013).Parallelisation of the hybrid QG-3D code relies on the Message Passing Interface (MPI) library.
The azimuthal (respectively radial) expansion involves adding zeros if m3D > m (Nr > Ns) and truncating fields at m = m3D = ℓ (Ns = Nr) if m3D < m (Nr < Ns).The same 3D and QG grids have been used when possible but because we explored cases with both P m < 1 and P r < 1, on some occasions we allowed the grid size to vary between fields.In these cases, nothing in particular has been done to bridge the two grids although we have used hyperdiffusion to mitigate this closure problem (Schaeffer 2005).

Hyperdiffusion
We included in our implementation an option to use hyperdiffusion.Following the formalism of Nataf & Schaeffer (2015) and Aubert et al. (2017), the diffusion operators entering Eqs. ( 9)-( 11)-( 16) are multiplied by hyperdiffusivity functions that solely depend on the azimuthal wavenumber m or on the spherical harmonic degree ℓ, such that on the velocity field, and on the magnetic field, where ℓH and mH are the cut-off degrees and azimuthal orders below which the hyperdiffusion has no effect.qH,u and qH,B are the strength of the hyperdiffusive effect on u ⊥ or B respectively and have been varied in the range 1.01 ≤ qH ≤ 1.08 (we do not apply any hyperdiffusion to the temperature field).
The values for mH and ℓH for the runs using hyperdiffusion on either u ⊥ or B are summarised in Table A1.
We have employed hyperdiffusion for two main reasons: (i) to mitigate closure issues when the 2D and 3D grids were different ; and (ii) to avoid the numerical problems in our most demanding runs -i.e. with the highest Rm -arising because of the tangent cylinder discontinuity and the approximations involved in the interpolation schemes.We have successfully removed the hyperdiffusion in several cases without observing any significant changes in the average properties.

Diagnostics
We now introduce the following notations for our various integral and average operators.For any quantity x, the hat x corresponds to its time average, that is where t0 is chosen such that any transient has been passed and with τ the averaging time interval long enough to reach a statistical equilibrium.The brackets {x}3D corresponds to the full spherical shell average, the brackets {x}S to a spherical surface average and the brackets {x}QG to the average over the equatorial annulus, such that, respectively with V3D the volume of the full spherical shell and VQG the volume outside of the tangent cylinder.
The dimensionless kinetic energy per unit volume Ekin, is defined by We similarly define the dimensionless magnetic energy per unit volume Emag, as The magnetic to kinetic energy ratio is then From these expressions, we define a diagnostic for the fluid velocity which characterises the average flow speed, based on the root-mean-square (r.m.s.) of the velocity, and which is denoted by the Reynolds number Then the magnetic Reynolds number is simply Rm = P m Re and we can finally define the Elsasser number, a non-dimensional measure of the magnetic field strength, which reads

RESULTS
We present here results of experiments conducted at control parameters of Ek = 3 × 10 −5 , P r = 10 −1 , varying P m, and focusing on a regime well above the onset of convection, i.e.Ra ⩾ 5 Rac.Time-integrating the nonlinear equations in absence of a magnetic field (Barrois et al. 2022), we estimated the critical Rayleigh number for this configuration to be Rac = 1.03 × 10 6 (which is thus the thermal convection critical value).In total 33 numerical simulations have been carried out, a list of their key diagnostics is given in Table A1 in Appendix A. Our first successful dynamo case -at parameters Ek = 3 × 10 −5 , P r = 0.1, P m = 0.9, Ra = 1.66 × 10 7 ∼ 16 Rac -was started from a motionless fluid with a strong axial dipole with Λ ∼ O(1) and a random perturbation in temperature.Subsequent experiments were initialized starting from a previously converged simulation.Attempts to restart the configurations with the largest Elsasser numbers from a strong axial dipole state were found to again yield the same final weak field multipolar solution.

Dynamo regime diagrams
A dynamo regime diagram as a function of P m and Ra/Rac for all our runs conducted at Ek = 3 × 10 −5 , P r = 0.1 is presented in Figure 1.The crosses correspond to simulations which failed to produce a self-sustained dynamo while the stars represent the growing dynamos.Similar to the case for 3D numerical dynamos we find that passing the onset for the dynamo action requires increasingly large P m on decreasing supercriticality.The shape of the dynamo threshold found in this diagram is qualitatively similar to that found by (Christensen & Aubert 2006, see their Fig. 1) despite the fact that they considered P r = 1 while we consider P r = 0.1.An important difference to note though is that we only find multipolar dynamos with the hybrid QG-3D formalism.
All the dynamos we have found so far have a low magnetic to kinetic energy ratio M < 1 and most of them have M < 10 −1 (with the exception of some of the points, e.g. at P m = 2.0, Ra/Rac ∼ 16.1 that reached a moderate M = 0.15 despite having the highest Λ = 8.02).They therefore fall into the weak-field dynamo regime, characterised by M ≪ 1 (e.g., Schaeffer et al. 2017;Aubert et al. 2017;Schwaiger et al. 2019).
Plotting the magnetic Reynolds number Rm against the supercriticality Ra/Rac in Figure 2, we observe that the minimum Rm required to obtain a self sustained dynamo in this setup is at least 500, in agreement with previous QG magnetic studies (e.g., Guervilly 2010).This value is about one order of magnitude higher than that found for 3D simulations, with for example, a critical value of Rmc ∼ 50 reported by Petitdemange (2018).We found that Λ > 0.05 for all the cases when a self-sustained dynamo was identified -with the highest value Λ ∼ 8 reached for the case at Ek = 3 × 10 −5 , P r = 0.1, P m = 2.0, Ra = 1.66 × 10 7 ∼ 16.1 Rac -and we generally observe an increase of Λ with increasing Rm.

An example weak field dynamo
To illustrate the typical features of our data set, we look at an example dynamo with parameters Ek = 3 × 10 −5 , P r = 0.1, P m = 0.9, Ra = 1.66×10 7 ∼ 16 Rac which was computed with a resolution of (Ns, Nm)/(Nr, ℓmax) = (385, 768)/(149, 148) for approximately ∼ 5τν = 4.5τ λ .This case is rather close to the onset of the dynamo action and decreasing P m by a factor 2 or decreasing Ra by a factor 1.5 was sufficient to lose dynamo action (see Fig. 1).The average magnetic Reynolds number of this simulation is Rm ≃ 900 and M = 0.04.
Figure 3 displays the time-averaged magnetic and kinetic energy spectra (a), a snapshot of the z-averaged axial vorticity ωz (b), a meridional section of the longitudinally-averaged azimuthal 3D velocity u ϕ3D (c), a snapshot of the radial magnetic field at the outer boundary of the dynamo region Br(ro) (d), and a meridional section of the ϕ-averaged azimuthal magnetic field B ϕ (e).
The power spectra shown in Fig. 3 (a) confirms that the dynamo is multipolar (the magnetic field is dominated by the degree ℓ = 2) and that the kinetic energy dominates over the magnetic energy at the largest length-scales by about 2 orders of magnitude (both for ℓ and m + 1 < 20, where m is here the QG azimuthal wavenumber).Additionally, we can observe that the velocity field is dominated by the azimuthal velocity at large length scales (m < 8) and that the magnetic field is strongly dominated by its toroidal part, this latter component being greater than the poloidal field by one order of magnitude at almost all degrees.This prevalence of the kinetic energy over the magnetic energy and of the toroidal component over the poloidal component of the magnetic field (dominated by the degree ℓ = 2) is typical of all our dynamos.
We can see in Fig. 3 (b-c) that the vorticity field (b) and the azimuthal velocity fields (c) are similar to what can be observed for a non-magnetic simulation (see, e.g.Fig. 3-4 in Barrois et al. 2022) with azimuthally elongated structures on the scale of the container in the case of the vorticity and with a strong zonal flow in the case of u ϕ3D .This suggests that the magnetic field does not strongly influence the velocity field, consistent with the low value of M for this case.
Turning to the magnetic field components, we can see in Fig. 3 (d-e) that Br(r = ro) (d) and B ϕ (e) have fairly low amplitude, and are mostly non-dipolar despite the clear equatorial antisymmetry.Br is dominated by small length-scales and the location of the strongest azimuthal magnetic field suggests that dynamo action mostly occurs at mid-latitudes and in the outer third of the shell.One could suspect that Ekman pumping is an important contributor to dynamo action, since at Ek = 3 × 10 −5 it is expected to have some impact close to the outer boundary.How-ever, although Ekman pumping has been shown to contribute to dynamo action close to the onset of dynamo (Busse 1975) it has been found in similar, but mechanically-forced, QG-models that removing the Ekman pumping flow does not significantly modify the dynamo onset (Schaeffer & Cardin 2006).We investigated the role of Ekman pumping in this dynamo by removing the Ekman pumping contribution to u3D -second part of Eq. ( 2) -used in the magnetic induction equation, and did not observe any major modifications in the resulting fields.We also conducted a test removing the thermal wind contribution to u ϕ3D -see Eq. ( 14) -and again did not observe any major change in the dynamo action, in agreement with Guervilly (2010) who also found that the thermal wind does not seem to have a strong impact on the dynamo onset.Our results are therefore consistent with QG flows influenced by a β-effect in spherical geometry -the first term of Eq. ( 2) -being sufficient to sustain a dynamo, in agreement with the earlier findings of Schaeffer & Cardin (2006) and Guervilly (2010).
The magnetic field structure reported in Fig. 3 is typical of our results.All our simulations display similar equatoriallyantisymmetric, predominantly degree 2, magnetic fields strongest at mid-to-high latitudes in the outer part of the shell.Note that the symmetry (or more correctly the anti-symmetry) of our magnetic field remains the same as that of our initial field, as there is, by construction, no ingredient to break the symmetry in our QG flows.

Dynamo mechanism
To illustrate the main mechanism underlying our dynamos, Figure 4 shows, from left to right respectively, meridional sections of the ϕ-averaged 3D-helicity H ≡ u3D • (∇ × u3D), the azimuthal component of the mean electromotive force (u ′ 3D × B ′ ) • e ϕ (hereafter EMF) and an estimate of the (rescaled) mean α-effect for the same dynamo presented in the previous section -where the prime x ′ denotes fluctuations about the azimuthal average of any quantity x, i.e. x ′ ≡ x − x and where τc is a typical convective turnover time (Brandenburg & Subramanian 2005) which can be approximated by τc ∼ d/{u ′ 3D }S following Brown et al. (2010) and Gastine et al. (2012).We find that the mean helicity, often thought to be a key ingredient in the magnetic induction via the so-called α-effect (e.g., Moffatt 1978;Jones 2008), changes sign between the North and South hemispheres and is stronger at mid latitudes and towards the outer boundary, where concentrations of B ϕ and strong zonal wind u ϕ3D have already been observed in the Fig. 3.Most of the 3D helicity is contained in its zcomponent (uz • (∇ × u3D) • ez) with only weak contributions from the sand ϕ-components (not shown).There is clearly a strong correlation between the mean EMF and the estimated α-effect despite the discrepancy in amplitudes (Fig. 4 b-c) with both quantities located where the maximum helicity and toroidal field are found, indicating that (i) our model produces an α-effect sufficient to support the dynamo action, (ii) kinetic helicity provides a reasonable mean-field approximation of the actual EMF.• e ϕ with superimposed poloidal field lines (bottom row).A strong correlation is found between both the Ω-effect and the toroidal field lines and between the EMF and the poloidal field lines, and both effects are concentrated in the same region localised in the upper part of the shell.The Ω-effect is stronger than the αeffect by 1-to-2 orders of magnitude and reaches its maxima in the region near to the outer boundary where both strong zonal flow and helicity are found.This is not unexpected as strong zonal winds sustained by Reynolds stresses -i.e. the correlations between the azimuthal motions u ′ ϕ3D and the cylindrically radial velocity u ′ sare expected to produce an important Ω-effect leading to multipolar dynamos dominated by the toroidal magnetic field.This has been found in a number of previous studies: for classical Boussinesq models (Sheyko et al. 2016), models with stress-free boundary conditions (Grote & Busse 2000;Goudard & Dormy 2008;Schrinner et al. 2012), models driven by strongly heterogeneous boundary heat flux (Dietrich et al. 2013) and anelastic dynamo models (Gastine et al. 2012).The main mechanism of our dynamos can thus be understood as powered by a strong Ω-effect -generated by the strong shear of of the zonal flow -with the toroidal field being converted into a poloidal field via an α-effect (Parker 1955).Such multipolar dynamos are often classified as αΩ or α 2 Ω type, following the mean field nomenclature of e.g.Steenbeck et al. (1966); Steenbeck & Krause (1969) or Malkus & Proctor (1975).

Dynamo waves
Figure 5 also shows that within a small fraction of a viscous diffusion time the toroidal and poloidal field lines have been mostly pushed towards the poles and replaced by field lines of the reversed polarity while deeper in the bulk the underlying Ω-effect and the EMF have changed sign.This process is reminiscent of dynamo waves that have been found in a 3D dynamo models with a range of different geometries, boundary conditions and driving mechanisms (see, e.g.Goudard & Dormy 2008;Schrinner et al. 2011;Simitev & Busse 2012;Sheyko et al. 2016) as well as in mechanically-forced QG dynamos (Schaeffer & Cardin 2006).
In Figure 6 we show the temporal evolution of B ϕ (θ) just beneath the outer boundary at r ∼ 0.92 ro -a so-called "butterflydiagram" -for a case at Ek = 3 × 10 −5 , P r = 0.1, P m = 0.9, Ra = 1.66 × 10 7 ∼ 16 Rac with a Rm of 909 (a) -same dynamo as presented in §3.2 -along with a more strongly driven case at Ek = 3×10 −5 , P r = 0.1, P m = 0.9, Ra = 3.0×10 7 ∼ 30 Rac with a higher Rm of 1284 (b).In both cases, we observe that the flux patches of the toroidal magnetic field appear symmetrically in both hemispheres at low-to-mid latitudes and move towards the poles until they reach the tangent cylinder which inhibits further motions.This is similar to what has been observed by e.g.Schrinner et al. ( 2011); Gastine et al. ( 2012) and Dietrich et al. (2013).We find that the cycle starts again after the flux patches have travelled all the way to the tangent cylinder and we can also see that some of the toroidal field travels equatorwards and quickly vanishes.The direction of migration (poleward or equatorward) of dynamo waves is controlled by gradients in the zonal flow (Yoshimura 1975).Waves are expected to travel polewards when the zonal flow increases with cylindrical radius (e.g., Sheyko et al. 2016), consistent with the zonal flow patterns in our dynamos (see Fig. 3 c).On increasing Rm between cases (a) and (b) the cycles become less coherent and less periodic with an oscillation period of 1 cycle per 4.5 × 10 −3 τ λ and 1 cycle per 2.7 × 10 −3 τ λ respectively.This is consistent with the expectation that the reversal frequency should increase with increasing forcing or increasing field strength; we have increased Ra/Rac leading to an increase in both Rm and the magnetic field strength -between cases (a) and (b).
In order to better characterise these oscillations, we test the approximate dispersion relation for (Parker) dynamo waves derived in Schrinner et al. (2011), which reads For the same case as in Fig. 3. lation periods extracted from our dynamos and those predicted for the Parker dynamo wave frequencies (above expression).We find an overall agreement similar to what has already been found in 3D dynamo studies (Busse & Simitev 2006;Schrinner et al. 2011;Gastine et al. 2012) even though some of the individual retrieved frequencies can be offset by up to a factor two compared to the theory.This is not unreasonable given all the approximations underlying the derivation of Eq. ( 32) so we can conclude that the observed oscillations in our dynamos can indeed be attributed to Parker dynamo waves.

Comparison with 3D dynamos
In order to compare our results with more conventional 3D simulations, we have computed a series of cases varying P m from P m = 0.5 to P m = 0.05 at Ek = 3 × 10 −5 , P r = 10 −1 , Ra = 2 × 10 7 ∼ 20 Rac using a full 3D method (MagIC, Wicht 2002).Here we compare these results with a hybrid QG-3D case at Ek = 3 × 10 −5 , P m = 0.5, P r = 10 −1 , Ra = 2 × 10 7 ∼ 20 Rac. Figure 8 displays results for one of these 3D simulations, computed at P m = 0.1 (left column) and for the hybrid QG-3D simulation computed at P m = 0.5 (right column).Snapshots of the z-averaged vorticity (panels a and e), meridional sections of the ϕ-averaged azimuthal 3D velocity (panels b and f), meridional sections of the ϕ-averaged kinetic helicity u3D • ∇ × u3D (panels c and g), and time series of the axisymmetric azimuthal magnetic field at r ∼ 0.92 ro (panels d and h) are presented for both the 3D and hybrid cases.
We stress that the comparison involves different P m for both setups.This is necessary because the onset for dynamo action happens at a lower P m with the 3D method (see e.g., Petitdemange 2018).Hence, the chosen parameters compare results at a similar level of supercriticality with respect to the dynamo onset.If one instead considers 3D and hybrid QG-3D results at the same control parameters, qualitative differences are found -the solution is already bi-stable in the 3D case with a strong dipolar field, while a multipolar solution is only found when starting from a very small seed magnetic field.
Despite the difference in P m, qualitatively similar solutions are then obtained when we are close to the onset of dynamo action for both configurations.Figure 8 (a-e) show convective patterns that are close to a non-magnetic case with large length-scale and azimuthally elongated z-averaged vorticity patterns.Both the 3D and hybrid QG-3D cases belong to the weak-field multipolarbranch dynamo with the kinetic energy that is greater than the magnetic energy, and the toroidal magnetic field that is much larger than the poloidal magnetic field.
Figure 8 (b-f) shows that the zonal flow patterns are similar in the two simulations, with a comparable level of geostrophy and similar radial gradients outside of the tangent cylinder, although the zonal flow is slightly stronger in the hybrid QG-3D case.The largest differences between the two cases can be observed in their respective helicities (Fig. 8 c-g), with a maximum helicity that is about one order of magnitude lower in the hybrid QG-3D case and a change of sign towards low latitudes at the outer boundary that is not visible in the 3D case.The spatial distribution of the helicity nevertheless remains comparable in the two cases with a segregation between mostly negative helicity in the northern hemisphere and mostly positive helicity in the southern hemisphere, similar to results previously reported for 3D simulations (see, for example Davidson & Ranjan 2018).
Finally, we can see in Fig. 8 (d-h) that both the 3D and hybrid QG-3D simulations display similar Parker-wave oscillations although the magnetic field is equatorially symmetric rather than antisymmetric and contains more small scale features in the 3D case.There is also much more magnetic energy in the 3D case with maxB 3D ϕ > 10 × maxB Hyb ϕ .A rough estimation of the Parker cy- • e ϕ with poloidal field lines superimposed (bottom row) during half a cycle of a dynamo wave (from a to e).Note that the time difference between two snapshots is not constant and that time is expressed in units of the viscous diffusion time.
cles gives an oscillation period of 1 cycle per 5 × 10 −3 τν for the 3D case and 1 cycle per 8.3 × 10 −3 τν for the hybrid QG-3D case, consistent with the fact that Rec and Rez are similar in the two simulations but H3D ≫ HHyb (see Eq. 32).
At P m = 0.5, Ra ∼ 20 Rac, we have computed the total r.m.s.kinetic helicity -as a proxy of the amplitude of the α-effect -for our weak field hybrid QG-3D case, as well as for a multipolar weak field 3D case, and a dipolar strong field 3D case (at these parameters the 3D run is bi-stable).Both 3D dynamos feature an average rms helicity roughly one order of magnitude larger than the hybrid QG-3D case (not shown).There is clearly a global lack of helicity in the hybrid QG-3D case compared with the 3D case, even if the spatial segregation of helicity remains as expected.
Overall, we conclude that although the comparison of 3D and hybrid QG-3D dynamos is not straightforward, because they give different solutions when run at the same control parameters, the hybrid QG-3D and the 3D methods do produce qualitatively similar results when considered at the same distance from the onset of dynamo action, and provided the 3D case is started from a weak seed.The major difference is that the hybrid QG-3D model involves much lower levels of kinetic heliticy.

DISCUSSION AND CONCLUSIONS
The results above demonstrate that it is possible to obtain dynamo action in a hybrid QG-3D magnetohydrodynamic model of rapidlyrotating convection in a spherical shell geometry, despite the strong assumptions described in section 2 concerning the velocity field.We have found several self-sustained multipolar dynamos in the parameter range Ek = 3×10 −5 , P r = 0.1−1, Ra ∼ (5−100) Rac and P m = 0.1 − 2.0 and have performed a detailed benchmarking of the Lorentz force in our setup as described in Appendix B.
Focusing on simulations conducted at Ek = 3 × 10 −5 , P r = 0.1 and a range of P m and Ra/Rac, we found QG dynamos characterised by a low magnetic to kinetic energy ratio M, a multipolar magnetic field, dominated by a toroidal field that is produced by an Ω-effect sustained by strong zonal winds and with only a relatively weak poloidal field.We have presented evidence for time-dependence in the form of dynamo waves in these solutions, similar to 3D models where the zonal flow plays an important role (e.g., Schrinner et al. 2011;Simitev & Busse 2012;Dietrich et al. 2013;Sheyko et al. 2016).A similar weak multipolar dynamo branch has been found in 3D models in a variety of setups (e.g., Christensen & Aubert 2006;Schrinner et al. 2012) and it seems that it is difficult to have both a strong dipolar magnetic field and a strong zonal flow (Gastine et al. 2012;Duarte et al. 2013).Although, we have not conducted a thorough analysis of the pa-Figure 6. Axisymmetric azimuthal magnetic field B ϕ (θ) at r ∼ 0.92 ro as a function of time -expressed in viscous units τν -(so called butterfly diagram) for (a) dynamo with control parameters Ek = 3 × 10 −5 , P r = 0.1, P m = 0.9, Ra = 1.66 × 10 7 and, (b) a second case, more strongly forced, with control parameters Ek = 3 × 10 −5 , P r = 0.1, P m = 0.9, Ra = 3.0 × 10 7 .In both panels, dashed lines mark the location of the tangent cylinder.
rameter space at P r = 1, for the cases tested we found similar results with weak and multipolar magnetic field solutions, suggesting that having P r = 0.1 rather than 1 is not crucial to the form of the dynamos reported here.A more detailed study focused on systematically varying the Prandtl number would be needed to fully quantify its role in the dynamo mechanism.
An initial motivation of this study was to test whether Earth-like dynamos could be achieved within a hybrid QG-3D convection-driven dynamo setup.We have found no examples of dipole-dominated dynamos.In contrast to 3D dynamos the poloidal field in our dynamos always remains weak compared with the toroidal field.This remains true, even at other values of P r (e.g.runs at P r = 1, not shown), which suggests that our model may lack some important ingredient in the field generation cycle that operates in dipole-dominated 3D dynamos.Schaeffer et al. (2016) previously found it was necessary to add an extra source of induction (magnetic pumping) in order to produce kinematic dynamos from observation-based QG flows.We suspect that a lack of αeffect associated with the lower level of helicity found in our dynamos compared to the 3D method is the main reason for their observed higher dynamo thresholds and possibly for the absence of a dipole-dominated branch in our configuration.On the contrary, we expect that if the zonal component of the flow was removed or damped, we would simply lose the dynamo action because of the loss of Ω-effect due to the strong zonal flow in the hybrid QG-3D set-up.The absence of equatorially antisymmetric axial flows and the associated missing correlations between uz and the temperature T3D have been shown to play a rather significant role in a lack of convective power already observed in the non-magnetic configuration investigated in (Barrois et al. 2022).These components are likely to also play a role in the helicity production (Ranjan et al. 2020).
A number of avenues can be envisaged for extending our model in order to enhance poloidal field generation.One obvious option is to add a simple α-effect term in the induction equation (e.g., Chan et al. 2001), but such simple functional forms are difficult to justify in terms of the underlying convection.Another option would be to follow Schaeffer et al. (2016) and implement a magnetic pumping whereby the velocity field is modified such that its helicity is enhanced based on an assumed toroidal magnetic field geometry (Sreenivasan & Jones 2011;Sreenivasan & Kar 2018).Slow MAC or MC waves might also be an important source of helicity for producing a dipolar dynamo (Varma & Sreenivasan 2022), but is unclear at the moment how best to parameterize their effect on the helicity, especially in the regions where the magnetic field is strong and heterogeneous.A final possibility could be to include a simple form of α-effect associated with helical waves propagating in the axial direction away from the equatorial plane where they are forced by turbulent convection.Davidson and co-workers have explored the hypothesis that such waves, forced by convection, can play a role in generating dipolar magnetic fields (for an overview see Davidson & Ranjan 2018).The axial averaging applied in our setup removes the dynamo action associated with such helical waves; Davidson & Ranjan (2015) set out how such an αeffect can be parameterized based on the kinetic energy of the flow in the equatorial plane.
An advantage of our hybrid QG-3D approach for the low magnetic Prandtl number regime of planetary cores is that it can treat the small scale velocity field efficiently within a QG framework while retaining a correct description of the 3D magnetic field and its boundary conditions.However further numerical work is needed before our model can be applied to this regime.So far, all our dynamos involved relatively weak Lorentz force and the energy of the magnetic field is much less than that associated with the velocity field.Moreover, in the present implementation, because of challenges associated with the tangent cylinder discontinuity and because of the rather crude interpolation schemes used to move between the QG and the 3D grids, a very large number of points was needed for accurate and stable computations.To take better advantage of the hybrid QG-3D approach with very different 3D and QG grid sizes it may be necessary to more carefully account for the action of the large length-scale magnetic field on the small length-scale velocity field, for example along the lines suggested by Schaeffer & Cardin (2006).
Returning to the geophysical context, we conclude that our hybrid QG-3D model seems incapable of producing Earth-like (strong-field, dipole-dominated) dynamos.This suggests, in agreement with the earlier findings of Schaeffer et al. (2016), that something important for geodynamo action is lost in moving between 3D flows and the simplified QG flows considered here.If hybrid QG-3D models are to be used to study the long-term behaviour of the geodynamo it will be necessary to find a principled scheme for parameterizing these missing effects, which may be related to structures in the axial flow component and their helicity.Hybrid QG-3D models could however already prove to be a valuable tool for studying the short-term behaviour of the geodynamo, on timescales shorter than the convective timescales when the dynamics is dominated by QG hydromagnetic waves (Aubert 2018;Aubert et al. 2022;Gillet et al. 2022).On these timescales, the dynamo-generated field can be considered steady and could be imposed, for example, based on results from a 3D simulation producing an Earth-like field (e.g., Aubert 2023).The hybrid QG-3D model is capable of efficiently representing both QG wave flows and related 3D magnetic field perturbations and has the potential to be significantly faster than full 3D simulations for studying such waves.APPENDIX A: RESULTS OF NUMERICAL SIMULATIONS Table A1: Summary of the hybrid QG-3D numerical simulations computed in this study at Ek = 3 × 10 −5 , P r = 0.1, using η = ri/ro = 0.35.Ra is the Rayleigh number (the supercriticality Sc = Ra/Rac is also provided, where Rac is the thermal convection critical value), P m is the magnetic Prandtl number, N u is the Nusselt number, Rm is the magnetic Reynolds number, Λ is the Elsasser number, M = Emag/ Ekin is the magnetic to kinetic energy ratio, ftor = Etor/ Emag is the toroidal to the total magnetic field energy ratio, ℓH and mH are the cut-off degree and azimuthal wavenumber below which the hyperdiffusion has no effect, and (Ns, Nm)/(Nr, ℓmax) are the grid-sizes of the run.The simulations with a ⋆ symbol in the second column are the growing dynamos we found.The determination of Eq. ( 21) relies on the computation of (j × B) which depends on ∇ × B and B. These quantities are computed using SHTns † functions, already validated in (Schaeffer 2013) and widely used in many codes.Given any fields (Bpol, Btor) in spectral space, and their radial-derivatives, the SHTns routines directly provide (Br, B θ , B ϕ ) and (jr, j θ , j ϕ ) on the 3D-physical-grid.The remaining part of the process involves the computation of the non-linear products (j × B), and the z-averaging of ∂s [s(j × B) ϕ ], (j × B) ϕ and (j × B)s which are computed on the physical grid.These quantities are then sent to the spectral space using fft functions where the ϕ-derivative is performed before assembling all the terms to obtain Eq. ( 21).Where the wide tildes x refers to the Fourier transform of any quantity x, that is the quantity x in the spectral space.

Ra
B1.2 Analytical benchmark from an artificial j × B field As steps 1) and 2) in (B.1) are respectively handled by SHTns routines and simply involve a linear product, we analytically validate our method from step 3).We need a field that is easily differentiable in s, easily integrated in z, periodic in ϕ, cancels at the outer boundary, but is not trivial.Following these constraints, we thus chose a completely artificial field, that reads

B1.3 Relative error definition
In order to discuss the validation and the accuracy of our numerical schemes, we define a relative error estimate, erel, as where the brackets in the above equation correspond to an average over the annulus as in Eq.( 26).

B2 Results
Figure A1 displays the convergence of the relative error for the Lorentz force terms as a function of the resolution.We can see that the computation of the Lorentz force term acting on the zonal flow FL , u ϕ (green curve) -that only involves a z-averaging of (j × B) ϕ -has an accuracy of order 4. Compared with Barrois et al. (2022), for all the z-integration steps we have used in this work a Simpson rule of integration, with an order 4 accuracy, and we satisfactorily retrieve the expected precision.On the other hand, we find a global accuracy of order 3 for the Lorentz force term acting on the vorticity FL , ωz (yellow curve) which involves additional operations (a ϕand s-derivatives).

B3 Conclusion
We have been able to retrieve the correct Lorentz force terms acting on both the z-averaged vorticity and the zonal-flow equations starting from an imposed artificial (j × B)-field and we have additionally verified that these results were consistent with an independent method (not shown).
The overall accuracy of the computation of the Lorentz force is limited by a number of interpolating schemes -i.e. a s-derivative scheme, a z-averaging scheme, and an ifft.We find an accuracy of order 4 converging toward an average relative error of 10 −16 (compared with an analytical solution) for the Lorentz force term acting on the zonal flow equation.And we find a global accuracy of order 3 converging toward an average relative error of 10 −6 (compared with an analytical solution) for the Lorentz force term acting on the vorticity equation.We thus consider the computation of Eq. ( 21) in our code to be validated.

Figure 1 .
Figure 1.Magnetic Prandtl number, P m, as a function of the supercriticality, Ra/Rac.Dynamo regime diagram computed for a series of cases at Ek = 3 × 10 −5 and P r = 0.1.Experiments that failed to produce a self-sustain dynamo are marked with a cross, those with a self-sustained multipolar dynamo are marked with a star and their symbol size is proportional to the Elsasser number, Λ.The different colors correspond to different P m.The dashed line marks the tentative limit between failed and growing dynamos.

Figure 2 .
Figure 2. Magnetic Reynolds Rm as a function of the supercriticality Ra/Rac.The symbols and the line have the same signification as in Fig. 1.

Further
insight into the dynamo mechanism at work in this dynamo is provided by Fig. 5 which displays a sequence of longitudinally-averaged snapshots of the Ω-effect sB • ∇ u ϕ3D s (e.g., Roberts & King 2013) with superimposed toroidal field lines (top row) and of the azimuthal component of the mean EMF (u ′ 3D × B ′ )

Figure 7
presents the resulting comparison between the oscil-

Figure 5 .
Figure 5. Temporal evolution of ϕ-averaged snapshots of the Ω-effect sB • ∇ u ϕ3D s with superimposed toroidal field lines (top row) and of the azimuthal component of the mean EMF (u ′ 3D × B ′ )• e ϕ with poloidal field lines superimposed (bottom row) during half a cycle of a dynamo wave (from a to e).Note that the time difference between two snapshots is not constant and that time is expressed in units of the viscous diffusion time.

Figure 7 .
Figure 7. Measured frequencies of oscillatory dynamos as a function of the frequency predicted by the Parker dynamo wave dispersion relation (32), for all of our dynamos.Different colors correspond to different magnetic Prandtl numbers and the dashed line corresponds to the expected relation (32).

Figure A1 .
Figure A1.Convergence of the relative error of the Lorentz force terms tested against an analytical solution.