Determination of the instantaneous geostrophic flow within the three-dimensional magnetostrophic regime

1. Introduction

Earth's magnetic field is generated by a self-excited dynamo process through the flow of electrically conducting fluid in the outer core. Although the set of equations that govern this process are known, their numerical solution is challenging because of the extreme dynamical conditions [1]. Of particular note is the extreme smallness of the core's estimated viscosity, and the large disparity between the daily timescale associated with Earth's rotation and the thousand-year timescale that governs the long-term geomagnetic evolution. Represented in terms of non-dimensional numbers, this means that the Rossby number (also known as the magnetic Ekman number, E_η, measuring the ratio of rotational to magnetic timescales) is R_o∼10⁻⁹ and the Ekman number (measuring the ratio of rotational to viscous effects) is E∼10⁻¹⁵. The smallness of these parameters means that rapid (sub-year) timescales associated with inertial effects (e.g. torsional waves) and extremely thin boundary layers (of depth about 1 m) must be resolved in any Earth-like numerical model, even though neither likely plays an important role in the long-term evolution of the geodynamo.

Over the past decades, modellers of the long-term geomagnetic field have followed one of two largely independent strategies in order to circumvent these problems. First, beginning with the work of [2,3], it was noted that by artificially increasing these two parameters by many orders of magnitude to now typical values of R_o = 10⁻³, E = 10⁻⁷ [4], the numerically difficult rapid timescales and short length scales are smoothed, allowing larger time steps, and therefore ultimately permitting a longer time period to be studied for a given finite computer resource. Although such (now mainstream) models can reproduce many characteristics of Earth's geomagnetic field, several studies have cast doubt as to whether they obey the correct force balance within the core [1,5,6], although some evidence points to models being on the cusp of faithfully representing Earth's core [7,8].

In the second strategy, which we consider here in this paper, the values of R_o and E are both set to zero [9]. By entirely neglecting inertia and viscosity, the challenging aspects of rapid timescales and very short viscous lengthscales are removed and this approximation will likely lead to a computationally less demanding set of equations to solve. The resulting dimensionless magnetostrophic regime then involves an exact balance between the Coriolis force, pressure, buoyancy and the Lorentz force associated with the magnetic field B itself:

where F_B is a buoyancy term that acts in the unit radial direction $\hat{\mathit{r}}$ and $\hat{\mathit{z}}$ is the unit vector parallel to the rotation axis [10]. In a full sphere (neglecting the solid inner core), a complete description of the geodynamo requires a solution of (1.1) alongside equations describing the evolution of B and F_B within the core, whose boundary conditions derive from the surrounding electrically insulating impenetrable overlying mantle. Denoting (s, ϕ, z) as cylindrical coordinates, Taylor [9] showed that, as a consequence of this magnetostrophic balance, the magnetic field must obey at all times t the well-known condition for any geostrophic cylinder C(s) of radius s coaxial with the rotation axis.

Taylor also showed that it is expedient to partition the magnetostrophic flow of (1.1), using a cylindrical average, into geostrophic and ageostrophic parts:

in which the ageostrophic flow u_a has an azimuthal component with zero cylindrical average. Provided equation (1.2) is satisfied, equation (1.1) can be used to find u_a directly (for example, by using Taylor's constructive method or the integral method of [1]), although the geostrophic flow remains formally unspecified by (1.1). As Taylor further showed however, the geostrophic flow can be constrained by insisting that equation (1.2) is not just satisfied instantaneously but for all time. The task of u_g is then to keep the magnetic field on the manifold of Taylor states [11]. It is noteworthy that, in such a model, at all times the flow is enslaved to B and F_B.

In his 1963 paper, Taylor showed that (for a fully three-dimensional system) the geostrophic flow was at every instant the solution of a certain second-order differential equation (ODE) whose coefficients depend on B and F_B. His elegant and succinct analysis has been reproduced many times in the literature. It may then come at some surprise that in the intervening five decades there have been no published implementations of his method (that the authors are aware of). Very likely, this is due to a subtle issue concerning the treatment of the magnetic boundary conditions. As we shall show, rather than being applicable to a general (Taylor state) B, Taylor's method is only valid for a small subset of Taylor states. Of crucial importance is that this subset does not include those states likely to be realized in any analytical example or in any practical numerical scheme to solve the magnetostrophic equations. The main goal of this paper is to describe why this happens, and to modify Taylor's method in order that it can apply more generally.

Despite the lack of headway using a direct application of Taylor's ODE, some alternative methods to evolve the magnetostrophic equation have shown success. By treating a version of the Taylor integral (1.2) that is specific to axisymmetry [12,13], Wu & Roberts [14] demonstrated that they could evolve the magnetostrophic system by solving a first-order differential equation for the geostrophic flow, rendering the Taylor integral zero to first order, and went on to apply it to a variety of examples. In an independent line of investigation [15] showed that, by using control theory, it is possible to find u_g implicitly such that the Taylor integral is zero at the end of any finite timestep. As we show later explicitly by example, their method is fundamentally three dimensional, although in their paper they only applied it to the axisymmetric case. The generalized version of Taylor's method that we present in this paper is also fully three dimensional and provides an alternative means to that of Li et al. [15] of calculating the geostrophic flow. Either of these methods may provide a route to create a fully three-dimensional magnetostrophic alternative to the mainstream numerical models with weak viscosity and inertia. We note however, that the methods we describe within this paper are restricted to the full sphere, and we do not attempt to incorporate the inner core or any of its dynamical effects.

An alternative route to finding a magnetostrophic dynamo is to reinstate viscosity and inertia and investigate the limit as both E and R_o become small [13]. In such models, it is important that the Lehnert number, λ (estimated to be 10⁻⁴ in Earth's core) is small in order that inertial modes separate from magneto-Coriolis waves and can be filtered out [16]. It is also worth noting that since a significant part of the Coriolis term may be balanced by the pressure gradient (e.g. [17]), the simple estimates of the Rossby number reported above may be too small. A different estimate of importance of inertia is the Alfvén number (measuring the square root of the ratio of kinetic to magnetic energies), whose small value of A∼10⁻² still supports neglecting the inertial term although with weaker justification [8]. Arguably retaining inertia and viscosity would result in models closer to geophysical reality than those that are purely magnetostrophic as this is precisely the regime of the Earth's core. A variety of studies reported evidence of behaviour independent of E in the inviscid Taylor-state limit, either from a direct solution [18], or from solving the equations assuming asymptotically small E [19,20]. To date, all models of this type have been axisymmetric and there have been no attempts at a general three-dimensional implementation of these ideas. One difficulty with treating asymptotically small E is that the resulting equation for u_g is an extremely delicate ratio of two small terms, whose form is dependent on the specific choice of mechanical boundary conditions [21]. The convergence of magnetostrophic and asymptotically low-E models remains an outstanding question.

The remainder of this paper is structured as follows. Before we can explain why Taylor's method of determining the geostrophic flow fails in general, we need to set out some general background and review other alternative schemes: this is accomplished in §§2–5. In §6, we discuss the importance of a key boundary term and why it restricts the validity of Taylor's method; we then show explicitly in a simple case that Taylor's method fails. In §§7–9, we generalize Taylor's method and give some examples, discussing the existence of weak singularities in §10 we end with a discussion in §11.

2. General considerations

(a) Non-dimensionalization

In the non-dimensionalization considered in this paper, length is scaled by L, the outer core radius 3.5 × 10⁶ m, time by the Ohmic diffusion time τ (250–540 kyr) [22] and speed by $\mathcal{U}=L{\tau }^{-1}\approx 5\times {10}^{-7}$ m s⁻¹. The scale used for the magnetic field is $\mathcal{B}={(2{\mathrm{\Omega }}_0{\mu }_0{\rho }_0\eta )}^{1/2}$ [10], where for Earth the physical parameters take the following values: angular velocity Ω₀ = 7.272 × 10⁻⁵ s⁻¹, permeability μ₀ = 4π × 10⁻⁷ NA⁻², density ρ₀ = 10⁴ kg m⁻³ and magnetic diffusivity η = 0.6–1.6 m² s⁻¹. These parameters lead to the non-dimensional parameters R_o = η/(2ΩL²) ≈ 10⁻⁹ and E = ν/(2ΩL²) ≈ 10⁻¹⁵, whose small values motivate neglecting the terms they multiply.

The value of $\mathcal{B}\approx 1.7\hspace{0.17em}\mathrm{mT}$ is close to the estimate of the geomagnetic field strength of Gillet et al. [23], and so we use dimensionless magnetic fields with toroidal or poloidal components of RMS (root mean squared) strength of unity. This corresponds to a dimensional RMS magnitude of 1.7 mT for purely toroidal or purely poloidal fields and $1.7\sqrt{2}\approx 2.4\hspace{0.17em}\mathrm{mT}$ for mixed states. Using $\mathcal{U}$, this choice enables the immediate interpretation of the dimensional scale of any flow that we show.

(b) Magnetic field representation and the initial state

In our full sphere of unit radius, the position r is naturally described in spherical coordinates (r, θ, ϕ), although the importance of the rotation axis also leads us to use cylindrical coordinates (s, ϕ, z). The magnetic field B can be written using a toroidal (T)-poloidal (S) framework

with S and T expanded as where Y^m_l is a spherical harmonic of degree l and order m. The functions $\mathcal{S}$ and $\mathcal{T}$ must be chosen to satisfy both Taylor's condition (1.2), along with the electrically insulating boundary conditions at r = 1 that can be written The fluid is assumed to be incompressible and hence the flow u can also be written in a comparable form, and due to the absence of viscosity only satisfies an impenetrability condition: u_r = 0 on r = 1. We cannot impose no-slip or stress-free conditions, there being no boundary layer to accommodate any adjustment from the free-stream inviscid structure.

All time-dependent magnetostrophic models, axisymmetric or three dimensional, require an initial state from which the system evolves. Because the flow is defined completely by the magnetic field and F_B, only the initial structure of the magnetic field B(0) and F_B(0) are needed: there is no need to specify the initial flow. A general scheme for finding an exact initial Taylor state using a poloidal–toroidal representation was described by Livermore et al. [24]; in general, it requires a highly specialized magnetic field to render its integrated azimuthal Lorentz force zero over all geostrophic cylinders. However, in a full sphere such cancellation can be achieved in a simple way by exploiting reflectional symmetry in the equator [25]. Using the Galerkin basis of single-spherical-harmonic modes that satisfy the boundary conditions (see appendix Aa), suitable simple modal expansions are automatically Taylor states.

(c) Overview of time evolution

Because of the absence of inertia, at each instant the magnetostrophic flow is entirely determined by B and F_B from equation (1.1): therefore the system, as a whole, only evolves through time-evolution of the quantities F_B and B. The evolution of F_B is assumed to be tractable and lies outside the scope of this study: for simplicity we shall henceforth assume that F_B = 0, although we note that all the methods nevertheless apply in the case of non-zero F_B. The evolution of the magnetic field is described by the induction equation:

where η≠0 is the magnetic diffusivity (assumed constant) and ∂_t = ∂/∂t. Assuming that we can evolve B and F_B (using standard methods), the major outstanding task is then to determine the flow at any instant given B and F_B.

The ageostrophic component of the flow, containing all the (possibly complex) axially asymmetric structure turns out to be straight-forward to calculate, as it can be determined either through the integral method of Roberts & King [1], the constructive method of Taylor [9] or the spectral method as described in appendix Ac. By contrast, the more elementary geostrophic flow, depending only on s, is surprisingly difficult to compute, owing to its key role of maintaining Taylor's constraint.

There are two ways in which the geostrophic flow may be found, which differ in philosophy. In the first, we may undertake an instantaneous analysis to find the geostrophic flow that gives zero rate of change of Taylor's constraint: ∂_tT(s, t) = 0 [9]. Because of the resulting closed-form analytic description, such methods can be useful in computing snapshot solutions that elucidate the mathematical structure of the geostrophic flow, for example, the presence of any singularities. However, as a practical time-evolution tool, their utility is not so obvious. For example, the simple explicit time-evolution scheme, defined by assuming an instantaneous solution is constant over a finite time interval, would lead to a rapid divergence from the Taylor manifold (see [11] for an example).

In the second type of method, we may consider taking a time step (of size h), determining the geostrophic flow implicitly by the condition that the magnetic field B(t + h) satisfies Taylor's constraint [14,15]. In general, implicit and instantaneous methods will only produce the same geostrophic flow in a steady state, or for a time-dependent state for infinitesimally small h.

All methods to determine the geostrophic flow do so up to an arbitrary solid body rotation: u_g = as. The constant a can be found through requiring zero global angular momentum

where $Z_T=\sqrt{1-s^2}$ is the half-height of C(s). We also assume the geostrophic flow is everywhere finite (although we permit singularities in the higher-order derivatives), which is implemented by additional conditions where necessary.

3. Braginsky's formulation

Before discussing the determination of the geostrophic flow in more detail, we briefly review a crucial alternative formulation of Taylor's constraint due to [12], which laid the foundations of many subsequent works on the subject (e.g. [14,26–29]). As an identity the Taylor integral (1.2) can be equivalently written as

where N and S are the northern and southern end caps, respectively, of the cylinder C(s) at the intersection with the spherical boundary at r = 1.

It is also useful to consider the net magnetic torque on all fluid enclosed within C(s), Γ_z, defined by

In our full-sphere geometry, it is clear that Γ_z(s, t) is zero if and only if T(s, t) is zero, although in a spherical shell it is possible that a piecewise (non-zero) solution exists for Γ_z. The condition Γ_z = 0 defines what we refer to as the Braginsky constraint: which is equivalent to Taylor's constraint, and simplifies for specific classes of magnetic fields that cause the boundary term to vanish. There are two such classes of magnetic fields which have no sources in the exterior of r = 1: fields with no radial component on r = 1 (e.g. toroidal fields) and fields that have a vanishing azimuthal component on r = 1 (e.g. axisymmetric fields).

It is important to note the significant difference in the mathematical structure between the constraints of Braginsky (3.2) and Taylor (1.2). In (3.2), there is a clear partition between the two surface integral terms on the right-hand side: the first term is an integral defined over C(s) that is independent of the magnetic field values on the end caps (these being a set of measure zero); the second end-cap term depends only on the boundary values of the magnetic field. By contrast, although ostensibly Taylor's integral (1.2) is an integral over the surface C(s), the integrand involves a spatial derivative (the curl of B) leading to a dependence on the boundary values of the magnetic field. As we will see later, this hidden dependence on the boundary conditions has a deep consequence on Taylor's method for determining the geostrophic flow.

4. Existing methods to determine the geostrophic flow

Our modification of Taylor's method described in §7 determines the instantaneous geostrophic flow in a fully three-dimensional geometry. In this section, we briefly review other methods available to calculate the geostrophic flow whose working assumptions are different: either they are axisymmetric or designed to take a finite time step and are not instantaneous. Where there is overlap in applicability, we will use these methods to numerically confirm our solutions.

(a) An axisymmetric first-order implicit method

As noted above, under axisymmetry Braginsky's condition collapses to

This simple form was exploited by Wu & Roberts [14] who considered taking a single timestep of duration h, after which they required

The left-hand side here approximates Γ_z(s, t + h), so this ensures that (4.1) is satisfied to first order. To find an equation for the geostrophic flow they differentiated equation (4.1) with respect to time and used the fact that the geostrophic term in the induction equation reduces to

They obtained the following first-order ordinary differential equation describing the geostrophic flow where and The subscripts of zero denote a restriction to axisymmetry of (more general) three-dimensional quantities that are defined subsequently. Wu & Roberts [14] implemented this method by solving equation (4.4) using a finite difference scheme. It is worth remarking that this scheme allows small numerical deviations from a Taylor state (since (4.2) is only approximate). Because the method depends upon (4.1) which is tied to axisymmetry, their method is not extendable to three dimensions.

(b) A three-dimensional fully implicit scheme

An alternative implicit scheme proposed by Li et al. [15] was to seek a geostrophic flow that ensured Taylor's constraint is satisfied (without error) in a numerical scheme after taking a single timestep h. By extending to multiple timesteps, this method is suitable to describe fully three-dimensional time-dependent dynamics. Although the authors only demonstrated its utility on axisymmetric examples, in this paper we will show how the method applies to three dimensions with a single short time-step.

The key idea is to minimize (hopefully to zero) the target function

by optimizing over all possible choices of u_g, assumed constant throughout the interval 0≤t≤h. Although [15] set out a sophisticated algorithm to do this in general based on control theory, here we describe a simplification of the method which is suitable for h≪1, which we can use to benchmark our instantaneous solutions of the generalized three-dimensional Taylor methodology.

Like Li et al. [15] we adopt a modal expansion of u_g, of which a general form is

where T_i(2s² − 1) are even Chebyshev polynomials of the first kind, and we allow a weak logarithmic singularity at the origin as required by our analytic results in §6a; see also §10.

Because we plan to take only a single time step of size h≪1, we adopt a very simple first-order explicit Euler time evolution scheme

which is then substituted into (4.6). For simplicity, we assume that the ageostrophic flow, calculated at t = 0, is also constant over the time-step. As a representation of the magnetic field (and its rate of change), we use a Galerkin scheme (see appendix Aa), which satisfies the boundary conditions (2.1) automatically. Practically, this means that we use $\overline{\mathcal{I}}$ (see equation (2.2)) in place of ∂_tB, where the overbar denotes the projection onto the Galerkin basis. The coefficients A_i and B are then found through minimizing Φ. We note that since B(t + h) is formally linear in u_g(s), T(s, t + h) is then quadratic and hence Φ quartic in the coefficients A_i and B. Li et al. [15] found the minimum using an iterative scheme, although we note that, in general (and without a good starting approximation), finding such a minimum may be problematic.

It is noteworthy, however, that in the axisymmetric case this analysis is greatly simplified. Through equation (4.3) only the azimuthal component of B(t + h) depends on u_g, and equation (3.1) shows that T(s) is now linear and Φ quadratic in u_g, hence finding the minimum of Φ is more straightforward.

(c) An instantaneous axisymmetric method

Wu & Roberts [14] also presented a method for finding an instantaneous solution for the geostrophic flow in axisymmetry. Through differentiating with respect to time equation (4.1) they arrive at the following first-order ODE, here referred to as the BWR¹ (Braginsky–Wu–Roberts) equation:

which is the same as (4.4) without the final term. This gives u_g(s) explicitly as In all the cases we consider, (4.9) can be solved analytically (with the assistance of computer algebra). A further property of this equation is that, for a purely poloidal axisymmetric magnetic field, the solution u_g is independent of the magnetic diffusivity η. This is because ∇²B is also purely poloidal and a purely poloidal field has no azimuthal component. Thus, and the diffusion term (within S₀) then never appears in (4.8). This differs from the case of a more general magnetic field with both toroidal and poloidal components, where u_g depends upon η.

We also observe that for an axisymmetric purely toroidal field, since B_s = 0 everywhere equation (4.8) is null because α₀ = S₀ = 0 reducing it to the tautology 0 = 0 and hence placing no constraint on the geostrophic flow.

(d) Taylor's three-dimensional instantaneous method

We end this section by discussing the well-known (instantaneous) method of Taylor, who determined the unknown geostrophic flow by differentiating with respect to time (denoted by the over-dot shorthand) the Taylor integral in equation (1.2) to produce:

On substituting directly for $\dot{\mathbf{B}}$ from equation (2.2) in addition to its curl (describing $\mathbf{\nabla }\times \dot{\mathbf{B}}$), Taylor showed that for fully three-dimensional Taylor states B the resulting equation for the geostrophic flow can be written in a remarkably succinct form as the second-order ordinary differential equation In the above, the coefficients are and G(s) is a function describing the interaction of u_a and the magnetic field defined as Note the mistake in Taylor [9] where a factor of s is omitted within the coefficient β. The functions α₀ and S₀, previously defined, are simply axisymmetric variants of α given above and S(s) defined as where C^a is as defined in equation (4.5). The fact that the coefficients α(s) and β(s) are spatially dependent means that analytic solutions to (4.11) are very rare and in general only numerical solutions are possible. Of crucial note is that the boundary conditions played no part in the derivation above.

5. Technical aside: higher-order boundary conditions

(a) Higher-order boundary conditions in the heat equation

Taylor's method is based on the instantaneous evolution (which we can take to be at time t = 0) of the magnetostrophic system whose magnetic field is prescribed and must satisfy Taylor's constraint. Here, we discuss higher-order boundary conditions, the importance of which has so far been overlooked. We start by introducing this concept in a simple PDE, then we discuss the relevance for Taylor's equation.

Suppose we are interested in finding f(x, t) on x∈[0, 1], whose evolution is described by the heat equation in the interior of the domain

to be solved with the boundary conditions f(0, t) = f(1, t) = 0. For this simple equation, the general solution can be written in the form Let us now suppose we have an initial state: which satisfies the boundary conditions. Its future evolution would be given by the projection onto the normal modes as above.

In Taylor's analysis, part of the integral in (1.2) could be converted to a boundary term. Here we consider an analogy which is exactly integrable:

using the boundary conditions. In Taylor's derivation, he differentiated under the integral sign and substituted directly for ∂f/∂t, in order to find the equation that u_g must satisfy using an instantaneous initial magnetic field. In our example, this produces

At t = 0, we evaluate the above expression as −6 (note that f_xxx(x, 0) =  − 6) resulting in an apparent contradiction with (5.1) and illustrating that this approach is not generally valid.

The problem arises because the initial state does not satisfy the condition f_xx(0, t) = f_xx(1, t) = 0, which arises from differentiating f(0, t) = f(1, t) = 0 with respect to time and substituting the PDE. The condition f_xx(0, t) = f_xx(1, t) = 0 is called the first-order boundary condition [30]. The consequence of the initial state not satisfying the first-order boundary condition is that the solution is not smooth at the boundary at t = 0. Specifically, the derivatives in (5.2) do not exist and thus the above derivation is not valid. As a simple illustration of the issue, note that the general solution implies that $f_{xxx}(x,0)=-\sum _n{n}^3{\pi }^3{A}_n\cos (n\pi x)$, which cannot represent the constant function f_xxx(x, 0) =  − 6 associated with the initial state. This lack of smoothness only occurs at the initial time t = 0. At any later time (t > 0), the solution is infinitely smooth; this is the smoothing property of the heat equation.

In the very special case that the initial state satisfies the first-order boundary conditions (e.g. f(x, 0) = x³(1 − x)³) then there is no contradiction and (5.1) and (5.2) are consistent. However, for a general initial condition, the procedure adopted is not valid.

(b) The relevance for Taylor's equation

We now discuss the relevance of the above discussion of higher-order boundary conditions in the context of the Earth's magnetic field. In the derivation of Taylor's second-order ODE (4.11), it is implicitly assumed that B and all its time derivatives are (initially) smooth everywhere. Although it is somewhat hidden in Taylor's original derivation, taking the time-derivative of the equivalent form of (3.1) makes this explicit:

Taylor substituted everywhere the induction equation (2.2), ${\mathrm{\partial }}_t\mathbf{B}=\mathcal{I}(\mathbf{u},\mathbf{B})$, but in view of the above discussion, we need to take care, particularly for the boundary terms.

We appeal to a reduced version of the magnetostrophic equations in order to probe what can be said about the behaviour of B(t) on the boundary at t = 0. Assuming that u(t) is given and is independent of B, the induction equation (2.2) is of standard parabolic form (like the heat equation), so its solution is smooth for all t > 0. If the initial condition B(0) is also smooth and satisfies the boundary condition (2.1), then the solution is smooth also at t = 0, except possibly at r = 1. For the solution to be smooth everywhere, including at r = 1, and for Taylor's substitution to be valid, we need the initial condition to satisfy not only the usual boundary condition (also termed the zero-order boundary conditions) but also the first-order boundary conditions: that ∂_tB, given by $\mathcal{I}(\mathbf{B},\mathbf{u})$ of (2.2) satisfies the boundary condition (2.1). Higher-order variants of the boundary conditions pertain to higher-order time derivatives. Assuming that this analysis extends to the full magnetostrophic equations, it provides strong constraints on the form of the initial condition that produces a solution that is smooth for t≥0 and all r≥0.

This issue of lack of smoothness of B occurs only instantaneously at t = 0. One may ask if it is possible to specify an initial field that satisfies Taylor's constraint and higher-order boundary conditions, making it possible to use equation (4.11) directly. Although in principle the answer is yes, it would be practically impossible because an evaluation of the first-order boundary condition requires knowledge of ∂_tB and therefore u_g. The logic is therefore circular: we need to know u_g in order to check the method that enables us to find u_g in the first place. It would seem that some additional insight or good fortune would be required to find a geostrophic flow that is self-consistently satisfies the boundary conditions. The complication compounds the already difficult task of finding an initial condition that satisfies the necessary condition of being a Taylor state.

It is worth noting, however, that once the system has evolved past the initial condition many of these problems vanish. For t > 0, solutions to parabolic systems are smooth and so automatically satisfy all higher-order boundary conditions. It follows that equation (4.11) is valid for t > 0, although this does not help find the geostrophic flow at t = 0.

(c) Schemes in which the boundary information is included

These concerns described above regarding boundary conditions do not carry over to the axisymmetric case, the plane layer situation nor the three-dimensional implicit schemes described. In the axisymmetric and Cartesian cases (e.g. [31]), the boundary conditions are evaluated to zero and the boundary value of the magnetic field or any of its time derivatives never enter any subsequent calculations. In the three-dimensional implicit scheme, because of the representation of all quantities (including B and any of its time derivatives) in terms of a Galerkin basis, boundary conditions to all orders are satisfied.

Thus in the axisymmetric and Cartesian cases, equation (4.8) and equation (4.4) are correct irrespective of the initial choice of Taylor state, as is the fully implicit method of §4b for the three-dimensional case. This is to be contrasted with (4.11) that is valid only for the subset of Taylor states satisfying zero- and first-order boundary conditions.

6. An appraisal of Taylor's method

(a) An illustration of when Taylor's method fails

We are now in a position to provide a first explicit demonstration that Taylor's ODE equation (4.11) fails when using an initial Taylor state that does not satisfy first-order boundary conditions. We show this in two parts. Firstly, within axisymmetry, we demonstrate that Taylor's equation (4.11) is formally inconsistent with the BWR equation (4.8); secondly, we plot an explicit solution of Taylor's equation and show that it does not agree with those derived from other methods known to be correct. In §§8–10, we will show that our generalized version of Taylor's method shows agreement among all methods.

We consider the simple case of the dipolar, single spherical harmonic l = 1 axisymmetric poloidal magnetic field

where $A=\sqrt{231/20584}$ is a scaling constant (see §2a). We note that B satisfies the electrically insulating boundary conditions (2.1), and is an exact Taylor state owing to its simple symmetry.

The ageostrophic flow (determined for example by the method described in appendix Ac) has only an azimuthal component given by

For this choice of B, equation (4.8) then provides an exact expression for the first derivative of u_g. Substituting this into Taylor's second-order equation (4.11) renders it unbalanced. Thus none of the solutions of the first-order ODE (4.8), satisfy the second-order ODE (4.11). Full details of this are given in the electronic supplementary material.

This specific case (which is illustrative of the general case) shows that equations (4.11) and (4.8) are inconsistent: in particular, the first-order equation (4.8) is not simply the first integral of the second-order equation (4.11). The reason why they are not consistent is that although the ODEs are derived from the equivalent forms (3.2) and (1.2), the boundary terms are used to derive (4.8) but not (4.11). Thus, the two equations embody different information. In this example, Taylor's method is equivalent to the erroneous replacement of ∂_tB_ϕ (which is zero) in the boundary term of (5.3), by ${\mathcal{I}}_{\varphi }\ne 0$. This can be seen in the expression for the coefficient β given in equation (6.10), where the boundary term is such that it does not vanish in the axisymmetric case. While the initial magnetic field has been chosen such that it satisfies the boundary condition (2.1), through computing ∂_tB we can show that, based on Taylor's solution, the initial rate of change of the magnetic field violates this boundary condition.

To confirm that Taylor's method is not generally valid, we now directly compare solutions from various methods. Integrating equation (4.8) analytically gives the solution

We note that the solution is a sum of odd polynomials, an s ln(s) term and additional (and non-singular) ln and arctan terms. The constant C is determined through enforcing zero solid body rotation (equation (2.3)). The solution for u_g is everywhere continuous and finite, only at s = 0 is there a weak singularity: ∂_s(u_g/s)∼1/s. We also observe that there is no singularity at s = 1. A comparable analytic solution but for a quadrupolar axisymmetric magnetic field was given in Li et al. [15], which is also regular everywhere except for a weak sln(s) singularity at s = 0. That the analytic expression (6.2) is indeed the true solution is confirmed by figure 1a which compares it to the geostrophic flow given by the independent three-dimensional implicit scheme of §4b; the two solutions over-plot. A contour plot of the total azimuthal flow is shown in §9 (figure 6a).

• Download figure
• Open in new tab
• Download PowerPoint

We now directly compare this solution with that obtained by solving Taylor's equation (4.11), shown as the dotted line of figure 1a. This solution is found by adopting the expansion (4.7) and minimizing the integrated squared residual

with respect to the spectral coefficients, whose truncation is increased until the solution converges.

Although all solutions agree at small s, Taylor's solution shows significant differences from the others for s > 0.8.

It is also of interest to assess numerical convergence to solutions of equations (4.11) and (4.8). Although we have an analytic solution to (4.8), we use the same numerical method as given above but now applied to (4.8) by minimizing

Figure 2 demonstrates that convergence of the solution is faster for the correct, first-order equation (4.8) than for Taylor's equation (4.11). Therefore, aside from Taylor's equation being generally inapplicable, it seems that converged solutions are also relatively more difficult to find.

• Download figure
• Open in new tab
• Download PowerPoint

(b) Specific cases when Taylor's method succeeds

For arbitrary purely toroidal Taylor states bounded by an electrical insulator, B vanishes on r = 1 and in this special case Taylor's methodology is correct. This is because the boundary term involving B_ϕB_r (see equation (3.2)) has a ‘double zero’ and so, when considering its time derivative, erroneous substitution for ∂_tB leaves it invariant as zero.

Taking the time derivative of (3.2), noting that the boundary term is zero, we obtain

Using the three-dimensional extension of (4.3) where ∂₁/∂ϕ is a derivative with respect to ϕ that leaves invariant the unit vectors (e.g. [29]), the term involving u_g in (6.5) becomes Noting that the last integral is zero, we obtain an equation (that holds in three dimensions) that is of the same form as the axisymmetric BWR equation (4.8) As an illustration, we consider the non-axisymmetric l = 1, m = 1 toroidal magnetic field where A is a scaling constant which takes the value $\frac{3}{4}\sqrt{105}$. The ageostrophic flow is and, solving (6.7), the geostrophic flow is where C₁ is determined through considerations of angular momentum. Note the absence of singularities in this solution.

This geostrophic flow is shown in figure 1b, and we note that the three-dimensional implicit method and Taylor's method give the same solution (not shown).

It is in fact simple for us to show analytically that for any purely toroidal field, Taylor's equation (4.11) and equation (6.7) are equivalent, up to the requirement of a further boundary condition for the second-order differential equation (4.11).

We note that via integration by parts, (4.12) can be written in the following way, which allows identification of the boundary term present within Taylor's method, as the second term in the following expression for the coefficient β:

We observe that since B_r = 0 for a purely toroidal field, then the boundary term within equation (6.10) will always vanish in this case, reducing Taylor's equation (4.11) to the BWR equation (6.7).

7. A generalization of Taylor's analysis

To modify the method of Taylor so that it applies to a magnetic field that does not satisfy the first-order boundary conditions, we use (5.3) to impose stationarity of the Taylor constraint. Equally, we could impose stationarity of the equivalent equation (3.2) but it is simpler to avoid the additional integral in s. Bearing in mind our discussion in §5, we take particular care to ensure correct handling of the boundary term.

The magnetic field matches continuously (since η≠0) with an external potential field within the mantle r≥1. Note that our assumption of a globally continuous solution differs from the case when η = 0, for which horizontal components of B may be discontinuous on r = 1 [32]. In our setting where η≠0, the field matches continuously but not necessarily smoothly across r = 1. We note however that owing to ∇ · B = 0, the radial component of B (and all its time derivatives) are always smooth at r = 1 (e.g. [33]): thus only the horizontal components B_θ and B_ϕ are not in general smooth.

Thus, in the first term of equation (5.3), we may substitute at t = 0

For the second (boundary) term, we may substitute ${\mathrm{\partial }}_t{B}_r={\mathcal{I}}_r(\mathbf{u},\mathbf{B})$ but the initial value of ∂_tB_ϕ at r = 1 is not specified by ${\mathcal{I}}_{\varphi }$ alone, as assumed in Taylor's derivation.

The key remaining issue is then to find the initial boundary value of ${\dot{B}}_{\varphi }$, for which we present three methods below. Having done this, all terms are defined and (5.3) provides an implicit determination of u_g up to the usual considerations of solid body rotation and regularity.

We observe that the form of equation (5.3) differs markedly from equations (4.11) and (4.8): in addition to the spatial derivatives of u_g (in the leftmost term), there is an explicit boundary term. For the general case, this boundary term must be retained, although it may be neglected under certain circumstances: e.g. those of §§6b and 9.

We remark that the above instantaneous method can be amended to a first-order implicit scheme (akin to equation (4.4)) by considering

As before, this equation is applicable even when Γ_z≠0, that is, if the solution is close but not exactly on the Taylor manifold.

(a) A potential-based spherical transform method

One way to find ${\dot{B}}_{\varphi }$ on r = 1 is to note that it is the azimuthal component of the potential field in r≥1

The potential $\dot{V}$ is itself determined through continuity of the radial component ${\dot{B}}_r$ at r = 1 and thus depends upon u_g. This method of determining ${\dot{B}}_{\varphi }$ has been introduced in the study of torsional waves by Jault [29], but is implemented here for the evaluation of the geostrophic flow.

The time derivative of the potential $\dot{V}\hspace{0.17em}$ can be written in terms of orthonormal spherical harmonics Y_lm with unknown coefficients a_lm as

where 0≤l≤L_max and −l≤m≤l and where Ω is an element of solid angle. It follows then that on r = 1 Key to the implementation of this method here is a spectral expansion of u_g, for example (4.7), because it allows ${\dot{B}}_r$ (which depends on the I + 2 spectral coefficients of u_g) to be evaluated everywhere on the boundary, as required in the above spherical transform. This is to be contrasted, for example, with a finite difference representation of u_g where no such evaluation is possible.

To find u_g, we note that all time-derivative terms in the left-hand side of (5.3), including those evaluated on the boundary, are linear in the unknown coefficients (A₀, A₁, …, A_I, B), and hence the residual is of the form

for some functions a_i, b and c that depend on B and u_a. We formulate a single equation for the coefficients defining u_g by minimizing the quantity ${\int }_0^1{R}^2\hspace{0.17em}\mathrm{d}s$ (which is quadratic in the coefficients that we seek). Note that the solution is approximate and depends on two parameters I and $L_{max}$, which represent the truncation of the expansion used and care must be taken to ensure we achieve convergence in each.

(b) A potential-based Green's function method

An alternative method for determining the potential $\dot{V}\hspace{0.17em}$ at the core mantle boundary is through the use of a Green's function convolved with ${\dot{B}}_r$ on r = 1. Following [34,35], the relevant Green's function associated with the Laplace equation in the exterior of a sphere with Neumann boundary conditions is

where x = 1/r, f = (1 − 2xμ + x²)^1/2, μ = cosθcosθ′ + sinθsinθ′cos(ϕ − ϕ′). This can be expressed as N(x, μ) = N(1/r, θ, θ′, ϕ − ϕ′), which is the potential at location (r, θ, ϕ) in r≥1 due to a singularity of unit strength in the radial field at (θ′, ϕ′) on the core–mantle boundary. Making use of the periodicity of ϕ, the magnetic potential in the region r≥1 can then be written as and so Like the previous method, this procedure of evaluating ${\dot{B}}_{\varphi }$ on r = 1 requires an integral over all solid angle. Using again our spectral expansion (4.7), this results in ${\dot{B}}_{\varphi }$ being a linear function of the unknown spectral coefficients; thus using equation (5.3) the geostrophic flow can then be determined as in §7a.

(c) A modal projection

A further alternative method to find ${\dot{B}}_{\varphi }$ on r = 1, which does not rely on a magnetic potential, is to employ a modal basis set for the magnetic field that is complete and satisfies the required boundary conditions. Here, we adopt a numerically expedient Galerkin basis set (see appendix Aa for details), whose orthonormal poloidal and toroidal modes are written, respectively, as ${\mathcal{S}}_{(l,n)}^m$ and ${\mathcal{T}}_{(l,n)}^m$.

By using such a representation, boundary conditions to all orders are automatically satisfied and therefore a direct substitution of the projected representation of $\mathcal{I}$,

for ∂_tB in all three components for the whole sphere r≤1 is justified. In the above, l is bounded by $L_{max}$, $0\le n\le N_{max}$ and $\overline{\mathbf{x}}$ indicates the modal projection of x (see appendix Ab).

As before, key to the method here is the spectral representation (4.7) for u_g; the coefficients c_l,m,n and d_l,m,n, found by integration (see appendix Ab), then depend linearly on the unknown coefficients A_i and B.

Equation (5.3) can be then written as the following, in which u_g appears explicitly

and where the modal projection onto the Galerkin basis is required for an accurate representation of ${\dot{B}}_{\varphi }$ within the boundary term. Note that it is not necessary to project the term representing ${\dot{B}}_r$, due to the fact that the radial component of a divergence-free field must be smooth at the boundary.

This approach may be considered as the most direct generalization of the BWR equation (4.8) to three dimensions. We note that under the assumption of axisymmetry, equation (7.4) can be directly integrated to obtain the BWR equation (4.8).

Although on one level, a simpler method than those previously presented because we do not need to calculate $\dot{V}\hspace{0.17em}$, in fact the method is more computationally expensive for two reasons. First, we need to check convergence in three parameters: $I,L_{max},N_{max}$, rather than just the first two; second, because the orthonormality requires an integration over radius, in addition to the integration over solid angle required by both methods.

8. Examples of the geostrophic flow in three dimensions

We now give some examples to illustrate our generalized methodology for computing the instantaneous geostrophic flow associated with three-dimensional Taylor states, using our spherical-transform method. These will be compared with the solution obtained using the fully implicit three-dimensional method with a very small timestep of h = 10⁻⁹; in all cases, the solutions overplot. In none of the cases is an analytic solution available for comparison. For further comparison, we plot also the solution of Taylor's ODE (see equation (6.3)).

We consider firstly an example of a non-axisymmetric l = 2, m = 2 poloidal magnetic field

where $A=1/(6\sqrt{390})$. Figure 3 shows that the implicit and instantaneous solutions agree, whereas similar to the axisymmetric case of figure 1a we can see that Taylor's solution differs significantly particularly near s = 1.

• Download figure
• Open in new tab
• Download PowerPoint

For all our three-dimensional solutions, the expansion for u_g differs from that in axisymmetry given in equation (4.7). We now do not include a logarithmic term. As discussed in §10a, the logarithmic behaviour is not expected outside of axisymmetry and would violate the assumed regularity of the magnetic field.

The approximate polynomial solution, with coefficients rounded to five significant figures, is

where the expansion has been truncated at s¹⁹ and convergence achieved with parameters $I=L_{max}=20$.

We secondly consider a more complex example of a non-axisymmetric magnetic field, which contains both l = 2, m = 1 toroidal and poloidal components

where $A_t=\frac{5}{4}\sqrt{21}$ and $A_p=\sqrt{7/262440}$. Figure 4 shows that again the solution using the instantaneous method is validated by the implicit method, whereas Taylor's solution deviates as s → 1. The figure also shows the geostrophic flow generated separately by either the purely toroidal or purely poloidal magnetic field component, each individually a Taylor state. As anticipated by the structure of the equation for u_g (nonlinear in B), the geostrophic flow driven by the total field does not equal the sum of the individually driven geostrophic flows.

• Download figure
• Open in new tab
• Download PowerPoint

9. Analytic approximation for an Earth-like field

Based on the present structure of the geomagnetic field, various studies show that it is reasonable to neglect the boundary term in equation (3.2) in an Earth-like context [1,36]. This is because not only is the magnetic field likely much stronger inside the core than on r = 1, but also because only the non-axisymmetric field contributes to the boundary term and it is relatively weak. The estimated strength of the magnetic inside the core is 5 mT, and that of the non-axisymmetric field on r = 1 is 0.5 mT; therefore, the relative magnitude of the boundary to the interior terms is about 1/10² or 1%. The negligible effect of the boundary term has been verified in the case of related studies of torsional waves [26,37].

Should we neglect the boundary term entirely, then the geostrophic flow is described by the same equation (6.7) that pertains to a purely toroidal field, whose solution is

If α(s) > 0, then this equation is integrable. A continuous solution for u_g does not exist, however, if B²_s is everywhere zero on a geostrophic cylinder C(s*) (rendering α(s*) = 0). Physically, this would mean that the magnetic field fails to couple cylinders on either side of s = s*, leading to a discontinuity in the geostrophic flow.

In the Taylor states we use, B is of polynomial form and it then follows that S and α are also polynomial (up to a square root factor arising from the geometry) and therefore u_g can (in general) be found in closed form. We note that, in general, S/α is O(1) and so u_g behaves as sln(s) as s → 0.

As an example of this approach, here we construct an Earth-like Taylor state comprising an axisymmetric poloidal mode and a non-axisymmetric toroidal mode, scaled such that the magnitude of the asymmetric part is 20% of the magnitude of the axisymmetric part, but that the total RMS field strength is unity:

where $A_t=\sqrt{28875}/4$ and $A_p=1/\sqrt{966}$. The analytic solution of (9.1) is which is shown in figure 5 and compared to our solution by the method in §7a in which full account is taken of the boundary terms. As anticipated, the two solutions are very similar and diverge only close to s = 1 (where the boundary term has most effect), with an RMS difference of about 1%, all of which occurs very close to the outer boundary. This validates the neglect of the boundary term for this example, and indicates the significance of equation (9.1) which can be used with confidence to analytically approximate the geostrophic flow generated by an Earth-like field. However, we note the presence of a logarithmic singularity that (in view of an earlier comment) that we do not expect in a non-axisymmetric case; this is discussed in the following section.

• Download figure
• Open in new tab
• Download PowerPoint

Finally, figure 6b shows contours of the total azimuthal component of the flow. Of note is the much higher amplitude of flow associated with the increased complexity of the magnetic field compared with the single-mode magnetic field example of figure 6a. The scale of this flow is as would be expected geophysically: maximum dimensionless velocities are of order 100, corresponding to dimensional velocities of order 10⁻⁴ ms⁻¹ consistent with large-scale core flows inferred by secular variation [38].

• Download figure
• Open in new tab
• Download PowerPoint

10. Singularities of u_g

A key benefit of having an instantaneous description of the geostrophic flow is to make explicit its analytic structure, which then motivates spectral expansions such as (4.7) for use with other methods. Assuming α(s) > 0, because the equation describing u_g is smooth and regular, u_g is expected to be an odd [39] finite function on 0 < s < 1. There are three places however where the solution may be singular: (i) s = 0; (ii) s = 1 and (iii) in the complex plane s = x + iy, away from the real axis (y≠0). We discuss each in turn.

(a) Singularities at s = 0

Firstly, we consider the presence of a singularity at s = 0. In axisymmetry, it is well established that u_g∼sln(s) as s → 0, resulting in a s⁻¹ singularity in ∂_s(u_g/s) [13,14,40], reproduced in our example (6.2). However, it has not been quite clear whether the logarithmic singularity pertains to a general asymmetric Taylor state: in particular, in axisymmetry s = 0 is a singular line of the coordinate system, whereas in three-dimensional spherical coordinates the only singular point is the origin r = 0. Roberts & Wu [36] showed that either by neglecting the boundary term (their (25a)) or considering Taylor's ODE directly, which we have shown to be of limited validity, (see their appendix B) leads to a general logarithmic behaviour.

At first inspection, it appears that the boundary term is negligible as s → 0. For a general three-dimensional field, both B and $\dot{\mathbf{B}}$ are $\mathcal{O}(1)$ on s = 0, suggesting that the interior term in equation (5.3) is $\mathcal{O}(1)$, whereas the boundary term is $\mathcal{O}(s)$ as $s\to 0$. Motivated by the example in §9, this suggests that a full treatment (including the boundary term) retains the singularity in three dimensions—however, we do not find this to be the case. Significant cancellation in the interior term occurs and while the integrand is O(1), the integral itself is O(s), as expected since we know that the interior term and boundary term must sum to zero for all s. Therefore, there is no evidence that the three-dimensional case has a logarithmic singularity at s = 0, and indeed all our numerical solutions and analytic solutions are regular there. In the purely toroidal field explored in §6b, the analytic solution given in equation (6.9) is purely polynomial, with no singular behaviour at the origin.

Theorem

This assertion can be strengthened into a theorem. The assumption of a magnetic field that is regular initially and remains so for all time places a restriction on the permitted behaviour of the geostrophic flow. In axisymmetry, the space of solutions allows a weak singularity in the geostrophic flow at s = 0. However, in three dimensions it is required that the geostrophic flow is regular at the origin in order to maintain regularity of the magnetic field.

Proof.

This result directly follows from the form of the geostrophic term in the induction equation. In axisymmetry this is given by equation (4.3), from which it is clear that it is permissible for u_g to contain a weak logarithmic singularity while maintaining a regular B. In three dimensions, the geostrophic term in the induction equation is given by equation (6.6). In the presence of a non-axisymmetric magnetic field, any logarithmic singularity in u_g would render ∂_tB non-regular. Hence the assumption of regular B(t) is incompatible with such a singular solution. ▪

While the analytic approximation in §9 is shown to produce accurate geostrophic flows for Earth-like magnetic fields, it should be used with caution, since the analytic structure of the solution will contain an slns dependence, that does not persist when the full balance including the boundary term is considered. For axisymmetric magnetic fields, this weak logarithmic singularity is not a significant concern since the geostrophic flow only enters the induction equation through ∂_s(u_g/s) and so the magnetic field remains regular everywhere. By contrast, in three dimensions the structure of the geostrophic term in the induction equation (given in equation (6.6)) means that the logarithmic singularity is imparted to the magnetic field itself, causing the magnetic field to diverge at the rotation axis and violating the standard assumption of a regular field. Thus, in a practical implementation, such singular behaviour must be filtered out of u_g.

(b) Singularities at s = 1

We also address the possible existence of a singularity at s = 1. For the specific case of an axisymmetric dipolar magnetic field, Roberts & Wu [36] presented an argument that ∂_su_g∼(1 − s²)^−1/2, although they conceded that this was not supported by their numerical examples. The same form of singular behaviour for u_g has been predicted for torsional waves [41,42], perturbations to Taylor states, whose eventual steady state at t = ∞ would be exactly magnetostrophic (if indeed steady Taylor states exist). However, there is no reason why the analytic structure of the oscillations should mirror that of the underlying background state, particularly as the manner of how the limit t → ∞ is reached at the end points where the wave speed may vanish is unclear [15,43].

Although we are not in a position to prove one way or the other the existence of singular behaviour at s = 1, we demonstrate by example that it is not generally present.

We find no singularity at s = 1 in the non-axisymmetric example of §8. A similar regular behaviour is shown in figure 7 (solid curve) for an axisymmetric example. Interestingly, for this latter case, the application of Taylor's ODE (which is invalid for this example) gives a solution that does show a singularity at s = 1 (dotted curve). In this instance, singular behaviour is simply an artefact of applying Taylor's ODE when it is not valid, and we have found no cases where a solution to our more general analysis behaves singularly at s = 1.

• Download figure
• Open in new tab
• Download PowerPoint

This observation may help explain why the prediction of a singularity at s = 1 [36] is not borne out in any numerical examples. They themselves discussed this discrepancy and hypothesized that a key issue is the lack of boundary information contained within Taylor's equation. We speculate that should their magnetic field satisfy not only Taylor's constraint and the boundary conditions but also crucially the first-order boundary conditions, that this singular behaviour will vanish and the geostrophic flow will remain regular at s = 1. We note however that certain magnetic forcing terms can render the geostrophic flow singular at s = 1: for example, that of a non-polynomial mean-field α-effect described in appendix F of Li et al. [15].

Finally, we remark that for a dipolar axisymmetric Taylor state, both Li et al. [15] and Roberts & Wu [43] showed evidence of non-singular but abrupt boundary-layer like behaviour close to s = 1, possibly because the equation describing the geostrophic flow is null at the equator (i.e. α = S = 0). A similar result was also found by Fearn & Proctor [40] who abandoned constraining their geostrophic flows near s = 1 due to anomalous behaviour. We note, however, in our analytical solutions, we find no evidence of such behaviour: for example, figure 1a shows a smooth solution at s = 1.

(c) Singularities off the s-axis

Finally, inspecting an example solution (6.2) shows that there can be either branch cuts or logarithmic singularities away from the real line. These do not affect the solution itself (defined on the real interval 0≤s≤1) but can influence convergence of the numerical method used to find u_g [44]. The closer the singularities lie to the real interval [0, 1] the slower the convergence. In general, we speculate that such singularities can lie arbitrarily close to the real line, possibly being associated with the breakdown of the magnetostrophic balance, for example, torsional waves.

11. Discussion

In this paper, we have discussed in some detail how the geostrophic flow, a fundamental part of any magnetostrophic dynamo, might be determined. Of particular note is that we have shown why the method introduced by Taylor [9] fails in most cases, because of its intrinsic (and, to date, unrecognized) assumption that the initial magnetic field structure must satisfy a higher-order boundary condition (that is, both the magnetic field and its time derivative must satisfy matching conditions pertaining to an exterior electrical insulator). We presented a generalized version of Taylor's method valid for an arbitrary initial magnetic Taylor state that is not subject to higher-order boundary conditions. In many of our examples, the magnetic fields of dimensional scale 1.7 mT drive flows of magnitude about 10⁻⁴ ms⁻¹, comparable to large-scale flows inferred for the core [38]. Thus, in concert with weakly viscous models, inviscid models also produce Earth-like solutions.

A broader point of note is the extent to which the restriction on the validity of Taylor's approach impacts the related derivation of the equation describing torsional waves [26]. A general treatment of torsional waves includes boundary terms, whose proper evaluation would require a method such as described in Jault [29]. However, the troublesome boundary terms are usually neglected, either because of axisymmetry or because of arguments based on the relative size of the asymmetric magnetic field [1]. Either way, these approaches remain unconstrained by any consideration of higher-order boundary conditions on the magnetic field and the theoretical description remains correct. However, in §10a we describe the danger of neglecting the boundary term, this leading to a logarithmic singularity not present in solutions of the full equation. This has potential implications for analysis of torsional waves, for which the avoidance of a logarithmic singularity may require the full boundary term.

It is worth noting that the weak logarithmic singularity u_g∼sln(s) as s → 0 in axisymmetric magnetostrophic models stands in contrast with weakly viscous models which are anticipated to be regular everywhere. For example, the asymptotic structure is u_g = O(s) in axisymmetry for both no-slip and stress-free boundary conditions (using the formulae summarized in eqns 8 and 9 of Livermore et al. [21] and the fact that B_s, B_ϕ∼s as s → 0) and u_g = O(1) in non-axisymmetry (using the formulae in eqn (33) of Hollerbach [45] and the fact that ([∇ × B] × B)_ϕ∼s through the properties of general vectors described by Lewis & Bellan [39]). The presence of a weak logarithmic singularity is therefore a feature unique to the axisymmetric inviscid case, and serves to distinguish the exact magnetostrophic balance (with zero viscosity) from models with arbitrarily small but non-zero viscosity. However, in three dimensions, there is no such distinction between the structure of u_g between E = 0 and E≪1: in both cases u_g is regular.

Given that the geometry of the outer core of the Earth is a spherical shell rather than a full sphere, a natural question to ask is how we would calculate the flow within this domain. The method for determining the ageostrophic flow would remain comparable although it could be discontinuous or singular across the tangent cylinder $\mathcal{C}$, the geostrophic cylinder tangent to the solid inner core [46]. As for the geostrophic flow, in the absence of viscosity, there is no reason why it must be continuous across $\mathcal{C}$; there are no known matching conditions that it must satisfy and such an analysis lies far beyond the scope of this work.

Although supplying an analytic structure of the evolving magnetostrophic flow, an instantaneous determination of the geostrophic component is not itself of practical use within a numerical method using finite timesteps of size h, as the solution will immediately diverge from the solution manifold [11]. However, as for the axisymmetric-specific method of Wu & Roberts [14], our three-dimensional instantaneous methods generalize simply to schemes that are accurate to first order in h, thus presenting a viable method for numerically evolving a three-dimensional magnetostrophic dynamo. A direct comparison of this method with the fully implicit (three-dimensional) method of Li et al. [15] would be an interesting study. Indeed, our three-dimensional first-order-accurate solutions could be used as a starting guess for their nonlinear iterative scheme, enabling much larger timesteps to be taken for which the geostrophic flow does not need to be close to its structure at the previous step.

Lastly, there is mounting evidence that rapid dynamics within the core is governed by quasi-geostrophic (QG) dynamics, in which the flow is quasi-invariant along the axis of rotation [47,48]. We briefly comment on whether the slowly evolving background magnetostrophic state is also likely to show such a structure. Both Li et al. [15] and Wu & Roberts [14] show axisymmetric magnetostrophic solutions that have largely z-invariant zonal flows. Here, in our three-dimensional cases, we also find that the geostrophic flow is comparable in magnitude to the ageostrophic zonal flow. In our Earth-like example, comparing figures 5 and 6b, the maximum value of the geostrophic flow is about one-quarter of that of the total zonal flow. Furthermore, our (large-scale) magnetostrophic solutions contain a significant z-invariant component, a finding that is consistent with [49] who have recently suggested the existence of a threshold lengthscale, below which the geodynamo is magnetostrophic and above which the dynamics are QG.

Data accessibility

Maple worksheets that reproduce our examples are included as electronic supplementary material.

Author's contributions

P.L. is responsible for the original ideas behind the work with input from K.L. C.H. wrote the Maple worksheets with input from J.N. and calculated the example solutions, which were validated by J.L. using Mathematica. C.H. and P.L. wrote the majority of the paper with input from J.N. and J.L. All authors read and commented on the manuscript.

Competing interests

We declare we have no competing interests.

Funding

This work was supported by the Engineering and Physical Sciences Research Council (EPSRC) Centre for Doctoral Training in Fluid Dynamics at the University of Leeds under grant no. EP/L01615X/1. P.W.L. was partially supported by NERC grant no. NE/G014043/1.

Footnotes

Electronic supplementary material is available online at http://dx.doi.org/10.6084/m9.figshare.c.4238159.