2 The magnetostrophic approximation

Magnetostrophy (abbreviated MS) is a term we shall use to mean the absence of viscous friction in Earth’s core. This approximation, used throughout this paper, is suggested by the extreme smallness of the ratio of the viscous and Coriolis forces, quantified by the Ekman number $E=\unicode[STIX]{x1D708}/2\unicode[STIX]{x1D6FA}r_{o}^{2}$ , where $\unicode[STIX]{x1D708}$ is the core’s kinematic viscosity, $\unicode[STIX]{x1D6FA}$ ( ${\approx}7.3\times 10^{-5}~\text{s}^{-1}$ ) is the angular speed of Earth’s rotation and $r_{o}$ ( ${\approx}3.48\times 10^{6}~\text{m}$ ) is the radius of the core, assumed spherical. For $\unicode[STIX]{x1D708}\approx 10^{-6}~\text{m}^{2}~\text{s}^{-1}$ , a value widely believed to be realistic for Earth’s core, $E\approx 10^{-15}$ ; e.g. de Wijs et al. (Reference de Wijs, Kresse, Voĉado, Dobson, Alfè, Gillan and Price1998), Voĉadlo et al. (Reference Voĉadlo, Alfè, Price and Gillan2000).

It is generally recognized that core motions are turbulent and, lacking a proper deductive theory of turbulence, several ad hoc representations of turbulent stresses have been made, including ones that simply replace the kinematic viscosity $\unicode[STIX]{x1D708}$ by a turbulent viscosity, $\unicode[STIX]{x1D708}_{\text{turb}}$ . Although this certainly oversimplifies turbulent dynamics (see e.g. Roberts & King Reference Roberts and King2013, § 9.2), it has been used to make the estimates $E_{\text{turb}}\approx 6\times 10^{-11}$ (Buffett Reference Buffett1992) and $E_{\text{turb}}\approx 2\times 10^{-12}$ (Buffett & Christensen Reference Buffett and Christensen2007). We believe that these turbulent Ekman numbers are still so small that the MS approximation is justified.

The fully three-dimensional MHD models that today dominate dynamo studies first appeared in the mid-1990s. We shall call them ‘today’s models’. Their impact has been so great that interest in the MS alternative has waned, even though demands for numerical stability condemn investigators to limit themselves to $E$ vastly greater than the geophysically realistic $10^{-15}$ . When $E$ has to exceed ${\sim}10^{-6}$ in simulations, Roberts & King (Reference Roberts and King2013) called them ‘viscously dominated’. They pointed out that, even if the Earth’s core were made of honey, $E$ would still be less than $10^{-11}$ !

A main preoccupation of dynamo simulation over the past 2 decades has been to reduce the gap between the $E$ of today’s models and realistic $E$ , using the ever-increasing capability of computer hardware. Recently Schaeffer et al. (Reference Schaeffer, Jault, Nataf and Fournier2017) were able to reach $E=3\times 10^{-7}$ in buoyantly driven three-dimensional (3-D) dynamos. Livermore, Jackson & Hollerbach (Reference Livermore, Jackson and Hollerbach2013) reached $E=3\times 10^{-9}$ in a limited simulation that tracked only the toroidal part of the magnetic field. Although it is remarkable what has been achieved in 2 decades, the drive for smaller $E$ becomes increasingly difficult as $E$ is decreased; see Davies, Gubbins & Jimack (Reference Davies, Gubbins and Jimack2011). Using Moore’s law, Roberts & King (Reference Roberts and King2013, § 9.1) estimated that advances in computer technology will be unable before AD 2038 to simulate buoyantly driven 3-D MHD dynamos for $E\approx 10^{-15}$ . The ever-increasing difficulty of attaining ever-smaller $E$ is caused principally by the need to resolve numerically the Ekman layers at the core boundaries, whose thicknesses are dimensionally $\unicode[STIX]{x1D6FF}^{E}\propto E^{1/2}r_{o}$ . For the numerically attainable $E=10^{-6}$ , $\unicode[STIX]{x1D6FF}^{E}$ is approximately 3 km; for the unattainable $E=10^{-15}$ , $\unicode[STIX]{x1D6FF}^{E}$ is about 10 cm. The rate ${\mathcal{Q}}_{\unicode[STIX]{x1D708}}$ of energy dissipation through viscosity is (dimensionally) the product of the volume of the boundary layer $O(r_{o}^{2}\unicode[STIX]{x1D6FF}^{E})$ and the dissipation rate per unit volume $O(\unicode[STIX]{x1D70C}\unicode[STIX]{x1D708}({\mathcal{U}}/\unicode[STIX]{x1D6FF}^{E})^{2})$ , where ${\mathcal{U}}$ is the characteristic speed of core motions; therefore ${\mathcal{Q}}_{\unicode[STIX]{x1D708}}\propto 2\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D70C}{\mathcal{U}}^{2}r_{o}^{2}\unicode[STIX]{x1D6FF}^{E}$ . In other words, as $E$ is reduced, the influence of the boundary layers on the mainstream solution (the solution at distances greater than $\unicode[STIX]{x1D6FF}^{E}$ from boundaries) diminishes until plausibly it is effectively non-existent. The boundary layers themselves become vanishingly thin. Why then continually press today’s models to deliver solutions for smaller and smaller $E$ ? Why not follow the example of every other modeller in theoretical physics: if a physical process is unimportant, exclude it from the model?! Although this may motivate magnetostrophic modelling strongly, it is insufficient; more should be said about the limitations of MS.

Today’s computer models have one significant attribute lacked by MS: they include the inertial force in full. MS solutions must, however, approximate the inertial force to avoid the large Ekman numbers that make today’s computer models irrelevant to the geodynamo and to most cosmic dynamos. This subject is the main topic of § 4.1 below. The MS approximation is often viewed as the logical way to find solutions when viscosity is unimportant, but in reality its success also depends indirectly on the smallness of the inertial force, as measured by the Rossby number, $Ro$ , defined in § 3.1 below.

It is not appropriate here to prolong discussion by including the inner core and the shear layers that plausibly exist at the tangent cylinder, which is the cylinder parallel to $Oz$ that touches the inner core at its equator. There is an extensive literature on these shear layers, going back to Proudman (Reference Proudman1956) and Stewartson (Reference Stewartson1966); see Dormy & Soward (Reference Dormy and Soward2007). Most MS models, including ours, do not recognize the existence of the inner core, but see Hollerbach (Reference Hollerbach1994) and Livermore & Hollerbach (Reference Livermore and Hollerbach2012).

It is natural to hope that today’s dynamos for $E\lesssim 10^{-6}$ will differ little from those operating at $E\approx 10^{-15}$ , but there are plausible reasons for thinking otherwise. These doubts followed the analysis of Chandrasekhar (Reference Chandrasekhar1954, Reference Chandrasekhar1956) of the marginal stability of a rotating Bénard layer in a vertical magnetic field, and were first articulated by Eltayeb & Roberts (Reference Eltayeb and Roberts1970) who argued that, when $E$ is small enough, any sufficiently strong magnetic field would facilitate (magneto-)convection so much that magnetic field generation would even be possible at smaller Rayleigh numbers $Ra$ than when the magnetic field is smaller or absent, i.e. they postulated two types of buoyantly driven dynamos, a weak-field, viscously dominated branch and a more easily excited, strong-field, magnetically dominated branch. (The Rayleigh number is a non-dimensional measure of the buoyancy force when temperature decreases with increasing height. Chandrasekhar was the first to give it a name. In chapter VI of his celebrated book (Chandrasekhar Reference Chandrasekhar1961), he defined $Ra$ for a sphere of radius $r_{o}$ uniformly heated by dissolved heat sources. Assuming a constant thermal diffusivity, $\unicode[STIX]{x1D705}$ , the ambient temperature gradient is $\unicode[STIX]{x1D737}=-\unicode[STIX]{x1D6FD}_{o}\boldsymbol{r}/r_{o}$ , where $\unicode[STIX]{x1D6FD}_{o}$ is the gradient at $r=r_{o}$ . For uniform fluid density, the gravitational force is $\boldsymbol{g}(\boldsymbol{x})=-g_{o}\boldsymbol{r}/r_{o}$ , where $g_{o}$ is gravity at $r=r_{o}$ . Chandrasekhar defined the Rayleigh number as $Ra=g_{o}\unicode[STIX]{x1D6FC}_{T}\unicode[STIX]{x1D6FD}_{o}r_{o}^{4}/\unicode[STIX]{x1D708}\unicode[STIX]{x1D705}$ , where $\unicode[STIX]{x1D6FC}_{T}$ is the thermal expansion coefficient. For a highly rotating system, particularly in the MS approximation, a more convenient definition is $Ra=g_{o}\unicode[STIX]{x1D6FC}_{T}\unicode[STIX]{x1D6FD}_{o}r_{o}^{2}/\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D702}$ .)

The speculation of Eltayeb & Roberts (Reference Eltayeb and Roberts1970) was soon supported by Soward (Reference Soward1974), who derived weak-field solutions for magneto-convection in the Bénard layer for $E\ll 1$ . He found that, as the Rayleigh number is increased, these terminated in an asymptote at finite $B$ . It was conjectured that this joined to the strong-field branch; see his figure 5. Dormy (Reference Dormy2016) recently found that both branches exist when $E=3\times 10^{-4}$ ; see his figure 3. Although this encourages further study of MS magneto-convection, our interest below is with mean-field dynamos driven by the $\unicode[STIX]{x1D6FC}$ -effect alone. There is no obvious physical reason why there should be two branches of $\unicode[STIX]{x1D6FC}^{2}$ dynamos, and neither this paper nor papers 1 and 2 indicate that two mean-field branches exist.

Because we ignore viscosity, $\boldsymbol{U}$ can satisfy only one condition on $S$ . This is $U_{r}(1,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719},t)=0$ , which is applied separately on the two hemispheres of $S$ :

(2.1a,b ) $$\begin{eqnarray}\displaystyle U_{r}^{{\mathcal{N}}}(1,\unicode[STIX]{x1D703},t)=0,\quad U_{r}^{{\mathcal{S}}}(1,\unicode[STIX]{x1D703},t)=0. & & \displaystyle\end{eqnarray}$$

Generally $U_{\unicode[STIX]{x1D703}}(1,\unicode[STIX]{x1D703},t)$ and $U_{\unicode[STIX]{x1D719}}(1,\unicode[STIX]{x1D703},t)$ are non-zero.

There is neither viscous coupling between core and mantle in the MS approximation nor magnetic coupling, in the assumed absence of sources of magnetic field in the mantle. We could easily include other coupling mechanisms (topographic, gravitational, viscous, $\ldots$ ; see e.g. Roberts & Aurnou Reference Roberts and Aurnou2012) but do not do so because we want to avoid complicating further an already complicated message. Therefore $\unicode[STIX]{x1D6E4}_{z}^{\text{total}}$ , the $z$ -component of the total torque on v, is zero so that the angular momentum, $m_{z}^{\text{total}}$ , of core motion about $Oz$ cannot change. Our reference frame is the rotating frame in which

(2.2) $$\begin{eqnarray}\displaystyle m_{z}^{\text{total}}\equiv \int _{\text{v}}sU_{\unicode[STIX]{x1D719}}\,\text{d}^{3}x=0. & & \displaystyle\end{eqnarray}$$

We start our numerical integrations from the state

(2.3a,b ) $$\begin{eqnarray}\displaystyle \boldsymbol{B}(\boldsymbol{x},0)=B_{0}[3sz\mathbf{1}_{s}+(5-6s^{2}-3z^{2})\mathbf{1}_{z}],\quad \boldsymbol{U}(\boldsymbol{x},0)=\mathbf{0}, & & \displaystyle\end{eqnarray}$$

where $B_{0}$ is an arbitrary constant representing magnetic field strength. This magnetic field is purely dipolar and is created by a zonal current density $J_{\unicode[STIX]{x1D719}}$ proportional to $s$ , so that the Lorentz force is $\boldsymbol{J}\times \boldsymbol{B}=\unicode[STIX]{x1D735}[(1/2)(3B_{0}s)^{2}(5-3r^{2})]$ . It is therefore in hydrostatic balance with constant $\unicode[STIX]{x1D6F1}+(1/2)(3B_{0}s)^{2}(5-3r^{2})$ , making it what in § 3.1 will be called a ‘Taylor state’. Clearly (2.1a,b ) are satisfied by (2.3b ), which also implies $m_{z}^{\text{total}}=0$ .

A short digression on equatorial symmetry may be welcome here. Although zero in (2.3a ), a non-zero $B_{\unicode[STIX]{x1D719}}(s,z,t)$ is generated in the evolution of $\boldsymbol{B}$ and $\boldsymbol{U}$ that we study below. This, like $B_{s}(s,z,t)$ , is antisymmetric in $z$ : i.e. $B_{\unicode[STIX]{x1D719}}(s,-z,t)\equiv -B_{\unicode[STIX]{x1D719}}(s,z,t)$ . The antisymmetry of $B_{s}$ and $B_{\unicode[STIX]{x1D719}}$ is preserved as $\boldsymbol{B}$ evolves, as is the $z$ -symmetry of $B_{z}$ : $B_{z}(s,-z,t)\equiv B_{z}(s,z,t)$ . This type of $z$ -symmetry is termed ‘dipole symmetry’. Also, starting from (2.3b ), $\boldsymbol{U}$ becomes non-zero and quadrupolar, with $z$ -symmetric $s$ - and $\unicode[STIX]{x1D719}$ -components and a $z$ -antisymmetric $z$ -component. In addition to this family of solutions in which $\boldsymbol{B}$ is permanently dipolar and $\boldsymbol{U}$ permanently quadrupolar, there is another family in which both $\boldsymbol{B}$ and $\boldsymbol{U}$ are permanently quadrupolar. Neither that family nor mixed dipolar/quadrupolar states are considered here.

The adjectives geostrophic and ageostrophic often appear in papers 1 & 2. There, as in almost all theoretical work in this area including this paper, the system is homogeneous and v is axisymmetric with respect to the rotation axis. (The axisymmetry of $\boldsymbol{B}$ , assumed in paper 2 and assumed in most mean-field models results in further simplification.) Axisymmetry of v simplifies the definition of geostrophy; see Greenspan (Reference Greenspan1968). The geostrophic cylinder ${\mathcal{C}}(s)$ of radius $s$ is defined by $0\leqslant \unicode[STIX]{x1D719}<2\unicode[STIX]{x03C0}$ and by $-z_{1}(s)\leqslant z\leqslant z_{1}(s)$ where $z_{1}(s)=\surd (1-s^{2})$ . Its area is ${\mathcal{A}}(s)=4\unicode[STIX]{x03C0}sz_{1}(s)$ . The cylinder ${\mathcal{C}}(s)$ is part of the boundary of a volume ${\mathcal{V}}(s)$ consisting of the interior of v together with northern and southern extensions bounded by spherical caps ${\mathcal{N}}(s)$ and ${\mathcal{S}}(s)$ which are segments of S that complete the boundary of ${\mathcal{V}}(s)$ . The geostrophic average of a scalar function $Q(\boldsymbol{x})$ is

(2.4a ) $$\begin{eqnarray}\displaystyle Q^{\text{G}}(s)=\frac{1}{2z_{1}(s)}\int _{-z_{1}(s)}^{z_{1}(s)}Q(s,z)\,\text{d}z. & & \displaystyle\end{eqnarray}$$

As $s\rightarrow 1$ , $Q(s,z)\rightarrow Q(1,0)$ so that $Q^{\text{G}}(s)\rightarrow Q(1,0)$ . If $Q(s,z)\propto z^{2}$ as $s\rightarrow 1$ , then $Q^{\text{G}}(s)\propto z_{1}^{2}$ , etc.

The ageostrophic (or non-geostrophic) part of $Q(\boldsymbol{x})$ is what remains after the geostrophic part $Q^{\text{G}}(s)$ has been removed:

(2.4b ) $$\begin{eqnarray}\displaystyle Q^{\text{A}}(\boldsymbol{x})=Q(\boldsymbol{x})-Q^{\text{G}}(s). & & \displaystyle\end{eqnarray}$$

The geostrophic and ageostrophic parts of a vector such as $\boldsymbol{U}(\boldsymbol{x},t)$ singles out its $\unicode[STIX]{x1D719}$ -component:

(2.4c ) $$\begin{eqnarray}\displaystyle \boldsymbol{Q}^{\text{A}}(\boldsymbol{x})=\boldsymbol{Q}(\boldsymbol{x})-\boldsymbol{Q}^{\text{G}}(s)=\boldsymbol{Q}(s,z)-Q_{\unicode[STIX]{x1D719}}^{\text{G}}(s)\mathbf{1}_{\unicode[STIX]{x1D719}}, & & \displaystyle\end{eqnarray}$$

so that, for example,

(2.4d-f ) $$\begin{eqnarray}\displaystyle U_{s,z}^{\text{A}}(\boldsymbol{x},t)=U_{s,z}(s,z,t),\quad U_{\unicode[STIX]{x1D719}}^{\text{A}}(\boldsymbol{x},t)=U_{\unicode[STIX]{x1D719}}(s,z,t)-U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t),\quad [U_{\unicode[STIX]{x1D719}}^{\text{A}}(\boldsymbol{x},t)]^{\text{G}}(s,t)=0. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

3 Basic equations and boundary conditions

Finding the magnetic field $\boldsymbol{B}$ created by the fluid velocity $\boldsymbol{U}$ is a standard kinematic dynamo problem that has been so often solved that we shall not dwell on it here. We shall merely, for reference purposes, state that when the $\unicode[STIX]{x1D6FC}$ -effect is included, $\boldsymbol{B}$ and $\boldsymbol{E}$ satisfy the pre-Maxwell equations. In v

(3.1a-d ) $$\begin{eqnarray}\displaystyle \boldsymbol{J}=\unicode[STIX]{x1D735}\times \boldsymbol{B},\quad \unicode[STIX]{x1D735}\times \boldsymbol{E}=-\unicode[STIX]{x2202}_{t}\boldsymbol{B},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}=0,\quad \boldsymbol{J}=\boldsymbol{E}+\boldsymbol{U}\times \boldsymbol{B}+\unicode[STIX]{x1D6FC}\boldsymbol{B}. & & \displaystyle\end{eqnarray}$$

The induction equation that follows from (3.1a–d ) is

(3.1e ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{B}}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D735}\times (\boldsymbol{U}\times \boldsymbol{B}+\unicode[STIX]{x1D6FC}\boldsymbol{B})+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{B}. & & \displaystyle\end{eqnarray}$$

In the electrically insulating exterior $\text{v}^{\ast }$ of v, (3.1a–c ) hold but (3.1d ) is replaced by $\boldsymbol{J}^{\ast }\equiv \mathbf{0}$ .

The solution to (3.1e ) has to obey the continuity condition,

(3.1f,g ) $$\begin{eqnarray}\displaystyle \boldsymbol{B}=\boldsymbol{B}^{\ast }\quad \text{on }r=1,\quad \text{and}\quad r^{2}\boldsymbol{B}^{\ast }\rightarrow \mathbf{0}\quad \text{as }r\rightarrow \infty , & & \displaystyle\end{eqnarray}$$

the latter excluding sources of $\boldsymbol{B}$ at infinity. Equations (3.1f,g ) lead to the so-called ‘dynamo conditions’, which for axisymmetric $\boldsymbol{B}$ , are

(3.1h,i ) $$\begin{eqnarray}\displaystyle B_{\unicode[STIX]{x1D719}}(1,\unicode[STIX]{x1D703},t)=0,\quad \left[\frac{\unicode[STIX]{x2202}B_{r,n}}{\unicode[STIX]{x2202}r}+(n+2)B_{r,n}\right]_{(1,t)}=0\quad \forall n\geqslant 0, & & \displaystyle\end{eqnarray}$$

where $B_{r,n}(r,t)$ is the $n$ th spherical harmonic components of $B_{r}(\boldsymbol{x},t)$ .

In § 4, the geostrophic and ageostrophic parts of the induction (3.1e ) will be needed. These are defined by

(3.2a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{B}(\boldsymbol{x},t)=\boldsymbol{C}^{\text{G}}(\boldsymbol{x},t)+\boldsymbol{C}^{\text{A}}(\boldsymbol{x},t),\quad \text{where} & & \displaystyle\end{eqnarray}$$

(3.2b,c ) $$\begin{eqnarray}\displaystyle \boldsymbol{C}^{\text{G}}(\boldsymbol{x},t)=\unicode[STIX]{x1D735}\times (\boldsymbol{U}^{\text{G}}\times \boldsymbol{B}),\quad \boldsymbol{C}^{\text{A}}(\boldsymbol{x},t)=\unicode[STIX]{x1D735}\times (\boldsymbol{U}^{\text{A}}\times \boldsymbol{B}+\unicode[STIX]{x1D6FC}\boldsymbol{B})+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{B}. & & \displaystyle\end{eqnarray}$$

In the axisymmetric case, it was shown in (2: equation (18d)) that

(3.2d-f ) $$\begin{eqnarray}\displaystyle \boldsymbol{C}^{\text{G}}(\boldsymbol{x},t)=sZ^{\text{G}}(s,t)B_{s}(\boldsymbol{x},t)\mathbf{1}_{\unicode[STIX]{x1D719}},\quad \text{where }Z^{\text{G}}(s,t)=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}(s,t)}{\unicode[STIX]{x2202}s}\text{ and }\unicode[STIX]{x1D701}(s,t)=\frac{U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)}{s} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

is the geostrophic shear.

In MS theory, the Boussinesq equations governing $\boldsymbol{U}(\boldsymbol{x},t)$ in mean-field models satisfies

(3.3a-c ) $$\begin{eqnarray}\displaystyle E_{\unicode[STIX]{x1D702}}D_{t}\boldsymbol{U}+\mathbf{1}_{z}\times \boldsymbol{U}=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F1}+\boldsymbol{F},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{U}=0,\quad \text{where }\boldsymbol{F}=\boldsymbol{J}\times \boldsymbol{B} & & \displaystyle\end{eqnarray}$$

is the Lorentz force, which provides the power that maintains $\boldsymbol{U}(\boldsymbol{x},t)$ and, when there is dynamo action, maintains $\boldsymbol{B}(\boldsymbol{x},t)$ too; in (3.3a ), $D_{t}\;(=\unicode[STIX]{x2202}_{t}+\boldsymbol{U}\boldsymbol{\cdot })$ is the material derivative. For simplicity, the possibility that $\boldsymbol{F}$ contains a radial buoyancy force has been excluded in paper 2 and in this paper (§ 2). The $\unicode[STIX]{x1D6FC}$ -effect alone powers our mean-field $\unicode[STIX]{x1D6FC}^{2}$ dynamos, but our $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D714}$ models are powered too by the rate of working of the imposed shearing velocity against the Lorentz force; see below.

Mass conservation (3.3b ) and the boundary conditions (2.1a,b ) imply (Taylor Reference Taylor1963)

(3.3d ) $$\begin{eqnarray}\displaystyle U_{s}^{\text{G}}(s,t)=0\quad \text{for all }0\leqslant s\leqslant 1. & & \displaystyle\end{eqnarray}$$

Apart from the pressure gradient, only the Lorentz and Coriolis forces affect the motion of fluid in the volume ${\mathcal{V}}(s)$ defined in § 2. Its equation of motion is, by (3.3a ),

(3.4a ) $$\begin{eqnarray}\displaystyle E_{\unicode[STIX]{x1D702}}[\unicode[STIX]{x2202}_{t}m_{z}(s,t)+f_{s}(s,t)]=\unicode[STIX]{x1D6E4}_{z}(s,t), & & \displaystyle\end{eqnarray}$$

where $m_{z}(s,t)$ is the $z$ -component of the angular momentum of ${\mathcal{V}}(s)$ ,

(3.4b ) $$\begin{eqnarray}\displaystyle m_{z}(s,t)\equiv \mathbf{1}_{z}\boldsymbol{\cdot }\int _{{\mathcal{V}}(s)}\boldsymbol{x}\times \boldsymbol{U}\,\text{d}^{3}x=\int _{0}^{s}s\left[\int _{{\mathcal{C}}(s)}U_{\unicode[STIX]{x1D719}}(\boldsymbol{x},t)\,\text{d}^{2}x\right]\text{d}s=\int _{0}^{s}s{\mathcal{A}}(s)U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)\,\text{d}s, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

$f_{s}(s,t)$ is the flux of angular momentum from ${\mathcal{V}}(s)$ across ${\mathcal{C}}(s)$ ,

(3.4c ) $$\begin{eqnarray}\displaystyle f_{s}(s,t)=\int _{{\mathcal{C}}(s)}sU_{s}(\boldsymbol{x},t)U_{\unicode[STIX]{x1D719}}(\boldsymbol{x},t)\,\text{d}^{2}x=s{\mathcal{A}}(s)(U_{s}U_{\unicode[STIX]{x1D719}})^{\text{G}}(s,t), & & \displaystyle\end{eqnarray}$$

and $\unicode[STIX]{x1D6E4}_{z}(s,t)$ is the $z$ -component of the torque on ${\mathcal{V}}(s)$ ,

(3.4d ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E4}_{z}(s,t)\equiv \mathbf{1}_{z}\boldsymbol{\cdot }\int _{{\mathcal{V}}(s)}\boldsymbol{x}\times \boldsymbol{F}\,\text{d}^{3}x=\int _{0}^{s}s\left[\int _{{\mathcal{C}}(s)}F_{\unicode[STIX]{x1D719}}(\boldsymbol{x},t)\,\text{d}^{2}x\right]\text{d}s=\int _{0}^{s}s{\mathcal{A}}(s)F_{\unicode[STIX]{x1D719}}^{G}(s,t)\,\text{d}s. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

By writing the Lorentz force (3.3c ) as the divergence of the magnetic stress tensor, and applying the dynamo condition (3.1h ), $\unicode[STIX]{x1D6E4}_{z}(s,t)$ can be expressed as a surface integral:

(3.5a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E4}_{z}(s,t)=E_{\unicode[STIX]{x1D702}}^{-1}\int _{{\mathcal{C}}(s)}sB_{s}B_{\unicode[STIX]{x1D719}}\,\text{d}^{2}x. & & \displaystyle\end{eqnarray}$$

Differentiating (3.5a ) by $t$ , applying (3.2a,d ) and adding a star to Taylor’s $\unicode[STIX]{x1D6FC}$ to avoid confusion with the $\unicode[STIX]{x1D6FC}$ -effect, we find

(3.5b,c ) $$\begin{eqnarray}\displaystyle E_{\unicode[STIX]{x1D702}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}_{z}(s,t)}{\unicode[STIX]{x2202}t}=s[s\unicode[STIX]{x1D6FC}^{\ast }Z(s,t)-S(s,t)],\quad \text{where }\unicode[STIX]{x1D6FC}^{\ast }(s,t)=\int _{{\mathcal{C}}(s)}B_{s}^{2}\,\text{d}^{2}x\quad \text{and}\qquad & & \displaystyle\end{eqnarray}$$

(3.5d ) $$\begin{eqnarray}\displaystyle S(s,t)=-\int _{{\mathcal{C}}(s)}(B_{s}C_{\unicode[STIX]{x1D719}}^{\text{A}}+B_{\unicode[STIX]{x1D719}}C_{s}^{\text{A}})\,\text{d}^{2}x. & & \displaystyle\end{eqnarray}$$

Equation (3.5b ) determines $Z(s,t)$ from $\unicode[STIX]{x1D6E4}_{z}(s,t)$ , $S(s,t)$ and $\unicode[STIX]{x1D6FC}^{\ast }(s,t)$ . For dipolar $\boldsymbol{B}$ , these are each $O(s^{4})$ as $s\rightarrow 0$ ; for quadrupolar $\boldsymbol{B}$ , they are each $O(s^{2})$ as $s\rightarrow 0$ . In both these cases, $Z(s,t)=O(s^{-1})$ , $\unicode[STIX]{x1D701}(s,t)=O(\log s)$ and $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)=O(s\log s)$ for $s\rightarrow 0$ . This weak singularity of $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)$ may be an endemic shortcoming of axisymmetric MS models and of the axisymmetric part of more general, three-dimensional MS models.

By (3.1h ), (3.4d ), (3.5a ), and $B_{\unicode[STIX]{x1D719}}(1,t)=0$ , by

(3.5e-g ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E4}_{z}(1,t)(=\unicode[STIX]{x1D6E4}_{z}^{\text{total}})=0,\quad \unicode[STIX]{x2202}_{s}\unicode[STIX]{x1D6E4}_{z}(1,t)=0,\quad F_{\unicode[STIX]{x1D719}}^{\text{G}}(1,t)=0. & & \displaystyle\end{eqnarray}$$

Finally, to show that $\unicode[STIX]{x1D6FC}$ is the energy source driving the $\unicode[STIX]{x1D6FC}^{2}$ system, we develop the energy equation. The total kinetic and magnetic energies are

(3.6a,b ) $$\begin{eqnarray}\displaystyle {\mathcal{K}}=\frac{1}{2}\int _{\text{v}}\boldsymbol{U}^{2}\,\text{d}^{3}x,\quad {\mathcal{M}}=\frac{1}{2}E_{\unicode[STIX]{x1D702}}^{-1}\int _{\text{v}_{\infty }}\boldsymbol{B}^{2}\,\text{d}^{3}x, & & \displaystyle\end{eqnarray}$$

where $\text{v}_{\infty }=\text{v}+\text{v}^{\ast }$ is all space. By taking the scalar product of (3.3a ) with $\boldsymbol{U}$ and the scalar product of (3.1e ) with $\boldsymbol{B}$ , integrating over v and applying the associated boundary conditions, it is found that

(3.6c,d ) $$\begin{eqnarray}\displaystyle E_{\unicode[STIX]{x1D702}}\frac{\unicode[STIX]{x2202}{\mathcal{K}}}{\unicode[STIX]{x2202}t}={\mathcal{L}},\quad \frac{\unicode[STIX]{x2202}{\mathcal{M}}}{\unicode[STIX]{x2202}t}={\mathcal{P}}-{\mathcal{L}}-{\mathcal{Q}}, & & \displaystyle\end{eqnarray}$$

where ${\mathcal{L}}$ is the rate of working of the Lorentz force, ${\mathcal{P}}$ is the power input and ${\mathcal{Q}}(t)$ is the dissipation, which is entirely ohmic in the MS approximation:

(3.6e-g ) $$\begin{eqnarray}\displaystyle {\mathcal{L}}=\int _{\text{v}}\boldsymbol{U}\boldsymbol{\cdot }(\boldsymbol{J}\times \boldsymbol{B})\,\text{d}^{3}x,\quad {\mathcal{P}}=\int _{\text{v}}\unicode[STIX]{x1D6FC}\boldsymbol{J}\boldsymbol{\cdot }\boldsymbol{B}\,\text{d}^{3}x,\quad {\mathcal{Q}}=\int _{\text{v}}J^{2}\,\text{d}^{3}x. & & \displaystyle\end{eqnarray}$$

Clearly $\unicode[STIX]{x2202}_{t}(E_{\unicode[STIX]{x1D702}}{\mathcal{K}}+{\mathcal{M}})={\mathcal{P}}-{\mathcal{Q}}$ . As stated above, the power input, ${\mathcal{P}}$ , is created by $\unicode[STIX]{x1D6FC}$ . Because the integrand of (3.6f ) is the product of two helicities, the kinetic helicity $\unicode[STIX]{x1D6FC}$ and ‘the current helicity’ $\boldsymbol{J}\boldsymbol{\cdot }\boldsymbol{B}$ , which is a pseudoscalar too, ${\mathcal{P}}$ is correctly independent of the right- or left-handedness of the coordinate frame.

In a steadily operating $\unicode[STIX]{x1D6FC}^{2}$ dynamo, (3.6c,d ) show that $\langle {\mathcal{L}}\rangle =0$ and $\langle {\mathcal{P}}\rangle =\langle {\mathcal{Q}}\rangle$ , where $\langle Q\rangle$ is the average of $Q(t)$ over a sufficiently long time, e.g. the magnetic diffusion time, which is $O(1)$ in our non-dimensionalization. The positivity of ${\mathcal{P}}$ implies that, in a steadily working $\unicode[STIX]{x1D6FC}^{2}$ dynamo, the kinetic helicity on average creates a current helicity of the same sign.

The kinetic energy ${\mathcal{K}}$ in (3.6a ) may be written as ${\mathcal{K}}=(1/2)\int _{\text{v}}[U_{\text{M}}^{2}+s^{2}\unicode[STIX]{x1D701}^{2}]\,\text{d}^{3}x$ where $\boldsymbol{U}_{\text{M}}=U_{s}\mathbf{1}_{s}+U_{z}\mathbf{1}_{z}$ is the meridional part of $\boldsymbol{U}$ . For an $\unicode[STIX]{x1D6FC}^{2}\unicode[STIX]{x1D714}$ dynamo, defined in § 1, this is replaced by ${\mathcal{K}}=(1/2)\int _{\text{v}}[U_{\text{M}}^{2}+s^{2}(\unicode[STIX]{x1D701}+\unicode[STIX]{x1D714})^{2}]\,\text{d}^{3}x$ which, in the extreme case of an $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D714}$ dynamo $(|\unicode[STIX]{x1D6FC}|\ll r_{o}|\unicode[STIX]{x1D714}|)$ , is ${\mathcal{K}}=(1/2)\int _{\text{v}}[U_{\text{M}}^{2}+s^{2}\unicode[STIX]{x1D714}^{2}]\,\text{d}^{3}x$ .

3.1 Original Taylor theory (OTT); Taylor’s constraint

This subsection summarizes what is needed later in this paper about OTT, our abbreviation for ‘original Taylor theory’ (Taylor Reference Taylor1963).

The study of magnetostrophy was initiated by the plasma physicist Bryan Taylor, whose paper had a major impact on the dynamo community. As well as neglecting viscosity, Taylor (Reference Taylor1963) discarded the inertial force in the rotating frame. Inertia is generally quantified by the Rossby number $Ro\equiv {\mathcal{U}}/2\unicode[STIX]{x1D6FA}r_{o}$ . For an estimated flow speed ${\mathcal{U}}=2\times 10^{-4}~\text{m}~\text{s}^{-1}$ , $Ro\approx 4\times 10^{-7}$ . Our choice of non-dimensional variables makes it more convenient to use instead the magnetic Ekman number (1.5). For $\unicode[STIX]{x1D702}\approx 0.7~\text{m}^{2}~\text{s}^{-1}$ (Pozzo et al. Reference Pozzo, Davies, Gubbins and Alfé2012), $E_{\unicode[STIX]{x1D702}}\approx 4\times 10^{-10}$ . The ratio of these 2 inertial measures is the magnetic Reynolds number: $Rm\equiv {\mathcal{U}}r_{o}/\unicode[STIX]{x1D702}=Ro/E_{\unicode[STIX]{x1D702}}$ , which is approximately 1000. The omission of the inertial force may therefore be stated either as $Ro=0$ or as $E_{\unicode[STIX]{x1D702}}=0$ . In writing $E_{\unicode[STIX]{x1D702}}=0$ , there is no implication that $\unicode[STIX]{x1D702}$ is as negligible as $\unicode[STIX]{x1D708}$ ; $E_{\unicode[STIX]{x1D702}}$ is the ratio of the dynamo time scale, $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}=r_{o}^{2}/\unicode[STIX]{x1D702}$ , and the other time scales governing geodynamo (identified in § 4.1 below) and is much greater than any of them; $E_{\unicode[STIX]{x1D702}}=0$ is an expression of that fact.

The smallness of $Ro$ encouraged Taylor (Reference Taylor1963) to ignore inertia entirely in the rotating frame. For $E_{\unicode[STIX]{x1D702}}=0$ , the Boussinesq MS equations (3.3a,b ) governing $\boldsymbol{U}$ become

(3.7a,b ) $$\begin{eqnarray}\displaystyle \mathbf{1}_{z}\times \boldsymbol{U}=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F1}+\boldsymbol{F},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{U}=0. & & \displaystyle\end{eqnarray}$$

The derivation of energy conservation given in § 3 remains valid for ${\mathcal{K}}=0$ . The equation of motion (3.4a ) of ${\mathcal{V}}(s)$ and its $s$ -differential are also simplified, becoming

(3.7c,d ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E4}_{z}(s,t)\equiv 0,\quad F_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)=0,\quad \text{for all }0\leqslant s\leqslant 1\text{ and all }t. & & \displaystyle\end{eqnarray}$$

Equivalent to (3.7d ) is

(3.7e ) $$\begin{eqnarray}\displaystyle \int _{{\mathcal{C}}(s)}(\boldsymbol{J}\times \boldsymbol{B})_{\unicode[STIX]{x1D719}}\,\text{d}^{2}x=0,\quad \text{for all }0\leqslant s\leqslant 1\text{ and all }t. & & \displaystyle\end{eqnarray}$$

Either of (3.7d,e ), but usually (3.7e ), is termed ‘Taylor’s constraint’ (aka ‘Taylor’s condition’). Provided it is satisfied, the solution of (3.7a,b ) can almost be reduced to quadratures, as Taylor (Reference Taylor1963) demonstrated.

To maintain his constraint, Taylor (Reference Taylor1963, equation (4.5)) derived a second-order equation for $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)$ . The first integral of that equation is

(3.8a,b ) $$\begin{eqnarray}\displaystyle Z(s,t)=\frac{q(s,t)}{s},\quad \text{where }q(s,t)=\frac{S(s,t)}{\unicode[STIX]{x1D6FC}^{\ast }(s,t)}. & & \displaystyle\end{eqnarray}$$

The function $q(s,t)$ defined by (3.8b ) is regular for $0\leqslant s\leqslant 1$ and non-zero at $s=0$ , as shown below (3.5d ). This means that $\unicode[STIX]{x1D701}(s,t)$ cannot be obtained by integrating $Z(s,t)$ directly because of the logarithmic divergence found at $s=0$ in § 3. Instead we introduce $\unicode[STIX]{x1D716}\;({>}0)$ and define $\unicode[STIX]{x1D701}^{\unicode[STIX]{x1D716}}(s,t)$ by excluding $s=0$ from the range of integration, using (3.8a ) to obtain

(3.8c,d ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D701}^{\unicode[STIX]{x1D716}}(s,t)=\int _{\unicode[STIX]{x1D716}}^{s}\frac{S(s,t)}{s\unicode[STIX]{x1D6FC}^{\ast }(s,t)}\,\text{d}s+D_{0}(t),~~\text{so that }\unicode[STIX]{x1D701}^{\unicode[STIX]{x1D716}}(s,t)\sim q(0,t)\ln \!\left(\frac{s}{\unicode[STIX]{x1D716}}\right)+D_{0}(t)\text{ for }s\rightarrow 0, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $D_{0}(t)$ is the integration ‘constant’. Clearly $\unicode[STIX]{x1D701}^{\unicode[STIX]{x1D716}}(s,t)$ is finite, but diverges as $s\rightarrow 0$ or $\unicode[STIX]{x1D716}\rightarrow 0$ . We approximated $\unicode[STIX]{x1D701}(s,t)$ numerically by

(3.8e,f ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D701}(s,t)\approx \unicode[STIX]{x1D701}^{\unicode[STIX]{x1D716}}(s,t),\quad \text{where }\unicode[STIX]{x1D716}=s_{1} & & \displaystyle\end{eqnarray}$$

is the smallest positive grid point (at $s=0$ , we arbitrarily take $\unicode[STIX]{x1D701}^{\unicode[STIX]{x1D716}}(0,t)=0$ ). Equation (2.2) shows that

(3.8g ) $$\begin{eqnarray}\displaystyle D_{0}(t)=-\frac{15}{4\unicode[STIX]{x03C0}}\int _{0}^{1}s^{2}z_{1}\unicode[STIX]{x1D701}(s,t)\,\text{d}^{2}x. & & \displaystyle\end{eqnarray}$$

(Even though this result was derived from (3.8c ), it features $\unicode[STIX]{x1D701}$ and not $\unicode[STIX]{x1D701}^{\unicode[STIX]{x1D716}}$ . This is because $\lim _{\unicode[STIX]{x1D716}\rightarrow 0}m_{z}^{\unicode[STIX]{x1D716},\text{total}}=m_{z}^{\text{total}}$ .)

Most modellers, including Taylor (Reference Taylor1963) and ourselves in this paper, assume the core is fluid throughout, ignoring the inner core entirely. Livermore, Ierley & Jackson (Reference Livermore, Ierley and Jackson2009) provide general examples of full-sphere Taylor states. Recognition of the inner core introduces complications touched on in § 2 above. Nevertheless, the behaviour of solutions as $s\rightarrow 0$ is not much affected; the logarithmic divergence of $\unicode[STIX]{x1D701}$ persists as does the endemic shortcoming of $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)$ identified in § 3.

In the present axisymmetric case, the cylindrical components of (3.7a ) are

(3.9a-c ) $$\begin{eqnarray}\displaystyle U_{s}(\boldsymbol{x},t)=F_{\unicode[STIX]{x1D719}}(\boldsymbol{x},t),\quad U_{\unicode[STIX]{x1D719}}(\boldsymbol{x},t)=-F_{s}(\boldsymbol{x},t)+\unicode[STIX]{x2202}_{s}\unicode[STIX]{x1D6F1}(\boldsymbol{x},t),\quad U_{z}(\boldsymbol{x},t)=-s^{-1}\unicode[STIX]{x2202}_{s}(sz_{1}\widetilde{F}_{\unicode[STIX]{x1D719}}), & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where

(3.9d,e ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F1}(\boldsymbol{x},t)=z_{1}\widetilde{F}_{z}(\boldsymbol{x},t)+w(s,t),\quad \frac{\unicode[STIX]{x2202}w(s,t)}{\unicode[STIX]{x2202}s}=F_{s}^{\text{G}}(s,t)+\frac{1}{z_{1}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}(s[zF_{z}]^{\text{G}}(s,t)),\quad & & \displaystyle\end{eqnarray}$$

(3.9f ) $$\begin{eqnarray}\displaystyle \widetilde{Q}(\boldsymbol{x},t)=\frac{1}{2z_{1}}\left(\int _{z_{1}}^{z}Q(s,z,t)+\int _{-z_{1}}^{z}Q(s,z,t)\right). & & \displaystyle\end{eqnarray}$$

Solutions of these equations were given in (1: § 4), where we recognized that the initial $\boldsymbol{B}(\boldsymbol{x},0)$ , given for example by (2.1a ), determines the initial $\boldsymbol{U}_{\unicode[STIX]{x1D719}}^{\text{G}}(s,0)$ .

Finding OTT solutions promises to be much easier than solving (3.1e ) and the unapproximated partial differential equations governing $\boldsymbol{U}$ that include viscous and inertial forces. This explains why Taylor’s paper made a major impact on the dynamo community with its publication in 1963. Nevertheless, the task proved to be more onerous than it appeared to be at first sight, and it took 52 years for Taylor’s idea to be successfully implemented. See (2: §§ 5,6), which presents solutions for mean-field $\unicode[STIX]{x1D6FC}^{2}$ and $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D714}$ -dynamos.

4 Modified Taylor theory (MTT); torsional waves

In this section, a generalization of OTT will be described which, as in paper 1, is termed ‘modified Taylor theory’ (abbreviated as ‘MTT’). Numerical consequences are presented in § 5.3 below. MTT recognizes the existence of torsional waves which OTT does not. Much has been written about torsional waves since their importance was first recognized by Braginsky (Reference Braginsky1970). See, for example, Wicht & Christensen (Reference Wicht and Christensen2010), Teed, Jones & Tobias (Reference Teed, Jones and Tobias2013, Reference Teed, Jones and Tobias2015), Schaeffer et al. (Reference Schaeffer, Jault, Nataf and Fournier2017). Interpretations of observations that appear to confirm that torsional waves are travelling across the core are due to Abarca del Rio, Gambis & Salstein (Reference Abarca del Rio, Gambis and Salstein2000), Gillet et al. (Reference Gillet, Jault, Canet and Fournier2010) and Holme & De Viron (Reference Holme and De Viron2013), who link torsional oscillations to observed changes in the length of day. But we should emphasize here, however, that the incorporation of torsional waves into our MS dynamo model is not our main objective. We have studied MTT primarily because it not only generalizes OTT but also because it is, as we see it, the best MS approximation possible, without re-introducing the shortcomings inherent in today’s viscously dominated models. We defend that assertion in § 4.1 below.

MTT ignores the nonlinearity of the inertial force; the $f_{s}$ defined by (3.4c ) is now zero, so that (3.4a,b,d ) and (3.3d ) imply

(4.1) $$\begin{eqnarray}\displaystyle E_{\unicode[STIX]{x1D702}}\frac{\unicode[STIX]{x2202}U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)}{\unicode[STIX]{x2202}t}=\frac{E_{\unicode[STIX]{x1D702}}}{s{\mathcal{A}}(s)}\frac{\unicode[STIX]{x2202}^{2}m_{z}(s,t)}{\unicode[STIX]{x2202}s\unicode[STIX]{x2202}t}=\frac{1}{s{\mathcal{A}}(s)}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}_{z}(s,t)}{\unicode[STIX]{x2202}s}=F_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t). & & \displaystyle\end{eqnarray}$$

The left-hand side of (4.1) is the only part of the inertial force $E_{\unicode[STIX]{x1D702}}(\unicode[STIX]{x2202}_{t}\boldsymbol{U}+\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{U})$ retained in MTT; the momentum equations (3.3a–c ) and (3.7a–c ) are replaced by

(4.2a-c ) $$\begin{eqnarray}\displaystyle E_{\unicode[STIX]{x1D702}}\unicode[STIX]{x2202}_{t}U_{\unicode[STIX]{x1D719}}^{\text{G}}\mathbf{1}_{\unicode[STIX]{x1D719}}+\mathbf{1}_{z}\times \boldsymbol{U}=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F1}+\boldsymbol{F},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{U}=0,\quad \boldsymbol{F}=\boldsymbol{J}\times \boldsymbol{B}. & & \displaystyle\end{eqnarray}$$

The derivation of energy conservation given in § 3 continues to be valid provided the kinetic energy ${\mathcal{K}}$ in (3.6a ) is redefined with $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)$ replacing $\boldsymbol{U}(\boldsymbol{x},t)$ .

Although (3.3d ) remains true ( $U_{s}^{\text{G}}(s,t)=0$ ), (4.2a–c ) no longer imply Taylor’s constraint (3.7d,e ). Because the predictive (4.1) replaces the diagnostic (3.7a ), it does not constrain $\boldsymbol{F}$ . Instead

(4.2d ) $$\begin{eqnarray}\displaystyle E_{\unicode[STIX]{x1D702}}\frac{\unicode[STIX]{x2202}U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)}{\unicode[STIX]{x2202}t}=\frac{1}{{\mathcal{A}}(s)}\int _{{\mathcal{C}}(s)}(\boldsymbol{J}\times \boldsymbol{B})_{\unicode[STIX]{x1D719}}\,\text{d}^{2}x,\quad \text{for all }0\leqslant s\leqslant 1. & & \displaystyle\end{eqnarray}$$

When integrating (4.2d ) from an initial state such as (2.3a,b ), the resulting $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)$ evidently depends on $E_{\unicode[STIX]{x1D702}}$ , and (4.1) shows that

(4.3a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}(E_{\unicode[STIX]{x1D702}}s{\mathcal{A}}(s)U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t))=\unicode[STIX]{x2202}_{s}\unicode[STIX]{x1D6E4}_{z}. & & \displaystyle\end{eqnarray}$$

By writing (4.2d ) in this conservative form, the condition (2.2) on the total angular momentum ( $m_{z}^{\text{total}}(t)=0$ ), is automatically satisfied as $\boldsymbol{U}(\boldsymbol{x},t)$ evolves (provided that $U_{\unicode[STIX]{x1D719}}^{\text{G}}(1,0)=0$ so that $m_{z}^{\text{total}}$ is initially zero). Equations (3.5e ) and (4.2d ) imply

(4.3b,c ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}U_{\unicode[STIX]{x1D719}}^{\text{G}}(1,t)=0,\quad \text{so that }U_{\unicode[STIX]{x1D719}}^{\text{G}}(1,t)=U_{\unicode[STIX]{x1D719}}^{\text{G}}(1,0)=\text{const.}=K,\text{ say}. & & \displaystyle\end{eqnarray}$$

In MTT, as for OTT, $\unicode[STIX]{x1D701}(s,t)$ is $O(\ln s)$ for $s\rightarrow 0$ . See figure 1, which was derived by integrating (4.3a ). Except for the first half-grid point, the computed values of $\unicode[STIX]{x1D701}$ fit well onto the curve $\unicode[STIX]{x1D701}(s,t_{0})=-20.5\ln s-50$ shown, where $t_{0}=54.8$ is the time at which the data were taken.

Figure 1. The function $\unicode[STIX]{x1D701}(s,t)$ at $t=54.8$ as a function of $s$ for the modified Taylor theory $\unicode[STIX]{x1D6FC}^{2}$ model with $\unicode[STIX]{x1D6FC}_{0}=14.5$ , $E_{\unicode[STIX]{x1D702}}=0.01$ and $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,0)=0$ is shown with blue curves marked by circles: (a) in the range $0<s\leqslant 1$ , and (b) in $0<s\leqslant 0.2$ . (The function $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)$ at the same time is shown in figure 3 below.) Also shown in red in (b) is the curve $\unicode[STIX]{x1D701}(s,t)=-20.5\ln (s)-50$ .

In OTT, (3.4a ) is the tautology 0 = 0; its left-hand side vanishes because $E_{\unicode[STIX]{x1D702}}=0$ and its right-hand side because of Taylor’s constraint (3.7d,e ). An important part of Taylor’s paper (1963, § 4) is devoted to finding his $u_{\unicode[STIX]{x1D719}}(s,t)$ , i.e. our $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)$ . We used his method in deriving OTT results in paper 2 but cannot apply it here because $F_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)$ is generally non-zero in MTT.

In MTT, $\unicode[STIX]{x1D6E4}_{z}(s,t)$ and $m_{z}(s,t)$ are generally non-zero and are related by (3.4a ) with $f_{s}(s,t)=0$ . They are responsible for coupling geostrophic cylinders together and for time variation of the angular momentum, both of which are associated with torsional waves and both of which are lacking in OTT. An equation governing $\unicode[STIX]{x1D701}(s,t)(=U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)/s)$ alone can by found by differentiating (3.5b ) by $s$ and then applying (4.2d ), to give

(4.4a ) $$\begin{eqnarray}\displaystyle E_{\unicode[STIX]{x1D702}}s^{2}{\mathcal{A}}(s)\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D701}}{\unicode[STIX]{x2202}t^{2}}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left[s\left(s\unicode[STIX]{x1D6FC}^{\ast }(s,t)\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}(s,t)}{\unicode[STIX]{x2202}s}-S(s,t)\right)\right]. & & \displaystyle\end{eqnarray}$$

Equation (4.4a ) governs a one-dimensional type of Alfvén wave that travels in the $\pm s$ -direction. This torsional wave equation may be re-written as (see (1: equation (22a)))

(4.4b ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D701}}{\unicode[STIX]{x2202}t^{2}}=\frac{1}{s^{2}{\mathcal{A}}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left(s^{2}{\mathcal{A}}{\mathcal{V}}_{\text{tort}}^{2}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}{\unicode[STIX]{x2202}s}\right)-\frac{1}{s^{2}{\mathcal{A}}E_{\unicode[STIX]{x1D702}}}\frac{\unicode[STIX]{x2202}(sS)}{\unicode[STIX]{x2202}s}, & & \displaystyle\end{eqnarray}$$

where ${\mathcal{V}}_{\text{tort}}(s,t)$ is the torsional wave speed (Braginsky Reference Braginsky1970):

(4.4c,d ) $$\begin{eqnarray}\displaystyle {\mathcal{V}}_{\text{tort}}(s,t)=\left(\frac{1}{E_{\unicode[STIX]{x1D702}}{\mathcal{A}}(s)}\int _{{\mathcal{C}}(s)}B_{s}^{2}(\boldsymbol{x},t)\,\text{d}^{2}x\right)^{1/2}\quad \text{or}\quad {\mathcal{V}}_{\text{tort}}(s,t)=E_{\unicode[STIX]{x1D702}}^{-1/2}\left(\frac{\unicode[STIX]{x1D6FC}^{\ast }(s,t)}{{\mathcal{A}}(s)}\right)^{1/2}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

It is generally true that $\unicode[STIX]{x1D6FC}^{\ast }(s,t)\propto z_{1}^{3}$ as $s\rightarrow 1$ and, since ${\mathcal{A}}(s)\propto z_{1}$ , ${\mathcal{V}}_{\text{tort}}(s,t)\propto (1-s)^{1/2}$ for $s\rightarrow 1$ by (4.4d ). For dipole parity, ${\mathcal{V}}_{\text{tort}}(s,t)\propto s$ as $s\rightarrow 0$ . A good approximation to the torsional wave speed is ${\mathcal{V}}_{\text{tort}}(s,t)={\mathcal{V}}_{0}(t)s(1-s)^{1/2}$ , where ${\mathcal{V}}_{0}(0)=(3/E_{\unicode[STIX]{x1D702}})^{1/2}B_{0}$ in the special case (2.3a ). For magnetic fields of quadrupolar parity, $\unicode[STIX]{x1D6FC}^{\ast }(s,t)\propto s$ for $s\rightarrow 0$ so that ${\mathcal{V}}_{\text{tort}}(0,t)=O(1)$ .

As $E_{\unicode[STIX]{x1D702}}$ is so small (of order $4\times 10^{-10}$ , according to § 3.1) and ${\mathcal{V}}_{\text{tort}}(s,t)=O(E_{\unicode[STIX]{x1D702}}^{-1/2}{\mathcal{B}})$ , torsional waves introduce a new short time scale, $\unicode[STIX]{x1D70F}_{\text{tort}}={\mathcal{V}}_{\text{tort}}^{-1}$ , into core dynamics. In dimensional units for Earth’s core, ${\mathcal{V}}_{\text{tort}}(s,t)\approx 2\times 10^{-2}~\text{m}~\text{s}^{-1}$ and is about $100{\mathcal{U}}$ ; for ${\mathcal{B}}=25$ gauss, $\unicode[STIX]{x1D70F}_{\text{tort}}\approx 5~\text{yr}$ so that $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}/\unicode[STIX]{x1D70F}_{\text{tort}}\approx 10^{5}$ , assuming the dynamo time scale $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}=r_{o}^{2}/\unicode[STIX]{x1D702}\approx 5\times 10^{5}~\text{yr}$ . So there are many torsional oscillations of period $\unicode[STIX]{x1D70F}_{\text{tort}}$ during the transition to the steady state of our $\unicode[STIX]{x1D6FC}^{2}$ model, but the duration of that transition is $O(\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}})$ , this being the time scale over which ohmic resistance damps out the waves. All our numerical integrations were carried out on a workstation; this made values of $E_{\unicode[STIX]{x1D702}}$ smaller than $10^{-4}$ inaccessible to us. Figure 2 below therefore shows only one oscillation convincingly because $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}/\unicode[STIX]{x1D70F}_{\text{tort}}\approx E_{\unicode[STIX]{x1D702}}^{-1/2}{\mathcal{B}}$ , is not large enough to show many. For the $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D714}$ models of paper 2, the final OTT solutions were periodic, with periods that were unrelated to $\unicode[STIX]{x1D70F}_{\text{tort}}$ , but the MTT evolution to these solutions is accompanied by torsional waves.

Figure 2. Internal magnetic energy, $M(t)$ , as a function of time for the $\unicode[STIX]{x1D6FC}^{2}$ model for $\unicode[STIX]{x1D6FC}_{0}=14.5$ . The solid curve is for MTT with $U_{\unicode[STIX]{x1D719}}^{\text{G}}(1,0)=0$ and $E_{\unicode[STIX]{x1D702}}=0.01$ ; the dashed curve is for $E_{\unicode[STIX]{x1D702}}=0.001$ . The dotted curve is from an OTT calculation with the same $\unicode[STIX]{x1D6FC}_{0}$ . During the linear phase of the growth, $M(t)\propto \exp (2\unicode[STIX]{x1D706}t)$ , where $\unicode[STIX]{x1D706}=2.55$ for OTT; $\unicode[STIX]{x1D706}=5.50$ for both MTT calculations.

The significance of the short time scale of torsional waves was not lost on Taylor (Reference Taylor1963, § 5), and it seems to us worth quoting a short extract from his paper where he recognizes the limitations of OTT that are now removed in MTT. We make minor changes to his text so that it conforms with our notation:

There is, however, a perplexing paradox. It appears, from the numerical results of § 5.3 for the $\unicode[STIX]{x1D6FC}^{2}$ model, that there is only one value of the constant $K$ defined in (4.3c ) for which $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,\infty )$ for MTT coincides for all $s$ with $[U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,\infty )]_{\text{O}\text{T}\text{T}}$ for OTT. The single exception is $K=[U_{\unicode[STIX]{x1D719}}^{\text{G}}(1,\infty )]_{\text{O}\text{T}\text{T}}$ for which the MTT and OTT solutions are the same at infinite $t$ . The point is simple: except for that special value of $K$ , one cannot require the MTT flow to satisfy $U_{\unicode[STIX]{x1D719}}^{\text{G}}(1,\infty )=K$ and $U_{\unicode[STIX]{x1D719}}^{\text{G}}(1,\infty )=[U_{\unicode[STIX]{x1D719}}^{\text{G}}(1,\infty )]_{\text{O}\text{T}\text{T}}$ simultaneously! Section 5.3 below points the way MTT solutions solve the paradox.

At long last, we return from finding $\boldsymbol{U}^{\text{G}}(\boldsymbol{x},t)\;(=U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)\mathbf{1}_{\unicode[STIX]{x1D719}})$ to determining $\boldsymbol{U}^{\text{A}}(\boldsymbol{x},t)$ . It is helpful to cast the defining (4.2a ) of MTT in the same form as (3.7a ) for OTT by replacing $\boldsymbol{F}(\boldsymbol{x},t)$ in (3.7a ) by

(4.5) $$\begin{eqnarray}\displaystyle \boldsymbol{H}=\boldsymbol{F}-E_{\unicode[STIX]{x1D702}}\unicode[STIX]{x2202}_{t}U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)\mathbf{1}_{\unicode[STIX]{x1D719}}. & & \displaystyle\end{eqnarray}$$

Note particularly that $H_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)=0$ , which is analogous to $F_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)=0$ for OTT; see (3.7d ). This means that the process of solving (4.2a,b ) is identical to that of solving (3.7a,b ), which was achieved by (3.9a–e ). Both processes determine the flow $\boldsymbol{U}^{\text{A}}(\boldsymbol{x},t)$ at the same time $t$ as they advance $\boldsymbol{B}(\boldsymbol{x},t)$ (in the case of OTT) or $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)$ (in the case of MTT). This point is elaborated in the next subsection.

5 Numerical results

5.1 Task for MTT

As already mentioned in § 2, we seek only MTT mean-field dynamos in this paper. Although this rescues the torsional waves abandoned by OTT, the inclusion of these waves by MTT is only incidental to our main objective. As we stressed in the last subsection, MTT is the best way of approximating the inertial force while avoiding the large Ekman numbers that are so devastating to the geophysical relevance of today’s computer models; see § 2. We now present numerical results for MTT dynamos; these may be compared with the OTT dynamos of paper 2.

We start by summarizing §§ 3 and 4. In one respect MTT solutions are simpler to obtain than OTT solution. There are no constraints and $\boldsymbol{U}$ is obtained by time stepping predictive equations. But one disadvantage of MTT is that one has to take more time steps in integrating the governing equations because the torsional wave time scale $\unicode[STIX]{x1D70F}_{\text{tort}}=1/V_{\text{tort}}\propto E_{\unicode[STIX]{x1D702}}^{1/2}$ has to be resolved, and when $E_{\unicode[STIX]{x1D702}}$ is small this is much less than the dynamo time scale $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}=r_{o}^{2}/\unicode[STIX]{x1D702}$ .

MTT requires the numerical solution of 4 predictive equations, namely the 3 scalar equations (3.1e ) and the modified Taylor condition (4.2a ):

(5.1a,b ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{B}}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D735}\times (\boldsymbol{U}\times \boldsymbol{B}+\unicode[STIX]{x1D6FC}\boldsymbol{B})+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{B},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}=0, & & \displaystyle\end{eqnarray}$$

(5.1c ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)}{\unicode[STIX]{x2202}t}=\frac{1}{E_{\unicode[STIX]{x1D702}}{\mathcal{A}}(s)}\int _{{\mathcal{C}}(s)}(\boldsymbol{J}\times \boldsymbol{B})_{\unicode[STIX]{x1D719}}\,\text{d}^{2}x. & & \displaystyle\end{eqnarray}$$

Equations (5.1a,b ) and the dynamo conditions (3.1h,i ) decide the evolution of $\boldsymbol{B}$ ; equation (5.1c ) and condition (2.2) determine how $U_{\unicode[STIX]{x1D719}}^{\text{G}}$ evolves. At each time step, $\boldsymbol{U}(\boldsymbol{x},t)$ must be found from $\boldsymbol{F}(\boldsymbol{x},t)$ and $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)$ by the process described in § 4.

We integrated (5.1a,c ), starting from

(5.2a,b ) $$\begin{eqnarray}\displaystyle \boldsymbol{B}(\boldsymbol{x},0)=B_{0}[3sz\mathbf{1}_{s}+(5-6s^{2}-3z^{2})\mathbf{1}_{z}],\quad U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,0)=0, & & \displaystyle\end{eqnarray}$$

for which $K\;(=U_{\unicode[STIX]{x1D719}}^{\text{G}}(1,0))$ is zero; see (4.3c ). Other $K$ are considered in § 5.3.

As in paper 2, we made Braginsky’s (Braginsky (Reference Braginsky1964)) choice of

(5.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}(s,z)=\frac{729}{16}\unicode[STIX]{x1D6FC}_{0}r^{7}z(1-r^{2})^{2}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}_{0}$ is the maximum of $\unicode[STIX]{x1D6FC}$ . For $\unicode[STIX]{x1D6FC}_{0}\geqslant \unicode[STIX]{x1D6FC}_{T}\approx 13.80$ , the solutions equilibrate to non-trivial steady states; for $\unicode[STIX]{x1D6FC}_{0}<\unicode[STIX]{x1D6FC}_{T}$ , $\boldsymbol{B}$ eventually decays to zero. For the $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D714}$ dynamos of paper 2, supercritical solutions equilibrate to steadily oscillating dynamos. In § 5.3, we shall compare MTT solutions with the OTT dynamos of paper 2 for $\unicode[STIX]{x1D6FC}^{2}$ mean-field models driven by (5.3), and shall also display results for MTT $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D714}$ dynamos.

5.2 Numerical methods

We used a semi-discrete method to solve the MTT model. We first discretized the equations in space and obtained a system of ordinary differential equations in time. These were solved using a Runge–Kutta scheme of Shu (Reference Shu1988).

The spatial discretization was carried out on three overlapping grids, which were described in detail in appendix B of Wu & Roberts (Reference Wu and Roberts2014). They involved

1. (i) a ‘square’ grid covering the region, $-\unicode[STIX]{x0394}s\leqslant s\leqslant 0.65$ , $-0.65\leqslant z\leqslant 0.65$ , where $\unicode[STIX]{x0394}s$ is the grid spacing in $s$ ;

2. (ii) a polar grid covering $0.5\leqslant r\leqslant 1$ and $-\unicode[STIX]{x0394}\unicode[STIX]{x1D703}\leqslant \unicode[STIX]{x1D703}\leqslant \unicode[STIX]{x03C0}+\unicode[STIX]{x0394}\unicode[STIX]{x1D703}$ , where $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}$ is the grid spacing in $\unicode[STIX]{x1D703}$ ; and

3. (iii) a cylindrical grid, as specified by Nakajima & Roberts (Reference Nakajima and Roberts1995a ,Reference Nakajima and Roberts b ), covering the entire computational domain.

There are 4 primary variables in the MTT calculations: $U_{\unicode[STIX]{x1D719}}^{\text{G}}$ and three components of $\boldsymbol{B}$ . Other quantities involving $\boldsymbol{U}$ are auxiliary variables that can be obtained from the 4 primary variables; the cylindrical grid was used mainly for that task. The polar grid and the square grid were used to solve the induction equation while enforcing the dynamo conditions (3.1h,i ). (If the polar grid alone had been used, a grid singularity would have occurred at $r=0$ . By supplementing the polar grid by the square grid, this difficulty is avoided.)

To preserve the integrity (2.2) of our reference frame, we used the discrete expression of (4.3a ) which, for our cylindrical grid, is

(5.4a ) $$\begin{eqnarray}\displaystyle E_{\unicode[STIX]{x1D702}}s_{i+1/2}{\mathcal{A}}(s_{i+1/2})\frac{\unicode[STIX]{x2202}U_{\unicode[STIX]{x1D719}}^{\text{G}}(s_{i+1/2},t)}{\unicode[STIX]{x2202}t}=\frac{\unicode[STIX]{x1D6E4}_{z}(s_{i+1},t)-\unicode[STIX]{x1D6E4}_{z}(s_{i},t)}{\unicode[STIX]{x0394}s}, & & \displaystyle\end{eqnarray}$$

where $s_{i}$ and $s_{i+1/2}$ refer to the integer and half-integer grid values, respectively. Thus the change in the total angular momentum in time becomes

(5.4b ) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}\left(\mathop{\sum }_{i}s_{i+1/2}{\mathcal{A}}(s_{i+1/2})U_{\unicode[STIX]{x1D719}}^{\text{G}}(s_{i+1/2},t)\unicode[STIX]{x0394}s\right)=\frac{\unicode[STIX]{x1D6E4}_{z}(1,t)-\unicode[STIX]{x1D6E4}_{z}(0,t)}{E_{\unicode[STIX]{x1D702}}}=0. & & \displaystyle\end{eqnarray}$$

The results presented below used the same grid system and the same number of grid points as Wu & Roberts (Reference Wu and Roberts2014). The calculations are second-order accurate in both space and time. Numerical stability limits the magnitude of the attainable $E_{\unicode[STIX]{x1D702}}$ which, for the grid just specified, is $E_{\unicode[STIX]{x1D702}}\geqslant 10^{-4}$ for the $\unicode[STIX]{x1D6FC}^{2}$ model and $E_{\unicode[STIX]{x1D702}}\geqslant 10^{-5}$ for the $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D714}$ model. To reach smaller values of $E_{\unicode[STIX]{x1D702}}$ , the number of grid points would have to be increased. We did not analyse the system to determine how the time step is constrained. The time step was found by trial and error guided by the Courant condition and by the constraint imposed by magnetic field diffusion.

5.3 MTT results for dipolar models

Most of the figures shown below are for $t$ large enough for the solutions to be close to their final steady or statistically steady states; we term these ‘saturated states’. We consider first dipolar $\unicode[STIX]{x1D6FC}^{2}$ models and compare them with the OTT dynamos of (2: § 5).

Because this paper uses MTT instead of the OTT of paper 2, the evolution of $\boldsymbol{B}$ from (2.1b,c ) differs from that of paper 2: evolution in MTT depends on $E_{\unicode[STIX]{x1D702}}$ but that in OTT does not. This is apparent in figure 2 which shows the internal magnetic energy, $M(t)$ , for $\unicode[STIX]{x1D6FC}_{0}=14.5$ for 2 values of $E_{\unicode[STIX]{x1D702}}$ . Also shown by the dotted curve is the corresponding $M(t)$ from paper 2 for OTT. The 3 curves approach one another as $t$ increases, and each approaches the same steady-state value $M(\infty )\approx 0.128$ . The existence of the torsional wave in MTT is evident; $M(t)$ for both MTT solutions overshoots $M(\infty )$ before recovering. In comparison, the $M(t)$ of the OTT solution increases monotonically to $M(\infty )$ (apart from an initial decrease whose origin was explained in paper 2). The OTT solution shown was extracted from paper 2.

Paper 1 pinpointed the way ohmic resistance damps out torsional waves; if there is no source to maintain a torsional wave, resistivity will remove it on the $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ time scale ( $r_{o}^{2}/\unicode[STIX]{x1D702}$ ). Our initial non-zero $U_{\unicode[STIX]{x1D719}}^{\text{G}}$ launches a torsional wave but does not maintain it, so its ‘half-life’ is only approximately $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ . Over this period, there will be approximately $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}/\unicode[STIX]{x1D70F}_{\text{tort}}$ maxima and minima of $M(t)$ , where $\unicode[STIX]{x1D70F}_{\text{tort}}\propto r_{o}/V_{\text{tort}}$ , the period of the torsional wave, is proportional to $E_{\unicode[STIX]{x1D702}}^{1/2}$ (see § 4). In the geophysical case, $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}/\unicode[STIX]{x1D70F}_{\text{tort}}\sim 10^{5}$ according to § 4.1, and the waves traverse the core so many times before being attenuated that the mechanism of their excitation becomes the primary target of theory. Such small $E_{\unicode[STIX]{x1D702}}$ and large $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}/\unicode[STIX]{x1D70F}_{\text{tort}}$ are not practical for dynamo simulations, and therefore far fewer extrema of $M$ can be anticipated. Numerical difficulties increase as $E_{\unicode[STIX]{x1D702}}$ diminishes, and our limited computer resources made it impossible for us to obtain reliable results for $E_{\unicode[STIX]{x1D702}}$ smaller than $10^{-4}$ . Figure 2 shows only one maximum of $M$ convincingly. Future MTT models at smaller $E_{\unicode[STIX]{x1D702}}$ should see more.

While $\boldsymbol{B}$ is small, it grows exponentially: $\boldsymbol{B}\propto \text{e}^{\unicode[STIX]{x1D706}t}$ and the $M(t)\propto \text{e}^{2\unicode[STIX]{x1D706}t}$ as shown in figure 2. This may be seen in the ‘straight line’ segments of the three $M(t)$ of figure 2. During this ‘linear phase’, $\unicode[STIX]{x1D706}=2.55$ for the OTT calculation and $\unicode[STIX]{x1D706}=5.50$ for both MTT cases. The different slopes of the straight line segments of MTT and OTT are apparent in figure 2. They arises because OTT solutions must satisfy Taylor’s constraint as $\boldsymbol{B}$ and $M(t)$ grow, but MTT solutions do not. Therefore the $\unicode[STIX]{x1D706}$ of the MTT solutions coincides with the $\unicode[STIX]{x1D706}$ of the kinematic dynamo for the assigned $\unicode[STIX]{x1D6FC}_{0}=14.5$ , but the $\unicode[STIX]{x1D706}$ of the OTT solution is smaller and vanishes for the marginal $\unicode[STIX]{x1D6FC}_{T}\approx 13.8$ .

The remaining three figures for the $\unicode[STIX]{x1D6FC}^{2}$ dynamos show only the geostrophic part, $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)$ , of the MTT flows.

Figure 3 is for $\unicode[STIX]{x1D6FC}_{0}=14.5$ , $E_{\unicode[STIX]{x1D702}}=0.01$ and $K\;(=U_{\unicode[STIX]{x1D719}}^{\text{G}}(1,0))=0$ ; see (4.3c ). It displays $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)$ for $t=34.2$ , 41.1, 47.9 and 54.8. There are 2 minima of $U_{\unicode[STIX]{x1D719}}^{\text{G}}$ . The one nearer to $s=1$ is denoted by $(U_{\unicode[STIX]{x1D719}}^{\text{G}})_{\text{min}}$ and its location by $s_{\text{min}}$ . For $t=34.2$ , 41.1, 47.9 and 54.8, $(s_{\text{min}},(U_{\unicode[STIX]{x1D719}}^{\text{G}})_{\text{min}})=(0.830,-1.635)$ , $(0.835,-1.777)$ , $(0.842,-1.898)$ and $(0.845,-2.000)$ , respectively. In the interval $s_{\text{min}}\leqslant s\leqslant 1$ , $U_{\unicode[STIX]{x1D719}}^{\text{G}}$ increases dramatically to 0. As $t$ increases, so does the depth $-(U_{\unicode[STIX]{x1D719}}^{\text{G}})_{\text{min}}$ of the minimum and its location $s_{\text{min}}$ . This trend suggests that, as $t\rightarrow \infty$ , $s_{\text{min}}$ will tend to 1 and $-(U_{\unicode[STIX]{x1D719}}^{\text{G}})_{\text{min}}$ will tend to $-(U_{\unicode[STIX]{x1D719}}^{\text{G}})_{\text{O}\text{T}\text{T}}(1,\infty )\approx 4.8$ .

This figure and the next indicate how the paradox of § 4 may be resolved: Because the rate at which $U_{\unicode[STIX]{x1D719}}^{\text{G}}(1,\infty )$ approaches $K$ is proportional to ${\mathcal{V}}_{\text{tort}}(s,t)$ and vanishes with it as $s\rightarrow 1$ , the differences $\boldsymbol{U}(\boldsymbol{x},t)-[\boldsymbol{U}(\boldsymbol{x},t)]_{\text{OTT}}$ and $\boldsymbol{B}(\boldsymbol{x},t)-[\boldsymbol{B}(\boldsymbol{x},t)]_{\text{OTT}}$ remain non-zero, becoming simple discontinuities on the equator, $(s,z)=(1,0)$ , at $t=\infty$ . There they have no dynamic or inductive effect, and do not affect $m_{z}^{\text{total}}$ .

Figure 3. The geostrophic flow, $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)$ , as a function of $s$ for the MTT $\unicode[STIX]{x1D6FC}^{2}$ model at different times; $\unicode[STIX]{x1D6FC}_{0}=14.5$ , $E_{\unicode[STIX]{x1D702}}=0.01$ and $U_{\unicode[STIX]{x1D719}}^{\text{G}}(1,t)=0$ . The red dashed-dotted curve is for $t=34.2$ , the magenta dotted curve is for $t=41.1$ , the green dashed curve is for $t=47.9$ , and the blue solid curve is for $t=54.8$ .

Figure 4. The geostrophic flow, $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,41.1)$ , as a function of $s$ at different $E_{\unicode[STIX]{x1D702}}$ for the MTT $\unicode[STIX]{x1D6FC}^{2}$ model; $\unicode[STIX]{x1D6FC}_{0}=14.5$ and $U_{\unicode[STIX]{x1D719}}^{\text{G}}(1,t)=0$ . The blue dotted curve is for $E_{\unicode[STIX]{x1D702}}=0.01$ , the green dashed curve is for $E_{\unicode[STIX]{x1D702}}=0.001$ and the red solid curve is for $E_{\unicode[STIX]{x1D702}}=0.0001$ . In addition, the black solid curve is for the OTT model for the same $\unicode[STIX]{x1D6FC}_{0}$ .

Figure 5. The geostrophic flow, $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,41.1)$ , as a function of $s$ at different $E_{\unicode[STIX]{x1D702}}$ for the MTT $\unicode[STIX]{x1D6FC}^{2}$ model; $\unicode[STIX]{x1D6FC}_{0}=14.5$ and $E_{\unicode[STIX]{x1D702}}=0.001$ . The magenta solid curve is for $U_{\unicode[STIX]{x1D719}}^{\text{G}}(1,0)=0$ , the curve of solid green squares is for $U_{\unicode[STIX]{x1D719}}^{\text{G}}(1,0)=-2.0$ , the red dashed-dotted curve is for $U_{\unicode[STIX]{x1D719}}^{\text{G}}(1,0)=-4.5$ and the black dotted curve is for $U_{\unicode[STIX]{x1D719}}^{\text{G}}(1,0)=-6.0$ . In addition, the blue dashed curve is for the OTT model for the same $\unicode[STIX]{x1D6FC}_{0}$ ; for this OTT model, $U_{\unicode[STIX]{x1D719}}^{\text{G}}\approx -4.8$ at $s=1$ .

Figure 6. The averaged internal magnetic energy, $\langle M\rangle$ , plotted in circles as a function of $E_{\unicode[STIX]{x1D702}}$ , for saturated MTT $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D714}$ states; $\unicode[STIX]{x1D6FC}_{0}(=-\unicode[STIX]{x1D714}_{0}^{\prime })=425$ . The $\langle M\rangle$ value based on the OTT calculation with the same $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D714}^{\prime }$ is indicated by the asterisk.

Figure 7. Time dependence of the internal magnetic energy, $M(t)$ , for the MTT $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D714}$ model with $\unicode[STIX]{x1D6FC}_{0}(=-\unicode[STIX]{x1D714}^{\prime })=425$ : (a) at an early phase, and (b) in the final saturated state. The blue curves are for $E_{\unicode[STIX]{x1D702}}=10^{-4}$ and the red curves for $E_{\unicode[STIX]{x1D702}}=10^{-5}$ .

Figure 8. The power spectral density (PSD) of the internal magnetic energy, $M(t)$ , as a function of frequency for saturated states of the $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D714}$ model with $\unicode[STIX]{x1D6FC}_{0}(=-\unicode[STIX]{x1D714}^{\prime })=425$ : (a) in red, for MTT and $E_{\unicode[STIX]{x1D702}}=10^{-4}$ ; (b) in green, for MTT and $E_{\unicode[STIX]{x1D702}}=10^{5}$ ; (c) in blue for OTT. In case (a) the PSD has been divided by 40.

Figure 4 re-enforces this interpretation. This figure displays $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)$ as $E_{\unicode[STIX]{x1D702}}$ is reduced for fixed $t=41.1$ ; here again $\unicode[STIX]{x1D6FC}_{0}=14.5$ and $K=0$ . Displayed are plots for $E_{\unicode[STIX]{x1D702}}=10^{-2},10^{-3}$ and $10^{-4}$ . Also shown is the OTT solution for the same $\unicode[STIX]{x1D6FC}_{0}$ . Employing the same terminology used when describing figure 3, $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)$ attains $(U_{\unicode[STIX]{x1D719}}^{\text{G}})_{\text{min}}$ , which becomes deeper and moves to larger $s$ as $E_{\unicode[STIX]{x1D702}}$ is reduced, while hugging the OTT solution until it leaves it near $s=s_{\text{m}\text{i}\text{n}}$ to climb rapidly to the assigned $U_{\unicode[STIX]{x1D719}}^{\text{G}}(1,t)=K=0$ . The similarity of the progression of solutions for fixed $t$ and decreasing $E$ is related to the progression of solutions for fixed $E$ and increasing $t$ (figure 3) through the $E^{1/2}$ time scale of torsional waves identified in § 4. (The slight irregularity of $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,t)$ for $E_{\unicode[STIX]{x1D702}}=10^{-4}$ is due to the same numerical difficulty as mentioned above in connection with figure 2.)

Figure 5 shows how $U_{\unicode[STIX]{x1D719}}^{\text{G}}(s,41.1)$ for MTT solutions changes when $t$ and $E_{\unicode[STIX]{x1D702}}$ are fixed and only $K$ defined in (4.3c ) is varied; as before $\unicode[STIX]{x1D6FC}_{0}=14.5$ and $E_{\unicode[STIX]{x1D702}}=0.001$ ). Displayed are plots for $K=0$ , $-2.0$ , $-4.5$ and $-6.0$ . That there are substantial differences near $s=1$ where the boundary condition (4.3c ) is enforced is no surprise. The OTT solution, for which $[U_{\unicode[STIX]{x1D719}}^{\text{G}}(1)]_{\text{O}\text{T}\text{T}}\approx -4.8$ is also shown. The MTT plots for $K=-4.5$ and $-6.0$ straddle the OTT solution. This demonstrates consistency between the OTT and MTT solutions.

Our final 3 figures are for MTT $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D714}$ models but we also compare them with OTT $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D714}$ dynamos of (2: § 6). In all 3 cases, $\unicode[STIX]{x1D6FC}_{0}\;(=-\unicode[STIX]{x1D714}_{0}^{\prime })=425$ .

Figure 6 shows the mean $\langle M(t)\rangle$ of the evolved internal magnetic energy, $M(t)$ , as a function of $E_{\unicode[STIX]{x1D702}}$ . The inferred value of $\langle M(t)\rangle \approx 0.004$ for $E_{\unicode[STIX]{x1D702}}=0$ may be compared with the OTT value of $\langle M(t)\rangle \approx 0.0043$ : see (2: § 6). Interestingly, the amplitude of $\langle M(t)\rangle$ increases with $E_{\unicode[STIX]{x1D702}}$ , as the torsional wave speed (4.4d ) deceases.

Figure 7 shows $M(t)$ for $E_{\unicode[STIX]{x1D702}}=10^{-4}$ and $E_{\unicode[STIX]{x1D702}}=10^{-5}$ at an early stage in the evolution of the solution from its starting point (5.2a,b ) and in its final evolved state. The increase in the amplitude of $M(t)$ with time is evident. As it is impossible to extract quantitative information from this figure, or from (2: figure 14a) for OTT, we performed some spectral analyses:

Figure 8 shows the power spectral density (PSD) of the three $M(t)$ of figure 7 in their saturated states. The spectral peaks for $E_{\unicode[STIX]{x1D702}}=10^{-5}$ are sharper than for OTT and are significantly displaced from them.