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

O. Barrois, T. Gastine, C. C. Finlay

► To cite this version:

O. Barrois, T. Gastine, C. C. Finlay. Quasi-geostrophic convection-driven dynamos in a thick spherical shell. Geophysical Journal International, 2024, 236, pp.1373-1389. 10.1093/gji/ggad478. insu-04462197

HAL Id: insu-04462197

https://insu.hal.science/insu-04462197

Submitted on 16 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Quasi-geostrophic convection-driven dynamos in a thick spherical shell

O. Barrois,<fnref n="1" /> T. Gastine<fnref n="1" />,<fnref n="2" /> and C.C. Finlay<fnref n="1" />,<fnref n="1" />

Accepted 2023 December 14. Received 2023 December 11; in original form 2023 July 3

## SUMMARY

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 3-D 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 antisymmetric 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 3-D simulations where zonal flows and the Ω-effect dominate. The time-dependence of these dynamos is characterized 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 α-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 3-D magnetic and temperature fields and a finer grid for the QG velocity field.

**Key words:** Dynamo: theories and simulations; Numerical modelling; Earth core.

## 1 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 characterizes the ratio between the viscous and the Coriolis forces, $Ek = v/\Omega d^2 \sim 10^{-15}$, where v is the kinematic viscosity, Ω is the Earth's rotation rate and d is the thickness of the outer core, and its Rossby number, which characterizes the ratio between the inertia and the Coriolis forces, $Ro = U/\Omega d \sim 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 0th order, that is 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, that is an equilibrium 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 QG-balance 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 \lesssim 10^{-7}$ and $Re \lesssim 5 \times 10^3$ to be compared with $Ek \sim 10^{-15}$ and $Re \sim 10^9$, where $Re = Ud/v$ 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; Aubert et al. 2017; Sheyko et al. 2018; 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 towards the

<footnote n="1">¹Division of Geomagnetism, DTU Space, Technical University of Denmark, Lyngby DK-2800, Denmark. E-mail: oli.barrois@gmail.com, ovrbr@dtu.dk</footnote>

<footnote n="2">²Institut de Physique du Globe de Paris, Université Paris Cité, CNRS, F-75005 Paris, France</footnote>
# Quasi-geostrophic convection-driven dynamos in a thick spherical shell

O. Barrois,<fnref n="1" /> T. Gastine<fnref n="1" />,<fnref n="2" /> and C.C. Finlay<fnref n="1" />,<fnref n="1" />

Accepted 2023 December 14. Received 2023 December 11; in original form 2023 July 3

## SUMMARY

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 3-D 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 antisymmetric 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 3-D simulations where zonal flows and the Ω-effect dominate. The time-dependence of these dynamos is characterized 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 α-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 3-D magnetic and temperature fields and a finer grid for the QG velocity field.

**Key words:** Dynamo: theories and simulations; Numerical modelling; Earth core.

## 1 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 characterizes the ratio between the viscous and the Coriolis forces, $Ek = v/\Omega d^2 \sim 10^{-15}$, where v is the kinematic viscosity, Ω is the Earth's rotation rate and d is the thickness of the outer core, and its Rossby number, which characterizes the ratio between the inertia and the Coriolis forces, $Ro = U/\Omega d \sim 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 0th order, that is 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, that is an equilibrium 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 QG-balance 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 \lesssim 10^{-7}$ and $Re \lesssim 5 \times 10^3$ to be compared with $Ek \sim 10^{-15}$ and $Re \sim 10^9$, where $Re = Ud/v$ 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; Aubert et al. 2017; Sheyko et al. 2018; 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 towards the

<footnote n="1">¹Division of Geomagnetism, DTU Space, Technical University of Denmark, Lyngby DK-2800, Denmark. E-mail: oli.barrois@gmail.com, ovrbr@dtu.dk</footnote>

<footnote n="2">²Institut de Physique du Globe de Paris, Université Paris Cité, CNRS, F-75005 Paris, France</footnote>
desired values, Aubert (2019, 2023) has been able to approach the conditions of Earth’s core, by focusing on large length-scales and employing hyperdiffusion. Such 3-D 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 2-D projection of the 3-D 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 2-D and contributions from the thermal wind are not included (Gillet & Jones 2006), this approach has proven effective at mimicking the behaviour of the full 3-D 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 3-D 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 = Ud/λ$, which characterizes the ratio of magnetic diffusion over the convection overturn timescales, where $λ$ is the magnetic diffusivity—was about 5–10 times higher than that expected for 3-D dynamos (e.g. Petitdemange 2018), with critical values of $Rm_c \sim 500$ in the QG compared with $Rm_c \sim 50$ in the full 3-D 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 3-D magnetic field and projecting the 3-D quantities on a QG basis (Gerick et al. 2021). Lately Jackson & Maffei (2020) have described a more complete 2-D model QG-MHD based on quadratic magnetic quantities.

Using mean-field electrodynamics theory (Steenbeck et al. 1966; Krause & Radler 1980), it is possible to characterize 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 3-D convection-driven models of the geodynamo that produce strong dipolar fields is usually classified as being of α² 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 convection 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 3-D 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 3-D, and exploring the type of dynamos that are possible in this configuration. We compute the Lorentz force in 3-D then z-average to obtain the impact on the QG flows. Building on previous QG-dynamos studies (Schaeffer & Cardin 2006; Guervilly 2010; 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.

# 2 METHODOLOGY

## 2.1 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, \phi, z$)—with unit vectors ($e_s, e_\phi, e_z$)—in a spherical shell between the inner and outer radii, $s_i$ and $s_o$ respectively, that rotates about the z-axis at a constant angular velocity $\Omega$. We take $\eta = s_i/s_o = 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 $\mathbf{u}$, the magnetic field $\mathbf{B}$ and the temperature field $T_{3D} \equiv T_{3D}^{\text{cond}} + \vartheta_{3D}$. The last two fields are fully treated in 3-D using the spherical coordinates system ($r, \theta, \phi_{3D}$), with unit vectors ($e_r, e_\theta, e_{\phi_{3D}}$). Both boundaries of the spherical shells are considered as electrically insulating, mechanically rigid and we impose a fixed temperature contrast $\Delta T = T_i - T_o = T_{3D}(r_i) - T_{3D}(r_o)$ which drives convection.

In order to non-dimensionalize our variables, we use the shell thickness $d = s_o - s_i$ as the reference length-scale, the viscous diffusion time $d^2/\nu$ as the reference timescale, the temperature contrast between the boundaries $\Delta T$ as the reference for temperature, and $\sqrt{\rho\mu_0\lambda\Omega}$—where $\rho$ and $\lambda$ are, respectively, the density and the
desired values, Aubert (2019, 2023) has been able to approach the conditions of Earth’s core, by focusing on large length-scales and employing hyperdiffusion. Such 3-D 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 2-D projection of the 3-D 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 2-D and contributions from the thermal wind are not included (Gillet & Jones 2006), this approach has proven effective at mimicking the behaviour of the full 3-D 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 3-D 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 = Ud/λ$, which characterizes the ratio of magnetic diffusion over the convection overturn timescales, where $λ$ is the magnetic diffusivity—was about 5–10 times higher than that expected for 3-D dynamos (e.g. Petitdemange 2018), with critical values of $Rm_c \sim 500$ in the QG compared with $Rm_c \sim 50$ in the full 3-D 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 3-D magnetic field and projecting the 3-D quantities on a QG basis (Gerick et al. 2021). Lately Jackson & Maffei (2020) have described a more complete 2-D model QG-MHD based on quadratic magnetic quantities.

Using mean-field electrodynamics theory (Steenbeck et al. 1966; Krause & Radler 1980), it is possible to characterize 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 3-D convection-driven models of the geodynamo that produce strong dipolar fields is usually classified as being of α² 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 convection 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 3-D 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 3-D, and exploring the type of dynamos that are possible in this configuration. We compute the Lorentz force in 3-D then z-average to obtain the impact on the QG flows. Building on previous QG-dynamos studies (Schaeffer & Cardin 2006; Guervilly 2010; 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.

# 2 METHODOLOGY

## 2.1 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, \phi, z$)—with unit vectors ($e_s, e_\phi, e_z$)—in a spherical shell between the inner and outer radii, $s_i$ and $s_o$ respectively, that rotates about the z-axis at a constant angular velocity $\Omega$. We take $\eta = s_i/s_o = 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 $\mathbf{u}$, the magnetic field $\mathbf{B}$ and the temperature field $T_{3D} \equiv T_{3D}^{\text{cond}} + \vartheta_{3D}$. The last two fields are fully treated in 3-D using the spherical coordinates system ($r, \theta, \phi_{3D}$), with unit vectors ($e_r, e_\theta, e_{\phi_{3D}}$). Both boundaries of the spherical shells are considered as electrically insulating, mechanically rigid and we impose a fixed temperature contrast $\Delta T = T_i - T_o = T_{3D}(r_i) - T_{3D}(r_o)$ which drives convection.

In order to non-dimensionalize our variables, we use the shell thickness $d = s_o - s_i$ as the reference length-scale, the viscous diffusion time $d^2/\nu$ as the reference timescale, the temperature contrast between the boundaries $\Delta T$ as the reference for temperature, and $\sqrt{\rho\mu_0\lambda\Omega}$—where $\rho$ and $\lambda$ are, respectively, the density and the
magnetic diffusivity of the fluid and $\mu_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

$$Ek = \frac{\nu}{\Omega d^2}, \quad Ra = \frac{\alpha_T g_o \Delta T d^3}{\kappa \nu}, \quad Pr = \frac{\nu}{\kappa}, \quad Pm = \frac{\nu}{\lambda}, \quad (1)$$

where $\alpha_T$ is the thermal expansion coefficient, $g_o = g(r_o)$ is the gravity at the outer boundary, and $\kappa$ is the thermal diffusivity. Note that the magnetic Prandtl number can equally be thought of as the ratio between the magnetic diffusion time $\tau_\lambda$ and the viscous diffusion time $\tau_v$, that is $Pm = \tau_\lambda/\tau_v$.

We further assume that the dynamics is well described by the evolution of the axial vorticity averaged in the z direction $\omega_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 $\mathbf{u}_\perp$, perpendicular to the rotation axis, are assumed to be mostly invariant along the rotation axis, that is $\mathbf{u}_\perp \sim (u_s, u_\phi, 0)$, where $u_s$ and $u_\phi$ are, respectively, the radial and azimuthal velocities. The axial velocity $u_z$ 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 $\mathcal{P}$ (Schaeffer & Cardin 2005; Gastine 2019), yields

$$u_z(s, \phi, z) = z \left[ \beta u_s + \frac{Ek}{2} \mathcal{P}(Ek, \mathbf{u}_\perp, \omega_z) \right], \quad (2)$$

where the Ekman pumping term $\mathcal{P}(Ek, \mathbf{u}_\perp, \omega_z)$ is deduced from Greenspan's formula (Greenspan et al. 1968) in a rigid sphere, that is

$$\mathcal{P}(Ek, \mathbf{u}_\perp, \omega_z) = - \left( \frac{s_o}{Ek} \right)^{1/2} \frac{1}{h^{3/2}} \left[ \omega_z - \frac{\beta}{2} u_\phi + \beta \frac{\partial u_s}{\partial \phi} - \frac{5s_o}{2h} u_s \right], \quad (3)$$

with $\beta = \frac{1}{h} \frac{dh}{ds} = -\frac{s}{h^2}$, and $h \equiv \sqrt{s_o^2 - s^2}$, the half-height of a cylinder aligned with the rotation axis at a radius $s$. Note that the singularity of $\beta$ at $s = s_o$ is not an issue since mechanical boundary conditions enforce $\mathbf{u} = \mathbf{0}$ there.

In this framework, the continuity equation $\nabla \cdot \mathbf{u} = 0$ reads

$$\frac{\partial(su_s)}{\partial s} + \frac{\partial u_\phi}{\partial \phi} + \beta s u_s = 0, \quad (4)$$

from which it follows that there is a stream function $\psi$, which satisfies

$$u_s = \frac{1}{s} \frac{\partial \psi}{\partial \phi}, \quad u_\phi = \overline{u_\phi} - \frac{\partial \psi}{\partial s} - \beta \psi, \quad (5)$$

which accounts for the non-axisymmetric QG-velocity. $\overline{u_\phi}$ is the remaining axisymmetric zonal flow component, with the overbar $\bar{x}$ denoting the azimuthal average of any quantity $x$, that is

$$\bar{x} = \frac{1}{2\pi} \int_{0}^{2\pi} x \, d\phi. \quad (6)$$

The dynamics of the non-axisymmetric motions are then described by the time-evolution of the z-averaged axial vorticity $\omega_z \equiv \langle (\nabla \times \mathbf{u}) \cdot e_z \rangle$, where the angular brackets $\langle x \rangle$ refer to the axial average of any quantity $x$, such that

$$\langle x \rangle = \frac{1}{2h} \int_{-h}^{h} x \, dz. \quad (7)$$

The axial vorticity can be expressed in our framework as

$$\omega_z = \frac{1}{s} \frac{\partial(s\overline{u_\phi})}{\partial s} - \nabla^2 \psi - \frac{1}{s} \frac{\partial(\beta s \psi)}{\partial s}, \quad (8)$$

and its time evolution reads

$$\begin{aligned} \frac{\partial \omega_z}{\partial t} + \nabla_{\perp} \cdot (\mathbf{u}_{\perp} \omega_z) - \frac{2}{Ek} \beta u_s = \nabla_{\perp}^2 \omega_z - \frac{Ra}{Pr} \left\langle \frac{1}{r_o} \frac{\partial \vartheta_{3D}}{\partial \phi_{3D}} \right\rangle \\ + \frac{1}{Ek Pm} \left\langle \nabla \times (\mathbf{j} \times \mathbf{B}) \cdot \mathbf{e}_z \right\rangle \\ + \mathcal{P}(Ek, \mathbf{u}_{\perp}, \omega_z), \end{aligned} \quad (9)$$

where the subscript $\perp$ corresponds to the horizontal part of the operators and $\mathbf{j} \equiv \nabla \times \mathbf{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 3-D temperature and magnetic fields. The above eq. (9) is thus coupled with the 3-D temperature equation

$$\frac{\partial \vartheta_{3D}}{\partial t} + \mathbf{u}_{3D} \cdot \nabla \vartheta_{3D} + u_r \frac{dT_{3D}^{\text{cond}}}{dr} = \frac{1}{Pr} \nabla^2 \vartheta_{3D}, \quad (10)$$

and with the 3-D magnetic induction equation

$$\frac{\partial \mathbf{B}}{\partial t} = \nabla \times (\mathbf{u}_{3D} \times \mathbf{B}) + \frac{1}{Pm} \nabla^2 \mathbf{B}, \quad (11)$$

where $\mathbf{u}_{3D} = (u_r, u_\theta, u_{\phi 3D})$ is the 3-D velocity in spherical coordinates.

In the above equations $T_{3D}^{\text{cond}}$ is the conducting temperature profile, a solution of $\nabla^2 T_{3D}^{\text{cond}} = 0$. For a fixed temperature contrast between $r_i$ and $r_o$ without internal heating this takes the form

$$T_{3D}^{\text{cond}}(r) = \frac{r_o r_i}{r} - r_i, \quad \frac{dT_{3D}^{\text{cond}}}{dr} = - \frac{r_i r_o}{r^2}. \quad (12)$$

We reconstruct the 3-D velocity field $\mathbf{u}_{3D}$ from the QG velocity field $\mathbf{u}_\perp$ using the conversion between cylindrical and spherical coordinate systems, where the cylindrical quantities are cast onto the 3-D grid using a bi-linear extrapolation (see Barrois et al. 2022, Appendix D for more details), such that

$$\left\{ \begin{aligned} u_r(r, \theta, \phi_{3D}) &= \sin \theta u_s(s, \phi) + \cos \theta u_z(s, \phi, z), \\ u_\theta(r, \theta, \phi_{3D}) &= \cos \theta u_s(s, \phi) - \sin \theta u_z(s, \phi, z), \\ u_{\phi 3D}(r, \theta, \phi_{3D}) &= u_\phi(s, \phi) + T_w(r, \theta), \end{aligned} \right. \quad (13)$$

where $u_z$ is obtained from eq. (2), and with an additional contribution to the axisymmetric azimuthal motions $\overline{u_{\phi 3D}}$ of the thermal wind, $T_w(r, \theta)$, which satisfies the relation

$$\frac{\partial \overline{u_{\phi 3D}}}{\partial z} = \frac{Ra Ek}{2Pr} \frac{g(r)}{r} \frac{\partial \overline{\vartheta_{3D}}}{\partial \theta}, \quad (14)$$

and is integrated between the position $z$ and the height of the column above the equator $h$, that is

$$T_w(r, \theta) = \frac{RaEk}{2Pr} \int_h^z \frac{g(r)}{r} \frac{\partial \overline{\vartheta_{3D}}}{\partial \theta} dz', \quad (15)$$

where $g(r) = r/r_o$ is the dimensionless 3-D gravity field. Note that inside the tangent cylinder, apart from the thermal wind contribution, the 3-D velocity components are set to zero and thus mainly temperature or magnetic diffusion occurs in that region.

Finally, the z-averaged axial vorticity eq. (9) has to be supplemented by an equation to account for the axisymmetric motions, that is

$$\begin{aligned} \frac{\partial \overline{u_\phi}}{\partial t} + \overline{u_s} \frac{\partial u_\phi}{\partial s} + \frac{\overline{u_s u_\phi}}{s} = \nabla_{\perp}^2 \overline{u_\phi} - \frac{1}{s^2} \overline{u_\phi} \\ + \frac{1}{Ek Pm} \left\langle (\mathbf{j} \times \mathbf{B}) \cdot \mathbf{e}_\phi \right\rangle \\ - \left( \frac{s_o}{Ek} \right)^{1/2} \frac{1}{h^{3/2}} \overline{u_\phi}, \end{aligned} \quad (16)$$
magnetic diffusivity of the fluid and $\mu_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

$$Ek = \frac{\nu}{\Omega d^2}, \quad Ra = \frac{\alpha_T g_o \Delta T d^3}{\kappa \nu}, \quad Pr = \frac{\nu}{\kappa}, \quad Pm = \frac{\nu}{\lambda}, \quad (1)$$

where $\alpha_T$ is the thermal expansion coefficient, $g_o = g(r_o)$ is the gravity at the outer boundary, and $\kappa$ is the thermal diffusivity. Note that the magnetic Prandtl number can equally be thought of as the ratio between the magnetic diffusion time $\tau_\lambda$ and the viscous diffusion time $\tau_v$, that is $Pm = \tau_\lambda/\tau_v$.

We further assume that the dynamics is well described by the evolution of the axial vorticity averaged in the z direction $\omega_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 $\mathbf{u}_\perp$, perpendicular to the rotation axis, are assumed to be mostly invariant along the rotation axis, that is $\mathbf{u}_\perp \sim (u_s, u_\phi, 0)$, where $u_s$ and $u_\phi$ are, respectively, the radial and azimuthal velocities. The axial velocity $u_z$ 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 $\mathcal{P}$ (Schaeffer & Cardin 2005; Gastine 2019), yields

$$u_z(s, \phi, z) = z \left[ \beta u_s + \frac{Ek}{2} \mathcal{P}(Ek, \mathbf{u}_\perp, \omega_z) \right], \quad (2)$$

where the Ekman pumping term $\mathcal{P}(Ek, \mathbf{u}_\perp, \omega_z)$ is deduced from Greenspan's formula (Greenspan et al. 1968) in a rigid sphere, that is

$$\mathcal{P}(Ek, \mathbf{u}_\perp, \omega_z) = - \left( \frac{s_o}{Ek} \right)^{1/2} \frac{1}{h^{3/2}} \left[ \omega_z - \frac{\beta}{2} u_\phi + \beta \frac{\partial u_s}{\partial \phi} - \frac{5s_o}{2h} u_s \right], \quad (3)$$

with $\beta = \frac{1}{h} \frac{dh}{ds} = -\frac{s}{h^2}$, and $h = \sqrt{s_o^2 - s^2}$, the half-height of a cylinder aligned with the rotation axis at a radius $s$. Note that the singularity of $\beta$ at $s = s_o$ is not an issue since mechanical boundary conditions enforce $\mathbf{u} = \mathbf{0}$ there.

In this framework, the continuity equation $\nabla \cdot \mathbf{u} = 0$ reads

$$\frac{\partial(su_s)}{\partial s} + \frac{\partial u_\phi}{\partial \phi} + \beta s u_s = 0, \quad (4)$$

from which it follows that there is a stream function $\psi$, which satisfies

$$u_s = \frac{1}{s} \frac{\partial \psi}{\partial \phi}, \quad u_\phi = \overline{u_\phi} - \frac{\partial \psi}{\partial s} - \beta \psi, \quad (5)$$

which accounts for the non-axisymmetric QG-velocity. $\overline{u_\phi}$ is the remaining axisymmetric zonal flow component, with the overbar $\bar{x}$ denoting the azimuthal average of any quantity $x$, that is

$$\bar{x} = \frac{1}{2\pi} \int_{0}^{2\pi} x \, d\phi. \quad (6)$$

The dynamics of the non-axisymmetric motions are then described by the time-evolution of the z-averaged axial vorticity $\omega_z \equiv \langle (\nabla \times \mathbf{u}) \cdot e_z \rangle$, where the angular brackets $\langle x \rangle$ refer to the axial average of any quantity $x$, such that

$$\langle x \rangle = \frac{1}{2h} \int_{-h}^{h} x \, dz. \quad (7)$$

The axial vorticity can be expressed in our framework as

$$\omega_z = \frac{1}{s} \frac{\partial(s \overline{u_\phi})}{\partial s} - \nabla^2 \psi - \frac{1}{s} \frac{\partial(\beta s \psi)}{\partial s}, \quad (8)$$

and its time evolution reads

$$\begin{aligned} \frac{\partial \omega_z}{\partial t} + \nabla_{\perp} \cdot (\mathbf{u}_{\perp} \omega_z) - \frac{2}{Ek} \beta u_s = & \nabla_{\perp}^2 \omega_z - \frac{Ra}{Pr} \left\langle \frac{1}{r_o} \frac{\partial \vartheta_{3D}}{\partial \phi_{3D}} \right\rangle \\ & + \frac{1}{Ek Pm} \left\langle \nabla \times (\mathbf{j} \times \mathbf{B}) \cdot \mathbf{e}_z \right\rangle \\ & + \mathcal{P}(Ek, \mathbf{u}_{\perp}, \omega_z), \end{aligned} \quad (9)$$

where the subscript $\perp$ corresponds to the horizontal part of the operators and $\mathbf{j} \equiv \nabla \times \mathbf{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 3-D temperature and magnetic fields. The above eq. (9) is thus coupled with the 3-D temperature equation

$$\frac{\partial \vartheta_{3D}}{\partial t} + \mathbf{u}_{3D} \cdot \nabla \vartheta_{3D} + u_r \frac{dT_{3D}^{\text{cond}}}{dr} = \frac{1}{Pr} \nabla^2 \vartheta_{3D}, \quad (10)$$

and with the 3-D magnetic induction equation

$$\frac{\partial \mathbf{B}}{\partial t} = \nabla \times (\mathbf{u}_{3D} \times \mathbf{B}) + \frac{1}{Pm} \nabla^2 \mathbf{B}, \quad (11)$$

where $\mathbf{u}_{3D} = (u_r, u_\theta, u_{\phi 3D})$ is the 3-D velocity in spherical coordinates.

In the above equations $T_{3D}^{\text{cond}}$ is the conducting temperature profile, a solution of $\nabla^2 T_{3D}^{\text{cond}} = 0$. For a fixed temperature contrast between $r_i$ and $r_o$ without internal heating this takes the form

$$T_{3D}^{\text{cond}}(r) = \frac{r_o r_i}{r} - r_i, \quad \frac{dT_{3D}^{\text{cond}}}{dr} = - \frac{r_i r_o}{r^2}. \quad (12)$$

We reconstruct the 3-D velocity field $\mathbf{u}_{3D}$ from the QG velocity field $\mathbf{u}_\perp$ using the conversion between cylindrical and spherical coordinate systems, where the cylindrical quantities are cast onto the 3-D grid using a bi-linear extrapolation (see Barrois et al. 2022, Appendix D for more details), such that

$$\left\{ \begin{array}{l} u_r(r, \theta, \phi_{3D}) = \sin \theta u_s(s, \phi) + \cos \theta u_z(s, \phi, z), \\ u_\theta(r, \theta, \phi_{3D}) = \cos \theta u_s(s, \phi) - \sin \theta u_z(s, \phi, z), \\ u_{\phi 3D}(r, \theta, \phi_{3D}) = u_\phi(s, \phi) + T_w(r, \theta), \end{array} \right. \quad (13)$$

where $u_z$ is obtained from eq. (2), and with an additional contribution to the axisymmetric azimuthal motions $\overline{u_{\phi 3D}}$ of the thermal wind, $T_w(r, \theta)$, which satisfies the relation

$$\frac{\partial \overline{u_{\phi 3D}}}{\partial z} = \frac{Ra Ek}{2Pr} \frac{g(r)}{r} \frac{\partial \overline{\vartheta_{3D}}}{\partial \theta}, \quad (14)$$

and is integrated between the position $z$ and the height of the column above the equator $h$, that is

$$T_w(r, \theta) = \frac{RaEk}{2Pr} \int_h^z \frac{g(r)}{r} \frac{\partial \overline{\vartheta_{3D}}}{\partial \theta} dz', \quad (15)$$

where $g(r) = r/r_o$ is the dimensionless 3-D gravity field. Note that inside the tangent cylinder, apart from the thermal wind contribution, the 3-D velocity components are set to zero and thus mainly temperature or magnetic diffusion occurs in that region.

Finally, the z-averaged axial vorticity eq. (9) has to be supplemented by an equation to account for the axisymmetric motions, that is

$$\begin{aligned} \frac{\partial \overline{u_\phi}}{\partial t} + u_s \frac{\partial u_\phi}{\partial s} + \frac{\overline{u_s u_\phi}}{s} = & \nabla_{\perp}^2 \overline{u_\phi} - \frac{1}{s^2} \overline{u_\phi} \\ & + \frac{1}{Ek Pm} \left\langle (\mathbf{j} \times \mathbf{B}) \cdot \mathbf{e}_\phi \right\rangle \\ & - \left( \frac{s_o}{Ek} \right)^{1/2} \frac{1}{h^{3/2}} \overline{u_\phi}, \end{aligned} \quad (16)$$
where the last term of the right-hand side is the Ekman-pumping contribution to the axisymmetric motions.

## 2.2 Computation of the QG-Lorentz force

As seen in eqs (9) and (16), our method requires that we compute the following two quantities related to the 3-D Lorentz force averaged over axial direction

$$ \mathbf{F}_{\mathcal{L}, \omega_z} = \frac{1}{Ek Pm} \langle \nabla \times (\mathbf{j} \times \mathbf{B}) \cdot \mathbf{e}_z \rangle, \quad (17) $$

$$ \mathbf{F}_{\mathcal{L}, \overline{u}_{\phi}} = \frac{1}{Ek Pm} \langle (\mathbf{j} \times \mathbf{B}) \cdot \mathbf{e}_{\phi} \rangle . \quad (18) $$

We can expand eq. (17) as

$$ \mathbf{F}_{\mathcal{L}, \omega_z} = \frac{1}{Ek Pm} \left\langle \frac{1}{s} \left( \frac{\partial}{\partial s} [s(\mathbf{j} \times \mathbf{B}) \cdot \mathbf{e}_{\phi}] - \frac{\partial}{\partial \phi} [(\mathbf{j} \times \mathbf{B}) \cdot \mathbf{e}_s] \right) \right\rangle, \quad (19) $$

and using the identity $\mathbf{e}_s = \sin \theta \mathbf{e}_r + \cos \theta \mathbf{e}_\theta$ this becomes

$$ \mathbf{F}_{\mathcal{L}, \omega_z} = \frac{1}{Ek Pm} \left\langle \frac{1}{s} \left( \frac{\partial}{\partial s} [s(\mathbf{j} \times \mathbf{B}) \cdot \mathbf{e}_{\phi}] - \frac{\partial}{\partial \phi} [\sin \theta (\mathbf{j} \times \mathbf{B}) \cdot \mathbf{e}_r + \cos \theta (\mathbf{j} \times \mathbf{B}) \cdot \mathbf{e}_{\theta}] \right) \right\rangle . \quad (20) $$

Recalling that the $s$, $\phi$ and $z$ components are orthogonal, that the z-averaging operator is linear so $<u + v> = <u> + <v>$, and adopting the more compact notation $\mathbf{f} \cdot \mathbf{e}_x \equiv \mathbf{f}_x$ for any field $\mathbf{f}$ and coordinate $x$, the Lorentz force terms become

$$ \mathbf{F}_{\mathcal{L}} = \left\{ \begin{array}{ll} \frac{1}{Ek Pm} \left( \frac{1}{s} \left\langle \frac{\partial}{\partial s} [s(\mathbf{j} \times \mathbf{B})_{\phi}] \right\rangle \right) & \\ -\frac{1}{s} \frac{\partial}{\partial \phi} [(\sin \theta (\mathbf{j} \times \mathbf{B})_r + \cos \theta (\mathbf{j} \times \mathbf{B})_{\theta})] & . \quad (21) \\ \frac{1}{Ek Pm} \overline{\langle (\mathbf{j} \times \mathbf{B})_{\phi} \rangle} & \end{array} \right. $$

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),

$$ \mathbf{F}_{\mathcal{L}, \omega_z} = \frac{1}{Ek Pm} \left( \frac{1}{s} \left[ \frac{\partial}{\partial s} [s(\mathbf{j} \times \mathbf{B})_{\phi}] \right. \right. \\ \left. + \beta s \left( (\mathbf{j} \times \mathbf{B})_{\phi} - \frac{1}{2} [(\mathbf{j} \times \mathbf{B})_{\phi}] (\pm h) \right) \right] \\ \left. - \frac{1}{s} \frac{\partial}{\partial \phi} [(\sin \theta (\mathbf{j} \times \mathbf{B})_r + \cos \theta (\mathbf{j} \times \mathbf{B})_{\theta})] \right), \quad (22) $$

where the surface term $(\mathbf{j} \times \mathbf{B})_{\phi}(\pm h)$ cancels when $\mathbf{B}$ matches to a potential field at the outer boundary. This avoids the need to compute the s-derivative on the 3-D 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.

## 2.3 Numerics

The calculations presented in this study have been carried out using an extension of the open-source pseudospectral 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 2-D quantities are expanded in Fourier series up to the degree $N_m$ in the azimuthal

direction and in Chebyshev polynomials up to degree $N_s$ in the radial direction. The 3-D fields are expanded in Spherical Harmonics up to the degree and order $\ell_{max}$ in the angular $(\theta, \phi_{3D})$ directions and in Chebyshev polynomials with $N_r$ collocation gridpoints in the radial direction. The open-source SHTns<fnref n="1" /> library is used to handle the Spherical Harmonic Transforms (Schaeffer 2013). Parallelization of the hybrid QG-3D code relies on the Message Passing Interface (MPI) library.

The azimuthal (respectively radial) expansion involves adding zeros if $m_{3D} > m$ ($N_r > N_s$) and truncating fields at $m = m_{3D} = \ell$ ($N_s = N_r$) if $m_{3D} < m$ ($N_r < N_s$). The same 3-D and QG grids have been used when possible but because we explored cases with both $Pm < 1$ and $Pr < 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).

## 2.4 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 $\ell$, such that

$$ f_{H,u}(m) = \begin{cases} 1 & \text{for } m < m_H, \\ q_{H,u}^{m-m_H} & \text{for } m \ge m_H, \end{cases}. \quad (23) $$

on the velocity field, and

$$ f_{H,B}(\ell) = \begin{cases} 1 & \text{for } \ell < \ell_H, \\ q_{H,B}^{\ell-\ell_H} & \text{for } \ell \ge \ell_H, \end{cases}. \quad (24) $$

on the magnetic field, where $\ell_H$ and $m_H$ are the cut-off degrees and azimuthal orders below which the hyperdiffusion has no effect. $q_{H,u}$ and $q_{H,B}$ are the strength of the hyperdiffusive effect on $\mathbf{u}_{\perp}$ or $\mathbf{B}$, respectively and have been varied in the range $1.01 \le q_H \le 1.08$ (we do not apply any hyperdiffusion to the temperature field). The values for $m_H$ and $\ell_H$ for the runs using hyperdiffusion on either $\mathbf{u}_{\perp}$ or $\mathbf{B}$ are summarized in Table A1.

We have employed hyperdiffusion for two main reasons: (i) to mitigate closure issues when the 2-D and 3-D grids were different and (ii) to avoid the numerical problems in our most demanding runs—that is 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.

## 2.5 Diagnostics

We now introduce the following notations for our various integral and average operators. For any quantity $x$, the hat $\hat{x}$ corresponds to its time average, that is

$$ \hat{x} = \frac{1}{\tau} \int_{t_0}^{t_0+\tau} x \, dt, \quad (25) $$

where $t_0$ is chosen such that any transient has been passed and with $\tau$ the averaging time interval long enough to reach a statistical

<footnote n="1">¹https://bitbucket.org/nschaeff/shtns</footnote>
where the last term of the right-hand side is the Ekman-pumping contribution to the axisymmetric motions.

## 2.2 Computation of the QG-Lorentz force

As seen in eqs (9) and (16), our method requires that we compute the following two quantities related to the 3-D Lorentz force averaged over axial direction

$$ \mathbf{F}_{\mathcal{L}, \omega_z} = \frac{1}{Ek Pm} \langle \nabla \times (\mathbf{j} \times \mathbf{B}) \cdot \mathbf{e}_z \rangle, \quad (17) $$

$$ \mathbf{F}_{\mathcal{L}, \overline{u}_{\phi}} = \frac{1}{Ek Pm} \langle (\mathbf{j} \times \mathbf{B}) \cdot \mathbf{e}_{\phi} \rangle . \quad (18) $$

We can expand eq. (17) as

$$ \mathbf{F}_{\mathcal{L}, \omega_z} = \frac{1}{Ek Pm} \left\langle \frac{1}{s} \left( \frac{\partial}{\partial s} [s(\mathbf{j} \times \mathbf{B}) \cdot \mathbf{e}_{\phi}] - \frac{\partial}{\partial \phi} [(\mathbf{j} \times \mathbf{B}) \cdot \mathbf{e}_s] \right) \right\rangle, \quad (19) $$

and using the identity $\mathbf{e}_s = \sin \theta \mathbf{e}_r + \cos \theta \mathbf{e}_\theta$ this becomes

$$ \mathbf{F}_{\mathcal{L}, \omega_z} = \frac{1}{Ek Pm} \left\langle \frac{1}{s} \left( \frac{\partial}{\partial s} [s(\mathbf{j} \times \mathbf{B}) \cdot \mathbf{e}_{\phi}] - \frac{\partial}{\partial \phi} [\sin \theta (\mathbf{j} \times \mathbf{B}) \cdot \mathbf{e}_r + \cos \theta (\mathbf{j} \times \mathbf{B}) \cdot \mathbf{e}_{\theta}] \right) \right\rangle . \quad (20) $$

Recalling that the $s$, $\phi$ and $z$ components are orthogonal, that the z-averaging operator is linear so $<u + v> = <u> + <v>$, and adopting the more compact notation $\mathbf{f} \cdot \mathbf{e}_x \equiv \mathbf{f}_x$ for any field $\mathbf{f}$ and coordinate $x$, the Lorentz force terms become

$$ \mathbf{F}_{\mathcal{L}} = \left\{ \begin{array}{l} \frac{1}{Ek Pm} \left( \frac{1}{s} \left\langle \frac{\partial}{\partial s} [s(\mathbf{j} \times \mathbf{B})_{\phi}] \right\rangle \right) \\ -\frac{1}{s} \frac{\partial}{\partial \phi} [(\sin \theta (\mathbf{j} \times \mathbf{B})_r + \cos \theta (\mathbf{j} \times \mathbf{B})_{\theta})] \\ \frac{1}{Ek Pm} \left\langle (\mathbf{j} \times \mathbf{B})_{\phi} \right\rangle \end{array} \right. . \quad (21) $$

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),

$$ \mathbf{F}_{\mathcal{L}, \omega_z} = \frac{1}{Ek Pm} \left( \frac{1}{s} \left[ \frac{\partial}{\partial s} [s(\mathbf{j} \times \mathbf{B})_{\phi}] \right. \right. \\ \left. \left. + \beta s \left( (\mathbf{j} \times \mathbf{B})_{\phi} - \frac{1}{2} [(\mathbf{j} \times \mathbf{B})_{\phi}] (\pm h) \right) \right] \right. \\ \left. - \frac{1}{s} \frac{\partial}{\partial \phi} [(\sin \theta (\mathbf{j} \times \mathbf{B})_r + \cos \theta (\mathbf{j} \times \mathbf{B})_{\theta})] \right), \quad (22) $$

where the surface term $(\mathbf{j} \times \mathbf{B})_\phi(\pm h)$ cancels when $\mathbf{B}$ matches to a potential field at the outer boundary. This avoids the need to compute the $s$-derivative on the 3-D 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.

## 2.3 Numerics

The calculations presented in this study have been carried out using an extension of the open-source pseudospectral 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 2-D quantities are expanded in Fourier series up to the degree $N_m$ in the azimuthal

direction and in Chebyshev polynomials up to degree $N_s$ in the radial direction. The 3-D fields are expanded in Spherical Harmonics up to the degree and order $\ell_{max}$ in the angular $(\theta, \phi_{3D})$ directions and in Chebyshev polynomials with $N_r$ collocation gridpoints in the radial direction. The open-source SHTns<fnref n="1" /> library is used to handle the Spherical Harmonic Transforms (Schaeffer 2013). Parallelization of the hybrid QG-3D code relies on the Message Passing Interface (MPI) library.

The azimuthal (respectively radial) expansion involves adding zeros if $m_{3D} > m$ ($N_r > N_s$) and truncating fields at $m = m_{3D} = \ell$ ($N_s = N_r$) if $m_{3D} < m$ ($N_r < N_s$). The same 3-D and QG grids have been used when possible but because we explored cases with both $Pm < 1$ and $Pr < 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).

## 2.4 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 $\ell$, such that

$$ f_{H,u}(m) = \begin{cases} 1 & \text{for } m < m_H, \\ q_{H,u}^{m-m_H} & \text{for } m \ge m_H, \end{cases}. \quad (23) $$

on the velocity field, and

$$ f_{H,B}(\ell) = \begin{cases} 1 & \text{for } \ell < \ell_H, \\ q_{H,B}^{\ell-\ell_H} & \text{for } \ell \ge \ell_H, \end{cases}. \quad (24) $$

on the magnetic field, where $\ell_H$ and $m_H$ are the cut-off degrees and azimuthal orders below which the hyperdiffusion has no effect. $q_{H,u}$ and $q_{H,B}$ are the strength of the hyperdiffusive effect on $\mathbf{u}_\perp$ or $\mathbf{B}$, respectively and have been varied in the range $1.01 \le q_H \le 1.08$ (we do not apply any hyperdiffusion to the temperature field). The values for $m_H$ and $\ell_H$ for the runs using hyperdiffusion on either $\mathbf{u}_\perp$ or $\mathbf{B}$ are summarized in Table A1.

We have employed hyperdiffusion for two main reasons: (i) to mitigate closure issues when the 2-D and 3-D grids were different and (ii) to avoid the numerical problems in our most demanding runs—that is 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.

## 2.5 Diagnostics

We now introduce the following notations for our various integral and average operators. For any quantity $x$, the hat $\hat{x}$ corresponds to its time average, that is

$$ \hat{x} = \frac{1}{\tau} \int_{t_0}^{t_0+\tau} x \, dt, \quad (25) $$

where $t_0$ is chosen such that any transient has been passed and with $\tau$ the averaging time interval long enough to reach a statistical

<footnote n="1">¹https://bitbucket.org/nschaeff/shtns</footnote>
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

$$
\begin{align}
\{x\}_{3D} &\equiv \frac{1}{V_{3D}} \int_{V_{3D}} x \, dV, \quad \{x\}_{S} &\equiv \frac{1}{4\pi} \int_{0}^{\pi} \int_{0}^{2\pi} x \sin\theta \, d\theta \, d\phi_{3D}, \nonumber \\
\{x\}_{QG} &\equiv \frac{2}{V_{QG}} \int_{V_{QG}} x h(s) \, ds \, d\phi, \tag{26}
\end{align}
$$

with $V_{3D}$ the volume of the full spherical shell and $V_{QG}$ the volume outside of the tangent cylinder.

The dimensionless kinetic energy per unit volume $E_{\text{kin}}$, is defined by

$$E_{\text{kin}} = \frac{1}{2} \{\mathbf{u}_{\perp}^2\}_{\text{QG}}. \qquad (27)$$

We similarly define the dimensionless magnetic energy per unit volume $E_{\text{mag}}$, as

$$E_{\text{mag}} = \frac{1}{2} \{\mathbf{B}^2\}_{3D}. \qquad (28)$$

The magnetic to kinetic energy ratio is then

$$\mathcal{M} = \frac{\hat{E}_{\text{mag}}}{\hat{E}_{\text{kin}}}. \qquad (29)$$

From these expressions, we define a diagnostic for the fluid velocity which characterizes the average flow speed, based on the root-mean-square (r.m.s.) of the velocity, and which is denoted by the Reynolds number

$$Re = \sqrt{2E_{\text{kin}}}. \qquad (30)$$

Then the magnetic Reynolds number is simply $Rm = Pm Re$ and we can finally define the Elsasser number, a non-dimensional measure of the magnetic field strength, which reads

$$\Lambda = \frac{Ek}{Pm} Rm^2 \mathcal{M} = 2Ek Pm \hat{E}_{\text{mag}}. \qquad (31)$$

# 3 RESULTS

We present here results of experiments conducted at control parameters of $Ek = 3 \times 10^{-5}$, $Pr = 10^{-1}$, varying $Pm$, and focusing on a regime well above the onset of convection, that is $Ra \ge 5 Ra_c$. Time-integrating the non-linear equations in absence of a magnetic field (Barrois et al. 2022), we estimated the critical Rayleigh number for this configuration to be $Ra_c = 1.03 \times 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 \times 10^{-5}$, $Pr = 0.1$, $Pm = 0.9$, $Ra = 1.66 \times 10^7 \sim 16 Ra_c$— was started from a motionless fluid with a strong axial dipole with $\Lambda \sim 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.

## 3.1 Dynamo regime diagrams

A dynamo regime diagram as a function of $Pm$ and $Ra/Ra_c$ for all our runs conducted at $Ek = 3 \times 10^{-5}$, $Pr = 0.1$ is presented in Fig. 1. The crosses correspond to simulations which failed to pro-

duce a self-sustained dynamo while the stars represent the growing dynamos. Similar to the case for 3-D numerical dynamos we find that passing the onset for the dynamo action requires increasingly large $Pm$ 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 $Pr = 1$ while we consider $Pr = 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 $\mathcal{M} < 1$ and most of them have $\mathcal{M} < 10^{-1}$ (with the exception of some of the points, for example at $Pm = 2.0$, $Ra/Ra_c \sim 16.1$ that reached a moderate $\mathcal{M} = 0.15$ despite having the highest $\Lambda = 8.02$). They therefore fall into the weak-field dynamo regime, characterized by $\mathcal{M} \ll 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/Ra_c$ in Fig. 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 3-D simulations, with for example, a critical value of $Rm_c \sim 50$ reported by Petitdemange (2018). We found that $\Lambda > 0.05$ for all the cases when a self-sustained dynamo was identified—with the highest value $\Lambda \sim 8$ reached for the case at $Ek = 3 \times 10^{-5}$, $Pr = 0.1$, $Pm = 2.0$, $Ra = 1.66 \times 10^7 \sim 16.1 Ra_c$—and we generally observe an increase of $\Lambda$ with increasing $Rm$.

## 3.2 An example weak field dynamo

To illustrate the typical features of our data set, we look at an example dynamo with parameters $Ek = 3 \times 10^{-5}$, $Pr = 0.1$, $Pm = 0.9$, $Ra = 1.66 \times 10^7 \sim 16 Ra_c$ which was computed with a resolution of $(N_s, N_m)/(N_r, \ell_{\text{max}}) = (385, 768)/(149, 148)$ for approximately $\sim 5 \tau_v = 4.5 \tau_\lambda$. This case is rather close to the onset of the dynamo action and decreasing $Pm$ 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 \simeq 900$ and $\mathcal{M} = 0.04$.

Fig. 3 displays the time-averaged magnetic and kinetic energy spectra (a), a snapshot of the z-averaged axial vorticity $\omega_z$ (b), a meridional section of the longitudinally averaged azimuthal 3-D velocity $\overline{u}_{\phi 3D}$ (c), a snapshot of the radial magnetic field at the outer boundary of the dynamo region $B_r(r_o)$ (d), and a meridional section of the $\phi$-averaged azimuthal magnetic field $\overline{B}_\phi$ (e).

The power spectra shown in Fig. 3(a) confirms that the dynamo is multipolar (the magnetic field is dominated by the degree $\ell = 2$) and that the kinetic energy dominates over the magnetic energy at the largest length-scales by about 2 orders of magnitude (both for $\ell$ 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 $\ell = 2$) is typical of all our dynamos.

We can see in Figs 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. figs 3 and 4 in Barrois
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

$$
\begin{align}
\{x\}_{3D} &\equiv \frac{1}{V_{3D}} \int_{V_{3D}} x \, dV, \quad \{x\}_{S} &\equiv \frac{1}{4\pi} \int_{0}^{\pi} \int_{0}^{2\pi} x \sin\theta \, d\theta \, d\phi_{3D}, \nonumber \\
\{x\}_{QG} &\equiv \frac{2}{V_{QG}} \int_{V_{QG}} x h(s) \, ds \, d\phi, \tag{26}
\end{align}
$$

with $V_{3D}$ the volume of the full spherical shell and $V_{QG}$ the volume outside of the tangent cylinder.

The dimensionless kinetic energy per unit volume $E_{\text{kin}}$, is defined by

$$E_{\text{kin}} = \frac{1}{2} \{\mathbf{u}_{\perp}^2\}_{\text{QG}}. \qquad (27)$$

We similarly define the dimensionless magnetic energy per unit volume $E_{\text{mag}}$, as

$$E_{\text{mag}} = \frac{1}{2} \{\mathbf{B}^2\}_{3D}. \qquad (28)$$

The magnetic to kinetic energy ratio is then

$$M = \frac{\hat{E}_{\text{mag}}}{\hat{E}_{\text{kin}}}. \qquad (29)$$

From these expressions, we define a diagnostic for the fluid velocity which characterizes the average flow speed, based on the root-mean-square (r.m.s.) of the velocity, and which is denoted by the Reynolds number

$$Re = \sqrt{2E_{\text{kin}}}. \qquad (30)$$

Then the magnetic Reynolds number is simply $Rm = Pm Re$ and we can finally define the Elsasser number, a non-dimensional measure of the magnetic field strength, which reads

$$\Lambda = \frac{Ek}{Pm} Rm^2 M = 2Ek Pm \hat{E}_{\text{mag}}. \qquad (31)$$

# 3 RESULTS

We present here results of experiments conducted at control parameters of $Ek = 3 \times 10^{-5}$, $Pr = 10^{-1}$, varying $Pm$, and focusing on a regime well above the onset of convection, that is $Ra \ge 5 Ra_c$. Time-integrating the non-linear equations in absence of a magnetic field (Barrois et al. 2022), we estimated the critical Rayleigh number for this configuration to be $Ra_c = 1.03 \times 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 \times 10^{-5}$, $Pr = 0.1$, $Pm = 0.9$, $Ra = 1.66 \times 10^7 \sim 16 Ra_c$— was started from a motionless fluid with a strong axial dipole with $\Lambda \sim 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.

## 3.1 Dynamo regime diagrams

A dynamo regime diagram as a function of $Pm$ and $Ra/Ra_c$ for all our runs conducted at $Ek = 3 \times 10^{-5}$, $Pr = 0.1$ is presented in Fig. 1. The crosses correspond to simulations which failed to pro-

duce a self-sustained dynamo while the stars represent the growing dynamos. Similar to the case for 3-D numerical dynamos we find that passing the onset for the dynamo action requires increasingly large $Pm$ 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 $Pr = 1$ while we consider $Pr = 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, for example at $Pm = 2.0$, $Ra/Ra_c \sim 16.1$ that reached a moderate $M = 0.15$ despite having the highest $\Lambda = 8.02$). They therefore fall into the weak-field dynamo regime, characterized by $M \ll 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/Ra_c$ in Fig. 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 3-D simulations, with for example, a critical value of $Rm_c \sim 50$ reported by Petitdemange (2018). We found that $\Lambda > 0.05$ for all the cases when a self-sustained dynamo was identified—with the highest value $\Lambda \sim 8$ reached for the case at $Ek = 3 \times 10^{-5}$, $Pr = 0.1$, $Pm = 2.0$, $Ra = 1.66 \times 10^7 \sim 16.1 Ra_c$—and we generally observe an increase of $\Lambda$ with increasing $Rm$.

## 3.2 An example weak field dynamo

To illustrate the typical features of our data set, we look at an example dynamo with parameters $Ek = 3 \times 10^{-5}$, $Pr = 0.1$, $Pm = 0.9$, $Ra = 1.66 \times 10^7 \sim 16 Ra_c$ which was computed with a resolution of $(N_s, N_m)/(N_r, \ell_{\text{max}}) = (385, 768)/(149, 148)$ for approximately $\sim 5 \tau_v = 4.5 \tau_\lambda$. This case is rather close to the onset of the dynamo action and decreasing $Pm$ 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 \simeq 900$ and $M = 0.04$.

Fig. 3 displays the time-averaged magnetic and kinetic energy spectra (a), a snapshot of the z-averaged axial vorticity $\omega_z$ (b), a meridional section of the longitudinally averaged azimuthal 3-D velocity $\overline{u}_{\phi 3D}$ (c), a snapshot of the radial magnetic field at the outer boundary of the dynamo region $B_r(r_o)$ (d), and a meridional section of the $\phi$-averaged azimuthal magnetic field $\overline{B}_\phi$ (e).

The power spectra shown in Fig. 3(a) confirms that the dynamo is multipolar (the magnetic field is dominated by the degree $\ell = 2$) and that the kinetic energy dominates over the magnetic energy at the largest length-scales by about 2 orders of magnitude (both for $\ell$ 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 $\ell = 2$) is typical of all our dynamos.

We can see in Figs 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. figs 3 and 4 in Barrois
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Magnetic Prandtl number, *Pm*, as a function of the supercriticality, *Ra/Ra*<sub>c</sub>. Dynamo regime diagram computed for a series of cases at *Ek* = 3 × 10<sup>-5</sup> and *Pr* = 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 *Pm*. The dashed line marks the tentative limit between failed and growing dynamos.</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>**Figure 2.** Magnetic Reynolds *Rm* as a function of the supercriticality *Ra/Ra*<sub>c</sub>. The symbols and the line have the same signification as in Fig. 1.</figcaption></figure>

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_{\phi 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 Figs 3(d) and (e) that $B_r(r=r_o)$ (d) and $\bar{B}_\phi$ (e) have fairly low amplitude, and are mostly non-dipolar despite the clear equatorial anti-symmetry. $B_r$ 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 \times 10^{-5}$ it is expected to have some impact close to the outer boundary. However, 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 $u_{3D}$—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 $\bar{u}_{\phi 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 $\beta$-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 equatorially antisymmetric, 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 antisymmetry) 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.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Magnetic Prandtl number, *Pm*, as a function of the supercriticality, *Ra/Ra*<sub>c</sub>. Dynamo regime diagram computed for a series of cases at *Ek* = 3 × 10<sup>-5</sup> and *Pr* = 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 *Pm*. The dashed line marks the tentative limit between failed and growing dynamos.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Magnetic Reynolds *Rm* as a function of the supercriticality *Ra/Ra*<sub>c</sub>. The symbols and the line have the same signification as in Fig. 1.</figcaption></figure>

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_{\phi 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 Figs 3(d) and (e) that $B_r(r=r_o)$ (d) and $\bar{B}_\phi$ (e) have fairly low amplitude, and are mostly non-dipolar despite the clear equatorial anti-symmetry. $B_r$ 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 \times 10^{-5}$ it is expected to have some impact close to the outer boundary. However, 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 $u_{3D}$—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 $\bar{u}_{\phi 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 $\beta$-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 equatorially antisymmetric, 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 antisymmetry) 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.
<figure><img src="image_3.png" /><figcaption>**Figure 3.** (a) Time-averaged magnetic and kinetic energy spectra. (b) Snapshot of the z-averaged vorticity $\omega_z$. (c) Meridional section of the $\phi$-averaged azimuthal 3-D velocity $\bar{u}_{\phi 3D}$. (d) Snapshot of the radial magnetic field at the outer boundary $B_r(r_o)$. (e) Meridional section of the $\phi$-averaged azimuthal magnetic field $\bar{B}_{\phi}$. For a dynamo with control parameters $Ek = 3 \times 10^{-5}$, $Pr = 0.1$, $Pm = 0.9$, $Ra = 1.66 \times 10^7 \sim 16 Ra_c$.</figcaption></figure>

### 3.3 Dynamo mechanism

To illustrate the main mechanism underlying our dynamos, Fig. 4 shows, from left to right, respectively, meridional sections of the $\phi$-averaged 3-D helicity $\mathcal{H} = \overline{u}_{3D} \cdot (\nabla \times \mathbf{u}_{3D})$, the azimuthal component of the mean electromotive force ($\mathbf{u}'_{3D} \times \mathbf{B}'$) $\cdot \mathbf{e}_\phi$ (hereafter EMF) and an estimate of the (rescaled) mean $\alpha$-effect $-\frac{\tau_c}{3} B_\phi \cdot (\mathbf{u}'_{3D} \cdot \nabla \times \mathbf{u}'_{3D})$ for the same dynamo presented in the previous section—where the prime $x'$ denotes fluctuations about the azimuthal average of any quantity $x$, that is $x' \equiv x - \bar{x}$ and, where $\tau_c$ is a typical convective turnover time (Brandenburg & Subramanian 2005) which can be approximated by $\tau_c \sim d/\{\mathbf{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 $\alpha$-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 $\bar{B}_\phi$ and strong zonal wind $\overline{u}_{\phi 3D}$ have already been observed in the Fig. 3. Most of the 3-D helicity is contained in its z-component ($\overline{u}_z \cdot (\nabla \times \mathbf{u}_{3D}) \cdot \mathbf{e}_z$) with only weak contributions from the s- and $\phi$-components (not shown). There is clearly a strong correlation between the mean EMF and the estimated $\alpha$-effect despite the discrepancy in amplitudes (Figs 4b and c) with both quantities located where the maximum helicity and toroidal field are found,

indicating that (i) our model produces an $\alpha$-effect sufficient to support the dynamo action and (ii) kinetic helicity provides a reasonable mean-field approximation of the actual EMF.

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 $\Omega$-effect $s\bar{\mathbf{B}} \cdot \nabla \left(\frac{\overline{u}_{\phi 3D}}{s}\right)$ (e.g. Roberts & King 2013) with superimposed toroidal field lines (top row) and of the azimuthal component of the mean EMF ($\mathbf{u}'_{3D} \times \mathbf{B}'$) $\cdot \mathbf{e}_\phi$ with superimposed poloidal field lines (bottom row). A strong correlation is found between both the $\Omega$-effect and the toroidal field lines and between the EMF and the poloidal field lines, and both effects are concentrated in the same region localized in the upper part of the shell. The $\Omega$-effect is stronger than the $\alpha$-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—that is the correlations between the azimuthal motions $\mathbf{u}'_{\phi 3D}$ and the cylindrically radial velocity $\mathbf{u}'_s$—are expected to produce an important $\Omega$-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
<figure><img src="image_4.png" /><figcaption>**Figure 3.** (a) Time-averaged magnetic and kinetic energy spectra. (b) Snapshot of the z-averaged vorticity $\omega_z$. (c) Meridional section of the $\phi$-averaged azimuthal 3-D velocity $\bar{u}_{\phi 3D}$. (d) Snapshot of the radial magnetic field at the outer boundary $B_r(r_o)$. (e) Meridional section of the $\phi$-averaged azimuthal magnetic field $\bar{B}_{\phi}$. For a dynamo with control parameters $Ek = 3 \times 10^{-5}$, $Pr = 0.1$, $Pm = 0.9$, $Ra = 1.66 \times 10^7 \sim 16 Ra_c$.</figcaption></figure>

### 3.3 Dynamo mechanism

To illustrate the main mechanism underlying our dynamos, Fig. 4 shows, from left to right, respectively, meridional sections of the $\phi$-averaged 3-D helicity $\mathcal{H} \equiv \bar{u}_{3D} \cdot (\nabla \times \bar{u}_{3D})$, the azimuthal component of the mean electromotive force ($\boldsymbol{u}'_{3D} \times \boldsymbol{B}'$) $\cdot \boldsymbol{e}_\phi$ (hereafter EMF) and an estimate of the (rescaled) mean $\alpha$-effect $-\frac{\tau_c}{3} B_\phi \cdot (\boldsymbol{u}'_{3D} \cdot \nabla \times \boldsymbol{u}'_{3D})$ for the same dynamo presented in the previous section—where the prime $x'$ denotes fluctuations about the azimuthal average of any quantity $x$, that is $x' \equiv x - \bar{x}$ and, where $\tau_c$ is a typical convective turnover time (Brandenburg & Subramanian 2005) which can be approximated by $\tau_c \sim d/\{\boldsymbol{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 $\alpha$-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 $\bar{B}_\phi$ and strong zonal wind $\bar{u}_{\phi 3D}$ have already been observed in the Fig. 3. Most of the 3-D helicity is contained in its z-component ($\boldsymbol{u}_z \cdot (\nabla \times \boldsymbol{u}_{3D}) \cdot \boldsymbol{e}_z$) with only weak contributions from the s- and $\phi$-components (not shown). There is clearly a strong correlation between the mean EMF and the estimated $\alpha$-effect despite the discrepancy in amplitudes (Figs 4b and c) with both quantities located where the maximum helicity and toroidal field are found,

indicating that (i) our model produces an $\alpha$-effect sufficient to support the dynamo action and (ii) kinetic helicity provides a reasonable mean-field approximation of the actual EMF.

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 $\Omega$-effect $s\bar{\mathbf{B}} \cdot \nabla \left(\frac{\boldsymbol{u}_{\phi 3D}}{s}\right)$ (e.g. Roberts & King 2013) with superimposed toroidal field lines (top row) and of the azimuthal component of the mean EMF ($\boldsymbol{u}'_{3D} \times \boldsymbol{B}'$) $\cdot \boldsymbol{e}_\phi$ with superimposed poloidal field lines (bottom row). A strong correlation is found between both the $\Omega$-effect and the toroidal field lines and between the EMF and the poloidal field lines, and both effects are concentrated in the same region localized in the upper part of the shell. The $\Omega$-effect is stronger than the $\alpha$-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—that is the correlations between the azimuthal motions $\boldsymbol{u}'_{\phi 3D}$ and the cylindrically radial velocity $\boldsymbol{u}'_s$—are expected to produce an important $\Omega$-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
<figure><img src="image_4.png" /><figcaption>**Figure 4.** (a) Meridional section of the longitudinally averaged 3-D helicity $\overline{u'_{3D} \cdot \nabla \times u_{3D}}$. (b) Azimuthal component of the mean electromotive force $\overline{(u'_{3D} \times B') \cdot e_{\phi}}$. (c) Estimated mean $\alpha$-effect $-\frac{d}{3\{u'_{3D}\}_S} B_{\phi} \cdot (u'_{3D} \cdot \nabla \times u'_{3D})$. For the same case as in Fig. 3.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 5.** Temporal evolution of $\phi$-averaged snapshots of the $\Omega$-effect $s\bar{B} \cdot \nabla \left( \frac{\mu_{\phi 3D}}{s} \right)$ with superimposed toroidal field lines (top row) and of the azimuthal component of the mean EMF $\overline{(u'_{3D} \times B') \cdot e_{\phi}}$ 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.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** (a) Meridional section of the longitudinally averaged 3-D helicity $\overline{u'_{3D} \cdot \nabla \times u_{3D}}$. (b) Azimuthal component of the mean electromotive force $\overline{(u'_{3D} \times B') \cdot e_{\phi}}$. (c) Estimated mean $\alpha$-effect $-\frac{d}{3\{u'_{3D}\}_S} B_{\phi} \cdot (u'_{3D} \cdot \nabla \times u'_{3D})$. For the same case as in Fig. 3.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Temporal evolution of $\phi$-averaged snapshots of the $\Omega$-effect $s\bar{B} \cdot \nabla \left( \frac{\mu_{\phi 3D}}{s} \right)$ with superimposed toroidal field lines (top row) and of the azimuthal component of the mean EMF $\overline{(u'_{3D} \times B') \cdot e_{\phi}}$ 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.</figcaption></figure>
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 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 α²Ω type, following the mean field nomenclature of, for example Steenbeck et al. (1966), Steenbeck & Krause (1969) or Malkus & Proctor (1975).

## 3.4 Dynamo waves

Fig. 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 3-D 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 Fig. 6 we show the temporal evolution of $\overline{B_\phi}(\theta)$ just beneath the outer boundary at $r \sim 0.92 r_o$—a so-called ‘butterfly-diagram’— for a case at $Ek = 3 \times 10^{-5}$, $Pr = 0.1$, $Pm = 0.9$, $Ra = 1.66 \times 10^7 \sim 16 Ra_c$ with a $Rm$ of 909 (a)—same dynamo as presented in Section 3.2—along with a more strongly driven case at $Ek = 3 \times 10^{-5}$, $Pr = 0.1$, $Pm = 0.9$, $Ra = 3.0 \times 10^7 \sim 30 Ra_c$ 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, for example 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. 3c). 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 \times 10^{-3} \tau_\lambda$ and 1 cycle per $2.7 \times 10^{-3} \tau_\lambda$, respectively. This is consistent with the expectation that the reversal frequency should increase with increasing forcing or increasing field strength; we have increased $Ra/Ra_c$ leading to an increase in both $Rm$ and the magnetic field strength—between cases (a) and (b).

In order to better characterize these oscillations, we test the approximate dispersion relation for (Parker) dynamo waves derived in Schrinner et al. (2011), which reads

$$ \omega \sim \left( \frac{H'_{\text{Int}}}{2r_o} \frac{Re_z}{Re_c} \right)^{1/2}, \qquad (32) $$

$$ \text{with } Re_z = \sqrt{\langle u_{\phi 3D}^2 \rangle_{3D}}, \quad Re_c = \sqrt{\langle (u'_{3D})^2 \rangle_{3D}} \text{ and } H'_{\text{Int}} = -\frac{1}{3} \{ |u'_{3D}| \cdot \nabla \times u'_{3D} | \}_{3D}. $$

Fig. 7 presents the resulting comparison between the oscillation 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 3-D 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.

## 3.5 Comparison with 3-D dynamos

In order to compare our results with more conventional 3-D simulations, we have computed a series of cases varying $Pm$ from $Pm = 0.5$ to $Pm = 0.05$ at $Ek = 3 \times 10^{-5}$, $Pr = 10^{-1}$, $Ra = 2 \times 10^7 \sim 20 Ra_c$ using a full 3-D method (MagIC, Wicht 2002). Here we compare these results with a hybrid QG-3D case at $Ek = 3 \times 10^{-5}$, $Pm = 0.5$, $Pr = 10^{-1}$, $Ra = 2 \times 10^7 \sim 20 Ra_c$.

Fig. 8 displays results for one of these 3-D simulations, computed at $Pm = 0.1$ (left-hand column) and for the hybrid QG-3D simulation computed at $Pm = 0.5$ (right-hand column). Snapshots of the z-averaged vorticity (panels a and e), meridional sections of the φ-averaged azimuthal 3-D velocity (panels b and f), meridional sections of the φ-averaged kinetic helicity $\dot{u}_{3D} \cdot \nabla \times \dot{u}_{3D}$ (panels c and g), and time series of the axisymmetric azimuthal magnetic field at $r \sim 0.92 r_o$ (panels d and h) are presented for both the 3-D and hybrid cases.

We stress that the comparison involves different $Pm$ for both setups. This is necessary because the onset for dynamo action happens at a lower $Pm$ with the 3-D 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 3-D and hybrid QG-3D results at the same control parameters, qualitative differences are found—the solution is already bi-stable in the 3-D 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 $Pm$, qualitatively similar solutions are then obtained when we are close to the onset of dynamo action for both configurations. Figs 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 3-D and hybrid QG-3D cases belong to the weak-field multipolar-branch 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.

Figs 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 (Figs 8c–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 3-D 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 3-D simulations (see, for example Davidson & Ranjan 2018).

Finally, we can see in Figs 8(d)–(h) that both the 3-D and hybrid QG-3D simulations display similar Parker-wave oscillations
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 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 α²Ω type, following the mean field nomenclature of, for example Steenbeck et al. (1966), Steenbeck & Krause (1969) or Malkus & Proctor (1975).

## 3.4 Dynamo waves

Fig. 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 3-D 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 Fig. 6 we show the temporal evolution of $\overline{B_\phi}(\theta)$ just beneath the outer boundary at $r \sim 0.92 r_o$—a so-called ‘butterfly-diagram’— for a case at $Ek = 3 \times 10^{-5}$, $Pr = 0.1$, $Pm = 0.9$, $Ra = 1.66 \times 10^7 \sim 16 Ra_c$ with a $Rm$ of 909 (a)—same dynamo as presented in Section 3.2—along with a more strongly driven case at $Ek = 3 \times 10^{-5}$, $Pr = 0.1$, $Pm = 0.9$, $Ra = 3.0 \times 10^7 \sim 30 Ra_c$ 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, for example 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. 3c). 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 \times 10^{-3} \tau_\lambda$ and 1 cycle per $2.7 \times 10^{-3} \tau_\lambda$, respectively. This is consistent with the expectation that the reversal frequency should increase with increasing forcing or increasing field strength; we have increased $Ra/Ra_c$ leading to an increase in both $Rm$ and the magnetic field strength—between cases (a) and (b).

In order to better characterize these oscillations, we test the approximate dispersion relation for (Parker) dynamo waves derived in Schrinner et al. (2011), which reads

$$ \omega \sim \left( \frac{H'_{\text{Int}}}{2r_o} \frac{Re_z}{Re_c} \right)^{1/2}, \qquad (32) $$

$$ \text{with } Re_z = \sqrt{\langle u_{\phi 3D}^2 \rangle_{3D}}, \quad Re_c = \sqrt{\langle (u'_{3D})^2 \rangle_{3D}} \text{ and } H'_{\text{Int}} = -\frac{1}{3} \{ |u'_{3D}| \cdot \nabla \times u'_{3D} | \}_{3D}. $$

Fig. 7 presents the resulting comparison between the oscillation 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 3-D 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.

## 3.5 Comparison with 3-D dynamos

In order to compare our results with more conventional 3-D simulations, we have computed a series of cases varying $Pm$ from $Pm = 0.5$ to $Pm = 0.05$ at $Ek = 3 \times 10^{-5}$, $Pr = 10^{-1}$, $Ra = 2 \times 10^7 \sim 20 Ra_c$ using a full 3-D method (MagIC, Wicht 2002). Here we compare these results with a hybrid QG-3D case at $Ek = 3 \times 10^{-5}$, $Pm = 0.5$, $Pr = 10^{-1}$, $Ra = 2 \times 10^7 \sim 20 Ra_c$.

Fig. 8 displays results for one of these 3-D simulations, computed at $Pm = 0.1$ (left-hand column) and for the hybrid QG-3D simulation computed at $Pm = 0.5$ (right-hand column). Snapshots of the z-averaged vorticity (panels a and e), meridional sections of the φ-averaged azimuthal 3-D velocity (panels b and f), meridional sections of the φ-averaged kinetic helicity $\dot{u}_{3D} \cdot \nabla \times \dot{u}_{3D}$ (panels c and g), and time series of the axisymmetric azimuthal magnetic field at $r \sim 0.92 r_o$ (panels d and h) are presented for both the 3-D and hybrid cases.

We stress that the comparison involves different $Pm$ for both setups. This is necessary because the onset for dynamo action happens at a lower $Pm$ with the 3-D 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 3-D and hybrid QG-3D results at the same control parameters, qualitative differences are found—the solution is already bi-stable in the 3-D 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 $Pm$, qualitatively similar solutions are then obtained when we are close to the onset of dynamo action for both configurations. Figs 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 3-D and hybrid QG-3D cases belong to the weak-field multipolar-branch 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.

Figs 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 (Figs 8c–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 3-D 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 3-D simulations (see, for example Davidson & Ranjan 2018).

Finally, we can see in Figs 8(d)–(h) that both the 3-D and hybrid QG-3D simulations display similar Parker-wave oscillations
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Axisymmetric azimuthal magnetic field $\overline{B}_{\phi}(\theta)$ at $r \sim 0.92 r_o$ as a function of time—expressed in viscous units $\tau_v$—(so called butterfly diagram) for (a) dynamo with control parameters $Ek = 3 \times 10^{-5}$, $Pr = 0.1$, $Pm = 0.9$, $Ra = 1.66 \times 10^7$ and (b) a second case, more strongly forced, with control parameters $Ek = 3 \times 10^{-5}$, $Pr = 0.1$, $Pm = 0.9$, $Ra = 3.0 \times 10^7$. In both panels, dashed lines mark the location of the tangent cylinder.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 7.** Measured frequencies of oscillatory dynamos as a function of the frequency predicted by the Parker dynamo wave dispersion relation (eq. 32), for all of our dynamos. Different colors correspond to different magnetic Prandtl numbers and the dashed line corresponds to the expected relation (eq. 32).</figcaption></figure>

although the magnetic field is equatorially symmetric rather than antisymmetric and contains more small scale features in the 3-D case. There is also much more magnetic energy in the 3-D case with $\max \overline{B}_{\phi}^{3D} > 10 \times \max \overline{B}_{\phi}^{Hyb}$. A rough estimation of the Parker cycles gives an oscillation period of 1 cycle per $5 \times 10^{-3} \tau_v$ for the 3-D case and 1 cycle per $8.3 \times 10^{-3} \tau_v$ for the hybrid QG-3D case, consistent with the fact that $Re_c$ and $Re_z$ are similar in the two simulations but $\mathcal{H}_{3D} \gg \mathcal{H}_{Hyb}$ (see eq. 32).

At $Pm = 0.5$, $Ra \sim 20 Ra_c$, we have computed the total rms kinetic helicity—as a proxy of the amplitude of the $\alpha$-effect— for our weak field hybrid QG-3D case, as well as for a multipolar weak field 3-D case, and a dipolar strong field 3-D case (at these parameters the 3-D run is bi-stable). Both 3-D 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
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Axisymmetric azimuthal magnetic field $\overline{B}_{\phi}(\theta)$ at $r \sim 0.92 r_o$ as a function of time—expressed in viscous units $\tau_v$—(so called butterfly diagram) for (a) dynamo with control parameters $Ek = 3 \times 10^{-5}$, $Pr = 0.1$, $Pm = 0.9$, $Ra = 1.66 \times 10^7$ and (b) a second case, more strongly forced, with control parameters $Ek = 3 \times 10^{-5}$, $Pr = 0.1$, $Pm = 0.9$, $Ra = 3.0 \times 10^7$. In both panels, dashed lines mark the location of the tangent cylinder.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Measured frequencies of oscillatory dynamos as a function of the frequency predicted by the Parker dynamo wave dispersion relation (eq. 32), for all of our dynamos. Different colors correspond to different magnetic Prandtl numbers and the dashed line corresponds to the expected relation (eq. 32).</figcaption></figure>

although the magnetic field is equatorially symmetric rather than antisymmetric and contains more small scale features in the 3-D case. There is also much more magnetic energy in the 3-D case with $\max \overline{B}_{\phi}^{3D} > 10 \times \max \overline{B}_{\phi}^{Hyb}$. A rough estimation of the Parker cycles gives an oscillation period of 1 cycle per $5 \times 10^{-3} \tau_v$ for the 3-D case and 1 cycle per $8.3 \times 10^{-3} \tau_v$ for the hybrid QG-3D case, consistent with the fact that $Re_c$ and $Re_z$ are similar in the two simulations but $\mathcal{H}_{3D} \gg \mathcal{H}_{Hyb}$ (see eq. 32).

At $Pm = 0.5$, $Ra \sim 20 Ra_c$, we have computed the total rms kinetic helicity—as a proxy of the amplitude of the $\alpha$-effect— for our weak field hybrid QG-3D case, as well as for a multipolar weak field 3-D case, and a dipolar strong field 3-D case (at these parameters the 3-D run is bi-stable). Both 3-D 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
<figure><img src="image_8.png" /><figcaption>**Figure 8.** Comparison of results from a 3-D dynamo simulation computed at $Ek = 3 \times 10^{-5}$, $Pr = 0.1$, $Pm = 0.1$, $Ra = 20 \times Ra_c$ (panels a to d) and a hybrid QG-3D case at $Ek = 3 \times 10^{-5}$, $Pr = 0.1$, $Pm = 0.5$, $Ra = 20 \times Ra_c$ (panels e to h). Panels (a) and (e) show snapshot of the z-averaged vorticity, (b) and (f) meridional section of the φ-averaged azimuthal 3-D velocity, (c) and (g) meridional section of the φ-averaged 3-D helicity, (d) and (h) time series of the axisymmetric azimuthal magnetic field at $r \sim 0.92 r_o$.</figcaption></figure>

the 3-D case, even if the spatial segregation of helicity remains as expected.

Overall, we conclude that although the comparison of 3-D 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 3-D methods do produce qualitatively similar results when considered at the same distance from the onset of dynamo action, and provided the 3-D case is started from a weak seed. The major difference is that the hybrid QG-3D model involves much lower levels of kinetic helicity.

# 4 DISCUSSION AND CONCLUSIONS

The results above demonstrate that it is possible to obtain dynamo action in a hybrid QG-3D magnetohydrodynamic model of rapidly rotating 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 \times 10^{-5}$, $Pr = 0.1 - 1$, $Ra \sim (5 - 100) Ra_c$ and $Pm = 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 \times 10^{-5}$, $Pr = 0.1$ and a range of $Pm$ and $Ra/Ra_c$, we found QG dynamos characterized by a low magnetic to kinetic energy ratio $M$, a multipolar magnetic
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Comparison of results from a 3-D dynamo simulation computed at $Ek = 3 \times 10^{-5}$, $Pr = 0.1$, $Pm = 0.1$, $Ra = 20 \times Ra_c$ (panels a to d) and a hybrid QG-3D case at $Ek = 3 \times 10^{-5}$, $Pr = 0.1$, $Pm = 0.5$, $Ra = 20 \times Ra_c$ (panels e to h). Panels (a) and (e) show snapshot of the z-averaged vorticity, (b) and (f) meridional section of the φ-averaged azimuthal 3-D velocity, (c) and (g) meridional section of the φ-averaged 3-D helicity, (d) and (h) time series of the axisymmetric azimuthal magnetic field at $r \sim 0.92 r_o$.</figcaption></figure>

the 3-D case, even if the spatial segregation of helicity remains as expected.

Overall, we conclude that although the comparison of 3-D 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 3-D methods do produce qualitatively similar results when considered at the same distance from the onset of dynamo action, and provided the 3-D case is started from a weak seed. The major difference is that the hybrid QG-3D model involves much lower levels of kinetic helicity.

# 4 DISCUSSION AND CONCLUSIONS

The results above demonstrate that it is possible to obtain dynamo action in a hybrid QG-3D magnetohydrodynamic model of rapidly rotating 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 \times 10^{-5}$, $Pr = 0.1 - 1$, $Ra \sim (5 - 100) Ra_c$ and $Pm = 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 \times 10^{-5}$, $Pr = 0.1$ and a range of $Pm$ and $Ra/Ra_c$, we found QG dynamos characterized 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 3-D 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 3-D models in a variety of set-ups (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 parameter space at $Pr = 1$, for the cases tested we found similar results with weak and multipolar magnetic field solutions, suggesting that having $Pr = 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 3-D dynamos the poloidal field in our dynamos always remains weak compared with the toroidal field. This remains true, even at other values of $Pr$ (e.g. runs at $Pr = 1$, not shown), which suggests that our model may lack some important ingredient in the field generation cycle that operates in dipole-dominated 3-D 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 3-D 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 $u_z$ and the temperature $T_{3D}$ 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 parametrize 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 coworkers 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 parametrized 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 3-D 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 3-D 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 3-D 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 3-D 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 parametrizing 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 3-D 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 3-D magnetic field perturbations and has the potential to be significantly faster than full 3-D simulations for studying such waves.

## ACKNOWLEDGMENTS

We thank two anonymous reviewers and the editor Richard Holme for their constructive comments which have helped to improve the manuscript. Nathanaël Schaeffer and Céline Guervilly are thanked for helpful discussions and for sharing their code for the z-averaging functions and for the z-integration of the thermal wind. QG and hybrid QG-3D numerical computations using the pizza code were performed at DTU Space, using the Humboldt and Larmor CPU clusters. OB and CCF were supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 772561).
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 3-D 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 3-D models in a variety of set-ups (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 parameter space at $Pr = 1$, for the cases tested we found similar results with weak and multipolar magnetic field solutions, suggesting that having $Pr = 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 3-D dynamos the poloidal field in our dynamos always remains weak compared with the toroidal field. This remains true, even at other values of $Pr$ (e.g. runs at $Pr = 1$, not shown), which suggests that our model may lack some important ingredient in the field generation cycle that operates in dipole-dominated 3-D 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 3-D 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 $u_z$ and the temperature $T_{3D}$ 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 parametrize 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 coworkers 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 parametrized 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 3-D 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 3-D 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 3-D 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 3-D 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 parametrizing 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 3-D 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 3-D magnetic field perturbations and has the potential to be significantly faster than full 3-D simulations for studying such waves.

## ACKNOWLEDGMENTS

We thank two anonymous reviewers and the editor Richard Holme for their constructive comments which have helped to improve the manuscript. Nathanaël Schaeffer and Céline Guervilly are thanked for helpful discussions and for sharing their code for the z-averaging functions and for the z-integration of the thermal wind. QG and hybrid QG-3D numerical computations using the pizza code were performed at DTU Space, using the Humboldt and Larmor CPU clusters. OB and CCF were supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 772561).
DATA AVAILABILITY

All our data are available upon reasonable request to the corre-
sponding author and some key parameters from our whole data-
set are already included in this paper. The codes used (pizza) is
freely available (at http://www.github.com/magic-sph/pizza/tree/hy
brid_QG-3D) under the GNU GPL v3 license.

REFERENCES

Aubert, J., 2018. Geomagnetic acceleration and rapid hydromagnetic wave dynamics in advanced numerical simulations of the geodynamo, *Geophys. J. Int.*, **214**(1), 531–547.

Aubert, J., 2019. Approaching earth's core conditions in high-resolution geodynamo simulations, *Geophys. J. Int.*, **219**(Suppl 1), S137-S151.

Aubert, J., 2023. State and evolution of the geodynamo from numerical models reaching the physical conditions of Earth's core, *Geophys. J. Int.*, **235**(1), 468–487.

Aubert, J., Gillet, N. & Cardin, P., 2003. Quasigeostrophic models of convection in rotating spherical shells, *Geochem. Geophys. Geosyst.*, **4**(7), 1–19.

Aubert, J., Finlay, C. C. & Fournier, A., 2013. Bottom-up control of geomagnetic secular variation by the Earth's inner core, *Nature*, **502**(7470), 219–223.

Aubert, J., Gastine, T. & Fournier, A., 2017. Spherical convective dynamos in the rapidly rotating asymptotic regime, *J. Fluid Mech.*, **813**, 558–593.

Aubert, J., Livermore, P.W., Finlay, C.C., Fournier, A. & Gillet, N., 2022. A taxonomy of simulated geomagnetic jerks, *Geophys. J. Int.*, **231**(1), 650–672.

Aurnou, J. & King, E., 2017. The cross-over to magnetostrophic convection in planetary dynamo systems, *Proc. R. Soc. A*, **473**(2199), doi:10.1098/rspa.2016.0731.

Barrois, O., Gastine, T. & Finlay, C., 2022. Comparison of quasi-geostrophic, hybrid and 3-D models of planetary core convection, *Geophys. J. Int.*, **231**(1), 129–158.

Brandenburg, A. & Subramanian, K., 2005. Astrophysical magnetic fields and nonlinear dynamo theory, *Phys. Rep.*, **417**(1-4), 1-209.

Brown, B.P., Browning, M.K., Brun, A.S., Miesch, M.S. & Toomre, J., 2010. Persistent magnetic wreaths in a rapidly rotating sun, *Astrophys. J.*, **711**(1), 424–438.

Busse, F.H., 1970. Thermal instabilities in rapidly rotating systems, *J. Fluid Mech.*, **44**(3), 441–460.

Busse, F.H., 1975. A model of the geodynamo, *Geophys. J. Int.*, **42**(2), 437–459.

Busse, F.H. & Simitev, R.D., 2006. Parameter dependences of convection-
driven dynamos in rotating spherical fluid shells, *Geophys. astrophys.
Fluid. Dyn.*, **100**(4–5), 341–361.

Calkins, M.A., 2018. Quasi-geostrophic dynamo theory, *Phys. Earth planet. Inter.*, **276**, 182–189.

Canet, E., Finlay, C.C. & Fournier, A., 2014. Hydromagnetic quasi-
geostrophic modes in rapidly rotating planetary cores, *Phys. Earth planet.
Inter.*, **229**, 1–15.

Cardin, P. & Olson, P., 1994. Chaotic thermal convection in a rapidly rotating spherical shell: consequences for flow in the outer core, *Phys. Earth planet. Inter.*, **82**(3-4), 235-259.

Chan, K.H., Zhang, K., Zou, J. & Schubert, G., 2001. A non-linear, 3-D spherical α2 dynamo using a finite element method, *Phys. Earth planet. Inter.*, **128**(1-4), 35–50.

Christensen, U.R. & Aubert, J., 2006. Scaling properties of convection-
driven dynamos in rotating spherical shells and application to planetary
magnetic fields, *Geophys. J. Int.*, **166**(1), 97–114.

Davidson, P., 2013. Scaling laws for planetary dynamos, *Geophys. J. Int.*, **195**(1), 67–74.

Davidson, P. & Ranjan, A., 2015. Planetary dynamos driven by helical waves
— II, *Geophys. J. Int.*, **202**(3), 1646–1662.

Davidson, P.A. & Ranjan, A., 2018. Are planetary dynamos driven by helical waves?, *J. Plasma Phys.*, **84**(3), doi:10.1017/S0022377818000466.

Dietrich, W., Schmitt, D. & Wicht, J., 2013. Hemispherical Parker waves driven by thermal shear in planetary dynamos, *Europhys. Lett.*, **104**(4), doi:10.1209/0295-5075/104/49001.

Dormy, E., 2016. Strong-field spherical dynamos, *J. Fluid Mech.*, **789**, 500–513.

Duarte, L.D., Gastine, T. & Wicht, J., 2013. Anelastic dynamo models with variable electrical conductivity: an application to gas giants, *Phys. Earth planet. Inter.*, **222**, 22–34.

Gastine, T., 2019. pizza: an open-source pseudo-spectral code for spherical quasi-geostrophic convection, *Geophys. J. Int.*, **217**(3), 1558–1576.

Gastine, T., Duarte, L. & Wicht, J., 2012. Dipolar versus multipolar dynamos: the influence of the background density stratification, *Astron. Astrophys.*, **546**, A19, doi:10.1051/0004-6361/201219799.

Gerick, F., Jault, D., Noir, J. & Vidal, J., 2020. Pressure torque of torsional Alfvén modes acting on an ellipsoidal mantle, *Geophys. J. Int.*, **222**(1), 338–351.

Gerick, F., Jault, D. & Noir, J., 2021. Fast quasi-geostrophic magneto-coriolis modes in the Earth's core, *Geophys. Res. Lett.*, **48**(4), e2020GL090803, doi:10.1029/2020GL090803.

Gillet, N. & Jones, C., 2006. The quasi-geostrophic model for rapidly rotating spherical convection outside the tangent cylinder, *J. Fluid Mech.*, **554**, 343–369.

Gillet, N., Gerick, F., Jault, D., Schwaiger, T., Aubert, J. & Istas, M., 2022. Satellite magnetic data reveal interannual waves in Earth's core, *Proc. Natl. Acad. Sci.*, **119**, e2115258119, doi:10.1073/pnas.2115258119.

Glatzmaier, G.A. & Roberts, P.H., 1995. A three-dimensional convective dynamo solution with rotating and finitely conducting inner core and mantle, *Phys. Earth planet. Int.*, **91**(1-3), 63–75.

Goudard, L. & Dormy, E., 2008. Relations between the dynamo region geometry and the magnetic behavior of stars and planets, *Europhysics Letters*, **83**(5), doi:10.1209/0295-5075/83/59001.

Greenspan, H.P.G. et al., 1968. *The Theory of Rotating Fluids*, CUP Archive.

Grote, E. & Busse, F.H., 2000. Hemispherical dynamos generated by convection in rotating spherical shells, *Phys. Rev. E*, **62**(3), 4457–4460.

Guervilly, C., 2010. Dynamos numériques planétaires générées par cisaille-
ment en surface ou chauffage interne, *Thesis*, Université de Grenoble,
France.

Guervilly, C. & Cardin, P., 2016. Subcritical convection of liquid metals in a rotating sphere using a quasi-geostrophic model, *J. Fluid Mech.*, **808**, 61–89.

Guervilly, C. & Cardin, P., 2017. Multiple zonal jets and convective heat transport barriers in a quasi-geostrophic model of planetary cores, *Geophys. J. Int.*, **211**(1), 455–471.

Guervilly, C., Cardin, P. & Schaeffer, N., 2019. Turbulent convective length scale in planetary cores, *Nature*, **570**(7761), 368–371.

Hollerbach, R., 1996. On the theory of the geodynamo, *Phys. Earth planet. Inter.*, **98**(3-4), 163–185.

Jackson, A. & Maffei, S., 2020. Plesio-geostrophy for Earth's core: I. basic equations, inertial modes and induction, *Proc. R. Soc. A*, **476**(2243), doi:10.1098/rspa.2020.0513.

Jones, C.A., 2008. Course 2 dynamo theory, in *Les Houches*, Vol. 88, pp. 45–135, Elsevier.

Krause, F. & Radler, K., 1980. *Mean-Field Magnetohydrodynamics and Dynamo Theory*, Akadernie-Verlag.

Labbé, F., Jault, D. & Gillet, N., 2015. On magnetostrophic inertia-less waves in quasi-geostrophic models of planetary cores, *Geophys. Astrophys. Fluid Dyn.*, **109**(6), 587–610.

Maffei, S., Jackson, A. & Livermore, P.W., 2017. Characterization of colum-
nar inertial modes in rapidly rotating spheres and spheroids, *Proc. R. Soc.
A*, **473**(2204), doi:10.1098/rspa.2017.0181.

Malkus, W. & Proctor, M., 1975. The macrodynamics of α-effect dynamos in rotating fluids, *J. Fluid Mech.*, **67**(3), 417–443.

Moffatt, H.K., 1978. Magnetic Field Generation in Electrically Conducting Fluids: Cambridge Monographs on Mechanics and Applied Mathematics, Cambridge Univ. Press.

Nataf, H.-C. & Schaeffer, N., 2015. Turbulence in the core, in *Treatise on Geophysics*, 2nd edn, Vol. 8, pp. 161–181, ed. Schubert, G., Elsevier.
DATA AVAILABILITY

All our data are available upon reasonable request to the corre-
sponding author and some key parameters from our whole data-
set are already included in this paper. The codes used (pizza) is
freely available (at http://www.github.com/magic-sph/pizza/tree/hy
brid_QG-3D) under the GNU GPL v3 license.

REFERENCES

Aubert, J., 2018. Geomagnetic acceleration and rapid hydromagnetic wave dynamics in advanced numerical simulations of the geodynamo, *Geophys. J. Int.*, **214**(1), 531–547.

Aubert, J., 2019. Approaching earth's core conditions in high-resolution geodynamo simulations, *Geophys. J. Int.*, **219**(Suppl 1), S137-S151.

Aubert, J., 2023. State and evolution of the geodynamo from numerical models reaching the physical conditions of Earth's core, *Geophys. J. Int.*, **235**(1), 468–487.

Aubert, J., Gillet, N. & Cardin, P., 2003. Quasigeostrophic models of convection in rotating spherical shells, *Geochem. Geophys. Geosyst.*, **4**(7), 1–19.

Aubert, J., Finlay, C. C. & Fournier, A., 2013. Bottom-up control of geomagnetic secular variation by the Earth's inner core, *Nature*, **502**(7470), 219–223.

Aubert, J., Gastine, T. & Fournier, A., 2017. Spherical convective dynamos in the rapidly rotating asymptotic regime, *J. Fluid Mech.*, **813**, 558–593.

Aubert, J., Livermore, P.W., Finlay, C.C., Fournier, A. & Gillet, N., 2022. A taxonomy of simulated geomagnetic jerks, *Geophys. J. Int.*, **231**(1), 650–672.

Aurnou, J. & King, E., 2017. The cross-over to magnetostrophic convection in planetary dynamo systems, *Proc. R. Soc. A*, **473**(2199), doi:10.1098/rspa.2016.0731.

Barrois, O., Gastine, T. & Finlay, C., 2022. Comparison of quasi-geostrophic, hybrid and 3-D models of planetary core convection, *Geophys. J. Int.*, **231**(1), 129–158.

Brandenburg, A. & Subramanian, K., 2005. Astrophysical magnetic fields and nonlinear dynamo theory, *Phys. Rep.*, **417**(1-4), 1-209.

Brown, B.P., Browning, M.K., Brun, A.S., Miesch, M.S. & Toomre, J., 2010. Persistent magnetic wreaths in a rapidly rotating sun, *Astrophys. J.*, **711**(1), 424–438.

Busse, F.H., 1970. Thermal instabilities in rapidly rotating systems, *J. Fluid Mech.*, **44**(3), 441–460.

Busse, F.H., 1975. A model of the geodynamo, *Geophys. J. Int.*, **42**(2), 437–459.

Busse, F.H. & Simitev, R.D., 2006. Parameter dependences of convection-
driven dynamos in rotating spherical fluid shells, *Geophys. astrophys.
Fluid. Dyn.*, **100**(4–5), 341–361.

Calkins, M.A., 2018. Quasi-geostrophic dynamo theory, *Phys. Earth planet. Inter.*, **276**, 182–189.

Canet, E., Finlay, C.C. & Fournier, A., 2014. Hydromagnetic quasi-
geostrophic modes in rapidly rotating planetary cores, *Phys. Earth planet.
Inter.*, **229**, 1–15.

Cardin, P. & Olson, P., 1994. Chaotic thermal convection in a rapidly rotating spherical shell: consequences for flow in the outer core, *Phys. Earth planet. Inter.*, **82**(3-4), 235-259.

Chan, K.H., Zhang, K., Zou, J. & Schubert, G., 2001. A non-linear, 3-D spherical α2 dynamo using a finite element method, *Phys. Earth planet. Inter.*, **128**(1-4), 35–50.

Christensen, U.R. & Aubert, J., 2006. Scaling properties of convection-
driven dynamos in rotating spherical shells and application to planetary
magnetic fields, *Geophys. J. Int.*, **166**(1), 97–114.

Davidson, P., 2013. Scaling laws for planetary dynamos, *Geophys. J. Int.*, **195**(1), 67–74.

Davidson, P. & Ranjan, A., 2015. Planetary dynamos driven by helical waves
— II, *Geophys. J. Int.*, **202**(3), 1646–1662.

Davidson, P.A. & Ranjan, A., 2018. Are planetary dynamos driven by helical waves?, *J. Plasma Phys.*, **84**(3), doi:10.1017/S0022377818000466.

Dietrich, W., Schmitt, D. & Wicht, J., 2013. Hemispherical Parker waves driven by thermal shear in planetary dynamos, *Europhys. Lett.*, **104**(4), doi:10.1209/0295-5075/104/49001.

Dormy, E., 2016. Strong-field spherical dynamos, *J. Fluid Mech.*, **789**, 500–513.

Duarte, L.D., Gastine, T. & Wicht, J., 2013. Anelastic dynamo models with variable electrical conductivity: an application to gas giants, *Phys. Earth planet. Inter.*, **222**, 22–34.

Gastine, T., 2019. pizza: an open-source pseudo-spectral code for spherical quasi-geostrophic convection, *Geophys. J. Int.*, **217**(3), 1558–1576.

Gastine, T., Duarte, L. & Wicht, J., 2012. Dipolar versus multipolar dynamos: the influence of the background density stratification, *Astron. Astrophys.*, **546**, A19, doi:10.1051/0004-6361/201219799.

Gerick, F., Jault, D., Noir, J. & Vidal, J., 2020. Pressure torque of torsional Alfvén modes acting on an ellipsoidal mantle, *Geophys. J. Int.*, **222**(1), 338–351.

Gerick, F., Jault, D. & Noir, J., 2021. Fast quasi-geostrophic magneto-coriolis modes in the Earth's core, *Geophys. Res. Lett.*, **48**(4), e2020GL090803, doi:10.1029/2020GL090803.

Gillet, N. & Jones, C., 2006. The quasi-geostrophic model for rapidly rotating spherical convection outside the tangent cylinder, *J. Fluid Mech.*, **554**, 343–369.

Gillet, N., Gerick, F., Jault, D., Schwaiger, T., Aubert, J. & Istas, M., 2022. Satellite magnetic data reveal interannual waves in Earth's core, *Proc. Natl. Acad. Sci.*, **119**, e2115258119, doi:10.1073/pnas.2115258119.

Glatzmaier, G.A. & Roberts, P.H., 1995. A three-dimensional convective dynamo solution with rotating and finitely conducting inner core and mantle, *Phys. Earth planet. Int.*, **91**(1-3), 63–75.

Goudard, L. & Dormy, E., 2008. Relations between the dynamo region geometry and the magnetic behavior of stars and planets, *Europhysics Letters*, **83**(5), doi:10.1209/0295-5075/83/59001.

Greenspan, H.P.G. et al., 1968. *The Theory of Rotating Fluids*, CUP Archive.

Grote, E. & Busse, F.H., 2000. Hemispherical dynamos generated by convection in rotating spherical shells, *Phys. Rev. E*, **62**(3), 4457–4460.

Guervilly, C., 2010. Dynamos numériques planétaires générées par cisaille-
ment en surface ou chauffage interne, *Thesis*, Université de Grenoble,
France.

Guervilly, C. & Cardin, P., 2016. Subcritical convection of liquid metals in a rotating sphere using a quasi-geostrophic model, *J. Fluid Mech.*, **808**, 61–89.

Guervilly, C. & Cardin, P., 2017. Multiple zonal jets and convective heat transport barriers in a quasi-geostrophic model of planetary cores, *Geophys. J. Int.*, **211**(1), 455–471.

Guervilly, C., Cardin, P. & Schaeffer, N., 2019. Turbulent convective length scale in planetary cores, *Nature*, **570**(7761), 368–371.

Hollerbach, R., 1996. On the theory of the geodynamo, *Phys. Earth planet. Inter.*, **98**(3-4), 163–185.

Jackson, A. & Maffei, S., 2020. Plesio-geostrophy for Earth's core: I. basic equations, inertial modes and induction, *Proc. R. Soc. A*, **476**(2243), doi:10.1098/rspa.2020.0513.

Jones, C.A., 2008. Course 2 dynamo theory, in *Les Houches*, Vol. **88**, pp. 45–135, Elsevier.

Krause, F. & Radler, K., 1980. *Mean-Field Magnetohydrodynamics and Dynamo Theory*, Akadernie-Verlag.

Labbé, F., Jault, D. & Gillet, N., 2015. On magnetostrophic inertia-less waves in quasi-geostrophic models of planetary cores, *Geophys. Astrophys. Fluid Dyn.*, **109**(6), 587–610.

Maffei, S., Jackson, A. & Livermore, P.W., 2017. Characterization of colum-
nar inertial modes in rapidly rotating spheres and spheroids, *Proc. R. Soc.
A*, **473**(2204), doi:10.1098/rspa.2017.0181.

Malkus, W. & Proctor, M., 1975. The macrodynamics of α-effect dynamos in rotating fluids, *J. Fluid Mech.*, **67**(3), 417–443.

Moffatt, H.K., 1978. Magnetic Field Generation in Electrically Conducting Fluids: Cambridge Monographs on Mechanics and Applied Mathematics, Cambridge Univ. Press.

Nataf, H.-C. & Schaeffer, N., 2015. Turbulence in the core, in *Treatise on Geophysics*, 2nd edn, Vol. 8, pp. 161–181, ed. Schubert, G., Elsevier.
Olson, P., Christensen, U. & Glatzmaier, G.A., 1999. Numerical modeling of the geodynamo: mechanisms of field generation and equilibration, *J. geophys. Res.*, **104**(B5), 10 383–10 404.

Parker, E.N., 1955. Hydromagnetic dynamo models, *Astrophys. J.*, **122**, 293, doi:10.1086/146087.

Petitdemange, L., 2018. Systematic parameter study of dynamo bifurcations in geodynamo simulations, *Phys. Earth planet. Inter.*, **277**, 113–132.

Ranjan, A., Davidson, P., Christensen, U. & Wicht, J., 2020. On the gener-
ation and segregation of helicity in geodynamo simulations, *Geophys. J.
Int.*, **221**, 741–757.

Roberts, P. & Scott, S., 1965. On analysis of the secular variation, *J. Geo-
magnet. Geoelectr.* **17**(2), 137–151.

Roberts, P.H. & King, E.M., 2013. On the genesis of the Earth’s magnetism,
*Rep. Prog. Phys.*, **76**(9), doi:10.1088/0034-4885/76/9/096801.

Schaeffer, N., 2005. *Instabilités, turbulence et dynamo dans une couche de fluide cisaillée en rotation rapide. Importance de l'aspect ondulatoire*, Université de Grenoble, France.

Schaeffer, N., 2013. Efficient spherical harmonic transforms aimed at pseu-
dospectral numerical simulations, *Geochem. Geophys. Geosyst.*, **14**(3),
751–758.

Schaeffer, N. & Cardin, P., 2005. Quasigeostrophic model of the instabilities of the Stewartson layer in flat and depth-varying containers, *Phys. Fluids*, **17**(10), doi:10.1063/1.2073547.

Schaeffer, N. & Cardin, P., 2006. Quasi-geostrophic kinematic dynamos at low magnetic Prandtl number, *Earth planet. Sci. Lett.*, **245**(3-4), 595-604.

Schaeffer, N., Silva, E.L. & Pais, M.A., 2016. Can core flows inferred from geomagnetic field models explain the Earth's dynamo?, *Geophys. J. Int.*, **204**(2), 868–877.

Schaeffer, N., Jault, D., Nataf, H.-C. & Fournier, A., 2017. Turbulent geody-
namo simulations: a leap towards Earth's core, *Geophys. J. Int.*, **211**(1),
1–29.

Schrinner, M., Rädler, K.-H., Schmitt, D., Rheinhardt, M. & Christensen, U.R., 2007. Mean-field concept and direct numerical simulations of rotating magnetoconvection and the geodynamo, *Geophys. Astro Fluid Dyn.*, **101**(2), 81–116.

Schrinner, M., Petitdemange, L. & Dormy, E., 2011. Oscillatory dy-
namos and their induction mechanisms, *Astron. Astrophys.*, **530**,
doi:10.1051/0004-6361/201016372.

Schrinner, M., Petitdemange, L. & Dormy, E., 2012. Dipole collapse and dynamo waves in global direct numerical simulations, *Astrophys. J.*, **752**(2), 121.

Schwaiger, T., Gastine, T. & Aubert, J., 2019. Force balance in numerical geodynamo simulations: a systematic study, *Geophys. J. Int.*, **219**(Suppl 1), S101–S114.

Sheyko, A., Finlay, C.C. & Jackson, A., 2016. Magnetic reversals from planetary dynamo waves, *Nature*, **539**(7630), 551–554.

Sheyko, A., Finlay, C., Favre, J. & Jackson, A., 2018. Scale separated low viscosity dynamos and dissipation within the earth's core, *Nat. Sci. Rep.*, **8**(1), 1–7.

Simitev, R.D. & Busse, F.H., 2012. How far can minimal models explain the solar cycle?, *Astrophys. J.* **749**(1), doi:10.1088/0004-637X/749/1/9.

Sreenivasan, B. & Jones, C.A., 2011. Helicity generation and subcritical behaviour in rapidly rotating dynamos, *J Fluid Mech.*, **688**, 5–30.

Sreenivasan, B. & Kar, S., 2018. Scale dependence of kinetic helicity and selection of the axial dipole in rapidly rotating dynamos, *Phys. Rev. Fluids*, **3**(9), doi:10.1103/PhysRevFluids.3.093801.

Steenbeck, M. & Krause, F., 1969. On the dynamo theory of stellar and planetary magnetic fields. I. AC dynamos of solar type, *Astronomische Nachrichten*, **291**, 49–84.

Steenbeck, M., Krause, F. & Rädler, K.-H., 1966. Berechnung der mit-
tleren lorentz-feldstärke für ein elektrisch leitendes medium in turbulen-
ter, durch coriolis-kräfte beeinflußter bewegung, *Zeitschrift für Natur-
forschung A*, **21**(4), 369–376.

Varma, A. & Sreenivasan, B., 2022. The role of slow magnetostrophic waves in the formation of the axial dipole in planetary dynamos, *Phys. Earth planet. Inter.*, **333**, doi:10.1016/j.pepi.2022.106944.

Wicht, J., 2002. Inner-core conductivity in numerical dynamo simulations,
*Phys. Earth planet. Inter.*, **132**(4), 281–302.

Yoshimura, H., 1975. Solar-cycle dynamo wave propagation, *Astrophys. J.*,
**201**, 740–748.
Olson, P., Christensen, U. & Glatzmaier, G.A., 1999. Numerical modeling of the geodynamo: mechanisms of field generation and equilibration, *J. geophys. Res.*, **104**(B5), 10 383–10 404.

Parker, E.N., 1955. Hydromagnetic dynamo models, *Astrophys. J.*, **122**, 293, doi:10.1086/146087.

Petitdemange, L., 2018. Systematic parameter study of dynamo bifurcations in geodynamo simulations, *Phys. Earth planet. Inter.*, **277**, 113–132.

Ranjan, A., Davidson, P., Christensen, U. & Wicht, J., 2020. On the gener-
ation and segregation of helicity in geodynamo simulations, *Geophys. J.*
*Int.*, **221**, 741–757.

Roberts, P. & Scott, S., 1965. On analysis of the secular variation, *J. Geo-magnet. Geoelectr.* **17**(2), 137–151.

Roberts, P.H. & King, E.M., 2013. On the genesis of the Earth's magnetism,
*Rep. Prog. Phys.*, **76**(9), doi:10.1088/0034-4885/76/9/096801.

Schaeffer, N., 2005. *Instabilités, turbulence et dynamo dans une couche de fluide cisaillée en rotation rapide. Importance de l'aspect ondulatoire*, Université de Grenoble, France.

Schaeffer, N., 2013. Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations, *Geochem. Geophys. Geosyst.*, **14**(3), 751–758.

Schaeffer, N. & Cardin, P., 2005. Quasigeostrophic model of the instabilities of the Stewartson layer in flat and depth-varying containers, *Phys. Fluids*, **17**(10), doi:10.1063/1.2073547.

Schaeffer, N. & Cardin, P., 2006. Quasi-geostrophic kinematic dynamos at low magnetic Prandtl number, *Earth planet. Sci. Lett.*, **245**(3-4), 595-604.

Schaeffer, N., Silva, E.L. & Pais, M.A., 2016. Can core flows inferred from geomagnetic field models explain the Earth's dynamo?, *Geophys. J. Int.*, **204**(2), 868–877.

Schaeffer, N., Jault, D., Nataf, H.-C. & Fournier, A., 2017. Turbulent geody-
namo simulations: a leap towards Earth's core, *Geophys. J. Int.*, **211**(1),
1–29.

Schrinner, M., Rädler, K.-H., Schmitt, D., Rheinhardt, M. & Christensen, U.R., 2007. Mean-field concept and direct numerical simulations of rotating magnetoconvection and the geodynamo, *Geophys. Astro Fluid Dyn.*, **101**(2), 81–116.

Schrinner, M., Petitdemange, L. & Dormy, E., 2011. Oscillatory dynamos and their induction mechanisms, *Astron. Astrophys.*, **530**, doi:10.1051/0004-6361/201016372.

Schrinner, M., Petitdemange, L. & Dormy, E., 2012. Dipole collapse and dynamo waves in global direct numerical simulations, *Astrophys. J.*, **752**(2), 121.

Schwaiger, T., Gastine, T. & Aubert, J., 2019. Force balance in numerical geodynamo simulations: a systematic study, *Geophys. J. Int.*, **219**(Suppl 1), S101–S114.

Sheyko, A., Finlay, C.C. & Jackson, A., 2016. Magnetic reversals from planetary dynamo waves, *Nature*, **539**(7630), 551–554.

Sheyko, A., Finlay, C., Favre, J. & Jackson, A., 2018. Scale separated low viscosity dynamos and dissipation within the earth's core, *Nat. Sci. Rep.*, **8**(1), 1–7.

Simitev, R.D. & Busse, F.H., 2012. How far can minimal models explain the solar cycle?, *Astrophys. J.* **749**(1), doi:10.1088/0004-637X/749/1/9.

Sreenivasan, B. & Jones, C.A., 2011. Helicity generation and subcritical behaviour in rapidly rotating dynamos, *J Fluid Mech.*, **688**, 5–30.

Sreenivasan, B. & Kar, S., 2018. Scale dependence of kinetic helicity and selection of the axial dipole in rapidly rotating dynamos, *Phys. Rev. Fluids*, **3**(9), doi:10.1103/PhysRevFluids.3.093801.

Steenbeck, M. & Krause, F., 1969. On the dynamo theory of stellar and planetary magnetic fields. I. AC dynamos of solar type, *Astronomische Nachrichten*, **291**, 49–84.

Steenbeck, M., Krause, F. & Rädler, K.-H., 1966. Berechnung der mittleren lorentz-feldstärke für ein elektrisch leitendes medium in turbulenter, durch coriolis-kräfte beeinflußter bewegung, *Zeitschrift für Naturforschung A*, **21**(4), 369–376.

Varma, A. & Sreenivasan, B., 2022. The role of slow magnetostrophic waves in the formation of the axial dipole in planetary dynamos, *Phys. Earth planet. Inter.*, **333**, doi:10.1016/j.pepi.2022.106944.

Wicht, J., 2002. Inner-core conductivity in numerical dynamo simulations,
*Phys. Earth planet. Inter.*, **132**(4), 281–302.

Yoshimura, H., 1975. Solar-cycle dynamo wave propagation, *Astrophys. J.*,
**201**, 740–748.
# APPENDIX A: RESULTS OF NUMERICAL SIMULATIONS

<table><caption>**Table A1.** Summary of the hybrid QG-3D numerical simulations computed in this study at $Ek = 3 \times 10^{-5}$, $Pr = 0.1$, using $\eta = r_i/r_0 = 0.35$. $Ra$ is the Rayleigh number (the supercriticality $Sc = Ra/Ra_c$ is also provided, where $Ra_c$ is the thermal convection critical value), $Pm$ is the magnetic Prandtl number, $Nu$ is the Nusselt number, $Rm$ is the magnetic Reynolds number, $\Lambda$ is the Elsasser number, $M = \hat{E}_{mag}/\hat{E}_{kin}$ is the magnetic to kinetic energy ratio, $f_{tor} = \hat{E}_{tor}/\hat{E}_{mag}$ is the toroidal to the total magnetic field energy ratio, $\ell_H$ and $m_H$ are the cut-off degree and azimuthal wavenumber below which the hyperdiffusion has no effect, and $(N_s, N_m)/(N_r, \ell_{max})$ are the grid sizes of the run. The simulations with a * symbol in the second column are the growing dynamos we found.</caption><thead><tr><th>Ra</th><th>Sc</th><th>Pm</th><th>Nu</th><th>Rm</th><th>Λ</th><th>M</th><th>f<sub>tor</sub></th><th>m<sub>H</sub>/ℓ<sub>H</sub></th><th>(N<sub>s</sub>, N<sub>m</sub>)/(N<sub>r</sub>, ℓ<sub>max</sub>)</th></tr></thead><tbody><tr><td>5.00 × 10<sup>6</sup></td><td>4.9</td><td>0.9</td><td>1.19</td><td>334.0</td><td>2.08 × 10<sup>-6</sup></td><td>4.86 × 10<sup>-7</sup></td><td>0.81</td><td>–/–</td><td>(193, 192)/(145, 144)</td></tr><tr><td>5.00 × 10<sup>6</sup></td><td>4.9</td><td>1.3</td><td>1.19</td><td>482.9</td><td>8.61 × 10<sup>-9</sup></td><td>1.39 × 10<sup>-9</sup></td><td>0.79</td><td>144/–</td><td>(193, 192)/(145, 144)</td></tr><tr><td>5.00 × 10<sup>6</sup></td><td>4.9</td><td>2.0</td><td>1.19</td><td>740.4</td><td>1.66 × 10<sup>-6</sup></td><td>1.76 × 10<sup>-7</sup></td><td>0.82</td><td>144/–</td><td>(385, 768)/(149, 148)</td></tr><tr><td>1.00 × 10<sup>7</sup></td><td>9.7</td><td>0.5</td><td>1.32</td><td>374.7</td><td>1.09 × 10<sup>-6</sup></td><td>1.13 × 10<sup>-7</sup></td><td>0.77</td><td>144/–</td><td>(385, 768)/(149, 148)</td></tr><tr><td>1.00 × 10<sup>7</sup></td><td>9.7</td><td>0.9</td><td>1.32</td><td>670.9</td><td>5.41 × 10<sup>-6</sup></td><td>3.15 × 10<sup>-7</sup></td><td>0.89</td><td>144/–</td><td>(385, 768)/(149, 148)</td></tr><tr><td>1.00 × 10<sup>7</sup></td><td>9.7</td><td>1.3</td><td>1.32</td><td>976.1</td><td>8.97 × 10<sup>-3</sup></td><td>3.56 × 10<sup>-4</sup></td><td>0.97</td><td>144/–</td><td>(385, 768)/(149, 148)</td></tr><tr><td>1.00 × 10<sup>7</sup></td><td>9.7*</td><td>2.0</td><td>1.32</td><td>1456.4</td><td>2.48 × 10<sup>-1</sup></td><td>6.79 × 10<sup>-3</sup></td><td>0.98</td><td>–/–</td><td>(145, 144)/(145, 144)</td></tr><tr><td>1.65 × 10<sup>7</sup></td><td>16.0*</td><td>0.9</td><td>1.40</td><td>894.3</td><td>1.31 × 10<sup>0</sup></td><td>4.28 × 10<sup>-2</sup></td><td>0.96</td><td>–/–</td><td>(129, 192)/(129, 192)</td></tr><tr><td>1.66 × 10<sup>7</sup></td><td>16.1</td><td>0.3</td><td>1.42</td><td>347.0</td><td>4.80 × 10<sup>-7</sup></td><td>3.48 × 10<sup>-8</sup></td><td>0.78</td><td>144/–</td><td>(385, 768)/(149, 148)</td></tr><tr><td>1.66 × 10<sup>7</sup></td><td>16.1</td><td>0.5</td><td>1.42</td><td>578.2</td><td>7.31 × 10<sup>-5</sup></td><td>3.18 × 10<sup>-6</sup></td><td>0.95</td><td>144/–</td><td>(385, 768)/(149, 148)</td></tr><tr><td>1.66 × 10<sup>7</sup></td><td>16.1*</td><td>0.9</td><td>1.41</td><td>909.1</td><td>1.27 × 10<sup>0</sup></td><td>4.03 × 10<sup>-2</sup></td><td>0.96</td><td>144/–</td><td>(385, 768)/(149, 148)</td></tr><tr><td>1.66 × 10<sup>7</sup></td><td>16.1*</td><td>1.1</td><td>1.41</td><td>1029.7</td><td>2.51 × 10<sup>0</sup></td><td>7.55 × 10<sup>-2</sup></td><td>0.96</td><td>–/–</td><td>(145, 192)/(145, 192)</td></tr><tr><td>1.66 × 10<sup>7</sup></td><td>16.1*</td><td>1.3</td><td>1.41</td><td>1174.0</td><td>3.68 × 10<sup>0</sup></td><td>1.00 × 10<sup>-1</sup></td><td>0.96</td><td>–/–</td><td>(145, 192)/(145, 192)</td></tr><tr><td>1.66 × 10<sup>7</sup></td><td>16.1*</td><td>2.0</td><td>1.40</td><td>1747.4</td><td>8.02 × 10<sup>0</sup></td><td>1.51 × 10<sup>-1</sup></td><td>0.95</td><td>288/–</td><td>(289, 416)/(193, 288)</td></tr><tr><td>2.00 × 10<sup>7</sup></td><td>19.4</td><td>0.1</td><td>1.46</td><td>135.6</td><td>3.58 × 10<sup>-9</sup></td><td>5.67 × 10<sup>-10</sup></td><td>0.69</td><td>144/–</td><td>(385, 768)/(149, 148)</td></tr><tr><td>2.00 × 10<sup>7</sup></td><td>19.4</td><td>0.3</td><td>1.46</td><td>407.1</td><td>8.20 × 10<sup>-6</sup></td><td>4.32 × 10<sup>-7</sup></td><td>0.89</td><td>144/–</td><td>(385, 768)/(149, 148)</td></tr><tr><td>2.00 × 10<sup>7</sup></td><td>19.4*</td><td>0.5</td><td>1.45</td><td>665.8</td><td>7.15 × 10<sup>-2</sup></td><td>2.35 × 10<sup>-3</sup></td><td>0.95</td><td>–/–</td><td>(97, 128)/(97, 128)</td></tr><tr><td>2.00 × 10<sup>7</sup></td><td>19.4*</td><td>0.9</td><td>1.44</td><td>913.0</td><td>3.46 × 10<sup>0</sup></td><td>1.08 × 10<sup>-1</sup></td><td>0.95</td><td>144/–</td><td>(145, 180)/(145, 144)</td></tr><tr><td>2.00 × 10<sup>7</sup></td><td>19.4*</td><td>1.3</td><td>1.45</td><td>1274.1</td><td>6.14 × 10<sup>0</sup></td><td>1.42 × 10<sup>-1</sup></td><td>0.96</td><td>–/–</td><td>(145, 192)/(145, 192)</td></tr><tr><td>3.00 × 10<sup>7</sup></td><td>29.1</td><td>0.1</td><td>1.55</td><td>188.9</td><td>3.86 × 10<sup>-7</sup></td><td>3.15 × 10<sup>-8</sup></td><td>0.73</td><td>144/–</td><td>(385, 768)/(149, 148)</td></tr><tr><td>2.99 × 10<sup>7</sup></td><td>29.0</td><td>0.2</td><td>1.54</td><td>375.3</td><td>2.51 × 10<sup>-7</sup></td><td>1.04 × 10<sup>-8</sup></td><td>0.77</td><td>–/–</td><td>(385, 960)/(385, 320)</td></tr><tr><td>3.00 × 10<sup>7</sup></td><td>29.1*</td><td>0.3</td><td>1.54</td><td>557.8</td><td>9.75 × 10<sup>-2</sup></td><td>2.74 × 10<sup>-3</sup></td><td>0.93</td><td>–/–</td><td>(97, 128)/(97, 128)</td></tr><tr><td>3.00 × 10<sup>7</sup></td><td>29.1*</td><td>0.5</td><td>1.52</td><td>691.2</td><td>3.16 × 10<sup>0</sup></td><td>9.58 × 10<sup>-2</sup></td><td>0.94</td><td>256/–</td><td>(433, 864)/(193, 256)</td></tr><tr><td>3.00 × 10<sup>7</sup></td><td>29.1*</td><td>0.9</td><td>1.57</td><td>1283.8</td><td>7.61 × 10<sup>0</sup></td><td>1.21 × 10<sup>-1</sup></td><td>0.92</td><td>224/224</td><td>(433, 480)/(193, 256)</td></tr><tr><td>5.00 × 10<sup>7</sup></td><td>48.5</td><td>0.1</td><td>1.66</td><td>275.7</td><td>3.84 × 10<sup>-7</sup></td><td>1.47 × 10<sup>-8</sup></td><td>0.79</td><td>–/–</td><td>(433, 800)/(149, 148)</td></tr><tr><td>5.00 × 10<sup>7</sup></td><td>48.5</td><td>0.2</td><td>1.66</td><td>550.9</td><td>6.91 × 10<sup>-6</sup></td><td>1.33 × 10<sup>-7</sup></td><td>0.84</td><td>816/–</td><td>(433, 864)/(149, 148)</td></tr><tr><td>4.99 × 10<sup>7</sup></td><td>48.4*</td><td>0.3</td><td>1.62</td><td>640.5</td><td>2.31 × 10<sup>0</sup></td><td>4.91 × 10<sup>-2</sup></td><td>0.91</td><td>576/144</td><td>(481, 912)/(193, 256)</td></tr><tr><td>8.00 × 10<sup>7</sup></td><td>77.7</td><td>0.1</td><td>1.73</td><td>376.5</td><td>7.80 × 10<sup>-7</sup></td><td>1.60 × 10<sup>-8</sup></td><td>0.83</td><td>144/–</td><td>(577, 864)/(193, 192)</td></tr><tr><td>8.00 × 10<sup>7</sup></td><td>77.7</td><td>0.2</td><td>1.73</td><td>754.4</td><td>1.07 × 10<sup>-4</sup></td><td>1.10 × 10<sup>-6</sup></td><td>0.86</td><td>144/–</td><td>(577, 864)/(193, 192)</td></tr><tr><td>8.00 × 10<sup>7</sup></td><td>77.7*</td><td>0.3</td><td>1.66</td><td>836.0</td><td>3.65 × 10<sup>0</sup></td><td>4.55 × 10<sup>-2</sup></td><td>0.92</td><td>192/–</td><td>(289, 384)/(193, 192)</td></tr><tr><td>1.03 × 10<sup>8</sup></td><td>100</td><td>0.1</td><td>1.76</td><td>439.9</td><td>4.42 × 10<sup>-5</sup></td><td>6.65 × 10<sup>-7</sup></td><td>0.84</td><td>768/–</td><td>(577, 864)/(193, 192)</td></tr><tr><td>1.03 × 10<sup>8</sup></td><td>100</td><td>0.2</td><td>1.74</td><td>861.7</td><td>9.89 × 10<sup>-3</sup></td><td>1.94 × 10<sup>-4</sup></td><td>0.94</td><td>768/–</td><td>(577, 912)/(193, 256)</td></tr><tr><td>1.03 × 10<sup>8</sup></td><td>100*</td><td>0.3</td><td>1.69</td><td>982.1</td><td>5.44 × 10<sup>0</sup></td><td>4.92 × 10<sup>-2</sup></td><td>0.94</td><td>192/–</td><td>(289, 384)/(193, 192)</td></tr></tbody></table>

# APPENDIX B: CODE VALIDATION AND BENCHMARKS

## B1 Benchmark of the Lorentz force computation

The determination of eq. (21) relies on the computation of $(\mathbf{j} \times \mathbf{B})$ which depends on $\nabla \times \mathbf{B}$ and $\mathbf{B}$. These quantities are computed using SHTns<fnref n="2" /> functions, already validated in (Schaeffer 2013) and widely used in many codes. Given any fields $(B_{pol}, B_{tor})$ in spectral space, and their radial-derivatives, the SHTns routines directly provide $(B_r, B_{\theta}, B_{\phi})$ and $(j_r, j_{\theta}, j_{\phi})$ on the 3D-physical-grid. The remaining part of the process involves the computation of the non-linear products $(\mathbf{j} \times \mathbf{B})$, and the z-averaging of $\partial_s[s(\mathbf{j} \times \mathbf{B})_\phi]$, $(\mathbf{j} \times \mathbf{B})_\phi$ and $(\mathbf{j} \times \mathbf{B})_s$ which are computed on the physical grid. These quantities are then sent to the spectral space using `fft` functions where the $\phi$-derivative is performed before assembling all the terms to obtain eq. (21).

<footnote n="2">²https://bitbucket.org/nschaeff/shtns</footnote>
# APPENDIX A: RESULTS OF NUMERICAL SIMULATIONS

<table><caption>**Table A1.** Summary of the hybrid QG-3D numerical simulations computed in this study at $Ek = 3 \times 10^{-5}$, $Pr = 0.1$, using $\eta = r_i/r_0 = 0.35$. $Ra$ is the Rayleigh number (the supercriticality $Sc = Ra/Ra_c$ is also provided, where $Ra_c$ is the thermal convection critical value), $Pm$ is the magnetic Prandtl number, $Nu$ is the Nusselt number, $Rm$ is the magnetic Reynolds number, $\Lambda$ is the Elsasser number, $M = \hat{E}_{mag}/\hat{E}_{kin}$ is the magnetic to kinetic energy ratio, $f_{tor} = \hat{E}_{tor}/\hat{E}_{mag}$ is the toroidal to the total magnetic field energy ratio, $\ell_H$ and $m_H$ are the cut-off degree and azimuthal wavenumber below which the hyperdiffusion has no effect, and $(N_s, N_m)/(N_r, \ell_{max})$ are the grid sizes of the run. The simulations with a * symbol in the second column are the growing dynamos we found.</caption><thead><tr><th>Ra</th><th>Sc</th><th>Pm</th><th>Nu</th><th>Rm</th><th>Λ</th><th>M</th><th>f<sub>tor</sub></th><th>m<sub>H</sub>/ℓ<sub>H</sub></th><th>(N<sub>s</sub>, N<sub>m</sub>)/(N<sub>r</sub>, ℓ<sub>max</sub>)</th></tr></thead><tbody><tr><td>5.00 × 10<sup>6</sup></td><td>4.9</td><td>0.9</td><td>1.19</td><td>334.0</td><td>2.08 × 10<sup>-6</sup></td><td>4.86 × 10<sup>-7</sup></td><td>0.81</td><td>–/–</td><td>(193, 192)/(145, 144)</td></tr><tr><td>5.00 × 10<sup>6</sup></td><td>4.9</td><td>1.3</td><td>1.19</td><td>482.9</td><td>8.61 × 10<sup>-9</sup></td><td>1.39 × 10<sup>-9</sup></td><td>0.79</td><td>144/–</td><td>(193, 192)/(145, 144)</td></tr><tr><td>5.00 × 10<sup>6</sup></td><td>4.9</td><td>2.0</td><td>1.19</td><td>740.4</td><td>1.66 × 10<sup>-6</sup></td><td>1.76 × 10<sup>-7</sup></td><td>0.82</td><td>144/–</td><td>(385, 768)/(149, 148)</td></tr><tr><td>1.00 × 10<sup>7</sup></td><td>9.7</td><td>0.5</td><td>1.32</td><td>374.7</td><td>1.09 × 10<sup>-6</sup></td><td>1.13 × 10<sup>-7</sup></td><td>0.77</td><td>144/–</td><td>(385, 768)/(149, 148)</td></tr><tr><td>1.00 × 10<sup>7</sup></td><td>9.7</td><td>0.9</td><td>1.32</td><td>670.9</td><td>5.41 × 10<sup>-6</sup></td><td>3.15 × 10<sup>-7</sup></td><td>0.89</td><td>144/–</td><td>(385, 768)/(149, 148)</td></tr><tr><td>1.00 × 10<sup>7</sup></td><td>9.7</td><td>1.3</td><td>1.32</td><td>976.1</td><td>8.97 × 10<sup>-3</sup></td><td>3.56 × 10<sup>-4</sup></td><td>0.97</td><td>144/–</td><td>(385, 768)/(149, 148)</td></tr><tr><td>1.00 × 10<sup>7</sup></td><td>9.7*</td><td>2.0</td><td>1.32</td><td>1456.4</td><td>2.48 × 10<sup>-1</sup></td><td>6.79 × 10<sup>-3</sup></td><td>0.98</td><td>–/–</td><td>(145, 144)/(145, 144)</td></tr><tr><td>1.65 × 10<sup>7</sup></td><td>16.0*</td><td>0.9</td><td>1.40</td><td>894.3</td><td>1.31 × 10<sup>0</sup></td><td>4.28 × 10<sup>-2</sup></td><td>0.96</td><td>–/–</td><td>(129, 192)/(129, 192)</td></tr><tr><td>1.66 × 10<sup>7</sup></td><td>16.1</td><td>0.3</td><td>1.42</td><td>347.0</td><td>4.80 × 10<sup>-7</sup></td><td>3.48 × 10<sup>-8</sup></td><td>0.78</td><td>144/–</td><td>(385, 768)/(149, 148)</td></tr><tr><td>1.66 × 10<sup>7</sup></td><td>16.1</td><td>0.5</td><td>1.42</td><td>578.2</td><td>7.31 × 10<sup>-5</sup></td><td>3.18 × 10<sup>-6</sup></td><td>0.95</td><td>144/–</td><td>(385, 768)/(149, 148)</td></tr><tr><td>1.66 × 10<sup>7</sup></td><td>16.1*</td><td>0.9</td><td>1.41</td><td>909.1</td><td>1.27 × 10<sup>0</sup></td><td>4.03 × 10<sup>-2</sup></td><td>0.96</td><td>144/–</td><td>(385, 768)/(149, 148)</td></tr><tr><td>1.66 × 10<sup>7</sup></td><td>16.1*</td><td>1.1</td><td>1.41</td><td>1029.7</td><td>2.51 × 10<sup>0</sup></td><td>7.55 × 10<sup>-2</sup></td><td>0.96</td><td>–/–</td><td>(145, 192)/(145, 192)</td></tr><tr><td>1.66 × 10<sup>7</sup></td><td>16.1*</td><td>1.3</td><td>1.41</td><td>1174.0</td><td>3.68 × 10<sup>0</sup></td><td>1.00 × 10<sup>-1</sup></td><td>0.96</td><td>–/–</td><td>(145, 192)/(145, 192)</td></tr><tr><td>1.66 × 10<sup>7</sup></td><td>16.1*</td><td>2.0</td><td>1.40</td><td>1747.4</td><td>8.02 × 10<sup>0</sup></td><td>1.51 × 10<sup>-1</sup></td><td>0.95</td><td>288/–</td><td>(289, 416)/(193, 288)</td></tr><tr><td>2.00 × 10<sup>7</sup></td><td>19.4</td><td>0.1</td><td>1.46</td><td>135.6</td><td>3.58 × 10<sup>-9</sup></td><td>5.67 × 10<sup>-10</sup></td><td>0.69</td><td>144/–</td><td>(385, 768)/(149, 148)</td></tr><tr><td>2.00 × 10<sup>7</sup></td><td>19.4</td><td>0.3</td><td>1.46</td><td>407.1</td><td>8.20 × 10<sup>-6</sup></td><td>4.32 × 10<sup>-7</sup></td><td>0.89</td><td>144/–</td><td>(385, 768)/(149, 148)</td></tr><tr><td>2.00 × 10<sup>7</sup></td><td>19.4*</td><td>0.5</td><td>1.45</td><td>665.8</td><td>7.15 × 10<sup>-2</sup></td><td>2.35 × 10<sup>-3</sup></td><td>0.95</td><td>–/–</td><td>(97, 128)/(97, 128)</td></tr><tr><td>2.00 × 10<sup>7</sup></td><td>19.4*</td><td>0.9</td><td>1.44</td><td>913.0</td><td>3.46 × 10<sup>0</sup></td><td>1.08 × 10<sup>-1</sup></td><td>0.95</td><td>144/–</td><td>(145, 180)/(145, 144)</td></tr><tr><td>2.00 × 10<sup>7</sup></td><td>19.4*</td><td>1.3</td><td>1.45</td><td>1274.1</td><td>6.14 × 10<sup>0</sup></td><td>1.42 × 10<sup>-1</sup></td><td>0.96</td><td>–/–</td><td>(145, 192)/(145, 192)</td></tr><tr><td>3.00 × 10<sup>7</sup></td><td>29.1</td><td>0.1</td><td>1.55</td><td>188.9</td><td>3.86 × 10<sup>-7</sup></td><td>3.15 × 10<sup>-8</sup></td><td>0.73</td><td>144/–</td><td>(385, 768)/(149, 148)</td></tr><tr><td>2.99 × 10<sup>7</sup></td><td>29.0</td><td>0.2</td><td>1.54</td><td>375.3</td><td>2.51 × 10<sup>-7</sup></td><td>1.04 × 10<sup>-8</sup></td><td>0.77</td><td>–/–</td><td>(385, 960)/(385, 320)</td></tr><tr><td>3.00 × 10<sup>7</sup></td><td>29.1*</td><td>0.3</td><td>1.54</td><td>557.8</td><td>9.75 × 10<sup>-2</sup></td><td>2.74 × 10<sup>-3</sup></td><td>0.93</td><td>–/–</td><td>(97, 128)/(97, 128)</td></tr><tr><td>3.00 × 10<sup>7</sup></td><td>29.1*</td><td>0.5</td><td>1.52</td><td>691.2</td><td>3.16 × 10<sup>0</sup></td><td>9.58 × 10<sup>-2</sup></td><td>0.94</td><td>256/–</td><td>(433, 864)/(193, 256)</td></tr><tr><td>3.00 × 10<sup>7</sup></td><td>29.1*</td><td>0.9</td><td>1.57</td><td>1283.8</td><td>7.61 × 10<sup>0</sup></td><td>1.21 × 10<sup>-1</sup></td><td>0.92</td><td>224/224</td><td>(433, 480)/(193, 256)</td></tr><tr><td>5.00 × 10<sup>7</sup></td><td>48.5</td><td>0.1</td><td>1.66</td><td>275.7</td><td>3.84 × 10<sup>-7</sup></td><td>1.47 × 10<sup>-8</sup></td><td>0.79</td><td>–/–</td><td>(433, 800)/(149, 148)</td></tr><tr><td>5.00 × 10<sup>7</sup></td><td>48.5</td><td>0.2</td><td>1.66</td><td>550.9</td><td>6.91 × 10<sup>-6</sup></td><td>1.33 × 10<sup>-7</sup></td><td>0.84</td><td>816/–</td><td>(433, 864)/(149, 148)</td></tr><tr><td>4.99 × 10<sup>7</sup></td><td>48.4*</td><td>0.3</td><td>1.62</td><td>640.5</td><td>2.31 × 10<sup>0</sup></td><td>4.91 × 10<sup>-2</sup></td><td>0.91</td><td>576/144</td><td>(481, 912)/(193, 256)</td></tr><tr><td>8.00 × 10<sup>7</sup></td><td>77.7</td><td>0.1</td><td>1.73</td><td>376.5</td><td>7.80 × 10<sup>-7</sup></td><td>1.60 × 10<sup>-8</sup></td><td>0.83</td><td>144/–</td><td>(577, 864)/(193, 192)</td></tr><tr><td>8.00 × 10<sup>7</sup></td><td>77.7</td><td>0.2</td><td>1.73</td><td>754.4</td><td>1.07 × 10<sup>-4</sup></td><td>1.10 × 10<sup>-6</sup></td><td>0.86</td><td>144/–</td><td>(577, 864)/(193, 192)</td></tr><tr><td>8.00 × 10<sup>7</sup></td><td>77.7*</td><td>0.3</td><td>1.66</td><td>836.0</td><td>3.65 × 10<sup>0</sup></td><td>4.55 × 10<sup>-2</sup></td><td>0.92</td><td>192/–</td><td>(289, 384)/(193, 192)</td></tr><tr><td>1.03 × 10<sup>8</sup></td><td>100</td><td>0.1</td><td>1.76</td><td>439.9</td><td>4.42 × 10<sup>-5</sup></td><td>6.65 × 10<sup>-7</sup></td><td>0.84</td><td>768/–</td><td>(577, 864)/(193, 192)</td></tr><tr><td>1.03 × 10<sup>8</sup></td><td>100</td><td>0.2</td><td>1.74</td><td>861.7</td><td>9.89 × 10<sup>-3</sup></td><td>1.94 × 10<sup>-4</sup></td><td>0.94</td><td>768/–</td><td>(577, 912)/(193, 256)</td></tr><tr><td>1.03 × 10<sup>8</sup></td><td>100*</td><td>0.3</td><td>1.69</td><td>982.1</td><td>5.44 × 10<sup>0</sup></td><td>4.92 × 10<sup>-2</sup></td><td>0.94</td><td>192/–</td><td>(289, 384)/(193, 192)</td></tr></tbody></table>

# APPENDIX B: CODE VALIDATION AND BENCHMARKS

## B1 Benchmark of the Lorentz force computation

The determination of eq. (21) relies on the computation of $(\mathbf{j} \times \mathbf{B})$ which depends on $\nabla \times \mathbf{B}$ and $\mathbf{B}$. These quantities are computed using SHTns<fnref n="2" /> functions, already validated in (Schaeffer 2013) and widely used in many codes. Given any fields $(B_{pol}, B_{tor})$ in spectral space, and their radial-derivatives, the SHTns routines directly provide $(B_r, B_{\theta}, B_{\phi})$ and $(j_r, j_{\theta}, j_{\phi})$ on the 3D-physical-grid. The remaining part of the process involves the computation of the non-linear products $(\mathbf{j} \times \mathbf{B})$, and the z-averaging of $\partial_s[s(\mathbf{j} \times \mathbf{B})_\phi]$, $(\mathbf{j} \times \mathbf{B})_\phi$ and $(\mathbf{j} \times \mathbf{B})_s$ which are computed on the physical grid. These quantities are then sent to the spectral space using `fft` functions where the $\phi$-derivative is performed before assembling all the terms to obtain eq. (21).

<footnote n="2">²https://bitbucket.org/nschaeff/shtns</footnote>
## B1.1 Summary of the steps involved in the Lorentz force computation

As a summary, the steps involved in our computation of the QG-Lorentz force are presented in the following steps

$$ \text{1) } (B_r, B_\theta, B_\phi) \times (j_r, j_\theta, j_\phi) \rightarrow (\mathbf{j} \times \mathbf{B}) \cdot (\mathbf{e}_r, \mathbf{e}_\theta, \mathbf{e}_\phi); \text{ non-linear products,} $$

$$ \text{2) } (\mathbf{j} \times \mathbf{B})_r, (\mathbf{j} \times \mathbf{B})_\theta \rightarrow (\mathbf{j} \times \mathbf{B})_s; \text{ linear operation,} $$

$$ \text{3) } (\mathbf{j} \times \mathbf{B})_\phi \rightarrow \frac{\partial}{\partial s} [s(\mathbf{j} \times \mathbf{B})_\phi]; \text{ s-derivative,} $$

$$ \text{4) } \left\{ \begin{array}{l} (\mathbf{j} \times \mathbf{B})_\phi \rightarrow \langle (\mathbf{j} \times \mathbf{B})_\phi \rangle \\ (\mathbf{j} \times \mathbf{B})_s \rightarrow \langle (\mathbf{j} \times \mathbf{B})_s \rangle; \text{ z-averages,} \\ \frac{\partial}{\partial s} [s(\mathbf{j} \times \mathbf{B})_\phi] \rightarrow \langle \partial_s [s(\mathbf{j} \times \mathbf{B})_\phi] \rangle \end{array} \right. $$

**physical-to-spectral space**

$$ \text{5) } \left\{ \begin{array}{l} \langle (\mathbf{j} \times \mathbf{B})_\phi \rangle \rightarrow \langle (\widetilde{\mathbf{j} \times \mathbf{B}})_\phi \rangle \\ \langle (\mathbf{j} \times \mathbf{B})_s \rangle \rightarrow \langle (\widetilde{\mathbf{j} \times \mathbf{B}})_s \rangle; \text{ Fourier transforms,} \\ \langle \partial_s [s(\mathbf{j} \times \mathbf{B})_\phi] \rangle \rightarrow \langle \partial_s [s(\widetilde{\mathbf{j} \times \mathbf{B}})_\phi] \rangle \end{array} \right. $$

**spectral space**

$$ \text{6) } \langle (\widetilde{\mathbf{j} \times \mathbf{B}})_s \rangle \rightarrow im(\langle \widetilde{\mathbf{j} \times \mathbf{B}} \rangle_s); \text{ } \phi\text{-derivative,} \\ \text{7) } \left\{ \begin{array}{l} \langle (\widetilde{\mathbf{j} \times \mathbf{B}})_\phi \rangle \rightarrow \frac{1}{Ek Pm} \langle (\widetilde{\mathbf{j} \times \mathbf{B}})_\phi \rangle_{m=0} = \widetilde{F_{L, \bar{u}\bar{\phi}}}; \text{ assemble Lorentz force for the zonal-flow eq.,} \\ \langle (\partial_s [s(\widetilde{\mathbf{j} \times \mathbf{B}})_\phi]), im(\langle \widetilde{\mathbf{j} \times \mathbf{B}} \rangle_s) \rangle \rightarrow \frac{1}{Ek Pm} \left( \frac{1}{s} \langle \partial_s [s(\widetilde{\mathbf{j} \times \mathbf{B}})_\phi] \rangle - \frac{i m}{s} \langle (\widetilde{\mathbf{j} \times \mathbf{B}})_s \rangle \right) = \widetilde{F_{L, \omega_z}} \end{array} \right. \quad \text{(B1)} \\ \text{assemble Lorentz force for the vorticity eq.,} $$

where the wide tilde $\tilde{x}$ refers to the Fourier transform of any quantity $x$, that is the quantity $x$ in the spectral space.

## B1.2 Analytical benchmark from an artificial $\mathbf{j} \times \mathbf{B}$ field

As steps (1) and (2) in eq. (B1) 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 $\phi$, cancels at the outer boundary, but is not trivial. Following these constraints, we thus chose a completely artificial field, that reads

$$ \left\{ \begin{array}{l} (\mathbf{j} \times \mathbf{B})_{\phi}^{\text{ref}}(r, \theta, \phi) = -\left(\frac{1}{2} + \cos \phi + \sin(4\phi)\right) \frac{\pi h^3}{s} \cos\left(\pi \frac{z}{2h}\right) \\ (\mathbf{j} \times \mathbf{B})_{s}^{\text{ref}}(r, \theta, \phi) = \sin^2 \phi s(s_o - s)(s - s_i) z^2 \\ (\mathbf{j} \times \mathbf{B})_{z}^{\text{ref}}(r, \theta, \phi) = \sin^2 \phi s(s_o - s)(s - s_i) z \end{array} \right. \quad \text{(B2)} $$

leading to

$$ F_{L}^{\text{ref}} = \begin{cases} [1/2 + \cos \phi + \sin(4\phi)] 8hs - \frac{2}{3} \cos \phi \sin \phi (s_o - s)(s - s_i) h^2 & . \\ -\frac{h^3}{s} & . \end{cases} \quad \text{(B3)} $$

## B1.3 Relative error definition

In order to discuss the validation and the accuracy of our numerical schemes, we define a relative error estimate, $e_{\text{rel}}$, as

$$ e_{\text{rel}}(f) = \left[ \frac{\{(f_{\text{ref}} - f)^2\}_{\text{QG}}}{\{f_{\text{ref}}^2\}_{\text{QG}}} \right]^{1/2}, \quad \text{(B4)} $$

where the brackets in the above equation correspond to an average over the annulus as in eq. (26).
## B1.1 Summary of the steps involved in the Lorentz force computation

As a summary, the steps involved in our computation of the QG-Lorentz force are presented in the following steps

### physical grid

$$1) (B_r, B_\theta, B_\phi) \times (j_r, j_\theta, j_\phi) \rightarrow (j \times B) \cdot (e_r, e_\theta, e_\phi); \text{ non-linear products,}$$

$$2) (j \times B)_r, (j \times B)_\theta \rightarrow (j \times B)_s; \text{ linear operation,}$$

$$3) (j \times B)_\phi \rightarrow \frac{\partial}{\partial s} [s(j \times B)_\phi]; \text{ s-derivative,}$$

$$4) \left\{ \begin{array}{l} (j \times B)_\phi \rightarrow \langle (j \times B)_\phi \rangle \\ (j \times B)_s \rightarrow \langle (j \times B)_s \rangle; \text{ z-averages,} \\ \frac{\partial}{\partial s} [s(j \times B)_\phi] \rightarrow \langle \partial_s [s(j \times B)_\phi] \rangle \end{array} \right. $$

### physical-to-spectral space

$$5) \left\{ \begin{array}{l} \langle (j \times B)_\phi \rangle \rightarrow \langle \widetilde{(j \times B)_\phi} \rangle \\ \langle (j \times B)_s \rangle \rightarrow \langle \widetilde{(j \times B)_s} \rangle; \text{ Fourier transforms,} \\ \langle \partial_s [s(j \times B)_\phi] \rangle \rightarrow \langle \partial_s [s(\widetilde{j \times B})_\phi] \rangle \end{array} \right. $$

### spectral space

$$7) \left\{ \begin{array}{l} \langle \widetilde{(j \times B)_\phi} \rangle \rightarrow \frac{1}{Ek Pm} \langle \widetilde{(j \times B)_\phi} \rangle_{m=0} = \widetilde{F_{L, \bar{u}\bar{\phi}}}; \text{ assemble Lorentz force for the zonal-flow eq.,} \\ \langle (\partial_s [s(\widetilde{j \times B})_\phi]), i m(\widetilde{(j \times B)_s}) \rangle \rightarrow \frac{1}{Ek Pm} \left( \frac{1}{s} \langle \partial_s [s(\widetilde{j \times B})_\phi] \rangle - \frac{i m}{s} \langle (\widetilde{j \times B})_s \rangle \right) = \widetilde{F_{L, \omega_z}} \quad (B1) \\ \text{assemble Lorentz force for the vorticity eq.,} \end{array} \right. $$

where the wide tilde $\tilde{x}$ refers to the Fourier transform of any quantity $x$, that is the quantity $x$ in the spectral space.

## B1.2 Analytical benchmark from an artificial $j \times B$ field

As steps (1) and (2) in eq. (B1) 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 $\phi$, cancels at the outer boundary, but is not trivial. Following these constraints, we thus chose a completely artificial field, that reads

$$ \left\{ \begin{array}{l} (j \times B)_\phi^{\text{ref}}(r, \theta, \phi) = -\left(\frac{1}{2} + \cos \phi + \sin(4\phi)\right) \frac{\pi h^3}{s} \cos\left(\pi \frac{z}{2h}\right) \\ (j \times B)_s^{\text{ref}}(r, \theta, \phi) = \sin^2 \phi s(s_o - s)(s - s_i) z^2 \\ (j \times B)_z^{\text{ref}}(r, \theta, \phi) = \sin^2 \phi s(s_o - s)(s - s_i) z \end{array} \right. \quad (B2) $$

leading to

$$ F_{L}^{\text{ref}} = \begin{cases} [1/2 + \cos \phi + \sin(4\phi)] 8hs - \frac{2}{3} \cos \phi \sin \phi (s_o - s)(s - s_i) h^2 & . \\ -\frac{h^3}{s} & . \end{cases} \quad (B3) $$

## B1.3 Relative error definition

In order to discuss the validation and the accuracy of our numerical schemes, we define a relative error estimate, $e_{rel}$, as

$$ e_{\text{rel}}(f) = \left[ \frac{\{(f_{\text{ref}} - f)^2\}_{\text{QG}}}{\{f_{\text{ref}}^2\}_{\text{QG}}} \right]^{1/2}, \quad (B4) $$

where the brackets in the above equation correspond to an average over the annulus as in eq. (26).
<figure><img src="image_9.png" /><figcaption>**Figure B1.** Convergence of the relative error of the Lorentz force terms tested against an analytical solution.</figcaption></figure>

## B2 Results

Fig. B1 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 $F_{L, \vec{u}_\phi}$ (green curve)—that only involves a z-averaging of $(\mathbf{j} \times \mathbf{B})_\phi$—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 $F_{L, \omega_z}$ (yellow curve) which involves additional operations (a $\phi$- 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 $(\mathbf{j} \times \mathbf{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, that is 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 towards 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><img src="image_10.png" /><figcaption>**Figure B1.** Convergence of the relative error of the Lorentz force terms tested against an analytical solution.</figcaption></figure>

## B2 Results

Fig. B1 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 $F_{L, \vec{u}_\phi}$ (green curve)—that only involves a z-averaging of $(\mathbf{j} \times \mathbf{B})_\phi$—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 $F_{L, \omega_z}$ (yellow curve) which involves additional operations (a $\phi$- 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 $(\mathbf{j} \times \mathbf{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, that is 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 towards 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.