2.1 The governing equations

The geometry of our model is a whole sphere of radius $R$ containing fluid of electrical conductivity $\unicode[STIX]{x1D70E}$ , density $\unicode[STIX]{x1D70C}$ , kinematic viscosity $\unicode[STIX]{x1D708}$ , thermal diffusivity $\unicode[STIX]{x1D705}$ , rotating about the $z$ -axis with angular velocity $\unicode[STIX]{x1D6FA}\,\hat{\boldsymbol{z}}$ , surrounded by an impenetrable electrical insulator. The body has a gravity profile $\boldsymbol{g}=g_{0}\boldsymbol{r}$ and thermal expansivity $\unicode[STIX]{x1D6FC}$ , which together supply thermal buoyancy. We present the equations pertinent to a buoyancy-driven flow under the Boussinesq approximation with future applications in mind, although in the calculations that follow we neglect the temperature equation and buoyancy effects. Since our intention in this work is to model the dynamics of Earth’s core over long time scales, it is essential that we non-dimensionalize appropriately. In particular, we scale length by the spherical radius $R$ , time by the Ohmic dissipation time scale $R^{2}/\unicode[STIX]{x1D702}$ , where $\unicode[STIX]{x1D702}=(\unicode[STIX]{x1D70E}\unicode[STIX]{x1D707}_{0})^{-1}$ and $\unicode[STIX]{x1D707}_{0}$ is the permeability of free space, magnetic field strength by $(2\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D707}_{0}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D702})^{1/2}$ and temperature by $\unicode[STIX]{x1D6FD}R$ , where $\unicode[STIX]{x1D6FD}$ is a typical temperature gradient. Then the non-dimensional equations governing the dynamics of the Earth’s core are (Fearn Reference Fearn1998)

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle R_{o}\left(\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}\right)+\hat{\boldsymbol{z}}\times \boldsymbol{u}=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D70B}+(\unicode[STIX]{x1D735}\times \boldsymbol{B})\times \boldsymbol{B}+q\,R_{a}\unicode[STIX]{x1D717}\boldsymbol{r}+E\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{B}}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D735}\times (\boldsymbol{u}\times \boldsymbol{B})+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{B}, & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D717}}{\unicode[STIX]{x2202}t}=-\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}+q\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D717}+h, & \displaystyle\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}=0, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}$ is the flow, $\boldsymbol{B}$ the magnetic field, $\unicode[STIX]{x1D717}$ the temperature, $\unicode[STIX]{x1D70B}$ the pressure and $h$ is an internal heat source. The non-dimensional parameters appearing are the modified Rayleigh number $R_{a}$ , the magnetic Rossby number $R_{o}$ , the Ekman number $E$ and Roberts number $q$ defined as

(2.5a-d ) $$\begin{eqnarray}\displaystyle R_{a}=\frac{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}g_{0}R}{2\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D705}},\quad R_{o}=\frac{\unicode[STIX]{x1D702}}{2\unicode[STIX]{x1D6FA}R^{2}},\quad E=\frac{\unicode[STIX]{x1D708}}{2\unicode[STIX]{x1D6FA}R^{2}},\quad q=\frac{\unicode[STIX]{x1D705}}{\unicode[STIX]{x1D702}}, & & \displaystyle\end{eqnarray}$$

the latter three of which take typical values $10^{-9}$ , $10^{-15}$ and $10^{-5}$ , respectively. Taylor (Reference Taylor1963) recognized the smallness of both $R_{o}$ and $E$ , and a dominant balance between Coriolis, pressure, Lorentz and buoyancy leads to the magnetostrophic approximation, here specialized to the case $R_{a}=0$ and described by

(2.6) $$\begin{eqnarray}\displaystyle \hat{\boldsymbol{z}}\times \boldsymbol{u}=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D70B}+(\unicode[STIX]{x1D735}\times \boldsymbol{B})\times \boldsymbol{B}, & & \displaystyle\end{eqnarray}$$

which replaces (2.1). We note parenthetically that the equations given by the $R_{o}=0$ limit of (2.1) were those solved in the seminal work of Glatzmaier & Roberts (Reference Glatzmaier and Roberts1995a ,Reference Glatzmaier and Roberts b ). In the case of axisymmetry, equation (2.2) has no self-sustained solutions. However, in this case, mean-field effects are often invoked and the induction equation takes the more general form

(2.7) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{B}}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D735}\times (\boldsymbol{u}\times \boldsymbol{B}+\boldsymbol{f}(\boldsymbol{B},\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D714}))+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{B}, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{f}$ encapsulates the effects of helical twisting (the $\unicode[STIX]{x1D6FC}$ -effect) and large-scale differential rotation (the $\unicode[STIX]{x1D714}$ -effect). One important assignment for $\boldsymbol{f}$ is as the classic isotropic $\unicode[STIX]{x1D6FC}$ -effect, which takes the form of an electromotive force (EMF) ${\mathcal{E}}$ , of the form

(2.8) $$\begin{eqnarray}\displaystyle {\mathcal{E}}=\unicode[STIX]{x1D6FC}\boldsymbol{B}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}$ is a spatially dependent prescribed scalar.

Equations (2.6) and (2.7) are solved within the unit sphere, $V$ , augmented by the boundary conditions $\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{r}=0$ ,  $[\boldsymbol{B}]=\mathbf{0}$  and $[\boldsymbol{r}\times \boldsymbol{E}]=\mathbf{0}$ on $r=1$ where $[~]$ indicates the jump across the boundary, $\boldsymbol{E}$ is the electric field and the magnetic field is potential in the exterior region $\hat{V}$ described by $r>1$ .

Equation (2.6) only has solutions when (1.1) is satisfied (Taylor Reference Taylor1963), and only gives the ageostrophic component of flow, $\boldsymbol{u}_{m}$ , also termed the magnetic wind (as signified by the subscript $m$ ). To this flow can be added an arbitrary geostrophic flow $\boldsymbol{u}_{g}=u_{g}(s)\hat{\unicode[STIX]{x1D753}}$ that is made specific by the following requirements: (i)  ${\mathcal{T}}$ remains zero at all times; (ii)  $u_{g}(s)\hat{\unicode[STIX]{x1D753}}$ is finite on the axis; and (iii) angular momentum is conserved. The problem as envisaged by Taylor now has a complete and unique solution.

It is noteworthy that, within a Taylor state, the geostrophic component directly contributes zero to the change in magnetic energy. This can be seen by taking the dot product of the induction equation (2.7) with $\boldsymbol{B}$ and integrating over all space:

(2.9) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}\int _{V+\hat{V}}\frac{1}{2}\boldsymbol{B}^{2}\,\text{d}V=\int _{V+\hat{V}}\boldsymbol{B}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\boldsymbol{B}}{\unicode[STIX]{x2202}t}\,\text{d}V=\int _{V+\hat{V}}\boldsymbol{B}\boldsymbol{\cdot }[\unicode[STIX]{x1D735}\times (\boldsymbol{u}\times \boldsymbol{B}+{\mathcal{E}})+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{B}]\,\text{d}V. & & \displaystyle\end{eqnarray}$$

The contribution from the geostrophic flow is

(2.10) $$\begin{eqnarray}\displaystyle \int _{V}\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times [\boldsymbol{u}_{g}\times \boldsymbol{B}]\,\text{d}V=\int _{V}\unicode[STIX]{x1D735}\boldsymbol{\cdot }[(\boldsymbol{u}_{g}\times \boldsymbol{B})\times \boldsymbol{B}]\,\text{d}V-\int _{V}\boldsymbol{u}_{g}\boldsymbol{\cdot }[(\unicode[STIX]{x1D735}\times \boldsymbol{B})\times \boldsymbol{B}]\,\text{d}V. & & \displaystyle\end{eqnarray}$$

The first term on the right-hand side can be written as a surface integral:

(2.11) $$\begin{eqnarray}\displaystyle \int _{r=1}[(\boldsymbol{u}_{g}\times \boldsymbol{B})\times \boldsymbol{B}]\boldsymbol{\cdot }\text{d}\boldsymbol{S}=\int _{r=1}[(\boldsymbol{u}_{g}\boldsymbol{\cdot }\boldsymbol{B})\,\boldsymbol{B}-|\boldsymbol{B}|^{2}\,\boldsymbol{u}_{g}]\boldsymbol{\cdot }\text{d}\boldsymbol{S}=\int _{r=1}[\boldsymbol{u}_{g}\boldsymbol{\cdot }\boldsymbol{B}]B_{r}\,\text{d}S, & & \displaystyle\end{eqnarray}$$

since $u_{r}=0$ on $r=1$ . Owing to the form of the magnetic boundary conditions and the assumed axisymmetry, $\boldsymbol{B}$ has no azimuthal component at $r=1$ , which renders the whole term zero. Thus the contribution of $u_{g}$ to the evolution of magnetic energy is

(2.12) $$\begin{eqnarray}\displaystyle -\int _{V}\boldsymbol{u}_{g}\boldsymbol{\cdot }[(\unicode[STIX]{x1D735}\times \boldsymbol{B})\times \boldsymbol{B}]\,\text{d}V=-\int _{0}^{1}A(s){\mathcal{T}}\,\boldsymbol{u}_{g}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D753}}\,\text{d}s=0, & & \displaystyle\end{eqnarray}$$

since the Taylor integral ${\mathcal{T}}$ is zero, where $A(s)=4\unicode[STIX]{x03C0}s\sqrt{1-s^{2}}$ is the area of a cylinder. The special form of the magnetic field (a Taylor state) together with the form of the geostrophic flow serves to annul the usual $\unicode[STIX]{x1D714}$ -effect. This result was shown by Fearn, Roberts & Soward (Reference Fearn, Roberts, Soward, Galdi and Straughan1988, pp. 185–186) and contrasts with the observation that for weakly viscous models the effect of the geostrophic flow on the magnetic energy is negative-semidefinite. The corresponding analysis with $u_{m}$ is not possible, as this flow is not purely azimuthal, and the boundary term remains.

2.2 Maintaining a Taylor state using optimal control

Key to our implementation is the ability to maintain ${\mathcal{T}}=0$ at all times, which we do by drawing on the theory of optimal control. Beginning with the foundations laid down by Athans & Falb (Reference Athans and Falb1966) and Lions (Reference Lions1971), developments of this technique are now widely used in numerous fields (e.g. Tarantola Reference Tarantola1984; Talagrand & Courtier Reference Talagrand and Courtier1987; Tromp, Tape & Lui Reference Tromp, Tape and Lui2005; Pringle & Kerswell Reference Pringle and Kerswell2010); a recent review is given by Kerswell, Pringle & Willis (Reference Kerswell, Pringle and Willis2014). In this section we introduce the theory in a continuous setting in order to illustrate the basic concepts. Actual implementation details, which depend on the choice of a numerical time-stepping algorithm, differ slightly and are described in appendix E. Our intention here is to illustrate the general concepts in an accessible manner.

We consider the problem of taking a single time step, from $t=0$ , at which ${\mathcal{T}}=0$ , to $t=\unicode[STIX]{x0394}t$ . At time $t=0$ the ageostrophic component $\boldsymbol{u}_{m}$ is known from (2.6), which uses the structure of the magnetic field at time $t=0$ and is assumed constant over the time step; the geostrophic component $u_{g}$ , required to remain on the Taylor manifold, is similarly assumed constant and is to be determined.

The variational calculus we employ (Roberts Reference Roberts1960) will be based on integrals taken over all space, $V+\hat{V}$ , denoted by $\langle \;\rangle$ , anticipating the simple form of the adjoint equations that arise in this case (Chen et al. Reference Chen, Herreman, Li, Livermore, Luo and Jackson2018, and appendix C):

(2.13) $$\begin{eqnarray}\displaystyle \langle \bullet \rangle =\int _{V\cup \hat{V}}\bullet \,\text{d}V. & & \displaystyle\end{eqnarray}$$

Since the Lorentz force $\boldsymbol{J}\times \boldsymbol{B}$ of (1.1) vanishes in $\hat{V}$ , we note that we can extend the definition of Taylor’s constraint ${\mathcal{T}}$ both in $z$ and in $s$ to all space:

(2.14) $$\begin{eqnarray}\displaystyle {\mathcal{T}}(s,t)=\int _{0}^{2\unicode[STIX]{x03C0}}\!\int _{-\infty }^{\infty }[\boldsymbol{J}\times \boldsymbol{B}]_{\unicode[STIX]{x1D719}}\,\text{d}z\,s\,\text{d}\unicode[STIX]{x1D719}=\{\boldsymbol{J}\times \boldsymbol{B}\},\quad 0\leqslant s<\infty , & & \displaystyle\end{eqnarray}$$

where ${\mathcal{T}}$ vanishes at $s=1$ and joins continuously to its value of zero in $s>1$ . In the above, we have introduced a notation $\{~\}$ to signify integration of the $\unicode[STIX]{x1D719}$ component over an infinite cylinder, $C_{\infty }(s)$ , the extension of $C(s)$ at radius $s$ :

(2.15) $$\begin{eqnarray}\displaystyle \{\boldsymbol{q}\}=\int _{C_{\infty }(s)}\hat{\unicode[STIX]{x1D753}}\boldsymbol{\cdot }\boldsymbol{q}\,\text{d}z\,s\,\text{d}\unicode[STIX]{x1D719}. & & \displaystyle\end{eqnarray}$$

We are now in a position to introduce the core of our methodology, involving the augmented objective function or Lagrangian, $\unicode[STIX]{x1D712}^{2}$ (redefined from (1.2)), which includes the constraints that must be obeyed by the magnetic field: those of the solenoidal condition and the equation of magnetic induction. For the sake of clarity of exposition, we specialize to the case $\boldsymbol{f}=\unicode[STIX]{x1D6FC}\boldsymbol{B}$ , the more general case proceeding along similar lines, and revert to dimensional equations. Our augmented Lagrangian is

(2.16) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D712}^{2} & = & \displaystyle \frac{1}{2}\int _{0}^{\infty }w(s){\mathcal{T}}(s,\unicode[STIX]{x0394}t)^{2}\,\text{d}s+\int _{0}^{\unicode[STIX]{x0394}t}\langle p^{\dagger }\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}\rangle \,\text{d}t\nonumber\\ \displaystyle & & \displaystyle +\,\int _{0}^{\unicode[STIX]{x0394}t}\left\langle \boldsymbol{B}^{\dagger }\boldsymbol{\cdot }\left[\frac{\unicode[STIX]{x2202}\boldsymbol{B}}{\unicode[STIX]{x2202}t}-\unicode[STIX]{x1D735}\times [(\boldsymbol{u}_{m}+\boldsymbol{u}_{g})\times \boldsymbol{B}+\unicode[STIX]{x1D6FC}\boldsymbol{B}]+\unicode[STIX]{x1D702}^{\prime }\unicode[STIX]{x1D735}\times \unicode[STIX]{x1D735}\times \boldsymbol{B}\right]\right\rangle \,\text{d}t,\end{eqnarray}$$

where ${\mathcal{T}}(s,\unicode[STIX]{x0394}t)$ is Taylor’s integral at $t=\unicode[STIX]{x0394}t$ , $\boldsymbol{B}^{\dagger }$ is a time-dependent Lagrange multiplier (also known as the adjoint magnetic field) and $p^{\dagger }$ is an adjoint pressure term that constrains $\boldsymbol{B}$ to be divergence-free. In the above, $\unicode[STIX]{x1D702}^{\prime }=\unicode[STIX]{x1D702}$ (a constant) in $V$ , and $\hat{\unicode[STIX]{x1D702}}$ in $\hat{V}$ , as defined in appendix C, and serves to allow the treatment of the magnetic field in the mantle in the limit of vanishingly small electrical conductivity (Namikawa & Matsushita Reference Namikawa and Matsushita1970).

Using the analysis of appendix C to rewrite (2.16) in modified form, we then take the variation, using the designation $D$ to signify the Gateaux derivative (see appendix B) and noting $D\boldsymbol{B}(0)=0$ , to obtain

(2.17) $$\begin{eqnarray}\displaystyle D\unicode[STIX]{x1D712}^{2} & = & \displaystyle \underbrace{D\int _{0}^{\infty }\frac{1}{2}w(s){\mathcal{T}}(s,\unicode[STIX]{x0394}t)^{2}\,\text{d}s}_{I_{1}}-\int _{0}^{\unicode[STIX]{x0394}t}\underbrace{\langle D\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{A}\rangle }_{I_{2}}\,\text{d}t-\int _{0}^{\unicode[STIX]{x0394}t}\underbrace{\langle \boldsymbol{B}^{\dagger }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times (Du_{g}\hat{\unicode[STIX]{x1D753}}\times \boldsymbol{B})\rangle }_{I_{3}}\,\text{d}t\nonumber\\ \displaystyle & & \displaystyle +\,\langle D\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{B}^{\dagger }(\unicode[STIX]{x0394}t)\rangle ,\end{eqnarray}$$

where $\boldsymbol{A}$ defines the following adjoint induction equation derived in appendix C,

(2.18) $$\begin{eqnarray}\displaystyle \boldsymbol{A}\equiv \frac{\unicode[STIX]{x2202}\boldsymbol{B}^{\dagger }}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\times \boldsymbol{B}^{\dagger }\times \boldsymbol{u}+\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D735}\times \boldsymbol{B}^{\dagger }+\unicode[STIX]{x1D702}^{\prime }\unicode[STIX]{x1D6FB}^{2}\boldsymbol{B}^{\dagger }+\unicode[STIX]{x1D735}p^{\dagger }=0, & & \displaystyle\end{eqnarray}$$

similar to that used in an allied context by Chen et al. (Reference Chen, Herreman, Li, Livermore, Luo and Jackson2018). The adjoint field $\boldsymbol{B}^{\dagger }$ is divergenceless and satisfies the same boundary conditions as $\boldsymbol{B}$ ; it is governed by the adjoint induction equation that operates in reverse time and is driven by initial conditions given at $t=\unicode[STIX]{x0394}t$ , arising from the term $I_{1}$ in (2.17), which we will analyse now.

The variation $I_{1}$ becomes

(2.19) $$\begin{eqnarray}\displaystyle D\int _{0}^{\infty }\frac{1}{2}w(s){\mathcal{T}}(s,\unicode[STIX]{x0394}t)^{2}\,\text{d}s=\int _{0}^{\infty }w(s){\mathcal{T}}\,D{\mathcal{T}}\,\text{d}s. & & \displaystyle\end{eqnarray}$$

Dropping the explicit time $\unicode[STIX]{x0394}t$ for convenience, the integrand is

(2.20) $$\begin{eqnarray}\displaystyle w(s){\mathcal{T}}\,D{\mathcal{T}}=w(s){\mathcal{T}}(s)\{\unicode[STIX]{x1D735}\times D\boldsymbol{B}\times \boldsymbol{B}+\unicode[STIX]{x1D735}\times \boldsymbol{B}\times D\boldsymbol{B}\}. & & \displaystyle\end{eqnarray}$$

Including the integral over $s$ , considering the first term on the right-hand side, we have

(2.21) $$\begin{eqnarray}\displaystyle & & \displaystyle \langle w(s){\mathcal{T}}(s)\hat{\unicode[STIX]{x1D753}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times D\boldsymbol{B}\times \boldsymbol{B}\rangle \nonumber\\ \displaystyle & & \displaystyle \quad =\langle \unicode[STIX]{x1D735}\times D\boldsymbol{B}\boldsymbol{\cdot }[\boldsymbol{B}\times w(s){\mathcal{T}}(s)\hat{\unicode[STIX]{x1D753}}]\rangle \nonumber\\ \displaystyle & & \displaystyle \quad =\langle \unicode[STIX]{x1D735}\boldsymbol{\cdot }(D\boldsymbol{B}\times [\boldsymbol{B}\times w(s){\mathcal{T}}(s)\hat{\unicode[STIX]{x1D753}}])\rangle +\langle D\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times [\boldsymbol{B}\times w(s){\mathcal{T}}(s)\hat{\unicode[STIX]{x1D753}}]\rangle \nonumber\\ \displaystyle & & \displaystyle \quad =\oint \hat{\boldsymbol{n}}\boldsymbol{\cdot }D\boldsymbol{B}\times [\boldsymbol{B}\times w(s){\mathcal{T}}(s)\hat{\unicode[STIX]{x1D753}}]\,\text{d}S+\langle D\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times [\boldsymbol{B}\times w(s){\mathcal{T}}(s)\hat{\unicode[STIX]{x1D753}}]\rangle .\end{eqnarray}$$

Looking at the first of the terms above (the surface term) we find it is

(2.22) $$\begin{eqnarray}\displaystyle \oint \hat{\boldsymbol{n}}\boldsymbol{\cdot }D\boldsymbol{B}\times [\boldsymbol{B}\times w(s){\mathcal{T}}(s)\hat{\unicode[STIX]{x1D753}}]\,\text{d}S & = & \displaystyle \oint D\boldsymbol{B}\boldsymbol{\cdot }(\boldsymbol{B}\times w(s){\mathcal{T}}(s)\hat{\unicode[STIX]{x1D753}})\times \hat{\boldsymbol{n}}\,\text{d}S\nonumber\\ \displaystyle & = & \displaystyle \oint D\boldsymbol{B}\boldsymbol{\cdot }w(s){\mathcal{T}}(s)B_{r}\hat{\unicode[STIX]{x1D753}}\,\text{d}S=0,\end{eqnarray}$$

since we evaluate the surface integral at infinity where $B_{r}$ goes to zero at least as fast as $r^{-3}$ . We note that there is no contribution from the surface $r=1$ , since all quantities in the integrand (in particular, both $B_{r}$ and ${\mathcal{T}}(s)$ ) are continuous.

The second term arising in (2.20) is much simpler:

(2.23) $$\begin{eqnarray}\displaystyle \langle w(s){\mathcal{T}}(s)\hat{\unicode[STIX]{x1D753}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{B}\times D\boldsymbol{B}\rangle =\langle D\boldsymbol{B}\boldsymbol{\cdot }w(s){\mathcal{T}}(s)\hat{\unicode[STIX]{x1D753}}\times (\unicode[STIX]{x1D735}\times \boldsymbol{B})\rangle . & & \displaystyle\end{eqnarray}$$

Collecting together the term at $t=\unicode[STIX]{x0394}t$ from (2.17), the last term of (2.21) and (2.23), this now gives us the initial value for $\boldsymbol{B}^{\dagger }$ at time $\unicode[STIX]{x0394}t$ , from which we integrate backwards:

(2.24) $$\begin{eqnarray}\displaystyle -\boldsymbol{B}^{\dagger }(\unicode[STIX]{x0394}t)=\unicode[STIX]{x1D735}\times (\boldsymbol{B}\times w(s){\mathcal{T}}(s)\hat{\unicode[STIX]{x1D753}})+w(s){\mathcal{T}}(s)\hat{\unicode[STIX]{x1D753}}\times (\unicode[STIX]{x1D735}\times \boldsymbol{B}). & & \displaystyle\end{eqnarray}$$

Since $\boldsymbol{B}$ and $\boldsymbol{B}^{\dagger }$ are known over the interval $[0,\unicode[STIX]{x0394}t]$ , we are now in a position to use $I_{3}$ to determine the derivative of the objective function with respect to $u_{g}$ , our ultimate goal.

Manipulation of $I_{3}$ shows the following:

(2.25) $$\begin{eqnarray}\displaystyle & & \displaystyle D\int _{0}^{\unicode[STIX]{x0394}t}\!\langle \boldsymbol{B}^{\dagger }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times (u_{g}\hat{\unicode[STIX]{x1D753}}\times \boldsymbol{B})\rangle \,\text{d}t\nonumber\\ \displaystyle & & \displaystyle \quad =\int _{0}^{\unicode[STIX]{x0394}t}\langle \boldsymbol{B}^{\dagger }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times (Du_{g}\hat{\unicode[STIX]{x1D753}}\times \boldsymbol{B})\rangle \,\text{d}t\nonumber\\ \displaystyle & & \displaystyle \quad =\int _{0}^{\unicode[STIX]{x0394}t}\!\langle (Du_{g}\hat{\unicode[STIX]{x1D753}}\times \boldsymbol{B})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{B}^{\dagger }-\unicode[STIX]{x1D735}\boldsymbol{\cdot }[\boldsymbol{B}^{\dagger }\times (Du_{g}\hat{\unicode[STIX]{x1D753}}\times \boldsymbol{B})]\rangle \,\text{d}t\nonumber\\ \displaystyle & & \displaystyle \quad =\int _{0}^{\unicode[STIX]{x0394}t}\langle Du_{g}\hat{\unicode[STIX]{x1D753}}\boldsymbol{\cdot }\boldsymbol{B}\times \unicode[STIX]{x1D735}\times \boldsymbol{B}^{\dagger }\rangle \,\text{d}t-\int _{0}^{\unicode[STIX]{x0394}t}\oint \hat{\boldsymbol{r}}\boldsymbol{\cdot }[\boldsymbol{B}^{\dagger }\times (Du_{g}\hat{\unicode[STIX]{x1D753}}\times \boldsymbol{B})]\,\text{d}t\nonumber\\ \displaystyle & & \displaystyle \quad =\int _{0}^{\unicode[STIX]{x0394}t}\int _{0}^{1}\{\boldsymbol{B}\times \unicode[STIX]{x1D735}\times \boldsymbol{B}^{\dagger }\}Du_{g}\,\text{d}s\,\text{d}t,\end{eqnarray}$$

where, in the penultimate line the last term vanishes. The reason for this is that it is $D\boldsymbol{E}_{\Vert }\boldsymbol{\cdot }\hat{\boldsymbol{r}}\times \boldsymbol{B}^{\dagger }$ , where $D\boldsymbol{E}_{\Vert }$ represents the incremental horizontal electric field created by $Du_{g}$ , and the analysis of appendix C shows that this must match to the electric field outside the core.

We are now ready to extract the relevant derivative of $\unicode[STIX]{x1D712}^{2}$ from the volume integral by the following definition,

(2.26) $$\begin{eqnarray}\displaystyle D\unicode[STIX]{x1D712}^{2}=\langle \unicode[STIX]{x1D735}_{u_{g}}\unicode[STIX]{x1D712}^{2}Du_{g}(s)\rangle =\int _{0}^{1}A(s)\unicode[STIX]{x1D735}_{u_{g}}\unicode[STIX]{x1D712}^{2}Du_{g}(s)\,\text{d}s\,\text{d}t, & & \displaystyle\end{eqnarray}$$

where

(2.27) $$\begin{eqnarray}\displaystyle A(s)=4\unicode[STIX]{x03C0}s\sqrt{1-s^{2}} & & \displaystyle\end{eqnarray}$$

is the area of the cylinder $C(s)$ at radius $s$ . Equation (2.25) shows the relevant derivative to be

(2.28) $$\begin{eqnarray}\displaystyle A(s)\unicode[STIX]{x1D735}_{u_{g}}\unicode[STIX]{x1D712}^{2} & = & \displaystyle -\int _{0}^{\unicode[STIX]{x0394}t}\{\boldsymbol{B}\times \unicode[STIX]{x1D735}\times \boldsymbol{B}^{\dagger }\}\,\text{d}t,\end{eqnarray}$$

(2.29) $$\begin{eqnarray}\displaystyle & = & \displaystyle -\int _{0}^{\unicode[STIX]{x0394}t}\{\boldsymbol{B}\times \boldsymbol{J}^{\dagger }\}\,\text{d}t.\end{eqnarray}$$

The scheme that now ensues is classic and has found wide use in many areas of science; a recent review is given by Kerswell et al. (Reference Kerswell, Pringle and Willis2014). We take the opportunity to reiterate the relevant equations; these form the basis of an iterative scheme in which we gradually improve our estimate of $u_{g}$ whilst driving down the value of $\unicode[STIX]{x1D712}^{2}$ :

1. (i) The adjoint field $\boldsymbol{B}^{\dagger }$ is governed by the adjoint equation

   (2.30) $$\begin{eqnarray}\displaystyle -\frac{\unicode[STIX]{x2202}\boldsymbol{B}^{\dagger }}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D735}\times \boldsymbol{B}^{\dagger }\times \boldsymbol{u}+\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D735}\times \boldsymbol{B}^{\dagger }+\unicode[STIX]{x1D702}^{\prime }\unicode[STIX]{x1D6FB}^{2}\boldsymbol{B}^{\dagger }+\unicode[STIX]{x1D735}p^{\dagger }, & & \displaystyle\end{eqnarray}$$
   over the reverse time interval, $t\in [0,\unicode[STIX]{x0394}t]$ .

2. (ii) The terminal condition is

   (2.31) $$\begin{eqnarray}\displaystyle \boldsymbol{B}^{\dagger }(\unicode[STIX]{x0394}t)=-\unicode[STIX]{x1D735}\times (w\boldsymbol{B}\times {\mathcal{T}}\hat{\unicode[STIX]{x1D753}})-w{\mathcal{T}}\hat{\unicode[STIX]{x1D753}}\times (\unicode[STIX]{x1D735}\times \boldsymbol{B}). & & \displaystyle\end{eqnarray}$$

3. (iii) The downhill direction is

   (2.32) $$\begin{eqnarray}\displaystyle -\unicode[STIX]{x1D735}_{\boldsymbol{u}_{g}}\unicode[STIX]{x1D712}^{2}=\frac{1}{A(s)}\int _{0}^{\unicode[STIX]{x0394}t}\{\boldsymbol{B}\times (\unicode[STIX]{x1D735}\times \boldsymbol{B}^{\dagger })\}\,\text{d}t, & & \displaystyle\end{eqnarray}$$
   where the adjoint field $\boldsymbol{B}^{\dagger }$ is a divergence-free field and satisfies the same insulating boundary conditions as those of $\boldsymbol{B}$ .

4. (iv) Two additional conditions, namely a regularity condition, $\boldsymbol{u}_{g}=0$ at $s=0$ , and conservation of angular momentum

   (2.33) $$\begin{eqnarray}\displaystyle \int _{V}s\,\hat{\unicode[STIX]{x1D753}}\boldsymbol{\cdot }\boldsymbol{u}_{g}\,\text{d}V=0 & & \displaystyle\end{eqnarray}$$
   must be imposed and serve to uniquely determine $u_{g}$ .

The equations developed thus far capture the spirit of our method for the calculation of the sensitivity of $\unicode[STIX]{x1D712}^{2}$ with respect to $u_{g}$ , allowing the limit $\unicode[STIX]{x1D712}^{2}\rightarrow 0$ to be reached. Actual details of implementation differ somewhat, as they depend on the choice of a method for the numerical time stepping of the equations. We discuss these details in appendix E where we give an implementation specific to the Adams–Bashforth time-stepping scheme.

3.1 Spatial discretization

Our implementation relies on a fully spectral Galerkin method, which is expected to deliver rapid convergence for spatially smooth solutions. We build on the work of Livermore & Ierley (Reference Livermore and Ierley2009), Livermore (Reference Livermore2010) and Li, Jackson & Livermore (Reference Li, Jackson and Livermore2011). For the magnetic and velocity fields we use the classic Mie representation, which ensures that the solenoidality condition is satisfied. Using capital letters to represent quantities associated with the magnetic field and lower case for velocity field, we therefore expand the ageostrophic flow (the ‘magnetic wind’) and the magnetic field, $\boldsymbol{u}_{m}$ and $\boldsymbol{B}$ , using the poloidal–toroidal decomposition

(3.1) $$\begin{eqnarray}\displaystyle \boldsymbol{B}=\boldsymbol{S}+\boldsymbol{T}=\mathop{\sum }_{(n,l,m)}S_{(n,l,m)}\hat{\boldsymbol{S}}_{(n,l,m)}+T_{(n,l,m)}\hat{\boldsymbol{T}}_{(n,l,m)}, & & \displaystyle\end{eqnarray}$$

where

(3.2a,b ) $$\begin{eqnarray}\displaystyle \hat{\boldsymbol{S}}_{(n,l,m)}=\unicode[STIX]{x1D735}\times \unicode[STIX]{x1D735}\times [\unicode[STIX]{x1D6F7}_{n}^{l}(r)Y_{l}^{m}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})\hat{\boldsymbol{r}}],\quad \hat{\boldsymbol{T}}_{(n,l,m)}=\unicode[STIX]{x1D735}\times [\unicode[STIX]{x1D6F9}_{n}^{l}(r)Y_{l}^{m}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})\hat{\boldsymbol{r}}]. & & \displaystyle\end{eqnarray}$$

The basis functions have the property

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle \int _{V+\hat{V}}\hat{\boldsymbol{S}}_{(n_{1},l_{1},m_{1})}\boldsymbol{\cdot }\hat{\boldsymbol{S}}_{(n_{2},l_{2},m_{2})}\,\text{d}V=\int _{V}\hat{\boldsymbol{T}}_{(n_{1},l_{1},m_{1})}\boldsymbol{\cdot }\hat{\boldsymbol{T}}_{(n_{2},l_{2},m_{2})}\,\text{d}V=\unicode[STIX]{x1D6FF}_{n_{1},n_{2}}\unicode[STIX]{x1D6FF}_{l_{1},l_{2}}\unicode[STIX]{x1D6FF}_{m_{1},m_{2}}, & \displaystyle\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle \int _{V+\hat{V}}\hat{\boldsymbol{S}}_{(n_{1},l_{1},m_{1})}\boldsymbol{\cdot }\hat{\boldsymbol{T}}_{(n_{2},l_{2},m_{2})}\,\text{d}V=0. & \displaystyle\end{eqnarray}$$

In (3.2), $Y_{l}^{m}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ is a fully normalized spherical harmonic of degree $l$ and order $m$ whose squared integral over the sphere is unity, and $S_{(n,l,m)}$ and $T_{(n,l,m)}$ are the spectral coefficients for $\boldsymbol{B}$ where $n$ is a radial index. The Galerkin radial basis functions $\unicode[STIX]{x1D6F7}_{n}^{l}$ and $\unicode[STIX]{x1D6F9}_{n}^{l}$ are polynomials in which the required boundary conditions are encoded (see appendix D for details). Thus each $\hat{\boldsymbol{S}}$ and $\hat{\boldsymbol{T}}$ individually satisfies the boundary conditions. These basis functions have been constructed by using the fundamental Jones–Worland polynomials (Livermore, Jones & Worland Reference Livermore, Jones and Worland2007) as building blocks. Such functions have been analysed by Marti & Jackson (Reference Marti and Jackson2016), particularly with respect to their smooth differentiable nature at the origin of the coordinate system. A key feature of the basis when used in a time-stepping scheme is the fact that the Courant–Friedrichs–Lewy condition remains innocuous near the origin (Marti & Jackson Reference Marti and Jackson2016), something that is hard to achieve in other (e.g. finite difference) methods.

Similarly, the ageostrophic flow is represented as

(3.5) $$\begin{eqnarray}\displaystyle \boldsymbol{u}_{m}=\mathop{\sum }_{(n,l,m)}s_{(n,l,m)}\hat{\boldsymbol{s}}_{(n,l,m)}+t_{(n,l,m)}\hat{\boldsymbol{t}}_{(n,l,m)}, & & \displaystyle\end{eqnarray}$$

where

(3.6a,b ) $$\begin{eqnarray}\displaystyle \hat{\boldsymbol{s}}_{(n,l,m)}=\unicode[STIX]{x1D735}\times \unicode[STIX]{x1D735}\times [\unicode[STIX]{x1D6F1}_{n}^{l}(r)Y_{l}^{m}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})\hat{\boldsymbol{r}}],\quad \hat{\boldsymbol{t}}_{(n,l,m)}=\unicode[STIX]{x1D735}\times [\unicode[STIX]{x1D6EF}_{n}^{l}(r)Y_{l}^{m}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})\hat{\boldsymbol{r}}], & & \displaystyle\end{eqnarray}$$

and the radial Galerkin functions $\unicode[STIX]{x1D6F1}_{n}^{l}$ and $\unicode[STIX]{x1D6EF}_{n}^{l}$ are suitably defined in order that $\hat{\boldsymbol{t}}_{(n,l,m)}$ and $\hat{\boldsymbol{s}}_{(n,l,m)}$ satisfy both their boundary conditions and a similar orthonormality relation to (3.4) (see appendix D). It is worth noting that, of the four classes of basis function, only $\hat{\boldsymbol{S}}_{(n,l,m)}$ is non-zero in $\hat{V}$ .

The remaining geostrophic flow is represented as

(3.7) $$\begin{eqnarray}\displaystyle \boldsymbol{u}_{g}=\mathop{\sum }_{k}g_{k}\unicode[STIX]{x1D6EC}_{k}(s)\hat{\unicode[STIX]{x1D753}}, & & \displaystyle\end{eqnarray}$$

with coefficients $g_{k}$ , where the basis functions $\unicode[STIX]{x1D6EC}_{k}$ (see appendix D) are orthonormal when integrated over the volume of the unit sphere (recall $A(s)$ is defined in (2.27)),

(3.8) $$\begin{eqnarray}\displaystyle \int _{0}^{1}A(s)\unicode[STIX]{x1D6EC}_{j}(s)\unicode[STIX]{x1D6EC}_{k}(s)\,\text{d}s=\unicode[STIX]{x1D6FF}_{j,k}. & & \displaystyle\end{eqnarray}$$

As will be described in § 3.2, we note that the only component of flow that carries angular momentum is $\unicode[STIX]{x1D6EC}_{1}$ .

Because our basis functions $\hat{\boldsymbol{S}}$ , $\hat{\boldsymbol{T}}$ , $\hat{\boldsymbol{s}}$ and $\hat{\boldsymbol{t}}$ are defined in terms of spherical harmonics and regular polynomials in radius, they may be equivalently represented in Cartesian coordinates as polynomials. This not only allows us to calculate quantities such as (1.1) exactly, but also permits us to define a largest polynomial degree of the Cartesian components which defines a measure of spatial complexity. By using a triangular truncation for the radial index (Ivers, Jackson & Winch Reference Ivers, Jackson and Winch2015), we adopt a resolution uniform in each coordinate direction (see table 1). At each spherical harmonic degree $l$ , the maximum radial index is $(L_{max}+q)/2-\lfloor (l-(1-q))/2\rfloor$ with $q=L_{max}~\text{mod}\,2$ .

Our fields $\boldsymbol{B}$ , $\boldsymbol{u}_{m}$ and $u_{g}$ are truncated at $L_{B}$ , $L_{u_{m}}$ and $N$ , respectively. We sometimes report a truncation of $L_{max}$ , which means $L_{B}=L_{max}$ , $L_{u_{m}}=2L_{B}$ and $N=L_{max}$ , which is the natural band limit for the fields (recall that $u_{g}/s$ is an even polynomial so that the largest degree is then $2L_{max}-1$ and has the same complexity as $\boldsymbol{u}_{m}$ ), but some calculations have to resort to more aggressively truncated fields.

We work in the subspace of truncated polynomials for all our fields. Little is known about the behaviour of the geostrophic flow in general, and there are examples of $\unicode[STIX]{x1D6FC}$ distributions, albeit non-polynomial ones, for which $u_{g}$ exhibits singularities. We have received a demonstration from Professor G. Ierley, who has completed the analysis of Greenspan (Reference Greenspan1974), showing that the geostrophic flow does not necessarily remain finite at $s=1$ . This is presented in appendix F.

Table 1. The maximum radial index $N_{max}$ as a function of spherical degree $l$ , where $L_{max}$ is the maximum degree of the spherical harmonic expansion.

3.2 Algorithm for finding the ageostrophic flow

In the exact magnetostrophic balance (2.6), in which the magnetic field necessarily satisfies Taylor’s constraint, Taylor showed that $\boldsymbol{u}_{m}$ is uniquely determined up to the geostrophic flow for a given forcing term $\boldsymbol{F}$ (here, the Lorentz force $\unicode[STIX]{x1D735}\times \boldsymbol{B}\times \boldsymbol{B}$ ). Because of our fully spectral representation of $\boldsymbol{B}$ , the forcing term has Cartesian components that are polynomial (of the form $x^{a}y^{b}z^{c}$ ) and are of some maximum degree $P$ . We further observe that if the operator $\hat{\boldsymbol{z}}\,\times$ acts upon a vector with polynomial coefficients, it does not alter the maximum degree. Because the pressure term can be eliminated by taking the curl, it follows that the space of solutions for the ageostrophic flow is spanned by those fully spectral basis functions (defined in § 3.1) which have Cartesian components of degree no greater than $P$ . We may construct a closed system of equations for the unknown poloidal and toroidal modal coefficients by back-projecting onto the same set of velocity modes of degrees up to $P$ ; thus we find the spectral flow coefficients as $\unicode[STIX]{x1D648}^{-1}\,\boldsymbol{f}$ , where

(3.9a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D614}_{i,j}=M(\boldsymbol{u}_{i},\boldsymbol{u}_{j})=\int _{V}\boldsymbol{u}_{i}\boldsymbol{\cdot }(\hat{\boldsymbol{z}}\times \boldsymbol{u}_{j})\,\text{d}V,\quad \boldsymbol{f}_{i}=\int _{V}\boldsymbol{u}_{i}\boldsymbol{\cdot }\boldsymbol{F}\,\text{d}V. & & \displaystyle\end{eqnarray}$$

We note that the pressure term vanishes in the above projection procedure by virtue of the divergence theorem, along with exploiting incompressibility and impenetrability of the boundary.

The results we shall present in §§ 6–8 have very small values of $\unicode[STIX]{x1D712}^{2}$ (as measured by (1.2)), indicating a close adherence to Taylor’s constraint. However, numerically it is impossible to achieve an exactly zero ${\mathcal{T}}$ everywhere in (1.1). It is therefore worth considering the repercussions of this for our algorithm for finding $\boldsymbol{u}_{m}$ , in light of the fact that (3.9) was developed assuming an exact magnetostrophic balance. When ${\mathcal{T}}\neq 0$ for all $s$ , the band-limited property of $\boldsymbol{u}_{m}$ alluded to above is lost. However, the matrix $\unicode[STIX]{x1D648}$ remains invertible and we can nevertheless continue to solve for $\boldsymbol{u}_{m}$ using (3.9), and we do so using the same truncation as before.

To simplify the numerics, we note that the sparse matrix $\unicode[STIX]{x1D648}$ separates into azimuthal modes and by equatorial symmetry (see appendix A). The axisymmetric class containing the geostrophic component needs to be treated carefully, as the matrix $\unicode[STIX]{x1D648}$ is formally singular in this case. One way of eliminating the geostrophic degeneracy is to remove all lowest-radial-order odd- $l$ toroidal modes from the representation. The cylindrical average of the azimuthal component can then be removed from the resulting flow. In all symmetry classes, therefore, this procedure results in the ageostrophic flow with no geostrophic component.

We remark that this construction ensures that the ageostrophic flow carries no angular momentum. To show this we note that to extract the angular momentum from a general flow $\boldsymbol{U}$ one needs to take the inner product

(3.10) $$\begin{eqnarray}\displaystyle \langle \hat{\boldsymbol{t}}_{(1,1,0)}\boldsymbol{\cdot }\boldsymbol{U}\rangle . & & \displaystyle\end{eqnarray}$$

However, no mode in the ageostrophic flow other than $\hat{\boldsymbol{t}}_{(1,1,0)}$ can carry angular momentum. But having removed the geostrophic component from $\boldsymbol{u}_{m}$ , one can see that it carries no angular momentum at all. Only the geostrophic component then contributes to the angular momentum:

(3.11) $$\begin{eqnarray}\displaystyle \int A(s)su_{g}(s)\,\text{d}s & & \displaystyle\end{eqnarray}$$

and we note that in our basis the polynomial ‘ $s$ ’ appearing in (3.11) is just a scalar multiple of $\unicode[STIX]{x1D6EC}_{1}$ , the solid-body component of the flow. Thus we see from (3.8) that all angular momentum is carried by the coefficient $g_{1}$ and therefore it must be set to zero. This is not updated by use of the gradient, whilst all other $N-1$ components are.

3.3 A smoothness-inducing linear friction model

Computationally we face a cascade to smaller wavenumbers at each time step: suppose that at time $t$ the spatial complexity of the magnetic field is $L_{max}$ . The magnetostrophic equation (2.6) is quadratic in $\boldsymbol{B}$ , which creates a $\boldsymbol{u}_{m}$ that is of complexity $2L_{max}$ , whilst the induction term, $\unicode[STIX]{x1D735}\times (\boldsymbol{u}\times \boldsymbol{B})$ , further increases the spatial complexity of $\boldsymbol{B}$ to $3L_{max}$ after a time step. Hence at every time step, we have to truncate the magnetic field at $L_{max}$ and it is this truncated field that is forced to satisfy Taylor’s condition (i.e. the maximum degree of $\boldsymbol{B}$ in (1.1) is $L_{max}$ while the degree of (1.2) is twice this). This truncation error is inevitably injected into the computation of $\boldsymbol{u}_{g}$ . After a large number of computational time steps, the accumulated truncation error inversely cascades into the spectrum of $\boldsymbol{u}_{g}$ and eventually destroys its convergence.

To combat this effect we introduce a linear friction model previously employed by Hollerbach & Ierley (Reference Hollerbach and Ierley1991), drawing on the oceanographic implementation of Salmon (Reference Salmon1986). We shall show the adoption of this idea to be beneficial in suppressing small-scale structures that have little effect on the magnetic field at the time horizon, and that the linear friction will always be below the level of the Lorentz torque ${\mathcal{T}}$ . We define a damping term as ${\mathcal{L}}=-\{\boldsymbol{u}_{g}\}$ (noting that the geostrophic flow vanishes outside of $V$ ), which together with a pseudo-Ekman number $\unicode[STIX]{x1D716}^{2}$ satisfies

(3.12) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D716}\Vert {\mathcal{L}}(\boldsymbol{u}_{g})\Vert \ll \Vert {\mathcal{T}}\Vert \ll 1. & & \displaystyle\end{eqnarray}$$

Correspondingly, in our numerical model, Taylor’s integral ${\mathcal{T}}$ is replaced by ${\mathcal{T}}-\unicode[STIX]{x1D716}\{\boldsymbol{u}_{g}\}$ and a suitably modified version of $\unicode[STIX]{x1D712}^{2}$ is used. We note that there is a link between this approach and that of an Ekman state formulation, in which boundary friction is used to balance Taylor torques. The correspondence is not exact, as we omit the geometric factor containing a singularity at $s=1$ , but $\unicode[STIX]{x1D716}$ can be seen to be akin to $\sqrt{E}$ . This friction term appears in our objective function but none of the underlying physical equations are altered. Thus we still seek purely magnetostrophic states (that is, with zero viscosity); the introduction of the pseudo-Ekman number merely assists us in finding them numerically. In appendix E, we derive the discretized version of the continuous model (E 3) presented in § 2 using this modified version of the target function.

3.4 On the choice of weight function $w$

The choice of the positive weight function $w(s)$ has thus far remained arbitrary. We will report calculations with $w(s)=(1-s^{2})^{\unicode[STIX]{x1D6FE}}$ with three choices of $\unicode[STIX]{x1D6FE}$ , namely $(-1/2,0,1/2)$ . We will show in subsequent examples that the choice of $\unicode[STIX]{x1D6FE}$ has very little effect on the solutions obtained.

There is a need to compute a number of integrals in cylindrical coordinates, such as the value of the Taylor torque (1.1) and its downhill direction (2.32). With our polynomial basis for $\boldsymbol{B}$ and $\boldsymbol{B}^{\dagger }$ , all of these integrals can be computed exactly by using quadrature rules appropriate for the weighting function chosen. For example, $w=1$ dictates that we use Gauss–Legendre quadrature, and $w=(1-s^{2})^{-1/2}$ dictates that we use Gauss–Chebyshev quadrature.

3.5 The use of gradients in the update of $\boldsymbol{u}_{g}$

Note that the target $\unicode[STIX]{x1D712}^{2}$ as defined in (1.2) is a quartic functional with respect to $\boldsymbol{B}$ that can tend to zero as a result of both adherence to Taylor’s constraint, or if $\boldsymbol{B}\rightarrow 0$ . We define a normalized quantity $\mathbb{T}$ , here dubbed the ‘Taylicity’ (we coin the term ‘Taylicity’ to avoid confusion with the term ‘Taylorization’ (Rotvig & Jones Reference Rotvig and Jones2002), which has a different definition and is not a differentiable function of the magnetic field) to quantify the adherence to Taylor’s constraint:

(3.13) $$\begin{eqnarray}\displaystyle \mathbb{T}=\frac{\unicode[STIX]{x1D712}^{2}}{\langle \boldsymbol{B}^{2}\rangle ^{2}}. & & \displaystyle\end{eqnarray}$$

We seek the optimal solution for $\boldsymbol{u}_{g}$ by minimizing the misfit, $\unicode[STIX]{x1D712}^{2}$ , using the following procedure:

1. (i) At $t=0$ , we start with a $\boldsymbol{B}$ field in a Taylor state, calculate $\boldsymbol{u}_{m}$ (assumed constant for $t\in [0,\unicode[STIX]{x0394}t]$ ) by (3.9) and take $\boldsymbol{u}_{g}=0$ ; we then forward evolve (2.7) from $t=0$ to $t=\unicode[STIX]{x0394}t$ .

2. (ii) We compute the Taylicity, $\mathbb{T}$ , defined in (3.13).

3. (iii) We compute the downhill direction, $-\unicode[STIX]{x1D735}_{\boldsymbol{u}_{g}}\unicode[STIX]{x1D712}^{2}$ , with respect to $\boldsymbol{u}_{g}$ using the adjoint method.

4. (iv) We use the limited memory L-BFGS algorithm (Nocedal & Wright Reference Nocedal and Wright2006), a gradient-based optimization method, to calculate an updated $\boldsymbol{u}_{g}$ :

   (3.14) $$\begin{eqnarray}\displaystyle \boldsymbol{u}_{g}^{(i+1)}=\boldsymbol{u}_{g}^{(i)}-\unicode[STIX]{x1D70C}{\mathcal{H}}^{-1}\unicode[STIX]{x1D735}_{\boldsymbol{u}_{g}}\unicode[STIX]{x1D712}^{2}, & & \displaystyle\end{eqnarray}$$
   where $\boldsymbol{u}_{g}^{(i)}$ is the solution of the geostrophic flow at the $i$ th iteration, ${\mathcal{H}}^{-1}$ is the estimate of the inverse Hessian built up by BFGS and based on previous iterations, and $\unicode[STIX]{x1D70C}$ is the step length (Nocedal Reference Nocedal1980). Angular momentum conservation is respected as described in § 3.2.

5. (v) We keep iterating (ii)–(iv) until the normalized misfit is less than or equal to a target or threshold value, $\unicode[STIX]{x1D712}_{T}^{2}$ , i.e., $\mathbb{T}\leqslant \unicode[STIX]{x1D712}_{T}^{2}$ .

6. (vi) When condition (v) is satisfied, the solution of $\boldsymbol{B}$ , $\boldsymbol{u}_{m}$ and $\boldsymbol{u}_{g}$ at $t=\unicode[STIX]{x0394}t$ is accepted as the optimal Taylor state solution and we move on to the next time step.

4.1 Numerical verification of the gradient computation

In order to demonstrate the veracity of the derivations and their numerical implementation, we undertake a comparison of our gradient calculation against the so-called ‘brute-force’ method. We test using the most general structure of the initial magnetic field, which is not a Taylor state, shown in (figure 1 a). We also use an arbitrary ageostrophic flow and geostrophic flow, $\boldsymbol{u}_{g}$ (figure 1 b); this assignment means that (2.6) is not obeyed in this test. No aspect of this brute-force test relies on the adherence to a magnetostrophic balance.

Figure 1. (a) The poloidal and (b) toroidal components of the axisymmetric initial magnetic field used in the gradient benchmark: (a) the stream function is shown, such that field lines are parallel to contours, and (b) the $\unicode[STIX]{x1D719}$ component of the field is shown. (c) The initial geostrophic flow, $\boldsymbol{u}_{g}$ , and (d) Taylor’s integral, ${\mathcal{T}}$ , of the prescribed $\boldsymbol{B}$ as functions of  $s$ .

We forward evolve (2.7) and (2.30)–(2.31) for one time step $\unicode[STIX]{x0394}t=10^{-6}$ , and compute the gradient, $\unicode[STIX]{x1D735}_{\boldsymbol{u}_{g}}\unicode[STIX]{x1D712}^{2}$ , via

1. (i) the adjoint method (2.32) and

2. (ii) a brute-force method in which we perturb each coefficient, $g_{k}$ , of $\boldsymbol{u}_{g}$ , in sequence with $\unicode[STIX]{x1D6FF}g_{k}/g_{k}=10^{-3}$ , and compute the $k$ th entry of the gradient in terms of $\unicode[STIX]{x0394}\unicode[STIX]{x1D712}_{k}^{2}/\unicode[STIX]{x1D6FF}g_{k}$ .

Figure 2 illustrates the gradient comparison for this scenario, comparing $\unicode[STIX]{x1D735}_{\boldsymbol{u}_{g}}\unicode[STIX]{x1D712}^{2}$ of the adjoint method and the brute-force method for each of the spectral coefficients of $\boldsymbol{u}_{g}$ using different weight functions, $w$ . We set the spatial resolution as $L_{B}=30$ , $L_{u_{m}}=60$ and $N_{u_{g}}=30$ . The adjoint method is seen to be accurate for all weight functions.

Figure 2. Comparison of the gradient computations, $\unicode[STIX]{x1D735}_{\boldsymbol{u}_{g}}\unicode[STIX]{x1D712}^{2}$ , using the adjoint method (green dots) and brute-force method (red open squares) for different weight functions $w=(1-s^{2})^{\unicode[STIX]{x1D6FE}}$ for (a) $\unicode[STIX]{x1D6FE}=0$ , (b) $\unicode[STIX]{x1D6FE}=-1/2$ and (c) $\unicode[STIX]{x1D6FE}=1/2$ . In all cases, the green and red symbols overplot. The ordinate shows the index of the spectral coefficient $g_{n}$ of $\boldsymbol{u}_{g}$ . The time step taken is $\unicode[STIX]{x0394}t=10^{-6}$ , and truncation levels are $L_{B}=30$ , $L_{u_{m}}=60$ and $N_{u_{g}}=30$ .

4.2 Benchmark against a Taylor-type analytic solution for $\boldsymbol{u}_{g}$

A second comparison can be made between our methodology and that of Taylor when $\unicode[STIX]{x0394}t$ is small. We look for an instantaneous solution for $\boldsymbol{u}_{m}$ and $\boldsymbol{u}_{g}$ when the field $\boldsymbol{B}$ is prescribed, where

(4.1) $$\begin{eqnarray}\displaystyle \boldsymbol{B}=16\sqrt{\frac{7\unicode[STIX]{x03C0}}{3}}\hat{\boldsymbol{S}}_{(1,2,0)}=\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D735}\times \unicode[STIX]{x1D735}\times [r^{3}(5r^{2}-7)(3\cos ^{2}\unicode[STIX]{x1D703}-1)\hat{\boldsymbol{r}}], & r\in [0,1],\\ -2\unicode[STIX]{x1D735}r^{-3}(3\cos ^{2}\unicode[STIX]{x1D703}-1), & r\in [1,\infty ).\end{array}\right. & & \displaystyle\end{eqnarray}$$

Such a single-mode choice for $\boldsymbol{B}$ is automatically a Taylor state (Livermore, Ierley & Jackson Reference Livermore, Ierley and Jackson2008, Reference Livermore, Ierley and Jackson2009). For this choice of $\boldsymbol{B}$ , the velocity fields $\boldsymbol{u}_{m}$ and $\boldsymbol{u}_{g}$ have closed analytic forms. We solve the magnetostrophic equation (2.6) to obtain $\boldsymbol{u}_{m}$ . To obtain $u_{g}$ , we differentiate (1.3) with respect to time and demand that $\unicode[STIX]{x2202}{\mathcal{T}}/\unicode[STIX]{x2202}t=0$ . This leads to a first-order differential equation (Wu & Roberts Reference Wu and Roberts2015) for $d_{s}(u_{g}/s)$ , where $d_{s}=\text{d}/\text{d}s$ , that has an analytical solution in this case. Note that in the presence of viscosity $\unicode[STIX]{x1D708}$ , the term $\unicode[STIX]{x1D708}sd_{s}(u_{g}/s)$ would be the stress of the geostrophic flow in the tangential direction of the geostrophic cylinder. We make use of this in our diagnostics.

For our chosen single mode $\boldsymbol{B}$ , the instantaneous solution for $\boldsymbol{u}_{m}$ is

(4.2) $$\begin{eqnarray}\displaystyle \boldsymbol{u}_{m}=210r\sin (\unicode[STIX]{x1D703})[22-70r^{2}+39r^{4}-2(r^{2}+7)r^{2}\cos (2\unicode[STIX]{x1D703})-7r^{4}\cos (4\unicode[STIX]{x1D703})]\hat{\unicode[STIX]{x1D753}}, & & \displaystyle\end{eqnarray}$$

and the geostrophic flow has a closed-form solution,

(4.3) $$\begin{eqnarray}\displaystyle u_{g} & = & \displaystyle {\textstyle \frac{5}{4}}s\{\!48(81s^{2}+62)s^{2}-688\log (5(s^{2}-2)s^{2}+6)\nonumber\\ \displaystyle & & \displaystyle -\,2464\sqrt{5}\tan ^{-1}(\sqrt{5}(s^{2}-1))-896\log (s)+4557\!\}+\,c_{1}s.\end{eqnarray}$$

The constant $c_{1}\approx -14282.627018357707$ is the strength of the solid-body rotation component and is chosen so that the angular momentum of the flow in the $\hat{\unicode[STIX]{x1D753}}$ direction remains zero. Within our axisymmetric system, we note that the form of the geostrophic flow contains an $s\log s$ term, which means that the flow is not differentiable at $s=0$ (e.g. Fearn & Proctor Reference Fearn and Proctor1987; Jault Reference Jault1995; Wu & Roberts Reference Wu and Roberts2015). We have chosen not to make any special accommodation for this term, remaining instead with our polynomial description of $u_{g}$ . The polynomial description can, of course, adequately synthesize most continuous functions, as we shall see, although it does not capture the singular derivatives, being itself everywhere regular.

We compare Taylor’s solution to the results of our own method with time step $\unicode[STIX]{x0394}t=10^{-6}$ and spatial and temporal resolutions $L_{B}=30$ , $L_{u_{m}}=60$ and $L_{u_{g}}=30$ . We determine $\boldsymbol{u}_{m}$ using our matrix solver and, in this special case, the solution for $\boldsymbol{u}_{m}$ is entirely toroidal; we find the analytic and the numerical solutions for $\boldsymbol{u}_{m}$ are identical up to machine precision. We initially set the geostrophic flow, $\boldsymbol{u}_{g}$ , to zero as our starting value and forward evolve the induction equation from $t=0$ to $t=\unicode[STIX]{x0394}t$ and obtain a normalized misfit (Taylicity), $\unicode[STIX]{x1D712}_{0}^{2}$ , as a reference value. We perform a number of iterations to gradually refine our estimate of $u_{g}$ . The optimization stops when $\mathbb{T}\leqslant \unicode[STIX]{x1D712}_{T}^{2}$ is satisfied. We use a reduction factor, $R_{d}$ , to quantify the stringentness of the optimization, setting the threshold to be $\unicode[STIX]{x1D712}_{T}^{2}=R_{d}\unicode[STIX]{x1D712}_{0}^{2}$ . Figure 3 compares the numerical and analytic solutions for different weight functions, $w$ , and reduction factors, $R_{d}$ . An accurate solution is obtained if $\unicode[STIX]{x1D712}_{T}^{2}$ is made sufficiently small. For the case $R_{d}=10^{-6}$ , the numerical and the analytic solutions are identical, except for a very small deviation of $\boldsymbol{u}_{g}$ at $s=0$ , due to the $s\log s$ singularity. The three different weight functions, $w$ , appear to perform equally well. As noted in § 3.4, we prefer to use $w=(1-s^{2})^{-1/2}$ as a weight function in order to cancel the geometric factor $\sqrt{1-s^{2}}$ , rendering all important quantities polynomial (e.g. $\boldsymbol{B}_{\unicode[STIX]{x0394}t}^{\dagger }$ , see (E 5)) and adopt this for all subsequent calculations.

Figure 3. Plots of $u_{g}$ and $d_{s}(u_{g}/s)$ as a function of $s$ comparing the numerical and analytic instantaneous solution (4.3) for the single-mode assignment of $\boldsymbol{B}$ given in (4.1). (a–e) The numerical solution (in dashed green) versus the analytic solutions (in red) with different weight functions, $w$ , and reduction factors, $R_{d}$ . (f) Plot of $d_{s}(u_{g}/s)$ as a function of $s$ , where the red curve is the analytic solution and the green, blue, orange and purple ones (these latter three curves almost entirely overlaying one another) correspond to the numerical solutions of $\boldsymbol{u}_{g}$ in (a–d), respectively. Only in such a derivative plot does the $s\log s$ singularity show itself markedly.

6.1 The torsional wave (TW) versus the Taylor state (TS) dynamo model

We now extend our investigations to take into account some dynamics. We compare models in which Taylor’s constraint is adhered to (termed TS models) with models that allow for deviations from this state and consequently generate torsional oscillations or waves (termed TW models), as presaged in § 5 of Taylor (Reference Taylor1963) and rigorously studied by Braginsky (Reference Braginsky1970). Such an approach was first implemented by Jault (Reference Jault1995), although much of our discussion builds upon the study of Roberts & Wu (Reference Roberts and Wu2014), who set out the theoretical links between TS (what they term ordinary Taylor theory, OTT) and TW (what they term MTT, modified Taylor theory).

Our model is the following, in which the inertial force is assumed to be negligible except for the geostrophic component, i.e.

(6.1) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle R_{o}\frac{\unicode[STIX]{x2202}\boldsymbol{u}_{g}}{\unicode[STIX]{x2202}t}+\hat{\boldsymbol{z}}\times \boldsymbol{u}_{m}=(\unicode[STIX]{x1D735}\times \boldsymbol{B})\times \boldsymbol{B}-\unicode[STIX]{x1D735}\unicode[STIX]{x1D70B},\\ \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{B}}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D735}\times [\boldsymbol{u}\times \boldsymbol{B}+{\mathcal{E}}]+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{B},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where ${\mathcal{E}}$ is defined in (5.5), $\boldsymbol{u}=\boldsymbol{u}_{m}+\boldsymbol{u}_{g}$ and $R_{o}$ is the magnetic Rossby number defined in (2.5). We use the $\unicode[STIX]{x1D6FC}_{B}$ -effect in these calculations. The TS model is, of course, defined by $R_{o}=0$ , in which case $\boldsymbol{u}_{g}$ is determined by the methods of § 2.2, whereas for $R_{o}\neq 0$ , $\boldsymbol{u}_{g}$ is explicitly deterministic. We envisage that (6.1) evolves towards an equilibrium state, in which the torsional wave stops and the magnetic field becomes stationary in time. Note that $\unicode[STIX]{x2202}_{t}u_{g}=0$ is true if and only if Taylor’s condition is satisfied, ${\mathcal{T}}=0$ for all $s$ . Therefore, at the equilibrium state, the solutions of the TS and TW models are expected to be identical. It is worth noting, however, that, if the torsional wave speed (Roberts & King Reference Roberts and King2013) is locally relatively small at $s=s_{a}$ , convergence of the TW model as $t\rightarrow \infty$ to the TS model may be slow in the neighbourhood of $s_{a}$ . However, the pointwise vanishing of the Alfvén speed is not, in itself, sufficient to prevent the propagation of torsional waves (Maffei & Jackson Reference Maffei and Jackson2016).

In (6.1), in contrast to what has gone before, the geostrophic component of the flow $\boldsymbol{u}_{g}$ is now determined prognostically and is solved explicitly in time using a second-order Adams–Bashforth method. We solve the TS and TW models for different values of $\unicode[STIX]{x1D6FC}_{0}$ with spatial resolutions $L_{B}=50$ , $L_{u_{m}}=60$ and $N_{u_{g}}=30$ . The stable time step is $\unicode[STIX]{x0394}t=10^{-5.5}{-}10^{-6}$ for the Taylor state model and $\unicode[STIX]{x0394}t\leqslant 10^{-6.5}$ for $R_{o}\leqslant 10^{-3.5}$ , determined empirically through tests.

We set $\boldsymbol{B}=\hat{\boldsymbol{S}}_{(1,1,0)}$ and $\boldsymbol{B}=\hat{\boldsymbol{S}}_{(1,2,0)}$ as the initial condition of the $\unicode[STIX]{x1D6FC}_{B}^{2}$ TS and TW dynamo systems for the dipole and quadrupole symmetries, respectively. (Note that by dint of (3.3) this means $\langle \boldsymbol{B}^{2}\rangle =1$ .) Although we focus on single-parity solutions, it is possible that mixed-parity solutions have a lower onset threshold (Jault Reference Jault1996). This prescription serves to determine the initial $\boldsymbol{u}_{m}$ , after which the whole system is forward evolved in time to equilibrium, the action of the magnetic field on the flow eventually saturating the initial exponential growth. In the quadrupole symmetry, we note that the initial $\boldsymbol{u}_{m}$ is given in (4.2); whereas for the dipolar case, the initial $\boldsymbol{u}_{m}=\mathbf{0}$ because the associated Lorentz force is curl-free. Figure 8 illustrates the Taylicity and the energy evolution of the $\unicode[STIX]{x1D6FC}_{B}^{2}$ dynamo for different values of $\unicode[STIX]{x1D6FC}_{0}$ in the dipole symmetry.

Figure 8. The Braginsky $\unicode[STIX]{x1D6FC}_{B}^{2}$ dynamo in dipole symmetry computed under two paradigms, either Taylor state (TS) or torsional wave (TW). (a) Time evolution of the Taylicity. (b) Energy evolution. The dashed red and green curves show the TW solution with $\unicode[STIX]{x1D6FC}_{0}=18$ , $R_{o}=10^{-3}$ and $10^{-4}$ , respectively. The blue and black curves show the TS solutions for $\unicode[STIX]{x1D6FC}_{0}=18$ and 14.25 with $\unicode[STIX]{x1D712}_{T}^{2}=10^{-7}$ and $\unicode[STIX]{x1D716}=10^{-2}\sqrt{\unicode[STIX]{x1D712}_{T}^{2}}$ .

For $\unicode[STIX]{x1D6FC}_{0}=18$ , we compute the TS and TW model for the purposes of comparison, where $R_{o}=10^{-3}$ and $10^{-4}$ for the TW model. For the TW model, the initial condition of the geostrophic flow is zero. We draw the following conclusions from figure 8:

1. (i) At the initial stage, $t\in [0,0.01]$ (magnetic decay times), the flow experiences a very rapid acceleration due to the rapid variation of the magnetic field, which evolves away from the Taylor manifold. At $t\sim 0.01$ , $\mathbb{T}$ increases from $0$ to $10^{-2}$ . After the rapid acceleration, Taylor’s condition is gradually restored. At $t=1$ , the Taylicity is 100 times smaller than its peak value at $t\approx 0.01$ . The TW simulations are terminated at approximately $t=2.8$ , since the reduction rate of $\unicode[STIX]{x1D712}^{2}$ becomes smaller and smaller and the solution $\boldsymbol{u}_{g}$ is nearly identical to its state at $t=1$ .

2. (ii) As a twin experiment, we evolve the Taylor state model in time with the same initial condition as for the TW model. At each time step, we enforce Taylor’s constraint, i.e., $\mathbb{T}\leqslant \unicode[STIX]{x1D712}_{T}^{2}=10^{-7}$ . This twin experiment shows that the TS model is the asymptotic limit of the TW model as $R_{o}\rightarrow 0$ . As we decrease the inertia, there is a tendency for the Taylicity to be better preserved, i.e. the dashed green line in figure 8(a) lies below the red line, although they do not reach the low values of Taylicity achieved by the TS model. One notes that the equilibrium energies are identical (figure 8 b), while there is no discernible difference in the field structure (e.g. figure 9).

3. (iii) Evolving the TS dynamo system towards the equilibrium state is approximately 3–5 times faster than computing the TW model, for two reasons. Firstly, in the TW model, torsional waves with local frequency proportional to $\{B_{s}^{2}\}^{1/2}$ need to be resolved; these are entirely absent from the TS model. Secondly, we use the solution of $\boldsymbol{u}_{g}$ obtained from the previous time step as a starting approximation for the current time step: it sometimes happens that this is an optimal solution within the tolerance that we seek, i.e., the resulting Lorentz force of the estimated magnetic field satisfies $\mathbb{T}\leqslant \unicode[STIX]{x1D712}_{T}^{2}$ . Thus no iterative improvement for the geostrophic flow is required in such a case. We note that the typical number of iterations required within the TS model is less than three per time step, except for the computations at the initial stages, $t<0.01$ , when 10–20 iterations are required at each time step.

4. (iv) We focus on the behaviour in figure 8(b). At the initial stage, $|{\mathcal{E}}|\ll |\unicode[STIX]{x1D735}\times [(\boldsymbol{u}_{m}+\boldsymbol{u}_{g})\times \boldsymbol{B}]+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{B}|$ and the dynamo systems decay exponentially with the same decay rates for all $\unicode[STIX]{x1D6FC}_{0}$ . The transition stage occurs when the $\unicode[STIX]{x1D6FC}$ -effect becomes significant. For stronger $\unicode[STIX]{x1D6FC}_{0}$ , the transition stage occurs earlier with shorter duration. For example, for $\unicode[STIX]{x1D6FC}_{0}=18$ , the transition occurs at $t=0.15$ and, within one decay time, the dynamo reaches the equilibrium state; while for $\unicode[STIX]{x1D6FC}_{0}=14.25$ (black line), the energy of the dynamo system keeps reducing until $t=0.6$ and slowly increases towards the equilibrium state. The solution finally settles down at approximately $t=4$ .

6.2 Saturation of the $\unicode[STIX]{x1D6FC}_{B}^{2}$ model

Figure 9. Contour plots of the solutions at the equilibrium state of the torsional wave model (a) and the Taylor state model (b) for the $\unicode[STIX]{x1D6FC}_{B}^{2}$ dynamo for $\unicode[STIX]{x1D6FC}_{0}=18$ with the threshold $\unicode[STIX]{x1D712}_{T}^{2}=10^{-7}$ in dipole symmetry. Shown from left to right are meridional field lines ( $\boldsymbol{B}_{FL}$ ), toroidal field $B_{\unicode[STIX]{x1D719}}$ , poloidal (meridional) flow ( $\boldsymbol{u}_{SF}$ ) from the ageostrophic flow $\boldsymbol{u}_{m}$ and the zonal component of $\boldsymbol{u}_{m}$ .

Table 2. The poloidal, toroidal and total energy $\langle \boldsymbol{B}^{2}\rangle /2$ of the $\unicode[STIX]{x1D6FC}_{B}^{2}$ TS dynamo for different $\unicode[STIX]{x1D6FC}_{0}$ in dipole and quadrupole symmetries at equilibrium.

Table 2 shows the variation of the saturated energy of the TS model as we increase $\unicode[STIX]{x1D6FC}_{0}$ for both the dipole and quadrupole symmetry solutions. Note that, in this study, the poloidal magnetic field energy is measured over all space, hence the energy of the poloidal field is slightly larger (approximately 10 % larger) than its energy within the unit sphere. The critical value of the $\unicode[STIX]{x1D6FC}_{B}^{2}$ TS dynamo is slightly larger than the critical value of the kinematic $\unicode[STIX]{x1D6FC}_{B}^{2}$ model (whose $\unicode[STIX]{x1D6FC}_{c}\approx 13.1$ for both the dipole and quadrupole fields). A referee suggested that the magnetic energy bifurcates from $\unicode[STIX]{x1D6FC}_{T}$ as $(\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FC}_{T})^{1/2}$ where $\unicode[STIX]{x1D6FC}_{T}\approx 13.8$ is the critical $\unicode[STIX]{x1D6FC}_{0}$ for a TS model, a value corroborated by Wu & Roberts (Reference Wu and Roberts2015). Table 2 supports this scaling.

At the equilibrium state, the solutions of the TS and the TW model are very similar for the same $\unicode[STIX]{x1D6FC}_{0}$ , as demonstrated by figures 9–11, which show stationary solutions of the TW and TS models for $\unicode[STIX]{x1D6FC}_{0}=18$ in the dipole and quadrupole symmetries. This property is expected and reassuring, serving as an internal check of our methodology. For an axisymmetric dipole magnetic field, the torsional wave speed vanishes at both $s=0$ and $s=1$ , a feature that is likely to locally slow convergence of the TW solution as $t\rightarrow \infty$ ; this may then explain the small difference in the TW and TS solutions at the end points of the $s$ domain at large but finite time.

Figure 10. Contour plots of the solutions at the equilibrium state of the torsional wave model (a) and the Taylor state model (b) for the $\unicode[STIX]{x1D6FC}_{B}^{2}$ dynamo for $\unicode[STIX]{x1D6FC}_{0}=18$ with the threshold $\unicode[STIX]{x1D712}_{T}^{2}=10^{-7}$ in quadrupole symmetry. Shown from left to right are meridional field lines ( $\boldsymbol{B}_{FL}$ ), toroidal field $B_{\unicode[STIX]{x1D719}}$ , poloidal (meridional) flow ( $\boldsymbol{u}_{SF}$ ) from the ageostrophic flow $\boldsymbol{u}_{m}$ and the zonal component of $\boldsymbol{u}_{m}$ .

Figure 11. Plots of $\boldsymbol{u}_{g}$ as a function of cylindrical radius, $s$ , with $\unicode[STIX]{x1D6FC}_{0}=18$ for the $\unicode[STIX]{x1D6FC}_{B}^{2}$ dynamo, where the green curves are the torsional wave solution and the red curves are the solution in the Taylor state model for $\unicode[STIX]{x1D716}=10^{-2}\sqrt{\unicode[STIX]{x1D712}_{T}^{2}}$ with the threshold $\unicode[STIX]{x1D712}_{T}^{2}=10^{-7}$ .

We note that Roberts & Wu (Reference Roberts and Wu2014) (their appendix B) proposed the presence of a singularity in $u_{g}$ at $s=1$ : for the dipole field, $B_{s}$ is always zero at $s=1$ due to symmetry, resulting in Taylor’s condition always being satisfied at the equator and therefore a nugatory role for the geostrophic flow at this point. In figure 11 (and also in the forthcoming figure 14), there is the suggestion of a boundary layer type structure in $u_{g}$ . In concert with Wu & Roberts (Reference Wu and Roberts2015) our solutions remain finite.

6.3 Force and energy balances in the $\unicode[STIX]{x1D6FC}_{B}^{2}$ model

For the $\unicode[STIX]{x1D6FC}_{B}^{2}$ model we now examine the contribution of the different flow components in more detail. Following in the spirit of Taylor (Reference Taylor1963) temporarily, the time-differentiated Taylor integral

(6.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}{\mathcal{T}}=s\int _{C(s)}[\unicode[STIX]{x2202}_{t}\boldsymbol{J}\times \boldsymbol{B}+\boldsymbol{J}\times \unicode[STIX]{x2202}_{t}\boldsymbol{B}]\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D753}}\,\text{d}z\,\text{d}\unicode[STIX]{x1D719}=0 & & \displaystyle\end{eqnarray}$$

can be split into several components as

(6.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}{\mathcal{T}}_{1} & = & \displaystyle s\int _{C(s)}[(\unicode[STIX]{x1D735}\times \unicode[STIX]{x0394}{\mathcal{I}}_{1})\times \boldsymbol{B}+(\unicode[STIX]{x1D735}\times \boldsymbol{B})\times \unicode[STIX]{x0394}{\mathcal{I}}_{1}]\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D753}}\,\text{d}z\,\text{d}\unicode[STIX]{x1D719},\end{eqnarray}$$

(6.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}{\mathcal{T}}_{2} & = & \displaystyle s\int _{C(s)}[(\unicode[STIX]{x1D735}\times \unicode[STIX]{x0394}{\mathcal{I}}_{2})\times \boldsymbol{B}]\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D753}}\,\text{d}z\,\text{d}\unicode[STIX]{x1D719},\end{eqnarray}$$

(6.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}{\mathcal{T}}_{3} & = & \displaystyle s\int _{C(s)}[(\!\unicode[STIX]{x1D735}\times \unicode[STIX]{x1D6FB}^{2}\boldsymbol{B}\times \boldsymbol{B}]\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D753}}+[(\!\unicode[STIX]{x1D735}\times \boldsymbol{B}\times \unicode[STIX]{x1D6FB}^{2}\boldsymbol{B}]\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D753}}\,\text{d}z\,\text{d}\unicode[STIX]{x1D719},\end{eqnarray}$$

which sum to zero at equilibrium, where the inductive secular variation contributions are

(6.6) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\unicode[STIX]{x0394}{\mathcal{I}}_{1}=\unicode[STIX]{x1D735}\times [\boldsymbol{u}_{m}\times \boldsymbol{B}+{\mathcal{E}}],\\ \unicode[STIX]{x0394}{\mathcal{I}}_{2}=\unicode[STIX]{x1D735}\times [\boldsymbol{u}_{g}\times \boldsymbol{B}],\end{array}\right\} & & \displaystyle\end{eqnarray}$$

and we note that there is no second term in (6.4) arising for the geostrophic flow contribution by dint of the fact that it creates only azimuthal secular variation. We analyse these contributions in figure 12(a), which shows the leading-order Lorentz force balance in the $\hat{\unicode[STIX]{x1D753}}$ direction. There is a close balance between the $\boldsymbol{u}_{m}$ and $\boldsymbol{u}_{g}$ contributions, with the diffusion operating on a much longer time scale and contributing relatively little. Figure 12(b) illustrates the resulting Taylor integral of $\boldsymbol{B}$ versus the residue damping term, $\unicode[STIX]{x1D716}{\mathcal{L}}(\boldsymbol{u}_{g})=-\unicode[STIX]{x1D716}A(s)\boldsymbol{u}_{g}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D753}}$ . Here we see that the linear friction term is always small in comparison to the Taylor torque.

Figure 12. The contributions to the time-differentiated Taylor torque at $t=3$ in dipole symmetry for $\unicode[STIX]{x1D6FC}_{0}=18$ with $\unicode[STIX]{x1D712}_{T}^{2}=10^{-7}$ and $\unicode[STIX]{x1D716}=10^{-2}\sqrt{\unicode[STIX]{x1D712}_{T}^{2}}$ . (a) The leading-order balance in Taylor’s integral from the magnetic wind $\unicode[STIX]{x0394}{\mathcal{T}}_{1}$ (in red) and the geostrophic flow $\unicode[STIX]{x0394}{\mathcal{T}}_{2}$ (in green), with the diffusive term $\unicode[STIX]{x0394}{\mathcal{T}}_{3}$ (in black) making a minor contribution. (b) The normalized Taylor’s integral, ${\mathcal{T}}(s)/\langle \boldsymbol{B}^{2}\rangle$ , in blue alongside the damping term $\unicode[STIX]{x1D716}{\mathcal{L}}(\boldsymbol{u}_{g})$ in red.

Another useful diagnostic is the energy balance:

(6.7) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}\int _{V+\hat{V}}\frac{1}{2}\boldsymbol{B}^{2}\,\text{d}V=\int _{V+\hat{V}}\boldsymbol{B}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\boldsymbol{B}}{\unicode[STIX]{x2202}t}\,\text{d}V=\int _{V+\hat{V}}\boldsymbol{B}\boldsymbol{\cdot }[\unicode[STIX]{x1D735}\times (\boldsymbol{u}\times \boldsymbol{B}+{\mathcal{E}})-\unicode[STIX]{x1D702}^{\prime }\unicode[STIX]{x1D735}\times \unicode[STIX]{x1D735}\times \boldsymbol{B}]\,\text{d}V, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

and we define the following rate-of-work terms:

(6.8) $$\begin{eqnarray}\displaystyle W_{\unicode[STIX]{x1D6FC}_{B}} & = & \displaystyle \int _{V+\hat{V}}\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D6FC}_{0}\unicode[STIX]{x1D6FC}_{B}\boldsymbol{B})\,\text{d}V,\end{eqnarray}$$

(6.9) $$\begin{eqnarray}\displaystyle W_{\unicode[STIX]{x1D6FB}^{2}\boldsymbol{B}} & = & \displaystyle \int _{V+\hat{V}}\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D6FB}^{2}\boldsymbol{B}\,\text{d}V,\end{eqnarray}$$

(6.10) $$\begin{eqnarray}\displaystyle W_{\boldsymbol{u}_{m}} & = & \displaystyle \int _{V+\hat{V}}\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times (\boldsymbol{u}_{m}\times \boldsymbol{B})\,\text{d}V,\end{eqnarray}$$

(6.11) $$\begin{eqnarray}\displaystyle W_{\boldsymbol{u}_{g}} & = & \displaystyle \int _{V+\hat{V}}\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times (\boldsymbol{u}_{g}\times \boldsymbol{B})\,\text{d}V.\end{eqnarray}$$

An analysis of the equilibrium state for $\unicode[STIX]{x1D712}_{T}^{2}=10^{-7}$ , $\unicode[STIX]{x1D6FC}_{0}=18$ and $\unicode[STIX]{x1D716}=10^{-2}\sqrt{\unicode[STIX]{x1D712}_{T}^{2}}$ yields $W_{\unicode[STIX]{x1D6FC}_{B}}=-W_{\unicode[STIX]{x1D6FB}^{2}\boldsymbol{B}}\approx 80$ in dipole parity, whereas $W_{\unicode[STIX]{x1D6FC}_{B}}=-W_{\unicode[STIX]{x1D6FB}^{2}\boldsymbol{B}}\approx 90$ in quadrupole parity, while the contributions of $\boldsymbol{u}_{m}$ and $\boldsymbol{u}_{g}$ are expected to be zero and indeed are very small in the computation, $W_{\boldsymbol{u}_{m}}$ and $W_{\boldsymbol{u}_{g}}$ both being ${\sim}O(10^{-4}){-}O(10^{-5})$ .

6.4 Role of linear friction

In order to obtain a smooth and well-converged solution, a damping term has been applied in the computation. We compare the $\unicode[STIX]{x1D6FC}_{B}^{2}$ TS dynamo of $\unicode[STIX]{x1D6FC}_{0}=18$ with and without the damping term ${\mathcal{L}}(\boldsymbol{u}_{g})$ in a quadrupole symmetry class computation, where $\boldsymbol{u}_{g}$ is well defined everywhere. Both systems evolve with an identical energy trajectory and saturate at the same equilibrium state with the magnetic energy, $E=\langle \boldsymbol{B}^{2}\rangle /2=0.427$ (as in table 2); however there are differences in the solution of the geostrophic flow. Figure 13 illustrates the solution for $\boldsymbol{u}_{g}$ with (in red) and without (in green) the damping term, $\unicode[STIX]{x1D716}{\mathcal{L}}(\boldsymbol{u}_{g})=-\unicode[STIX]{x1D716}A(s)\boldsymbol{u}_{g}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D753}}$ . For the large-scale component of the geostrophic flow, the two solutions are in agreement; however, for the solution without the residue absorbing term, numerical error accumulates at the small scale and results an oscillatory solution. A very small residue damping term significantly improves the convergence and smoothness of the solution for the geostrophic flow.

Figure 13. The $\unicode[STIX]{x1D6FC}_{B}^{2}$ TS dynamo solutions of $\boldsymbol{u}_{g}$ as a function of $s$ for $\unicode[STIX]{x1D6FC}_{0}=18$ in Qp with (in red) and without (in green) the damping term $\unicode[STIX]{x1D716}{\mathcal{L}}(\boldsymbol{u}_{g})$ , where, for both solutions, the threshold value of the Taylor state is $\unicode[STIX]{x1D712}_{T}^{2}=10^{-7}$ and $\unicode[STIX]{x1D716}=10^{-2}\sqrt{\unicode[STIX]{x1D712}^{2}}$ .

Figure 14. The solution for $\boldsymbol{u}_{g}$ of the $\unicode[STIX]{x1D6FC}_{B}^{2}$ TS dynamo for different $\unicode[STIX]{x1D6FC}_{0}$ at equilibrium, where the black, red and green curves stand for the cases of $\unicode[STIX]{x1D6FC}_{0}=14.25$ , 16 and 18 in dipole symmetry (a) and for $\unicode[STIX]{x1D6FC}_{0}=14$ , 16 and 18 in quadrupole symmetry (b).

Figure 15. The total azimuthal flow for the $\unicode[STIX]{x1D6FC}_{B}^{2}$ TS dynamo at equilibrium for (a) dipole symmetry with $\unicode[STIX]{x1D6FC}_{0}=14.25$ , 16 and 18 and (b) for quadrupole symmetry with $\unicode[STIX]{x1D6FC}_{0}=14$ , 16 and 18.

6.5 Dominance of $u_{g}$

We conclude this section by focusing on the azimuthal flow in the steady saturated solutions, noting that the geostrophic component is overwhelmingly dominant. We plot the solutions for $\boldsymbol{u}_{g}$ as a function of $\unicode[STIX]{x1D6FC}_{0}$ in figure 14 and the total azimuthal flow $u_{\unicode[STIX]{x1D719}}=[u_{m}]_{\unicode[STIX]{x1D719}}+u_{g}$ in figure 15 in the equilibrium states for the cases of $\unicode[STIX]{x1D6FC}_{0}=14.25,~16$ and 18 in the dipole symmetry and for $\unicode[STIX]{x1D6FC}_{0}=14,~16$ and 18 in the quadrupole symmetry. For all cases, we set $\unicode[STIX]{x1D712}_{T}^{2}=10^{-7}$ and $\unicode[STIX]{x1D716}=10^{-2}\sqrt{\unicode[STIX]{x1D712}_{T}^{2}}$ . For the $\unicode[STIX]{x1D6FC}_{B}^{2}$ dynamo, the flow is dominated by the geostrophic component. As we increase the driving force, $\unicode[STIX]{x1D6FC}_{0}$ , the intensity of the magnetostrophic and the geostrophic flow increases, although the flow pattern remains essentially the same. Table 3 records the root-mean-square values of the poloidal, toroidal and total flows. All solutions demonstrate spectral convergence, as shown in figure 16 for the energy spectra of $\boldsymbol{B}$ and $\boldsymbol{u}=\boldsymbol{u}_{m}+\boldsymbol{u}_{g}$ ; more than eight decades of convergence in energy are found in the poloidal components, with slightly less in the toroidal components.

Figure 16. The poloidal and toroidal energy of $\boldsymbol{B}$ and $\boldsymbol{u}=\boldsymbol{u}_{m}+\boldsymbol{u}_{g}$ for the $\unicode[STIX]{x1D6FC}_{B}^{2}$ TS dynamo with $\unicode[STIX]{x1D6FC}_{0}=18$ in dipole (a–d) and quadrupole (e–h) symmetry classes. The triangles stand for the solution of $\unicode[STIX]{x1D712}_{T}^{2}=10^{-7}$ and $\unicode[STIX]{x1D716}=10^{-2}\sqrt{\unicode[STIX]{x1D712}_{T}^{2}}$ for both the Dp and Qp; and the diamonds are for $\unicode[STIX]{x1D712}_{T}^{2}=10^{-8}$ and $\unicode[STIX]{x1D716}=10^{-2}\sqrt{\unicode[STIX]{x1D712}_{T}^{2}}$ in dipole symmetry and for $\unicode[STIX]{x1D712}_{T}^{2}=10^{-9}$ and $\unicode[STIX]{x1D716}=10^{-2}\sqrt{\unicode[STIX]{x1D712}_{T}^{2}}$ in quadrupole symmetry.

Table 3. The $L_{2}$ -norm of $\boldsymbol{u}_{m}$ of the poloidal and the toroidal components, $\sqrt{\int _{V}[\boldsymbol{u}_{m}]_{P}^{2}\,\text{d}V}$ and $\sqrt{\int _{V}[\boldsymbol{u}_{m}]_{T}^{2}\,\text{d}V}$ , and the norm of $\boldsymbol{u}_{g}$ , $\sqrt{\int _{V}\boldsymbol{ u}_{g}^{2}\,\text{d}V}$ , of the $\unicode[STIX]{x1D6FC}_{B}^{2}$ TS dynamo for different $\unicode[STIX]{x1D6FC}_{0}$ in dipole and quadrupole symmetries.