1. Introduction
Rotating convection plays an important dynamical role in stars and planets, where it is believed to be one of the primary drivers of global scale magnetic fields (Busse Reference Busse1975; Glatzmaier & Roberts Reference Glatzmaier and Roberts1995; Kageyama & Sato Reference Kageyama and Sato1995; Stanley & Glatzmaier Reference Stanley and Glatzmaier2010; Jones Reference Jones2011; Aurnou et al. Reference Aurnou, Calkins, Cheng, Julien, King, Nieves, Soderlund and Stellmach2015), and possibly gives rise to coherent large-scale flows such as zonal jets and vortices, as observed on the giant planets (Heimpel et al. Reference Heimpel, Yadav, Featherstone and Aurnou2022; Siegelman et al. Reference Siegelman2022; Böning et al. Reference Böning, Wulff, Dietrich, Wicht and Christensen2023). Understanding the physics of turbulence driven by rotating convection remains challenging due to the vast range of spatiotemporal scales. As a result, the parameter space accessible to direct numerical simulation (DNS) and laboratory experiments remains distant from that which characterises natural convective systems. An approach that is often taken to overcome this limitation is to identify asymptotic scaling behaviour in model output so that extrapolation to the conditions of natural systems is possible (e.g. Christensen Reference Christensen2002; Aurnou Reference Aurnou2007). The development of asymptotically reduced equation sets is another strategy that has been particularly helpful for improving understanding of rotating convective turbulence and dynamos in simplified planar geometries since more extreme parameter regimes can be accessed in comparison to DNS (Sprague et al. Reference Sprague, Julien, Knobloch and Werne2006; Calkins, Julien & Marti Reference Calkins, Julien and Marti2013; Calkins et al. Reference Calkins, Julien, Tobias and Aurnou2015; Yan & Calkins Reference Yan and Calkins2022). In general, excellent agreement has been found between the results of plane layer asymptotic models and DNS when an overlap in parameter space is possible (Plumley et al. Reference Plumley, Julien, Marti and Stellmach2016; Yan & Calkins Reference Yan and Calkins2022). However, asymptotic models have not yet been developed for global spherical geometries. Towards this end, the present investigation utilises DNS in spherical shell geometries to determine the asymptotic scaling behaviour of various key quantities, including force balances, length scales, and convective and zonal flow speeds.
We consider Boussinesq convection in a rotating spherical shell with angular velocity 
$\varOmega$. The geometry is specified by the aspect ratio 
$\eta = r_i/r_o$, where 
$r_i$ is the radius of the inner sphere and 
$r_o$ is the radius of the outer sphere. The fluid is forced via a temperature contrast 
$\Delta T$ between the inner and outer boundaries, and gravity varies linearly with radius. For this system, the convection dynamics are determined by the sizes of several non-dimensional parameters, including the Rayleigh number and Ekman number, respectively defined by
\begin{equation} Ra = \frac{g_o \alpha \Delta T H^3}{\kappa \nu}, \quad Ek = \frac{\nu}{\varOmega H^2}, \end{equation}
where 
$g_o$ is the gravitational acceleration at the outer boundary, 
$\alpha$ is the coefficient of thermal expansion, 
$H = r_o - r_i$ is the depth of the fluid layer, 
$\kappa$ is the thermal diffusivity and 
$\nu$ is the kinematic viscosity. The most unstable state consists of convective Rossby waves that align with the rotation axis and drift in the prograde direction (Busse Reference Busse1970). Asymptotic theory, valid in the limit 
$Ek \rightarrow 0$, has shown that the critical Rayleigh number scales as 
$Ra_c = O(Ek^{-4/3})$, and the critical azimuthal (zonal) wavenumber scales as 
$m_c = O(Ek^{-1/3})$ (Roberts Reference Roberts1968; Busse Reference Busse1970; Jones, Soward & Mussa Reference Jones, Soward and Mussa2000; Dormy et al. Reference Dormy, Soward, Jones, Jault and Cardin2004). Thus, large Rayleigh numbers are needed to drive rotating convection and the subsequent motions become increasingly smaller scale as the Ekman number is reduced.
As the Rayleigh number is increased beyond critical, strong zonal flows develop in rotating spherical convection provided stress-free mechanical boundary conditions are applied on the inner and outer spherical surfaces (e.g. Gilman Reference Gilman1978; Aurnou & Olson Reference Aurnou and Olson2001; Christensen Reference Christensen2002). These zonal flows are characterised by alternating regions of prograde and retrograde motion that are approximately invariant in the direction of the rotation axis. For a fixed value of 
$Ek$, the number of jets is controlled both by the Rayleigh number and the shell aspect ratio, 
$\eta$. In full sphere and small aspect ratio geometries (
$\eta \lesssim 0.6$), simulations typically find a single prograde jet in the equatorial region and retrograde jets at higher latitudes (Aurnou & Olson Reference Aurnou and Olson2001; Christensen Reference Christensen2002; Lin & Jackson Reference Lin and Jackson2021). As 
$\eta$ and 
$Ra$ are increased, there is a tendency for multiple jets to form at higher latitudes, leading to a banded structure that is reminiscent of the flows observed on the gas giant planets (Christensen Reference Christensen2001; Heimpel, Aurnou & Wicht Reference Heimpel, Aurnou and Wicht2005; Jones & Kuzanyan Reference Jones and Kuzanyan2009; Gastine, Heimpel & Wicht Reference Gastine, Heimpel and Wicht2014; Heimpel et al. Reference Heimpel, Yadav, Featherstone and Aurnou2022). The number of zonal jets that appear can be related to the Rhines length scale (Rhines Reference Rhines1975; Heimpel et al. Reference Heimpel, Aurnou and Wicht2005; Gastine et al. Reference Gastine, Heimpel and Wicht2014). For sufficiently large Rayleigh numbers, there is an eventual transition to a retrograde equatorial jet and prograde high-latitude jets (Aurnou, Heimpel & Wicht Reference Aurnou, Heimpel and Wicht2007; Soderlund Reference Soderlund2019).
Steady zonal flows are driven by Reynolds stresses and damped by global scale viscous stresses. Thus, the scaling behaviour of the zonal flow is intrinsically linked to the scaling of the correlations of the convective flows (Christensen Reference Christensen2002). Such correlations are not known a priori, though various scaling theories have been presented in the literature. For perfectly correlated small-scale velocity components, as relevant near the onset of convection, the zonal flow amplitude exhibits a quadratic dependence on the small-scale velocity (Aubert et al. Reference Aubert, Brito, Nataf, Cardin and Masson2001). However, Christensen (Reference Christensen2002) found that correlations decrease with increasing supercriticality, which causes the quadratic scaling of the zonal flow to eventually break down. Later investigations have found that this loss of correlation in the small-scale velocity is a monotonically decreasing function of 
$Ra$ (Showman, Kaspi & Flierl Reference Showman, Kaspi and Flierl2011; Gastine & Wicht Reference Gastine and Wicht2012). For sufficiently large Rayleigh numbers, Christensen (Reference Christensen2002) and Lin & Jackson (Reference Lin and Jackson2021) find evidence that the zonal flow scaling approaches a ‘diffusion-free’ regime, so-called because of the lack of dependence on 
$\nu$ and 
$\kappa$, though within this regime, the dynamics is no longer geostrophic since inertia and the Coriolis force are of the same order of magnitude.
Aside from their correlations, it is also important to understand the scaling of the convective flow speeds themselves. One scaling theory that is often invoked to explain convective flow speed scaling behaviour in rotating convection is the so-called Coriolis-inertia-Archimedean (CIA) balance (e.g. Cardin & Olson Reference Cardin and Olson1994; Jones Reference Jones2015; Aurnou, Horn & Julien Reference Aurnou, Horn and Julien2020). This balance predicts that the dominant convective length scale behaves like 
$\ell \sim \widetilde {\textit {{Ra}}}^{1/2} Ek^{1/3}$ and the global scale Reynolds number should scale like 
$Re = U H/\nu \sim Ek Ra/Pr$, where 
$U$ is a characteristic flow speed, the reduced Rayleigh number is defined as 
$\widetilde {\textit {Ra}} = Ra Ek^{4/3}$ and the thermal Prandtl number is 
$Pr = \nu /\kappa$. King & Buffett (Reference King and Buffett2013) analysed length scales and flow speeds in a broad suite of numerical dynamo simulations in spherical geometries and concluded that viscous effects remained important in controlling these quantities. CIA theory is often contrasted with the scaling of the linear instability scale, 
$\ell \sim Ek^{1/3}$. However, it is important to emphasise that, from the point of view of asymptotics, both the linear ‘viscous’ length scale and the length scale predicted by CIA theory are of the size 
$O(Ek^{1/3})$ given that 
$\widetilde {\textit {Ra}}$ is an order unity asymptotic parameter.
Several previous investigations have tested these CIA scaling predictions in spherical geometries with no-slip boundary conditions (Gastine, Wicht & Aubert Reference Gastine, Wicht and Aubert2016; Guervilly, Cardin & Schaeffer Reference Guervilly, Cardin and Schaeffer2019; Long et al. Reference Long, Mound, Davies and Tobias2020), laboratory experiments in rotating cylindrical geometries (Madonia et al. Reference Madonia, Aguirre Guzmán, Clercx and Kunnen2021), as well as asymptotic models of plane layer convection with stress-free boundary conditions (Maffei et al. Reference Maffei, Krouss, Julien and Calkins2021; Oliver et al. Reference Oliver, Jacobi, Julien and Calkins2023). Gastine et al. (Reference Gastine, Wicht and Aubert2016) carried out a comprehensive survey of rotating spherical convection and found that the length scale for their smallest Ekman number cases (
$Ek=3\times 10^{-7}$) approached the Rhines scaling predicted by the CIA balance, and that the interior dissipation could also be approximated through a CIA balance. However, as we discuss in the present study, Gastine et al. (Reference Gastine, Wicht and Aubert2016) did not investigate the convection and zonal flow separately. A similar approach was taken by Long et al. (Reference Long, Mound, Davies and Tobias2020) in which constant heat flux thermal boundary conditions were used. Guervilly et al. (Reference Guervilly, Cardin and Schaeffer2019) simulated rotating spherical convection with no-slip boundary conditions and 
$Pr=0.01$ using both a two-dimensional quasi-geostrophic model and three-dimensional DNS, and found that the length scale increases with Rayleigh number for all parameter ranges used; they find evidence of length scales and flow speeds approaching the CIA scaling predictions as the Ekman number is reduced.
The rotating convection experiments of Madonia et al. (Reference Madonia, Aguirre Guzmán, Clercx and Kunnen2021) exhibited an increase of the integral length scale with Rayleigh number, but at a slower rate than that predicted by CIA theory. The experiments of Hawkins et al. (Reference Hawkins, Cheng, Abbate, Pilegard, Stellmach, Julien and Aurnou2023) and Abbate & Aurnou (Reference Abbate and Aurnou2023) suggest that the turbulent length scale remains comparable to the viscous length scale over the parameter space investigated. In asymptotic models, Maffei et al. (Reference Maffei, Krouss, Julien and Calkins2021) found flow speed scaling behaviour consistent with CIA theory over a finite range of Rayleigh numbers; the deviation from CIA theory at large Rayleigh numbers was attributed to the effects of the large scale vortex (LSV) that is generated in this system. Oliver et al. (Reference Oliver, Jacobi, Julien and Calkins2023) also found that certain measures of the convective length scales show an increase with 
$\widetilde {\textit {Ra}}$, but at a rate that is slower than the exponent of 
$1/2$. Importantly, however, the asymptotic models do not find a CIA force balance in the fluid interior. Instead, the buoyancy force is only comparable to the Coriolis and inertial forces within the thermal boundary layers, though viscous effects are equally important in these regions of the flow domain. Moreover, the ratio of the viscous force to the buoyancy force was found to be an increasing function of 
$\widetilde {\textit {Ra}}$, indicating that the CIA balance is never achieved in plane layer convection (Maffei et al. Reference Maffei, Krouss, Julien and Calkins2021; Oliver et al. Reference Oliver, Jacobi, Julien and Calkins2023). The conclusion from these asymptotic studies is that convective length scales remain viscously controlled. This same conclusion was reached by Yan & Calkins (Reference Yan and Calkins2022) who found similar behaviour in DNS of rapidly rotating convection driven dynamos in the plane layer geometry. An additional important finding in the studies of Madonia et al. (Reference Madonia, Aguirre Guzmán, Clercx and Kunnen2021), Yan & Calkins (Reference Yan and Calkins2022), Oliver et al. (Reference Oliver, Jacobi, Julien and Calkins2023) is that the viscous dissipation length scale remains approximately constant with increasing 
$\widetilde {\textit {Ra}}$ – this indicates that length scale evolution in rotating convective turbulence is fundamentally different than non-rotating convective turbulence where the dissipation length decreases strongly with increasing Rayleigh number (e.g. Yan, Tobias & Calkins Reference Yan, Tobias and Calkins2021). We observe similar behaviour for the length scales and force balances in the spherical simulations reported in the present investigation.
The scaling behaviour of key quantities such as flow speeds and length scales is linked to the force balances in the governing equations. To our knowledge, the asymptotic scaling behaviour of the force balances (and terms in the heat equation) in spherical convection simulations have not been reported to date, though several previous dynamo studies have computed these forces over a range of parameters. For an Ekman number of 
$Ek=10^{-4}$ in a spherical dynamo, Soderlund, King & Aurnou (Reference Soderlund, King and Aurnou2012) noted that the Lorentz force was smaller than the inertial term and that the Lorentz force did not significantly alter the convective flow as compared with the purely hydrodynamic model. However, Soderlund et al. (Reference Soderlund, King and Aurnou2012) also noted that at smaller Ekman numbers, the Lorentz force seemed to make larger changes to the convection. In general, simulations find that the force balance in the mean equations is thermal wind to leading order, with the Lorentz force entering at higher order in the zonal component of the mean momentum equation (Aubert Reference Aubert2005; Calkins, Orvedahl & Featherstone Reference Calkins, Orvedahl and Featherstone2021; Orvedahl, Featherstone & Calkins Reference Orvedahl, Featherstone and Calkins2021). For the small-scale convective dynamics, dynamo studies find that the zeroth-order force balance is geostrophic, the first-order force balance is between the ageostrophic Coriolis force, the buoyancy force and the Lorentz force, and inertial and viscous forces enter at the next order (Yadav et al. Reference Yadav, Gastine, Christensen, Wolk and Poppenhaeger2016) (note, however, these authors did not separate the mean and fluctuating dynamics). Here, we find a similar sequence of balances in the fluctuating momentum equation, though there does not appear to be an asymptotic difference between the buoyancy force, the viscous force and the inertial force, which is similar to plane layer rotating convection. Also in agreement with previous asymptotic studies, we find that the ratio of the viscous force to the buoyancy force is an increasing function of 
$\widetilde {\textit {Ra}}$. For the large-scale dynamics, we find that the flows are geostrophically balanced to leading order.
In this paper, we investigate the asymptotic behaviour of rapidly rotating convection and the associated zonal flows in a spherical shell with stress-free boundary conditions. Knowledge of such scaling behaviour is crucial for the development of asymptotic models. Overall, we find excellent agreement between the asymptotic predictions and the results of the nonlinear simulations. We develop a prediction for the asymptotic scaling of the amplitude of the zonal flow and conclude that the zonal flow must remain dependent on viscosity given that viscous stresses are the sole saturation mechanism for this component of the flow. The paper is organised as follows. In § 2, we describe the model and governing equations. We give a brief overview of the asymptotic theory in § 3 and numerical results are analysed in § 4. A discussion is provided in § 5.
2. Model
The governing equations consist of the conservation laws for momentum, thermal energy and mass. We non-dimensionalise these equations using the length 
$H$, the large-scale viscous diffusion time 
$H^2/\nu$ and the temperature scale 
$\Delta T$. With this non-dimensionalisation, the governing equations are given by
$$\begin{gather} {\partial_t} {\boldsymbol{u}} + {\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}} =-\frac{2}{Ek} \boldsymbol{\hat{z}} \times {\boldsymbol{u}} - \frac{1}{Ek}\boldsymbol{\nabla} P + \frac{Ra}{Pr} \left ( \frac{r}{r_o} \right ) T \, \hat{{\boldsymbol{r}}}+ \nabla^2 {\boldsymbol{u}}, \end{gather}$$
$$\begin{gather}{\partial_t} T + {\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} T = \frac{1}{Pr} \nabla^2T, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{u}} = 0, \end{gather}$$
where 
${\boldsymbol {u}}=\langle u_r, u_\theta, u_\phi \rangle$ is the fluid velocity, 
$T$ is the temperature, 
$P$ is the pressure and 
$r$ is radius. The ‘axial’ direction points in the direction of the rotation axis, and the axial and radial unit vectors are denoted by 
$\boldsymbol {\hat {z}}$ and 
$\boldsymbol {\hat {r}}$, respectively. In all of the simulations presented here, we fix 
$Pr = 1$.
The boundary conditions are impenetrable (
$u_r=0$), stress-free and fixed temperature. We use the code Rayleigh to numerically solve the governing equations (Featherstone et al. Reference Featherstone, Edelmann, Gassmoeller, Matilsky, Orvedahl and Wilson2022). Rayleigh is a pseudo-spectral code which uses spherical harmonics to represent data on spherical shells and Chebyshev polynomials to represent data in the radial direction. A 
$2/3$ de-aliasing is used for both the spherical harmonics and the Chebyshev polynomials. Rayleigh has been tested against the benchmark cases from Christensen et al. (Reference Christensen2001). We choose our spherical resolutions such that the kinetic energy spectra exhibits a minimum of four orders of magnitude in separation between the peak value and the value at the largest spherical harmonic degree. A similar method is used for resolving the radial direction where we compute the magnitudes of the corresponding Chebyshev coefficients. Moreover, the output from our simulations shows excellent agreement with previous studies (e.g. Christensen Reference Christensen2002). Time stepping is carried out using a second-order semi-implicit Crank–Nicolson method for the linear terms and a second-order Adams–Bashforth method for the nonlinear terms. The time step is chosen adaptively to maintain numerical stability.
We opt for stress-free mechanical boundary conditions in the present investigation since one of the primary purposes is to determine how the amplitude of zonal flows depends on the Ekman number and Rayleigh number. As shown in many previous investigations, zonal flows tend to be strongly damped when no-slip boundary conditions are used at the parameter values accessible in most simulations. Stress-free boundary conditions are also more applicable to natural systems such as stars and giant planets.
As discussed by Jones et al. (Reference Jones, Boronski, Brun, Glatzmaier, Gastine, Miesch and Wicht2011), the use of stress-free boundary conditions implies that angular momentum is a conserved quantity, yet numerical simulations can exhibit a spurious growth of this quantity due to the intrinsic error associated with numerical time-stepping discretisations. To deal with this problem, we monitor the angular momentum that builds up in the simulations. We found that this build up can lead to a relative error in our Ekman number of up to 
$0.1\,\%$, though it is typically many orders of magnitude smaller than this. To determine whether this build up has an influence on our output quantities, we performed a test at one particular parameter combination (
$Ek=10^{-5}, Ra=2\times 10^{8}$) that showed a higher growth of angular momentum than most other cases. One case used a strict conservation of angular momentum setting which changes the boundary condition of the 
$l=1$ spherical harmonic mode on the outer boundary to enforce zero total angular momentum in the simulation for all time (as discussed by Jones et al. Reference Jones, Boronski, Brun, Glatzmaier, Gastine, Miesch and Wicht2011). We ran this case from the same initial state that had a buildup of angular momentum both for strict conservation of angular momentum and with stress-free boundaries. Between these two runs, we found that the mean Reynolds number differed by less than one percent and the fluctuating Reynolds number differed by less than three percent. Some of the observed differences are likely the result of time averaging errors due to the large time variability associated with the presence of relaxation oscillations. Given the small change in values, we believe the amount of growth in the angular momentum is not affecting our results.
2.1. Notation and outputs
Due to the symmetry of the model set-up around the rotation axis, it is convenient to define mean and fluctuating components of some scalar quantity 
$X$ relative to an azimuthal or zonal average, i.e.
\begin{equation} \bar{X} = \int_0^{2{\rm \pi}} X \,{\rm d}\phi, \end{equation}
with 
$X'=X-\bar {X}$.
We define the Reynolds number as
\begin{equation} Re = {\overline{\sqrt{\langle {\boldsymbol{u}} \boldsymbol{\cdot} {\boldsymbol{u}} \rangle}},}^t \end{equation}
where 
$\langle {\cdot } \rangle$ denotes a volume average and 
${\overline {({\cdot })}}^t$ denotes a time average. We further define the mean and fluctuating (convective) Reynolds numbers respectively as
\begin{equation} Re_z = {\overline{\sqrt{\langle \bar{\boldsymbol{u}}\phantom{\prime^{\prime}}\! \boldsymbol{\cdot} \bar{\boldsymbol{u}} \rangle}}}^t, \quad Re_c = {\overline{\sqrt{\langle {\boldsymbol{u}}' \boldsymbol{\cdot} {\boldsymbol{u}}' \rangle}}}^t. \end{equation}In all of the simulations, we find that the zonal component of the mean flow dominates and we therefore refer to the mean Reynolds number as the ‘zonal’ Reynolds number.
We analyse several different length scales in this paper. Christensen & Aubert (Reference Christensen and Aubert2006) (see also Gastine et al. Reference Gastine, Wicht and Aubert2016; Long et al. Reference Long, Mound, Davies and Tobias2020) define the spherical harmonic length scale as
\begin{equation}
\left(\ell_{sh}\right)^{-1} =
{\overline{\left(\dfrac{\displaystyle \sum_{l=0}^{l_{max}} \sum_{m=0}^{l}
l\mathcal{E}_l^m}{\displaystyle \pi \sum_{l=0}^{l_{max}} \sum_{m=0}^{l}
\mathcal{E}_l^m}\right)}}^t,
\end{equation}
where 
$l_{max}$ is the maximum spherical harmonic degree in the simulation and 
$\mathcal {E}_l^m$ is the radially averaged kinetic energy density of spherical harmonic degree 
$l$ and order 
$m$. Thus, the volume averaged kinetic energy density is given by
\begin{equation} \frac{1}{2} \langle {\boldsymbol{u}} \boldsymbol{\cdot} {\boldsymbol{u}} \rangle = \sum_{l=0}^{l_{max}} \sum_{m=0}^{l} \mathcal{E}_l^m. \end{equation}Due to the strong influence of the zonal flow in our simulations, we separate the spherical harmonic length scale into a zonal and a non-zonal component. We define these length scales as
\begin{equation}
\left(\ell_{sh}'\right)^{-1}
= {\overline{\left(\dfrac{\displaystyle \sum_{l=1}^{l_{max}} \sum_{m=1}^{l}
l\mathcal{E}_l^m}{\displaystyle\pi \sum_{l=1}^{l_{max}} \sum_{m=1}^{l}
\mathcal{E}_l^m}\right)}}^t, \quad
\left(\bar{\ell}_{sh}\right)^{-1} =
{\overline{\left(\dfrac{\displaystyle\sum_{l=0}^{l_{max}}l\mathcal{E}_l^{m=0}}{\displaystyle\pi
\sum_{l=0}^{l_{max}} \mathcal{E}_l^{m=0}}\right)}}^t.
\end{equation}We define the fluctuating and mean Taylor microscales as
\begin{equation}
\ell'_{tm} =
{\overline{\sqrt{\dfrac{\langle {\boldsymbol{u}}'
\boldsymbol{\cdot} {\boldsymbol{u}}' \rangle }{\langle
(\boldsymbol{\nabla} \times {\boldsymbol{u}}')
\boldsymbol{\cdot} (\boldsymbol{\nabla} \times
{\boldsymbol{u}}') \rangle}}}}^t, \quad \bar{\ell}_{tm} =
{\overline{\sqrt{\dfrac{\vphantom{\frac{1^3}{1^2}}\langle \bar{\boldsymbol{u}}
\boldsymbol{\cdot} \bar{\boldsymbol{u}} \rangle }{\langle
(\boldsymbol{\nabla} \times \bar{\boldsymbol{u}})
\boldsymbol{\cdot}(\boldsymbol{\nabla} \times
\bar{\boldsymbol{u}}) \rangle}}}}^t,
\end{equation}respectively. The Taylor microscale can be considered a viscous dissipation length scale since it characterises the length scale at which viscous effects become important.
The Nusselt number is calculated according to
\begin{equation} Nu = \frac{ \partial_r \bar{T}^{\theta,\phi,t} \bigr\rvert_{r=r_{o}}}{\partial_r T_c\bigr\rvert_{r=r_{o}}}, \end{equation}
where 
${\overline {({\cdot })}}^{\theta,\phi,t}$ is a shell and time average, and 
$T_c$ is the conductive temperature profile, which satisfies
\begin{equation} \nabla^2 T_c = 0, \quad T_c(r_i) = 1, \quad T_c(r_o) = 0. \end{equation}We also define the viscous dissipation rates of the mean and fluctuating velocity fields according to
\begin{equation} \bar{\varepsilon} = {\overline{\langle \left(\boldsymbol{\nabla} \times \bar{{\boldsymbol{u}}} \right)^2 \rangle}}^t \quad \hbox{and} \quad \varepsilon' = {\overline{\langle \left( \boldsymbol{\nabla} \times {\boldsymbol{u}}' \right)^2 \rangle}}^t, \end{equation}respectively.
3. Theory
Here we provide arguments for the scaling behaviour of various quantities. We use the term asymptotic to mean rotationally constrained motions in which the Ekman number and Rossby number, 
$Ro = U/ ( 2 \varOmega H )$, are both small relative to unity, i.e. 
$(Ro, Ek) \ll 1$. The dimensional flow speed 
$U$ characterises the magnitude of the fluctuating velocity field. We expect many aspects of the asymptotic theory for linear rotating convection in spherical geometries presented by Jones et al. (Reference Jones, Soward and Mussa2000) and Dormy et al. (Reference Dormy, Soward, Jones, Jault and Cardin2004) to hold here, though some of these scalings must be modified to account for nonlinear terms in the governing equations. In particular, the scaling of the convective flow speeds and temperature perturbation need to be reduced by a factor of 
$Ek^{1/3}$ relative to those of Jones et al. (Reference Jones, Soward and Mussa2000) and Dormy et al. (Reference Dormy, Soward, Jones, Jault and Cardin2004), though this difference does not influence the leading order force balance in the fluctuating momentum equation when zonal flows are weak. The scaling of the large-scale zonal flow depends on the small-scale convective velocity, so we first consider theoretical scaling laws for the convective velocity and the corresponding convective length scale. Such scaling laws can be found by examining the fluctuating momentum equation and the fluctuating heat equation, which are respectively given by
\begin{gather} \begin{aligned}&\left[{\partial_t}
{\boldsymbol{u}}^{\prime} + \bar{\boldsymbol{u}}
\boldsymbol{\cdot}\boldsymbol{\nabla}
{\boldsymbol{u}}^{\prime}\right] +
{\boldsymbol{u}}^{\prime} \boldsymbol{\cdot}
\boldsymbol{\nabla} \bar{\boldsymbol{u}} +
{\boldsymbol{u}}^{\prime}
\boldsymbol{\cdot}\boldsymbol{\nabla}
{\boldsymbol{u}}^{\prime}-
\overline{{\boldsymbol{u}}^{\prime}
\boldsymbol{\cdot}\boldsymbol{\nabla}
{\boldsymbol{u}}^{\prime}}\cr &\qquad =-\frac{2}{Ek}
\boldsymbol{\hat{z}} \times {\boldsymbol{u}}^{\prime} -
\frac{1}{Ek}\boldsymbol{\nabla} P' + \frac{Ra}{Pr} \left (
\frac{r}{r_o} \right ) T' \, \hat{\boldsymbol{r}}+ \nabla^2
{\boldsymbol{u}}^{\prime}, \end{aligned}
\end{gather}
\begin{gather} \left[{\partial_t} T' + \bar{\boldsymbol{u}} \boldsymbol{\cdot} \nabla T'\right] + {\boldsymbol{u}}^{\prime} \boldsymbol{\cdot} \boldsymbol{\nabla} \bar{T} + {\boldsymbol{u}}^{\prime} \boldsymbol{\cdot} \boldsymbol{\nabla} T' - \overline{ {\boldsymbol{u}}^{\prime} \boldsymbol{\cdot} \boldsymbol{\nabla} T' } = \frac{1}{Pr} \nabla^2T'. \end{gather}
We find that the terms in brackets can be large compared to some terms due to the large amplitude of the zonal flow. However, summing the terms in brackets leads to results smaller than the individual terms in the brackets, which physically means that the Lagrangian time derivative is smaller than the Eulerian time derivative. We will therefore consider the sum of the terms in brackets rather than each individually. It is well known from linear theory that the Coriolis force and pressure gradient force are dominant terms in the limit 
$( Ro, Ek ) \rightarrow 0$. We will see below that the zonal flow can also modify the leading order force balance. To study the first-order effects, we eliminate the pressure gradient by taking the curl of the fluctuating momentum equation. This operation yields
\begin{align} &\boldsymbol{\nabla} \times\left[{\partial_t} {\boldsymbol{u}}^{\prime} + \bar{\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}}^{\prime}\right] + \boldsymbol{\nabla} \times ({\boldsymbol{u}}^{\prime} \boldsymbol{\cdot} \boldsymbol{\nabla} \bar{\boldsymbol{u}}) + \boldsymbol{\nabla}\times ({\boldsymbol{u}}^{\prime} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}}^{\prime} - \overline{{\boldsymbol{u}}^{\prime} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}}^{\prime}}) \nonumber\\ &\quad =-\frac{2}{Ek} \frac{\partial {\boldsymbol{u}}^{\prime}}{\partial z} +\frac{Ra}{Pr} \boldsymbol{\nabla} \times \left ( \frac{r}{r_o} T' \, \hat{\boldsymbol{r}} \right ) + \nabla^2 \left ( \boldsymbol{\nabla} \times {\boldsymbol{u}}^{\prime} \right ) . \end{align}
Assuming length scales of order one in the 
$z$-direction, length scales of order one for azimuthally averaged terms and length scales of order 
$\ell$ otherwise, the vorticity equation and heat equation can be approximately written as
$$\begin{gather} \boldsymbol{\nabla} \times\left[{\partial_t} {\boldsymbol{u}}^{\prime} + \bar{\boldsymbol{u}} \boldsymbol{\cdot}\boldsymbol{\nabla} {\boldsymbol{u}}^{\prime}\right] + \frac{\bar{\boldsymbol{u}} {\boldsymbol{u}}^{\prime}}{\ell} + \frac{\left({\boldsymbol{u}}^{\prime}\right)^2}{\ell^2} - \boldsymbol{\nabla}\times \overline{{\boldsymbol{u}}^{\prime} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}}^{\prime}} \sim-\frac{{\boldsymbol{u}}^{\prime}}{Ek} + \frac{Ra}{Pr}\frac{T'}{\ell} +\frac{{\boldsymbol{u}}^{\prime}}{\ell^3}, \end{gather}$$
$$\begin{gather}\left[{\partial_t} T' + \bar{\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} T'\right] + {\boldsymbol{u}}^{\prime} \bar{T} + \frac{{\boldsymbol{u}}^{\prime} T'}{\ell} - \overline{ {\boldsymbol{u}}^{\prime} \boldsymbol{\cdot} \boldsymbol{\nabla} T' } \sim \frac{1}{Pr}\frac{T'}{\ell^2}, \end{gather}$$
where factors of order one have been dropped. We now follow Aurnou et al. (Reference Aurnou, Horn and Julien2020) and assume a balance between 
${\boldsymbol {u}}^{\prime } \bar {T}$ and 
${\boldsymbol {u}}^{\prime } T'\ell ^{-1}$ in the temperature equation. Using 
$\bar {T}=O(1)$ yields 
$T' \sim \ell$. Plugging this relation in for 
$T'$ in the momentum equation and assuming a CIA balance where the fluctuating–fluctuating advection term is used yields
\begin{equation} \frac{\left({\boldsymbol{u}}^{\prime}\right)^2}{\ell^2} \sim-\frac{{\boldsymbol{u}}^{\prime}}{Ek} \sim \frac{Ra}{Pr}, \end{equation}which can be solved for the convective velocity and length scale to give
\begin{equation} {\boldsymbol{u}}^{\prime} \sim \frac{Ra Ek}{Pr}, \quad \ell \sim \sqrt{\frac{Ra Ek^2}{Pr}}. \end{equation}We can rewrite these expressions in terms of the reduced Rayleigh number as
\begin{equation} {\boldsymbol{u}}^{\prime} \sim Ek^{-1/3}\frac{\widetilde{\textit{Ra}}}{Pr}, \quad \ell \sim Ek^{1/3} \sqrt{\frac{\widetilde{\textit{Ra}}}{Pr}}. \end{equation}
Here again we note that both the CIA and viscous length scale have the same 
$Ek^{1/3}$ dependence. Thus, for a fixed value of 
$\widetilde {\textit {Ra}}$, the CIA length scale follows the same scaling as predicted for the viscous length scale. Indeed, one of the points that we stress in the present study is that all length scales in this system are viscously selected to leading order, with order one variations away from this viscous scale since 
$\widetilde {\textit {Ra}} = O(1)$ (e.g. Yan & Calkins Reference Yan and Calkins2022; Oliver et al. Reference Oliver, Jacobi, Julien and Calkins2023). We note that the statement 
$\widetilde {\textit {Ra}} = O(1)$ is in reference to the scaling of this quantity with respect to the Ekman number only, and it does not require that the reduced Rayleigh number is small, though the simulations presented here have a limited range.
We can understand why these different balances produce the same Ekman dependence by making a different set of assumptions than the assumptions used for a CIA balance. The first assumption we make is that the Ekman dependence for any term can be written in terms of a power law. We will assume that the Ekman dependence of the convective flow speeds, the length scale and the fluctuating temperature can be written as 
${\boldsymbol {u}}^{\prime } = O(Ek^{x_u})$, 
$\ell = O(Ek^{x_\ell })$ and 
$T' = O(Ek^{x_T})$, respectively. The second assumption is that ratios of certain terms in the momentum and heat equations do not change when the Ekman number is reduced. We will assume that the advection of fluctuating velocity by fluctuating velocity, viscosity and buoyancy follow the same Ekman number scaling in the fluctuating momentum equation. In the heat equation, we will assume that conduction and advection of the mean temperature by the fluctuating velocity follow the same Ekman number scaling. These assumptions produce the system of equations given by
\begin{equation} 2x_u-x_\ell = x_u-2x_\ell =-4/3 + x_T, \quad x_T-2x_\ell = x_u. \end{equation}
Note that we have assumed the mean temperature and the length scale of the mean temperature do not depend on the Ekman number. Solving this system of equations yields 
$x_u=-1/3$, 
$x_\ell = 1/3$ and 
$x_T = 1/3$, which implies that 
${\boldsymbol {u}}^{\prime } = O(Ek^{-1/3})$, 
$\ell = O(Ek^{1/3})$ and 
$T' = O(Ek^{1/3})$. Note that the Ekman number dependence derived here is the same as the Ekman number dependence derived using the CIA balance written in terms of the reduced Rayleigh number. Therefore, we can get the 
$Ek^{1/3}$ scaling for the length scale by assuming that various terms follow the same Ekman number scaling without actually assuming any balances a priori. These scalings are equivalent to the scalings used to derive the asymptotic model of rotating convection in a Cartesian geometry (e.g. Sprague et al. Reference Sprague, Julien, Knobloch and Werne2006).
An asymptotic constraint on the amplitude of the zonal flow can now be obtained upon examination of the mean momentum equation. The mean momentum equation is given by
\begin{equation} {\partial_t} \bar{\boldsymbol{u}} + \bar{\boldsymbol{u}} \boldsymbol{\cdot}\boldsymbol{\nabla} \bar{\boldsymbol{u}} + \overline{ {\boldsymbol{u}}^{\prime} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}}^{\prime} } + \frac{2}{Ek} \, \boldsymbol{\hat{z}} \times \bar{\boldsymbol{u}} = -\frac{1}{Ek}\boldsymbol{\nabla} \bar{P} + \frac{Ra}{Pr} \left ( \frac{r}{r_o} \right ) \bar{T} \, \hat{{\boldsymbol{r}}} + \nabla^2 \bar{\boldsymbol{u}} . \end{equation}The zonal component is then
\begin{equation} {\partial_t} \bar{u}_{\phi}+ \left [ \bar{\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \bar{\boldsymbol{u}} \right ] _{\phi} + \left [ \overline{ {\boldsymbol{u}}^{\prime} \boldsymbol{\cdot}\boldsymbol{\nabla} {\boldsymbol{u}}^{\prime} } \right ] _{\phi} + \frac{2}{Ek} \, \left [ \boldsymbol{\hat{z}} \times \bar{\boldsymbol{u}} \right ] _{\phi} = \left [ \nabla^2 \bar{\boldsymbol{u}} \right ] _{\phi} , \end{equation}
where the square brackets and corresponding subscript are used for brevity. The above equation allows for a straightforward interpretation of the zonal flow dynamics. Time dependence and advection of the zonal flow by mean meridional flows are captured by the first two terms on the left-hand side; for simplicity, we find it useful to refer to the latter quasi-linear term with the descriptor ‘mean-mean’. The term 
$[\overline { {\boldsymbol {u}}^{\prime } \boldsymbol {\cdot } \boldsymbol {\nabla } {\boldsymbol {u}}^{\prime } }]_{\phi }$ includes the divergence of the Reynolds stresses and acts as the primary source of the zonal flow. The zonal component of the mean Coriolis force can be written as 
$[\boldsymbol {\hat {z}} \times \bar {\boldsymbol {u}}]_{\phi } = \cos \theta \, \bar {u}_{\theta } + \sin \theta \, \bar {u}_r$ such that only meridional circulation appears in this term.
Averaging (3.11) along the axial direction and in time leaves only the two advective terms and the viscous term. For a careful derivation of the balance between the averaged advection term and the viscous term in spherical coordinates, see Dietrich, Gastine & Wicht (Reference Dietrich, Gastine and Wicht2017). Previous studies have found that the mean-mean term is small (e.g. Dietrich et al. Reference Dietrich, Gastine and Wicht2017), and we also find this to be the case in our simulations. While this advection term may be important for transporting angular momentum (McIntyre Reference McIntyre1998; Miesch & Hindman Reference Miesch and Hindman2011), it cannot play a role in setting the amplitude of the time-averaged (steady) zonal flow. Thus, we anticipate that there must be a balance between the fluctuating advection term and the viscous term when the flow is rapidly rotating. Letting 
$\overline {({\cdot })}^{\theta,z,t}$ denote an average over 
$\phi$, 
$z$ and time, the preceding argument implies that
\begin{equation} \left [ \overline{ {\boldsymbol{u}}^{\prime} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}}^{\prime} }^{\theta,z,t} \right ] _{\phi} \approx \left [ \nabla^2 \bar{{\boldsymbol{u}}}^{\theta,z,t} \right ] _{\phi} . \end{equation}
The above balance is expected to hold for all values of 
$Ek$ and 
$Ra$; the zonal flow is therefore intrinsically dependent on viscosity and we should not expect its scaling behaviour to be ‘diffusion-free’. We note that a similar balance holds for mean flows in planar geometries, allowing for constraints on the amplitude of analogous ‘zonal flows’ (e.g. Nicoski, Yan & Calkins Reference Nicoski, Yan and Calkins2022). Finally, since averaged quantities vary on order one length scales, this balance suggests
\begin{equation} \bar{\boldsymbol{u}}_\phi = O ( C_R u' u' ) , \end{equation}
where 
$C_R$ represents the correlation of the fluctuating velocity components. We might expect that this correlation gets weaker as the reduced Rayleigh number is increased and the flow becomes less constrained by rotation; this weakening was confirmed via direct computation by Christensen (Reference Christensen2002). However, if we restrict our analysis to the regime of asymptotically small Ekman number, then the dynamics should depend on 
$\widetilde {\textit {Ra}}$, rather than on 
$Ek$ and 
$Ra$ separately; in this sense, we do not expect 
$C_R$ to depend on the Ekman number and we can then use (3.13) to determine how the zonal flow scales with 
$Ek$. The asymptotic analysis predicts that the convective velocity scales as 
$\boldsymbol {u}' = O(Ek^{-1/3})$, which implies that
\begin{equation} \bar{\boldsymbol{u}}_\phi = O \left ( Ek^{-2/3} \right ) , \end{equation}thus indicating that the zonal flow is intrinsically dependent on viscosity, albeit in an asymptotic sense. For time dependence in the zonal flow, we require that the time derivative is comparable to the viscous force so that
\begin{equation} {\partial_t} \bar{\boldsymbol{u}}_\phi = O \left ( Ek^{-2/3} \right ) , \end{equation}
which, along with (3.14), indicates that the zonal flow time scale is 
$O(1)$ in our non-dimensional, large-scale viscous diffusion units. Thus, the zonal flow varies on a large-scale viscous diffusion time. This property is one of the reasons that computations of rotating spherical convection with stress free boundary conditions are so demanding – very long integration is necessary to saturate the amplitude of the zonal flow and reach a statistically stationary state.
The asymptotic constraint on the zonal flow amplitude allows for additional insight into the force balance by which it is constrained. In the limit 
$Ek \rightarrow 0$, the Rayleigh number must scale as 
$Ra = O (Ek^{-4/3} )$ to generate convection (Roberts Reference Roberts1968). Along with the fact that the magnitude of the mean temperature is independent of the Ekman number, this indicates that the mean buoyancy force scales as 
$O(Ek^{-4/3} )$. The radial and co-latitudinal components of the mean Coriolis force both contain 
$\bar {u}_{\phi }$, thus indicating that these components scale as 
$O ( Ek^{-5/3})$, which is larger than the mean buoyancy force by a factor of 
$O (Ek^{-1/3} )$. Thus, the zonal flow is geostrophically balanced to leading order, i.e.
\begin{equation} -2 \sin \theta \, \bar{u}_{\phi} \approx- \partial_{r} \bar{P}, \quad - 2 \cos \theta \, \bar{u}_{\phi} \approx- r^{-1} \partial_{\theta} \bar{P} . \end{equation}It is informative to compare the scaling behaviour of zonal flows that are geostrophically balanced with a zonal flow that is in thermal wind balance (e.g. Calkins et al. Reference Calkins, Orvedahl and Featherstone2021). An order of magnitude estimate for the scaling of the thermal wind component of the zonal flow can be obtained if we balance the mean Coriolis force with the mean buoyancy force,
\begin{equation} - \frac{2 \sin \theta }{Ek} \, \bar{u}_{\phi}^{tw} \approx \frac{Ra}{Pr} \left ( \frac{r}{r_o} \right ) \bar{T} \quad \Rightarrow \quad \bar{u}_{\phi}^{tw} \sim \frac{Ra Ek}{Pr} . \end{equation}Using the definition of the reduced Rayleigh number, this becomes
\begin{equation} \bar{u}_{\phi}^{tw} \sim \frac{\widetilde{\textit{Ra}}}{Pr} Ek^{-1/3}. \end{equation}
Additionally, since 
$\widetilde {\textit {Ra}} = O(1)$, this implies 
$\bar {u}_{\phi }^{tw} = O(Ek^{-1/3})$, which is of the same order as the convective flow speeds. Thus, zonal flows that satisfy a thermal wind balance are substantially weaker than those that are geostrophic. Interestingly, the above scaling also provides an estimate for the scaling of the thermal wind with Rayleigh number and represents a ‘diffusion-free’ scaling in the sense that it indicates that the thermal wind does not depend on either 
$\nu$ or 
$\kappa$. The thermal wind scaling seems to be loosely consistent with the zonal flows present in the dynamos of Calkins et al. (Reference Calkins, Orvedahl and Featherstone2021); when a magnetic field is present, the Lorentz force strongly damps the geostrophic component of the zonal flow.
4. Numerical results
4.1. Overview
To test the theoretical arguments given in the previous section, we perform a suite of direct numerical simulations of convection across a range of Ekman number and Rayleigh number. We expect the agreement between simulation output and the asymptotic trends to become better as the Ekman and Rossby numbers are made smaller. We consider the aspect ratios 
$\eta =0.35$ and 
$\eta =0.7$. The smallest Ekman numbers that were simulated are 
$Ek=10^{-6}$ and 
$Ek=3\times 10^{-5}$ for the 
$\eta =0.35$ and 
$\eta =0.7$ aspect ratios, respectively.
Details of the simulations are contained in tables 1–2 in the Appendix. The Rossby number, 
$Ro=EkRe_c$, is often used to determine the rotational constraint. We note that for all of our cases, the Rossby number is smaller than 
$0.1$, which would suggest that our simulations are rotationally constrained, and this is confirmed directly via a force balance analysis. Another common measure of rotational constraint in convection is the convective Rossby number, 
$Ro_c = \sqrt {RaEk^2/Pr}$. The largest convective Rossby number in our simulations is 
$Ro_c=0.94$ for the largest Rayleigh number case with 
$Ek=3\times 10^{-4}$. However, we note that for our lower Ekman number cases, the convective Rossby number tends to be smaller. For example, our largest convective Rossby number at 
$Ek=3\times 10^{-6}$ is 
$Ro_c=0.11$.
Our simulations share some overlap in parameter space with several previous investigations. In particular, Christensen (Reference Christensen2002) ran cases similar to our 
$\eta =0.35$ cases to Ekman numbers as low as 
$Ek=10^{-5}$. We reach Ekman numbers as low as 
$Ek=10^{-6}$ at that aspect ratio, although the highest Rayleigh numbers we reach are approximately a factor of three smaller than Christensen (Reference Christensen2002) reached because we are mainly interested in the rapidly rotating regime. As an example, we reach 
$Ra/Ra_c \sim 30$ for 
$Ek=10^{-5}$, while Christensen (Reference Christensen2002) reached 
$Ra/Ra_c \sim 105$ for this same Ekman number.
All of the simulations presented in this study generate zonal flows with an amplitude that is at least comparable to the amplitude of the underlying convection, though in many of the cases, the zonal flow is substantially larger in magnitude compared to the convection. Although we fix the Prandtl number to be unity, smaller Prandtl numbers tend to yield stronger zonal flows (e.g. Aubert et al. Reference Aubert, Brito, Nataf, Cardin and Masson2001). The qualitative nature of the zonal flows was similar for much of the parameter space covered, and representative cases for both shell thicknesses are shown in figure 1. The zonal flows we observe in this study are similar to zonal flows observed in previous works. These zonal flows are characterised by a nearly invariant structure in the axial (
$z$) direction, and consist of a single prograde jet at the equator with retrograde jets at higher latitudes. Three of the simulations that were performed with a thin shell (
$\eta =0.7$) developed high-latitude jets; one example is shown in figure 1(c). As discussed in previous work (e.g. Heimpel et al. Reference Heimpel, Aurnou and Wicht2005), the number of jets is related to the Rhines length scale; in general, high-latitude jets are more likely to form when the shell is made thinner, the Ekman number is made smaller and the Rayleigh number is larger. A subset of our cases exhibit relaxation oscillations in which the convection mainly occurs during short bursts; this behaviour was also observed in previous work (e.g. Christensen Reference Christensen2002). Movie 1 from the supplementary material available at https://doi.org/10.1017/jfm.2024.78 shows the radial velocity in the equatorial plane over the course of one relaxation oscillation with the convective Reynolds number shown for reference. From this movie, we see that during times of weak convection, the convection is strongest near the inner boundary in the equatorial plane. However, during times of strong convection, the convection fills the whole region in the equatorial plane.
Figure 1. Instantaneous visualisations of the zonal flow where red indicates prograde motion and blue indicates retrograde motion: (a) 
$\eta =0.35$, 
$Ek=10^{-5}$, 
$Ra=2\times 10^8$ (
$\widetilde {\textit {Ra}} \approx 43$); (b) 
$\eta =0.7$, 
$Ek=3\times 10^{-5}$, 
$Ra=3.2\times 10^7$ (
$\widetilde {\textit {Ra}} \approx 30$); (c) 
$\eta =0.7$, 
$Ek=3\times 10^{-5}$, 
$Ra=6.4\times 10^{7}$ (
$\widetilde {\textit {Ra}} \approx 60$).
We find that using 
$\widetilde {\textit {Ra}}$, as opposed to the supercriticality measure, 
$Ra/Ra_c$, results in improved collapse of our data when comparing with asymptotic predictions; Christensen (Reference Christensen2002) also found that using 
$\widetilde {\textit {Ra}}$ improved the collapse of some data. This effect likely arises from the slow rate of convergence of the critical Rayleigh number to the predicted asymptotic scaling of 
$Ra_c \sim Ek^{-4/3}$ (e.g. Dormy et al. Reference Dormy, Soward, Jones, Jault and Cardin2004; Barik et al. Reference Barik, Triana, Calkins, Stanley and Aurnou2023).
4.2. Flow speeds
Global root mean square (r.m.s.) values of both the fluctuating and mean velocity are computed to determine their scaling behaviour with respect to the Ekman number. Note that the mean velocity is dominated by the zonal flow (i.e. the 
$\phi$-component of the mean velocity). Figure 2(a) shows the fluctuating Reynolds number for both the 
$\eta =0.35$ cases and the 
$\eta =0.7$ cases. Figure 2(b) shows the asymptotically rescaled fluctuating Reynolds number, i.e. 
$\widetilde {Re}_c = Ek^{1/3} Re_c$, for the two different aspect ratios. We find that the rescaled data are order unity and collapse onto a single curve, which supports the 
${\boldsymbol {u}}^{\prime } = O(Ek^{-1/3})$ asymptotic scaling for the fluctuating velocity. However, we note that the Ekman number scaling of the convective Reynolds number might be time dependent. If we calculated the convective Reynolds number using data only during the convective peaks of the relaxation oscillations, we would obtain a steeper scaling closer to 
$Ek^{-1/2}$. This suggests that the time series for the convective velocity becomes more strongly peaked at lower Ekman number. Note that we expect deviation from this asymptotic scaling behaviour as the system loses rotational constraint; this deviation is particularly noticeable in figure 2(b) for the high-Rayleigh-number regime (
$\widetilde {\textit {Ra}} \gtrsim 100$) for the two largest Ekman numbers used in the 
$\eta =0.35$ simulations, 
$Ek=3\times 10^{-4}$ and 
$Ek=10^{-4}$. Figure 2(c,d) shows two versions of the compensated convective Reynolds number: 
$\widetilde {\textit {Re}}_c \widetilde {\textit {Ra}}{\vphantom {Ra}}^{-1}$ and 
$\widetilde {\textit {Re}}_c \widetilde {\textit {Ra}}{\vphantom {Ra}}^{-3/2}$. We observe in figure 2(c) that the compensated Reynolds number 
$\widetilde {\textit {Re}}_c \widetilde {\textit {Ra}}{\vphantom {Ra}}^{-1}$ becomes nearly horizontal for our large Rayleigh number cases at large Ekman number and small aspect ratio, which suggests these cases may be scaling as 
$Re_c \sim \widetilde {\textit {Ra}}$. However, this scaling behaviour may be localised in 
$\widetilde {\textit {Ra}}$ space. For sufficiently small Ekman number and large Rayleigh number, the compensated plot for 
$\widetilde {\textit {Re}}_c \widetilde {\textit {Ra}}{\vphantom {Ra}}^{-3/2}$ collapses the data well, though the scaling appears slightly weaker than 
$\widetilde {\textit {{Ra}}}^{3/2}$ which suggests that the convective Reynolds number scales approximately as 
$Re_c \sim \widetilde {\textit {Ra}}{\vphantom {Ra}}^{3/2}$ in this regime.
Figure 2. Reynolds number characterising the flow speeds of the fluctuating (convective) velocity versus the reduced Rayleigh number, 
$\widetilde {\textit {Ra}}$: (a) the convective Reynolds number 
$Re_c$; (b) the rescaled convective Reynolds number 
$\widetilde {\textit {Re}}_c = Ek^{1/3} Re_c$; (c) the compensated convective Reynolds number 
$\widetilde {\textit {Re}}_c \widetilde {\textit {Ra}}{\vphantom {Ra}}^{-1}$; (d) the compensated convective Reynolds number 
$\widetilde {\textit {Re}}_c \widetilde {\textit {Ra}}{\vphantom {Ra}}^{-3/2}$. The filled symbols represent 
$\eta =0.35$ cases and the hollow symbols represent 
$\eta =0.7$ cases.
Figure 3(a) shows the mean Reynolds number as a function of 
$\widetilde {\textit {Ra}}$, and figure 3(b) shows the corresponding asymptotically rescaled mean Reynolds number, 
$\widetilde {Re}_z = Ek^{2/3} Re_z$. As mentioned previously, the mean flow is dominated by the zonal component in all of our simulations. While there is clearly some spread in the rescaled data for the thick shell cases, there is an indication that the data collapse to an asymptotic state as 
$Ek \rightarrow 0$. Moreover, the rescaled values are order unity. There appears to be better collapse for the thin shell cases, indicating that the fluid depth may play an important role. Also note the excellent collapse for the three thin shell data points near 
$\widetilde {\textit {Ra}} \approx 60$ – the two lower Ekman number cases of these three develop prograde high-latitude jets as shown in figure 1(c). Taken together, the data seem to support that the zonal flow scales as 
$Ek^{-2/3}$ in the rapidly rotating regime. We note that Gastine et al. (Reference Gastine, Heimpel and Wicht2014) found an empirical scaling for the zonal Rossby number of 
$Ro_{zon} \sim Ra^{0.6} Ek^{0.99}$, which, when converted to a Reynolds number and using 
$Ra \sim Ek^{-4/3}$, gives 
$Re_z \sim Ek^{-0.81}$. This result broadly agrees with the scaling of 
$Re_z \sim Ek^{-2/3}$ derived in this paper.
Figure 3. Reynolds number characterising the flow speeds of the mean (zonal) velocity field versus 
$\widetilde {\textit {Ra}}$: (a) the zonal Reynolds number 
$Re_z$; (b) the rescaled zonal Reynolds number 
$\widetilde {Re}_z=Ek^{2/3} Re_z$. The solid line is a least squares fit of the data for the 
$\eta =0.35$, 
$Ek=10^{-5}$ cases to a line and is given by 
$Re_z =(0.022\widetilde {\textit {Ra}}-0.19)Ek^{-2/3}$; the dashed line is a least squares fit of all the 
$\eta =0.7$ data to a line and is given by 
$Re_z=(0.031\widetilde {\textit {Ra}}-0.016)Ek^{-2/3}$. The symbols are the same as defined in figure 2.
While the balance between viscosity and Reynolds stresses predicts a zonal flow scaling of 
$Ek^{-2/3}$, this balance is unable to predict how the zonal flow scales with the reduced Rayleigh number due to the fact that the correlation of the fluctuating velocity components is an unknown function of 
$\widetilde {\textit {Ra}}$. Thus, we make an empirical fit of 
$Re_z$ with respect to the reduced Rayleigh number, which is shown in figure 3(b). We find that a line does a good job of fitting our small Ekman number cases. However, our larger Ekman number cases at 
$\eta =0.35$ are run to larger values of the reduced Rayleigh number and are not well fit by a line at larger values of 
$\widetilde {\textit {Ra}}$. One possibility for this effect is that these larger Ekman number cases at large Rayleigh number are no longer in the rapidly rotating regime and therefore follow a different trend. Another possibility is that the affine dependence we have used to fit the data only holds near the onset of convection and that the zonal Reynolds number follows a different trend for large enough values of the reduced Rayleigh number, even as the Ekman number is decreased.
It is also interesting to consider the relationship between the zonal and convective Reynolds numbers. Equation (3.13) suggests that the relevant quantity to consider is 
$Re_z/Re_c^2$, which we show in figure 4. We predicted that this ratio is independent of the Ekman number, which we see is a good approximation for our range of parameters. Christensen (Reference Christensen2002) noted that the relation 
$Re_z \sim Re_c^2$ holds only when the correlation between the fluctuating velocity components is constant, which occurs near the onset of convection. Anelastic simulations of rotating spherical convection also find that 
$Re_z/Re_c^2$ is nearly constant near the onset of convection and that 
$Re_z/Re_c^2$ decreases with increasing 
$Ra$ for sufficiently large Rayleigh number (Gastine & Wicht Reference Gastine and Wicht2012). Our simulations are also consistent with this behaviour, showing that 
$Re_z/Re_c^2$ is approximately constant up to 
$\widetilde {\textit {Ra}} \approx 15$. For 
$\widetilde {\textit {Ra}} \gtrsim 15$, we observe that 
$Re_z/Re_c^2$ decreases with increasing 
$\widetilde {\textit {Ra}}$. This trend of 
$Re_z/Re_c^2$ decreasing as 
$\widetilde {\textit {Ra}}$ is increased does not show a systematic dependence on the Ekman number, since all data points appear to collapse well. Because 
$Re_z/Re_c^2$ is related to the correlation of the convective velocity components, this suggests that the correlation of the convective velocity components decreases as 
$\widetilde {\textit {Ra}}$ is increased, even at asymptotically small Ekman numbers.
Figure 4. Ratio of the zonal Reynolds number to the square of the convective Reynolds number. The least squares fit to a power law scaling calculated using the 
$\eta =0.35$, 
$Ek=10^{-4}$ cases is shown as a black line, and is given by 
$C_R = 14.21\widetilde {\textit {{Ra}}}^{-1.526}$. The symbols are the same as defined in figure 2.
4.3. Length scales
In this section, we compute length scales for both the mean flow and the small-scale convection. To provide a physical picture of how these length scales vary with Ekman number, we show snapshots of the radial velocity in the equatorial plane for three different Ekman number cases with 
$\widetilde {\textit {Ra}} \approx 50$ in figure 5. As expected, the typical length scale of the convection decreases with decreasing Ekman number.
Figure 5. Visualisations of the radial velocity for three thick shell cases during a time interval of strong convection. Equatorial slices are shown in panels (a–c) and meridional slices are shown in panels (d–f) for the cases (a,d) 
$Ek=3\times 10^{-4}$, 
$Ra=2.5\times 10^6$ (
$\widetilde {\textit {Ra}} \sim 50.2$); (b,e) 
$Ek=3\times 10^{-5}$, 
$Ra=5.4\times 10^7$ (
$\widetilde {\textit {Ra}} \sim 50.3$); (c,f) 
$Ek=3\times 10^{-6}$, 
$Ra=1.2\times 10^{9}$ (
$\widetilde {\textit {Ra}} \sim 51.9$).
One way to quantitatively study how the length scale varies with Ekman number is to consider how the kinetic energy power spectrum varies with Ekman number. We define the sum of the kinetic energy power spectrum over the order 
$m$ and the normalisation of the kinetic energy power spectrum summed over the order 
$m$ as
\begin{equation}
E_l =
{\overline{\displaystyle\sum\limits_{m=1}^{l} \mathcal{E}_l^m}^t}, \quad
\tilde{E}_l = {\overline{\left(\frac{\displaystyle\sum_{m=1}^{l}
\mathcal{E}_l^m}{\displaystyle\sum_{l=1}^{l_{max}} \sum_{m=1}^{l}
\mathcal{E}_l^m}\right)}^t}.
\end{equation}
Due to the influence of the zonal flow, we again only sum over 
$m\geqslant 1$ modes to remove the influence of the zonal flow. Figure 6 shows how these kinetic energy spectra vary with Ekman number. The critical azimuthal wavenumber 
$m_c$ is shown for reference, and 
$m_c$ takes the values 
$6$, 
$8$, 
$11$, 
$16$ and 
$23$ for the cases of 
$Ek=3\times 10^{-4}$, 
$Ek=10^{-4}$, 
$Ek=3\times 10^{-5}$, 
$Ek=10^{-5}$ and 
$Ek=3\times 10^{-6}$, respectively. These values were calculated using the eigenvalue solver Kore (Barik et al. Reference Barik, Triana, Calkins, Stanley and Aurnou2023). We see that the smaller Ekman number cases have more power in each mode and extend to higher values of 
$l$ in comparison to larger Ekman number cases. We attempt to collapse these data by multiplying the degree 
$l$ by 
$Ek^{1/3}$, the expected scaling for the length scale, and multiplying the amplitude of the data by 
$Ek^{1/3}$. We choose the amplitude scaling such that the area under the kinetic energy power spectrum, the volume averaged kinetic energy density, scales as 
$Ek^{-2/3}$, which is consistent with the scaling of our convective Reynolds number. This collapsed data for the kinetic energy power spectrum are shown in figure 6(b). We see that the collapse of the data is reasonable for values of 
$l$ greater than 
$m_c$, but not as good for values of 
$l$ less than 
$m_c$. This discrepancy might suggest that large length scales are following a different scaling than anticipated. We investigate the behaviour of the kinetic energy power spectrum for the small 
$l$ modes in figure 6(c,d). The power spectrum for these small 
$l$ modes appears to shift more slowly to the right with Ekman number than the large 
$l$ modes. Several papers on the linear onset of convection predict a longer length scale in the cylindrical radial direction than in the azimuthal direction. Dormy et al. (Reference Dormy, Soward, Jones, Jault and Cardin2004) predicts an 
$Ek^{2/9}$ length scale for cases with differential heating, and an 
$Ek^{1/6}$ length scale for cases with certain boundary conditions and heating conditions. We show the collapse of our data for the predicted 
$Ek^{2/9}$ length scale in figure 6(d), which does a reasonable job of collapsing our data. However, we do not have a large enough range in Ekman number to discern between an 
$Ek^{1/6}$ and an 
$Ek^{2/9}$ scaling.
Figure 6. Time- and radially averaged spherical harmonic kinetic energy spectra. Cases with 
$\eta =0.35$ and similar values of the reduced Rayleigh number are shown. Circles denote the critical azimuthal wavenumber (
$m_c$) at the onset of convection for each parameter set and squares denote the degree 
$l$ derived from the spherical harmonic length scale 
$\ell _{sh}^{\prime } = {\rm \pi}/l$. (a) Spherical harmonic spectra. (b) Rescaled spherical harmonic spectra. (c) Spherical harmonic spectra at small 
$l$. (d) Rescaled spherical harmonic spectra at small 
$l$.
Figure 7 shows the kinetic energy spectra for the 
$\eta =0.35$, 
$Ek=10^{-5}$ cases averaged in time and radius. Figure 7(b) shows the spherical harmonic spectra normalised to have area one. We find that the kinetic energy contained in high 
$l$ values increases slightly with increasing 
$\widetilde {\textit {Ra}}$, but that the overall shape of the spectra does not change significantly. This behaviour indicates that small-scale convection becomes more prominent as the Rayleigh number is increased. The spectra show that many different length scales are present within the flow, and these scales may exhibit different scaling behaviour with varying Ekman number. Thus, different methods for calculating length scales may lead to different conclusions. Here we use a few different methods to hopefully capture the scaling for both large and small scales.
Figure 7. Time- and radially averaged spherical harmonic spectra. Select cases from 
$\eta =0.35$, 
$Ek=10^{-5}$ are shown. Circles denote the critical azimuthal wavenumber at the onset of convection and squares denote the spherical harmonic degree calculated from 
$\ell _{sh}^{\prime }$, as in figure 6. (a) Spherical harmonic spectra. (b) Spherical harmonic spectra normalised to have an area of one.
The computed length scales are all shown in figure 8. The spherical harmonic length scale is shown in panels (a,b), and the Taylor microscale is shown in panels (c,d); the asymptotically rescaled lengths are shown in panels (b,d). The spherical harmonic length scale is observed to decrease rapidly for 
$\widetilde {\textit {Ra}} \lesssim 50$, then levels off for larger values of 
$\widetilde {\textit {Ra}}$; not surprisingly, this behaviour is consistent with the trends observed in the kinetic energy spectra. There may be a trend suggesting that this length scale converges to a nearly constant value at the largest values of 
$\widetilde {\textit {Ra}}$, though this behaviour may be occurring outside of the rapidly rotating regime. Both the decreasing and constant dependence of the length scale on Rayleigh number observed here is in contrast to the prediction made by the CIA balance where the length scale increases as the Rayleigh number increases. Figure 8(b) shows the asymptotically rescaled spherical harmonic length scale. As with the flow speeds, this rescaled quantity is order unity and we find significantly less scatter in the data when viewed in this rescaled coordinate, suggesting that this particular length scale does approximately scale as 
$Ek^{1/3}$. However, we also note that a subtle Ekman number dependence still exists in the rescaled length scale; larger Ekman number cases seem to level off at smaller rescaled length scales in comparison to the smaller Ekman number cases. This behaviour might indicate that the system is either converging slowly towards the 
$O(Ek^{1/3})$ length scale with decreasing Ekman number or that the asymptotic dependence of this length scale is slightly weaker than 
$Ek^{1/3}$. Our data for the kinetic energy power spectrum indicate that a second length scale may be present at small 
$l$ values which follows a different Ekman number scaling. Thus, the spherical harmonic length scale may have a composite asymptotic dependence on 
$Ek$ since it may be capturing both the 
$Ek^{2/9}$ and 
$Ek^{1/3}$ asymptotic scalings.
Figure 8. Convective length scales for the 
$\eta =0.35$ cases: (a,b) spherical harmonic length scale; (c,d) Taylor microscale. Panel (b,d) shows the asymptotically rescaled length scales. The symbols are the same as defined in figure 2.
Figure 8(c,d) shows the Taylor microscale, where we see that it follows a very similar trend in comparison to the spherical harmonic length scale. The collapse of the rescaled quantity indicates that the viscous dissipation length scale is consistent with an 
$Ek^{1/3}$ dependence across the entire range of parameters used. The Taylor microscale appears to decrease with 
$\widetilde {\textit {Ra}}$ over a wider range of 
$\widetilde {\textit {Ra}}$ than the spherical harmonic length scale, which becomes approximately constant with 
$\widetilde {\textit {Ra}}$. If we interpret the spherical harmonic length scale as the energy containing length scale (similar to an integral length scale), then our computed length scales indicate that the scale separation present in these simulations is relatively small across the entire range of parameters since 
$\ell _{sh}^{\prime }$ and 
$\ell _{tm}^{\prime }$ are not too dissimilar in value. For instance, at 
$\widetilde {\textit {Ra}} \approx 50$, 
$Ek^{-1/3} \ell_{sh}^{\prime} \approx 4$ and 
$Ek^{-1/3} \ell_{tm}^{\prime} \approx 0.9$. We note that the trends observed in both length scales contrast sharply with computed length scales in non-rotating Rayleigh–Bénard convection, where both the energy containing length scale and the dissipation length scale decrease rapidly with increasing Rayleigh number (e.g. Yan et al. Reference Yan, Tobias and Calkins2021).
To study the behaviour of the 
$Ek^{2/9}$ length scale found in the kinetic energy spectra, we define a length scale based on the degree 
$l_{peak}$ where the kinetic energy spectrum (again with the 
$m=0$ mode removed) reaches a maximum. The peak length scale is then given by 
$\ell ^{\prime }_{peak}={\rm \pi} /l_{peak}$. To smooth the resulting data set, we first fit each spectrum with a high-order polynomial; we found that 15th degree and higher worked well. Figure 9(a) shows the peak length scale and figure 9(c) shows the peak length scale rescaled by 
$Ek^{-2/9}$. We note that while the peak length scale approximately scales as 
$Ek^{2/9}$, we obtain a somewhat better fit for 
$Ek^{1/6}$ as shown in figure 9(d). Regardless, this peak length scale seems to scale differently than 
$\ell ^{\prime }_{sh}$, which suggests that 
$\ell ^{\prime }_{peak}$ does in fact correspond to a different length scale than the small-scale convective length scale. We also note that the length scale 
$\ell ^{\prime }_{peak}$ increases with increasing Rayleigh number, in contrast to 
$\ell ^{\prime }_{sh}$ and 
$\ell ^{\prime }_{tm}$, which either decrease or remain constant as the Rayleigh number is increased. Figure 9(b) shows the peak length scale as a function of the Rossby number characterising the convective flow speeds, 
$Ro=EkRe_c$, where it can be seen that 
$\ell ^{\prime }_{peak}$ approximately scales as 
$Ro^{1/4}$, which is in contrast to the 
$Ro^{1/2}$ scaling that is predicted from the CIA balance. To ensure the validity of our results with polynomial interpolation, we also estimated the peak length scale by replacing 
$\mathcal {E}_l^m$ in (2.9a,b) with 
$(\mathcal {E}_l^m)^{10}$, which weights the length scale more towards the peak. Using this estimate of the peak length scale, we found qualitatively similar results to what is shown in figure 9.
Figure 9. Length scale calculated from the degree 
$l_{peak}$ where the fluctuating kinetic energy power spectrum peaks: (a) peak length scale; (b) peak length scale versus the Rossby number characterising the convective flow speeds (two Rossby number scalings are shown for reference); (c) peak length scale rescaled by 
$Ek^{-2/9}$; (d) peak length scale rescaled by 
$Ek^{-1/6}$. The symbols are as defined in figure 2.
We now consider length scales of the zonal flow, which are shown in figure 10. Unlike the convective length scales, we expect these mean length scales to be order unity and approximately independent of the Ekman number; the data set confirms these expectations. The mean length scale reaches a maximum at approximately 
$\widetilde {\textit {Ra}} \sim 25$ and, for larger values of the reduced Rayleigh number, decays slowly over the range of Rayleigh numbers used in this study. We also computed a Taylor microscale for the mean flow, as shown in figure 10(b); although this length scale is slightly smaller than the spherical harmonic length scale, it shows a similar behaviour with 
$Ek$ and 
$\widetilde {\textit {Ra}}$ in comparison with the spherical harmonic length scale.
Figure 10. Length scales computed from the mean flow: (a) spherical harmonic length scale; (b) Taylor microscale. The symbols are the same as defined in figure 2.
4.4. Heat equation balances
The asymptotic scaling behaviour of terms in both the momentum and heat equations are linked. For this purpose, we study the scaling of various terms in both equations beginning with the heat equation in the present section. Figure 11 shows the r.m.s. of the fluctuating temperature for all of the thick shell cases. We observe a systematic decrease in the magnitude of the fluctuating temperature as the Ekman number is reduced. Figure 11(b) shows that the fluctuating temperature is well described by the relationship 
$T' = O(Ek^{1/3})$ for fixed 
$\widetilde {\textit {Ra}}$, as predicted by asymptotic theory. A line gives a decent approximation for how the fluctuating temperature varies with 
$\widetilde {\textit {Ra}}$, in contrast to the CIA balance, which predicts a reduced Rayleigh number dependence of 
$T'\sim \widetilde {\textit {Ra}}{\vphantom {Ra}}^{1/2}$. Figure 11(c) shows the compensated fluctuating temperature using the reduced Rayleigh number dependence predicted by the CIA balance: 
$\widetilde {\textit {Ra}}{\vphantom {Ra}}^{1/2}$. We see that the data are not very well described by this Rayleigh number dependence; for none of the low-Ekman-number cases does there exist a large range of 
$\widetilde {\textit {Ra}}$ where the compensated data are horizontal. Some of the larger Ekman number cases appear to be fit by this relation for 
$60 \lesssim \widetilde {\textit {Ra}} \lesssim 85$. However, this is a rather narrow range, so it is difficult to say whether these cases are actually following the predicted CIA scaling. We also note that for the 
$Ek=3\times 10^{-4}$ cases, the largest 
$\widetilde {\textit {Ra}}$ cases exhibit a slight decrease in the r.m.s. value of 
$T'$ as compared to the lower 
$\widetilde {\textit {Ra}}$ cases, so 
$T'$ is not strictly increasing with 
$\widetilde {\textit {Ra}}$ in our data.
Figure 11. Time-averaged volume r.m.s. of the fluctuating temperature; (a) raw data; (b) asymptotically rescaled data; (c) compensated data. The line 
$Ek^{-1/3}T'=0.01\widetilde {\textit {Ra}}+0.5$ is shown for reference in panel (b). The symbols are the same as defined in figure 2.
Figure 12 shows how the various terms from the fluctuating heat equation depend on Ekman number and reduced Rayleigh number. Figure 12(b,d, f) shows how well the asymptotic prediction of the scaling collapses the data. The diffusion term and the advection of the mean temperature by the fluctuating velocity term are well collapsed by the predicted 
$Ek^{-1/3}$ scaling, although the advection of the fluctuating temperature by the fluctuating velocity is not as well collapsed by the asymptotic scaling. The data suggest that a slightly stronger dependence on the Ekman number is present since the rescaled data from the smaller Ekman number cases lie above that of the larger Ekman number cases. However, we note that the magnitude of the rescaled terms are all comparable, thus providing support for the asymptotic theory.
Figure 12. Time-averaged volume r.m.s. values of various terms from the fluctuating heat equation. The terms shown are: (a) 
$\boldsymbol {u}' \boldsymbol {\cdot } \boldsymbol {\nabla } T'$; (c) 
$\boldsymbol {u}' \boldsymbol {\cdot } \boldsymbol {\nabla } \bar {T}$; (e) 
$\nabla ^2 T'$. Asymptotically rescaled data: (b) 
$Ek^{1/3}\boldsymbol {u}' \boldsymbol {\cdot } \boldsymbol {\nabla } T'$; (d) 
$Ek^{1/3}\boldsymbol {u}' \boldsymbol {\cdot }\boldsymbol {\nabla } \bar {T}$; (f) 
$Ek^{1/3}\nabla ^2 T'$. The symbols are the same as defined in figure 2.
Figure 13 shows how the various terms from the fluctuating heat equation vary with the reduced Rayleigh number for two fixed values of the Ekman number. The large magnitude of the time derivative term is a consequence of the zonal flow; advection by the zonal flow is not balanced and causes large accelerations. Therefore, the sum of the advection by the zonal flow and the time derivative is smaller by comparison and this sum is approximately balanced by the term 
$\boldsymbol {u}' \boldsymbol {\cdot } \boldsymbol {\nabla } T'$. Another interesting observation from figure 13 is that there is a near balance of the diffusion term and the advection of the mean temperature term for a wide range of 
$\widetilde {\textit {Ra}}$. This balance was also noted by Maffei et al. (Reference Maffei, Krouss, Julien and Calkins2021).
Figure 13. Time-averaged volume r.m.s. of several terms from the fluctuating heat equation averaged in time for thick shell simulations: (a) 
$Ek=3\times 10^{-4}$; (b) 
$Ek=10^{-5}$.
4.5. Force balances
In this section, we investigate the scaling behaviour of the forces present in the radial component of the fluctuating momentum equation (similar scaling was observed in the other two components). Figure 14 shows the time average and r.m.s. over the entire domain of the viscous force, buoyancy force and the fluctuating advective term. For the fluctuating advective term, we remove the spherically symmetric 
$l=0$ mode as this mode is not dynamically relevant. From the arguments in the theory section, we expect these three terms to all scale as 
$Ek^{-1}$ in our non-dimensional units. We test these scalings by multiplying the data in figure 14(a,c,e) by 
$Ek$, which is shown in figure 14(b,d,f). We see that the predicted Ekman number scalings are consistent with the data. Similar to the heat equation terms in the previous section, the collapse for the advective term is not as good as the collapse for the other terms. This difference either indicates that the advective term is converging slower than the other terms or that the advective term is converging to a slightly different scaling than predicted by asymptotic theory. We also note that the scaling of the advective term is time-dependent; removing the intervals of time in which strong convective bursts are occurring in the time series produces a scaling for the fluctuating advective term that more closely follows the predicted 
$Ek^{-1}$ scaling. This effect suggests that the relaxation oscillations are playing some role in the scaling, which may be due to the larger Rossby numbers that occur during these times. The Coriolis force is not shown since its scaling follows from the scaling of the fluctuating velocity given in figure 2 – in our units, it scales as 
$Ek^{-4/3}$.
Figure 14. Time-averaged volume r.m.s. of various terms from the radial component of the fluctuating momentum equation: (a) viscous force; (c) buoyancy force; (e) fluctuating advection term; (b) rescaled viscous force; (d) rescaled buoyancy force; (f) rescaled fluctuating advection term. The symbols are the same as defined in figure 2.
These scalings suggest that for a fixed value of the reduced Rayleigh number, the ratio of the viscous force, buoyancy force and fluctuating advection term remain approximately fixed as the Ekman number is decreased. However, the ratio of either the viscous force, the buoyancy force or the fluctuating advective term to the Coriolis force will scale as 
$Ek^{1/3}$. While we have only shown the radial component of the fluctuating force balance, we note that the 
$\theta$ and 
$\phi$ components of the fluctuating force balance show similar behaviour. We also tested removing the thermal boundary layers, and did not observe a qualitative change in any of the trends shown.
Some terms from the momentum equation follow a stronger Ekman number scaling due to the presence of the zonal flow, and these terms are shown in figure 15. The zonal flow strongly advects the fluctuating velocity, which is unbalanced and causes large accelerations in the small-scale fluid structures. Because the zonal flow scales as 
$Ek^{-2/3}$, the fluctuating velocity scales as 
$Ek^{-1/3}$ and the length scale of the fluctuating velocity scales as 
$Ek^{1/3}$, we expect the advection by the zonal flow to scale as 
$Ek^{-4/3}$, which we also expect for the scaling of the time derivative. Figure 15(b,d) shows the collapse of the mean advection term and the time derivative term for this scaling. We see that both the advection by the zonal flow and the time derivative follow an 
$Ek^{-4/3}$ scaling, which we note is the same scaling followed by the Coriolis force. Therefore, the advection by the zonal flow and the time derivative appear at leading order asymptotically, and the sum of these two terms is somewhat smaller than either individually, as shown in figure 15(c,e) and in figure 16.
Figure 15. Time-averaged volume r.m.s. of various terms from the radial component of the fluctuating momentum equation: (a) advection by the zonal flow; (c) time derivative (inertia); (e) the sum of the time derivative and advection by the zonal flow; (b) rescaled advection by the zonal flow; (d) rescaled time derivative; (f) rescaled sum of the time derivative and advection by the zonal flow. The symbols are the same as defined in figure 2.
Figure 16. Volume r.m.s. of several terms from the fluctuating radial momentum equation averaged in time for thick shell simulations: (a) 
$Ek=3\times 10^{-4}$; (b) 
$Ek=10^{-5}$.
Figure 16 shows how the different terms in the fluctuating radial momentum equation depend on the reduced Rayleigh number for two Ekman numbers. For both 
$Ek=3\times 10^{-4}$ and 
$Ek=10^{-5}$, the buoyancy force is approximately an order of magnitude larger than viscosity or advection at small 
$\widetilde {\textit {Ra}}$, and both viscosity and advection grow more rapidly with Rayleigh number as compared to the buoyancy force. For 
$Ek=3\times 10^{-4}$, advection becomes larger than buoyancy near 
$\widetilde {\textit {Ra}} \approx 140$ and continues growing larger than buoyancy. For 
$Ek=10^{-5}$, we do not reach large enough Rayleigh numbers for advection to become as large as buoyancy. Therefore, we do not observe a CIA balance in our simulations since the buoyancy force and advection scale differently with Rayleigh number for the parameter space we surveyed.
In figure 17, we plot the ratio of the viscous force to the buoyancy force, and the ratio of the fluctuating advective term to the buoyancy force. We find that the ratio of the viscous force to buoyancy is an increasing function of the Rayleigh number, which suggests that viscosity does not become negligible at more extreme parameters. This behaviour is in contrast to the diffusion-free scaling that is used in many previous studies, which assume that viscosity is negligible. Therefore, the trends we observe in our data do not support viscosity becoming negligible, although it is possible that this trend changes at higher values of the reduced Rayleigh number than we achieved in this study. However, Oliver et al. (Reference Oliver, Jacobi, Julien and Calkins2023) also investigated this ratio over a larger range of 
$\widetilde {\textit {Ra}}$ in the plane geometry and found similar behaviour, where it was attributed to the fact that the fluid interior acts as the dominant control on heat transport. In this sense, this measure may simply be a manifestation of the energetic constraints on the system in which work done by the viscous force must balance the work done by the buoyancy force when averaged over time and space (e.g. Siggia Reference Siggia1994).
Figure 17. Ratios of forces in the fluctuating radial momentum equation where 
$F_v=[\nabla ^2 \boldsymbol {u}']_r$, 
$F_a=[\boldsymbol {u}'\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {u}']_r$ and 
$F_b=(Ra/Pr)(r/r_0)T'$: (a) ratio of the viscous force to the buoyancy force; (b) ratio of the fluctuating advective term to the buoyancy force. The symbols are the same as defined in figure 2.
4.6. Heat flow and dissipation
Figure 18 shows the Nusselt number for all of the 
$\eta =0.35$ cases. The raw data are shown in panel (a); asymptotically rescaled data are shown in panel (b); and the Nusselt number compensated by the diffusion-free scaling is shown in panel (c). As for much of the data shown previously, in comparison to panel (a), the collapse of the data in panel (b) is indicative that the simulations are in the regime of rotationally constrained dynamics. Moreover, the force balances presented previously indicate that none of our simulations are within the buoyancy-dominated regime. The solid line in panel (b) shows the predicted scaling for diffusion-free heat transport, 
$Nu \sim \widetilde {\textit {Ra}}{\vphantom {Ra}}^{3/2}$, and the dashed line shows the scaling 
$Nu \sim \widetilde {\textit {Ra}}$. It appears that our data have not quite reached the predicted 
$Nu \sim \widetilde {\textit {Ra}}{\vphantom {Ra}}^{3/2}$ scaling, although some of our 
$Ek=3\times 10^{-6}$ cases appear to be close to this scaling. This might suggest that the 
$Nu \sim \widetilde {\textit {Ra}}{\vphantom {Ra}}^{3/2}$ scaling can be reached at lower Ekman number and higher Rayleigh number.
Figure 18. Time-averaged Nusselt number calculated at the outer boundary: (a) Nusselt number versus Rayleigh number; (b) Nusselt number versus reduced Rayleigh number; (c) compensated Nusselt number. The symbols are the same as defined in figure 2.
Figure 19(a,c) shows the dissipation by both the convective flow and the mean flow. The dissipation by the convective flow is greater than the dissipation by the mean flow for all of the cases we studied. Figure 19(b,d) shows that the dissipation follows the expected scaling of 
$Ek^{-4/3}$ and that the convective dissipation follows the scaling 
$Ek^{-4/3}$ slightly better than the mean flow. Figure 20 shows the ratio of the dissipation by the mean flow to the dissipation by the convective flow. We see that for all Ekman numbers, this ratio reaches a maximum when 
$\widetilde {\textit {Ra}} \sim 30$, and that for higher values of the reduced Rayleigh number, the ratio of mean to convective dissipation decreases. This trend might suggest that for large values of the reduced Rayleigh number, the dissipation is mainly controlled by the convection. There is also a trend where the maximum value of the ratio of the mean dissipation to the convective dissipation increases as the Ekman number is decreased. However, this trend is rather weak, with only approximately a factor of two change in the ratio while the change in Ekman number was two orders of magnitude.
Figure 19. Viscous dissipation rates for (a,b) the convective flow and (c,d) the zonal flow. The symbols are the same as defined in figure 2.
Figure 20. Ratio of the large-scale viscous dissipation rate to the small-scale viscous dissipation rate. The symbols are the same as defined in figure 2.
5. Conclusions
A systematic investigation of rotating convection in a spherical shell with stress-free boundary conditions was carried out. The scaling behaviour of various quantities was investigated for varying Rayleigh number and Ekman number, and two different aspect ratios were employed. An emphasis was placed on characterising the asymptotic nature of the system as the Ekman number is reduced. Here we have used the known scalings obtained from both the linear theory of rotating spherical convection (Jones et al. Reference Jones, Soward and Mussa2000; Dormy et al. Reference Dormy, Soward, Jones, Jault and Cardin2004) and the closely related nonlinear reduced model for rotating Rayleigh–Bénard convection (e.g. Sprague et al. Reference Sprague, Julien, Knobloch and Werne2006). Overall, we find good agreement between the asymptotic scalings and the DNS, though the parameter space explored with the simulations is limited, especially when compared with natural systems.
In general, we find that the asymptotic nature of the system can be demonstrated when plotting asymptotically rescaled quantities as a function of the reduced Rayleigh number, 
$\widetilde {\textit {Ra}}$. Using this approach, we find that the convective flow speeds, as measured by the large-scale Reynolds number, scale as 
$Re_c = O(Ek^{-1/3})$ for both thin and thick shell simulations. This scaling can be deduced from the governing equations by requiring the convection to be geostrophically balanced to leading order, with the buoyancy force entering the next asymptotic order. As discussed by previous work, the amplitude of the zonal flow is limited solely by large-scale viscous diffusion; thus, we do not anticipate a diffusion-free scaling for this component of the velocity field at any combination of Ekman or Rayleigh number in the rotationally constrained regime. By balancing the large-scale viscous diffusion of zonal momentum with the Reynolds stresses, and noting that 
$Re_c = O(Ek^{-1/3})$, we deduce that the zonal flow speeds scale as 
$Re_z = O(Ek^{-2/3})$. The numerical data support a trend towards this scaling as 
$Ek \rightarrow 0$. Moreover, we find that because of this diffusion-Reynolds-stress-balance, the correlations in the convective flow field are independent of the Ekman number. As noted in previous work (Christensen Reference Christensen2002; Gastine & Wicht Reference Gastine and Wicht2012), these correlations decrease rapidly with increasing Rayleigh number (approximately as 
$\sim \widetilde {\textit {Ra}}{\vphantom {Ra}}^{-3/2}$). It is also interesting to consider how the boundary conditions used would affect our scalings. Data from Gastine et al. (Reference Gastine, Wicht and Aubert2016) suggest that cases with no-slip boundaries follow the same 
$Re_c = O(Ek^{-1/3})$ scaling we observed in this paper. However, zonal flows with no-slip boundaries tend to be significantly suppressed as compared with stress-free boundaries. This suggests that cases with no-slip boundary conditions either need to reach smaller Ekman number to achieve the same scaling for the zonal Reynolds number as cases with stress-free boundary conditions, or it suggests that the zonal flows for no-slip boundary conditions are asymptotically weaker than the zonal flows with stress-free boundary conditions due to the presence of Ekman pumping (e.g. Soward Reference Soward1977).
Kinetic energy spectra exhibit a trend with increasing Rayleigh number in which energy builds at nearly all spatial scales, with no clear preference for scales significantly larger than the linear convective instability, as indicated by an 
$Ek^{1/3}$ or 
$Ek^{2/9}$ scaling. Using the length scale calculation approach of Christensen & Aubert (Reference Christensen and Aubert2006) and Gastine et al. (Reference Gastine, Wicht and Aubert2016), we find that the dominant convective length scale decreases with increasing 
$\widetilde {\textit {Ra}}$ and scales approximately as 
$Ek^{1/3}$. At the largest accessible Rayleigh numbers, but also the largest Ekman numbers, this length scale appears to saturate. The Taylor microscale also scales as 
$Ek^{1/3}$ and shows a very similar trend in comparison to the spherical-harmonic-based length scale; after initially decreasing with increasing 
$\widetilde {\textit {Ra}}$, it tended to level off at the largest accessible Rayleigh numbers. Moreover, both of these length scales remain comparable in value across the entire range of parameters that were simulated, and the collapse of the asymptotically rescaled length scales suggests that this trend continues for more extreme parameter values. The largest length scales present in the fluctuating kinetic energy spectra scale weaker than the primary convective instability and may be consistent with the 
$Ek^{2/9}$ radial scale from linear theory (Dormy et al. Reference Dormy, Soward, Jones, Jault and Cardin2004). This 
$Ek^{2/9}$ length scale increases with increasing Rayleigh number, and approximately follows the scaling 
$\ell ^{\prime }_{peak} \sim Ro_c^{1/4}$. Our results suggest that there are multiple length scales present in rotating spherical convection, and that these length scales exhibit different dependencies on the input parameters. Nevertheless, all of these length scales appear to depend on the Ekman number and are therefore viscously controlled to some degree.
Several previous investigations have used no-slip boundary conditions that yield different scaling behaviour of the length scales in comparison with the present study. Gastine et al. (Reference Gastine, Wicht and Aubert2016) computed length scales over a broad range of parameters with 
$\eta =0.6$. These authors used a similar definition of the spherical-harmonic-based length scale (2.7), but also included 
$m=0$ in their calculations. In general, they found non-monotonic behaviour in their values of 
$\ell _{sh}$, though they did find an increase in the convective length scale provided 
$Ra$ was sufficiently large and 
$Ek$ was sufficiently small. We have found that including the 
$m=0$ component in the length scale calculation causes the length scale to increase with increasing Rayleigh number over some range of 
$\widetilde {\textit {Ra}}$, suggesting it is necessary to remove this component when deducing length scale behaviour of the small-scale convection. Using constant heat flux thermal boundary conditions, Long et al. (Reference Long, Mound, Davies and Tobias2020) also used the full spectrum version of (2.7) (i.e. they included 
$m=0$) and found an increase in 
$\ell _{sh}$ at sufficiently small Ekman number.
The mechanical boundary conditions that are employed in the model might play an important role in the scaling behaviour of the convective length scales. No-slip boundary conditions result in the formation of Ekman layers, which tend to significantly damp the large-scale flows (Soward Reference Soward1977; Christensen Reference Christensen2002). One might expect that the strong radial shear generated by the zonal flow limits the size to which convective structures can grow, therefore yielding distinct scaling behaviour in comparison with simulations in which the zonal flow is suppressed. It is also possible that a transition to a regime in which the convective length scale increases occurs in a parameter regime outside of that used in the presented study. However, our data show no evidence towards such a transition. Another possibility for the observed differences in the scaling behaviour of the length scale could be that the axisymmetric (
$m=0$) component of the flow is influencing the evolution of the length scale even when no-slip conditions are used. While the zonal flow is smaller for no-slip simulations, zonal flows can still be present and could influence 
$l_{sh}$. It is also worth noting that Gastine & Aurnou (Reference Gastine and Aurnou2023) found that heat transport for convection in rotating spherical shells is strongly dependent on latitude. This suggests that other quantities, such as the length scale and Reynolds number, may also depend strongly on latitude. In this case, it may be illuminating to study the Rayleigh number dependence of the length scale as a function of latitude rather than as a globally averaged quantity.
The asymptotic scaling of various terms in the fluctuating momentum and heat equations was investigated. In particular, for the fluctuating momentum equation, we found that the viscous force and the buoyancy force scale as 
$Ek^{-1}$, and that the fluctuating–fluctuating advective term scales approximately as 
$Ek^{-1}$, or perhaps slightly stronger than 
$Ek^{-1}$ (with our large-scale viscous diffusion non-dimensionalisation), whereas the Coriolis force scales as 
$Ek^{-4/3}$. These scalings again agree with asymptotic theory up to a small difference in the advective term. We find that viscosity does not become smaller compared to buoyancy as the reduced Rayleigh number is increased, so we do not observe a trend in which viscosity becomes negligible. We also do not observe a balance between buoyancy and advection for any of our simulations; the advective term is always growing more rapidly with Rayleigh number than the buoyancy force. Indeed, for the largest Ekman number considered for 
$\eta =0.35$ (
$Ek=3\times 10^{-4}$), the magnitude of inertia surpasses that of the buoyancy force, though this behaviour is outside the regime of rapidly rotating convection. This finding suggests that, at least for the parameter space observed in this study, there is no CIA balance that occurs.
The force balances that were computed for our spherical simulations are similar to the results found by Aguirre Guzmán et al. (Reference Aguirre Guzmán, Madonia, Cheng, Ostilla-Mónico, Clercx and Kunnen2021) for rapidly rotating convection with no-slip boundaries in Cartesian coordinates. In particular, Aguirre Guzmán et al. (Reference Aguirre Guzmán, Madonia, Cheng, Ostilla-Mónico, Clercx and Kunnen2021) found that for large Rayleigh numbers, the viscous force was approximately as large as the buoyancy force and that the inertial term would be larger than buoyancy. Simulations of the asymptotic model for the plane layer also show that the ratio of the viscous force to the buoyancy force approaches unity in the high-Rayleigh-number regime (Maffei et al. Reference Maffei, Krouss, Julien and Calkins2021; Oliver et al. Reference Oliver, Jacobi, Julien and Calkins2023). This similarity indicates that the Rayleigh number dependence of the forces may not depend too sensitively on the geometry of the system. Some dynamo simulations also show similar scalings compared with the hydrodynamic cases. Yan & Calkins (Reference Yan and Calkins2022) carried out numerical simulations of a dynamo in a plane layer and also found that viscosity and buoyancy follow the same Ekman number dependence with no indication of a regime where viscosity becomes small compared to buoyancy. However, Yadav et al. (Reference Yadav, Gastine, Christensen, Wolk and Poppenhaeger2016) carried out numerical simulations of spherical dynamo cases with no-slip boundaries and found that the ratio of viscosity to buoyancy decreased as the Ekman number was decreased at constant 
$Ra/Ra_c$. Some other papers studied the spectral decomposition of the forces for spherical dynamo cases and also found that viscosity was much smaller than buoyancy at small values of the Ekman number (e.g. Schwaiger, Gastine & Aubert Reference Schwaiger, Gastine and Aubert2019, Reference Schwaiger, Gastine and Aubert2021). These dynamo simulations where viscosity seems to become small compared to buoyancy at small values of the Ekman number could be following a different asymptotic scaling than we observe in the hydrodynamic cases in this paper, although dynamo plane layer simulations seem to follow similar scaling laws as the hydrodynamic model (Yan & Calkins Reference Yan and Calkins2022).
Supplementary material and movie
Supplementary material and movie are available at https://doi.org/10.1017/jfm.2024.78.
Acknowledgements
The authors gratefully acknowledge feedback from three anonymous referees that helped to significantly improve the manuscript.
Funding
This research was supported with funding from the National Science Foundation through grant EAR-1945270. The Stampede2 and Anvil supercomputers at the Texas Advanced Computing Center and Purdue University, respectively, were made available through allocation PHY180013 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants #2138259, #2138286, #2138307, #2137603 and #2138296. Simulations were also conducted on the Summit and Alpine supercomputers at the University of Colorado, Boulder. Summit is a joint effort of the University of Colorado Boulder and Colorado State University, supported by NSF awards ACI-1532235 and ACI-1532236.
Declaration of interests
The authors report no conflict of interest.
Appendix. Numerical simulation data
Table 1. Summary of the cases with 
$\eta =0.35$. Here 
$N_r$ is the number of radial grid points used in the simulation and 
$l_{max}$ is the maximum spherical harmonic degree used in the simulation. The standard deviation in time of given quantities is shown after the ‘
$\pm$’.
Table 2. Summary of the cases with 
$\eta =0.7$. Here, 
$N_r$ is the number of radial grid points used in the simulation and 
$l_{max}$ is the maximum spherical harmonic degree used in the simulation. The standard deviation in time of given quantities is shown after the ‘
$\pm$’.
Open access
