Self-similarity of the dipole–multipole transition in rapidly rotating dynamos

Abstract

The dipole–multipole transition in rapidly rotating dynamos is investigated through the analysis of forced magnetohydrodynamic waves in an unstably stratified fluid. The focus of this study is on the inertia-free limit applicable to planetary cores, where the Rossby number is small not only on the core depth but also on the length scale of columnar convection. By progressively increasing the buoyant forcing in a linear magnetoconvection model, the slow magnetic–Archimedean–Coriolis (MAC) waves are significantly attenuated so that their kinetic helicity decreases to zero; the fast MAC wave helicity, on the other hand, is practically unaffected. In turn, polarity reversals in low-inertia spherical dynamos are shown to occur when the slow MAC waves disappear under strong forcing. Two dynamically similar regimes are identified – the suppression of slow waves in a strongly forced dynamo and the excitation of slow waves in a moderately forced dynamo starting from a small seed field. While the former regime results in polarity reversals, the latter regime produces the axial dipole from a chaotic multipolar state. For either polarity transition, a local Rayleigh number based on the mean wavenumber of the energy-containing scales bears the same linear relationship with the square of the peak magnetic field measured at the transition. The self-similarity of the dipole–multipole transition can place a constraint on the Rayleigh number for polarity reversals in the Earth.

1. Introduction

The dynamo operating in the Earth's outer core generates a predominantly north–south dipole magnetic field. Occasionally, the magnetic dipole axis flips its orientation and retains its approximate alignment with the Earth's rotation axis. The last such polarity reversal occurred nearly 0.78 million years ago (Merrill 2011). Geomagnetic excursions, the periods during which the magnetic axis wanders up to $45^\circ$ from the rotation axis before returning to its original state, have been more frequent in the Earth's past. As reversals and excursions are likely to result from similar convective states of the core (Gubbins 1999; Valet, Meynadier & Guyodo 2005), it is possible that the geodynamo operates in a strongly driven state below the threshold for reversals.

The first polarity reversals in numerical dynamo models were obtained by Glatzmaier & Roberts (1995a,b). Since then, several other studies reported reversals in comparable parameter regimes (Sarson & Jones 1999; Kutzner & Christensen 2002; Wicht & Olson 2004; Olson, Driscoll & Amit 2009; Sreenivasan, Sahoo & Dhama 2014), and it is now well accepted that dynamo reversals occur under strong buoyancy-driven convection. Any explanation for polarity reversals must follow from an explanation for the preference for the axial dipole in planetary dynamos. The kinetic helicity $\boldsymbol {u}\boldsymbol {\cdot} \boldsymbol {\zeta }$ (where $\boldsymbol {u}$ is the velocity and $\boldsymbol {\zeta }$ is the vorticity) generated in convection columns that arise in rapid rotation is thought to be essential for dipolar dynamo action (Moffatt 1978; Olson, Christensen & Glatzmaier 1999). Consequently, the loss of columnar helicity can lead to collapse of the dipole (e.g. Soderlund, King & Aurnou 2012). A number of numerical dynamo models show a dipole–multipole transition for $Ro_\ell \approx 0.1$, where $Ro_\ell$ is a ‘local’ Rossby number that measures the ratio of nonlinear inertial to Coriolis forces on the length scale of columnar convection (Christensen & Aubert 2006). Since it has been hypothesized that the columnar vortices in rotating turbulence are formed by the propagation of linear inertial waves (Davidson, Staplehurst & Dalziel 2006), the above polarity transition may result from the suppression of inertial waves in the dynamo above a critical value of $Ro_\ell$ (McDermott & Davidson 2019). Multipolar solutions may also be found when the ratio of nonlinear inertia to Lorentz forces exceeds a critical value of $O(1)$ (Tassin, Gastine & Fournier 2021; Zaire et al. 2022). That having been said, in the rapidly rotating limit of zero nonlinear inertia, the helicity of columnar vortices aligned with the rotation axis is considerably enhanced by the magnetic field in the dynamo (Sreenivasan & Jones 2011). It is thought that this field-induced helicity is essential for the formation of the axial dipole in planetary dynamos. Under strong buoyancy-driven convection, the magnetically enhanced columnar flow breaks down, likely causing collapse of the axial dipole (Sreenivasan & Jones 2011). Given the spatial inhomogeneity of the magnetic field (e.g. Schaeffer et al. 2017), an analysis of volume-averaged forces in the dynamo cannot possibly reveal how buoyancy offsets the Lorentz–Coriolis force balance in isolated columns. The analysis of helical dynamo waves at progressively increasing forcing can, on the other hand, provide an insight into the role of buoyancy in dipole collapse. This study distinguishes between polarity-reversing and multipolar dynamos in the inertia-free limit, where the Rossby number is small not only on the depth of the planetary core but also on the columnar length scale transverse to the rotation axis. The analysis of dynamo waves in this limit in turn places a constraint on the convective state of the dynamo that admits polarity reversals.

In convectively driven planetary cores, isolated density disturbances generate fast and slow magnetic–Archimedean–Coriolis (MAC) waves under the combined influence of background rotation, magnetic field and unstable stratification. The fast waves are inertial waves weakly modified by the magnetic field and buoyancy while the slow waves are magnetostrophic waves produced by localized balances between the magnetic, Coriolis and buoyancy forces (Braginsky 1967; Busse et al. 2007, pp. 165–168). The MAC waves in the stably stratified regions of the Earth's core, where gravity acts as a restoring force, have been an active area of research (e.g. Nicolas & Buffett 2023). The present study, on the other hand, focuses on MAC waves in an unstably stratified medium that supports convection in a planetary dynamo. In a rapidly rotating dynamo, the fast MAC waves of frequency $\sim \omega _C$, the frequency of linear inertial waves, exist even in the strong field state where the square of the scaled mean magnetic field, $B^2/2 \varOmega \rho \mu \eta = O(1)$, where $\varOmega$ is the angular velocity of rotation, $\rho$ is the fluid density, $\eta$ is the magnetic diffusivity and $\mu$ is the magnetic permeability. The intensity of slow MAC wave motions becomes comparable to that of the fast waves when the ratio of Alfvén to inertial wave frequencies

(1.1)\begin{equation} \left(\frac{\omega_M}{\omega_C}\right)_{\!0} \sim \frac{V_M}{2 \varOmega \delta} \sim 10^{{-}2}, \end{equation}

where $V_M$ is the Alfvén wave velocity based on the peak magnetic field, $\delta$ is the length scale of the buoyancy disturbance and the subscript ‘0’ refers to the initial state of the disturbance as it is released into the flow (Varma & Sreenivasan 2022). Because of the anisotropy of the convection, the instantaneous value of $\omega _M/\omega _C$ for parity between fast and slow wave motions would be higher than its initial value, and inferred to be ${\sim }0.1$ from the spherical shell dynamo models of Varma & Sreenivasan. In the energy-containing scales, given by the range of spherical harmonic degrees lesser than the mean value at energy injection, the absolute kinetic helicity in the nonlinear dynamo is nearly twice that in the equivalent non-magnetic convection problem. This result implies that the helicity produced by the slow MAC waves in the nonlinear dynamo would be of the same magnitude as that produced by the inertial waves in non-magnetic convection. Further, we infer that the slow waves are essential for dipole formation since the hydrodynamic dynamo at the same parameters, where the Lorentz force is zero, does not generate the axial dipole (see also Sreenivasan & Kar 2018). Although the magnetic diffusion frequency $\omega _\eta$ is the lowest frequency in the dynamo, small but finite magnetic diffusion can place a lower bound on the length scale that supports slow MAC waves in the energy-containing scales of the dynamo. Linear magnetoconvection analysis of a forced damped system (Sreenivasan & Maurya 2021) indicates that, for

(1.2)\begin{equation} \left(\frac{\omega_\eta}{\omega_C}\right)_{\!0} \sim \frac{\eta}{2 \varOmega \delta^2} \lesssim 10^{{-}5}, \end{equation}

the energies of the fast and slow wave motions are comparable. (Here, $\omega _M/\omega _C \sim 0.1$.) Therefore, slow wave motions at length scales $\delta \sim 10\,{\rm km}$ can be influential in the generation of the dipole field. At higher values of the ratio $(\omega _\eta /\omega _C)_0$, the slow wave energy falls below the fast wave energy, implying that scales smaller than $\sim$10 km are rapidly damped by magnetic diffusion.

The role of buoyancy in a rapidly rotating dynamo is not merely in the excitation of MAC waves. In a convective dynamo that evolves from a small seed magnetic field, the slow MAC waves are first excited when $|\omega _M| > |\omega _A|$, where $\omega _A$ is a measure of the strength of buoyancy in an unstably stratified fluid and has the same magnitude as that of the frequency of internal gravity waves in a stably stratified fluid (Varma & Sreenivasan 2022). As the dynamo field intensity increases from its starting seed value, $|\omega _M|$ progressively increases and exceeds $|\omega _A|$, after which the axial dipole eventually forms from a chaotic multipolar state. Thus, the inequality $|\omega _C| > |\omega _M| >|\omega _A|> |\omega _\eta |$ represents a large region of the parameter space where dipole-dominated dynamos exist. If a rotating dynamo is subject to progressively increasing buoyant forcing, the larger self-generated fields would result in progressively larger $\omega _M$ until a state is reached where $|\omega _A| \sim |\omega _M|$ as the field attains its highest intensity for a given rotation rate. Here, the slow MAC waves disappear, likely causing collapse of the dipole. This transition from dipole to polarity-reversing states is dynamically similar to the transition from a multipolar state to the dipole that occurs in the growth phase of a dynamo starting from a seed field. Since $\omega _C$ remains the dominant frequency while forcing is increased, the dynamo reverses polarity in a rotationally dominant regime where slow magnetostrophic waves are suppressed. Although increased forcing may result in enhanced nonlinear inertia in numerical dynamo models, we consider it unlikely that inertia has any role in polarity transitions in the Earth even for the smallest buoyancy disturbances that support MAC waves. For disturbances of size $L_\perp \approx 15$ km transverse to the rotation axis, the actual ratio of nonlinear inertia to Coriolis forces

(1.3)\begin{equation} \dfrac{|\boldsymbol{\nabla} \times (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{u}|} {|2 (\boldsymbol{\varOmega} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{u}|} \sim \dfrac{u_\star L}{2 \varOmega \, L_\perp^2} \approx 0.03, \end{equation}

taking $L= 2260$ km and $u_\star = 5 \times 10^{-4}$ m s$^{-1}$ (Starchenko & Jones 2002). Low-inertia numerical dynamos must have a small Rossby number based on the length scale of convection. The use of magnetic Prandtl number $Pm =\nu /\eta \sim 1\unicode{x2013}10$, where $\nu$ is the kinematic viscosity, is useful in realizing strong magnetic fields in numerical simulations (Sreenivasan & Jones 2006; Willis, Sreenivasan & Gubbins 2007; Teed, Jones & Tobias 2015; Dormy 2016) in the inertialess regime relevant to rotating planetary cores. Although the choice of a large $Pm$ at moderate Ekman number $E=\nu /2\varOmega L^2 \sim 10^{-4}$–$10^{-6}$ has the unphysical consequence of the viscous dissipation being at least as high as the Ohmic dissipation, the advantage derived in terms of reproducing the MAC force balance in the energy-containing scales is crucial in the understanding of wave motions in both dipole-dominated and reversing regimes. The present study analyses polarity transitions in strongly driven, low-inertia dynamos.

In § 2, we consider the evolution of a buoyancy disturbance in an unstably stratified rotating fluid subject to a magnetic field. In a Cartesian geometry, the axes parallel to gravity, rotation and the magnetic field are chosen to be mutually orthogonal, which is referred to as the ‘equatorial toroidal configuration’ by Loper, Chulliat & Shimizu (2003). At times much shorter than the time scale for exponential increase of the perturbation, the relative intensities of the fast and slow MAC wave motions generated by the perturbation are studied. Apart from the dipole-dominated regime given by $|\omega _C|> |\omega _M| > |\omega _A| > |\omega _\eta |$, the regime thought to be relevant to polarity transitions, $|\omega _C|> |\omega _M| \approx |\omega _A| > |\omega _\eta |$, is analysed. Since the relative magnitudes of the frequencies in these inequalities do not depend on the precise orientation of the gravity, rotation and magnetic field axes, the Cartesian linear model serves as the basis for the study of the role of wave motions in polarity transitions in nonlinear dynamo models, given in § 3. Here, we find that the formation of the axial dipole from a chaotic multipolar state and the collapse of the axial dipole into a polarity-reversing state are dynamically similar phenomena, in that they both occur at $|\omega _A/\omega _M| \sim 1$. While dipole formation requires the excitation of slow MAC waves as a dynamo evolves from a small seed magnetic field, polarity reversals occur when the slow waves are suppressed in a strongly driven dynamo. The self-similarity of the dipole–multipole transition in the inertia-free regime places a constraint on the Rayleigh number for reversals. In § 4, we discuss the implications of our results for planetary cores and future work.

2. Evolution of a density disturbance under rapid rotation and a magnetic field

2.1. Problem set-up and governing equations

A localized density disturbance $\rho ^\prime$ that occurs in an unstably stratified rotating fluid layer threaded by a uniform magnetic field is considered. Since $\rho ^\prime$ is related to a temperature perturbation $\varTheta$ by $\rho ^\prime = -\rho \alpha \varTheta$, where $\rho$ is the ambient density and $\alpha$ is the coefficient of thermal expansion, an initial temperature perturbation is chosen in the form

(2.1)\begin{equation} \varTheta_0=C \exp [-(x^2+y^2+z^2)/\delta^2], \end{equation}

where $C$ is a constant and $\delta$ is the length scale of the perturbation. Figure 1 shows the initial perturbation which subsequently evolves under gravity $\boldsymbol {g} = -g \hat {\boldsymbol {e}}_y$, background rotation $\boldsymbol {\varOmega } = \varOmega \hat {\boldsymbol {e}}_z$ and a uniform magnetic field $\boldsymbol {B} = B \hat {\boldsymbol {e}}_x$ in Cartesian coordinates $(x,y,z)$. In an otherwise quiescent medium, the initial temperature perturbation (2.1) gives rise to a velocity field $\boldsymbol {u}$, which in turn interacts with $\boldsymbol {B}$ to generate the induced magnetic field $\boldsymbol {b}$. The initial velocity $\boldsymbol {u}_0$ and induced field $\boldsymbol {b}_0$ are both zero. In the Boussinesq approximation, the following magnetohydrodynamic (MHD) equations give the evolution of $\boldsymbol {u}$, $\boldsymbol {b}$ and $\varTheta$:

(2.2)$$\begin{gather} \frac{\partial{\boldsymbol{u}}}{\partial{t}} ={-}\frac{1}{\rho}\boldsymbol{\nabla}{p^*} -2 \boldsymbol{\varOmega} \times \boldsymbol{u} + \frac{1}{\mu \rho}\left(\boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} \right) \boldsymbol{b} - \boldsymbol{g} \alpha\varTheta +\nu\nabla^2\boldsymbol{u}, \end{gather}$$

(2.3)$$\begin{gather}\frac{\partial{\boldsymbol{b}}}{\partial{t}} =\left(\boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} \right) {\boldsymbol{u}} +\eta \nabla^2 \boldsymbol{b}, \end{gather}$$

(2.4)$$\begin{gather}\frac{\partial{\varTheta}}{\partial{t}}={-}\gamma \boldsymbol{\hat{e}}_y \boldsymbol{\cdot} \boldsymbol{u}+ \kappa\nabla^2\varTheta, \end{gather}$$

(2.5)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{b}=0, \end{gather}$$

where $\nu$ is the kinematic viscosity, $\kappa$ is the thermal diffusivity, $\eta$ is the magnetic diffusivity, $\mu$ is the magnetic permeability, $\boldsymbol {\varOmega }= \varOmega \hat {\boldsymbol {e}}_z$, $p^*=p -(\rho /2) |{\boldsymbol {\varOmega }} \times {\boldsymbol {x}}|^2+ (\boldsymbol {B} \boldsymbol {\cdot} \boldsymbol {b})/\mu$ and $\gamma =\partial T_0/\partial y <0$ is the mean temperature gradient in the unstably stratified fluid.

Figure 1. The initial state of a density perturbation $\rho ^{\prime }$ that evolves in an unstably stratified fluid subject to a uniform magnetic field $\boldsymbol {B} = B \hat {\boldsymbol {e}}_x$, background rotation $\boldsymbol {\varOmega } = \varOmega \hat {\boldsymbol {e}}_z$ and gravity $\boldsymbol {g} = -g \hat {\boldsymbol {e}}_y$ in Cartesian coordinates $(x,y,z)$.

2.2. Solutions for the velocity field

Taking the curl of (2.2) and (2.3) and eliminating the electric current density between them, we obtain for the velocity field

(2.6)\begin{align} &\left[\left(\frac{\partial}{\partial{t}} -\nu \nabla^2\right)\left(\frac{\partial}{\partial{t}}-\eta \nabla^2\right) -V_M^2 \frac{\partial^2{}}{\partial{x^2}}\right]^2 (-\nabla^2 \boldsymbol{u})\nonumber\\ &\quad =4\varOmega^2\left(\frac{\partial}{\partial{t}} -\eta\nabla^2\right)^2\frac{\partial^2{\boldsymbol{u}}}{\partial{z^2}} -2\varOmega\alpha\frac{\partial}{\partial{z}}\left(\frac{\partial}{\partial{t}} -\eta\nabla^2\right)^2(\boldsymbol{\nabla}\times\varTheta\boldsymbol{g})\nonumber\\ &\qquad -\alpha\left(\frac{\partial}{\partial{t}} -\eta \nabla^2\right)\left[\left(\frac{\partial}{\partial{t}} -\nu \nabla^2\right)\left(\frac{\partial}{\partial{t}}-\eta \nabla^2\right) - V_M^2 \frac{\partial^2{}}{\partial{x^2}}\right] (\boldsymbol{\nabla}\times\boldsymbol{\nabla}\times\varTheta\boldsymbol{g}), \end{align}

where $V_M= B/\sqrt {\mu \rho }$ is the Alfvén velocity. The assumption of an unbounded domain facilitates the use of Fourier transforms along each coordinate direction,

(2.7) \begin{equation} \mathscr{F}(\boldsymbol{A}) = \hat{\!\!\boldsymbol{A}}= \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\boldsymbol{A} {{\rm e}}^{- \mathrm{i} \, \boldsymbol{k}\ \boldsymbol{\cdot}\ \boldsymbol{x}}\,\mathrm{d}\kern0.7pt x\, \mathrm{d}y\,\mathrm{d}z, \end{equation}

and

(2.8) \begin{equation} \mathscr{F}^{{-}1}\left(\,\,\hat{\!\!\boldsymbol{A}}\right)=\boldsymbol{A}= {\frac{1}{(2{\rm \pi})^3}} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\hat{\!\!\boldsymbol{A}} \mathrm{e}^{\mathrm{i} \,\boldsymbol{k}\ \boldsymbol{\cdot}\ \boldsymbol{x}}\, \mathrm{d}k_x\, \mathrm{d}k_y\,\mathrm{d}k_z, \end{equation}

where $\boldsymbol {k}=(k_x,k_y,k_z)$ represents the wave vector such that $|\boldsymbol {k}|=k=\sqrt {k_x^2+k_y^2+k_z^2}$. Application of the Fourier transform (2.7) to the Cartesian components of (2.6) gives

(2.9)\begin{align} &\left[\left(\left(\frac{\partial}{\partial{t}} +\omega_\nu\right)\left(\frac{\partial}{\partial{t}} +\omega_\eta\right)+\omega_M^2\right)^2\left(\frac{\partial}{\partial{t}} +\omega_\kappa\right)+\omega_C^2\left(\frac{\partial}{\partial{t}} +\omega_\eta\right)^2\left(\frac{\partial}{\partial{t}} +\omega_\kappa\right)\right]\hat{u}_x\nonumber\\ &\quad=\left[-\omega_C\omega_A^2\frac{k_zk}{k_x^2+k_z^2} \left(\frac{\partial}{\partial{t}}+\omega_\eta\right)^2\right.\nonumber\\ &\left.\qquad+\,\omega_A^2 \frac{k_xk_y}{k_x^2+k_z^2}\left(\frac{\partial}{\partial{t}} +\omega_\eta\right) \left(\left(\frac{\partial}{\partial{t}} +\omega_\nu\right)\left(\frac{\partial}{\partial{t}} +\omega_\eta\right)+\omega_M^2\right)\right]\hat{u}_y, \end{align}

(2.10)\begin{align} &\left[\left(\left(\frac{\partial}{\partial{t}} +\omega_\nu\right)\left(\frac{\partial}{\partial{t}} +\omega_\eta \right)+\omega_M^2\right)^2\left(\frac{\partial}{\partial{t}} +\omega_\kappa\right)+\omega_C^2\left(\frac{\partial}{\partial{t}} +\omega_\eta\right)^2\left(\frac{\partial}{\partial{t}} +\omega_\kappa\right)\right.\nonumber\\ &\left.\qquad +\,\omega_A^2\left(\frac{\partial}{\partial{t}} +\omega_\eta\right) \left(\left(\frac{\partial}{\partial{t}} +\omega_\nu\right)\left(\frac{\partial}{\partial{t}} +\omega_\eta\right)+\omega_M^2\right)\right]\hat{u}_y=0, \end{align}

(2.11)\begin{align} &\left[\left(\left(\frac{\partial}{\partial{t}} +\omega_\nu\right)\left(\frac{\partial}{\partial{t}} +\omega_\eta\right)+\omega_M^2\right)^2\left(\frac{\partial}{\partial{t}} +\omega_\kappa\right)+\omega_C^2\left(\frac{\partial}{\partial{t}} +\omega_\eta\right)^2\left(\frac{\partial}{\partial{t}} +\omega_\kappa\right)\right]\hat{u}_z\nonumber\\ &\quad=\left[\omega_C\omega_A^2\frac{k_x k}{k_x^2+k_z^2} \left(\frac{\partial}{\partial{t}} +\omega_\eta\right)^2\right.\nonumber\\ &\left.\qquad+\,\omega_A^2\frac{k_yk_z}{k_z^2+k_x^2}\left(\frac{\partial}{\partial{t}} +\omega_\eta\right) \left(\left(\frac{\partial}{\partial{t}} +\omega_\nu\right)\left(\frac{\partial}{\partial{t}} +\omega_\eta\right)+\omega_M^2\right)\right]\hat{u}_y, \end{align}

where we have combined the transformed temperature equation (2.4) with the transform of (2.6). In (2.9)–(2.11),

(2.12a–c)\begin{equation} \omega_C^2=4 \varOmega^2 k_z^2/k^2, \quad \omega_A^2=g\alpha \gamma (k_x^2+k_z^2)/k^2 , \quad \omega_M^2=V_M^2 k_x^2, \end{equation}

represent the squares of the frequencies of linear inertial, buoyancy and Alfvén waves, respectively (Braginsky 1967; Busse et al. 2007). In an unstably stratified medium, wherein $\omega _A^2 <0$, $|\omega _A|$ is a measure of the strength of buoyancy. As the present study focuses on a system where the viscous and thermal diffusion are much smaller than magnetic diffusion, the frequencies $\omega _\nu = \nu k^2$ and $\omega _\kappa = \kappa k^2$ in (2.9)–(2.11) are small compared with $\omega _\eta = \eta k^2$. In this limit, the solution of the form $\hat {u}_y \sim {{\rm e}}^{\mathrm {i} \lambda t}$ for the homogeneous equation (2.10) gives the following quintic equation in $\lambda$:

(2.13)$$\begin{gather} \lambda^5- 2 \mathrm{i}\omega_\eta \lambda^4 - (\omega_A^2+\omega_C^2+\omega_\eta^2+2\omega_M^2) \lambda^3 +2 \mathrm{i} \omega_\eta( \omega_A^2+\omega_C^2+\omega_M^2) \lambda^2\nonumber\\ +\,(\omega_A^2\omega_\eta^2+\omega_C^2\omega_\eta^2+ \omega_A^2\omega_M^2+\omega_M^4) \lambda- \mathrm{i} \omega_A^2\omega_\eta\omega_M^2=0, \end{gather}$$

the approximate roots of which were discussed in Sreenivasan & Maurya (2021) for known relative orders of magnitudes of $\omega _M$, $\omega _A$, $\omega _C$ and $\omega _\eta$. In the solution for (2.10)

(2.14)\begin{equation} \hat{u}_y = \sum_{m=1}^{5}D_{m} \mbox{e}^{\mathrm{i} \lambda_{m}t}, \end{equation}

the coefficients $D_{m}$ are determined using the initial conditions for $\hat {u}_y$ and its time derivatives (§ 2.3). Of the five terms in the expansion on the right-hand side of (2.14), two terms represent oppositely travelling fast MAC waves, two other terms represent oppositely travelling slow MAC waves and the fifth term represents the overall growth of the velocity perturbation.

By substituting (2.14) in (2.9) and (2.11), the following solutions for (2.9) and (2.11) are obtained:

(2.15)\begin{gather} \hat{u}_{x} =\hat{u}^H_{x}+\hat{u}^P_{x}=\sum_{m=1}^{5}A_{m} \mathrm{e}^{\mathrm{i}\lambda_m^H t}+\sum_{m=1}^{5}M_m{\rm e}^{\mathrm{i} \lambda_m^P t}, \end{gather}

(2.16)\begin{gather} \hat{u}_{z} =\hat{u}^H_{z}+\hat{u}^P_{z}=\sum_{m=1}^{5}C_{m} \mathrm{e}^{\mathrm{i}\lambda_m^H t} +\sum_{m=1}^{5}N_m{\rm e}^{\mathrm{i}\lambda_m^P t} .\end{gather}

In (2.15) and (2.16), $\hat {u}_x^H$ and $\hat {u}_z^H$ are the homogeneous solutions of (2.9) and (2.11), respectively, $\hat {u}_x^P$ and $\hat {u}_z^P$ are the particular solutions, $\lambda _m^P$ are the roots of (2.13) and $\lambda _m^H$ are the roots of (2.13) with $\omega _A=0$

(2.17)\begin{equation} \left. \begin{aligned} \lambda^H_1 & =\frac{1}{2}\left(\omega_C+ \mathrm{i} \omega_\eta+\sqrt{\omega_C^2-2 \mathrm{i} \omega_C\omega_\eta-\omega_\eta^2+4\omega_M^2}\right),\\ \lambda^H_2 & = \frac{1}{2} \left(- \omega_C+ \mathrm{i} \omega_\eta - \sqrt{\omega_C^2 + 2 \mathrm{i} \omega_C\omega_\eta-\omega_\eta^2+4\omega_M^2}\right),\\ \lambda^H_3 & =\frac{1}{2}\left(\omega_C+ \mathrm{i} \omega_\eta-\sqrt{\omega_C^2-2 \mathrm{i}\omega_C\omega_\eta-\omega_\eta^2+4\omega_M^2}\right),\\ \lambda^H_4 & =\frac{1}{2} \left(- \omega_C+ \mathrm{i} \omega_\eta + \sqrt{\omega_C^2 +2 \mathrm{i}\omega_C\omega_\eta -\omega_\eta^2+4\omega_M^2}\right),\\ \lambda^H_5 & =0, \end{aligned} \right\} \end{equation}

which give the frequencies of the fast and slow magneto–Coriolis (MC) waves in the absence of buoyancy (Sreenivasan & Narasimhan 2017). The coefficients $A_m$, $C_m$, $M_m$ and $N_m$ in (2.15) and (2.16) are evaluated as in § 2.3 below.

The solutions for the induced magnetic field transforms $\hat {b}_x$, $\hat {b}_y$ and $\hat {b}_z$ are obtained following a similar approach.

2.3. Evaluation of spectral coefficients

From (2.14), the initial conditions for $\hat {u}_y$ and its time derivatives are given by

(2.18)\begin{equation} \mathrm{i}^n \sum_{m=1}^{5}D_m \lambda_{m}^n =\left(\frac{\partial^n\hat{u}_y}{\partial t^n}\right)_{t=0}=a_{n+1}, \quad n=0,1,2,3,4. \end{equation}

Algebraic simplifications give the right-hand sides of (2.18) in the limit of $\nu =\kappa =0$, as follows:

(2.19)\begin{equation} \left. \begin{aligned} a_1 & =\hat{u}_y|_{t=0}=0,\\ a_2 & =\frac{\partial{\hat{u}_y}}{\partial{t}}|_{t=0}= \alpha g \left(\frac{k_z^2+k_x^2}{k^2}\right)\hat{\varTheta}_0,\\ a_3 & =\frac{\partial^2 \hat{u}_y}{\partial{t^2}}|_{t=0}=0,\\ a_4 & =\frac{\partial^3{\hat{u}_y}}{\partial{t^3}}|_{t=0}={-}(\omega_M^2+\omega_C^2+\omega_A^2) a_2,\\ a_5 & =\frac{\partial^4{\hat{u}_y}}{\partial{t^4}}|_{t=0}= \omega_M^2\omega_\eta a_2. \end{aligned} \right\} \end{equation}

The coefficients $D_m$ may now be obtained using the roots of (2.13). For example, we obtain

(2.20)$$\begin{gather} D_{1}=\dfrac{a_{5}- \mathrm{i} \, a_{4}(\lambda_2+\lambda_3+\lambda_4+\lambda_5) +\mathrm{i} \, a_{2}(\lambda_2\lambda_4\lambda_5+\lambda_3\lambda_4\lambda_5 +\lambda_2\lambda_3\lambda_4+\lambda_2\lambda_3\lambda_5)}{(\lambda_1-\lambda_2) (\lambda_1-\lambda_3)(\lambda_1-\lambda_4)(\lambda_1-\lambda_5)}, \end{gather}$$

(2.21)$$\begin{gather}D_{3}=\dfrac{a_{5}- \mathrm{i} \, a_{4}(\lambda_1+\lambda_2+\lambda_4+\lambda_5) + \mathrm{i} \,a_{2}(\lambda_1\lambda_4\lambda_5+\lambda_2\lambda_4\lambda_5 +\lambda_1\lambda_2\lambda_4+\lambda_1\lambda_2\lambda_5)}{(\lambda_3-\lambda_1) (\lambda_3-\lambda_2)(\lambda_3-\lambda_4)(\lambda_3-\lambda_5)}, \end{gather}$$

for the forward-travelling fast and slow wave solutions, respectively. The coefficients $A_m$ and $C_m$ are determined in a similar way. The coefficients $M_m$ and $N_m$ in equations (2.15) and (2.16) are obtained using the method of undetermined coefficients, as follows:

(2.22)\begin{align} M_m &= \frac{D_{m}}{T_m}\left[-\omega_C\omega_A^2 \left(\frac{kk_z}{k_x^2+k_z^2}\right) \big(-(\lambda^P_m)^2+\omega_\eta^2 +2 \mathrm{i} \omega_\eta\lambda^P_m\big)\right.\nonumber\\ &\quad \left.+\,\omega_A^2\frac{k_xk_y}{k_x^2+k_z^2} \big(- \mathrm{i}(\lambda^P_m)^3-2\omega_\eta(\lambda^P_m)^2 + \mathrm{i}(\omega_\eta^2+\omega_M^2)\lambda^P_m +\omega_M^2\omega_\eta\big)\right], \end{align}

and

(2.23)\begin{align} N_m& = \frac{D_{m}}{T_m}\left[\omega_C\omega_A^2 \left(\frac{kk_x}{k_x^2+k_z^2}\right) \big(-(\lambda^P_m)^2+\omega_\eta^2 +2 \mathrm{i}\omega_\eta\lambda^P_m\big)\right.\nonumber\\ &\left.\quad+\,\omega_A^2\frac{k_zk_y}{k_x^2+k_z^2} \big(- \mathrm{i}(\lambda^P)^3-2\omega_\eta(\lambda^P_m)^2 + \mathrm{i}(\omega_\eta^2+\omega_M^2)\lambda^P_m +\omega_M^2\omega_\eta\big)\right], \end{align}

with

(2.24)$$\begin{gather} T_m = \mathrm{i} (\lambda^P_m)^5+2(\lambda^P_m)^4\omega_\eta - \mathrm{i}(\lambda^P_m)^3(\omega_C^2+\omega_\eta^2 +2\omega_M^2)-2\omega_\eta(\lambda^P_m)^2(2\omega_C^2+\omega_M^2)\nonumber\\ +\, \mathrm{i}\lambda^P_m(\omega_C^2\omega_\eta^2+\omega_M^4), \quad m=1,2,\ldots,5. \end{gather}$$

2.4. Fast and slow MAC waves in unstable stratification

From the solution of the initial value problem, the velocity and induced magnetic field can be obtained at discrete points in time. The analysis of the solutions is limited to times much shorter than the time scale for the exponential increase of the perturbations. When the buoyancy force is small compared with the Lorentz force ($|\omega _A/\omega _M| \ll 1$), the parameter regime is determined by the Lehnert number $Le$ and the magnetic Ekman number $E_\eta$

(2.25a,b)\begin{equation} Le= \dfrac{V_M}{2 \varOmega \delta}, \quad E_\eta= \dfrac{\eta}{2 \varOmega \delta^2}, \end{equation}

both based on the length scale $\delta$ of the initial buoyancy perturbation (2.1).

Figure 2 shows the evolution of the kinetic helicity $\boldsymbol {u} \boldsymbol {\cdot} \boldsymbol {\zeta }$ at $y=0$. The real-space fields are obtained from the transforms $\hat {\boldsymbol {u}}$ and $\hat {\boldsymbol {\zeta }}$ via the inverse Fourier transform (2.8). Here, a truncation value of $\pm 10/\delta$ is used for the three wavenumbers in the integrals since the initial wavenumber $k_0= \sqrt {3}/\delta$ (Appendix A). Apart from the segregation of oppositely signed helicity between the two halves about the mid-plane $z=0$, the evolution of blobs into columnar structures through the propagation of damped waves is evident.

Figure 2. Evolution of the kinetic helicity on the $x$–$z$ plane at $y=0$ with time (measured in units of the magnetic diffusion time $t_\eta$) for $Le=0.03$ and $E_\eta =2\times 10^{-5}$. The snapshots are at (a) $t/t_\eta =1\times 10^{-4}$, (b) $t/t_\eta =2.5\times 10^{-3}$ and (c) $t/t_\eta = 1\times 10^{-2}$. The ratio $|\omega _A/\omega _M| =0.05$ at times after the formation of the waves.

As we seek to understand the role of MAC waves in the dipolar and multipolar dynamo regimes, we may separate the fast and slow MAC wave parts of the general solution, which is a linear superposition of the two wave solutions (see Sreenivasan & Maurya 2021). For example

(2.26)\begin{equation} \left. \begin{aligned} \hat{u}_{x,f} & =M_1 \mathrm{e}^{\mathrm{i}\lambda^P_1 t}+ M_2 \mathrm{e}^{\mathrm{i}\lambda^P_2 t} + A_{1} \mathrm{e}^{\mathrm{i}{\lambda}_1^H t}+ A_{2} \mathrm{e}^{\mathrm{i} {\lambda}_2^H t},\\ \hat{u}_{y,f} & = D_1 \mathrm{e}^{\mathrm{i}\lambda_1 t}+ D_2 \mathrm{e}^{\mathrm{i}\lambda_2 t},\\ \hat{u}_{z,f} & =N_1 \mathrm{e}^{\mathrm{i} \lambda^P_1 t}+ N_2 \mathrm{e}^{\mathrm{i} \lambda^P_2 t} +C_{1} \mathrm{e}^{\mathrm{i} {\lambda}_1^H t}+C_{2} \mathrm{e}^{\mathrm{i}{\lambda}_2^H t}, \end{aligned} \right\} \end{equation}

and

(2.27)\begin{equation} \left. \begin{aligned} \hat{u}_{x,s} & = M_3 \mathrm{e}^{\mathrm{i} \lambda^P_3 t} + M_4 \mathrm{e}^{\mathrm{i}\lambda^P_4 t} +A_{3} \mathrm{e}^{\mathrm{i}{\lambda}_3^H t} + A_{4} \mathrm{e}^{\mathrm{i}{\lambda}_4^H t},\\ \hat{u}_{y,s} & = D_3 \mathrm{e}^{\mathrm{i}\lambda_3 t}+D_4 \mathrm{e}^{\mathrm{i}\lambda_4 t},\\ \hat{u}_{z,s} & =N_3 \mathrm{e}^{\mathrm{i}\lambda^P_3 t}+ N_4 \mathrm{e}^{\mathrm{i}\lambda^P_4 t}+C_{3} \mathrm{e}^{\mathrm{i}{\lambda}_3^H t}+C_{4} \mathrm{e}^{\mathrm{i} {\lambda}_4^H t}, \end{aligned} \right\} \end{equation}

where the subscripts $f$ and $s$ on the left-hand sides of (2.26) and (2.27) denote the fast and slow wave parts of the solution. Figure 3(a) shows the variation of fundamental frequencies with $Le$. To compute the frequencies, the mean wavenumbers are first calculated through ratios of $L^2$ norms; e.g.

(2.28a,b)\begin{equation} \bar{k}_x = \frac{\|k_x \hat{u}\| }{\| \hat{u}\|}, \quad \bar{k} = \frac{\| k \hat{u}\|}{\| \hat{u}\|}, \end{equation}

which are based on the velocity field. The values of $|\omega _M/\omega _C|$, given in brackets below the horizontal axis of figure 3(b), are systematically higher than that of $Le$, which is essentially the initial value of this ratio. The enhanced instantaneous value of $|\omega _M/\omega _C|$ is due to the anisotropy of the columnar flow, and would not be evident if $\boldsymbol {\varOmega }$ were aligned with $\boldsymbol {B}$ (Sreenivasan & Maurya 2021). For $E_\eta = 2 \times 10^{-5}$, all calculations for $Le > 2 \times 10^{-3}$ satisfy the inequality $|\omega _C| > |\omega _M| > |\omega _A|>|\omega _\eta |$, thought to be essential for axial dipole formation in convective dynamos. Figure 3(b) shows the variation of the dimensionless helicity $h^*$ of the fast and slow MAC waves, obtained by summing the helicity at all points in $(x,z)$ for $z<0$ and $y=0$ and then normalizing this value by the non-magnetic helicity. For $Le > 2 \times 10^{-3}$ ($|\omega _M/\omega _C| > 0.045$), the slow wave helicity increases dramatically, and for $Le \sim 10^{-2}$ ($|\omega _M/\omega _C| \sim 0.1$), the slow wave helicity is greater than the fast wave helicity. This result is consistent with the higher intensity of the slow wave motions relative to that of the fast wave motions for $|\omega _M/\omega _C| \sim 0.1$, inferred from the fast Fourier transform (FFT) spectra in dynamo simulations (Varma & Sreenivasan 2022).

Figure 3. (a) Variation of absolute values of frequencies with $Le$. (b) Variation with $Le$ of the total kinetic helicity in the region $z<0$ normalized by the non-magnetic helicity. The values of $|\omega _M/\omega _C|$ that correspond to the values of $Le$ are given in brackets below the horizontal axis. All calculations are performed for $E_\eta = 2 \times 10^{-5}$. The slow wave frequency, $\omega _s$, takes non-zero values for $Le > 2 \times 10^{-3}$, when $|\omega _M| > |\omega _A|$.

In figure 4, the contours of the fast and slow MAC wave helicities are shown at two times for $|\omega _M/\omega _C|=0.13$, which lies in the region of slow wave dominance in figure 3(b). The fast waves split in two and propagate rapidly along $z$. The slow waves do not propagate as far as the fast waves at the same time due to their lower group velocity. Yet, as indicated by the colour bars, the slow waves are markedly more intense than the fast waves. Both fast and slow wave columns propagate along $x$ at the Alfvén velocity (§ 2.5 below).

Figure 4. Helicity on the section $y=0$ of fast (a,b) and slow (c,d) MAC waves at two times, $t/t_\eta = 2.5\times 10^{-3}$ (a,c) and $t/t_\eta =1 \times 10^{-2}$ (b,d). Here, $E_\eta = 2 \times 10^{-5}$ and $Le=0.03$ ($|\omega _M/\omega _C| =0.13$).

The effect of progressively increasing the buoyant forcing on the fast and slow waves is examined in the following section.

2.5. Effect of progressively increasing buoyancy

For the regime given by $|\omega _C| \gg |\omega _M| \gg |\omega _A| \gg |\omega _\eta |$, the roots of the homogeneous equation (2.13) are approximated by (Sreenivasan & Maurya 2021)

(2.29)$$\begin{gather} \lambda_{1,2} \approx{\pm} \left(\omega_C+ \frac{\omega_M^2}{\omega_C} \right) + \mathrm{i} \, \frac{\omega_M^2\omega_\eta}{\omega_C^2}, \end{gather}$$

(2.30)$$\begin{gather}\lambda_{3,4} \approx{\pm} \left(\frac{\omega_M^2}{\omega_C}+ \frac{\omega_A^2}{2\omega_C}\right) + \mathrm{i} \, \omega_\eta \, \left(1-\frac{\omega_A^2}{2\omega_M^2}\right), \end{gather}$$

(2.31)$$\begin{gather}\lambda_{5} \approx \mathrm{i} \, \frac{\omega_A^2\omega_\eta}{\omega_M^2}. \end{gather}$$

As $|\omega _A|$ increases relative to $|\omega _M|$, the fast MAC waves, given by frequencies $\lambda _{1,2}$, are nearly unaffected. However, the slow waves, whose real frequencies are approximated by

(2.32)\begin{equation} \mbox{Re}(\lambda_{3,4}) \approx{\pm} \left(\frac{\omega_M^2}{\omega_C}+ \frac{\omega_A^2}{2\omega_C}\right) \approx{\pm} \frac{\omega_M^2}{\omega_C} \, \left(1+ \frac{\omega_A^2}{\omega_M^2} \right)^{1/2}, \end{equation}

for $|\omega _C| \gg |\omega _M|, |\omega _A|$ (Braginsky 1967), would be significantly attenuated in an unstably stratified fluid with $\omega _{A}^2<0$ as $|\omega _A|$ nears $|\omega _M|$. We see below that the decrease of the slow wave frequency translates into the marked decrease of the slow wave helicity relative to the fast wave helicity.

Figure 5(a) indicates that both fast and slow MAC waves propagate along the mean-field direction $x$ such that $x/\delta = t/t_a$, where $t_a$ is the Alfvén wave travel time. Contrary to that in the inertial–Alfvén wave system discussed in Bardsley & Davidson (2017), the Alfvén waves propagating in the direction orthogonal to the rotation axis are simply the degenerate form of the MAC waves (see also Varma & Sreenivasan 2022). For small $|\omega _A/\omega _M|$, the helicity of slow wave motions is greater than that of fast waves, but as $|\omega _A/\omega _M|$ approaches unity, the slow wave helicity weakens considerably and falls below that of the fast wave. The effect of increasing buoyancy forcing on the fast and slow wave helicity is shown graphically in figure 6. The fast waves are practically unaffected by the strength of forcing as their intensity and $z$ propagation rate are nearly invariant for $|\omega _A/\omega _M|$ in the range 0.1–1 (figure 6a–c). The slow wave helicity, on the other hand, is substantially weakened as $|\omega _A/\omega _M|$ increases in the same range (figure 6d–f). For $|\omega _A/\omega _M| \approx 1$, a state is reached where the slow wave helicity is nearly zero. The induced magnetic field $b_z$ also weakens considerably with increasing $|\omega _A/\omega _M|$ (figure 6g–i), which indicates that only the slow waves have a direct bearing on field generation.

Figure 5. (a) Variation of the fast wave (dashed line) and slow wave (solid line) helicity for $z<0$, normalized by the non-magnetic helicity, with $x/\delta$ and $|\omega _A/\omega _M|$ at different times (measured in units of the Alfvén wave travel time $t_a$). Here, $E_\eta =2\times 10^{-5}$ and $Le=0.03$. (b) Variation of the normalized slow wave helicity with the local Rayleigh number $Ra_\ell$ (defined in (2.33)) and $|\omega _A/\omega _M|$ for different $Le$.

Figure 6. Fast MAC wave helicity (a–c), slow MAC wave helicity (d–f) and $z$-component of the induced magnetic field, $b_z$ (g–i) for $|\omega _A/\omega _M| = 0.1$ (a,d,g), 0.6 (b,e,h), 0.95 (c,f,i). The plots are generated at time $t/t_\eta =0.01$ for the parameters $Le=0.03$ and $E_\eta =2 \times 10^{-5}$.

Figure 5(b) shows the normalized slow wave helicity against a ‘local’ Rayleigh number based on the length scale of the initial perturbation

(2.33)\begin{equation} Ra_\ell = \dfrac{g \alpha |\gamma| \delta^2}{2 \varOmega \eta}, \end{equation}

as well as $|\omega _A/\omega _M|$. Evidently, the forcing needed to suppress the slow waves increases with $Le$, although the total suppression of these waves occurs universally at $|\omega _A/\omega _M| \approx 1$. This result prompts us to look at the condition for vanishing slow wave helicity through a relation between $Ra_\ell$ and a parameter $\varLambda$ defined by

(2.34)\begin{equation} \varLambda=\left(\dfrac{\omega_M^2}{\omega_C \omega_\eta}\right)_{\!0} \sim \dfrac{V_M^2}{2 \varOmega \eta}, \end{equation}

which measures the initial ratio of the slow MC wave frequency to the magnetic diffusion frequency. Figure 7 shows that the same linear relation between $Ra_\ell$ and $\varLambda$, whose values are tabulated in table 1, holds for any $E_\eta$. Both $Ra_\ell$ and $\varLambda$ are measurable in dynamo models, the latter being of the same order of magnitude as the Elsasser number – the square of the scaled magnetic field – in many models. If the state of vanishing slow wave helicity is taken as a proxy for polarity transitions, then we may expect a self-similar relationship between $Ra_\ell$ and $\varLambda$ in low-inertia dynamos, where the nonlinear inertial force is small compared with the Coriolis force. This idea is explored further in § 3.

Figure 7. Variation of the local Rayleigh number $Ra_\ell$, defined in (2.33), shown against $\varLambda$, defined in (2.34), for the state of approximately zero slow MAC wave helicity. Three values of $E_\eta$ are considered – the diamonds represent $E_\eta =6\times 10^{-7}$, circles represent $E_\eta =6\times 10^{-6}$ and triangles represent $E_\eta =2\times 10^{-5}$.

Table 1. Values of $Ra_\ell$, defined in (2.33), at different $\varLambda$, defined in (2.34), for suppression of slow MAC waves in the linear magnetoconvection calculations. The dimensionless parameters $Le$ and $E_\eta$ are defined in (2.25a,b).

3. Nonlinear dynamo simulations

We consider a convection-driven dynamo operating in a spherical shell, the boundaries of which correspond to the inner core boundary and the core–mantle boundary. The ratio of inner to outer radius is 0.35. Fluid motion is driven by thermal buoyancy, although our formulation can also study thermochemical buoyancy using the codensity formulation (Braginsky & Roberts 1995). The other body forces acting on the fluid are the electromagnetic Lorentz force and the Coriolis force. Lengths are scaled by the thickness of the spherical shell $L$ and time is scaled by magnetic diffusion time $L^2/\eta$. The velocity $\boldsymbol {u}$ and magnetic field $\boldsymbol {B}$ are scaled by $\eta /L$ and $(2\varOmega \rho \mu \eta )^{1/2}$, respectively. The temperature is scaled by $\beta L$, where $\beta$ is the radial temperature gradient at the outer boundary. In the Boussinesq approximation, the non-dimensional MHD equations for the velocity, magnetic field and temperature are given by

(3.1)$$\begin{gather} E Pm^{{-}1} \left(\frac{\partial {\boldsymbol u}}{\partial t} + (\boldsymbol{\nabla} \times {\boldsymbol u}) \times {\boldsymbol u} \right)+ {\hat{\boldsymbol{z}}} \times {\boldsymbol u} ={-} \boldsymbol{\nabla} p^\star + Ra Pm Pr^{{-}1} T {\boldsymbol r} \nonumber\\ + \,(\boldsymbol{\nabla} \times {\boldsymbol B})\times {\boldsymbol B} + E\nabla^2 {\boldsymbol u}, \end{gather}$$

(3.2)$$\begin{gather}\frac{\partial {\boldsymbol B}}{\partial t} = \boldsymbol{\nabla} \times ({\boldsymbol u} \times {\boldsymbol B})+ \nabla^2 {\boldsymbol B}, \end{gather}$$

(3.3)$$\begin{gather}\frac{\partial T}{\partial t} +({\boldsymbol u} \boldsymbol{\cdot} \boldsymbol{\nabla}) T = Pm Pr^{{-}1} \nabla^2 T, \end{gather}$$

(3.4)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol u} = \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol B} = 0. \end{gather}$$

The modified pressure $p^*$ in (3.1) is given by $p+ \tfrac {1}{2} E Pm^{-1} |\boldsymbol {u}|^2$. The dimensionless parameters in the above equations are the Ekman number $E=\nu /2\varOmega L^2$, the Prandtl number, $Pr=\nu /\kappa$, the magnetic Prandtl number, $Pm=\nu /\eta$ and the modified Rayleigh number $g \alpha \beta L^2/2 \varOmega \kappa$. Here, $g$ is the gravitational acceleration, $\nu$ is the kinematic viscosity, $\kappa$ is the thermal diffusivity and $\alpha$ is the coefficient of thermal expansion.

The basic-state temperature profile represents a basal heating given by $T_0(r) = r_i r_o/r$, where $r_i$ and $r_o$ are the inner and outer radii of the spherical shell. The velocity and magnetic fields satisfy the no-slip and electrically insulating conditions respectively at the two boundaries. The inner boundary is isothermal while the outer boundary has constant heat flux. The calculations are performed by a pseudospectral code that uses spherical harmonic expansions in the angular coordinates $(\theta,\phi )$ and finite differences in radius $r$ (Willis et al. 2007).

As in recent studies (Varma & Sreenivasan 2022), the dynamo simulations begin from a dipole seed magnetic field of small volume-averaged intensity $\bar {B} =$ 0.01. The runs are performed for at least two magnetic diffusion times, well into the saturated state of the dynamo. The main output parameters of the dynamo simulations, given in table 2, are time-averaged values in the saturated state. For three values of the Ekman number $E$, a series of simulations at progressively increasing Rayleigh number $Ra$ are performed, spanning the dipole-dominated regime up to the onset of polarity transitions. The mean spherical harmonic degrees for convection and energy injection are defined by

(3.5a,b)\begin{equation} l_{c}= \dfrac{\varSigma l E_{k}(l)} {\varSigma E_{k}(l)}; \quad l_{E}= \dfrac{\varSigma l E_{T}(l)}{\varSigma E_{T}(l)}, \end{equation}

where $E_{k}(l)$ is the kinetic energy spectrum and $E_{T}(l)$ is the spectrum obtained from the product of the transform of $u_{r}T$ and its conjugate. The total kinetic and magnetic energies in the saturated dynamo are given by the volume integrals

(3.6a,b)\begin{equation} E_k= \dfrac{1}{2} \int \boldsymbol{u}^2 \,\text{d}V; \quad E_m= \dfrac{Pm}{2E} \int \boldsymbol{B}^2 \,\text{d}V. \end{equation}

The relative dipole field strength $f_{dip}$, which is the ratio of the mean dipole field strength to the field strength in harmonic degrees $l =$ 1–12 at the outer boundary (Christensen & Aubert 2006), takes values $>0.5$ in all the dipole-dominated runs (table 2.) The distinction between multipolar and polarity-reversing runs, however, is well understood from the evolution of the dipole colatitude, presented in § 3.1.

Table 2. Summary of the main input and output parameters in the dynamo simulations considered in this study. Here, $Ra$ is the modified Rayleigh number, $Ra_c$ is the modified critical Rayleigh number for onset of non-magnetic convection, $N_r$ is the number of radial grid points, $l_{max}$ is the maximum spherical harmonic degree, $Rm$ is the magnetic Reynolds number, $Ro_\ell$ is the local Rossby number, $l_C$ and $l_E$ are the mean spherical harmonic degrees of convection and energy injection respectively (defined in (3.5a,b)), $\bar {m}$ is the mean spherical harmonic order in the range $l \leq l_E$, $\bar {k}_s$ and $\bar {k}_z$ are the mean $s$ and $z$ wavenumbers in the range $l \leq l_E$, $E_k$ and $E_m$ are the time-averaged total kinetic and magnetic energies defined in (3.6a,b), $f_{dip}$ is the relative dipole field strength, $Ra_\ell$ is the local Rayleigh number defined in (3.16) and $\varLambda$ is the peak Elsasser number obtained from the square of the measured peak field at the earliest time of excitation of slow MAC waves in the dynamo run starting from a small seed field.

$^{a}$The last run in each Ekman number series is a polarity-reversing dynamo, for which $\varLambda$ is the square of the measured peak field when slow MAC waves cease to exist in the run starting from the saturated state of the penultimate run in that series (see also § 3.3).

For each $E$, the value of $Pm=Pr$ is chosen such that the local Rossby number $Ro_\ell$, which gives the ratio of the inertial to Coriolis forces on the characteristic length scale of convection (Christensen & Aubert 2006) is $<0.1$ (table 2). Therefore, our dynamo simulations lie in the rotationally dominant, or low-inertia, regime. While low-$Pm$ simulations starting from a seed field can result in multipolar fields even for $Ro_\ell <0.1$ (Petitdemange 2018), our choice of $Pm$ ensures that the low-inertia dynamos are dipole dominated regardless of whether one begins from a seed field or a strong field. The dipole dominance exists for a wide range of forcing, up to $Ra/Ra_c \sim 10^3$, where $Ra_c$ is the critical Rayleigh number at the onset of non-magnetic convection. The only exception here is the last run in each Ekman number series (marked by the superscript a in table 2), where a multipolar or polarity-reversing dynamo is obtained depending on the initial field strength (see § 3.1 below).

The role of the magnetic field in helicity generation is well understood by comparing a dynamo calculation with its equivalent non-magnetic convection calculation. For $Pm=Pr$, the dynamo obtained by solving (3.1)–(3.4) is compared with its non-magnetic counterpart, obtained by solving (B1)–(B3), Appendix B.

3.1. Dipolar, multipolar and polarity-reversing dynamos

Figure 8 shows the magnetic colatitude of the dipole field, $\theta$ at the upper boundary obtained from spherical harmonic Gauss coefficients, as follows:

(3.7a,b)\begin{equation} \cos\theta = g_1^0/|\boldsymbol{m}|, \quad \boldsymbol{m}= (g_1^0, g_1^1, h_1^1), \end{equation}

where $g_1^0$, $g_1^1$ and $h_1^1$ are obtained from the Schmidt-normalized expansion for the scalar potential of the field (Glatzmaier 2013, pp. 142–143). For $E=6 \times 10^{-5}$ and $Pm=Pr=5$, the evolution of $\theta$ is shown for runs at three closely spaced values of $Ra$ in strongly supercritical convection. For $Ra=20\,000$, the run that begins from a dipole seed field enters a chaotic multipolar state and subsequently regains dipolarity (figure 8a), a behaviour first noted by Sreenivasan & Kar (2018). For the same $Ra$, the run that begins from the saturated (strong field of mean square intensity $B^2 = O(1)$) state of $Ra=18\,000$ gives a stable dipole field throughout (figure 8d). Both these runs saturate at nearly identical dipolar states. For $Ra=21\,000$, the run that begins from a seed field is multipolar throughout (figure 8b). At the same $Ra$, the run that begins from a strong field produces occasional polarity reversals separated by well-defined periods of normal and reversed polarity (figure 8e), as in earlier reversing simulations (Glatzmaier et al. 1995; Kutzner & Christensen 2002) and experiments (Monchaux et al. 2009). The radial magnetic fields at the outer boundary before and after a reversal are shown in figure 9. For still higher forcing ($Ra=24\,000$), both seed field and strong field runs give multipolar solutions (figure 8c,f). Therefore, polarity reversals happen in a narrow range of $Ra$ that lies between dipolar and multipolar regimes, possibly due to variations in outer boundary heat flux.

Figure 8. Evolution of dipole colatitude with magnetic diffusion time for three closely spaced values of $Ra$ in strongly driven dynamo simulations. The top panels (a)–(c) are for simulations starting from a seed magnetic field while the bottom panels (d)–(f) are for simulations starting from a strong magnetic field. The dynamo parameters are $E = 6 \times 10^{-5},Pm=Pr=5$; (a) $Ra=20\,000$, (b) $Ra=21\,000^{{{\dagger}} }$, (c) $Ra=24\,000$, (d) $Ra=20\,000$, (e) $Ra=21\,000^{\ddagger}$, (f) $Ra=24\,000$. For $Ra=21\,000$, the superscript ${{\dagger}}$ denotes a seed field start whereas the superscript $\ddagger$ denotes a strong field start.

Figure 9. The contours of the radial magnetic field at the outer boundary for $Ra=21\,000$ at magnetic diffusion times (a) $t_d=1.1$ and (b) $t_d=2.3$ shown in figure 8(e). The other dynamo parameters are $E = 6 \times 10^{-5}$, $Pm=Pr=5$.

The range of spherical harmonic degrees $l \leq l_E$ is of particular interest since kinetic helicity is known to be generated in the nonlinear dynamo in this range of energy-containing scales (see § 1). A relative helicity is defined that measures the augmentation of lower-hemisphere helicity in the nonlinear dynamo (magnetic) run relative to that in the equivalent non-magnetic simulation

(3.8)\begin{equation} H_r = \dfrac{h_{m} - h_{nm}}{h_{nm}}, \end{equation}

for $l \leq l_E$. Here, the subscripts $m$ and $nm$ denote the magnetic (dynamo) and non-magnetic values, respectively. The variation of $H_r$ with time in dynamos at $E=6 \times 10^{-5}$ and progressively increasing $Ra$ is given in figure 10. The runs that begin from a seed field and produce an axial dipole show an approximately twofold increase in helicity (figure 10a,b). For $Ra=21\,000$, the run that begins from a seed field and produces a multipolar solution shows no noticeable increase in helicity over the non-magnetic state (figure 10c). At the same $Ra$, the run that begins from a strong field and produces polarity reversals undergoes a decrease in helicity from its initial value to the non-magnetic value (figure 10d). We emphasize that the effect of the magnetic field on helicity generation can be experienced only in the energy-containing range of harmonic degrees $l \leq l_E$. If the entire spectrum is considered (Ranjan et al. 2020), the dipolar dynamo helicity would be lower than the non-magnetic helicity.

Figure 10. Relative helicity, $H_r$, as defined in (3.8), in dynamo simulations for (a) $Ra=3000$ (seed field start), (b) $Ra=18\,000$ (seed field start), (c) $Ra=21\,000^{{\dagger}}$ (seed field start), (d) $Ra=21\,000^{\ddagger}$ (strong field start). The superposed red lines are polynomial fits showing the trend of the evolution. The dotted vertical lines in (a,b) mark the approximate times of formation of the axial dipole from a multipolar field at the outer boundary in the run. The other dynamo parameters are $E = 6 \times 10^{-5}$, $Pm=Pr=5$.

The results in figure 8 and figure 10 prompt us to examine the nature of wave motions in dipolar and reversing dynamos. An important aim of this study is to show through the analysis of MAC waves that the formation of the dipole from a multipolar state and the onset of polarity reversals lie in similar regimes. Such an argument is essential to place a constraint on the Rayleigh number that admits polarity reversals (see § 3.3).

3.2. The suppression of slow MAC waves in polarity-reversing dynamos

Isolated density disturbances in a rotating stratified fluid layer excite MAC waves whose frequencies depend on the fundamental frequencies $\omega _M$, $\omega _A$ and $\omega _C$, representing Alfvén waves, internal gravity waves and linear inertial waves, respectively. In unstable stratification that drives planetary core convection, $\omega _A^2 <0$, where $|\omega _A|$ is simply a measure of the strength of buoyancy. Since the dimensional frequencies $\omega _M^2$, $-\omega _A^2$ and $\omega _C^2$ in the dynamo are given by (Varma & Sreenivasan 2022)

(3.9a–c)\begin{equation} \omega_{M}^2 =\frac{(\boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{k})^2}{\mu \rho},\quad -\omega_{A}^2= g\alpha\beta\left(\frac{ k_{z}^2+ k_{\phi}^2}{k^2}\right),\quad \omega_{C}^2 = \frac{4(\boldsymbol{\varOmega} \boldsymbol{\cdot} \boldsymbol{k})^2}{ k^2}, \end{equation}

and scaling the frequencies by $\eta /L^2$, we obtain, in dimensionless units,

(3.10a–c)\begin{equation} \omega_M^2 = \frac{Pm}{E} (\boldsymbol{B} \boldsymbol{\cdot} {\boldsymbol k})^2, \quad -\omega_A^2 = \frac{Pm^2 Ra}{Pr \,E} \, \left(\frac{{k_z}^2 + k_{\phi}^2}{k^2}\right), \quad \omega_C^2 = \frac{Pm^2}{E^2} \frac{k_z^2}{k^2}, \end{equation}

where $k_s$, $k_{\phi }$ and $k_z$ are the radial, azimuthal and axial wavenumbers in cylindrical coordinates $(s,\phi,z)$, $k_\phi =m/s$, where $m$ is the spherical harmonic order, and $k^2=k_s^2+k_\phi ^2+k_z^2$. Here, $\omega _A$ is evaluated on the equatorial plane where the buoyancy force is maximum. The magnetic (Alfvén) wave frequency $\omega _M$ is based on the three components of the measured magnetic field at the peak-field location. The wavenumber $k_\phi$ is evaluated at $s=1$, approximately mid-radius of the spherical shell.

In the limit of zero magnetic diffusion ($\omega _\eta =0$), the characteristic equation (2.13) reduces to

(3.11)\begin{equation} \lambda^4-\lambda^2(\omega_A^2+\omega_C^2+2\omega_M^2)+\omega_A^2\omega_M^2+\omega_M^4=0, \end{equation}

the roots of which are (Sreenivasan & Maurya 2021)

(3.12)$$\begin{gather} \lambda_{1,2} ={\pm}\frac{1}{\sqrt{2}} \sqrt{\omega_A^2+\omega_C^2+2\omega_M^2+\sqrt{\omega_A^4+2\omega_A^2\omega_C^2+4\omega_M^2\omega_C^2+\omega_C^4}}, \end{gather}$$

(3.13)$$\begin{gather}\lambda_{3,4} ={\pm}\frac{1}{\sqrt{2}} \sqrt{\omega_A^2+\omega_C^2+2\omega_M^2-\sqrt{\omega_A^4+2\omega_A^2\omega_C^2+4\omega_M^2\omega_C^2+\omega_C^4}}. \end{gather}$$

For the inequality $|\omega _C| > |\omega _M| > |\omega _A|$, $\lambda _{1,2}$ and $\lambda _{3,4}$ give the frequencies of the fast ($\omega _f$) and slow ($\omega _s$) MAC waves, respectively. While the fast waves are linear inertial waves weakly modified by the magnetic field and buoyancy, the slow waves are magnetostrophic.

In figure 11, the magnitudes of the fundamental frequencies are shown as a function of the spherical harmonic order $m$ in the saturated state of the dynamo run at $E=6 \times 10^{-5}$, and $Pr=Pm=5$. The frequencies are computed from (3.10a–c) using the mean values of the $s$ and $z$ wavenumbers. For example, real-space integration over $(s,\phi )$ gives the kinetic energy as a function of $z$, the Fourier transform of which gives the one-dimensional spectrum $\hat {u}^2 (k_z)$. Subsequently, we obtain

(3.14)\begin{equation} \bar{k}_z = \frac{\varSigma k_z \hat{u}^2 (k_z)} {\varSigma \hat{u}^2 (k_z)}. \end{equation}

A similar approach gives $\bar {k}_s$. The computed frequencies in figure 11(a–c), shown for dynamos with $Ra= 6000\unicode{x2013} 18\,000$, satisfy the inequality $|\omega _C| > |\omega _M| > |\omega _A|$ in a range of the spherical harmonic order $m$. The dashed vertical lines show the value of $m$ below which the helicity in the nonlinear dynamo is greater than that in the non-magnetic run at the same parameters. The frequency root $\lambda _3$ in (3.13), shown by the black lines, represents the slow MAC waves of frequency $\omega _s$ when the inequality $|\omega _C| > |\omega _M| > |\omega _A|$ holds true. Evidently, the scales of helicity generation in the nonlinear dynamo overlaps with the scales where the slow MAC waves are generated. The range of $m$ over which the above frequency inequality holds narrows down with $Ra$, and for the polarity-reversing dynamo with $Ra=21\,000$, this inequality does not exist at any $m$ (figure 11d).

Figure 11. Panels (a i,b i,c i,d i) show the absolute values of wave frequencies plotted for the saturated state of the dynamo. The dashed vertical lines show the upper boundary of the range of wavenumbers $m$ for which the helicity of the dynamo run is greater than that of the equivalent non-magnetic run. Panels of (a ii,b ii,c ii,d ii) show the spectral distribution of the power supplied to the axial dipole, defined in (3.15). The dynamo parameters are $E=6 \times 10^{-5}$ and $Pm=Pr=5$; (a) $Ra=6000$, (b) $Ra=18\,000$, (c) $Ra=20\,000$, (d) $Ra=21\,000$.

Figure 11(a–d) also shows the spectral distribution of the power supplied to the poloidal part of the axial dipole field $B^P_{\mathrm {10}}$, given by (e.g. Buffett & Bloxham 2002)

(3.15)\begin{equation} P_{\mathrm{10}}= \int_{V} \boldsymbol {B}^{P}_{\mathrm{10}} \boldsymbol{\cdot} [\boldsymbol{\nabla} \times(\boldsymbol {u} \times \boldsymbol{ B})_{m}] \, {\rm d}V, \end{equation}

where $\boldsymbol {u}$ and $\boldsymbol {B}$ share the same value of $m$ (Bullard & Gellman 1954). The axial dipole is predominantly generated in the wavenumbers where the MAC waves are generated, except in the dynamo close to reversals (figure 11c) where the peak of $P_{10}$ lies outside the MAC wave window. In the reversing dynamo without the MAC wave window, the power supplied to the dipole is small.

In figure 12, the upper panels show the absolute values of the frequencies for a range of spherical harmonic order $m$ at three different times during the evolution of the dynamo from a seed field. The lower panels show the $m$-distribution of the power supplied to the axial dipole. As the magnetic field increases from a seed, the slow MAC waves of frequency $\omega _s$ are absent at early times but are subsequently excited in the scales where the inequality $|\omega _{C}|>|\omega _{M}|>|\omega _{A}|$ is satisfied. Kinetic helicity is generated in these scales, where a peak of the axial dipole power $P_{10}$ is also noted. In the growth phase of the dynamo where the field is multipolar, the dominant contribution to $P_{10}$ comes from higher wavenumbers of $\boldsymbol {u}$ and $\boldsymbol {B}$ within the MAC wave window (figure 12b); here, helicity is generated via slow MAC waves at the peak-field locations (figures S1 and S2 in the supplementary material available at https://doi.org/10.1017/jfm.2024.1). Once the axial dipole is established, it can hold itself up (Sreenivasan & Jones 2011) by inducing the generation of slow waves (figure 12c).

Figure 12. Panels of (a i–c i) show the absolute values of wave frequencies plotted for three different times in the growth phase of a dynamo starting from a seed field. The shaded grey area shows the range of scales where the helicity of the dynamo run is greater than that of the equivalent non-magnetic run. Panels (a ii–c ii) show the spectral distribution of the power supplied to the axial dipole, defined in (3.15). The dynamo parameters are $Ra= 18\,000, E=6 \times 10^{-5}$ and $Pm=Pr=5$. The plots are shown at times (a) $t_d=0.073$, (b) $t_d=0.168$, (c) $t_d=0.3$.

In figure 13, the dynamo frequencies are computed from (3.10a–c) using the mean values of the $s$, $\phi$ and $z$ wavenumbers in the saturated state of the dynamo in the range $l \le l_E$, which represents the energy-containing scales. The mean spherical harmonic order $\bar {m}$ is evaluated through a weighted average as in (3.14), but over the range of $m$ within $l \le l_E$. As the field increases from a small seed value in the dipolar dynamo run at $Ra=3000$ (figure 13a), slow MAC waves of frequency $\omega _s$ are first excited when $|\omega _M| > |\omega _A|$. The formation of the axial dipole from a chaotic field, marked by the dotted vertical line, follows slow wave excitation. In the run at $Ra=21\,000$ that begins from a seed field and produces a multipolar solution, $|\omega _M|$ remains lower than $|\omega _A|$ throughout (figure 13b), so the slow waves are never excited. Since polarity reversals are found in strongly driven dynamos starting from a strong field (see figure 8e), the variation of the frequencies with strength of forcing in dynamos starting from a strong field is studied in figure 13(c). Here, we note that $|\omega _M|$ falls below $|\omega _A|$ at $Ra \approx 21\,000$, which indicates that polarity reversals would indeed onset in the regime $|\omega _M| \approx |\omega _A|$ when slow MAC waves disappear. Further increase in forcing places the dynamo in a multipolar regime, the transition to which is marked by the dotted vertical line in figure 13(c). Thus, polarity reversals take place in a narrow range of $Ra$ situated between the dipolar and multipolar regimes. The volume-averaged mean square value of the axial dipole field is much smaller in the multipolar regime of $Ra=21\,000$ than in the stable dipole regime of $Ra=3000$ (figure 13d), which suggests that the slow MAC waves have an important role in the formation of the axial dipole. In the run at $Ra=21\,000$ that begins from a strong field and produces polarity reversals, the dipole intensity decreases and eventually reaches the same order of magnitude as that in the multipolar run. This is consistent with the suppression of the slow waves at this $Ra$, and reflected in the decrease in kinetic helicity to the non-magnetic value (figure 10d).

Figure 13. (a,b) Absolute values of the measured frequencies $\omega _M$, $\omega _A$, $\omega _C$ and $\omega _s$ plotted against time (measured in units of the magnetic diffusion time, $t_d$) for $Ra=3000$ (a) and $Ra=21\ 000$ (b). These simulations study the evolution of the dynamo starting from a small seed magnetic field. The dotted vertical line in (a) marks the time of formation of the axial dipole from a multipolar field. (c) Frequencies shown against $Ra$ for saturated dynamos starting from a strong field. The dashed vertical line indicates the value of $Ra$ at which the slow MAC wave frequency $\omega _s$ goes to zero as $\omega _M \approx \omega _A$. This marks the transition from the dipolar regime to the polarity-reversing regime. The dotted vertical line marks the beginning of the multipolar regime. (d) The Elsasser number of the axial dipole field component, based on its root mean square value, for $Ra=3000$, $Ra=21\,000^{{\dagger}}$ (seed field start) and $Ra=21\,000^{\ddagger}$ (strong field start). The dynamo parameters are $E = 6 \times 10^{-5}$, $Pm=Pr=5$. The colour codes in (a) are also used in (b,c).

Figure 14 shows the measurement of wave motion in the saturated state of dynamos at $E=6 \times 10^{-5}$ and $Pr=Pm=5$ and three values of $Ra$ spanning the dipole-dominated regime and reversals. Contours of $\dot {u}_z$ at cylindrical radius $s=1$ are plotted over small time windows in which the ambient magnetic field and wavenumbers are approximately constant. These contours show the propagation paths of the fluctuating $z$ velocity. In line with the discussion so far, the measurement of axial motions is limited to the energy-containing scales $l \leq l_E$, with no restriction on the wavenumber. The measured axial group velocity of the waves, $U_{g,z}$ – obtained from the slope of the black lines in figure 14 – is compared with the estimated fast ($U_f$) and slow ($U_s$) group velocities obtained by taking the derivatives of the respective frequencies in (3.12) and (3.13) with respect to $k_z$ (table 3). The theoretical frequencies $\omega _f$ and $\omega _s$ are estimated using the three components of the magnetic field at the peak-field location and the mean values of $k_s$, $k_z$ and $m$ over the range of energy-containing scales, $l \leq l_E$. In the dipole-dominated run at $Ra=6000$, the slow and fast waves coexist, but the intensity of slow wave motions measured at the peak-field locations is at least as high as that of the fast wave motions. At $Ra=20\,000$, the increasing occurrence of the fast waves alongside the slow waves is noted from the nearly vertical propagation paths (figure 14b). The measured slow wave group velocity at $Ra=20\,000$ is greater than that at $Ra=6000$, which reflects the larger self-generated field at the higher Rayleigh number. In the reversing dynamo at $Ra=21\,000$, slow waves are totally absent while the fast waves are abundant (figure 14c).

Figure 14. Contour plots of $\partial {u}_z/\partial t$ at cylindrical radius $s=1$ for $l\leq l_E$ and small intervals of time in the saturated state of three dynamo simulations. The parallel black lines indicate the predominant direction of travel of the waves and their slope gives the group velocity $U_{g,z}$. The dynamo parameters are $E=6 \times 10^{-5}$, $Pr=Pm=5$; (a) $Ra=6000$, (b) $Ra=20\ 000$, (c) $Ra= 21\ 000$. The estimated group velocity of the fast and slow MAC waves ($U_f$ and $U_s$ respectively) and $U_{g,z}$ are given in table 3.

Table 3. Summary of the data for MAC wave measurement in the dynamo models. The sampling frequency $\omega _n$ is chosen to ensure that the fast MAC waves are not missed in the measurement of group velocity. The values of $\omega _M^2$, $-\omega _A^2$ and $\omega _C^2$ are calculated from (3.10a–c) using the mean values of $m$, $k_s$ and $k_z$ over the range of energy-containing scales, $l \leq l_E$. The measured group velocity in the $z$ direction ($U_{g,z}$) is compared with the estimated fast ($U_f$) or slow ($U_s$) MAC wave velocity.

The dominance of the slow MAC waves in the dipole-dominated dynamo and the fast MAC waves in the reversing dynamo is further evident in figure 15, where the FFT of $\dot {u}_z$ is shown. The flow largely consists of waves of frequency $\omega \sim \omega _s$ in the dipolar dynamo (figure 15a), whereas in the reversing dynamo, waves of much higher frequency $\omega \sim \omega _f$ are dominant (figure 15b). The coloured vertical lines in figure 15(a,b) give the estimated magnitudes of $\omega _C$ and $\omega _M$ normalized by $\omega _f$, where $\omega _M$ is based on the peak magnetic field.

Figure 15. (a) The FFT spectrum of $\partial u_z/\partial t$ at cylindrical radius $s = 1$ for the scales $l \le l_E$. The spectra are computed at discrete $\phi$ points and then averaged azimuthally. The solid black vertical lines in (a) and (b) correspond to $\omega /\omega _f=1$ and the dotted vertical line in (a) corresponds to $\omega _s^\star =\omega _s/\omega _f$, where $\omega _f$ and $\omega _s$ are the estimated fast and slow MAC wave frequencies. The coloured vertical lines give $\omega _C$ and $\omega _M$, both normalized by $\omega _f$. The dynamo parameters are $E = 6 \times 10^{-5}$ and $Pm=Pr=\,$5; (a) $Ra=6000$, (b) $Ra=21\,000$.

3.3. Self-similarity of the dipole–multipole transition

In the presence of small but finite magnetic diffusion, the slow MAC waves are known to disappear in an unstably stratified medium for $|\omega _A/\omega _M| \approx 1$. The same condition must hold for the appearance of slow waves in a medium where the magnetic field progressively increases from a small value. In simulations starting from a small seed field, the earliest time of excitation of the slow MAC waves is noted from group velocity measurements at closely spaced times during the growth phase of the dynamo. The peak Elsasser number $\varLambda =B_{peak}^2$ at this time is obtained from the three components of the field at the peak-field location, and presented in the last column of table 2. In each of the three dynamo series considered in this study, the last run is a polarity-reversing dynamo, for which $\varLambda$ is the measured peak Elsasser number when slow MAC waves cease to exist in the run starting from the saturated (strong field) state of the penultimate run in that series. Figure 16(a) shows that the variation of $Ra$ with $\varLambda$ in the three dynamo series is nearly linear. The Rayleigh number corresponding to reversals (at which slow MAC waves disappear) also lies on this line, indicating that the appearance and disappearance of MAC waves are in similar regimes.

Figure 16. (a) Variation of the modified Rayleigh number $Ra$ with the peak Elsasser number $\varLambda$ (square of the peak magnetic field) at both excitation and suppression of slow MAC waves. The hollow symbols represent the states where slow waves are first excited as the dynamos evolve from a small seed field; the filled symbols represent the states where slow waves are suppressed in the polarity-reversing dynamos. The parameters of the three dynamo series and their symbolic representations are as follows: $E=3\times 10^{-4},Pm=Pr=20$ (circles), $E=6\times 10^{-5},Pm=Pr=5$ (triangles), $E=1.2\times 10^{-5},Pm=Pr=1$ (diamonds). (b) Variation of the local Rayleigh number $Ra_\ell$, defined in (3.16), with $\varLambda$. The values of $Ra$, $Ra_\ell$ and $\varLambda$ in the plots are given in table 2. The sections 1 and 2 marked on the self-similar line are analysed further in figures 17 and 18 below.

Following the analysis in figure 7, where the Rayleigh number based on the length scale of the buoyant perturbation was studied, we define a local Rayleigh number in the dynamo

(3.16)\begin{equation} Ra_\ell = \dfrac{g \alpha \beta}{2 \varOmega \eta}\left(\dfrac{2 {\rm \pi}}{\bar{m}}\right)^2, \end{equation}

where $\bar {m}$ is the mean spherical harmonic order evaluated over $m$ within the energy-containing scales $l \leq l_E$. The behaviour of $Ra_\ell$, which is defined for the scales where the MAC waves are excited by buoyancy, is approximately self-similar (figure 16b). The values of $Ra_\ell$ at the onset of polarity reversals, where the slow MAC waves disappear, also lie on the same self-similar branch. While the conventional Rayleigh numbers $Ra$ at the onset of reversals in the three dynamo series lie far apart (see the filled symbols in figure 16a), the respective local Rayleigh numbers $Ra_\ell$ are remarkably close, and $\sim 10^4$ (see the filled symbols in figure 16(b) and table 2). The magnetic Ekman number based on $\bar {m}$ in the three dynamo series takes values of $E_\eta \sim 10^{-5}$ for a wide range of $Ra$, which indicates an energy-containing length scale $\sim$ 10 km for Earth (see § 1).

The fact that the linear variation of $Ra_\ell$ with $\varLambda$ demarcates the boundary between dipolar and multipolar states is evident by traversing the sections 1 and 2 marked on the $Ra_\ell$ line from left to right (figure 16b). In practice, this is done by following the evolution of the dynamo from a small seed field. Figure 17(a) shows the section 1 within dashed vertical lines, where $|\omega _M|$ crosses $|\omega _A|$. The variation of the dipole colatitude $\theta$, shown in figure 17(b), indicates a multipolar field until this crossing, and a stable dipole thereafter. While the flow is predominantly made up of fast MAC waves of frequency $\omega \sim \omega _f$ before the transition, the slow waves of frequency $\omega \sim \omega _s$ are dominant after the transition (figure 17c,d). The multipole–dipole transitions are further evident in the contour plots of the radial magnetic field at the outer boundary, given in figure 18 for sections 1 and 2 marked in figure 16(b).

Figure 17. (a) Evolution the dynamo frequencies in a simulation beginning from a small seed magnetic field at $E = 1.2 \times 10^{-5}$, $Pm=Pr=1$, $Ra=4000$. The two dashed vertical lines at $t_d=0.55$ and $t_d=0.66$ represent the end points of section 1 marked in figure 16(b) with peak Elsasser numbers $\varLambda = 32$ and 56, respectively. (b) Evolution of the dipole colatitude in the above simulation. (c,d) FFT spectra of $\partial u_z/\partial t$ at cylindrical radius $s = 1$ for the scales $l \le l_E$ at $\varLambda = 32$ and 56, respectively. The spectra are computed at discrete $\phi$ points and then averaged azimuthally. In (d), $\omega _s^* = \omega _s/\omega _f$, where $\omega _f$ and $\omega _s$ are the estimated fast and slow MAC wave frequencies.

Figure 18. Shaded contours of the radial magnetic field at the outer boundary in two dynamo simulations: (a,b) $Ra_\ell =2505$, corresponding to section 1 in figure 16(b), at $\varLambda = 32$ and 56, respectively. The dynamo parameters are $E=1.2\times 10^{-5}$, $Ra=4000$, $Pm=Pr=1$; (c,d) $Ra_\ell =11\,740$, corresponding to section 2 in figure 16(b), at $\varLambda =$ 187 and 216, respectively. The dynamo parameters are $E=6\times 10^{-5}$, $Ra=18\,000$, $Pm=Pr=5$.

4. Concluding remarks

The present study investigates the dipole–multipole transition in rotating dynamos through the analysis of MHD wave motions. The limit of small Rossby number, based not only on the planetary core depth but also on the length scale of core convection, is considered. In this inertia-free limit, the dynamo polarity depends on the relative magnitudes of $\omega _M$ and $\omega _A$, which in turn depend on the peak intensity of the self-generated field and the strength of buoyant forcing in the unstably stratified fluid layer. The linear magnetoconvection model shows that the slow MAC (magnetostrophic) waves have greater helicity than the fast MAC waves for $|\omega _C| > |\omega _M|> |\omega _A| > |\omega _\eta |$, the regime thought to be relevant to dipole-dominated dynamos. Although buoyancy-induced helicity can be spatially segregated about the equator by linear inertial waves in the absence of the magnetic field (Davidson & Ranjan 2018), this purely hydrodynamic process is not likely to generate the axial dipole in rotating dynamos, where the helicity of slow MAC waves is at least as high as that of the fast waves of frequency $\sim \omega _C$, which already exist in the multipolar state (figure S2 in the supplementary material). Moreover, the time scale in which the slow waves establish the dipole field can be a significant fraction of the magnetic diffusion time (see figure 10a,b), whereas the fast waves would have to form the dipole on a much shorter time scale. Since the slow MAC waves are produced by magnetostrophic balances at the peak-field locations, the analysis of wave motions with $\omega _M$ based on the peak-field intensity is performed here rather than a scale-dependent analysis of volume-averaged forces (e.g. Schwaiger, Gastine & Aubert 2019). The present study shows that the condition of vanishing slow waves defines polarity transitions in rapidly rotating dynamos. The variation of the local Rayleigh number $Ra_\ell$ with the Elsasser number $\varLambda$ at the transition is very similar in the linear magnetoconvection and dynamo models, which is remarkable given that the Alfvén frequency $\omega _M$ is determined by the imposed field in the linear model and by the self-generated field in the dynamo model.

While the onset of slow magnetostrophic waves for $|\omega _C| > |\omega _M| \approx |\omega _A| > |\omega _\eta |$ is known to produce the axial dipole from a chaotic multipolar state (see figures 17a,b and 18), the suppression of slow waves for the same condition leads to polarity reversals in strongly driven dynamos starting from a strong field state of mean square intensity $O(1)$ (figure 13c). This study differentiates between polarity-reversing and multipolar solutions, both of which onset at the same Rayleigh number. Polarity reversals result from the suppression of slow waves that already exist in a strong field state, whereas the multipolar solution forms from a seed field when the forcing is sufficiently strong that the slow waves are never excited. Multipolar solutions are also obtained at stronger forcing, regardless of whether one begins from a seed field or a strong field.

The critical ratios of forces and energies in the dynamo have been instrumental in setting criteria for polarity transitions. Earlier studies (Christensen & Aubert 2006; Olson & Christensen 2006) have proposed a critical value $\approx 0.1$ for the local Rossby number $Ro_\ell$, above which the dipole–multipole transition occurs. However, this condition relies on the finiteness of nonlinear inertia on the scale of the buoyancy perturbations, which we consider unlikely in rapidly rotating planetary cores. While the condition of the energy ratio $E_k/E_m >1$ (table 2) by itself is admissible as a criterion for the dipole–multipole transition (Kutzner & Christensen 2002; Tassin et al. 2021), the increase of this ratio is likely a consequence of the transition rather than a cause for it. Polarity reversals are never obtained for $E_k/E_m <1$, so it is not certain that the reversals in the present study are Earth like.

The self-similarity of polarity reversals in the inertia-free regime can place a useful constraint on the Rayleigh number that admits reversals in the Earth. For $Ra_\ell \sim 10^4$, the classical Rayleigh number $R=Ra/E \sim 10^{17}$, taking a turbulent viscosity $\nu \approx \eta$ and a plausible ratio of core depth to the convective length scale, $L/\delta \sim 10^2$. That having been said, the Earth's core is thought to convect in response to large lateral variations in lower-mantle heat flux (e.g. Olson et al. 2015; Mound et al. 2019), which can potentially reduce the value of $Ra$ for reversals. Our understanding of the convective regime for polarity reversals is not complete, but the idea that polarity transitions are self-similar would eventually lead to improved parameter constraints for reversals.