Tensor calculus in spherical coordinates using Jacobi polynomials, Part-II: Implementation and Examples

We present a simulation code which can solve broad ranges of partial differential equations in a full sphere. The code expands tensorial variables in a spectral series of spin-weighted spherical harmonics in the angular directions and a scaled Jacobi polynomial basis in the radial direction, as described in Part-I. Nonlinear terms are calculated by transforming from the coefficients in the spectral series to the value of each quantity on the physical grid, where it is easy to calculate products and perform other local operations. The expansion makes it straightforward to solve equations in tensor form (i.e., without decomposition into scalars). We propose and study several unit tests which demonstrate the code can accurately solve linear problems, implement boundary conditions, and transform between spectral and physical space. We then run a series of benchmark problems proposed in Marti et al (2014), implementing the hydrodynamic and magnetohydrodynamic equations. We are able to calculate more accurate solutions than reported in Marti et al 2014 by running at higher spatial resolution and using a higher-order timestepping scheme. We find the rotating convection and convective dynamo benchmark problems depend sensitively on details of timestepping and data analysis. We also demonstrate that in low resolution simulations of the dynamo problem, small changes in a numerical scheme can lead to large changes in the solution. To aid future comparison to these benchmarks, we include the source code used to generate the data, as well as the data and analysis scripts used to generate the figures.


Introduction
Stars and planets are spherical to an excellent approximation. This makes spherical coordinates a natural choice for solving problems in astrophysical and geophysical fluid dynamics (e.g., [15,6], but see [10] for an alternative view). The spherical coordinates (r, θ, φ) have two types of coordinates singularities: at θ = 0, π; and at r = 0. Functions written in spherical coordinates must satisfy regularity conditions near these coordinate singularities [e.g., 11].
In [19, hereafter, Part-I], we discuss a strategy for computing general tensor-calculus operations on functions in the three-dimensional ball. This naturally leads to methods for solving a wide class of partial differential equations (PDEs) in spherical coordinates. We expand each of the PDEs' dependent variables in a spectral series using spin-weighted spherical harmonics for the θ and φ dependence [e.g., 16] and a scaled class of Jacobi polynomials for the r dependence [similar to 14]. Each basis function satisfies the regularity conditions at the coordinate singularities, so their sum automatically does as well. This is a similar to our approach to simulations in cylindrical geometry [18]. The results from the disk provide an introduction to the more complex geometry of the full three-dimensional ball.
Previous researchers [e.g., 11,14] have derived similar radial basis functions for scalar variables. Part-I provides a more thorough overview of the numerous different methods that have been developed to accurately cope with the large dynamic range associated with polar coordinate singularities. [3] also provide an excellent introduction to the topic in general. Although vectors and higher order tensors can be decomposed into their scalar components (e.g., toroidal-poloidal decomposition for divergence-free vectors), this becomes tedious for high rank tensors. In contrast, Part-I derives basis functions for arbitrary tensorial variables. Tensors of different ranks are linked by sparse derivative operators, expressing various tensorial relations. For example, the gradient of a vector is a rank-2 tensor, the divergence of a vector is a scalar, etc. This makes it possible to solve tensorial equations in primitive form (e.g., without decomposition into scalars), making this method applicable to wide classes of PDEs. [9] derived a similar basis for tensorial quantities, but mostly focused on the cartesian components (i.e., x, y, z components of a tensor). This paper contains a series of tests which demonstrate the utility of our method. We run unit tests: eigenvalue problems (section 4) and boundary value problems (section 5). These tests demonstrate that we can accurately solve linear problems and perform transforms from physical space to spectral space. They also test our implementation of a wide variety of boundary conditions that are used in hydrodynamics and magnetohydrodynamics.
We also run full-code tests. In section 6 we simulate all three of the full-sphere benchmark problems described in [13, hereafter, M14]. Section 6 of this paper can be thought of as a follow-up to M14 as we calculate converged solutions to higher precision by running at higher resolution and with higher-order time steppers. We include details of our simulations and data analysis that are necessary to compare between codes (e.g., timestepping scheme, timestep size, etc.). The supplementary materials also include the full source code so the interested reader could confirm any details of the simulations. Reduced data outputs and analysis scripts are public and in the repository https://github.com/lecoanet/dedalus_ sphere and at https://princeton.edu/~lecoanet/data.

Summary of the Algorithm
In this paper we solve initial value problems, boundary value problems, and eigenvalue problems using the algorithms derived in Part-I. As an example, consider the initial value problem, where X is a state vector consisting of a list of tensorial fields. Examples from fluid dynamics include: scalar fields (rank 0), e.g., density, temperature, pressure, the divergence of the velocity; vector fields (rank 1), e.g., velocities, magnetic fields, temperature gradient; rank 2 tensor fields, e.g., the strain rate, Maxwell stress; and higher order tensors. M and L are linear operators, possibly including derivative operators such as gradients, divergence, curl, etc. F is a general nonlinear function.
We solve equation (1) in spherical polar coordinates (r, θ, φ), with r ∈ [0, 1], θ ∈ [0, π], and φ ∈ [0, 2π), subject to the boundary conditions where B is a linear operator and E is a nonlinear function, and initial conditions There are analogous formulations of boundary value and eigenvalues problems that we discuss below.
Our approach is to expand X in the bases described in Part-I, and then rewrite the problem in terms of the coefficients of the basis elements. Here we briefly summarize some important results of Part-I. Consider a rank-r tensor T. Then T has 3 r components, corresponding to a linear combination of tensor products of the coordinate unit vectors e r , e θ , and e φ .
The element e(i) represents a single tensorial component using multi-index notation. For example for rank-3 tensors, i = {0, 0, 1} corresponds to e(i) = e r ⊗ e r ⊗ e θ . See the appendix in Part-I for a discussion of the multi-index notation. Therefore, T(r, θ, φ) = i,σ,a m, nT a m, ,n Q α, +ā n (r) Q (σ, a) Yσ ,m (θ, φ) U (σ, i) e(i), where the Q α,l+a n (r) are related to a set of Jacobi polynomials, Q l is an -dependent orthogonal transformation, Yσ ,m (θ, φ) are the spin-weighted spherical harmonics, U (σ, i) is a unitary transformation, and overbars denote a sum over multi-indices. Thus,T a m, ,n is a (complexvalued) coefficient of T using this basis. Our choice of basis ensures the solutions satisfy regularity conditions at the poles (θ = 0 and π) and the origin (r = 0), and these bases ensure that derivative operators are maximally sparse (Part-I).
For calculations, we must truncate the sums in equation (4). We pick a value of L max and N max . The maximum spherical harmonic degree is L max , and the azimuthal order m ranges 0 ≤ m ≤ L max . Note that we do not need negative values of m because we assume the tensor field is real. The truncation for modes depends on the spin of the component in question, σ. For each azimuthal order m, we have max(m, |σ|) ≤ ≤ L max . We use· to denote the sum of the elements of a given spin or regularity multi-index, i.e., if σ = {+1, −1}, then σ = (+1) + (−1) = 0 (see Part-I for more details). The truncation for tensors is similar to the familiar triangular truncation for scalar spherical harmonics. The truncation requires that the degree of sin(θ) in the spin-weighted spherical harmonic is no greater than L max .
The Q polynomials are where P n is a Jacobi polynomial of degree n and the proportionality is determined by a normalization factor. Thus, this Q polynomial has an r degree of +ā + 2n. Similar to spherical harmonics, we also impose a triangular truncation, requiring that the radial degree of each Q polynomial is bounded. For each problem, we determine the highest tensor rank we are interested in, R max , soā ≤ R max . We then require that 0 ≤ 2n ≤ 2N max − + R max . The range of n's depends on the value of . Our use of R max means there are the same number of n for tensor components with different a, which greatly simplifies our analysis.
The sum over i, σ, and τ are, respectively, sums over the spherical components of the tensor T, the spin indices, and the regularity indices. Each multi-index of a rank-r tensor has 3 r elements. There is also one additional index, α. All fields start off with α = 0, but operations like differentiation increase the value of α by one. We use different values of α to keep the differentiation matrices maximally sparse. This is equivalent to the sparse derivative relation between Chebyshev-T and Chebyshev-U polynomials. We use conversion matrices to ensure all variables in an equation have the same value of α.
After the truncation of equation (4), the state vector consists of O(L 2 max N max ) expansion coefficients for each tensor component. In this paper we only study problems in which the linear operators (e.g., M , L, B) do not contain any explicit dependence on θ or φ (but we allow coupling in all directions through gradient operators). In this case, the linear operators only couple different radial modes together. Thus, we can consider equation (1) as O(L 2 max ) different equations for X m, , the state vector corresponding to spherical harmonic order m and . Each tensor component has O(N max ) components in X m, . The linear operators acting on X m, are sparse with O(N max ) elements. They can be easily inverted with off-the-shelf sparse linear algebra packages. Coupling in only the radial direction and between different field variables allows the parallelization of the code across both m and . With this restriction on linear operators, we cannot treat terms like the Coriolis force with implicit time stepping (though we can treat it as a part of the nonlinear operator F (X)).
To calculate the nonlinear terms F (X) and E(X)| r=1 , we transform the coefficientsT a m, ,n into the tensor field T(r, θ, φ) in physical space, and then perform any local operations (e.g., products) in physical space. The transform requires O(L 2 max ) matrix-multiply transforms for the radial basis (i.e., multiplication by a dense, O(N 2 max ) matrix), and O(L max N max ) matrix-multiply transforms for the angular basis (i.e., multiplication by a dense, O(L 2 max ) matrix). Thus, the transformations require O(N max L 2 max max(N max , L max )) operations and are expected to be the slowest part of the calculation when L max and N max become large. Practically speaking, these transformations are reliant on the speed of the linear algebra library, in particular the speed of matrix and vector dot products; these are typically welloptimized numerical operations.
In section 3, we describe how we transform between data in physical space (T(r, θ, φ)), and the coefficient expansion of equation (4) (T a m, ,n ). In the subsequent sections, we describe the implementation of this formulation for eigenvalue problems (section 4), a boundary value problem (section 5), and the three initial value problems described in M14 (section 6).

Transforms
Here we describe how we transform data from physical space (T(φ, θ, r)), back and forth from the coefficient expansion (T a m, ,n ) in terms of scaled Jacobi polynomials in the radial direction and spin-weighted spherical harmonics in the angular directions. This is crucial for efficiently calculating nonlinear products, as well as visualizing our data.
Consider a rank-r tensor T(φ, θ, r). We represent T with N c = 3 r components, each of which have N r radial points, N θ latitudinal points and N φ longitudinal points. The data are initially on the quadrature nodes of the spin-weighted spherical harmonics in φ, θ and the Jacobi polynomials in r. We assume the components of T are real. We represent each component of T as a Field in Dedalus 1 , as it wraps FFTW's Fourier transforms and parallel transposes. The data for T(φ, θ, r) are stored as a N c × N φ × N θ × N r numpy array.
The grid points and transform matrices are calculated using Gaussian quadrature. We generate a guess for the quadrature grid and weights using the Golub-Welsch algorithm. After this we polish the results using a Newton iteration. We use Legendre quadrature for the latitudinal direction. We use Jacobi quadrature with parameters (0, 1/2) for the radial direction. After obtaining the grid, we compute all higher-order spin-weighted spherical harmonics and generalised spherical Zernike polynomials out of the three-term recursion for Jacobi polynomials. Each basis function comprises a spatial envelope (e.g., r l ) multiplied by a Jacobi polynomial of some kind (e.g., P (0,l+1/2) n (2r 2 − 1) ) We initialise with the appropriate spatial envelope and recurse up to the desired degree from there. In some cases the spatial envelope contains an extreme dynamic range (e.g., r l for l 1). This can lead the initialization to underflow to zero. Grid points where this occurs can never return to finite values, even though the eventual basis function should end up moderate at such points. To avoid underflow problems, we used 128-bit precision for the construction phase of the transform grid, weights and matrices. We cast the results to 64-bit precision after the initial construction. There are more sophisticated methods available to handle the same problems, e.g., the recursions presented in [17], but our simple method works for all polynomial degrees up to roughly ≈ 10, 000. The speed of the construction is fast enough considering we store the transform matrices for later use.
We can transform data from physical space to coefficient space, and back, either in serial or parallelized across cores (using MPI). We can parallelize in either one or two directions for 3D calculations. If the data size is O(N ) in each direction, parallelization across two directions allows a calculation to be run efficiently on O(N 2 ) cores. We describe the transform assuming parallelization across two directions, but note how the calculation differs if parallelized in a single direction.
In physical space, each core has the data for all φ points, but for only a subset of points in θ and r (or r only for parallelization across one direction). We say that the data are local in φ, but distributed across θ and r. First we use Dedalus to perform a real to complex Fourier transform in φ, so we have T m (θ, r). If the data is distributed across cores in θ, we use Dedalus to perform a parallel transpose across m and θ so the data are local in θ and distributed across m. Each processor has data for a sequential subset of m values.
Next we multiply by the unitary matrix U † to transform from the components i to spins σ. For each m and spin σ, we use a matrix multiplication transform to calculate the spinweighted spherical harmonic coefficientŝ where the matrices Sσ l (θ i ) represent whichever spin-weighted transform matrix is appropriate at the time; each matrix has size (L max − L min + 1) × N θ , where L min = max(|m|, |σ|).
Although there is less data for higher m, we pad with zeros soT σ m, (r) is a rectangular N c × (L max + 1) × (L max + 1) × N r array. The S matrix is the product of a spin-weighted spherical harmonic function and Gaussian quadrature weights. The inverse transform matrix is simply a spin-weighted spherical harmonic function, properly transposed [2]. We next use Dedalus to perform a parallel transpose across and r so the data are local in r and distributed across m and (or only if parallelized across one direction). Each processor has data for a sequential subset of the values.
We next multiply by the -dependent orthogonal matrix Q to transform from the spins σ to regularities a. This givesT a m, (r). Then for each we use a matrix multiplication transform to calculate the Jacobi polynomial coefficients, where w(r i ) is the Gaussian quadrature weight, and Q l+ā n (r i ) represents whichever Jacobipolynomial transform matrix is appropriate at the time; each matrix has size (N max − N min + 1) × N r , where N min = | − R max |/2 . The data are then stored as a list of arrays with N C (N max −N min +1)×N m elements for each value of . The forward transform is the product of a Zernike polynomial and Gaussian quadrature weight. The inverse radial transform matrix is simply a Zernike polynomial, properly transposed.
This gives the coefficient expansion,T a m, ,n associated with the tensor T(r, θ, φ). As each step is invertible, the algorithm can be inverted to calculate the grid values T(r, θ, φ) which correspond to coefficientsT a m, ,n .

Eigenvalue Problems
In this section, we solve eigenvalue problems of the form, subject to boundary conditions where M , L, and B are linear operators, and λ is the eigenvalue.

Surface Rossby Waves
Before discussing problems in the full sphere, we briefly mention an example of two-dimensional flow on the surface of a sphere. In this case, we expand the solution only in spin-weighted spherical harmonics. We consider Rossby waves. We solve the Laplace-tidal equations on the sphere r = 1, where we normalize the problem such the Coriolis parameter is equal to unity. Lamb's parameter γ = 4Ω 2 a 2 /gH, which we take to be zero. The statevector is where u − and u + are the two spin components of the angular velocities. The linear operators M and L are where again we take γ = 0.
We expand p, and u ± in 512 spin-weighted spherical harmonics each. See Part-I for the definitions of the cosine operators C and the angular derivatives k ± . The problem is coupled in but not in m. Thus, we can solve each m independently. We pick m = 50 for illustrative purposes. The generalized eigenvalue problem is solved using scipy's eig routine. The analytical eigenvalues are This is a non-trivial problem because eigen-solutions cannot be represented in terms of a pure spin-weighted spherical harmonic; multiplication by cos(θ) complicates the situation. However, the solution is expressible in a simple finite combination of spherical harmonics. We demonstrate in figure 1 that we correctly find all the eigenvalues.

Spherical Bessel Equation
We next solve the scalar-Laplacian eigenvalue equation with the boundary condition, This is an eigenvalue problem with eigenvalue κ 2 . We posed this problem as an example of our approach in Part-I, here we provide the full numerical solution. The equation is separable into a radial and an angular component. The scalar function f can be expanded in scalar spherical harmonics, R ,m (r)Y 0 ,m (θ, φ). Then R ,m satisfies the spherical Bessel equation The solutions are spherical Bessel functions of the first kind, j (κr). The boundary condition at r = 1 requires κ to be a zero, i.e., j (κ) = 0.
To solve this numerically, we take the state vector X = f , L = D − 1, +1 D + 0, , and the eigenvalue λ = κ 2 . Because L.X is in the α = 2 function space, owing to the two derivative operators, in equation (8) the matrix M = C 1, C 0, , where the C α, are conversion matrices which increment α → α + 1. Paper-I discusses the details of how we construct a matrix system of equations from the original PDE (15) We expand f in a basis of 512 polynomials. To implement the boundary condition, the last row of L is replaced with the r = 1 restriction operator Q α, (r = 1), a row vector of the Q polynomials evaluated at r = 1, and the last row of M is replaced by zeros. The generalized eigenvalue problem is solved using scipy's eig routine. The eigenvectors are transformed to the grid to compare to the spherical Bessel function.
In figure 2, we plot the 100 th eigenmode solution to equation (15) with = 50 (top panel), along with the error (|f − j (κr)|; bottom panel). The inset shows that near the origin, f ∼ r 50 , as required by the regularity condition at r = 0. This regularity condition is satisfied automatically by our choice of radial basis. Figure 1 plots the fractional error in the eigenvalue κ, where the analytic eigenvalues are the zeros j (κ a ) = 0. As expected, about the first half of the eigenvalues are very accurate. Eigenvalues corresponding to eigenmodes with high radial wavenumbers tend to have higher errors; these eigenvalues can be computed to machine precision by increasing the number of polynomials in the basis.

Linear Diffusion Equation & Boundary Conditions
In this section, we solve the linearized, diffusive hydrodynamics equations We assume ∂ t u = −κ 2 u, and solve for the damping rate κ 2 . Here we denote the pressure with p. This is completely equivalent to the linearized equation for the magnetic vector potential Here Φ is the scalar potential. The general solution can be derived analytically, so we can implement a wide range of boundary conditions and calculate exact solutions. This makes this problem very useful for insuring the proper implementation of boundary conditions in the code.
The hydrodynamics problem can have no-slip boundaries, or stress-free boundaries at r = 1.
where we have also assumed impenetrability and where is the rank-2 stress tensor.
There are also several choices for magnetic boundary conditions. Potential boundary conditions assume the magnetic field matches onto a harmonic field for r > 1. This is a non-local condition that is commonly specified by decomposing A into spherical harmonic degrees, A . A perfectly-conducting boundary has no normal magnetic field and no tangential electric fields. The pseudo-vacuum boundary condition is that the tangential magnetic field is zero.
We solve each of these problems analytically in appendix A. In each case, the eigenvalues are related to the zeros of spherical Bessel functions of different orders.
In this problem the statevector is where u − , u 0 , and u + are the components of u in regularity classes. The linear operators M and L are For no-slip boundary conditions, all three components of u are zero at the boundary. For stress-free boundary conditions, to impose no normal flow, we set where Q is the rank one orthogonal matrix and Q α, +a (r = 1), the restriction operator, is a row vector of the Q polynomials evaluated at r = 1. In equation (31), α = 0. The two other conditions are equivalent to Thus we impose +1 a,b=−1 Note that here we must take α = 1 because the D operator increases α from 0 to 1.
For the magnetic problem, the M and L matrices are identical, but the statevector changes to The magnetic boundary conditions have a simple form in terms of regularities: Although potential and pseudo-vacuum boundary conditions look similar, the pseudo-vacuum conditions can be expressed locally (equation (27)), whereas the potential conditions cannot (equation (25)).
We apply boundary conditions using the tau method (see Part-I & references within). Boundary conditions are imposed by adding a correction term (called τ ) to the equations. One arrives at (often subtly) different answers depending on the assumed form of the correction term. We assume τ takes the form of one of our basis polynomials, Q α, +ā n (r), where n is the highest order radial mode for the chosen values of , N max , and R max .
There is a remaining choice for what value of α to use. We call this value α BC . We use either α BC = 2 or α BC = 0. Using α BC = 2 is equivalent to replacing the last row of the L matrix with the boundary condition, as we did in section 4.2.
For α BC = 0, we add a single extra element to the state vector for each boundary condition, which corresponds to each τ correction. Then we must add extra rows to the matrices, which are the boundary conditions. We also must add extra columns to maintain square matrices. The column associated with a given τ is given by the final column of the C 1, C 0, matrix for the equation corresponding to the τ error (note it is dependent). The extra columns are only non-zero for the variables for which we are applying boundary conditions, e.g., the divergence conditions are imposed exactly. We have checked that if we instead use the final column of the identity matrix, we find the same results as for α BC = 2.
To solve the problem numerically, we fix = 50 and use 256 terms in the radial expansion of each variable in the statevector. We also apply boundary conditions using α BC = 2 and α BC = 0. Thus, for α BC = 2, the M and L matrices have size 1024 2 , whereas for α BC = 0 they have size 1027 2 because of the three extra rows and columns to incorporate the τ errors. We solve this generalized eigenvalue problem with scipy's eig routine.
For each of these five sets of boundary conditions, we have two sets of eigenvalues because each problem decouples into a problem for the toroidal component, and a problem for the In each case, we find that about half the eigenvalues are very accurate. When we use α BC = 0 to apply the boundary conditions, we find that the number of accurate eigenvalues is about 10% greater than when we use α BC = 2.
poloidal component (Appendix A). We calculate the analytic values of the eigenvalues, κ a , using the formulae in Appendix A. We sort the numerical and analytic eigenvalues and compare them. Figure 3 shows that about half of our eigenvalues are accurate. The errors using α BC = 2 and α BC = 0 are similar for all boundary conditions except perfectly-conducting, where α BC = 0 is more accurate. However, for each case, when we use α BC = 0 there are about ∼ 20 extra accurate eigenvalues compared to α BC = 2. This suggests that using α BC = 0, we are able to correctly resolve smaller scale features at a given resolution. Note that since our resolution is 256, we expect to have around 512 eigenvalues because there are two independent solutions.

Boundary Value Problems
We next discuss the solution of boundary value problems. We solve the equation subject to the boundary conditions where L and B are linear operators, R is a state vector, and E is a state vector restricted to r = 1.
In section 6.4, we initialize the magnetic vector potential from a specified magnetic field. This is a boundary value problem, and we use it as an example. Specifically, we are solving Thus, our state vector is If we set L equal to the curl operator, then we would have However, this gives two redundant equations for A 0 , and cannot uniquely determine A ± because we have not set the gauge. Thus, we replace the third row with the gauge condition ∇ · A = 0, The right hand side vector is given by B 0 , but we need to multiply by a conversion matrix C to increase the α index to 1, since L and L are both α = 1 (all terms carry a derivative operator): Finally, we must apply boundary conditions. We use potential boundary conditions, which we impose in the last rows of the three components of the L matrix (α BC = 2).
The dynamo problem in section 6.4 starts with an initial magnetic field To solve the boundary value problem numerically, we represent B 0 with N max = 31, L max = 31, R max = 2, and no dealiasing. We invert the L matrix for each using scipy's sparse solver (sparse.linalg.spsolve).
For this simple problem, we can solve the problem analytically and compare to the numerical solution. The cleanest way to write A in the Coulomb gauge is in terms of a poloidal function, where and One can check that ∇ × A analytic = B 0 and that A analytic satisfies the potential boundary conditions (equation (36)).
To validate our numerical solution, we calculate the error where i = r, θ, φ and the maximum is across all grid points. We find the error is 2.2 × 10 −14 , 2.7 × 10 −15 , and 2.9 × 10 −15 for the three components r, θ, and φ. Thus we have an accurate solution to this boundary value problem.

Initial Value Problems
We now discuss the three benchmark problems of M14. The three problems are posed as initial value problems in the form of equation (1). As the code described in this paper is an extension of the Dedalus code, we refer to it with D. We compare our results to the Marti & Jackson code [12], and the Hollerbach code [7,8]. The Marti & Jackson code (hereafter MJ) decomposes all variables into scalar functions, and then expands the scalar functions in scalar spherical harmonics in the angular directions, and in Jacobi polynomials weighted by r in the radial direction. The Hollerbach code (hereafter H) also decomposes all variables into scalar functions, and uses scalar spherical harmonics in the angular directions, but uses Chebyshev polynomials in the radial direction. Thus, H does not explicitly impose the regularity conditions at r = 0, unlike MJ and D.

Timestepping
We timestep equations of the form We use multistep implicit-explicit (IMEX) methods. Terms on the left hand side of equation (53) are treated by linearly-implicit methods and must be linear in the evolution variables X, while terms on the right hand side (F (X)) are treated explicitly and can include both linear or nonlinear terms. For a general multistep IMEX integrator, the new statevector at time n + 1 is related to the statevector at earlier times by We use two different time-stepping schemes: the second-order, two-step Crank-Nicolson-Adams-Bashforth scheme, CNAB2; and the fourth order, four-step semi-implicit backwards differencing formula scheme, SBDF4 [both described in 20]. Although M14 does not precisely state what time stepper MJ or H use, they likely used the second-order Runge-Kutta scheme described in [7] or [12, although it is unclear how many iterations were used]. We refer to this time-stepping scheme as RK2.
Two of the benchmark solutions are not stationary, and we find that the choice of timestepper and timestep size plays an important role in resolving the solutions. There is no discussion of how the timestep size is chosen in M14. Although we have run simulations with adaptive timestepping, to simplify our results and enhance reproducibility, we only report simulations with constant timesteps.

Resolution and Degrees of Freedom
It is not trivial to compare the resolutions in D simulations to resolutions in MJ or H simulations. This is because we use a triangular truncation in the radial direction, unlike MJ or H. We report our radial resolution in terms of N max . However, the number of radial modes averaged over is roughly For instance, if N max + 1 = 1 2 (L max + 1), then the average number of radial modes is about 1 2 (N max +1). This is complicated slightly by the regularity dependence of radial modes. Also, large 's are associated with fewer m modes than small 's.
In contrast, MJ do not appear to use a triangular truncation, and instead appear to use a constant number of radial modes, N r , for every . This means that the maximal radial order depends on . If N r = 1 2 (L max + 1), as is often the case in M14, then the highest mode is a polynomial with order 3N r . Thus, one would require ≈ 3N r grid points to prevent aliasing errors, rather than the ≈ 3N r /2 grid points required to dealias when using the triangular truncation.
To make a fair comparison between the codes, we report two quantities related to the number of radial modes. First, we report N max , which is half the highest radial order. This is analogous to reporting the angular resolution in terms of L max . For H, N max is equal to the number of radial modes. We also report the total number of spatial degrees of freedom, or DoF.

Energy Calculations
A main output of the benchmark problems are the energies of the equilibrated states. However, one must take care to accurately calculate volume integrals of quantities like the energy, lest error in the volume integral itself dominate the reported results. The weights of the scalar spherical harmonics is sin(θ), and since we use their quadrature nodes for the θ grid, we can use their quadrature weights to calculate angular integrals with spectral accuracy. Similarly, the weights of the Q polynomials is r 2 , and the quadrature weights can be used again to calculate radial integrals with spectral accuracy. Explicitly, we calculate the kinetic energy with where the sum is over each point on the r = (φ, θ, r) grid. The quadrature weights are w φ = 2π/N φ , and w r , w θ derived from their respective Gaussian quadrature.
In contrast, the polynomials used in H & MJ have an integral weight of (1 − r 2 ) −1/2 . Thus, to calculate integrals via quadrature, one must also include a factor of r 2 √ 1 − r 2 in the sum in equation (56). However, √ 1 − r 2 is not analytic at r = 1, so this reduces the accuracy of the integration scheme to second order. It is possible to have very accurate solutions, but to report inaccurate energies due to a low order integration scheme. In the hydrodynamics benchmark problem, H & D converge to the solutions with the same energy at fairly low resolution. This suggests H is not using this quadrature scheme to calculate the energy. However, MJ converges much more slowly. We hypothesize this is not due to inaccuracies in their solution, but instead due to inaccuracies in their integration scheme used for measuring KE.

Hydrodynamics Problem
The simplest problem in M14 solves the incompressible hydrodynamics problem with imposed velocity boundary conditions (benchmark 3), The terms on the left of the equals sign are timestepped implicitly, whereas the terms on the right of the equals sign are timestepped explicitly. The boundary conditions at r = 1 are u = u 0 , where Following M14, we take ν = 10 −2 , Ω = 10, and u 0 = 3/(2π).
We evolve the variables u − , u 0 , u + , and p, where u is written in terms of the three regularity classes. We report the matrices used for the problem in Appendix B.
The simulation is initialized with zero initial flow, but quickly reaches a stationary equilibrium state. To quantitatively describe this state, we calculate the volume-integrated kinetic energy, Figure 4 shows the kinetic energy as a function of time. We run simulations to t = 40 so we can minimize any transient effects from our initial condition. The simulations are run with N max = L max and R max = 3. The low-resolution simulations were dealiased, but the higher resolution simulations were not because they are already well-resolved. We ran with both α BC = 0 and α BC = 2. High resolution simulations gave identical energies for both values of α BC . For timestepping we use a CNAB2 with a constant timestep ∆t of 2 × 10 −2 or 10 −2 .
Even at late times, there are small changes in the kinetic energy. However, the kinetic energy is constant to ten decimal places between t = 35 and t = 40, so we report the values to ten  We find temporal and spatial convergence to ten decimal places, and find KE h = 0.06183074756. We can reach this converged solution at a resolution of N max = L max = 23 and timestep size of ∆t = 0.02. At high resolutions, we do not find any differences in the kinetic energy in simulations with or without dealiasing, or with α BC = 0 or α BC = 2. At low resolutions, we find that the energies were closer to KE h when we used α BC = 0 than when we used α BC = 2. The algorithm of MJ corresponds to α BC = −1/2, whereas there is no equivalent parameter for H.
We also report the kinetic energy of several MJ and H simulations reported in M14. They do not report their timestep size. Our simulations are consistent with the values reported by H. With similar numbers of degrees of freedom, H appears to be slightly more accurate than our D simulations. On the other hand, our simulations appear to be significantly more accurate than those of MJ. For this problem, we find no difference in the energy of the stationary state for simulations with different timesteps. This is not the case for the next two problems.

Convection Problem
Next we consider the rotating convection problem of M14 (benchmark 1). For this problem, the equations are The Ekman number E = 3 × 10 −4 , the Rayleigh number Ra = 95, the Prandtl number Pr = 1, and the temperature source term S = 3. The temperature has an equilibrium base state of T = 0.5(1 − r 2 ). The vector r ≡ re r represents the full radial vector, and the gravity is linearly increasing, as is appropriate to a self-gravitating incompressible sphere like the Earth's core. To reach the appropriate solution, we initialize the problem with the temperature initial condition specified in M14: The initial velocity is taken to be zero. The boundary conditions are impenetrable and stress-free for the velocity, and fixed temperature.
We evolve the variables u − , u 0 , u + , p, and T . As for the hydrodynamics benchmark, u is written in terms of the three regularity classes. We report the matrices used for this problem in Appendix C.
With this initial condition, the fluid is expected to evolve to a traveling wave state with constant kinetic energy. As in the hydrodynamics problem, we report the volume-integrated kinetic energy KE in Table 2. The simulations are run with N max = L max and R max = 3. The low-resolution simulations were dealiased, but the higher resolution simulations were Both images are created from the same vantage point. Volume and streamline renderings created using Vapor [4,5].
not because they are already well-resolved. We ran with both α BC = 0 and α BC = 2. For timestepping, we use either CNAB2 or SBDF4 with constant timesteps. We run all simulations for 20 diffusion times.
The structure of the equilibrated-flow, shown in Figure 5, is an m = 3 travelling wave. The structure and amplitude of the flow match those shown in M14. The radial u r and azimuthal u φ flow are symmetric across the equator while the latitudinal u θ flow is antisymmetric about the equator; as such we show a half-hemisphere, containing the equator and the south pole, in the volume and streamline renderings. The spiralling nature of the flow is visible in the streamline rendering of Figure 5, with the flow dominated by u r and u φ and with slower flows along the rotation axis. This snapshot is taken from the equilibrated state at t = 20 of our D simulation with L max = 31, N max = 31, α BC = 2, and ∆t = 10 −5 with the SBDF4 timestepper. Figure 6 shows the kinetic energy as a function of time. The left panel shows that after a few diffusion times, the kinetic energy becomes approximately constant, indicating that we have reached the traveling wave solution. However, the right panel shows that there is different secular behavior for different simulations. In our simulations with larger timestep size (∆t ≥ 2 × 10 −5 ) and α BC = 2, we find a secular energy growth of ∼ 10 −9 over the course of the simulation, similar to the red curve. The simulations with smaller timestep size or α BC = 0 show regular (blue curve) or irregular (green curve) low amplitude oscillations. Because of these oscillations and secular variation, we report the kinetic energy at t = 20 to ten decimal places.
Although the kinetic energy of the traveling wave is close to constant, the simulation must resolve the advection of the wave around the domain. We find that the choice of timestepper and timestep size plays an important role in determining the kinetic energy of the traveling wave state. In table 2 we report the kinetic energy of the traveling wave with different simulation parameters. We achieve spatial and temporal convergence to ten decimal places,    and find the kinetic energy to be KE c = 29.12045489. We can reach this level of accuracy with the fourth order timestepper SBDF4 with timestep size of 2 × 10 −5 and a spatial resolution of N max = L max = 31 with α BC = 2. At low resolution, we find that simulations with α BC = 0 reach convective states with energies closer to KE c than simulations with α BC = 2.
However, our value of the kinetic energy is inconsistent with the H and MJ values reported in M14. We believe the discrepancy is due to timestepping. Both use a second order scheme, RK2, which is less accurate than SBDF4. To test this, we ran simulations with the second order timestepper, CNAB2.
Simulations with the second order timestepper do not show temporal convergence (to ten decimal places) with timestep sizes greater than or equal to 10 −5 . We find that the kinetic energies are converging to KE c like ∆t 2 ( figure 7). This indicates that the dominant errors in the simulations are due to timestepping. The kinetic energies reported by Marti and Hollerbach are similar to the kinetic energies of our simulations with timestep size 10 −5 and 2 × 10 −5 . Thus, our results are consistent with those reported by the Marti and Hollerbach codes when we use a low order timestepper. Since the CNAB2 timestepper has some well known flaws [1], we also test the L-stable Modified CNAB2 (MCNAB2) timestepper of [20] at ∆t = 10 −5 . Evidently, the flaws are minor at worst.
This sensitivity of the kinetic energy to the details of timestepping make it difficult to determine the accuracy of the spatial discretization using this problem. Temporal convergence studies, as we have done here, are necessary to distinguish spatial from temporal errors.

Dynamo Problem
The last problem we discuss is the rotating convective dynamo problem of M14 (benchmark 2). This is the most challenging of the three benchmark problems. We solve the equations where A is the magnetic vector potential in the Coulomb gauge, and can be used to calculate the magnetic field B using ∇ × A = B. The scalar potential Φ enforces the Coulomb gauge constraint. We solve for u − , u 0 , and u + (the three regularity components of u), p and T , A − , A 0 , and A + (the three regularity components of A), and Φ. The exact implementation is described in Appendix D. As in M14, we set the magnetic Rossby number Ro = 5 7 ×10 −4 , the Ekman number E = 5×10 −4 , the Roberts number q = 7, the Rayleigh number Ra = 200, and the temperature source term to S = 3q = 21. For boundary conditions we use impenetrable, stress-free boundary conditions for the velocity, fixed temperature, and potential boundary conditions for the magnetic field (see section 4.3).
We use the temperature and velocity initial condition described in M14. Because we evolve the magnetic vector potential and not the magnetic field directly, we solve a boundary value problem to initialize the vector potential (see section 5). We run simulations with N max = L max and R max = 3. The low-resolutions simulations were dealiased, but the higher resolution simulations were not because they are already well-resolved. We ran with both α BC = 0 and α BC = 2, and again we found more accurate solutions at low resolutions using α BC = 0. For timestepping, we use either CNAB2 or SBDF4 with constant timesteps.
The system approaches an oscillating dynamo solution. The kinetic and magnetic energy, undergo regular variations over the oscillation period. We plot the kinetic and magnetic energy in figure 8. After an initial transient of ∼ 2 magnetic diffusion times, the system approaches an oscillating dynamo solution. The kinetic and magnetic energy oscillate rapidly in this state. For some of the highest resolution simulations, we restarted the simulation from a low-resolution simulation evolved beyond the initial transient.
A volume rendering of a late state of the oscillating dynamo simulation with N max = L max = 63, α BC = 2, no dealiasing, and timestep size of 2.5 × 10 −6 using the SBDF4 scheme is shown in Figure 9. Shown are the radial velocity field u r in half-domain rendering (as in Figure 5), and magnetic field line rendering in the full domain. In contrast to the rotating convection problem, the dynamo problem has asymmetries in the velocity field, with a region of strong up and downflow visible near the front of the rendering. The strong fields themselves are We use the SBDF4 timestepper with timesteps of 2.5 × 10 −6 . The kinetic and magnetic energies oscillate rapidly, which leads to an extended opaque region on the plot. In each inset, we zoom in onto the top or bottom of each oscillatory region, and plot the deviation of the energies from their limit supremathey might not go to zero because we round the limit suprema. For the kinetic energy, we use KE inf = 33681.31 and KE sup = 37444.32 (to within an accuracy of 10 −2 ). For the magnetic energy, we use ME inf = 867.7413 and ME sup = 943.4111 (to within an accuracy of 10 −4 ). These can be combined to calculate KE, ∆KE, etc. Figure 9: Volume rendering of radial velocity (left) and magnetic field lines (right) at a late time in the dynamo problem. In the radial velocity u r volume rendering (left), the upper half of the sphere has been cut away, showing an equatorial slice with columns descending below to the south pole and a marked asymmetry in flow structures. The transfer function is asymmetric, but with both color ranges diverging from zero flow. In the magnetic field line rendering (right), the full volume is shown from the same vantage point as the radial flow rendering. Field lines are seeded in the strongest regions of |B φ |, and field lines are colored by the phi component B φ . Here, in contrast to u r , the transfer function and the field B φ are both symmetric. The blue arrow is aligned with the rotation axis, pointing north. Volume and field rendering created using Vapor [4,5] from the same vantage point.
found in the region of strong flow. These volume and field line renderings were created using Vapor 2 . This snapshot is taken from the equilibrated oscillating state at t = 5.5 of our high resolution D simulation with L max = 127, N max = 63, α BC = 2, and ∆t = 1.25 × 10 −6 with the SBDF4 timestepper. At t = 5.5, KE(5.5) = 33860.4 and ME(5.5) = 933.404 and the solution is the decreasing magnetic energy phase of the oscillation.

Energy Diagnostics
Because the energy is not constant, there are different ways to characterize the system. M14 decomposes the energy into its two dominant temporal Fourier modes Most of the energy are in these modes, and the higher harmonics of f . The 2f harmonic contains a few percent of the energy of the time series. To perform the decomposition, we take each energy time series over the final half a diffusion time of the simulation. Then we identify the first and last energy maximum in the time series. We truncate the time series so it ranges from the time of the first maximum to the last output time before the last maximum. This makes the time series approximately periodic. We then take the Fourier transform. The amplitudes of the first and second peak give C and A.
Unfortunately, this approach is sensitive to various choices in the algorithm. For instance: using the energy time series for a full diffusion time or a half of a diffusion time; running the algorithm with a different output cadence; or, different truncation methods, all give different coefficients C and A. It is difficult to calculate the first two coefficients of an expansion to high precision (e.g., 10 −7 ), when the third coefficient has size 10 −2 .
We also consider a new, more robust metric to characterize the system. The simulation approaches an oscillatory state as t → ∞. The dynamo oscillation has well defined minima and maxima of kinetic and magnetic energy. Thus, we characterize the solution by calculating these minima and maxima. Formally these are the limit extrema as t → ∞. We define We refer to the limit extrema as KE inf , KE sup , and similar for the magnetic energy. ME ME α BC = 2 α BC = 3/2 α BC = 0 Figure 10: The magnetic energy as a function of time for three dynamo simulations with N max = L max = 23, but with different implementation of boundary conditions (α BC = 0, 3/2, or 2). Despite only very minor changes in the numerical algorithm, we find completely different solutions at these low, unresolved resolutions. Although α BC = 3/2 appears closest to the correct solution (black line), its magnetic energy is decaying secularly with time-we hypothesize its magnetic energy trends to zero as t → ∞.
To determine the limit extrema energies, we include insets in figure 8 in which we zoom in to energy scales close to the limit extrema. Here we can see a regular pattern in the energy extrema. This is because the, e.g., maximal energy might occur between time steps, so we must integrate for many oscillation periods before our integration reaches a time which, by coincidence, is very close to the time of an energy extremum. The limit extrema are rounded to 10 −2 for kinetic energy and 10 −4 for magnetic energy (7 digits reported in both cases). We show below that these metrics are more robust than the (C k , A k ) and (C m , A m ) decompositions described in M14.

Low-resolution Solutions
We found the low-resolution calculations (N max and L max less than 63) are extremely sensitive to numerical details. We perform a set of simulations with N max = L max = 23, using the CNAB2 timestepper with timestep size of 5 × 10 −6 . We use three different values of α BC : 0, 3/2, and 2. The magnetic energy of each solution is shown in figure 10. For comparison, we also show M E = 905.566 (the average magnetic energy at late times) as calculated in high resolution simulations (see table 4).
Although almost all simulation parameters are the same, we find completely different solutions with different values of α BC . It may appear that the most accurate solution uses α BC = 3/2. However, a more careful inspection of the data shows that the magnetic energy decays secularly by about 2.5 energy units every magnetic diffusion time. The energy in the simulation with α BC = 2 also decays at late times. We suspect that in both simulations, the magnetic energy decays to zero as t → ∞, i.e., they do not represent dynamo solutions. We also find that the magnetic energy decays slowly when using α BC = 2 at resolutions  N max = L max = 31 and 47, although the decay rate becomes smaller as the resolution increases. As shown in Figure 8, the energy asymptotes to a constant oscillate at the medium resolution of 63. In contrast, the magnetic energy in the simulation with α BC = 0 stays constant at late times at all resolutions we tried.
As described above, we found α BC = 0 is consistently more accurate than α BC = 2 at low resolutions in all our tests. For the other problems, the difference between α BC = 0 and α BC = 2 was at most a few extra digits of accuracy. However, for this problem, there is an order unity difference between the simulations with α BC = 0 and α BC = 2. For this reason, we think the low resolution (N max = L max = 23) dynamo problem is a good numerical test. It can show the limitations of a given numerical scheme and sensitivity to different methods in marginally resolved simulations.

Convergence Study
We report the kinetic and magnetic energies in our simulations in tables 3 & 4. Our high resolution simulations are run without dealiasing, and with a higher order timestepper. We find no significant difference between α BC = 0 and α BC = 2. Using the decomposition of the kinetic and magnetic energy into the first two Fourier modes, our simulations appear to show convergence to 5 decimal places in C k , 4 decimal places in A k , 5 decimal places in C m  Table 4: The normalized sums and differences of the limit extrema of the kinetic and magnetic energy (equations 74-77) for the convective dynamo test problem of M14. The correct digits for each solution are underlined. We find spatial and temporal convergence to six or seven digits for each quantity. and 4 decimal places in A m . Subsequent digits of each quantity vary with different spatial or temporal resolution. However, these differences are primarily due to the decomposition algorithm, rather than differences in the actual dynamo solution.
The low-resolution Dedalus simulations have similar accuracies as H and MJ. For the lowest resolution (N max = L max = 23), our solution is more accurate than MJ, but less accurate than H. Of course, all three simulations are far from the correct solution. At a resolution of N max = L max = 31, our solution has very similar energy to MJ-H appears to be more accurate for the magnetic energy, but less accurate for the kinetic energy. Higher resolutions cannot be easily compared to a reference solutions because the differences in energies between the different simulations are likely due to differences in the decomposition into C and A rather than real differences in the solutions. Table 4 shows the normalized sums and differences of the limit extrema of the kinetic and magnetic energy. We can only report Dedalus simulations as this quantity was not reported in M14. Note that each quantity is different from their corresponding values in table 3 by a few percent. This is due to the effect of higher order harmonics dropped in equations 72 & 73.
Our fiducial simulation has resolution of N max = L max = 63, α BC = 2, is run without dealiasing, and uses the SBDF4 timestepper with constant timestep size of 2.5 × 10 −6 . If we increase the radial resolution to N max = 127, none of the quantities change to the accuracy reported. This indicates that the simulation is radially well-resolved to this level of accuracy with N max = 63. We then fix N max = 63 and increase the angular resolution. We find the same value for KE and ∆ME to all digits reported between L max = 95 and L max = 127. However, ∆KE differs by 10 −2 and ME differs by 3 × 10 −4 . Finally, we checked temporal convergence by running our highest resolution simulation with time- Simulations with higher angular resolution (or run with dealiasing) could increase the accuracy of ∆KE and ME. The other quantities cannot be determined to higher accuracy without running the simulations for longer to minimize the effects of transients (see figure 8).

Conclusions
This paper describes the implementation of a new method for the solution of a wide class of partial differential equations in a full sphere, described in Part-I. We represent tensor variables using spin-weighted spherical harmonics in the angular direction, and a scaled class of Jacobi polynomials in the radial direction. This ensures that each quantity satisfies the appropriate regularity conditions, both at the poles, and at r = 0. We can calculate nonlinear quantities by transforming the solution from spectral space to physical space, and performing products or other operations in physical space.
To demonstrate the accuracy of this method, we first discuss a series of unit tests which test specific aspects of the code. The first is the Bessel's equation eigenvalue problem (section 4.2). This is a non-trivial problem, as Bessel functions are not polynomials. The accurate solution of Bessel's equation thus demonstrates the exponential convergence of our algorithm.
We next solve for the decaying eigenmodes of a diffusing, divergence-free vector field (section 4.3). Here we use the tensorial nature of our algorithm to rewrite the three components of the vector into three regularity classes, each of which has different behavior as r → 0. We are able to impose different boundary conditions on the eigenvalue problem: no slip or stress-free boundary conditions if the vector field is the velocity; potential, conducting, or pseudo-vacuum boundary conditions if the vector field is the magnetic vector potential associated with a magnetic field in the Coulomb gauge. In each case, we can solve the eigenvalue problem analytically in terms of Bessel functions. This allows us to precisely compare our eigenvalues to the analytical eigenvalues, and demonstrates that we correctly impose all boundary conditions.
We include an example of the solution of a boundary value problem (section 5). We solve for the magnetic vector potential which corresponds to a specified magnetic field. The magnetic field is a polynomial in r and only includes a few spherical harmonic components in the angular direction. This problem is a useful test of our transforms. These unit tests are invaluable for code verification.
The remainder of the paper describes our solutions to three full-code, initial value problems proposed in M14. The first problem we consider is a hydrodynamics problem (M14's benchmark 3; section 6.2). The system is forced with a velocity boundary condition, which leads to a stationary, nonlinear equilibrium. Crucially, the solution is independent of the details of timestepping, which is not the case for the other two problems. This makes the problem an excellent test of the spatial discretization. It also ensures the code can correctly impose the boundary conditions, and calculate the Coriolis force and the u · ∇u nonlinearity. We find excellent agreement with Hollerbach's simulations (M14). We show spatial and temporal convergence of the kinetic energy to 10 digits of precision: KE h = 0.06183074756.
The second problem is a rotating convection problem (M14's benchmark 1; section 6.3).
The system evolves toward a traveling wave solution, whose kinetic energy is constant in time. We find that the saturated value of the kinetic energy depends on the timestepping scheme. We show that for moderate spatial resolution, the error in the kinetic energy can be dominated by timestepping errors, and decreases like ∆t 2 for a second-order timestepper ( figure 7). Thus, we find that this problem is more of a test of a code's timestepper, rather than its spatial discretization. This makes it difficult to compare to previous results in M14, as they do not provide sufficient details about their timestepping. Nevertheless, we find similar results to the Hollerbach and Marti-Jackson codes when we run with a secondorder accurate timestepper. Switching to a fourth-order accurate timestepper, we are able to show spatial and temporal convergence of the kinetic energy to 10 digits of precision: The last problem is a convective dynamo problem (M14's benchmark 2; section 6.4). This is the most challenging problem in M14. This problem is very sensitive to the numerical method at low resolution, where we find that minor numerical choices can lead to completely different solutions (see figure 10).
The desired end state for this problem is an oscillating dynamo solution, for which both the kinetic and magnetic energy are variable. Thus, one must decide how to characterize the solutions. In M14, the kinetic and magnetic energy time series were decomposed in a series expansion, and they report the first two terms in the series. This is not very precise, as the third term in the expansion has a relative size of O(10 −2 ), and different methods of carrying out the decomposition may give systematically different results. Nevertheless, if we implement this decomposition, we recover solutions similar to those in M14. It is difficult to tell if differences between the codes are due to differences in the spatial discretization, the temporal discretization, or the algorithm used to preform the energy decomposition.
We also discuss a new diagnostic for this convective dynamo problem which is very precise. We calculate the limit superior and limit inferior of the kinetic and magnetic energies, and then calculate their sums and differences. This considers all possible terms in the series expansion considered in M14 and gives very precise solutions. We find temporally and spatially converged solutions to 5 or 6 digits of precision, as reported in equations 78-81. We hope this new diagnostic makes this problem more useful for quantitative code comparison.
We found it difficult to compare to the solutions of the Hollerbach and Marti-Jackson code because M14 does not include some important details of the simulations, e.g., timestepping scheme or timestep size, all the details of the algorithm for calculating the volume-integrated energy, etc. This makes it unclear if different results are due to important differences between spatial discretization schemes, or simply due to different timestep sizes. To aid future comparison to the solutions we describe here, we include the source code used to run the simulations, as well as the data and analysis scripts used to generate the plots. We believe that this information makes future code comparisons much more fruitful.
We have implemented the algorithms of Part-I using aspects of the Dedalus code. In the future, we will fully incorporate spherical geometry into the main Dedalus codebase. This will allow us to use the Dedalus equation parser to write equations out as strings, rather than manually constructing matrices and the nonlinear terms. This will also allow the user to specify complicated simulation outputs (e.g., enstrophy, Reynolds stresses, etc.) in string form. These features will make this algorithm straightforward to use for the solution of many different PDEs in spherical geometry.
for some coefficients c l and shift in the regularity class parameter, l + a. The spherical Bessel function are easy to evaluate with standard packages; e.g., scipy. This make it simple to generate good guesses for the roots of equation (90). We use Newton's method to polish the zeros to a high degree of accuracy. Table 5 shows the different dispersion (decay-rate) relations for different boundary conditions and modes.

B Matrices for Hydrodynamic Benchmark
In our formulation of this problem, we use the statevector The linear operators M and L are Thus, the part of the problem treated implicitly is identical to the linear diffusion equation described in section 4.3. The explicit terms are For boundary conditions, we either replace the last rows of the three components of u in the L and M matrices with the r = 1 operator, and then replace the corresponding entries of F (X) with the appropriate components of u 0 (α BC = 2), or we impose the boundary conditions with τ corrections (α BC = 0). We fix the = 0 component of all fields to be zero.

C Matrices for Convection Benchmark
We use the statevector The explicit terms are For the boundary conditions, we replace the bottom row of the three velocity components and temperature blocks with the impenetrable, stress-free, and fixed temperature conditions (see section 4.3) when we use α BC = 2, and implement the boundary conditions using τ corrections for α BC = 0. For the = 0 mode, we set the pressure and velocities to zero, but evolve the temperature equation normally.