Convection in a Rapidly Rotating Spherical Shell: Newton’s Method Using Implicit Coriolis Integration

1. Introduction

Convection in a rotating spherical annulus is a classic configuration in hydrodynamics [1,2,3,4]. At moderate parameter values, it provides beautiful and deep examples of symmetry-breaking bifurcations leading to complex patterns, e.g. [5,6,7,8,9,10,11,12,13]. At more extreme parameter values, the turbulent flow it engenders is among the fundamental processes in planets and stars, e.g. [14,15]. We do not attempt here to survey the considerable literature devoted to this configuration nor to derive new results in either of these well-studied regimes. Instead we describe the implementation and capabilities of a pair of numerical codes – two drivers relying on a common set of subroutines – that incorporate an implicit treatment of the Coriolis force to carry out both timestepping and also continuation via Newtons’s method.

We will illustrate, test, and compare the codes on the simplest convective states, which are rotating waves (RW). Their appearance at onset is a mathematical consequence of the invariance of the system under azimuthal rotation, but not under azimuthal reflection, which is broken by the rotation of the system [16,17,18]. We have carried out numerical continuation in both Rayleigh and Ekman numbers. We have found interesting examples of branches which contain saddle-node bifurcations separating plateaus in drift frequency.

Our numerical algorithm uses the basic framework described in [11,12,13,19] and the spherical harmonic library of [20]. The research which is closest to ours is that of [21,22,23,24], who, like us, have carried out implicit timestepping and continuation for convection in a rotating spherical shell and explored various interesting parameter regimes. There are a number of differences between our continuation method and theirs, some of which have been analyzed in [25] and which we will discuss after we have described our numerical methods. We mention here also that [14,15,26] have implemented implicit Coriolis timestepping.

2. Physical Description

2.1. Governing Equations

We study classical Rayleigh-Bénard convection in a spherical fluid shell rotating with constant angular velocity $\mathrm{\Omega }$ about the z axis. The shell is heated from within by imposing a temperature difference $\Delta T$ between the inner and outer spherical boundaries. Lengths are non-dimensionalized with the gap size d, so that the dimensionless outer and inner radii are $r_{\mathrm{out}}$ and $r_{\mathrm{in}}=r_{\mathrm{out}}-1$, respectively. The aspect ratio can be specified via $r_{\mathrm{out}}$ or via $\eta \equiv r_{\mathrm{in}}/r_{\mathrm{out}}=(r_{\mathrm{out}}-1)/r_{\mathrm{out}}$. Time is scaled by the viscous diffusion time $d^2/\nu$, where $\nu$ is the kinematic viscosity and the flow velocity $\mathit{u}$ is scaled by the viscous diffusion velocity $\nu /d$. Measuring temperature from a reference temperature $T_0$ at which the mass density is ${\rho }_0$, we nondimensionalize pressure by ${\rho }_0\nu$, and temperature by $\Delta T$.

Figure 1. Geometry. The domain is the spherical annulus between an inner and outer sphere of radius $r_{\mathrm{in}}$ and $r_{\mathrm{out}}$ whose temperatures are fixed at $T=1$ and $T=0$. Both spheres rotate around the z axis with angular velocity $\mathrm{\Omega }$.

Figure 1. Geometry. The domain is the spherical annulus between an inner and outer sphere of radius $r_{\mathrm{in}}$ and $r_{\mathrm{out}}$ whose temperatures are fixed at $T=1$ and $T=0$. Both spheres rotate around the z axis with angular velocity $\mathrm{\Omega }$.

The gravitational acceleration g is assumed to be proportional to the distance r from the center of the sphere (as is the case for the self-gravity in the interior of a spherical body with constant mass density) and is thus of the form $g=-(g_o/r_{\mathrm{out}})r$, where $g_o$ is the acceleration at radius $r_{\mathrm{out}}$. The resulting nondimensional Boussinesq equations read:

$$\begin{array}{cc}\hfill \mathrm{Ek}\left[{\displaystyle \frac{\partial \mathit{U}}{\partial t}}+\left(\mathit{U}·\nabla \right)\mathit{U}-{\nabla }^2\mathit{U}\right]& =-\nabla P+\mathrm{Ra}T{\displaystyle \frac{\mathit{r}}{r_{\mathrm{out}}}}-2{\mathit{e}}_z\times \mathit{U}\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial T}{\partial t}}+\mathit{U}·\nabla T& ={\displaystyle \frac{1}{\mathrm{Pr}}}{\nabla }^{2}T\hfill \end{array}$$

(1b)

$$\begin{array}{cc}\hfill \nabla ·\mathit{U}& =0\hfill \end{array}$$

(1c)

The centrifugal acceleration and a gradient part of the acceleration due to buoyancy have been included in the pressure gradient. The non-dimensional parameters used in (1a) above are the Ekman, Prandtl, and Rayleigh numbers:

$$\begin{array}{c}\hfill \mathrm{Ek}\equiv {\displaystyle \frac{\nu }{\mathrm{\Omega }{d}^2}},\mathrm{Pr}\equiv {\displaystyle \frac{\nu }{\kappa }},\mathrm{Ra}={\mathrm{Ra}}_{\mathrm{rot}}\equiv {\displaystyle \frac{d\alpha g_o\Delta T}{\mathrm{\Omega }\nu }}\end{array}$$

(2)

where $\kappa$ is the thermal diffusivity and $\alpha$ is the thermal expansion coefficient. The Rayleigh number $\mathrm{Ra}$ used here is adapted to the rotating case and is related to the conventional thermal Rayleigh number ${\mathrm{Ra}}_{\mathrm{therm}}$ by

$${\mathrm{Ra}}_{\mathrm{therm}}\equiv {\displaystyle \frac{d^3\alpha g_0\Delta T}{\kappa \nu }}={\displaystyle \frac{d\alpha g_o\Delta T}{\mathrm{\Omega }\nu }}{\displaystyle \frac{\mathrm{\Omega }{d}^2}{\nu }}{\displaystyle \frac{\nu }{\kappa }}={\mathrm{Ra}}_{\mathrm{rot}}{\displaystyle \frac{\mathrm{Pr}}{\mathrm{Ek}}}$$

(3)

In what follows, we will denote ${\mathrm{Ra}}_{\mathrm{rot}}$ merely by $\mathrm{Ra}$. No-slip and fixed temperature boundary conditions are applied at the inner and outer radii:

$$\begin{array}{cc}\hfill \mathit{U}\left(r_{\mathrm{in}}\right)=0,& T\left(r_{\mathrm{in}}\right)=1\hfill \end{array}$$

(4a)

$$\begin{array}{cc}\hfill \mathit{U}\left(r_{\mathrm{out}}\right)=0,& T\left(r_{\mathrm{out}}\right)=0\hfill \end{array}$$

(4b)

and the conductive solution is:

$$\begin{array}{cc}\hfill \mathit{U}=0,& T_{\mathrm{cond}}\left(r\right)={\displaystyle \frac{r_{\mathrm{out}}{r}_{\mathrm{in}}}{r}}-r_{\mathrm{in}}\hfill \end{array}$$

(5)

Following the benchmark study [27], we set the Prandtl number to $\mathrm{Pr}=1$ and the radii to $r_{\mathrm{out}}=20/13$ and $r_{\mathrm{in}}=7/13$ (so that $\eta =7/20=0.35$) throughout the investigation and we will vary $\mathrm{Ek}$ and $\mathrm{Ra}$.

2.2. Overview

The first states to appear at the onset of convection are rotating waves. We denote these by ${\mathrm{RW}}_M$ where M is the azimuthal wavenumber of the rotating wave. Figure 2(a) and Figure 2(b), modeled on those in [12], display properties of some of these rotating wave solutions for $\mathrm{Ek}={10}^{-3}$. Four traveling waves are shown, with azimuthal wavenumber from $M=2$ to $M=5$. At $\mathrm{Ek}={10}^{-3}$, ${\mathrm{RW}}_4$ is the first in Rayleigh number to bifurcate, and therefore it is the only one which is stable at onset. The drift frequency, i.e. the frequency relative to the imposed frequency $\mathrm{\Omega }$, decreases as $\mathrm{Ra}$ increases, from prograde (faster than $\mathrm{\Omega }$) to retrograde (slower than $\mathrm{\Omega }$). These results were obtained via timestepping simulations starting from initial conditions of the form $cos\left(M\phi \right)$.

Figure 3 shows the qualitative evolution of the rotating waves as $\mathrm{Ek}$ is decreased from ${10}^{-2}$ to ${10}^{-5}$: the flow is increasingly stronger near the inner sphere and the azimuthal wavenumber increases. These two properties are related. Each half-wave tends to occupy an equal length in radius ($r_{\max }-r_{\mathrm{in}}$) and in circumference ($\sim \pi r_{\mathrm{in}}/M$), leading to $M\approx \pi r_{\mathrm{in}}/2(r_{max}-r_{\mathrm{in}})$. With decreasing $\mathrm{Ek}$, the radius $r_{\max }$ decreases [28] and so M increases. The ${\mathrm{RW}}_8$ appeared naturally when performing timestepping simulations from an initial condition of the form $cos\left(4\phi \right)$ at $\mathrm{Ek}={10}^{-4}$, as did the ${\mathrm{RW}}_{12}$ at $\mathrm{Ek}=3\times {10}^{-5}$. This result was then used as input for a path-following computation to obtain the ${\mathrm{RW}}_{12}$ data at $\mathrm{Ek}={10}^{-5}$. The same techniques were applied to obtain the ${\mathrm{RW}}_4$ data at $\mathrm{Ek}={10}^{-2}$ starting from the existing results of the same rotating wave at $\mathrm{Ek}={10}^{-3}$.

3. Numerical Methods

3.1. Spatial Representation

The velocity is decomposed into toroidal and poloidal potentials e and f such that

$$\mathit{U}(r,\theta ,\phi ,t)=\nabla \times e{\mathit{e}}_r+\nabla \times \nabla \times f{\mathit{e}}_r$$

(6)

These and the temperature are expanded in spherical harmonics $Y_l^m(\theta ,\phi )$ and Chebyshev polynomials $T_k$. of the shifted radius $x\equiv (2r-(r_{\mathrm{o}}+r_{\mathrm{i}}))/(r_{\mathrm{o}}-r_{\mathrm{i}})$.

$$f(r,\theta ,\phi ,t)=\sum _{k,l,m}{f}_{k,l,m}\left(t\right)T_k\left(x\right)Y_l^m(\theta ,\phi )$$

(7)

To carry out transformations to and from the spherical harmonic space, we use the SHTns library [20]. This is a highly efficient code that uses the FFTW library to perform fast transformations to and from Fourier space, and the recursion relations presented in [29] for optimized computation of the associated Legendre polynomials. Although the library allows for multithreaded computations as well as GPU offloading, these capabilities were not used in the present work.

Dealiasing in the r direction is carried out according to the $3/2$ rule, whereas in the $\theta$ and $\phi$ directions, it is performed by SHTns. To achieve this, the library determines a suitable number of physical points based on the non-linear order of the equations to be solved (two in this case) with the aim of maximizing the efficiency of the FFT algorithm. Moreover, for solutions which have azimuthal periodicity $2\pi /M$, we restrict our calculations to a segment of the sphere $0\le \phi <2\pi /M$, or, equivalently, to those which have Fourier expansions containing only multiples of the wavenumber M; this is easily implemented using options in SHTns. Lastly, we will denote the physical points in $r,\theta$ and $\phi$ as $N_r,N_{\theta }$ and $N_{\phi }\times M$, respectively. The resolutions used will be given in the captions of each of the figures.

The radial component of the curl and of the double curl of equation (1a) are taken, leading to:

$$\begin{array}{cc}\hfill & \sum _{l,m}{\displaystyle \frac{l(l+1)}{r^2}}\mathrm{Ek}\left({\displaystyle \frac{\partial }{\partial t}}-L_l\right)e_{l,m}(r,t)Y_l^m(\theta ,\phi )={\mathit{e}}_{\mathit{r}}·\nabla \times \mathit{F}\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill -& \sum _{l,m}{\displaystyle \frac{l(l+1)}{r^2}}\mathrm{Ek}\left({\displaystyle \frac{\partial }{\partial t}}-L_l\right)L_l{f}_{l,m}(r,t)Y_l^m(\theta ,\phi )={\mathit{e}}_r·\nabla \times \nabla \times \mathit{F}\hfill \end{array}$$

(9)

where

$$\begin{array}{c}\hfill L_l\equiv {\displaystyle \frac{{\partial }^2}{\partial r^2}}-{\displaystyle \frac{l(l+1)}{r^2}}and\mathit{F}\equiv \mathrm{Ek}(\mathit{U}·\nabla )\mathit{U}-\nabla P+{\displaystyle \frac{Ra}{r_{\mathrm{out}}}}T\mathit{r}-2{\mathit{e}}_z\times \mathit{U}\end{array}$$

(10)

Equations (8) and (9) are decoupled in ℓ, m and in potential (e vs f), and so its discrete version involves a block-diagonal matrix with $3$ blocks that are each $N_r\times N_r$, as in Figure 5(a).

The temperature field is also expanded in Chebyshev polynomials and spherical harmonics like (7). Its evolution is governed by the discretized form of (1b).

3.2. Implicit Coriolis Integration

To carry out implicit integration of the Coriolis term, we include the Coriolis force in the left-hand-side and remove it from $\mathbf{F}$, leading to:

$$\begin{array}{cc}\hfill & \sum _{m,l}\left[{\displaystyle \frac{l(l+1)}{r^2}}\left({\displaystyle \frac{\partial }{\partial t}}-\mathrm{Ek}{L}_l\right)-{\displaystyle \frac{2im}{r^2}}\right]e_{lm}(r,t)Y_l^m(\theta ,\phi )\hfill \\ \hfill & +{\displaystyle \frac{2}{r^2}}\left({\displaystyle \frac{l(l+1)}{r}}-{\displaystyle \frac{\partial }{\partial r}}\right)f_{lm}(r,t)sin\theta {\displaystyle \frac{d}{d\theta }}{Y}_l^m(\theta ,\phi )\hfill \\ \hfill & +2{\displaystyle \frac{l(l+1)}{r^2}}\left({\displaystyle \frac{2}{r}}-{\displaystyle \frac{\partial }{\partial r}}\right)f_{lm}(r,t)cos\theta Y_l^m(\theta ,\phi )={\mathit{e}}_{\mathit{r}}·\nabla \times {\mathit{F}}^{\mathrm{implicit}}\hfill \end{array}$$

(11a)

$$\begin{array}{cc}\hfill & \sum _{m,l}-\left[{\displaystyle \frac{l(l+1)}{r^2}}\left({\displaystyle \frac{\partial }{\partial t}}-\mathrm{Ek}{L}_l\right)-{\displaystyle \frac{2im}{r^2}}\right]L_l{f}_{\mathrm{lm}}(r,t)Y_l^m(\theta ,\phi )\hfill \\ \hfill & +{\displaystyle \frac{2}{r^2}}\left({\displaystyle \frac{l(l+1)}{r}}-{\displaystyle \frac{\partial }{\partial r}}\right)e_{lm}(r,t)sin\theta {\displaystyle \frac{d}{d\theta }}{Y}_l^{\left|m\right|}(\theta ,\phi )\hfill \\ \hfill & +2{\displaystyle \frac{l(l+1)}{r^2}}\left({\displaystyle \frac{2}{r}}-{\displaystyle \frac{\partial }{\partial r}}\right)e_{lm}(r,t)cos\theta Y_l^m(\theta ,\phi )={\mathit{e}}_{\mathit{r}}·\nabla \times \nabla \times {\mathit{F}}^{\mathrm{implicit}}\hfill \end{array}$$

(11b)

where

$$\begin{array}{c}\hfill {\mathit{F}}^{\mathrm{implicit}}\equiv \mathrm{Ek}(\mathit{U}·\nabla )\mathit{U}-\nabla P+{\displaystyle \frac{\mathrm{Ra}}{r_{\mathrm{out}}}}T\mathit{r}\end{array}$$

(12)

The recursion relations of the Legendre polynomials

$$\begin{array}{cc}\hfill sin\theta {\displaystyle \frac{d}{d\theta }}{P}_l^{\left|m\right|}& ={\displaystyle \frac{l(l-|m|+1)}{2l+1}}{P}_{l+1}^{\left|m\right|}-{\displaystyle \frac{(l+|m\left|\right)(l+1)}{2l+1}}{P}_{l-1}^{\left|m\right|}\hfill \end{array}$$

(13a)

$$\begin{array}{cc}\hfill cos\theta P_l^{\left|m\right|}& ={\displaystyle \frac{l-\left|m\right|+1}{2l+1}}{P}_{l+1}^{\left|m\right|}+{\displaystyle \frac{(l+|m\left|\right)}{2l+1}}{P}_{l-1}^{\left|m\right|}\hfill \end{array}$$

(13b)

together with

$$\begin{array}{c}\hfill Y_l^m(\theta ,\phi )=N_l^m{P}_l^m(cos\theta )e^{im\phi }{N}_l^m=\sqrt{(2-{\delta }_{m0})(2l+1){\displaystyle \frac{(l-m)!}{l+m)!}}}\end{array}$$

(14)

can be used to transform (11a) and (11b) into:

$$\begin{array}{cc}\hfill & \sum _{m,l}{Y}_l^m(\theta ,\phi )\left\{\left[{\displaystyle \frac{l(l+1)}{r^2}}\left({\displaystyle \frac{\partial }{\partial t}}-\mathrm{Ek}{L}_l\right)-{\displaystyle \frac{2im}{r^2}}\right]e_{lm}(r,t)\right.\hfill \\ \hfill & +{\displaystyle \frac{2}{r^2}}{\displaystyle \frac{(l-1)(l-|m\left|\right)}{2l-1}}\left[\left({\displaystyle \frac{(l-1)l}{r}}-{\displaystyle \frac{\partial }{\partial r}}\right)+l\left({\displaystyle \frac{2}{r}}-{\displaystyle \frac{\partial }{\partial r}}\right)\right]{\displaystyle \frac{N_{l-1}^m}{N_l^m}}{f}_{l-1,m}(r,t)\hfill \\ \hfill & \left.+{\displaystyle \frac{2}{r^2}}{\displaystyle \frac{(l+1+|m\left|\right)(l+2)}{2l+3}}\left[-\left({\displaystyle \frac{(l+1)(l+2)}{r}}-{\displaystyle \frac{\partial }{\partial r}}\right)+(l+1)\left({\displaystyle \frac{2}{r}}-{\displaystyle \frac{\partial }{\partial r}}\right)\right]{\displaystyle \frac{N_{l+1}^m}{N_l^m}}{f}_{l+1,m}(r,t)\right\}\hfill \\ \hfill & ={\mathit{e}}_{\mathit{r}}·\nabla \times {\mathit{F}}^{\mathrm{implicit}}\hfill \end{array}$$

(15a)

$$\begin{array}{cc}\hfill & \sum _{m,l}{Y}_l^m(\theta ,\phi )\left\{-\left[{\displaystyle \frac{l(l+1)}{r^2}}\left({\displaystyle \frac{\partial }{\partial t}}-\mathrm{Ek}{L}_l\right)-{\displaystyle \frac{2im}{r^2}}\right]L_l{f}_{lm}(r,t)\right.\hfill \\ \hfill & +{\displaystyle \frac{2}{r^2}}{\displaystyle \frac{(l-1)(l-|m\left|\right)}{2l-1}}\left[\left({\displaystyle \frac{(l-1)l}{r}}-{\displaystyle \frac{\partial }{\partial r}}\right)+l\left({\displaystyle \frac{2}{r}}-{\displaystyle \frac{\partial }{\partial r}}\right)\right]{\displaystyle \frac{N_{l-1}^m}{N_l^m}}{e}_{l-1,m}(r,t)\hfill \\ \hfill & \left.+{\displaystyle \frac{2}{r^2}}{\displaystyle \frac{(l+1+|m\left|\right)(l+2)}{2l+3}}\left[-\left({\displaystyle \frac{(l+1)(l+2)}{r}}-{\displaystyle \frac{\partial }{\partial r}}\right)+(l+1)\left({\displaystyle \frac{2}{r}}-{\displaystyle \frac{\partial }{\partial r}}\right)\right]{\displaystyle \frac{N_{l+1}^m}{N_l^m}}{e}_{l+1,m}(r,t)\right\}\hfill \\ \hfill & ={\mathit{e}}_{\mathit{r}}·\nabla \times \nabla \times {\mathit{F}}^{\mathrm{implicit}}\hfill \end{array}$$

(15b)

In contrast with (8)-(9), equations (15a)-(15b) are coupled in several ways:

• The spectral coefficients $e_{\mathrm{lm}}$ and $f_{\mathrm{lm}}$ both appear in both equations.
• While each m can be treated independently, components ℓ, $\ell +1$, and $\ell -1$ are coupled.
• The real and imaginary parts of $e_{\mathrm{lm}}$ and $f_{\mathrm{lm}}$ are coupled via the imaginary coefficient $2im/r^2$.

However, two decoupled classes of coefficients appear, one class containing coefficients with odd l for $e_l$ and even l for $f_l$ and a second class with the opposite property. For example the real component of $e_{lm}$ is coupled to its imaginary component and to the real parts of $f_{l\pm 1,m}$. The sums in (15a) and (15b) are represented by block pentadiagonal matrices, as illustrated in Figure 5(b).

3.3. Newton Method

We represent our system of governing equations schematically by

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \mathit{U}}{\partial t}}=\mathcal{N}\left(\mathit{U}\right)+\mathcal{L}\mathit{U}\end{array}$$

(16)

For our problem, $\mathit{U}$ corresponds to the poloidal and toroidal fields e and f and the temperature T. For explicit or implicit Coriolis integration, $\mathcal{N}\left(\mathit{U}\right)$ corresponds to the right-hand-sides of (8)-(9) or of (15a)-(15b), respectively. $\mathcal{L}\mathit{U}$ corresponds to minus the left-hand-sides of the same equations with the time derivative omitted. (The correspondence is imperfect, since $\partial /\partial t$ does not act directly on $\mathit{U}$.)

We represent our timestepping code schematically by the following implicit-explicit Euler method:

$$\begin{array}{cc}\hfill \mathit{U}(t+\Delta t)& =\mathit{U}\left(t\right)+\Delta t\left(\mathcal{N}\right(\mathit{U}\left(t\right))+\mathcal{L}\mathit{U}(t+\Delta t\left)\right)\hfill \\ \hfill ⟹\mathit{U}(t+\Delta t)& =\mathcal{B}\left(\mathit{U}\left(t\right)\right)\equiv {(I-\Delta t\mathcal{L})}^{-1}\left(I+\Delta t\mathcal{N}\right)\mathit{U}\left(t\right)\hfill \end{array}$$

(17)

where $\mathcal{L}$ is linear (though $\mathcal{N}$ need not be nonlinear). Steady states are solutions to

$$\begin{array}{c}\hfill 0=\mathcal{A}\left(\mathit{U}\right)\equiv \mathcal{N}\left(\mathit{U}\right)+\mathcal{L}\mathit{U}\end{array}$$

(18)

or, alternatively

$$\begin{array}{c}\hfill 0=\mathit{U}(t+\Delta t)-\mathit{U}\left(t\right)\end{array}$$

(19)

where $\mathit{U}$ is a solution to the continuous-time differential equation (16). Surprisingly, (19) is also a criterion for stationarity when $\mathit{U}(t+\Delta t)$ is calculated via the time-stepping scheme (17) for any value of $\Delta t$, as can be seen from

$$\begin{array}{cc}\hfill (\mathcal{B}-I)\left(\mathit{U}\right(t\left)\right)=\mathit{U}(t+\Delta t)-\mathit{U}\left(t\right)& ={(I-\Delta t\mathcal{L})}^{-1}\left(I+\Delta t\mathcal{N}\right)\mathit{U}\left(t\right)-\mathit{U}\left(t\right)\hfill \\ \hfill & ={(I-\Delta t\mathcal{L})}^{-1}\left(I+\Delta t\mathcal{N}-(I-\Delta t\mathcal{L})\left)\mathit{U}\right(t)\right)\hfill \\ \hfill & =\underset{\mathcal{P}\left(\Delta t\right)}{\underbrace{{(I-\Delta t\mathcal{L})}^{-1}\Delta t}}\underset{\mathcal{A}}{\underbrace{(\mathcal{N}+\mathcal{L})}}\mathit{U}\left(t\right)\hfill \end{array}$$

(20)

Thus $\mathcal{B}-I=\mathcal{P}\left(\Delta t\right)\mathcal{A}$ and so roots of $\mathcal{A}$ of (18) can be found by computing roots of $\mathcal{B}-I$, which is simply the difference between two large consecutive timesteps. Note that (20) does not use Taylor series; $\Delta t$ can be of any size. Indeed, for $\Delta t$ large, $\mathcal{P}\left(\Delta t\right)\to -{\mathcal{L}}^{-1}$ acts as a preconditioner (approximate inverse) for $\mathcal{A}$, whose poor conditioning is due to that of $\mathcal{L}$. Typically, for Newton solving we use $\Delta t=200$.

Finding the roots of (20) via Newton’s method requires the solution of the linear system

$$\begin{array}{cc}\hfill ({\mathcal{B}}_{\mathit{U}}-I)\mathit{u}& =(\mathcal{B}-I)\mathit{U}\hfill \end{array}$$

(21)

where $\mathit{U}$ is the current estimate and ${\mathcal{B}}_{\mathit{U}}$ is the Jacobian of $\mathcal{B}$ evaluated at $\mathit{U}$, with ${\mathcal{N}}_{\mathit{U}}$ substituted for $\mathcal{N}$ in (17). For example, if $\mathcal{N}\left(\mathit{U}\right))$ contains $(\mathit{U}·\nabla )\mathit{U}$, this term becomes $(\mathit{U}·\nabla )\mathit{u}+(\mathit{u}·\nabla )\mathit{U}$ in ${\mathcal{N}}_{\mathit{U}}\mathit{u}$, as well as in the operators $\mathcal{A}$ and $\mathcal{B}$.

We solve (21) using a Krylov method such as GMRES [30]. In such matrix-free methods, one need only provide the right-hand-side $(\mathcal{B}-I)\mathit{U}$ for the current estimate, and a routine which carries out the action of $({\mathcal{B}}_{\mathit{U}}-I)$ on any vector $\mathit{u}$. The computational cost of a Newton step can be measured by the number of such actions (i.e. GMRES iterations) taken by the Krylov method, since the right-hand-side remains constant throughout the step. The number of GMRES iterations required is low if $({\mathcal{B}}_{\mathit{U}}-I)=\mathcal{P}\left(\Delta t\right)\mathcal{A}$ is well-conditioned; it is for this reason that we take $\Delta t$ large. Once the decrement $\mathit{u}$ is determined by solving (21), it is subtracted from $\mathit{U}$ to form an improved estimate. We accept $\mathit{u}$ as the decrement if the linear system (21) is solved by GMRES to relative accuracy ${10}^{-10}$. For Newton’s method, we accept $\mathit{U}$ as a steady state, i.e. a solution to (19), if $\left|\right|(\mathcal{B}-I)\mathit{U}\left|\right|<{10}^{-7}$.

This method [11,12,13,25,31,32] has been called Stokes preconditioning because the linear operator $\mathcal{L}$ which is integrated implicitly in (17) is usually the viscous term that occurs in Stokes equation. However, it is beneficial to include as many possible other terms in $\mathcal{L}$, as long as they are linear and can be inverted directly. It is for the reason that we have included the Coriolis force in $\mathcal{L}$, as described in Section 3.2. In Section 5, we will compare the cost of computations using Newton’s method with implicit vs. explicit Coriolis integration for various values of the Ekman number by comparing the number of GMRES iterations necessary for each.

3.4. Traveling Waves

Newton’s method can also be used to compute traveling waves in the same way. Azimuthal rotating waves satisfy $\mathit{U}(\phi ,t)=\tilde{\mathit{U}}(\phi -Ct)=\tilde{\mathit{U}}\left(\tilde{\phi }\right)$, where C is an unknown wavespeed. Thus, ${\partial }_t\mathit{U}=-C{\partial }_{\tilde{\phi }}\tilde{\mathit{U}}$. Substituting into (16) and dropping the tildes leads to

$$\begin{array}{cc}\hfill 0& =\mathcal{N}\left(\mathit{U}\right)+C{\partial }_{\phi }\mathit{U}+\mathcal{L}\mathit{U}\hfill \end{array}$$

(22)

The explicit portion of the timestep is augmented to include $C{\partial }_{\phi }$ as well as $\mathcal{N}$:

$$\begin{array}{cc}\hfill \mathit{U}(t+\Delta t)& ={(I-\Delta t\mathcal{L})}^{-1}\left(I+\Delta t(\mathcal{N}+C{\partial }_{\phi })\right)\mathit{U}\left(t\right)\hfill \end{array}$$

(23)

If $\mathit{U}$ is expressed in terms of spherical harmonics in which the $\phi$ dependence is trigonometric, the action of $C{\partial }_{\phi }$ on a Fourier component ${\mathit{U}}_m$ is merely $C{\partial }_{\phi }{\mathit{U}}_m=Cim{\mathit{U}}_m$.

To use Newton’s method to determine the unknown field $\mathit{U}$ and wavespeed C, we substitute $\mathit{U}\to \mathit{U}-\mathit{u}$ and $C\to C-c$ into (22):

$$\begin{array}{cc}\hfill 0& =\mathcal{N}(\mathit{U}-\mathit{u})(C-c)+{\partial }_{\phi }(\mathit{U}-\mathit{u})+\mathcal{L}(\mathit{U}-\mathit{u})\hfill \end{array}$$

(24)

Expanding around $(\mathit{U},C)$ and truncating at first order leads to the linear system that must be solved for the decrements $(\mathit{u},c)$:

$$\begin{array}{c}\hfill ({\mathcal{N}}_{\mathit{U}}+C{\partial }_{\phi }+\mathcal{L})\mathit{u}+c{\partial }_{\phi }\mathit{U}=(\mathcal{N}+C{\partial }_{\phi }+\mathcal{L})\left(\mathit{U}\right)\end{array}$$

(25)

The preconditioner remains ${\mathcal{P}}^{(}{\Delta t)=(I-\Delta t\mathcal{L})}^{-1}\Delta t$ with large $\Delta t$, since $\mathcal{L}$ continues to be responsible for the large condition number of (25).

$$\begin{array}{c}\hfill {\mathcal{P}}^{(}\Delta t)\left[\left({\mathcal{N}}_{\mathit{U}}+C{\partial }_{\phi }+\mathcal{L}\right)\mathit{u}+c{\partial }_{\phi }\mathit{U}\right]={\mathcal{P}}^{(}\Delta t)(\mathcal{N}+C{\partial }_{\phi }+\mathcal{L})\left(\mathit{U}\right).\end{array}$$

(26)

An additional equation must be added to the system to compensate for the additional unknown. We choose a phase condition, more specifically we require the imaginary part of a single component (temperature, toroidal or poloidal field; radial and angular and azimuthal mode) of $\mathit{u}$ be zero, i.e. that the corresponding phase component (whose index we shall call J) of $\mathit{U}$ remains unchanged.

This simple choice suggests a trick for retaining the size of the unknown $\mathit{u}$ rather than using augmented fields, since ${\mathit{u}}_J$ is no longer an unknown. The routine which acts on $(\mathit{u},c)$ is defined such that c is stored in ${\mathit{u}}_J$. At the beginning of an action, c is extracted, ${\mathit{u}}_J$ is set to zero, and the explicit part of the action on $\mathit{u}$ is carried out. Then we multiply the stored value of ${\partial }_{\phi }\mathit{U}$ by c and add the result to the explicit part of the action. When the Krylov method converges and returns $\mathit{u}$, we must again extract c from ${\mathit{u}}_J$, after which ${\mathit{u}}_J$ is set to zero. Effectively, although C and c have been added to the unknowns for the Newton method and for the linear equation, ${\mathit{U}}_J$ and ${\mathit{u}}_J$ are no longer unknowns.

3.5. Continuation

We compute solutions along a branch via Newton’s method as described above. We do not impose additional equations, such as requiring a new solution along the branch to be a certain distance from a prior solution or that the increment along the branch must be perpendicular to some direction. Therefore the only additional ingredients that must be introduced to discuss continuation are the choice of initial estimate for each solution along the branch and the parameters that are prescribed for the Newton iteration.

If the continuation is in Rayleigh number, we increment or decrement this parameter according to the number of Newton iterations required for convergence in the previous step as follows [33]:

$$R{a}^{(i+1)}=R{a}^{\left(i\right)}+\Delta Ra=R{a}^{\left(i\right)}+\alpha (R{a}^{\left(i\right)}-R{a}^{(i-1)})$$

(27)

where ${\mathrm{Ra}}^{\left(i\right)}$, ${\mathrm{Ra}}^{(i-1)}$ denote the Rayleigh numbers used in the two previous continuation steps, $\Delta Ra$ is the increment or decrement, and $\alpha$ is defined by:

$$\alpha ={\displaystyle \frac{N^{\mathrm{opt}}+1}{N^{\left(i\right)}+1}}$$

(28)

$N^{\mathrm{opt}}$ is the target number of Newton iterations. If $N^{\left(i\right)}=N^{\mathrm{opt}}$, then $\Delta Ra$ remains unchanged, whereas it $\Delta Ra$ is reduced (increased) if $N^{\left(i\right)}>(<)N^{\mathrm{opt}}$. The choice of $N^{\mathrm{opt}}$ is guided by two considerations. The first is the level of sampling desired along the branch. The second is economy: a smaller value of $\alpha$ will engender more values of $\mathrm{Ra}$, but each calculation will be faster. For our computations we fixed $N^{\mathrm{opt}}$ to be between 4 and 6.

To choose an initial estimate $\mathit{U}$ for the next solution along a branch at ${\mathrm{Ra}}^{(i+1)}$, first-order and even zero-th order extrapolation (i.e. using just the previously computed field) can be used to follow a smooth and monotonically varying branch. However, quadratic extrapolation is essential for going around turning points. Lagrange interpolation uses the three previous Rayleigh numbers ${\mathrm{Ra}}^{\left(i\right)}$, ${\mathrm{Ra}}^{(i-1)}$, ${\mathrm{Ra}}^{(i-2)}$ and the new Rayleigh number ${\mathrm{Ra}}^{(i+1)}$ from (27) to determine coefficients $a,b,c$ such that

$$\begin{array}{cc}\hfill R{a}^{(i+1)}& =aR{a}^{\left(i\right)}+bR{a}^{(i-1)}+cR{a}^{(i-2)}\hfill \end{array}$$

(29)

after which we set

$$\begin{array}{cc}\hfill 1-1{\mathit{U}}^{(i+1)}& =a{\mathit{U}}^{\left(i\right)}+b{\mathit{U}}^{(i-1)}+c{\mathit{U}}^{(i-2)}\hfill \end{array}$$

(30)

to form a new estimate of the solution. Newton iterations are then used to refine $\mathit{U}$ until the norm of (22) is ${10}^{-7}$. This procedure requires saving the three previous solution vectors and Rayleigh numbers.

3.6. Turning Points

We now turn to the more complicated matter of extrapolating near saddle-node bifurcations, at which $\mathit{U}$ ceases to be a single-valued function of $\mathrm{Ra}$. The code detects that the current step is in the vicinity of a turning point by comparing relative changes between the solution vector $\mathit{U}$, and the Rayleigh number, i.e. we compare

$$\begin{array}{c}\hfill \left|{\displaystyle \frac{\Delta U_J}{U_J}}\right|\equiv \left|{\displaystyle \frac{U_J^{\left(i\right)}-U_J^{(i-1)}}{U_J^{\left(i\right)}}}\right|\mathrm{with}\gamma \left|{\displaystyle \frac{\Delta Ra}{Ra}}\right|\equiv \gamma \left|{\displaystyle \frac{R{a}^{\left(i\right)}-R{a}^{(i-1)}}{R{a}_{\left(i\right)}}}\right|\end{array}$$

(31)

where J is the index denoting the element of the solution vector $U^{\left(i\right)}$ of highest absolute value and $\gamma$ is a constant factor that makes the two magnitudes comparable. For our computations we fixed this constant to be between 10 and 100.

Since the normal form of a saddle-node bifurcation at $(\mu ,x)=(0,0)$ is $\dot{x}=\mu -x^2$, dependent variables (x, $\mathit{U}$) vary like the square root of the control parameter ($\mu$, $\mathrm{Ra}-R{a}_c$) near a turning point. Consequently, as the current solution approaches a turning point, the relative difference $|\Delta U_J/U_J|$ increases until it is greater than that of $\gamma |\Delta Ra/Ra|$. It is at this stage that the code switches to extrapolating quadratically, using the component $U_J$ as the independent variable instead of the Rayleigh number. This allows $\Delta Ra$ to change sign and the continuation to go around the turning point. We replace (27) and (29) by analogous equations in ${\mathit{U}}_J$ to determine ${\mathit{U}}_J^{(i+1)}$ and $a,b,c$, and then replace (29) and (30) by

$$\begin{array}{cc}\hfill U_j^{(i+1)}& =a{U}_j^{\left(i\right)}+b{U}_j^{(i-1)}+c{U}_j^{(i-2)}forj\ne J,and\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill R{a}^{(i+1)}& =aR{a}^{\left(i\right)}+bR{a}^{(i-1)}+cR{a}^{(i-2)}\hfill \end{array}$$

(33)

to set the estimate of the new solution. Note that we have only changed the method of determining the initial guess. $\mathrm{Ra}$ remains fixed to ${\mathrm{Ra}}^{(i+1)}$ during the Newton step, while all the elements of $\mathit{U}$ are allowed to vary.

Continuation follows in this manner until eventually $|\Delta {\mathit{U}}_J/{\mathit{U}}_J|$ exceeds $\gamma |\Delta Ra/Ra|$, at which point the code switches back to using the Rayleigh number as the independent variable. The parameter $\gamma$ can be obtained by analyzing the behaviour of the continuation near a turning point. Indeed, if $\gamma$ is too large, then the code will continue to use $\mathrm{Ra}$ as the independent variable past the point at which there is no solution for the next continuation step, and Newton’s method will not converge. Should this situation arise, the continuation can be restarted at the last converged solution using a suitably reduced value of $\gamma$.

4. Branch Following

4.1. Continuation in Rayleigh Number

We used the continuation code to compute branches of rotating waves. As previously explained, to search for solutions having an M azimuthal periodicity we perform our computations to a fraction of the domain given by $2\pi /M$. The solutions are computed independently of their stability; most of the solutions we compute are unstable, either within the restricted domain $0\le \phi <2\pi /M$, or to instabilities which break this restriction such as rotating waves with a different M.

Fixing the Ekman number, we carried out continuation in Rayleigh number $\mathrm{Ra}$. Most of the branches that we computed showed no unusual features, varying smoothly and monotonically down to their threshold at a supercritical Hopf bifurcation. However, a few branches presented some interesting non-monotonic behavior, which we show in order to display the capacities of our code. Figure 6 tracks the Rayleigh-number dependence of the ${\mathrm{RW}}_4$ branch for $\mathrm{Ek}=3.53\times {10}^{-5}$. This branch consists of three smooth regions that are separated from one another by short intervals which each contain two saddle-node bifurcations and rapid changes in drift frequency. In order to focus in on these intervals, we have computed branches at slightly higher and lower Ekman numbers in Figure 7 and Figure 8. Each of these branches contains a single interval of rapid change. We also computed an ${\mathrm{RW}}_8$ branch for $\mathrm{Ek}=1.26\times {10}^{-5}$, shown in Figure 9. Indeed, for this value of $\mathrm{Ek}$, rotating waves with azimuthal wavenumber 8 are more appropriate, i.e. more likely to be stable. This branch also contains a plateau in drift frequency adjacent to a short interval at $\mathrm{Ra}\approx 159$ of rapid change delimited by two saddle-node bifurcations. These rapid rises and plateaus in drift frequency should have a physical or phenomenological explanation, but we have not been able to find one.

4.2. Continuation in Ekman Number and in Resolution

We adapted our method to follow branches in Ekman number, varying $\mathrm{Ek}$ on a logarithmic scale. We replace $\mathrm{Ra}$ with ${log}_{10}\left(\mathrm{Ek}\right)$ in (27) as follows:

$$\begin{array}{cc}\hfill {log}_{10}\left({\mathrm{Ek}}^{(i+1)}\right)& ={log}_{10}\left({\mathrm{Ek}}^{\left(i\right)}\right)+\Delta \left[{log}_{10}\left(\mathrm{Ek}\right)\right]\hfill \\ & ={log}_{10}\left({\mathrm{Ek}}^{\left(i\right)}\right)+\alpha ({log}_{10}\left({\mathrm{Ek}}^{\left(i\right)}\right)-{log}_{10}\left({\mathrm{Ek}}^{(i-1)}\right))\hfill \end{array}$$

(34)

where the computation of $\alpha$ remains as in equation (28). In this case, we used $N^{\mathrm{opt}}=3$ for all of the simulations. Furthermore, in order to go around turning points, we now compare $|\Delta {\mathit{U}}_j/{\mathit{U}}_j|$ and $\gamma |\Delta \left[{log}_{10}\left(\mathrm{Ek}\right)\right]/{log}_{10}\left(\mathrm{Ek}\right)|$, where we chose $\gamma$ between 400 and 500.

Apart from these minor changes, there are several major differences between performing continuation in Ekman and Rayleigh numbers. The first is that the matrices that represent the left-hand-sides of (8)-(9) or (11a)-(11b) must be recomputed every time a new value of $\mathrm{Ek}$ is chosen. Indeed, while the Rayleigh number only appears in the right hand side of these equations, $\mathrm{Ek}$ is also present in the diffusion terms, which are treated implicitly regardless of how the Coriolis term is treated. This means that extra work must be done at the beginning of every continuation step.

Secondly, as the Ekman number is decreased, the critical Rayleigh number for onset of convection increases. That is, rotation stabilizes the configuration against convection. A classic result is that the critical value of the usual thermal Rayleigh number ${\mathrm{Ra}}_{\mathrm{therm}}$ varies like ${\mathrm{Ek}}^{-4/3}$, so that the critical value of our Rayleigh number ${\mathrm{Ra}}_{\mathrm{rot}}={\mathrm{Ra}}_{\mathrm{therm}}\mathrm{Ek}/\mathrm{Pr}$ varies like ${\mathrm{Ek}}^{-1/3}/\mathrm{Pr}$. Therefore, instead of keeping $\mathrm{Ra}$ fixed as we vary $\mathrm{Ek}$, we have kept $\mathrm{Ra}{\mathrm{Ek}}^{1/3}$ fixed by increasing $\mathrm{Ra}$ every time we decrease $\mathrm{Ek}$.

Third, as the simulation ventures towards lower Ekman numbers, the fields require a higher resolution. This is clearly seen from the visualizations in Figure 3, which show that the radial extent and azimuthal wavelength decrease with $\mathrm{Ek}$. To achieve this, grid refinement is introduced in the code such that the spectral resolution is increased by 20% whenever under-resolution is detected in the Chebyshev or the spherical harmonic modes. This is done by introducing a threshold for the amplitude ratio between the mode of highest absolute value and of highest wavenumber. When this ratio exceeds the threshold, the code calls a grid refinement subroutine which first deallocates the previous grid, then creates a new one using the new resolution and finally represents the fields from the last continuation step on this new grid by filling in the new modes with zero. This last step is then recomputed with the new resolution and path following continues as intended. In our computations we chose the under-resolution threshold to be between ${10}^{-6}$ and ${10}^{-7}$.

Figure 10 shows the result of continuation in $\mathrm{Ek}$. Because the convective threshold obeys a scaling law ${\mathrm{Ra}}_{\mathrm{therm}}^{\mathrm{crit}}\sim {\mathrm{Ek}}^{-4/3}$[2], we fixed $\mathrm{Ra}$ to be proportional to ${\mathrm{Ek}}^{-1/3}$ so as to remain at approximately the same distance from the convective threshold. This branch showed three coexisting solutions over the range $4\times {10}^{-4}\lesssim \mathrm{Ek}\lesssim 6\times {10}^{-5}$. Our grid refinement algorithm increased the resolution in r by about 20% and by 50% in $\theta$ and $\phi$. We emphasize that most continuations showed monotonic behavior and we have chosen deliberately to present those that did not.

5. Timing Comparisons

Newton’s method is very fast, typically requiring 3-10 iterations. The bottleneck in applying Newton’s method to systems of partial differential equations is always the solution of the linear systems. In our case, this is quantified by the total number of GMRES iterations necessary to compute a new state along the branch. We find that the number of iterations depends on the Ekman number, but is fairly independent of the Rayleigh number. We therefore average the number of GMRES iterations over all the points computed in a branch for a fixed Ekman number. We do this for the Newton method with both implicit and explicit Coriolis to produce the curves in Figure 11, Figure 12, and Figure 13.

Measuring the economy realized in timestepping is more problematic. Because explicit integration effectively approximates a decaying exponential by a polynomial, it displays artificial temporal divergence if the timestep is too large. Because implicit integration instead approximates the exponential by a decaying rational function, the timestep is not constrained by stability. The timestep is, however, still constrained by accuracy. A particularly demanding criterion for judging accuracy is the wavespeed [34]. We compute the relative errors in wavespeeds $C_{\mathrm{error}}\equiv (C_{\Delta t}-C_{\mathrm{exact}})/C_{\mathrm{exact}}$, where $C_{\mathrm{exact}}$ is obtained from Newton’s method and hence has no timestepping error. Figure 14 presents the relative errors as a function of $\mathrm{Ek}$ and $\Delta t$, with the other parameters set to the values below the figure. The errors obtained from the explicit and implicit methods are the same for the same value of $\Delta t$. For each $\mathrm{Ek}$, there is a minimum $\Delta t$ above which the explicit simulation diverges in time, but the implicit simulation still yields a solution with an error that is less than 5%. This value of $\Delta t$ is indicated in Figure 14 by the left endpoint of an arrow for each Ekman number. But eventually, $\Delta t$ is so large that the result of the implicit method can even be qualitatively incorrect. As an upper limit for the implicit method, we choose that the relative error on $C_{\mathrm{error}}$ must be less than 5%. The limiting values of large $\Delta t$ are indicated as the right endpoints of the arrows in Figure 14 and in Table 1. Hence, for each $\mathrm{Ek}$, there is a range of $\Delta t$-values that can be used only for implicit timestepping and whose error does not exceed $5\%$. These ranges are given in Table 1. Although these results are specific to the parameter values and flows that we have simulated, they suggest that implicit Coriolis integration can be orders of magnitude faster, enabling the simulation of convection in a spherical annulus at low Ekman numbers.

6. Discussion

We have developed a pair of codes for simulating thermal convection in a rotating spherical fluid shell that rely on an implicit treatment of the Coriolis force. The numerical cost of this improvement is quite manageable: the block diagonal matrix systems which arise from implicit treatment of the diffusive terms must be replaced by block pentadiagonal matrix systems, which can still be solved by block banded LU decomposition and backsolving. Once implemented in a timestepping code, implicit integration with a very large pseudo-timestep ($\Delta t=200$) can be leveraged to precondition the large linear systems that are the core of Newton’s method. When only the diffusive terms are treated implicitly, this is known as Stokes preconditioning [11,12,13,25,31,32]; we could call the method developed here Stokes-and-Coriolis preconditioning.

We have demonstrated this method’s capabilities by carrying out continuation in Rayleigh number at various values of the Ekman number on the order of ${10}^{-5}$ for rotating waves with azimuthal wavenumbers 4, 8, and 12. We have found several intriguing examples of such branches contain plateaus in drift frequency, separated by intervals of rapid change delimited by pairs of saddle-node bifurcations. The physical, or at least phenomenological, reasons for these properties remain to be discovered. We have also implemented continuation in Ekman number, spaced logarithmically, in which we automatically measure and increase the resolution as needed. We are unaware of any previous examples of continuation in Ekman number or with automatic resolution adjustment in the literature.

We have measured the economy that is realized by implicit treatment of the Coriolis force. For Newton’s method, the advantage of the implicit over the explicit treatment is dramatic. For example, the explicit algorithm takes 20 times as many GMRES iterations (and hence 20 times as much CPU time) as the implicit method to compute states on the ${\mathrm{RW}}_{12}$ branch at $\mathrm{Ek}=3\times {10}^{-5}$. For lower Ekman numbers, the explicit method takes an unmanageable number of GMRES iterations, making it impossible to use in practice.

For timestepping, the advantage of the implicit over the explicit treatment is equally spectacular. For computing rotating waves at $\mathrm{Ek}={10}^{-5}$, the implicit algorithm can reasonably use timesteps that are 200 times past the temporal stability limit, meaning that simulations in this regime can be 200 times faster than if an explicit method were used.

We now describe the differences between Stokes preconditioning and the method used in [21,22,23,24] for convection in rotating spheres, as well as by [35,36,37,38] for other hydrodynamic problems. These authors seek steady states or traveling waves as roots of $\mathit{U}(t+T)-\mathit{U}\left(t\right)$, where $\mathit{U}(t+T)$ is computed from $\mathit{U}\left(t\right)$ in the usual manner by carrying out many timesteps, each with a small $\Delta t$; we call this the integration method. In contrast, although equations (19)-(20) would seem to imply that we too seek roots of $\mathit{U}(t+T)-\mathit{U}\left(t\right)$, this is not at all the case. In the Stokes preconditioning method, we compute $\mathit{U}(t+T)$ from $\mathit{U}\left(t\right)$ by carrying out a single very large implicit-explicit Euler timestep (17), with a timestep $\Delta t=T$ that is so large that the difference $\mathit{U}(t+\Delta t)-\mathit{U}\left(t\right)$ no longer approximates the time derivative. Via equation (20), it turns out that this difference is a preconditioned version of the operator whose roots are sought by our method, namely $\mathcal{L}+\mathcal{N}$, preconditioned by $\mathcal{L}$. Both methods solve the linear systems that are the core of Newton’s method via GMRES, which relies on repeated actions of the Jacobian. For the integration method, the action of the Jacobian consists of integrating the linearized equations via many small timesteps, while for us it consists taking a single large linearized timestep. Here, there is a trade-off. Despite the preconditioning displayed by (20), the Jacobian resulting from Stokes preconditioning is less well conditioned than that which results from taking many small linearized timesteps, and so more Jacobian actions are required for GMRES to converge to a solution of the linear system (21). However, each action costs less, since it consists of only a single timestep. This trade-off – number of actions vs. timesteps per action – can be quantified via the total number of timesteps required to compute a steady state or travelling wave. For wall-bounded shear flows in the transitional regime, we have found that the Stokes preconditioning method is faster than the integration method by a factor of 11 for plane Couette flow and by a factor of 35 for pipe flow; see [25] for details. We have not carried out a timing comparison of the integration and the Stokes preconditioning method for convection in rotating spheres.

Although the Stokes preconditioning method is faster the integration method is more general. The Stokes preconditioning method is only capable of computing steady states or traveling waves, while the integration method can compute periodic orbits of all kinds, including, for example modulated rotating waves or pulsing states.

A specific difference between our implementation of implicit Coriolis integration and that of [21,22,23,24,28,39] is that we solve the linear systems arising from implicit treatment of the diffusive and Coriolis terms directly via block pentadiagonal LU decomposition and backsolving, whereas [21,22,23,24,28,39] use Krylov methods to solve these systems. Time-integration using explicit, implicit and semi-implicit treatment of the Coriolis term is compared by [39]. Our strategy is to solve linear systems directly as much as possible and to resort to Krylov only to invert the full Jacobian.

One might wonder about the applicability of this achievement to geophysical and astrophysical research. Small Ekman numbers, like large Reynolds numbers, are associated with chaos and turbulence, not with the regularity and periodicity of traveling waves, which moreover are almost invariably unstable. The Boussinesq and Navier-Stokes equations generally contain a very large number of solution branches, most of them partly or completely unstable [e.g., 33,36–38,40–42]. This profusion of states, sometimes called Exact Coherent Structures, has become the focus of extensive research, motivated by the idea [43,44] that turbulence could be viewed from a dynamical systems perspective as a collection of trajectories ricocheting between periodic orbits along their unstable directions, a kind of deterministic analogue of ergodic theory.

We now turn to the relevance of implicit Coriolis timestepping for geophysics. Implicit Coriolis timestepping not only speeds up calculations of stable and realistic flows, it also allows us to take much larger timesteps that yield inaccurate results. Constraints on timesteps (whether for stability or accuracy) arise from different physical forces. In rapidly rotating fluids, the Coriolis force generates inertial waves, which have a large range of frequencies and are continually generated and damped. Using a timestep that greatly exceeds the constraint for resolving inertial waves can be considered analogous to using the incompressible Navier-Stokes equations, which deliberately neglect the high-frequency sound waves [14,15].