1. Introduction
Convection is commonplace in the interiors of stars and planets, in which rotation plays a dominant role (e.g. Zhang & Liao Reference Zhang and Liao2017). Despite a concerted study of this system, both numerically and theoretically, over many decades, our understanding of the process remains incomplete. It is of great interest to properly characterise both the morphology of fluid flow and the mechanism of heat transfer as a function of the dynamical regime as controlled by the relevant non-dimensional numbers. The highly supercritical regime, which is both the most challenging and the most relevant in planetary settings, remains poorly understood.
Three classes of problem have been tackled in order to make progress. By far the most studied is the problem of plane layer Rayleigh–Bénard convection, in both non-rotating and rotating settings (Plumley & Julien Reference Plumley and Julien2019). The second class is rotating spherical shell convection, naturally motivated by the geometry of planets that possess solid cores; such a geometry lends itself naturally to efficient numerical methods (Glatzmaier Reference Glatzmaier2014; Marti, Calkins & Julien Reference Marti, Calkins and Julien2016). This problem has been studied using laboratory experiments (Aubert et al. Reference Aubert, Brito, Nataf, Cardin and Masson2001) and numerical simulations in both Boussinesq (Christensen Reference Christensen2002; Aurnou, Heimpel & Wicht Reference Aurnou, Heimpel and Wicht2007) and anelastic settings (Gastine, Wicht & Aurnou Reference Gastine, Wicht and Aurnou2013), with a variety of gravity profiles (Gastine, Wicht & Aubert Reference Gastine, Wicht and Aubert2016) and boundary conditions (Mound & Davies Reference Mound and Davies2017; Long et al. Reference Long, Mound, Davies and Tobias2020). The results of the first two problems have been recently reviewed by Jones (Reference Jones2015), Aurnou et al. (Reference Aurnou, Calkins, Cheng, Julien, King, Nieves, Soderlund and Stellmach2015) and Plumley & Julien (Reference Plumley and Julien2019), and we refer the reader to these review articles.
The third class of study, which is the subject of the present research, is the canonical problem of rotating convection in a whole sphere in which an internal heat source is present; the strength of the heating is measured by the Rayleigh number 
$Ra$. Although early theoretical work concentrated on the whole-sphere geometry (Chandrasekhar Reference Chandrasekhar1961; Roberts Reference Roberts1968; Busse Reference Busse1970; Kida Reference Kida1994; Jones, Soward & Mussa Reference Jones, Soward and Mussa2000; Zhang & Liao Reference Zhang and Liao2004) and the determination of the critical Rayleigh number for convection (
$Ra_c$), there has been much less numerical work focusing on this area. Notwithstanding the studies that reduce the size of the inner core to negligible size without actually removing it, notable studies using a true full-sphere geometry are those of Jones et al. (Reference Jones, Soward and Mussa2000), Sánchez, Garcia & Net (Reference Sánchez, Garcia and Net2016) and Kaplan et al. (Reference Kaplan, Schaeffer, Vidal and Cardin2017), which mainly concerned themselves with the onset of convection or subcritical convection. More recently, Guervilly, Cardin & Schaeffer (Reference Guervilly, Cardin and Schaeffer2019) examined the convective length scales in a rotating sphere using both quasi-geostrophic and fully three-dimensional methods in the regime of very low Ekman number and small Prandtl number, albeit with a weak thermal forcing (around the onset of convection). Convection was studied experimentally in the full sphere by Chamberlain & Carrigan (Reference Chamberlain and Carrigan1986) using centrifugal gravity, though practical considerations precluded the exploration of very supercritical Rayleigh numbers: in that study, 
$Ra \leqslant 2.5Ra_c$. The highly supercritical regime of convection in a whole sphere is of great relevance to the dynamics of the Earth's core prior to inner core nucleation and to some fully convective stars, yet the problem remains largely unexplored.
In the present work we concentrate on the problem in which the diffusivities of heat and momentum are equal, leading to a so-called Prandtl number (
$Pr$) of unity. At fixed viscosity this leaves the rotation period (measured by the Ekman number 
$E$) and the strength of internal heating (characterised by 
$Ra$) as the only control parameters. We perform a set of direct numerical simulations with moderate Ekman numbers (
$E\sim 10^{-5}$), but over a wide range of Rayleigh numbers (from near the onset up to 
$10^{4}Ra_c$).
A primary focus of our work has been to delineate the regime boundaries that occur on increasing the Rayleigh number. We make use of the significant step forward made by Gilman (Reference Gilman1977) who identified the relevance of the convective Rossby number 
$Ro_c$ in supercritical convection. In terms of the aforementioned parameters, 
$Ro_c=\sqrt {RaE/Pr}$, and we find that 
$Ro_c$ is indeed the control parameter that defines the regime boundaries. A relaxation oscillation regime occurs just beyond the onset of convection and the boundary with the geostrophic turbulence (GT) regime that follows subsequently is well described by a critical value of 
$Ro_c$. Likewise, we find a second bifurcation to a regime with a large-scale vortex (LSV), which is perfectly described by a constant value of 
$Ro_c$. Thus our results support the significance of this parameter. Conversely, when looking quantitatively within these regimes, we find that this parameter alone is insufficient to fully describe the scaling of flow speed and heat transfer. In particular, we find that the formation of LSVs appears to reduce the efficiency of the heat transfer and the convective flow speed.
The identification of axisymmetric LSVs in spherical rotating convection is entirely new. We argue that the omission of an inner core removes the obstacle to such geostrophic structures that otherwise would be present. Large-scale coherent vortices have been previously observed in planar rotating convection in both the anelastic case (Chan Reference Chan2007; Käpylä, Mantere & Hackman Reference Käpylä, Mantere and Hackman2011; Chan & Mayr Reference Chan and Mayr2013) and Boussinesq case (Julien et al. Reference Julien, Rubio, Grooms and Knobloch2012b; Favier, Silvers & Proctor Reference Favier, Silvers and Proctor2014; Guervilly, Hughes & Jones Reference Guervilly, Hughes and Jones2014; Rubio et al. Reference Rubio, Julien, Knobloch and Weiss2014; Stellmach et al. Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014). The formation of LSV is attributed to an inverse cascade, i.e. upscale energy transfer, in rotating turbulence. In spherical containers, large-scale coherent structures may be associated with mean zonal flows (Christensen Reference Christensen2002). We find a reversal of the mean zonal flow when the LSVs are formed. Such a reversal has been observed in thin spherical shells for both Boussinesq convection (Aurnou et al. Reference Aurnou, Heimpel and Wicht2007) and anelastic convection (Gastine et al. Reference Gastine, Wicht and Aurnou2013) when the Rayleigh number is sufficiently large, but never before in a full sphere. Effects of boundary conditions on the formation of LSVs are also briefly discussed.
2. Numerical models
2.1. Governing equations
We consider Boussinesq convection in a full sphere with radius of 
$r_o$, which rotates at 
$\boldsymbol {\varOmega }=\varOmega \hat {\boldsymbol {z}}$. The sphere is filled with an incompressible fluid of density 
$\rho$, viscosity 
$\nu$ and thermal diffusivity 
$\kappa$. Convective motions are driven by a homogeneous internal heat source 
$S$, under a gravitational field 
$\boldsymbol {g}=-g_o\boldsymbol {r}/r_o$. In the absence of convection, the basic state temperature is determined by the conduction alone and is given as
\begin{equation} {\widetilde{T_b}}=\frac{\beta}{2}(r_o^{2}-r^{2}), \end{equation}
where 
$\beta =S/(3 k)$ in which 
$k$ is thermal conductivity.
Using the radius 
$r_o$ as the length scale, the viscous diffusion time 
$r_o^{2}/\nu$ as the time scale and 
$\beta r_o^{2}$ as the temperature scale, the system is governed by the following dimensionless equations in the rotating frame (Marti et al. Reference Marti2014):
\begin{gather} E \left(\frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}\right)+{\hat{\boldsymbol{z}}} \times \boldsymbol{u}={-}\boldsymbol{\nabla} p +\frac{Ra}{Pr}T\boldsymbol{r}+E\nabla^{2} \boldsymbol{u}, \end{gather}
\begin{gather}Pr \left(\frac{ \partial T}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} T \right)=\nabla^{2} (T-T_b), \end{gather}
\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}
where 
$\boldsymbol {u}$ is the fluid velocity, 
$p$ is the reduced pressure, 
$T$ is the total temperature and 
$T_b$ is the dimensionless basic temperature given as
\begin{equation} T_b=\tfrac{1}{2}(1-r^{2}). \end{equation} The dimensionless parameters are the Ekman number 
$E$, the rotationally modified Rayleigh number 
$Ra$ and the Prandtl number 
$Pr$:
\begin{equation} E=\frac{\nu}{2\varOmega r_o^{2}}, \quad Ra=\frac{\alpha \beta g_o r_o^{3}}{2 \varOmega \kappa},\quad Pr=\frac{\nu}{\kappa}, \end{equation}
where 
$\alpha$ is the thermal expansion coefficient. The Prandtl number is fixed to be unity (
$Pr=1$) in this study. Note that 
$Ra$ in a full sphere is naturally flux-based as 
$\beta$ is related to the total heat flux 
$\tilde {Q}$ through the boundary by
\begin{equation} \tilde{Q}=\frac{4{\rm \pi} r_0^{3}}{3}S=4 {\rm \pi}r_0^{3} k \beta, \end{equation}
regardless of thermal boundary conditions, and the non-dimensional heat flux 
$q$ through the boundary is 
$q=-\boldsymbol {\nabla } T$. Here 
$Ra$ is related to the convectional Rayleigh number 
$Ra_T=\alpha g_0\Delta Tr_o^{3}/\nu \kappa$ based on the temperature difference by 
$Ra=2NuERa_T$, where 
$Nu$ is the Nusselt number (see the definition given in (2.16) below).
As pointed out by Gilman (Reference Gilman1977), the convective Rossby number, which is a combination of 
$E$, 
$Ra$ and 
$Pr$, is a key control parameter in supercritical rotating convection. Based on our definition, the convective Rossby number is given as
\begin{equation} Ro_c=\sqrt{\frac{ERa}{Pr}}=\frac{\sqrt{\alpha \beta g_0 r_o}}{2\varOmega}, \end{equation}
which does not depend on the viscosity 
$\nu$ and the thermal diffusivity 
$\kappa$. The convective Rossby number is the ratio between the rotation time scale and the free-fall time scale, and characterises the global force balance between the buoyancy force with respect to the Coriolis force (Gilman Reference Gilman1977). Christensen (Reference Christensen2002) introduced a similar diffusivity-free parameter 
$Ra^{*}=ERa/Pr$, which is the square of the convective Rossby number, i.e. 
$Ra^{*}=Ro_c^{2}$.
We adopt the stress-free boundary condition for the velocity 
$\boldsymbol {u}$ and a fixed temperature (
$T=0$ for simplicity) at the boundary of the sphere in most of the numerical simulations. The no-slip boundary condition and a fixed heat flux, i.e. 
$\partial T/\partial r=-1.0$, boundary condition are also briefly considered.
2.2. Numerical method
The governing equations (2.2)–(2.4) subject to the boundary conditions are numerically solved using a spectral method in a whole sphere (Marti & Jackson Reference Marti and Jackson2016). As we consider an incompressible fluid, the velocity 
$\boldsymbol {u}$ is decomposed into toroidal and poloidal components:
\begin{equation} \boldsymbol{u}= \boldsymbol{\nabla} \boldsymbol{\times}(\mathcal{T} \boldsymbol{r}) +\boldsymbol{\nabla} \boldsymbol{\times}\boldsymbol{\nabla} \boldsymbol{\times}(\mathcal{P} \boldsymbol{r}). \end{equation}
The toroidal 
$\mathcal {T}$ and poloidal 
$\mathcal {P}$ scalar fields and the temperature field 
$T$ are expanded in terms of spherical harmonic expansion on spherical surfaces and the so-called Jones–Worland polynomials in the radial direction. The spectral expansion is truncated up to spherical harmonics of degree 
$L=255$ and order 
$M=255$, and up to 
$N=127$ for the Jones–Worland polynomials. A critical element of the technique is its graceful handling of the origin at the centre of the sphere, such that the Courant–Friedrichs–Lewy condition is not unduly stringent in this region. Time-stepping has been performed using a second-order Runge–Kutta scheme, with adaptive time steps. The typical time step is about 
$10^{-7}$ (viscous time scale) for the most demanding calculations. We refer the reader to Marti & Jackson (Reference Marti and Jackson2016) for more details regarding the numerical scheme. The numerical code has been benchmarked for several hydrodynamic and magnetohydrodynamic problems in a full sphere, including the rotating thermal convection problem (Marti et al. Reference Marti2014).
2.3. Diagnostics
In order to qualitatively describe numerical results, the following quantities derived from the velocity 
$\boldsymbol {u}$ and the temperature 
$T$ will be used.
The total kinetic energy is given by
\begin{equation} \mathcal{E} =\tfrac{1}{2}\int_v \boldsymbol{u}^{2} \,\mathrm{d}V. \end{equation}The development of strong zonal flow is a distinct feature in the nonlinear regime; hence we can decompose the total energy into zonal and non-zonal parts:
\begin{equation} \mathcal{E}_{{zonal}} =\tfrac{1}{2}\int_v \left( u_{\phi}^{0}\right)^{2} \,\mathrm{d}V, \quad \mathcal{E}_{{non}}=\mathcal{E}-\mathcal{E}_{{zonal}}, \end{equation}
where 
$u_{\phi }^{0}$ is the axisymmetric (
$m=0$) component of the azimuthal velocity 
$u_\phi$, i.e. zonal flow. We define the Rossby number as
\begin{equation} Ro=\frac{U_{rms}}{2\varOmega r_0}, \end{equation}
where 
$U_{rms}$ is the dimensional root mean square velocity. Base on the non-dimensionalisation we used, the Rossby number is related to the non-dimensional energy as
\begin{equation} Ro=E\sqrt{\frac{2\mathcal{E}}{V}}, \end{equation}
where 
$V=4 {\rm \pi}/3$ is the non-dimensional volume of the sphere. Accordingly, we have the zonal Rossby number
\begin{equation} Ro_{{zon}}=E\sqrt{\frac{2\mathcal{E}_{zonal}}{V}} \end{equation}and the non-zonal Rossby number
\begin{equation} Ro_{{non}}=E\sqrt{\frac{2\mathcal{E}_{{non}}}{V}}. \end{equation}Heat transport owing to convection can be measured by the Nusselt number. For the convection of internal heating in a full sphere, the Nusselt number can be estimated by the drop of the temperature with respect to the basic state at the centre (Guervilly & Cardin Reference Guervilly and Cardin2016; Kaplan et al. Reference Kaplan, Schaeffer, Vidal and Cardin2017):
\begin{equation} Nu=\frac{T_b(r=0)}{T(r=0)}.\end{equation}Apart from these global quantities, we often use the axial vorticity and the zonal velocity to visualise flows. The axial vorticity is calculated by
\begin{equation} \omega_z=2 E \hat{\boldsymbol{z}}\boldsymbol{\cdot} (\boldsymbol{\nabla} \boldsymbol{\times} \boldsymbol{u}). \end{equation}
Note that the prefactor of 
$2E$ leads to the vorticity in units of the rotation frequency 
$\varOmega$ given the non-dimensionalisation we have used. The zonal velocity is also rescaled, i.e. 
$U_\phi ^{0}=2E u_\phi ^{0}$, such that 
$U_\phi ^{0}$ is in the units of 
$\varOmega r_0$. The mean zonal flow 
$\overline {U_\phi ^{0}}$ is time-averaged 
$U_\phi ^{0}$.
3. Results
Numerical simulations in this study are listed in detail in table 1 in appendix A. We fix 
$Pr=1$ in this study and vary the control parameters 
$E$ and 
$Ra$, leading to a variety of dynamical regimes. Stress-free and fixed temperature boundary conditions are assumed unless stated otherwise. Figure 1 summarises simulations as a function of 
$(E, Ra)$ with symbols indicating different flow regimes, which are discussed in § 3.1.
Figure 1. Regime diagram in the plane of 
$(E, Ra)$. Different symbols represent different flow regimes. Open circles correspond to the steadily drifting regime near the onset; filled circles correspond to the relaxation oscillation regime; blue diamonds correspond to the GT regime; filled triangles correspond to the LSV regime; and the open triangle in the top right corresponds to the non-rotating regime.
Table 1. Summary of numerical simulations. Symbols # and * denote the cases using different boundary conditions (see text for details).
3.1. Flow regimes
Convection takes place only when 
$Ra$ is larger than a critical value 
$Ra_c$. At 
$Pr=1$ with the stress-free and the fixed temperature boundary conditions, the critical value is given as 
$Ra_c\approx 4.1 E^{-1/3}+17.8$ (black solid line in figure 1) according to the asymptotic theory of the linear onset (Jones et al. Reference Jones, Soward and Mussa2000). Our fully nonlinear numerical simulations are consistent with the linear theory regarding the onset of convection. Figure 2 shows a snapshot of the axial vorticity and the zonal flow when 
$Ra$ is marginally above the critical value (
$Ra=220\approx 1.05Ra_c$) at 
$E=10^{-5}$. The flow near the onset is in the form of the so-called thermal Rossby waves, which are quasi-geostrophic, i.e. nearly invariant along the rotation axis (figure 2b). The kinetic energy becomes quasi-steady after the initial growing phase. The azimuthal wavenumber 
$m=14$ is also in agreement with the linear theory. We also see that a weak zonal flow develops (figure 2c) with differential rotation (prograde near the equator and retrograde in the inner region) due to the nonlinear interactions of the linear onset mode (Zhang Reference Zhang1992).
Figure 2. Flows near the onset at 
$E=10^{-5}$, 
$Pr=1$, 
$Ra=220\approx 1.05Ra_c$. (a) Axial vorticity in the equatorial plane and (b) in the meridional plane; (c) zonal flow in the meridional plane.
When 
$Ra$ is a few times larger than the critical value, convective motions become quasi-periodic. Figure 3 shows an example of the relaxation oscillation regime at 
$E=10^{-5}$ and 
$Ra=1000\approx 4.8 Ra_c$. We can see that the time evolution of the energy exhibits a quasi-periodic manner (figure 3a). The growth of the onset mode (figure 3b) leads to subsequent instabilities, resulting in more vigorous and complex convection (figure 3c). A movie of the axial vorticity for this case is provided (see supplementary movie 1). As the non-zonal energy approaches the maxima, there is a burst of the zonal energy, followed by a relaminarisation process (decreasing of both the zonal and non-zonal energy). During these cycles, the non-zonal energy can vary over several orders of magnitude between peaks and troughs. Similar relaxation oscillations have been observed in numerical simulations of convection in spherical shells, and are also referred to as convective bursts (Grote & Busse Reference Grote and Busse2001; Busse Reference Busse2002; Christensen Reference Christensen2002; Heimpel & Aurnou Reference Heimpel and Aurnou2012). The relaxation oscillations can be attributed to the competition between the nonlinear effect and the viscous effect (Christensen Reference Christensen2002). The oscillatory regime exists only at intermediate Rayleigh numbers (filled circles in figure 1).
Figure 3. An example of the oscillatory flow regime at 
$E=10^{-5}$, 
$Ra=1000\approx 4.8 Ra_c$. (a) Time evolution of the zonal and non-zonal kinetic energy. Time is in units of the viscous time scale. (b,c) Axial vorticity in the equatorial plane and in the meridional plane at the instants indicated by the blue dot and red dot in (a), respectively. (A movie of this time series is provided; see supplementary movie 1 available at https://doi.org/10.1017/jfm.2020.1151.)
As we increase 
$Ra$, time evolution of the kinetic energy becomes less periodic and eventually becomes irregular fluctuations, but the flow structures are still elongated along the rotation axis due to the constraint of rotation. We refer to this regime as the GT regime. The transition from the relaxation oscillation to the GT regime seems to be controlled by the convective Rossby number, i.e. 
$Ro_c\approx 0.2$ (blue dashed line in figure 1). Figure 4 shows a typical case of the GT regime at 
$E=10^{-5}$ and 
$Ra=10^{5}$. In this regime, convective motions are charaterised by space-filling and sustained turbulence (figure 4a,b). The axial vorticity in the meridional plane (figure 4b) shows elongated structures along the rotation axis due to the constraint of rotation in this regime. This is a typical flow regime in rotating convection when 
$Ra$ is supercritical, which has been observed in a wide range of numerical simulations and laboratory experiments (Aurnou et al. Reference Aurnou, Calkins, Cheng, Julien, King, Nieves, Soderlund and Stellmach2015).
Figure 4. An example of the GT regime at 
$E=10^{-5}$ and 
$Ra=10^{5}$. (a) Axial vorticity in the equatorial plane. (b) Axial vorticity in the meridional plane. (c) Zonal flow in the meridional plane. (d,e) Energy spectra as a function of the azimuthal wavenumber 
$m$ and spherical harmonic degree 
$l$. Black dashed lines represent the 
$-5/3$ power law purely for reference.
Apart from small-scale turbulence, geostrophic zonal flows (figure 4c) are developed in this case, with a strong prograde jet near the equator and mainly retrograde flows in the bulk. The significance of the zonal flow is clearly evident in the energy spectra (figure 4d,e) as a function of the azimuthal wavenumber 
$m$ and spherical harmonic degree 
$l$. The kinetic energy is dominated by the 
$m=0$ toroidal component, which is basically the zonal flow. In the 
$l$ spectra, 
$l=0$ and 
$l=2$ toroidal components contribute the most kinetic energy in the flow.
Further increasing 
$Ra$ leads to a completely different flow regime (red triangles in figure 1), which is characterised by the formation of LSV. The transition from the GT regime to the LSV regime is also determined by the convective Rossby number, i.e. 
$Ro_c\approx 1.5$ (the red dashed line in figure 1). Figure 5 shows an example of the LSV at 
$E=10^{-5}$ and 
$Ra=3\times 10^{5}$. We can see a notable cyclonic vortex located at the centre in the equatorial plane (figure 5a) and extended along the rotation axis (figure 5b). The formation process of the LSV is further discussed in § 3.2. The vortex is accompanied by a strong geostrophic zonal flow (figure 5c). However, the direction of the zonal flow is reversed compared to the previous regime. The zonal flow is retrograde near the equator and prograde in the inner region. The kinetic energy is dominated by the 
$m=0$ toroidal component, which corresponds to the columnar vortex and the zonal flow. The sawtooth-shaped curves in the 
$l$ spectra at large 
$l$ reflect the equatorially symmetric nature of the LSV and the zonal flow.
Figure 5. An example of the LSV regime at 
$E=10^{-5}$ and 
$Ra=3\times 10^{5}$. (a) Axial vorticity in the equatorial plane. (b) Axial vorticity in the meridional plane. (c) Zonal flow in the meridional plane. (d,e) Energy spectra as a function of the azimuthal wavenumber 
$m$ and spherical harmonic degree 
$l$. Black dashed lines represent the 
$-5/3$ power law purely for reference.
At sufficiently large 
$Ra$, large-scale coherent structures break down and convection tends to behaviour similar to that in a non-rotating regime. Figure 6 shows an example in this regime at 
$E=6.0\times 10^{-5}$ and 
$Ra=2.0\times 10^{6}\approx 1.6\times 10^{4} Ra_c$ (corresponding to the open triangle in figure 1). Unlike previous cases, axial vorticity in the equatorial and meridional planes does not exhibit obvious anisotropy (figure 6a,b). The zonal flow becomes much less geostrophic compared with previous cases (figure 6c). Energy spectra in both 
$m$ and 
$l$ are in agreement with the Kolmogorov scaling of 
$k^{-5/3}$ over a wide range of wavenumbers (figure 6d,e). This non-rotating regime exists at very large 
$Ra$ and it is computationally demanding to simulate this regime at lower 
$E$. As our focus in this study is the LSV, we do not explore this regime further.
Figure 6. An example of LSV breakdown at 
$E=6.0\times 10^{-5}$ and 
$Ra=2.0\times 10^{6}$. (a) Axial vorticity in the equatorial plane. (b) Axial vorticity in the meridional plane. (c) Zonal flow in the meridional plane. (d,e) Energy spectra as a function of the azimuthal wavenumber 
$m$ and spherical harmonic degree 
$l$. Black dashed lines represent the 
$-5/3$ power law.
3.2. Formation of the LSV
We now turn to analyse the formation process of the LSV. Figure 7 shows the time evolution of the zonal energy and the non-zonal energy (figure 7a), as well as the kinetic energy carried by several different 
$m$ during the growth phase of the LSV (figure 7b) at 
$E=10^{-5}$ and 
$Ra=3.0\times 10^{5}$ (corresponding to the case in figure 5). The simulation started from a saturated state at lower 
$Ra$ in the GT regime. We can see that the zonal energy overtakes the non-zonal energy after the transient stage and then the non-zonal energy drops. The simulation saturates to a state of which the zonal energy is dominant over the non-zonal energy.
Figure 7. Time series of energy at 
$E=10^{-5}$ and 
$Ra=3.0\times 10^{5}$. (a) Zonal and non-zonal energy. (b) Energy contained in several different azimuthal wavenumbers 
$m$ during the formation stage of the LSV. Time is in units of the viscous time scale.
The energy transfer from a small-scale flow to a large-scale flow is probably best illustrated in figure 7(b). While the energy in the 
$m=0$ component keeps growing and becomes dominant, the energy in other components initially grows but subsequently drops. We note that the larger the azimuthal wavenumber 
$m$, the earlier the energy drop, suggesting a successive energy transfer from larger 
$m$ to smaller 
$m$. Similar upscale energy transfer has been reported in rotating Rayleigh–Bénard convection, leading to box-sized vortices (Guervilly et al. Reference Guervilly, Hughes and Jones2014; Rubio et al. Reference Rubio, Julien, Knobloch and Weiss2014). In the full sphere we consider here, the upscale energy transfer leads to a large cyclonic vortex located at the centre and strong zonal flow (figure 5).
Figure 8 shows a sequence of snapshots (corresponding to the black dots in figure 7b) of the axial vorticity and the zonal flow during the formation stage of the LSV at 
$E=10^{-5}$ and 
$Ra=3.0\times 10^{5}$. We clearly see the process of vortex merging from small scales (figure 8a) to a few prominent vortices (figure 8b,c), which further merge into one dominant cyclonic vortex (figure 8d). The eye of the cyclone spirals inwards and eventually settles at the centre (figure 8e,f). The vortex is elongated along the rotation axis as we have shown in figure 5(b). As the vortices merge and migrate inwards, the prograde zonal flow also shifts inwards, suggesting an inwards transport of angular momentum. We further discuss the mean zonal flow in § 3.3.
Figure 8. Snapshots (corresponding to the black dots in figure 7b) of the axial vorticity in the equatorial plane and zonal flow in the meridional plane at 
$E=10^{-5}$ and 
$Ra=3.0\times 10^{5}$. (A movie of this time series is provided; see supplementary movie 2.)
Figure 9 shows the same snapshots as in figure 8 but for the temperature field. We also see the merging process from small-scale to larger-scale structures of the temperature field, culminating with a relatively hot columnar aggregation around the rotation axis as the LSV is formed. The instantaneous Nusselt number 
$Nu$ is indicated at the bottom of each snapshot. We note that the formation of the LSV results in a significant drop of 
$Nu$, e.g. from 
$Nu=56.9$ (in figure 9e) to 
$Nu=25.11$ (in figure 9f). This is perhaps due to the fact that a hot column of fluid is trapped around the rotation axis by the LSV, which reduces the amount of heat transported by the convection. We further discuss effects of the LSV on the heat transport in § 3.4.
Figure 9. As for figure 8 but for the temperature in the equatorial plane and in the meridional plane. Instantaneous Nusselt number 
$Nu$ is also indicated.
3.3. Mean zonal flows
As we have shown, zonal flows are always developed once the convective motions set in. In this section, we compare the zonal flow in different regimes and investigate how the zonal flow varies depending on the control parameters.
Figure 10 shows the mean zonal flow, i.e. time-averaged zonal flow 
$\overline {U_\phi ^{0}}$, in the meridional plane at various 
$Ra$ but fixed 
$E=10^{-5}$. In all cases, the mean zonal flows are almost invariant along the rotation axis. However, the direction of the mean zonal flow reverses when the convective Rossby number is sufficiently large such that the LSV forms, i.e. 
$Ro_c \gtrsim 1.5$. Similar reversal of the zonal flow has been observed in thin shells for both Boussinesq convection (Aurnou et al. Reference Aurnou, Heimpel and Wicht2007) and anelastic convection (Gastine et al. Reference Gastine, Wicht and Aurnou2013). They also found that the reversal is controlled by the convective Rossby number 
$Ro_c$, or equivalently by 
$Ra^{*}$.
Figure 10. Mean zonal flows in the meridional plane at different 
$Ra$ for fixed 
$E=10^{-5}$: (a) 
$Ra=1.0\times 10^{3}$, 
$Ro_c=0.1$, in the relaxation oscillation regime; (b) 
$Ra=1.0\times 10^{5}$, 
$Ro_c=1.0$, in the GT regime; (c) 
$Ra=3.0\times 10^{5}$, 
$Ro_c=\sqrt {3}$, in the LSV regime.
As the mean zonal flow is essentially geostrophic, we use the zonal velocity profile in the equatorial plane to represent the zonal flow. Figure 11(a) shows the zonal profiles at various 
$E$ and 
$Ra$ in the GT regime. These profiles collapse into an invariant profile when the zonal velocity is divided by 
$Ro_c^{3/5}$. Note that the exponent 
$3/5$ is merely an empirical estimate and may be subject to small revisions. Nevertheless, it is clear that the zonal flow is controlled by the the convective Rossby number in the GT regime.
Figure 11. (a) Mean zonal flow on the surface as a function of the latitude at various 
$E$ and 
$Ra$, but all cases in the GT regime. (b) As for (a), but the mean zonal flow is rescaled by 
$Ro_c^{3/5}$. Solid (dashed) lines correspond to cases at 
$E=10^{-5}$ (
$E=3\times 10^{-5}$).
Figure 12(a) shows the zonal profiles in the LSV regime for various 
$E$ and 
$Ra$, but the same 
$Ro_c=\sqrt {3}$. The zonal flow profile in the LSV regime (figure 12a) is almost, but not exactly, the reverse of the profile in the GT regime (figure 11b). Aurnou et al. (Reference Aurnou, Heimpel and Wicht2007) proposed that the vigorous convection at large 
$Ra$ homogenises the angular momentum in the inertial frame, leading to a retrograde equatorial jet and a prograde zonal flow inside in the rotating frame. Figure 12(b) shows the specific angular momentum in the inertial frame, i.e. 
$M_I=s^{2}+s\overline {U_\phi ^{0}}$ (non-dimensional). We see an inward transfer of the angular momentum compared to solid body rotation, but the angular momentum profiles are far from a homogeneous distribution, even in the outer region. One may expect a more homogeneous distribution of the angular momentum on further increasing 
$Ra$, but the zonal flow appears to saturate in the LSV regime as we can see from figure 13(a).
Figure 12. (a) Mean zonal flow on the surface as a function of the latitude at various 
$E$ and 
$Ra$, but all cases in the LSV regime. (b) As for (a), but showing the specific angular momentum in the inertial frame.
Figure 13. Zonal Rossby number 
$Ro_{zon}$ (a) and rescaled non-zonal Rossby number 
$Ro_{non}/E^{2/5}$ (b) as a function of the convective Rossby number 
$Ro_c$. Different symbols represent different flow regimes as in figure 1. 
$Ro_{zon}$ and 
$Ro_{non}$ are time-averaged over several convective turnover time scales.
We have been unable to discover a principle that predicts this profile. Of interest is the fact that the angular momentum in the LSV regime is clearly prograde for 
$s\lessapprox 0.72$ and retrograde for 
$s\gtrapprox 0.72$, with the sum being zero, of course, in the rotating frame. But what determines this changeover point 
$s^{\ast }$? We merely remark that a fluid in solid body rotation contains half its angular momentum in 
$s \leqslant 0.77$ and half exterior to this radius. Any connection between 
$s^{\ast }$ and this value may be fortuitous.
The zonal Rossby number 
$Ro_{zon}$ increases with increasing 
$Ro_c$ with varying slopes depending on the flow regimes, and becomes saturated (
$Ro_{zon}\sim 0.15$) in the LSV regime (figure 13a). The plot is truncated to 
$Ro_{zon} \geqslant 10^{-3}$ to highlight the supercritical cases as 
$Ro_{zon}$ is very small near the onset. In the oscillatory regime 
$Ro_c<0.2$, 
$Ro_{zon}$ also depends on the Ekman number, suggesting that the viscosity plays a part in the zonal flow. Based on a balance between the Reynolds stress and the internal viscous stress, Christensen (Reference Christensen2002) proposed 
$Ro_{zon} \propto (Ro_{non})^{2}/E$. The non-zonal Rossby number as a function of 
$Ro_c$ is shown in figure 13(b), but the above relation does not hold for our numerical results as both 
$R_{zon}$ and 
$Ro_{non}$ follow the solid black line and roughly scale as 
$Ro_c^{2}$ in the oscillatory regime 
$Ro_c<0.2$.
In the GT regime 
$0.2 \leqslant Ro_c \leqslant 1.5$, 
$Ro_{zon}$ becomes independent of the viscosity and the empirical scaling 
$Ro_{zon}\sim Ro_c^{3/5}$ (black dashed line in figure 13a) is evident, in agreement with figure 11(b). In this regime, the non-zonal Rossby number is expected to follow the so-called inertial scaling, which is based on a force balance between Coriolis, inertia and Archimedean force, if the viscous effect is negligible (Aubert et al. Reference Aubert, Brito, Nataf, Cardin and Masson2001). The inertial scaling can be expressed as the following based on our definition of the control parameters (see appendix B):
\begin{equation} Ro_{non}\sim \left(\frac{E}{Pr}\right)^{2/5}\left(Ro_c\right)^{4/5}. \end{equation}
In figure 13(b), 
$Ro_{non}$ is divided by 
$E^{2/5}$ (
$Pr=1$) to compare with the inertial scaling (black dashed line). While the data points at different 
$E$ are well collapsed, we note that numerical results exhibit a steeper trend (larger exponent) than the inertial scaling. Previous numerical simulations in spherical shells also found steeper exponents (e.g. Gastine et al. Reference Gastine, Wicht and Aubert2016; Long et al. Reference Long, Mound, Davies and Tobias2020) and the deviation is usually attributed to the non-negligible effects of the viscosity in numerical simulations.
3.4. Heat transport
Figure 14 shows the Nusselt number as a function of the convective Rossby number 
$Ro_c$. Near the onset of the convection, the convective heat transfer is proportional to the supercriticality according to a weakly nonlinear theory (Gillet & Jones Reference Gillet and Jones2006), i.e.
\begin{equation} Nu-1 \propto \frac{Ra}{Ra_c}-1. \end{equation}
We can see from the lower left-hand corner of figure 14(a) that our numerical results can be well fitted by the aforementioned scaling when 
$Ra<6Ra_c$ for 
$E=10^{-5}$ (the blue solid line).
Figure 14. Nusselt number 
$Nu$ (a) and 
$NuE^{2/5}$ (b) as a function of the convective Rossby number 
$Ro_c$. Different symbols represent different flow regimes as in figure 1. The blue line depicting supercriticality is plotted in the lower left-hand corner of (a), for 
$Ro_c \leqslant 0.15$ and 
$E=10^{-5}$. Here 
$Nu$ is time-averaged over several convective turnover time scales.
In rotating turbulent convection, the inertial scaling has been promoted to describe heat transport (Stevenson Reference Stevenson1979; Aubert et al. Reference Aubert, Brito, Nataf, Cardin and Masson2001; Julien et al. Reference Julien, Knobloch, Rubio and Vasil2012a; Barker, Dempsey & Lithwick Reference Barker, Dempsey and Lithwick2014; Jones Reference Jones2015). Based on our definition of the non-dimensional parameters, this inertial scaling for the Nusselt number can be written as (appendix B)
\begin{equation} Nu \sim \left(\frac{E}{Pr}\right)^{{-}2/5}\left(Ro_c\right)^{6/5}. \end{equation} We can see that our numerical results in the GT regime (blue diamonds in figure 14) are in very good agreement with the inertial scaling. Note that we use the stress-free boundary condition in this study. Previous studies using the no-slip boundary found that the scaling exponent in 
$Ra$ depends on the Ekman number (Cheng et al. Reference Cheng, Stellmach, Ribeiro, Grannan, King and Aurnou2015; Gastine et al. Reference Gastine, Wicht and Aubert2016; Long et al. Reference Long, Mound, Davies and Tobias2020), i.e. 
$Nu\sim Ra^{\lambda (E)}$, though it may asymptotically approach the inertial scaling at low Ekman numbers (Gastine et al. Reference Gastine, Wicht and Aubert2016).
There is a sudden drop of 
$Nu$ on increasing 
$Ro_c$ when the flow transits from the GT regime to the LSV regime 
$Ro_c>1.5$. As we have mentioned, the reduction of the heat transport by the LSV is mainly due to the fact that a relatively hot column is formed around the rotation axis. Nusselt number 
$Nu$ continues to increase on further increase of 
$Ra$ in the LSV regime, but the slope seems to be less steep compared to the GT regime. It is difficult to deduce any definitive scaling tendency, as higher forcings would be required, and this is computationally prohibitive.
Figure 15(a) shows the heat flux 
$q$ on the surface in the LSV regime at 
$E=10^{-5}$ and 
$Ra=3.0\times 10^{5}$. The heat flux exhibits both small-scale perturbations and large-scale coherent variation in latitude, in line with the formation of LSV. The time and azimuthal average in figure 15(b) shows that the heat flux at high latitudes is about three times larger than that around the equator. The latitudal variations of the heat flux again suggest that the formation of LSV provides a barrier for the horizontal heat transport in the system.
Figure 15. Heat flux 
$q$ on the surface in the LSV regime at 
$E=10^{-5}$ and 
$Ra=3.0\times 10^{5}$. (a) Instantaneous heat flux on the surface (Mollweide projection). (b) Time and azimuthal average of the heat flux 
$q$ on the surface as a function of latitude.
3.5. Effects of boundary conditions
Thus far, we have made use of stress-free and fixed temperature boundary conditions. In this section, we briefly discuss effects of boundary conditions on the formation of LSV.
Previous studies (Guervilly et al. Reference Guervilly, Hughes and Jones2014; Stellmach et al. Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014), studying planar rotating convection, found that the LSV exists only when a stress-free boundary condition is imposed. In the present study, we briefly explore the possibility of formation of LSV with a no-slip boundary, but do not observe any LSV being formed, at least in the parameter range we could access. Figure 16 shows an example with a no-slip boundary condition at 
$E=10^{-5}$ and 
$Ra=5\times 10^{5}$, at which the LSV would form with the stress-free boundary condition. With the no-slip boundary condition, we see that the flow is characterised by small-scale structures elongated along the rotation axis (figure 16a,b). We note that a zonal flow is also developed with a prograde jet near the equator and a retrograde flow in the inner region (figure 16c), but the amplitude is much smaller compared to the zonal flow with a stress-free boundary condition. In fact, the zonal Rossby number is smaller than the non-zonal Rossby number when no-slip boundary conditions are employed (table 1), in contrast to the case with the stress-free boundary condition.
Figure 16. Axial vorticity and zonal flow with the no-slip boundary condition at 
$E=1.0\times 10^{-5}$ and 
$Ra=5.0\times 10^{5}$. (a,b) Axial vorticity in the equatorial plane and in the meridional plane. (c) Zonal flow in the meridional plane.
Figure 17 compares different thermal boundary conditions, i.e. fixed temperature (figure 17a,b) and fixed flux (figure 17c,d) conditions at 
$E=1.0\times 10^{-5}$ and 
$Ra=5.0\times 10^{5}$. Both cases use the stress-free boundary condition. It can be seen that in this regime the thermal boundary conditions play a minor role. A major difference is the temperature field in the polar regions, where the hot columnar structure associated with the LSV needs to adapt different thermal boundary conditions. The hot columnar region can extend to the polar surface when the fixed-flux condition is imposed.
Figure 17. Comparison of the fixed temperature boundary condition (a,b) and the fixed flux boundary condition (c,d) for 
$E=3.0\times 10^{-5}$ and 
$Ra=10^{5}$. (a,c) Axial vorticity in the equatorial plane and in the meridional plane. (b,d) Temperature in the equatorial plane and in the meridional plane.
4. Conclusions
We have observed the previously unseen presence of LSVs in full-sphere rotating convection, thus complementing the previous discoveries in Boussinesq rotating boxes (Favier et al. Reference Favier, Silvers and Proctor2014; Guervilly et al. Reference Guervilly, Hughes and Jones2014; Rubio et al. Reference Rubio, Julien, Knobloch and Weiss2014; Stellmach et al. Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014). Relatively simple large-scale zonal flow profiles are developed in this regime, with a prograde axial vortex and a retrograde equatorial zonal flow. While the geostrophic character of these flows is pervasive, the ubiquitous radial structure of the zonal flow in the LSV regime remains a theoretical challenge to comprehend. We have been able to fully document the transition between regimes as the Rayleigh number is increased from that required to first initiate convection, finding that the convective Rossby number is a completely prognostic control parameter: 
$Ro_c\approx 0.2$ and 
$Ro_c\approx 1.5$ delineate the boundaries between the oscillatory regime, the GT regime and the LSV regime.
In the GT regime, the zonal Rossby number 
$Ro_{zon}$ solely depends on the convective Rossby number and follows an empirical scaling law 
$Ro_{zon}\sim Ro_c^{3/5}$, while the non-zonal Rossby number and the Nusselt number can be described by the well-known inertial scaling. The direction of the zonal flow is reversed when the flow transitions from the GT regime to the LSV regime, while the non-zonal Rossby number and the Nusselt number are reduced by the formation of LSV. In the LSV regime we witness a saturation of the zonal flow speed, while the non-zonal Rossby number and the Nusselt number show an abrupt drop followed by a gradual increase, whose slopes are reduced compared to their systematics in the GT regime.
The effects of both the mechanical and thermal boundary conditions on the formation of LSVs are briefly considered. While the thermal boundary condition plays a minor role, the no-slip boundary condition prohibits the formation of LSVs, at least for the parameter regimes we could access, in line with previous studies in planar rotating convection (Guervilly et al. Reference Guervilly, Hughes and Jones2014; Stellmach et al. Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014). It is possible that the dynamics with the stress-free and the no-slip conditions would converge at very low Ekman numbers, but it is computationally demanding to study the highly supercritical convection at very low Ekman numbers via direct numerical simulations. We note the prospect of better experimental delineation of zonal flows in the future using the ZoRo experiment (Su et al. Reference Su, Cébron, Nataf, Cardin, Vidal, Solazzo and Do2020). This rapidly rotating spheroid has the capability of zonal flow determination using acoustic normal mode velocimetry, and will be subject to differential heating in the future.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2020.1151.
Acknowledgements
We thank P. Marti for the provision of the computer code EPMDynamo, and much implementation advice. Simulations were performed at the Swiss National Supercomputer Centre (CSCS) under account s872, on the Euler cluster of ETH Zurich and on the Taiyi cluster supported by the Center for Computational Science and Engineering of Southern University of Science and Technology.
Funding
The research was partially supported by a European Research Council Advanced grant 833848 (UEMHP) to A.J. under the European Union's Horizon 2020 research and innovation programme. Y.L. was supported by the B-type Strategic Priority Program of the Chinese Academy of Sciences (XDB41000000), the National Natural Science Foundation of China (grant no. 41904066), the pre-research project on Civil Aerospace Technologies of China National Space Administration (D020308) and the Macau Foundation. Y.L. was partly supported by a Swiss NSF Advanced PostDoc Mobility Fellowship when this study was initiated.
Declaration of interests
The authors report no conflict of interest.
Appendix A. List of numerical simulations
All of the numerical simulations are listed in table 1 including the control parameters, time-averaged diagnostic parameters, numerical resolutions (
$N,L,M$) and flow regimes. Here SD represents the steadily drifting regime near the onset; RO represents the relaxation oscillation regime; GT represents the geostrophic turbulence regime; LSV represents the large-scale vortices; NR represents the non-rotating regime. The case marked with # used the stress-free and fixed flux boundary condition, while the one marked with * used the no-slip and fixed temperature boundary conditions. All other cases used the stress-free and fixed temperature boundary conditions.
Appendix B. Inertial scaling expressed in our parameters
In rapidly rotating convection, the so-called inertial scaling is proposed to describe the convective flow speed, typical length scale and heat transport (Stevenson Reference Stevenson1979; Aubert et al. Reference Aubert, Brito, Nataf, Cardin and Masson2001). The inertial scaling is based on a balance of Coriolis–inertial–Archimedean terms, incorporated with the mixing length theory, leading to a scaling independent of the viscosity. Here we follow the review article of Jones (Reference Jones2015) and express the scaling in terms of dimensionless parameters defined in this study. Note that the Rayleigh number 
$Ra$ we defined is rotationally modified and based on heat flux, which is loosely related to 
$Ra_Q$ in Jones (Reference Jones2015) by 
$Ra\approx ERa_Q$. According to (119) in Jones (Reference Jones2015), the inertial scaling predicts the non-zonal Rossby number (representing the convective flow velocity) as
\begin{equation} Ro_{non}\sim \left(\frac{E}{Pr}\right)(EPr)^{1/5}Ra_Q^{2/5} \sim \left(\frac{E}{Pr}\right)^{4/5}Ra^{2/5}. \end{equation}
We can write the above scaling in the convective Rossby number 
$Ro_c=\sqrt {RaE/Pr}$ as
\begin{equation} Ro_{non}\sim \left(\frac{E}{Pr}\right)^{2/5}Ro_c^{4/5}. \end{equation}Note that the inertial scaling does not provide direct prediction of the zonal Rossby number.
The inertial scaling predicts the Nusselt number as ((119) in Jones (Reference Jones2015))
\begin{equation} Nu \sim E^{4/5}Pr^{{-}1/5}Ra_Q^{3/5} \sim \left(\frac{E}{Pr}\right)^{1/5}Ra^{3/5}. \end{equation}
Note that we have assumed 
$Nu-1 \sim Nu$ in the strongly supercritical regime. Using 
$Ra^{*}=RaE/Pr$, we obtain
\begin{equation} Nu \sim \left(\frac{E}{Pr}\right)^{{-}2/5}\left(Ra^{*}\right)^{3/5}. \end{equation}The above scaling is equivalent to
\begin{equation} Nu \sim E^{2} Pr^{{-}1/2} Ra_T^{3/2}, \end{equation}
where 
$Ra_T$ is the convectional Rayleigh number based on the temperature difference, i.e. 
$Ra_{T}=\alpha g_0 \Delta T r_o^{3}/{\nu \kappa }$, and is related to 
$Ra$ by
\begin{equation} Ra=2NuERa_{T}. \end{equation}
Equation (B4) can be also represented in terms of the convective Rossby number 
$Ro_c$ using 
$Ra^{*}=Ro_c^{2}$:
\begin{equation} Nu \sim \left(\frac{E}{Pr}\right)^{{-}2/5}\left(Ro_c\right)^{6/5}. \end{equation}
