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

Juan Cruz Gonzalez Sembla, Camille Rambert, Fred Feudel, Laurette S Tuckerman

► To cite this version:

Juan Cruz Gonzalez Sembla, Camille Rambert, Fred Feudel, Laurette S Tuckerman. Convection in a Rapidly Rotating Spherical Shell: Newton's Method Using Implicit Coriolis Integration. Mathematics , 2025, 13 (13), pp.2113. 10.3390/math13132113 . hal-05368782

HAL Id: hal-05368782

https://cnrs.hal.science/hal-05368782v1

Submitted on 17 Nov 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Cite as: J.C. Gonzalez Sembla, C. Rambert, F. Feudel, L.S. Tuckerman, Mathematics (MDPI) 13, 2113 (2025).

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

Juan Cruz Gonzalez Sembla¹, Camille Rambert ¹, Fred Feudel², and Laurette S. Tuckerman¹

¹PMMH, CNRS, ESPCI Paris, Université PSL, Sorbonne Université, Université de Paris, 75005 Paris, France; juan-cruz.gonzalez-sembla@espci.fr; camille.rambert@espci.fr; laurette@pmmh.espci.fr

²Institut für Physik und Astronomie, Universität Potsdam, 14476 Potsdam, Germany; ffeudel@uni-potsdam.de

Abstract

Geophysical flows are characterized by rapid rotation. Simulating these flows requires small timesteps to achieve stability and accuracy. Numerical stability can be greatly improved by the implicit integration of the terms that are most responsible for destabilizing the numerical scheme. We have implemented an implicit treatment of the Coriolis force in a rotating spherical shell driven by a radial thermal gradient. We modified the resulting timestepping code to carry out steady-state solving via Newton's method, which has no timestepping error. The implicit terms have the effect of preconditioning the linear systems, which can then be rapidly solved by a matrix-free Krylov method. We computed the branches of rotating waves with azimuthal wavenumbers ranging from 4 to 12. As the Ekman number (the non-dimensionalized inverse rotation rate) decreases, the flows are increasingly axially independent and localized near the inner cylinder, in keeping with well-known theoretical predictions and previous experimental and numerical results. The advantage of the implicit over the explicit treatment also increases dramatically with decreasing Ek, reducing the cost of computation by as much as a factor of 20 for Ekman numbers of order of $10^{-5}$. We carried out continuation for both the Rayleigh and Ekman numbers and obtained interesting branches in which the drift velocity remained unchanged between pairs of saddle-node bifurcations.
# I. INTRODUCTION

Convection in a rapidly rotating spherical shell subjected to differential heating or imposed heat flux and a radial gravity field is a standard model for planetary interiors. The onset of convection in this configuration has been extensively studied [1–8]. Due to the symmetry of the configuration, a primary Hopf bifurcation leads to traveling (rotating) waves [9–11], which take the form of fluid columns perpendicular to the axis of rotation, sometimes referred to as thermal Rossby waves [12–15]. Scaling laws have been derived for the critical Rayleigh number and for the azimuthal wavenumber for these patterns and their drifting frequencies [6] as a function of the Ekman number, which is the ratio of the viscous to the Coriolis force. A secondary Hopf bifurcation leads to modulated rotating waves, i.e., quasi-periodic patterns with an additional modulation frequency [16–20].

Although convection in astrophysics follows the same basic fluid-dynamical principles as convection in geophysics, our work is primarily relevant to geophysical applications, in which the density variations are sufficiently small that the fluid can be treated as Boussinesq, the convection is due to temperature (rather than concentration) gradients, the domain is a spherical shell, and the Coriolis force is dominant. Concerning the physics of the problem, the balance of forces can be used to define various regimes, each with a different convective length and velocity scales. In particular, the Coriolis-inertial-buoyancy balance can be contrasted to a regime called viscous balance in which the viscous force replaces the inertial force [21–24]. Experimental and numerical studies have analyzed the effect of rotation on convection and on the generation of a magnetic field; these are reviewed in e.g., [12, 15, 23, 25]. Some notable numerical papers incorporating these effects are [26], which describes a pioneering 3D code for solving the anelastic magnetohydrodynamic equations, and [27] and [28], which compare the performance of various codes on Boussinesq and anelastic convection-driven dynamo benchmark problems, respectively (the anelastic approximation lies between the Boussinesq approximation and the fully compressible equations, capturing some compressibility effects while filtering out fast sound waves). Laboratory experiments on radial gravity in a sphere are difficult to perform because they are affected by Earth’s vertical gravity. For this reason, microgravity experiments called Geoflow and Atmoflow, respectively, have been placed on a space station to mimic radial convection in rotating spherical annuli within the Earth [29, 30] and in the atmosphere [31].

Three-dimensional numerical simulations in a spherical geometry usually represent the velocity and magnetic fields in terms of poloidal and toroidal potentials, which reduce a 3D solenoidal field to two scalar fields. These potentials and the temperature are decomposed into spherical harmonics in the angular directions, with either Chebyshev polynomials or finite differences used in the radial direction. Chebyshev collocation methods lead to dense matrices, but sparse matrices can be obtained when using a Galerkin approach [32, 33], resulting in lower time and memory requirements.

Rapid rotation leads to a large Coriolis force. The Coriolis force introduces coupling (of the velocity components and/or the spherical harmonics) and so it is usually treated explicitly in numerical codes. This leads to a significant limitation on the timestep. A natural approach is then to treat this term implicitly, leading to a more stable scheme; its implementation for finite difference or Chebyshev collocation discretizations in the radial direction is discussed in [34, 35], respectively. The implicit treatment of the Coriolis term has been successfully implemented by [36–39] in a planar geometry and by [20, 33, 40, 41] in a spherical geometry.
Another form of numerical computation is continuation, i.e., the computation of steady states, traveling waves, or other periodic states independently of their stability via Newton's or a related method. Although unstable solutions cannot be observed experimentally, they are interesting because they form a framework for understanding the stable solutions, which in turn are accessible to experiment with.

In this work, 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. We will illustrate, test, and compare the codes on the thermally driven rotating waves. We carried out numerical continuation in both Rayleigh and Ekman numbers. We found interesting examples of branches which contain saddle-node bifurcations separating plateaus in drift frequency. Our numerical algorithm uses the basic framework described in [18, 19, 35, 42, 43] and the spherical harmonic library of [44]. The research which is closest to ours is that of [20, 45–47], 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 [48] and which we will discuss after we have described our numerical methods.

In Section II, we describe the problem, governing equations and non-dimensional parameters in addition to introducing rotating wave solutions of different azimuthal symmetries. The numerical methods, including spatial discretization, Newton, and path-following methods are detailed in Section III. Continuation results are described in Section IV, while timing comparisons and timestepping results are presented in Section V. We discuss and provide a conclusion to our results in Sections VI and VII.

## II. PHYSICAL DESCRIPTION

### A. Governing Equations

We study classical Rayleigh-Bénard convection in a spherical fluid shell rotating with constant angular velocity $\Omega$ about the $z$ axis, as illustrated in Figure 1. 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_{\text{out}}$ and $r_{\text{in}} = r_{\text{out}} - 1$, respectively. The aspect ratio can be specified via $r_{\text{out}}$ or via $\eta = r_{\text{in}}/r_{\text{out}} = (r_{\text{out}} - 1)/r_{\text{out}}$. Time is scaled by the viscous diffusion time $d^2/\nu$, where $\nu$ is the kinematic viscosity and the velocity 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 non-dimensionalize pressure by $\rho_0\nu$, and temperature by $\Delta T$.

The gravitational acceleration $g$ is assumed to be proportional to the distance $r$ from the center of the sphere (as is the case for self-gravity in the interior of a spherical body with constant mass density) and is thus of the form $g = -(g_o/r_{\text{out}})r$, where $g_o$ is the acceleration at radius $r_{\text{out}}$. The resulting non-dimensional Boussinesq equations read
<figure><img src="image_2.png" /><figcaption>FIG. 1: Geometry. The domain is the spherical annulus between an inner and outer sphere of radius $r_{in}$ and $r_{out}$ whose temperatures are fixed at $T = 1$ and $T = 0$. Both spheres rotate around the $z$ axis with angular velocity $\Omega$. Figure taken and adapted from [49].</figcaption></figure>

$$ \mathrm{Ek} \left[ \frac{\partial \mathbf{U}}{\partial t} + (\mathbf{U} \cdot \nabla) \mathbf{U} - \nabla^2 \mathbf{U} \right] = -\nabla P + \mathrm{Ra} T \frac{\mathbf{r}}{r_{\mathrm{out}}} - 2 \mathbf{e}_z \times \mathbf{U} \quad (1a) $$

$$ \frac{\partial T}{\partial t} + \mathbf{U} \cdot \nabla T = \frac{1}{\mathrm{Pr}} \nabla^2 T \quad (1b) $$

$$ \nabla \cdot \mathbf{U} = 0 \quad (1c) $$

When the density is constant, the gravitational and centrifugal forces can be written as gradients and included in the pressure gradient. Variable density leads to a non-gradient portion of these forces. The non-gradient part arising from gravity drives the convection. We make the common assumption in the geophysics of neglecting the non-gradient part arising from the centrifugal force [15, 47]; see [15, 50, 51] for examples of the effect of centrifugal buoyancy.

The non-dimensional parameters used in Equation (1a) above are the Ekman, Prandtl, and Rayleigh numbers:

$$ \mathrm{Ek} \equiv \frac{\nu}{\Omega d^2}, \quad \mathrm{Pr} \equiv \frac{\nu}{\kappa}, \quad \mathrm{Ra} = \mathrm{Ra}_{\mathrm{rot}} \equiv \frac{d\alpha g_o \Delta T}{\Omega \nu} \quad (2) $$

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

$$ Ra_{\mathrm{therm}} \equiv \frac{d^3 \alpha g_0 \Delta T}{\kappa \nu} = \frac{d\alpha g_o \Delta T}{\Omega \nu} \frac{\Omega d^2}{\nu} \frac{\nu}{\kappa} = Ra_{\mathrm{rot}} \frac{\mathrm{Pr}}{\mathrm{Ek}} \quad (3) $$
In what follows, we will denote Ra$_{rot}$ merely by Ra. No-slip and fixed-temperature boundary conditions are applied at the inner and outer radii:

$$ U(r_{\text{in}}) = 0 \quad U(r_{\text{out}}) = 0 \qquad (4a) $$

$$ T(r_{\text{in}}) = 1 \quad T(r_{\text{out}}) = 0 \qquad (4b) $$

and the conductive solution is

$$ U_{\text{cond}} = 0, \quad T_{\text{cond}}(r) = \frac{r_{\text{out}} r_{\text{in}}}{r} - r_{\text{in}} \qquad (5) $$

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

## B. Overview

The first states to appear at the onset of convection are rotating waves. We denote these by $RW_M$, where $M$ is the azimuthal wavenumber of the rotating wave. Figure 2a,b, modeled on those in [18], display properties of some of these rotating wave solutions for $Ek = 10^{-3}$. Four traveling waves are shown, with azimuthal wavenumbers from $M = 2$ to $M = 5$. At $Ek = 10^{-3}$, $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 $\Omega$, decreases as Ra increases, from prograde (faster than $\Omega$) to retrograde (slower than $\Omega$). These results were obtained via timestepping simulations starting from the initial conditions of the form $\cos(M\varphi)$.

Figure 3 shows the qualitative evolution of the rotating waves as 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: the radial interval over which the convection is most active decreases, and the azimuthal wavelength decreases as well [6, 41]. $RW_8$ appeared naturally when performing timestepping simulations from an initial condition of the form $\cos(4\varphi)$ at $Ek = 10^{-4}$, as did $RW_{12}$ at $Ek = 3 \times 10^{-5}$. This result was then used as input for a path-following computation to obtain $RW_{12}$ at $Ek = 10^{-5}$. The same techniques were applied to obtain $RW_4$ at $Ek = 10^{-2}$ starting from $RW_4$ at $Ek = 10^{-3}$. Figure 4 gives a more detailed visualization of an $RW_4$ state.

# III. NUMERICAL METHODS

## A. Spatial Representation

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

$$ U(r, \theta, \varphi, t) = \nabla \times e \, e_r + \nabla \times \nabla \times f \, e_r \qquad (6) $$

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

$$ f(r, \theta, \varphi, t) = \sum_{k,l,m} f_{k,l,m}(t) T_k(x) Y_l^m(\theta, \varphi) \qquad (7) $$
<figure><img src="image_3.png" /><figcaption>FIG. 2: (a) Bifurcation diagram of convection in a rotating spherical annulus for Ek = 10⁻³. Branches RW₂ (blue), RW₃ (black), RW₄ (red), and RW₅ (green) are shown. (b) Kinetic energy density. (c) Drift frequency. The drift frequencies for each rotating wave decreases from prograde (faster than imposed velocity Ω) to retrograde (slower than Ω) with increasing Ra with the same slope. The resolution used for these simulations was (N<sub>r</sub>, N<sub>θ</sub>, N<sub>φ</sub> × M) = (46, 72, 128 × 1).</figcaption></figure>

To carry out transformations to and from the spherical harmonic space, we use the SHTns library [44]. 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 [52] for an 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. This library has promising scalability prospects for computations in massively parallel clusters, as demonstrated in the XSHELLS and Magic codes [53, 54]. Furthermore, our code agrees with the non-magnetic benchmark proposed in [27] and the results presented in [18], validating our numerical method as well as our choice of using the SHTns library.

De-aliasing in the *r* direction is carried out according to the 3/2 rule, whereas in the *θ* and *φ* directions, it is performed by SHTns. To achieve this, the library determines a suitable number of physical points based on the nonlinear 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π/M, we restrict our calculations to a segment of the sphere 0 ≤ *φ* < 2π/M, or, equivalently, to those whose Fourier expansions contain only multiples of the wavenumber *M*; this is easily implemented using options in SHTns. Lastly, we will denote the number of gridpoints in *r*, *θ* and *φ* as N<sub>r</sub>, N<sub>θ</sub> and N<sub>φ</sub> × *M*, respectively. The resolutions chosen for the simulations are similar to those that have been shown to be adequate in the benchmark paper [27] and in more recent studies [7, 20, 55]. The spatial resolutions are 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,
<figure><img src="image_4.png" /><figcaption>FIG. 3: Velocity fields of rotating waves as the Ekman number is reduced. The above shows RW₄ for Ek = 10⁻² and Ra = 140; RW₄ for Ek = 10⁻³ and Ra = 77; RW₈ for Ek = 10⁻⁴ and Ra = 88; and RW₁₂ for Ek = 10⁻⁵ and Ra = 126. A state with azimuthal wavenumber M has M convection cells, forming M Taylor columns. All of these states have an approximate kinetic energy density of 25. For each field, the radial velocity Uₜ on the equatorial plane (top) and the azimuthal velocity Uφ on a meridional plane (bottom) are shown. As Ek decreases, the convective region becomes more confined near the inner sphere and the flow becomes almost axially independent, according to the Taylor-Proudman theorem. The resolutions used for these simulations were (Nₜ, Nφ, Nφ × M) = (46, 72, 32 × 4) for RW₄ at both Ek = 10⁻² and Ek = 10⁻³; (60, 100, 28 × 8) for RW₈; and (90, 184, 32 × 12) for RW₁₂.</figcaption></figure>

leading to

$$ \sum_{l,m} \frac{l(l+1)}{r^2} \mathrm{Ek} \left( \frac{\partial}{\partial t} - L_l \right) e_{l,m}(r,t) Y_l^m(\theta, \varphi) = e_r \cdot \nabla \times \mathbf{F} \quad (8a) $$

$$ -\sum_{l,m} \frac{l(l+1)}{r^2} \mathrm{Ek} \left( \frac{\partial}{\partial t} - L_l \right) L_l f_{l,m}(r,t) Y_l^m(\theta, \varphi) = e_r \cdot \nabla \times \nabla \times \mathbf{F} \quad (8b) $$

where

$$ L_l \equiv \frac{\partial^2}{\partial r^2} - \frac{l(l+1)}{r^2} \quad \text{and} \quad \mathbf{F} \equiv \mathrm{Ek} (\mathbf{U} \cdot \nabla)\mathbf{U} - \nabla P + \frac{Ra}{r_{\mathrm{out}}} T \mathbf{r} - 2e_z \times \mathbf{U} \quad (9) $$

Equations (8a) and (8b) are decoupled in l, m and in potential (e vs. f), and so their discrete versions involve block-diagonal matrices, with an Nₜ × Nₜ block for each (l, m) pair, as in Figure 5(a). Boundary conditions (4a) on U become

$$ e_{lm} = f_{lm} = df_{lm}/dr = 0 \quad \text{at } r_{\mathrm{in}}, r_{\mathrm{out}} \quad (10) $$
<figure><img src="image_5.png" /><figcaption>FIG. 4: Visualizations of the components of RW$_4$ at $Ek = 10^{-3}$ and $Ra = 140$. **Top:** equatorial plane. **Bottom:** meridional plane. The resolution used here was $(N_r, N_\theta, N_\varphi \times M) = (46, 72, 32 \times 4)$.</figcaption></figure>

The number of boundary conditions on $e$ and $f$ matches the fact that Equation (8a) for $e$ is of second order in $r$ while Equation (8b) for $f$ is of fourth order in $r$. The boundary conditions are imposed by replacing the last two rows of Equation (8a) or the last four rows of Equation (8b), corresponding to the highest Chebyshev polynomials. See, for example, [35, 56].

The temperature field is also expanded in Chebyshev polynomials and spherical harmonics like Equation (7). Its evolution is governed by the discretized form of Equation (1b) and the boundary conditions in Equation (4b) on the temperature are imposed straightforwardly on the rows corresponding to the last two Chebyshev polynomials.

## B. Implicit Coriolis Integration

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

$$ \sum_{m,l} \left[ \frac{l(l+1)}{r^2} \mathrm{Ek} \left( \frac{\partial}{\partial t} - L_l \right) - \frac{2im}{r^2} \right] e_{lm}(r,t) Y_l^m(\theta, \varphi) \\ + \frac{2}{r^2} \left( \frac{l(l+1)}{r} - \frac{\partial}{\partial r} \right) f_{lm}(r,t) \sin\theta \frac{d}{d\theta} Y_l^m(\theta, \varphi) \\ + 2 \frac{l(l+1)}{r^2} \left( \frac{2}{r} - \frac{\partial}{\partial r} \right) f_{lm}(r,t) \cos\theta Y_l^m(\theta, \varphi) = e_r \cdot \nabla \times \mathbf{F}^{\text{implicit}} \quad (11a) $$
$$
\begin{aligned}
& \sum_{m,l} -\left[ \frac{l(l+1)}{r^2} \text{Ek} \left( \frac{\partial}{\partial t} - L_l \right) - \frac{2im}{r^2} \right] L_l f_{lm}(r,t) Y_l^m(\theta, \varphi) \\
& + \frac{2}{r^2} \left( \frac{l(l+1)}{r} - \frac{\partial}{\partial r} \right) e_{lm}(r,t) \sin \theta \frac{d}{d\theta} Y_l^{|m|}(\theta, \varphi) \\
& + 2 \frac{l(l+1)}{r^2} \left( \frac{2}{r} - \frac{\partial}{\partial r} \right) e_{lm}(r,t) \cos \theta Y_l^m(\theta, \varphi) = e_r \cdot \nabla \times \nabla \times \mathbf{F}^{\text{implicit}}
\end{aligned}
\qquad (11b)
$$

where

$$
\mathbf{F}^{\text{implicit}} \equiv \text{Ek} (\mathbf{U} \cdot \nabla) \mathbf{U} - \nabla P + \frac{\text{Ra}}{r_{\text{out}}} T \mathbf{r} \quad (12)
$$

The recursion relations of the Legendre polynomials

$$
\sin \theta \frac{d}{d\theta} P_l^{|m|} = \frac{l(l - |m| + 1)}{2l + 1} P_{l+1}^{|m|} - \frac{(l + |m|)(l + 1)}{2l + 1} P_{l-1}^{|m|} \quad (13a)
$$

$$
\cos \theta P_l^{|m|} = \frac{l - |m| + 1}{2l + 1} P_{l+1}^{|m|} + \frac{(l + |m|)}{2l + 1} P_{l-1}^{|m|} \quad (13b)
$$

together with

$$
Y_l^m(\theta, \varphi) = N_l^m P_l^m(\cos \theta) e^{im\varphi} \qquad N_l^m = \sqrt{(2 - \delta_{m0})(2l + 1) \frac{(l-m)!}{l+m}!} \quad (14)
$$

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

$$
\begin{aligned}
& \sum_{m,l} Y_l^m(\theta, \varphi) \left\{ \left[ \frac{l(l+1)}{r^2} \text{Ek} \left( \frac{\partial}{\partial t} - L_l \right) - \frac{2im}{r^2} \right] e_{lm}(r,t) \right. \\
& \qquad \left. + \frac{2}{r^2} \frac{(l-1)(l+1)(l-|m|)}{2l-1} \left[ \frac{l}{r} - \frac{\partial}{\partial r} \right] \frac{N_{l-1}^m}{N_l^m} f_{l-1,m}(r,t) \right. \\
& \qquad \left. + \frac{2}{r^2} \frac{l(l+2)(l+1+|m|)}{2l+3} \left[ -\frac{(l+1)}{r} - \frac{\partial}{\partial r} \right] \frac{N_{l+1}^m}{N_l^m} f_{l+1,m}(r,t) \right\} = e_r \cdot \nabla \times \mathbf{F}^{\text{implicit}} \quad (15a)
\end{aligned}
$$

$$
\begin{aligned}
& \sum_{m,l} Y_l^m(\theta, \varphi) \left\{ -\left[ \frac{l(l+1)}{r^2} \text{Ek} \left( \frac{\partial}{\partial t} - L_l \right) - \frac{2im}{r^2} \right] L_l f_{lm}(r,t) \right. \\
& \qquad \left. + \frac{2}{r^2} \frac{(l-1)(l+1)(l-|m|)}{2l-1} \left[ \frac{l}{r} - \frac{\partial}{\partial r} \right] \frac{N_{l-1}^m}{N_l^m} e_{l-1,m}(r,t) \right. \\
& \qquad \left. + \frac{2}{r^2} \frac{l(l+2)(l+1+|m|)}{2l+3} \left[ -\frac{(l+1)}{r} - \frac{\partial}{\partial r} \right] \frac{N_{l+1}^m}{N_l^m} e_{l+1,m}(r,t) \right\} = e_r \cdot \nabla \times \nabla \times \mathbf{F}^{\text{implicit}} \quad (15b)
\end{aligned}
$$

Appendix A presents further details of these computations. In contrast with
Equations (8a) and (8b), Equations (15a) and (15b) are coupled in several ways:

* The spectral coefficients $e_{lm}$ and $f_{lm}$ both appear in both equations.

* While each *m* can be treated independently, components *l*, *l*+1, and *l*−1 are coupled.
<figure><img src="image_6.png" /><figcaption>FIG. 5: (a) Block-diagonal matrix. (b) Block-pentadiagonal matrix.</figcaption></figure>

* The real and imaginary parts of $e_{lm}$ and $f_{lm}$ are coupled via the imaginary coefficient $2im/r^2$.

However, two decoupled classes of coefficients appear, with 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, in Equation (15a), the real component of $e_{lm}$ is coupled to its imaginary component and to the real parts of $f_{l\pm 1,m}$ (with opposite parity to $l$) and in Equation (15b), the real component of $f_{lm}$ is coupled to its imaginary component and to the real parts of $e_{l\pm 1,m}$. For each class, the sums in Equations (15a) and (15b) are represented by a block pentadiagonal matrix, as illustrated in Figure 5b.

### C. Newton Method

We represent our system of governing equations schematically by

$$ \frac{\partial \mathbf{U}}{\partial t} = \mathcal{N}(\mathbf{U}) + \mathcal{L}\mathbf{U} \qquad (16) $$

For our problem, $\mathbf{U}$ corresponds to the poloidal and toroidal fields $e$ and $f$ and the temperature $T$. For explicit or implicit Coriolis integration, $\mathcal{N}(\mathbf{U})$ corresponds to the right-hand-sides of Equations (8a) and (8b) or of Equations (11a) and (11b), respectively. $\mathcal{L}\mathbf{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 $\mathbf{U}$.)

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

$$ \begin{aligned} \mathbf{U}(t + \Delta t) &= \mathbf{U}(t) + \Delta t(\mathcal{N}(\mathbf{U}(t)) + \mathcal{L}\mathbf{U}(t + \Delta t)) \\ \Rightarrow \mathbf{U}(t + \Delta t) &= \mathcal{B}(\mathbf{U}(t)) \equiv (I - \Delta t\mathcal{L})^{-1}(I + \Delta t\mathcal{N})\mathbf{U}(t) \end{aligned} \qquad (17) $$

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

$$ 0 = \mathcal{A}(\mathbf{U}) \equiv \mathcal{N}(\mathbf{U}) + \mathcal{L}\mathbf{U} \qquad (18) $$
or, alternatively

$$0 = U(t + \Delta t) - U(t) \qquad (19)$$

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

$$ \begin{aligned} (\mathcal{B} - I)(U(t)) = U(t + \Delta t) - U(t) &= (I - \Delta t\mathcal{L})^{-1}(I + \Delta t\mathcal{N})U(t) - U(t) \\ &= (I - \Delta t\mathcal{L})^{-1}(I + \Delta t\mathcal{N} - (I - \Delta t\mathcal{L}))U(t) \\ &= \underbrace{(I - \Delta t\mathcal{L})^{-1}\Delta t}_{\mathcal{P}(\Delta t)} \underbrace{(\mathcal{N} + \mathcal{L})}_{\mathcal{A}} U(t) \end{aligned} \qquad (20) $$

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

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

$$(\mathcal{B}_U - I)\mathbf{u} = (\mathcal{B} - I)\mathbf{U} \qquad (21)$$

where $U$ is the current estimate and $\mathcal{B}_U$ is the Jacobian of $\mathcal{B}$ evaluated at $U$. The linear Jacobian operator $\mathcal{B}_U$ acting on $\mathbf{u}$ is derived from $\mathcal{B}$ by modifying the nonlinear operator $\mathcal{N}$. We recall that $U$ represents $(U, T)$ and $\mathbf{u}$ represents $(\mathbf{u}, \tau)$ where $\tau$ is a temperature perturbation. To form $\mathcal{N}_U$, we replace the nonlinear terms in Equations (1a) and (1b) as follows:

$$(U \cdot \nabla)U \implies (U \cdot \nabla)\mathbf{u} + (\mathbf{u} \cdot \nabla)U \qquad (22a)$$

$$U \cdot \nabla T \implies U \cdot \nabla \tau + \mathbf{u} \cdot \nabla T \qquad (22b)$$

These substitutions produce $\mathcal{N}_U$, which is then used in Equation (17) to form $\mathcal{B}_U$. In addition, homogeneous boundary conditions must be applied to the perturbations $\mathbf{u}$ and $\tau$ instead of the inhomogeneous boundary conditions Equation (4). We emphasize that the operations on both sides of Equation (21) are just the ordinary and the linearized timestepping operators $\mathcal{B}$ and $\mathcal{B}_U$.

We solve Equation (21) using a Krylov method such as GMRES [57]. In such matrix-free methods, one only needs to provide the right-hand-side $(\mathcal{B} - I)U$ for the current estimate, and a routine which carries out the action of $(\mathcal{B}_U - I)$ on any vector $\mathbf{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}_U - I) = \mathcal{P}(\Delta t)\mathcal{A}$ is well conditioned; it is for this reason that we take $\Delta t$ large. Once the decrement $\mathbf{u}$ is determined by solving Equation (21), it is subtracted from $U$ to form an improved estimate. We accept $\mathbf{u}$ as the decrement if the linear system Equation (21) is solved by GMRES to relative accuracy $10^{-10}$. For Newton's method, we accept $U$ as a steady state, i.e., a solution to Equation (19), if $\| (\mathcal{B} - I)U \| < 10^{-7}$. Similar values of this tolerance have been used
in [20]. In addition, solutions found via Newton's method are far more accurate than those found by timestepping, which is usually halted when a solution remains constant to three to five significant digits.

The method described above is called Stokes preconditioning. It was first applied to calculate bifurcation diagrams in spherical Couette flow [58, 59] and subsequently to a wide variety of hydrodynamic problems, e.g., [60–67] and even to Bose-Einstein condensation [68, 69]. In Stokes preconditioning, $\Delta t$ no longer plays the role of a timestep but serves to better condition the system of linear equations in each Newton step. We investigated the effect of the choice of $\Delta t$ on Newton solving with Stokes preconditioning in [48, 58] and on the computation of eigenvalues by a similar method in [70]. In these articles, time has been non-dimensionalized by the advective time, and the usual value used in timestepping is $\Delta t \approx 0.01$. An increase in the value of $\Delta t$ from 0.01 to 0.1 and then to 10 leads to a decrease in the required CPU-time for convergence GMRES by 1–2 orders of magnitude per $\Delta t$ decade until an asymptotic limit of $\Delta t \approx 100$ or 1000 is reached. We emphasize that the final result is independent of the choice of $\Delta t$: only the time taken for convergence is affected. In the current study, we use the value $\Delta t = 200$ for Newton solving.

The name of the method comes from the fact that the linear operator $\mathcal{L}$ which is integrated implicitly in Equation (17) is usually the viscous term that occurs in Stokes equation. This preconditioning becomes less powerful as the Reynolds number increases and other terms begin to dominate the viscous term; this effect was demonstrated and measured quantitatively in [48]. Hence, 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 this reason that we have included the Coriolis force in $\mathcal{L}$, using the equations derived in Section III B. In Section V, 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.

## D. Traveling Waves

Newton's method can also be used to compute traveling waves in the same way. Azimuthal rotating waves satisfy $U(\varphi, t) = \tilde{U}(\varphi - Ct) = \tilde{U}(\tilde{\varphi})$, where $C$ is an unknown wavespeed. Thus, $\partial_t U = -C\partial_{\tilde{\varphi}}\tilde{U}$. Substituting into Equation (16) and dropping the tildes leads to

$$0 = \mathcal{N}(U) + C\partial_{\varphi}U + \mathcal{L}U \qquad (23)$$

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

$$U(t + \Delta t) = (I - \Delta t\mathcal{L})^{-1} (I + \Delta t(\mathcal{N} + C\partial_{\varphi})) U(t) \qquad (24)$$

If $U$ is expressed in terms of spherical harmonics in which the $\varphi$ dependence is trigonometric, the action of $C\partial_{\varphi}$ on a Fourier component $U_m$ is merely $C\partial_{\varphi}U_m = CimU_m$.

To use Newton's method to determine the unknown field $U$ and wavespeed $C$, we substitute $U \rightarrow U - u$ and $C \rightarrow C - c$ into Equation (23):

$$0 = \mathcal{N}(U - u) + (C - c)\partial_{\varphi}(U - u) + \mathcal{L}(U - u) \qquad (25)$$

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

$$(\mathcal{N}_U + C\partial_{\varphi} + \mathcal{L})u + c\partial_{\varphi}U = (\mathcal{N} + C\partial_{\varphi} + \mathcal{L})(U) \qquad (26)$$
The preconditioner remains $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 Equation (26).

$$P(\Delta t)[(N_U + C\partial_{\varphi} + \mathcal{L})u + c\partial_{\varphi}U] = P(\Delta t)(N + C\partial_{\varphi} + \mathcal{L})(U). \quad (27)$$

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 $u$ be zero, i.e., that the corresponding phase component (whose index we shall call $J$) of $U$ remains unchanged.

This simple choice suggests a trick for retaining the size of the unknown $u$ rather than using augmented fields, since $u_J$ is no longer an unknown. The routine which acts on $(u, c)$ is defined such that $c$ is stored in $u_J$. At the beginning of an action, $c$ is extracted, $u_J$ is set to zero, and the explicit part of the action on $u$ is carried out. Then, we multiply the stored value of $\partial_{\varphi}U$ by $c$ and add the result to the explicit part of the action. When the Krylov method converges and returns $u$, we must again extract $c$ from $u_J$, after which $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, $U_J$ and $u_J$ are no longer unknowns.

## E. 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 the Rayleigh number, we increment or decrement this parameter according to the number of Newton iterations required for convergence in the previous step as follows [64]:

$$Ra^{(i+1)} = Ra^{(i)} + \Delta Ra = Ra^{(i)} + \alpha(Ra^{(i)} - Ra^{(i-1)}) \quad (28)$$

where $Ra^{(i)}$, $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 = \frac{N^{\text{opt}} + 1}{N^{(i)} + 1} \quad (29)$$

$N^{\text{opt}}$ is the target number of Newton iterations. If $N^{(i)} = N^{\text{opt}}$, then $\Delta Ra$ remains unchanged, whereas $\Delta Ra$ is reduced (increased) if $N^{(i)} > (<)N^{\text{opt}}$. The choice of $N^{\text{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 Ra, but each calculation will be faster. For our computations, we fixed $N^{\text{opt}}$ to be between 4 and 6.

To choose an initial estimate $U$ for the next solution along a branch at $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 Ra⁽^i⁾, Ra⁽^i⁻¹⁾, Ra⁽^i⁻²⁾ and the new Rayleigh number Ra⁽^i+¹⁾ from Equation (28) to determine coefficients a, b, c such that

$$Ra^{(i+1)} = aRa^{(i)} + bRa^{(i-1)} + cRa^{(i-2)}. \qquad (30a)$$

We then use a, b, c to form a new estimate of the solution by setting

$$U^{(i+1)} = aU^{(i)} + bU^{(i-1)} + cU^{(i-2)} \qquad (30b)$$

Newton iterations are then used to refine U until the norm of Equation (23) is 10⁻⁷. This procedure requires saving the three previous solution vectors and Rayleigh numbers.

## F. Turning Points

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

$$\left| \frac{\Delta U_J}{U_J} \right| \equiv \left| \frac{U_J^{(i)} - U_J^{(i-1)}}{U_J^{(i)}} \right| \quad \text{with} \quad \gamma \left| \frac{\Delta Ra}{Ra} \right| \equiv \gamma \left| \frac{Ra^{(i)} - Ra^{(i-1)}}{Ra^{(i)}} \right| \qquad (31)$$

where J is the index denoting the element of the solution vector U⁽^i⁾ of the highest absolute value and γ 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 (μ, x) = (0, 0) is ˙x = μ − x², dependent variables (x, U) vary like the square root of the control parameter (μ, Ra − Raₑ) near a turning point. Consequently, as the current solution approaches a turning point, the relative difference |ΔU_J/U_J| increases until it is greater than γ|Δ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 ΔRa to change sign and the continuation to go around the turning point. We replace Equations (28) and (30a) by analogous equations in U_J to determine U_J⁽^i+¹⁾ and a, b, c, and then replace Equation (30) by

$$U_j^{(i+1)} = aU_j^{(i)} + bU_j^{(i-1)} + cU_j^{(i-2)} \quad \text{for } j \neq J, \text{ and} \qquad (32a)$$

$$Ra^{(i+1)} = aRa^{(i)} + bRa^{(i-1)} + cRa^{(i-2)} \qquad (32b)$$

to set the estimate of the new solution. Note that we have only changed the method of determining the initial guess. Ra remains fixed to Ra⁽^i+¹⁾ during the Newton step, while all the elements of U are allowed to vary.

Continuation follows in this manner until eventually |ΔU_J/U_J| exceeds γ|ΔRa/Ra|, at which point the code switches back to using the Rayleigh number as the independent variable. The parameter γ can be obtained by analyzing the behavior of the continuation near a turning point. Indeed, if γ is too large, then the code will continue to use 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$. Although the parameter $\gamma$ is necessary for comparing $|\Delta U_j/U_j|$ and $|\Delta Ra/Ra|$ and thus for determining when to switch continuation variables, its value has no influence on the location of the turning point or of the solution branch. The choice of $\gamma$ only affects the speed of continuation, as we determined by varying its value near a turning point. See also [64].

# IV. BRANCH FOLLOWING

## A. Continuation in Rayleigh Number

Our first computations with the implicit Coriolis method validated the new code by comparing the branches it produced for $Ek = 10^{-3}$ with those in [18], which were obtained using an entirely different continuation code. These branches were presented in Figure 2 above. As explained in Section III, to search for solutions having an $M$ azimuthal periodicity, we perform our computations in a fraction of the annulus of azimuthal width $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 \varphi < 2\pi/M$, or to instabilities which break this restriction such as rotating waves with a different $M$.

We then used the continuation code to compute branches of rotating waves at lower values of the Ekman number. Fixing the Ekman number, we carried out continuation in Rayleigh number $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 $RW_4$ branch for $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.

We also computed an $RW_8$ branch for $Ek = 1.26 \times 10^{-5}$, as shown in Figure 7. Indeed, for this value of $Ek$, rotating waves with an azimuthal wavenumber of 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 $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 at least a phenomenological explanation, but we have not yet been able to find one.

## B. Continuation in Ekman Number and in Resolution

We adapted our method to follow branches in Ekman number, varying $Ek$ on a logarithmic scale. We replace $Ra$ with $\log_{10}(Ek)$ in Equation (28) as follows:

$$ \begin{aligned} \log_{10}(Ek^{(i+1)}) &= \log_{10}(Ek^{(i)}) + \Delta[\log_{10}(Ek)] \\ &= \log_{10}(Ek^{(i)}) + \alpha(\log_{10}(Ek^{(i)}) - \log_{10}(Ek^{(i-1)})) \end{aligned} \quad (33) $$

where the computation of $\alpha$ remains as in Equation (29). In this case, we used $N^{\text{opt}} = 3$ for all of the simulations. Furthermore, in order to go around turning points, we now compare $|\Delta U_j/U_j|$ and $\gamma|\Delta[\log_{10}(Ek)]/\log_{10}(Ek)|$, 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
<figure><img src="image_7.png" /><figcaption>FIG. 6: Bifurcation diagram of rotating wave $RW_4$ for $Ek = 3.53 \times 10^{-5}$ as a function of Ra. (a) Kinetic energy. (b) Drift frequency. The branch contains three long smooth regions of almost constant drift frequency separated by two short intervals (at $Ra \approx 140$ and $Ra \approx 158$) of rapid change containing saddle-node bifurcations. The resolution used for this computation was $(N_r, N_\theta, N_\phi \times M) = (60, 80, 40 \times 4)$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 7: Bifurcation diagram of rotating wave $RW_8$ for $Ek = 1.26 \times 10^{-5}$ as a function of Ra. (a) Kinetic energy. (b) Drift frequency. The branch contains a single interval of rapid change at $Ra \approx 159$. The resolution used for this computation was $(N_r, N_\theta, N_\phi \times M) = (68, 100, 28 \times 8)$.</figcaption></figure>

the left-hand-sides of Equations (8a) and (8b) or Equations (15a) and (15b) must be re-computed every time a new value of $Ek$ is chosen. Indeed, while the Rayleigh number only appears in the right-hand side of these equations, $Ek$, is also present in the diffusion terms, which are treated implicitly regardless of how the Coriolis term is handled. This means that extra work must be performed at the beginning of every continuation step.

Secondly, as the Ekman number is decreased, the critical Rayleigh number for the on-
<figure><img src="image_9.png" /><figcaption>FIG. 8: Bifurcation diagram of the rotating wave $RW_4$ as a function of $Ek$, while fixing $Ra = 6.5 \times Ek^{-1/3}$. (a) Kinetic energy. (b) Drift frequency. Three solution branches exist over the range $4 \times 10^{-5} \lesssim Ek \lesssim 6 \times 10^{-5}$. The resolutions used for this computation ranged from $(N_r, N_\theta, N_\phi \times M) = (68, 92, 48 \times 4)$ to $(80, 136, 70 \times 4)$.</figcaption></figure>

set of convection increases. That is, rotation stabilizes the configuration against convection. A classic result [3] is that the critical value of the usual thermal Rayleigh number $Ra_{\text{therm}}$ varies like $Ek^{-4/3}$, so that the critical value of our rotational Rayleigh number $Ra_{\text{rot}} = Ra_{\text{therm}}Ek/Pr$ varies like $Ek^{-1/3}/Pr$. Therefore, instead of keeping $Ra$ fixed, we have increased $Ra$ every time we decrease $Ek$, so as to remain at approximately the same distance from the convective threshold.

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 decreasing $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 performed 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 procedure, called Fourier interpolation [71, 72], involves no loss of accuracy or stability. The most recent Newton step is then recomputed, using the representation with the finer resolution as an initial condition, and the path following continues as intended. Most of our runs were carried out with an under-resolution threshold (ratio of the highest-amplitude mode to that of the highest wavenumber) of $10^{-6}$. Decreasing this threshold to $10^{-7}$ led to minimal gains in accuracy at a high cost in CPU time.

Figure 8 shows the result of continuation in $Ek$. (Recall that $Ra$ is kept proportional to $Ek^{-1/3}$.) This branch showed three coexisting solutions over the range $4 \times 10^{-5} \lesssim Ek \lesssim 6 \times 10^{-5}$. Our grid refinement algorithm increased the resolution by about 20% in $r$ and by 50% in $\theta$ and $\phi$. We emphasize that most continuations showed monotonic behavior and we have deliberately chosen to present those that did not.
V. 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 Figures 9–11 for $RW_4$, $RW_8$, and $RW_{12}$. The ratio of the number of GMRES iterations between the explicit and implicit methods ranges from 2, for $RW_4$ at $Ek = 3 \times 10^{-2}$, to 20 for $RW_{12}$ at $Ek = 3 \times 10^{-5}$. This ratio would surely continue to increase as Ek decreases, but the explicit Coriolis calculation became prohibitively expensive below $Ek = 3 \times 10^{-5}$. Figure 9 also demonstrates that carrying out the calculation in the restricted domain $0 \le \varphi \le 2\pi/M$ requires fewer GMRES iterations (in addition to each iteration being at least $M$ times faster).

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, i.e., numerical instability, 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 [73]. We compute the relative errors in wavespeeds $C_{error} = (C_{\Delta t} - C_{exact})/C_{exact}$, where $C_{exact}$ is obtained from Newton's method and hence has no timestepping error. Figure 12 presents the relative errors as a function of $Ek$ and $\Delta t$, with the other parameters set to the values below the figure. $C_{error}$ obtained from the explicit and implicit methods are the same or very close for the same value of $\Delta t$. $C_{error}$ increases like $(\Delta t)^2$, indicating that both time-integration methods are second-order in time. Indeed, for timestepping, we have not used the first-order scheme given in Equation (17) and used for our Newton algorithm, but a second-order method combining Runge-Kutta and a predictor-corrector scheme, as described in [35].

For each $Ek$, there is a minimum $\Delta t$ above which the explicit Coriolis simulation diverges in time. This value of $\Delta t$ is indicated in Figure 12 for each Ekman number as the left endpoint of an arrow and also in Table I. Above this value of $\Delta t$, implicit Coriolis timestepping must be used. But eventually, $\Delta t$ is so large that the results of the implicit method are too inaccurate to be useful. As an upper limit for the implicit method, we make the arbitrary choice that the relative error on $C_{error}$ must be less than 2%. These limiting values of large $\Delta t$ are indicated as the right endpoints of the arrows in Figure 12 and in Table I. Hence, for each $Ek$, there is a large range of $\Delta t$-values that can be used only for implicit timestepping and whose error does not exceed 2%. For example, for $Ek = 10^{-5}$, the allowable timestep (i 2% error) for the implicit Coriolis method is almost 20 times that at which the explicit method diverges. These ranges are given in Table I. Fitting to the two smaller values of $Ek$, both endpoints of the range are approximately proportional to $Ek$: the explicit method diverges for $\Delta t \gtrsim 0.5 Ek^{1.1}$ while the error for the implicit method remains less than 2% for $\Delta t \lesssim 0.8 Ek^{0.9}$. If a 2% error in the wavespeed is considered to be too large, the data in Figure 12 can be used to determine the largest $\Delta t$ that achieves the desired error, by drawing a horizontal line at the desired level. Although these results are specific to the parameter values and flows that we have simulated, they imply that implicit Coriolis integration can
<figure><img src="image_10.png" /><figcaption>FIG. 9: Total number of matrix-vector actions required by nested Newton-GMRES algorithm to compute $RW_4$ as a function of $Ek$ with explicit (dashed) and implicit (solid) implementation of Coriolis force. An average is taken over a branch of Ra values. The number of actions required by the explicit algorithm is always greater than that required by the implicit algorithm, with the ratio between them increasing from approximately 2 at $Ek = 10^{-2}$ to approximately 9 at $Ek = 10^{-4}$. Using the full domain $0 \le \varphi < 2\pi$ (green) instead of the reduced domain $0 \le \varphi < \pi/2$ (red) increases the number of GMRES iterations, as well as the cost of each iteration. On average, the full domain requires 1.5 times more actions than the reduced domain counterpart. The explicit algorithm is prohibitively time-consuming for $Ek < 10^{-4}$. The resolutions used in this case were $(N_r, N_\theta, N_\varphi \times M) = (46, 72, 128 \times 1)$ for the computations in the full circle and from $(46, 72, 32 \times 4)$ to $(76, 80, 40 \times 4)$ for those in the quarter circle.</figcaption></figure>

be orders of magnitude faster as $Ek$ decreases, enabling the simulation of convection in a spherical annulus at low $Ek$-man numbers.
<figure><img src="image_11.png" /><figcaption>FIG. 10: Total number of matrix-vector actions required by the nested Newton-GMRES algorithm to compute $RW_4$ (red) and $RW_8$ (blue) as a function of $Ek$ with an explicit (dashed) and implicit (solid) implementation of Coriolis force. For clarity, the insets show the full domain, although the computations are carried out in reduced domains $0 \le \varphi < 2\pi/M$. An average is taken over many values of Ra. The number of actions required by the explicit algorithm is always greater than that required by the implicit algorithm, with the ratio between them increasing from approximately 4 at $Ek = 3 \times 10^{-3}$ to 15 at $Ek = 3 \times 10^{-5}$, below which the explicit algorithm is prohibitively time-consuming. The resolutions used for the $RW_8$ computations ranged from $(N_r, N_\theta, N_\varphi \times M) = (64, 100, 28 \times 8)$ to $(76, 124, 32 \times 8)$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIG. 11: Total number of matrix-vector actions required by the nested Newton-GMRES algorithm to compute $RW_{12}$ (black) and $RW_8$ (blue) as a function of $Ek$ with explicit (dashed) and implicit (solid) implementation of Coriolis force. For clarity, the insets show the full domain, although the computations are carried out in reduced domains $0 \le \varphi < 2\pi/M$. An average is taken over a branch of Ra values. The number of actions required by the explicit algorithm is always greater than that required by the implicit algorithm, with the ratio between them increasing from approximately 2.5 for $Ek = 3 \times 10^{-3}$ to 20 for $Ek = 3 \times 10^{-5}$, below which the explicit algorithm fails. The resolutions used for the $RW_8$ computations ranged from $(N_r, N_\theta, N_\varphi \times M) = (60, 112, 20 \times 12)$ to $(90, 184, 32 \times 12)$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIG. 12: Accuracy of the time-dependent calculation of rotating waves $RW_M$ as a function of timestep $\Delta t$ using implicit (circles) and explicit (crosses) timestepping of the Coriolis term. Shown here is the relative error $C_{error} = |(C_{\Delta t} - C_{exact})/C_{exact}|$ of the drift velocity as a function of the timestep for various sets of parameter values (see table) where $C_{exact}$ is calculated via Newton's method. Both the explicit and implicit methods are second order in time, with $C_{error} \sim \Delta t^2$. Vertical lines indicate the limiting $\Delta t$ above which explicit timestepping diverges and cannot be used. The horizontal dashed line indicates $C_{error} = 2\%$. The arrows designate the range of $\Delta t$ for which implicit timestepping is advantageous: it does not diverge and $C_{error}$ remains less than 2%. An approximate fit yields the range $0.5 \times Ek^{1.1} \lesssim \Delta t \lesssim 0.8 \times Ek^{0.9}$. See Table I.</figcaption></figure>

<table><caption>TABLE I: Parameters used for testing and useful range of $\Delta t$ for implicit timestepping.</caption><thead><tr><th rowspan="2">Ek</th><th rowspan="2">Ra</th><th rowspan="2">M</th><th rowspan="2">(N<sub>r</sub>, N<sub>θ</sub>, N<sub>φ</sub> × M)</th><th rowspan="2">C<sub>exact</sub></th><th colspan="3">Δt</th></tr><tr><th>Explicit</th><th>Diverges</th><th>2% Error</th><th>Ratio</th></tr></thead><tbody><tr><td>10<sup>-3</sup></td><td>120</td><td>4</td><td>(46, 72, 32 × 4)</td><td>-2.7647</td><td>4.0 × 10<sup>-4</sup></td><td>1.1 × 10<sup>-3</sup></td><td>3</td></tr><tr><td>10<sup>-4</sup></td><td>130</td><td>8</td><td>(60, 80, 40 × 4)</td><td>3.9998</td><td>2.0 × 10<sup>-5</sup></td><td>1.9 × 10<sup>-4</sup></td><td>9</td></tr><tr><td>10<sup>-5</sup></td><td>140</td><td>12</td><td>(76, 136, 24 × 4)</td><td>51.2817</td><td>1.6 × 10<sup>-6</sup></td><td>2.8 × 10<sup>-5</sup></td><td>17</td></tr></tbody></table>

We studied the transition to chaos for $Ek = 5 \times 10^{-6}$. Explicit Coriolis timestepping requires $\Delta t \lesssim 5 \times 10^{-7}$. Using the implicit method with $5 \Delta t$ and then $10 \Delta t$ allowed us to simulate five or 10 times longer and to tentatively confirm the temporally chaotic behavior we observed with the explicit method. If we wish to study this phenomenon further, we will decrease $\Delta t$ in order to obtain more accurate results.

As discussed in Section III B, treating the Coriolis term implicitly requires the use of block pentadiagonal matrices rather than block diagonal matrices. However, the additional cost is negligible. Part of the reason for this is that, in pseudospectral methods, the cost of carrying out the nonlinear terms (or the linearized nonlinear terms), more specifically the transforms to and from spectral space, is by far the most time-consuming portion of
the calculation. Thus, the cost of a timestep is barely increased by the implicit treatment
of the Coriolis term. In timestepping, if $\Delta t$ can be increased by a factor of 100 (assuming
the accuracy is still acceptable), then the simulation to the same final time costs (almost
exactly) 100 times less. In Newton's method, if 10 times fewer actions of the Jacobian (i.e.,
linearized timesteps) are required for a Newton step, then the computation will be (almost
exactly) 10 times faster.

# VI. DISCUSSION

## A. Differences with Other Methods

We now describe the differences between Stokes preconditioning and the method used
in [20, 45-47] for convection in rotating spheres, as well as in [74-78] for other hydrodynamic
problems. These authors seek steady states or traveling waves as roots of $U(t+T) - U(t)$,
where $U(t+T)$ is computed from $U(t)$ by carrying out many timesteps, each with a small $\Delta t$;
we therefore call this the integration method. In contrast, although Equations (19) and (20)
would seem to imply that we too seek roots of $U(t+T) - U(t)$, this is not at all the case.
In the Stokes preconditioning method, we compute $U(t+T)$ from $U(t)$ by carrying out a
single very large implicit-explicit Euler timestep Equation (17), with a timestep $\Delta t = T$
that is so large that the difference $U(t+\Delta t) - U(t)$ 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} + N$, preconditioned by $\mathcal{L}$.
Both methods solve the linear systems that are the core of Newton's method via GMRES,
which relies on the 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 of taking a single large linearized timestep. Here, there is a trade-off.
Despite the preconditioning displayed by Equation (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 our linear system Equation (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 a traveling 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 [48] 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 much faster (at least in the cases and
the parameter regimes that we have studied), the integration method has the considerable
advantage of being far 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, standing
waves, or pulsing states.

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

B. Relevance for Geophysics

One might wonder about the applicability of this achievement to geophysical research. We
begin by discussing the relevance of calculating traveling waves via Newton’s method. 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 which are partly or completely unstable, e.g.,
[64, 75–81]. This profusion of states, sometimes called exact coherent structures, has become
the focus of extensive research, motivated by the idea [82, 83] 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. Our Newton solver would enable the application of this line of research to geophysical
flows. More generally, unstable states can cast light on the origin and organization of
turbulent states; in one of the best known examples, ghostly Taylor vortices exist in even
highly turbulent Taylor–Couette flow [84–87] and the unstable underlying vortices may even
reproduce the mean properties of the turbulent flow [88].

We now turn to the relevance for geophysics of implicit Coriolis timestepping with large
timesteps. Constraints on timesteps (whether for stability or accuracy) arise from the differ-
ent physical forces in the equations. 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 (or the
anelastic approximation), which act to filter out the high-frequency sound waves [38, 39].

Using values of $\Delta t$ at which implicit Coriolis timestepping is possible but inaccurate can
also help to explore the large parameter space. The most economical procedure is to carry
out a preliminary large-scale survey with a coarse spatial resolution and a large timestep.
Interesting regions can then be accurately simulated using finer resolutions and smaller
timesteps.

Previous research by two of the authors [18, 19, 42, 43] originated in the context of the
Geoflow consortium, centered on microgravity experiments meant to mimic convective flows
within the Earth [29, 30]. These experiments were run in the Fluid Science Laboratory
facility located in the European Space Agency’s Columbus laboratory on the International
Space Station. The computations in [18, 19, 42, 43] used an explicit treatment of the Coriolis
force and could not be continued to Ekman numbers smaller than $10^{-4}$. Geoflow has been
succeeded by Atmoflow [31], whose purpose is to mimic convective flows in the atmosphere.
We hope to use the new code to carry out computations relevant to AtmoFlow and to be
able to achieve lower Ekman numbers.
VII. CONCLUSIONS

We have developed a pair of codes for simulating thermal convection in a rotating spherical
fluid shell that relies 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
the 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
(Δ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; the method developed here could be called Stokes-and-Coriolis
preconditioning.

We demonstrated this method’s capabilities by carrying out continuation in Rayleigh
number at various values of the Ekman number on the order of 10⁻⁵ for rotating waves with
azimuthal wavenumbers 4, 8, and 12. We found several intriguing examples of branches
containing 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 the
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 the
Ekman number or with automatic resolution adjustment in the literature.

We have measured the economy that is realized by the 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
RW₁₂ branch at Ek = 3 × 10⁻⁵. 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 Ek = 10⁻⁵, the implicit algorithm can rea-
sonably use timesteps that are almost 20 times past the temporal stability limit, meaning
that simulations in this regime can be 20 times faster than if an explicit method were used.
In conclusion, the implicit treatment of the Coriolis force greatly improves the efficiency of
computations in the low Ekman number regime.

**Code availability:** The code is in active development and fully available under the GNU General Public License v2.0 on the Github repository https://github.com/JuanSembla/spherical_code_pmmh.

**Funding:** The calculations for this work were performed using high performance computing resources provided by the Grand Equipement National de Calcul Intensif (GENCI) at the Institut du Développement et des Ressources en Informatique Scientifique (IDRIS, CNRS) through Grant No. A0162A01119.

**Acknowledgments:** We gratefully acknowledge Alan Riquier for programming assistance and we thank Adrian van Kan, Benjamin Favier, and Rainer Hollerbach for helpful discussions.
Appendix A

Here, we present the intermediate calculations leading from (11a) and (11b) to (15a) and (15b). Using (14), the recursion relations (13) for the associated Legendre polynomials can be rewritten in terms of the spherical harmonics as

$$ \sin \theta \frac{d}{d\theta} Y_l^{[m]} = \frac{l(l-|m|+1)}{2l+1} \frac{N_l^m}{N_{l+1}^m} Y_{l+1}^{[m]} - \frac{(l+|m|)(l+1)}{2l+1} \frac{N_l^m}{N_{l-1}^m} Y_{l-1}^{[m]} \quad (A1a) $$

$$ \cos \theta Y_l^{[m]} = \frac{l-|m|+1}{2l+1} \frac{N_l^m}{N_{l+1}^m} Y_{l+1}^{[m]} + \frac{(l+|m|)}{2l+1} \frac{N_l^m}{N_{l-1}^m} Y_{l-1}^{[m]} \quad (A1b) $$

Applying (A1) to (11a) leads to

$$ \begin{aligned} e_r \cdot \nabla \times \mathbf{F}^{\text{implicit}} = & \sum_{m,l} \left[ \frac{l(l+1)}{r^2} \text{Ek} \left( \frac{\partial}{\partial t} - L_l \right) - \frac{2im}{r^2} \right] e_{lm} Y_l^m \\ & + \frac{2}{r^2} \left( \frac{l(l+1)}{r} - \frac{\partial}{\partial r} \right) f_{lm} \left[ \frac{l(l-|m|+1)}{2l+1} \frac{N_l^m}{N_{l+1}^m} Y_{l+1}^{[m]} - \frac{(l+|m|)(l+1)}{2l+1} \frac{N_l^m}{N_{l-1}^m} Y_{l-1}^{[m]} \right] \\ & + 2 \frac{l(l+1)}{r^2} \left( \frac{2}{r} - \frac{\partial}{\partial r} \right) f_{lm} \left[ \frac{l-|m|+1}{2l+1} \frac{N_l^m}{N_{l+1}^m} Y_{l+1}^{[m]} + \frac{(l+|m|)}{2l+1} \frac{N_l^m}{N_{l-1}^m} Y_{l-1}^{[m]} \right] \end{aligned} $$

Grouping terms in $Y_{l+1}^{[m]}$ and $Y_{l-1}^{[m]}$, we obtain

$$ \begin{aligned} & = \sum_{m,l} \left[ \frac{l(l+1)}{r^2} \text{Ek} \left( \frac{\partial}{\partial t} - L_l \right) - \frac{2im}{r^2} \right] e_{lm} Y_l^m \\ & \quad + \frac{2}{r^2} \frac{l(l-|m|+1)}{2l+1} \left[ \left( \frac{l(l+1)}{r} - \frac{\partial}{\partial r} \right) + (l+1) \left( \frac{2}{r} - \frac{\partial}{\partial r} \right) \right] \frac{N_l^m}{N_{l+1}^m} f_{lm} Y_{l+1}^m \\ & \quad + \frac{2}{r^2} \frac{(l+|m|)(l+1)}{2l+1} \left[ -\left( \frac{l(l+1)}{r} - \frac{\partial}{\partial r} \right) + l \left( \frac{2}{r} - \frac{\partial}{\partial r} \right) \right] \frac{N_l^m}{N_{l-1}^m} f_{lm} Y_{l-1}^m \end{aligned} $$

$$ \begin{aligned} & = \sum_{m,l} \left[ \frac{l(l+1)}{r^2} \text{Ek} \left( \frac{\partial}{\partial t} - L_l \right) - \frac{2im}{r^2} \right] e_{lm} Y_l^m \\ & \quad + \frac{2}{r^2} \frac{(l-|m|+1)l(l+2)}{2l+1} \left[ \left( \frac{l+1}{r} - \frac{\partial}{\partial r} \right) \frac{N_l^m}{N_{l+1}^m} f_{lm} Y_{l+1}^m \right. \\ & \quad \left. + \frac{2}{r^2} \frac{(l+|m|)(l+1)(l-1)}{2l+1} \left[ -\frac{l}{r} - \frac{\partial}{\partial r} \right] \frac{N_l^m}{N_{l-1}^m} f_{lm} Y_{l-1}^m \right] \end{aligned} $$
Finally, changing the names of the indices so that all terms correspond to the same spherical harmonic $Y_l^m$, we obtain Equation (15a):

$$ = \sum_{m,l} Y_l^m \left\{ \left[ \frac{l(l+1)}{r^2} \text{Ek} \left( \frac{\partial}{\partial t} - L_l \right) - \frac{2im}{r^2} \right] e_{lm} \right. \\ \left. + \frac{2}{r^2} \frac{(l-1)(l+1)(l-|m|)}{2l-1} \left[ \frac{l}{r} - \frac{\partial}{\partial r} \right] \frac{N_{l-1}^m}{N_l^m} f_{l-1,m} \right. \\ \left. + \frac{2}{r^2} \frac{(l+1+|m|)l(l+2)}{2l+3} \left[ -\frac{(l+1)}{r} - \frac{\partial}{\partial r} \right] \frac{N_{l+1}^m}{N_l^m} f_{l+1,m} \right\} \quad (A2) $$

The same calculations can be applied to (11b) to obtain (15b):

$$ \begin{aligned} e_r \cdot \nabla \times \nabla \times F^{\text{implicit}} &= \sum_{m,l} - \left[ \frac{l(l+1)}{r^2} \text{Ek} \left( \frac{\partial}{\partial t} - L_l \right) - \frac{2im}{r^2} \right] L_l f_{lm} Y_l^m \\ &\quad + \frac{2}{r^2} \left( \frac{l(l+1)}{r} - \frac{\partial}{\partial r} \right) e_{lm} \left[ \frac{l(l-|m|+1)}{2l+1} \frac{N_l^m}{N_{l+1}^m} Y_{l+1}^{|m|} - \frac{(l+|m|)(l+1)}{2l+1} \frac{N_l^m}{N_{l-1}^m} Y_{l-1}^{|m|} \right] \\ &\quad + 2 \frac{l(l+1)}{r^2} \left( \frac{2}{r} - \frac{\partial}{\partial r} \right) e_{lm} \left[ \frac{l-|m|+1}{2l+1} \frac{N_l^m}{N_{l+1}^m} Y_{l+1}^{|m|} + \frac{(l+|m|)}{2l+1} \frac{N_l^m}{N_{l-1}^m} Y_{l-1}^{|m|} \right] \end{aligned} $$

$$ = \sum_{m,l} - \left[ \frac{l(l+1)}{r^2} \text{Ek} \left( \frac{\partial}{\partial t} - L_l \right) - \frac{2im}{r^2} \right] L_l f_{lm} Y_l^m \\ + \frac{2}{r^2} \frac{l(l-|m|+1)}{2l+1} \left[ \left( \frac{l(l+1)}{r} - \frac{\partial}{\partial r} \right) + (l+1) \left( \frac{2}{r} - \frac{\partial}{\partial r} \right) \right] \frac{N_l^m}{N_{l+1}^m} e_{lm} Y_{l+1}^m \\ + \frac{2}{r^2} \frac{(l+|m|)(l+1)}{2l+1} \left[ - \left( \frac{l(l+1)}{r} - \frac{\partial}{\partial r} \right) + l \left( \frac{2}{r} - \frac{\partial}{\partial r} \right) \right] \frac{N_l^m}{N_{l-1}^m} e_{lm} Y_{l-1}^m $$

$$ = \sum_{m,l} - \left[ \frac{l(l+1)}{r^2} \text{Ek} \left( \frac{\partial}{\partial t} - L_l \right) - \frac{2im}{r^2} \right] L_l f_{lm} Y_l^m \\ + \frac{2}{r^2} \frac{l(l-|m|+1)(l+2)}{2l+1} \left[ \frac{(l+1)}{r} - \frac{\partial}{\partial r} \right] \frac{N_l^m}{N_{l+1}^m} e_{lm} Y_{l+1}^m \\ + \frac{2}{r^2} \frac{(l+|m|)(l+1)(l-1)}{2l+1} \left[ -\frac{l}{r} - \frac{\partial}{\partial r} \right] \frac{N_l^m}{N_{l-1}^m} e_{lm} Y_{l-1}^m $$

$$ = \sum_{m,l} Y_l^m \left\{ - \left[ \frac{l(l+1)}{r^2} \text{Ek} \left( \frac{\partial}{\partial t} - L_l \right) - \frac{2im}{r^2} \right] L_l f_{lm} \right. \\ \left. + \frac{2}{r^2} \frac{(l-1)(l-|m|)(l+1)}{2l-1} \left[ \frac{l}{r} - \frac{\partial}{\partial r} \right] \frac{N_{l-1}^m}{N_l^m} e_{l-1,m} \right. \\ \left. + \frac{2}{r^2} \frac{(l+1+|m|)(l+2)l}{2l+3} \left[ -\frac{(l+1)}{r} - \frac{\partial}{\partial r} \right] \frac{N_{l+1}^m}{N_l^m} e_{l+1,m} \right\} \quad (A3) $$
[1] Chandrasekhar, S. *Hydrodynamic and Hydromagnetic Stability*; Oxford University Press: Oxford, UK, 1961.

[2] Roberts, P.H. On the thermal instability of a rotating-fluid sphere containing heat sources. *Phil. Trans. R. Soc. Lond. A* **1968**, *263*, 93-117.

[3] Busse, F. Thermal instabilities in rapidly rotating systems. *J. Fluid Mech.* **1970**, *44*, 441-460.

[4] Zhang, K.K.; Busse, F.H. On the onset of convection in rotating spherical shells. *Geophys. Astrophys. Fluid Dyn.* **1987**, *39*, 119-147.

[5] Jones, C.A.; Soward, A.M.; Mussa, A.I. The onset of thermal convection in a rapidly rotating sphere. *J. Fluid Mech.* **2000**, *405*, 157-179.

[6] Dormy, E.; Soward, A.; Jones, C.; Jault, D.; Cardin, P. The onset of thermal convection in rotating spherical shells. *J. Fluid Mech.* **2004**, *501*, 43-70.

[7] Net, M.; Garcia, F.; Sánchez, J. On the onset of low-Prandtl-number convection in rotating spherical shells: Non-slip boundary conditions. *J. Fluid Mech.* **2008**, *601*, 317-337.

[8] Skene, C.S.; Tobias, S.M. Weakly nonlinear analysis of the onset of convection in rotating spherical shells. *Geophys. Astrophys. Fluid Dyn.* **2024**, 1-22.

[9] Ruelle, D. Bifurcations in the presence of a symmetry group. *Arch. Rat. Mech. Anal.* **1973**, *51*, 136-152.

[10] Ecke, R.; Zhong, F.; Knobloch, E. Hopf bifurcation with broken reflection symmetry in rotating Rayleigh-Bénard convection. *Europhys. Lett.* **1992**, *19*, 177.

[11] Chossat, P.; Lauterbach, R. *Methods in Equivariant Bifurcations and Dynamical Systems*; World Scientific: Singapore, 2000.

[12] Busse, F. Convective flows in rapidly rotating spheres and their dynamo action. *Phys. Fluids* **2002**, *14*, 1301-1314.

[13] Hindman, B.W.; Jain, R. Radial trapping of thermal Rossby waves within the convection zones of low-mass stars. *Astrophys. J.* **2022**, *932*, 68.

[14] Dormy, E. Rapidly Rotating Magnetohydrodynamics and the Geodynamo. *Annu. Rev. Fluid Mech.* **2025**, *57*.

[15] Pothérat, A.; Horn, S. Seven decades of exploring planetary interiors with rotating convection experiments. *C. R. Phys.* **2025**, *25*, 233.

[16] Zhang, K. Convection in a rapidly rotating spherical shell at infinite Prandtl number: Transition to vacillating flows. *Phys. Earth Planet. Inter.* **1992**, *72*, 236-248.

[17] Simitev, R.; Busse, F. Patterns of convection in rotating spherical shells. *New J. Phys.* **2003**, *5*, 97. https://doi.org/10.1088/1367-2630/5/1/397.

[18] Feudel, F.; Seehafer, N.; Tuckerman, L.S.; Gellert, M. Multistability in rotating spherical shell convection. *Phys. Rev. E* **2013**, *87*, 023021.

[19] Feudel, F.; Tuckerman, L.S.; Gellert, M.; Seehafer, N. Bifurcations of rotating waves in rotating spherical shell convection. *Phys. Rev. E* **2015**, *92*, 053015.

[20] Garcia, F.; Net, M.; Sánchez, J. Continuation and stability of convective modulated rotating waves in spherical shells. *Phys. Rev. E* **2016**, *93*, 013119.

[21] Cardin, P.; Olson, P. Chaotic thermal convection in a rapidly rotating spherical shell: Consequences for flow in the outer core. *Phys. Earth Planet. Inter.* **1994**, *82*, 235-259.

[22] Aubert, J.; Brito, D.; Nataf, H.C.; Cardin, P.; Masson, J.P. A systematic experimental study of rapidly rotating spherical convection in water and liquid gallium. *Phys. Earth Planet. Inter.*
2001, 128, 51–74.

[23] Jones, C.A. Thermal and compositional convection in the outer core. *Core Dyn.* **2007**, *8*, 131–185.

[24] Aurnou, J.M.; Horn, S.; Julien, K. Connections between nonrotating, slowly rotating, and rapidly rotating turbulent convection transport scalings. *Phys. Rev. Res.* **2020**, *2*, 043115.

[25] Olson, P. Laboratory experiments on the dynamics of the core. *Phys. Earth Planet. Inter.* **2011**, *187*, 139–156.

[26] Glatzmaier, G.A. Numerical simulations of stellar convective dynamos. I. the model and method. *J. Comput. Phys.* **1984**, *55*, 461–484. https://doi.org/https://doi.org/10.1016/0021-9991(84)90033-0.

[27] Christensen, U.; Aubert, J.; Cardin, P.; Dormy, E.; Gibbons, S.; Glatzmaier, G.; Grote, E.; Honkura, Y.; Jones, C.; Kono, M.; et al. A numerical dynamo benchmark. *Phys. Earth Planet. Inter.* **2001**, *128*, 25–34.

[28] Jones, C.; Boronski, P.; Brun, A.; Glatzmaier, G.; Gastine, T.; Miesch, M.; Wicht, J. Anelastic convection-driven dynamo benchmarks. *Icarus* **2011**, *216*, 120–135.

[29] Futterer, B.; Krebs, A.; Plesa, A.C.; Zaussinger, F.; Hollerbach, R.; Breuer, D.; Egbers, C. Sheet-like and plume-like thermal flow in a spherical convection experiment performed under microgravity. *J. Fluid Mech.* **2013**, *735*, 647–683.

[30] Zaussinger, F.; Haun, P.; Szabo, P.S.B.; Travnikov, V.; Al Kawwas, M.; Egbers, C. Rotating spherical gap convection in the GeoFlow International Space Station (ISS) experiment. *Phys. Rev. Fluids* **2020**, *5*, 063502. https://doi.org/10.1103/PhysRevFluids.5.063502.

[31] Travnikov, V.; Egbers, C. Numerical investigation of atmospherelike flows in a spherical geometry. *Phys. Rev. E* **2021**, *104*, 065110. https://doi.org/10.1103/PhysRevE.104.065110.

[32] Tuckerman, L.S. Transformations of matrices into banded form. *J. Comput. Phys.* **1989**, *84*, 360–376.

[33] Marti, P.; Calkins, M.; Julien, K. A computationally efficient spectral method for modeling core dynamics. *Geochem. Geophys. Geosystems* **2016**, *17*, 3031–3053.

[34] Kuang, W.; Bloxham, J. Numerical modeling of magnetohydrodynamic convection in a rapidly rotating spherical shell: Weak and strong field dynamo action. *J. Comput. Phys.* **1999**, *153*, 51–81.

[35] Hollerbach, R. A spectral solution of the magneto-convection equations in spherical geometry. *Int. J. Numer. Methods Fluids* **2000**, *32*, 773–797.

[36] Sprague, M.; Julien, K.; Knobloch, E.; Werne, J. Numerical simulation of an asymptotically reduced system for rotationally constrained convection. *J. Fluid Mech.* **2006**, *551*, 141–174.

[37] Stellmach, S.; Hansen, U. An efficient spectral method for the simulation of dynamos in Cartesian geometry and its implementation on massively parallel computers. *Geochem. Geophys. Geosyst.* **2008**, *9*, https://doi.org/10.1029/2007GC001778.

[38] Julien, K.; van Kan, A.; Miquel, B.; Knobloch, E.; Vasil, G. Rescaled Equations of Fluid Motion for Well-Conditioned Direct Numerical Simulations of Rapidly Rotating Convection. *arXiv* **2024**, arXiv:2410.02702.

[39] van Kan, A.; Julien, K.; Miquel, B.; Knobloch, E. Bridging the Rossby number gap in rapidly rotating thermal convection. *J. Fluid Mech.* **2025**, *1010*, A42.

[40] García González, F.; Net Marcé, M.; Garcia-Archilla, B.; Sánchez Umbría, J. A comparison of high-order time integrators for thermal convection in rotating spherical shells. *J. Comput. Phys.* **2010**, *229*, 7997–8010.
[41] Garcia, F.; Sánchez, J.; Dormy, E.; Net, M. Oscillatory Convection in Rotating Spherical Shells: Low Prandtl Number and Non-Slip Boundary Conditions. *SIAM J. Appl. Dyn. Syst.* **2015**, *14*, 1787–1807.

[42] Feudel, F.; Bergemann, K.; Tuckerman, L.S.; Egbers, C.; Futterer, B.; Gellert, M.; Hollerbach, R. Convection patterns in a spherical fluid shell. *Phys. Rev. E* **2011**, *83*, 046304.

[43] Feudel, F.; Tuckerman, L.S.; Zaks, M.; Hollerbach, R. Hysteresis of dynamos in rotating spherical shell convection. *Phys. Rev. F* **2017**, *2*, 053902.

[44] Schaeffer, N. Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations. *Geochem. Geophys. Geosyst.* **2013**, *14*, 751–758.

[45] Sánchez, J.; Net, M.; Garcia-Archilla, B.; Simó, C. Newton–Krylov continuation of periodic orbits for Navier–Stokes flows. *J. Comput. Phys.* **2004**, *201*, 13–33.

[46] Garcia, F.; Sánchez, J.; Net, M. Antisymmetric Polar Modes of Thermal Convection in Rotating Spherical Fluid Shells at High Taylor Numbers. *Phys. Rev. Lett.* **2008**, *101*, 194501.

[47] Sánchez, J.; Garcia, F.; Net, M. Computation of azimuthal waves and their stability in thermal convection in rotating spherical shells with application to the study of a double-Hopf bifurcation. *Phys. Rev. E* **2013**, *87*, 033014.

[48] Tuckerman, L.S.; Langham, J.; Willis, A. Order-of-magnitude speedup for steady states and traveling waves via Stokes preconditioning in Channelflow and Openpipeflow. In *Computational Modelling of Bifurcations and Instabilities in Fluid Dynamics*; Springer: Berlin/Heidelberg, Germany, 2018; pp. 3–31.

[49] Gastine, T.; Wicht, J.; Barik, A.; Duarte, L. Formulation of the (magneto)-hydrodynamics problem. Available at https://magic-sph.github.io/ (accessed on 30 May 2025).

[50] Marques, F.; Mercader, I.; Batiste, O.; Lopez, J.M. Centrifugal effects in rotating convection: Axisymmetric states and three-dimensional instabilities. *J. Fluid Mech.* **2007**, *580*, 303–318.

[51] Satake, H.; Tagawa, T. Influence of Centrifugal Buoyancy in Thermal Convection within a Rotating Spherical Shell. *Symmetry* **2022**, *14*, 2021. https://doi.org/10.3390/sym14102021.

[52] Ishioka, K. A New Recurrence Formula for Efficient Computation of Spherical Harmonic Transform. *J. Meteorol. Soc. Jpn.* **2018**, *96*, 241–249. https://doi.org/10.2151/jmsj.2018-019.

[53] Schaeffer, N.; Jault, D.; Nataf, H.C.; Fournier, A. Turbulent geodynamo simulations: A leap towards Earth’s core. *Geophys. J. Int.* **2017**, *211*, 1–29.

[54] Lago, R.; Gastine, T.; Dannert, T.; Rampp, M.; Wicht, J. MagIC v5. 10: A two-dimensional message-passing interface (MPI) distribution for pseudo-spectral magnetohydrodynamics simulations in spherical geometry. *Geosci. Model Dev.* **2021**, *14*, 7477–7495.

[55] Garcia, F.; Sánchez, J.; Net, M. Numerical simulations of thermal convection in rotating spherical shells under laboratory conditions. *Phys. Earth Planet. Inter.* **2014**, *230*, 28–44. https://doi.org/https://doi.org/10.1016/j.pepi.2014.02.004.

[56] Gottlieb, D.; Orszag, S.A. *Numerical Analysis of Spectral Methods: Theory and Applications*; SIAM: Philadelphia, US 1977.

[57] Saad, Y.; Schultz, M.H. GMRES: A generalized minimal residual algorithm for solving non-symmetric linear systems. *SIAM J. Sci. Stat. Comput.* **1986**, *7*, 856–869.

[58] Tuckerman, L.S. Steady-state solving via Stokes preconditioning; recursion relations for elliptic operators. In Proceedings of the 11th International Conference on Numerical Methods in Fluid Dynamics, Williamsburg, VA, USA, 27 June 1988; Springer: Berlin/Heidelberg, Germany, 1989; pp. 573–577.

[59] Mamun, C.K.; Tuckerman, L.S. Asymmetry and Hopf bifurcation in spherical Couette flow.
Phys. Fluids **1995**, *7*, 80–91.

[60] Bergeon, A.; Henry, D.; Benhadid, H.; Tuckerman, L.S. Marangoni convection in binary mixtures with Soret effect. *J. Fluid Mech.* **1998**, *375*, 143–177.

[61] Nore, C.; Tuckerman, L.S.; Daube, O.; Xin, S. The 1:2 mode interaction in exactly counter-rotating von Kármán swirling flow. *J. Fluid Mech.* **2003**, *477*, 51–88.

[62] Meca, E.; Mercader, I.; Batiste, O.; Ramírez-Piscina, L. Blue sky catastrophe in double-diffusive convection. *Phys. Rev. Lett.* **2004**, *92*, 234501.

[63] Batiste, O.; Knobloch, E.; Alonso, A.; Mercader, I. Spatially localized binary-fluid convection. *J. Fluid Mech.* **2006**, *560*, 149–158.

[64] Borońska, K.; Tuckerman, L.S. Extreme multiplicity in cylindrical Rayleigh-Bénard convection. II. Bifurcation diagram and symmetry classification. *Phys. Rev. E* **2010**, *81*, 036321.

[65] Mercader, I.; Batiste, O.; Alonso, A.; Knobloch, E. Convectons, anticonvectons and multi-convectons in binary fluid convection. *J. Fluid Mech.* **2011**, *667*, 586–606.

[66] Beaume, C.; Bergeon, A.; Kao, H.C.; Knobloch, E. Convectons in a rotating fluid layer. *J. Fluid Mech.* **2013**, *717*, 417–448.

[67] Torres, J.F.; Henry, D.; Komiya, A.; Maruyama, S. Bifurcation analysis of steady natural convection in a tilted cubical cavity with adiabatic sidewalls. *J. Fluid Mech.* **2014**, *756*, 650–688.

[68] Huepe, C.; Metens, S.; Dewel, G.; Borckmans, P.; Brachet, M.E. Decay rates in attractive Bose-Einstein condensates. *Phys. Rev. Lett.* **1999**, *82*, 1616.

[69] Huepe, C.; Tuckerman, L.S.; Métens, S.; Brachet, M.E. Stability and decay rates of non-isotropic attractive Bose-Einstein condensates. *Phys. Rev. A* **2003**, *68*, 023609.

[70] Tuckerman, L.S. Laplacian preconditioning for the inverse Arnoldi method. *Commun. Comput. Phys.* **2015**, *18*, 1336–1351.

[71] Lanczos, C. Trigonometric interpolation of empirical and analytical functions. *J. Math. Phys.* **1938**, *17*, 123–199.

[72] Trefethen, L.N. *Spectral Methods in MATLAB*; SIAM, Philadelphia, US 2000.

[73] Marcus, P.S. Simulation of Taylor-Couette flow. Part 1. Numerical methods and comparison with experiment. *J. Fluid Mech.* **1984**, *146*, 45–64.

[74] Duguet, Y.; Pringle, C.C.; Kerswell, R.R. Relative periodic orbits in transitional pipe flow. *Phys. Fluids* **2008**, *20*, 114102.

[75] Gibson, J.F.; Halcrow, J.; Cvitanović, P. Equilibrium and travelling-wave solutions of plane Couette flow. *J. Fluid Mech.* **2009**, *638*, 243–266.

[76] Zheng, Z.; Tuckerman, L.S.; Schneider, T.M. Natural convection in a vertical channel. Part 1. Wavenumber interaction and Eckhaus instability in a narrow domain. *J. Fluid Mech.* **2024**, *1000*, A28.

[77] Zheng, Z.; Tuckerman, L.S.; Schneider, T.M. Natural convection in a vertical channel. Part 2. Oblique solutions and global bifurcations in a spanwise-extended domain. *J. Fluid Mech.* **2024**, *1000*, A29.

[78] Zheng, Z.; Tuckerman, L.S.; Schneider, T.M. Natural convection in a vertical channel. Part 3. Bifurcations of many (additional) unstable equilibria and periodic orbits. *arXiv* **2025**, arXiv:2504.05524.

[79] Schmiegel, A. Transition to Turbulence in Linearly Stable Shear Flows. PhD Thesis, Philipps-Universität Marburg, Marburg, Germany, 1999. Available online: http://archiv.ub.uni-marburg.de/diss/z2000/0062/ (accessed on 30 May 2025).

[80] Pringle, C.C.; Duguet, Y.; Kerswell, R.R. Highly symmetric travelling waves in pipe flow.
*Phil. Trans. R. Soc. Lond. A 2009*, **367**, 457–472.

[81] Borońska, K.; Tuckerman, L.S. Extreme multiplicity in cylindrical Rayleigh-Bénard convection. I. Time dependence and oscillations. *Phys. Rev. E* **2010**, *81*, 036320.

[82] Cvitanovic, P.; Eckhardt, B. Periodic orbit expansions for classical smooth flows. *J. Phys. A-Math. 1991*, **24**, L237.

[83] Kawahara, G.; Uhlmann, M.; Van Veen, L. The significance of simple invariant solutions in turbulent flows. *Annu. Rev. Fluid Mech. 2012*, **44**, 203–225.

[84] Lathrop, D.P.; Fineberg, J.; Swinney, H.L. Transition to shear-driven turbulence in Couette-Taylor flow. *Phys. Rev. A* **1992**, *46*, 6390.

[85] Dong, S. Direct numerical simulation of turbulent Taylor–Couette flow. *J. Fluid Mech. 2007*, **587**, 373–393.

[86] Huisman, S.G.; Van Der Veen, R.C.; Sun, C.; Lohse, D. Multiple states in highly turbulent Taylor–Couette flow. *Nat. Comm. 2014*, **5**, 3820.

[87] Grossmann, S.; Lohse, D.; Sun, C. High-Reynolds number Taylor-Couette turbulence. *Annu. Rev. Fluid Mech. 2016*, **48**, 53–80.

[88] Eckhardt, B.; Doering, C.R.; Whitehead, J.P. Exact relations between Rayleigh-Bénard and rotating plane Couette flow in two dimensions. *J. Fluid Mech. 2020*, **903**, R4.