1. Introduction
The spherical Couette system consists of two concentric spheres differentially rotating about a common axis, with the space in between filled with a viscous fluid. The differential rotation is considered ‘positive’ when the inner sphere rotates faster than the outer sphere and ‘negative’ when it rotates slower than or in the opposite direction to the outer sphere. Being the spherical analogue of the more well-known Taylor–Couette system (Chossat & Iooss Reference Chossat and Iooss1994), it is an interesting fluid dynamical system in its own right with very different instabilities. Applications to the interiors of astrophysical bodies (e.g. planetary interiors and stellar radiative zones) seem more obvious than in the Taylor–Couette geometry. The study of the spherical Couette system goes back to the analytical asymptotic formulation of Proudman (Reference Proudman1956) for an infinitely fast rotating outer sphere and an infinitesimal differential rotation. He showed that most of the fluid differential rotation remains confined within the cylinder tangent to the inner sphere equator, known as the tangent cylinder (TC), whereas the fluid outside the TC co-rotates with the outer boundary. A complex nested shear layer at the TC, known as the Stewartson layer (Stewartson Reference Stewartson1966), accommodates the jump in the fluid rotation rate and its derivatives. For a spherical Couette system with a wide gap, this shear layer is the source of the first flow instabilities for a rapidly rotating outer boundary. Note that when the gap becomes narrow, the flow instabilities resemble Taylor rolls similar to the Taylor–Couette system (Egbers & Rath Reference Egbers and Rath1995). Instabilities of a Stewartson layer driven by differential rotation were first studied experimentally by Hide & Titman (Reference Hide and Titman1967) for a cylindrical system with a differentially rotating disc and theoretically by Busse (Reference Busse1968). For the case of the spherical Couette system, Stewartson layer instabilities as well as other instabilities have been studied extensively using experiments (e.g. Sorokin, Khlebutin & Shaidurov Reference Sorokin, Khlebutin and Shaidurov1966; Munson & Menguturk Reference Munson and Menguturk1975; Egbers & Rath Reference Egbers and Rath1995; Kelley et al. Reference Kelley, Triana, Zimmerman, Tilgner and Lathrop2007, Reference Kelley, Triana, Zimmerman and Lathrop2010; Triana Reference Triana2011; Zimmerman et al. Reference Zimmerman, Triana, Nataf and Lathrop2014; Hoff, Harlander & Triana Reference Hoff, Harlander and Triana2016; Yoshikawa, Itano & Sugihara-Seki Reference Yoshikawa, Itano and Sugihara-Seki2023) and numerical computations (e.g. Munson & Joseph Reference Munson and Joseph1971a,Reference Munson and Josephb; Hollerbach Reference Hollerbach2003; Hollerbach et al. Reference Hollerbach, Futterer, More and Egbers2004; Hollerbach, Junk & Egbers Reference Hollerbach, Junk and Egbers2006; Wei & Hollerbach Reference Wei and Hollerbach2008; Matsui et al. Reference Matsui, Adams, Kelley, Triana, Zimmerman, Buffett and Lathrop2011; Rieutord et al. Reference Rieutord, Triana, Zimmerman and Lathrop2012; Wicht Reference Wicht2014). These studies have revealed a complex zoo of instabilities and have left many open questions.
Our previous study (Barik et al. Reference Barik, Triana, Hoff and Wicht2018, hereafter referred to as B18) and the present study are based on the experiments of Hoff et al. (Reference Hoff, Harlander and Triana2016) (hereafter referred to as H16) in a wide-gap spherical Couette set-up. Once the radius ratio of the two spheres is fixed, the system is characterised by two parameters, the Ekman number 
$E=\nu /\varOmega _o L^2$ and the differential rotation, 
$\Delta \varOmega /\varOmega = (\varOmega _i - \varOmega _o)/\varOmega _o$. Here, 
$\nu$ is the viscosity of the fluid, 
$L$ is the thickness of the spherical shell and 
$\varOmega _i$ and 
$\varOmega _o$ denote the rotation rates of the inner and outer sphere, respectively. H16 and B18 both focused on the case when the differential rotation was negative, i.e. when the inner sphere rotated slower than or in the opposite direction compared with the outer sphere. At intermediate or low Ekman numbers (
$3\times 10^{-6}\leq E \leq 10^{-4}$), as the differential rotation is made progressively more negative, the flow transitions through either four or five different hydrodynamic regimes.
(i) An axisymmetric flow described by Proudman (Reference Proudman1956).
(ii) The axisymmetric flow gives rise to a linear non-axisymmetric instability of the Stewartson shear layer with a fixed azimuthal wavenumber
$m$.(iii) The first instability gives way to a regime with a mode with
$m=1$. For a certain moderate to low range of Ekman numbers ($3\times 10^{-5}\leq E\leq 10^{-4}$), these two regimes may coincide and the first non-axisymmetric instability may occur in the form of $m=1$.(iv) The above regime gives way to equatorially antisymmetric (EA) wave-like ‘inertial modes’ which have been observed in several past studies (Kelley et al. Reference Kelley, Triana, Zimmerman, Tilgner and Lathrop2007, Reference Kelley, Triana, Zimmerman and Lathrop2010; Matsui et al. Reference Matsui, Adams, Kelley, Triana, Zimmerman, Buffett and Lathrop2011; Triana Reference Triana2011; Rieutord et al. Reference Rieutord, Triana, Zimmerman and Lathrop2012; Wicht Reference Wicht2014) and formed the focus of B18.
(v) Finally, a sharp and sudden transition to bulk turbulence takes place at a critical negative differential rotation.
These regimes have been observed in the simulations of Wicht (Reference Wicht2014) and B18 and the experiments of H16. More specifically, H16 observed that the transition to turbulence was characterised by a broadband temporal power spectra. Well-defined inertial mode peaks observed on top of these broadband spectra displayed an abrupt change in frequency right at the onset of turbulence. In addition, there was an increase in the spatial extent of the axisymmetric zonal flow and a decrease in the energy content of the inertial modes. They further observed a dependence of the critical differential rotation required for transition 
$|\Delta \varOmega /\varOmega |_c$ on the Ekman number as 
$|\Delta \varOmega /\varOmega |_c\sim E^{1/5}$. In the present study we concentrate on this transition to turbulence, addressing the following questions: how does the flow behave during and beyond the transition and what causes its onset?
There have been a few other studies on turbulence in spherical Couette flow, but not for the radius ratios and parameter ranges used in this study. Belyaev et al. (Reference Belyaev, Monakhov, Shcherbakov and Yavorskaya1979) experimentally analysed a wide-gap spherical Couette system for a stationary outer sphere and postulated that the transition to turbulence seems to follow the scenario of Ruelle & Takens (Reference Ruelle and Takens1971) but with several differences such as the existence of discrete peaks on top of a continuous background power spectrum. Yavorskaya & Belyaev (Reference Yavorskaya and Belyaev1986) experimentally investigated a thin-gap system (
$r_i/r_o = 0.9$) with both spheres rotating over a wide parameter range. They noted that the transition to turbulence involves an onset of ‘spatial intermittency’ in the form of small-scale structures on top of large-scale flow. Wulf, Egbers & Rath (Reference Wulf, Egbers and Rath1999) studied two different gap widths (
$r_i/r_o = 0.75$ and 
$0.67$) and found that the transition to turbulence was characterised by broadband temporal power spectra with some well-defined peaks.
The rest of the paper is arranged as follows. Section 2 provides details of the formulation of the problem, a brief description of the numerical methods used for simulation and the methods used to construct spectrograms and distinguish regimes (i)–(v) mentioned previously. Our results begin in § 3 with a discussion of the temporal and spatial spectra of flow and their variations. Section 4 discusses our results in the physical space with an analysis of the mean zonal flow, angular momentum transport and the effect of turbulence on inertial modes. Section 5 provides insight into the transition to turbulence by investigating force balances in the system. Section 6 investigates the instability of the boundary layer at the inner boundary and its effects and provides a heuristic explanation of the 
$E^{1/5}$ scaling law obtained by H16. Finally, § 7 discusses our main conclusions and open questions.
2. Numerical methods
2.1. Simulation set-up
Let us denote the radii and dimensional rotation rates of the two coaxial spheres as 
$r_i$ and 
$\varOmega _i$ for the inner sphere and 
$r_o$ and 
$\varOmega _o$ for the outer sphere, respectively. To simulate this system, we solve the Navier–Stokes and continuity equations in a reference frame rotating with the outer boundary. We use spherical coordinates 
$(r,\theta,\phi )$ denoting radial distance, colatitude and longitude, respectively. We also use 
$s=r\sin \theta$ to denote the cylindrical radius, perpendicular to the rotation axis. The equations are non-dimensionalised using 
$L=r_o-r_i$ as the length scale and the viscous diffusion time 
$\tau = L^2/\nu$ as the time scale, where 
$\nu$ is the kinematic viscosity of the fluid. This gives us,
$$\begin{gather} \frac{{\partial} \boldsymbol{u}}{{\partial} t} ={-} \boldsymbol{\nabla} p -\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} - \frac{2}{E}\hat{\boldsymbol{z}}\times \boldsymbol{u} + \boldsymbol{\nabla}^{2} \boldsymbol{u}, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0. \end{gather}$$
Here, 
$\boldsymbol {u}$ represents velocity, 
$p$ represents an effective pressure that includes centrifugal forces due to the outer boundary (system) rotation. The Ekman number 
$E = \nu /\varOmega _o L^2 = 1/\varOmega$, where 
$\varOmega$ is the non-dimensional outer boundary rotation rate. The inner sphere rotation rate (in the rotating frame) can also be similarly non-dimensionalised: 
$\Delta \varOmega = (\varOmega _i - \varOmega _o)L^2/\nu$. The system and the coordinate system is illustrated in figure 1.
Schematic of the spherical Couette system, indicating the rotation rates, the radii of the outer and inner boundaries and the spherical coordinate system 
$(r,\theta,\phi )$.
The flow problem is then characterised by three non-dimensional numbers: the Ekman number, 
$E$, the differential rotation 
$\Delta \varOmega /\varOmega$ and the radius ratio 
$\eta = r_i/r_o$, which is set to either 
$0.33$ or 
$0.35$. The former is the same as used in H16, whereas the latter is close to the ratio for Earth's core and has been used in B18 and other previous studies. No-slip boundary conditions allow for the viscous driving of the flow:
\begin{equation} \left.\begin{gathered} \boldsymbol{u}(r_o) = \boldsymbol{0},\\ \boldsymbol{u}(r_i) = (u_r,u_\theta,u_\phi) = (0,0,\Delta\varOmega r_i\sin\theta). \end{gathered}\right\} \end{equation}We numerically solve these equations using two independent pseudo-spectral codes: MagIC (Wicht Reference Wicht2002; Gastine et al. Reference Gastine, Barik, Rraynaud, Putigny, Thtassin, Johannes and Dintrans2021) (see https://github.com/magic-sph/magic) and XSHELLS (Figueroa et al. Reference Figueroa, Schaeffer, Nataf and Schmitt2013) (see https://bitbucket.org/nschaeff/xshells). The details of the numerical methods can be found in the respective publications. Both codes use the SHTns library (Schaeffer Reference Schaeffer2013) for spherical harmonic transforms.
As in H16 and B18, we concentrate on the case 
$\Delta \varOmega /\varOmega < 0$. The evolution of the flow is studied by keeping the outer boundary rotation (or Ekman number 
$E$) constant and running a simulation at a fixed 
$\Delta \varOmega /\varOmega$ and letting the kinetic energy reach a statistically stationary state. This state is then used as an initial condition to start the simulation for more negative 
$\Delta \varOmega /\varOmega$. The various parameters used in simulations and experiments along with the critical 
$\Delta \varOmega /\varOmega$ for transition to turbulence are listed in table 1, with each suite of experiments and simulations identified by a unique ID (first column). In B18 we already presented the simulation suites S1 and S3. In this study, we have run the rest of the suites, S2, S3a, S4 and S4a, with the suffix ‘a’ representing cases where the parameters are the same as the other case but the simulation is axisymmetric. The case S2 was run to confirm that numerical calculations yield the same critical 
$\Delta \varOmega /\varOmega$ for the turbulent transition as the experimental case E1. The ranges of 
$\Delta \varOmega /\varOmega$ in the table indicate how the differential rotation was changed within a suit, each using the previous simulation as starting condition (e.g. 
$-1.00$ to 
$-3.50$ means 
$\Delta \varOmega /\varOmega$ was made more negative starting from 
$-$1). This is important since the behaviour of the system has some amount of hysteresis (Egbers & Rath Reference Egbers and Rath1995). Through the rest of the paper, we mostly focus on simulation suites S3 and S4 with some comparisons with their axisymmetric counterparts S3a and S4a, respectively, and with experiments of H16 where appropriate. ‘Simulations’ will thus refer to simulations using MagIC unless otherwise specified. Figure 2 shows a diagram of the different regime transitions identified in simulations (filled circles) and experiments (open triangles). The suites S3 and S4 that are used throughout this paper clearly marked using squares (before transition to turbulence) and crosses (after transition). This does not show the suites S3a and S4a which would largely overlap with S3 and S4.
Parameters for different experiments at BTU C-S and simulations using MagIC and XSHELLS used in this study. Label ‘a’ denotes axisymmetric simulations and 
$\Delta \varOmega /\varOmega _c$ indicates the critical 
$\Delta \varOmega /\varOmega$ for transition to turbulence.
Physical regimes covered in H16, B18 and the present study. The horizontal axis shows Ekman number, 
$E$, whereas the vertical axis shows 
$|\Delta \varOmega /\varOmega |$. The open triangles and filled circles represent where various transitions to another regime have been found by experiments of H16 and simulations of B18 and Wicht (Reference Wicht2014), respectively. The brown squares and purple crosses show the simulation suites S3 and S4 from table 1, with squares indicating simulations in the EA inertial mode regime and crosses indicating simulations in the turbulent regime.
2.2. Spectrograms and identification of inertial modes
It has been shown in previous studies that inertial waves and modes are fundamental instabilities of the spherical Couette system. They obey the linear Euler equation,
\begin{equation} \frac{{\partial}\boldsymbol{u}}{{\partial} t} ={-}\boldsymbol{\nabla} p -2\varOmega\hat{\boldsymbol{z}}\times\boldsymbol{u}. \end{equation}This can be written as (see e.g. Greenspan Reference Greenspan1968; Tilgner Reference Tilgner2007)
\begin{equation} \frac{{\partial}^2}{{\partial} t^2}\nabla^2\boldsymbol{u} + 4\varOmega^2\frac{{\partial}^2}{{\partial} z^2}\boldsymbol{u} = 0, \end{equation}
which supports plane wave solutions (
$\boldsymbol {u}\propto \exp ({\mathrm {i}(\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol {r} - \omega t)})$), called ‘inertial waves’, in an unbounded fluid or bounded global oscillatory modes (
$\boldsymbol {u}\propto \exp ({\mathrm {i}(m\phi - \omega t)})$), called ‘inertial modes’, in a bounded container (Greenspan Reference Greenspan1968). In both cases, it can be shown that 
$|\omega |\leq 2\varOmega$ where 
$\omega$ is the frequency associated with a drift in azimuth 
$\phi$. Here, 
$\boldsymbol {k}$ and 
$m$ are the radial wavevector and the azimuthal wavenumber, respectively. For a spherical container, the solutions for inertial modes can be obtained analytically (Bryan Reference Bryan1889; Kudlick Reference Kudlick1966; Greenspan Reference Greenspan1968; Zhang et al. Reference Zhang, Earnshaw, Liao and Busse2001) and have the form of a spherical harmonic at the surface. Consequently, they are identified using indices 
$(l,m)$ corresponding to the spherical harmonic degree and order. These, together with the drift frequency 
$\omega$, uniquely determine a mode. Thus, as in our previous study (Barik et al. Reference Barik, Triana, Hoff and Wicht2018), we denote a mode using the notation 
$(l,m,\omega /\varOmega )$.
The different hydrodynamic regimes (i)–(v) mentioned in the introduction can be clearly identified using the spectrograms obtained from experimental data. The spectrograms are built by taking the single-sided fast Fourier transform (FFT) amplitude spectrum of the radially averaged azimuthal velocity 
$u_\phi$ at each 
$\Delta \varOmega /\varOmega$. The velocity measurements were performed via particle image velocimetry (PIV) techniques using a laser sheet perpendicular to the axis of rotation. The method is described in further detail in Hoff et al. (Reference Hoff, Harlander and Triana2016) and Hoff & Harlander (Reference Hoff and Harlander2019). Such spectrograms can also be constructed for simulations where we obtained data for the azimuthal component at eight different locations: 
$\theta = ({\rm \pi} /4,3{\rm \pi} /4)$ and 
$\phi = (0,{\rm \pi} /2,{\rm \pi},3{\rm \pi} /2)$, all on a single radial surface, 
$r = 0.7r_o$, which were stacked after correcting for their phase shift using cross-correlation of the time series. Thereafter, we performed a Fourier transform of this stacked time series to obtain spectra at each 
$\Delta \varOmega /\varOmega$ (see also Barik et al. Reference Barik, Triana, Hoff and Wicht2018). The spectrograms obtained from two suites of simulations are shown in figure 3(a,b), with identified inertial modes denoted by the indices 
$(l,m)$. Having access to the full three-dimensional (3-D) flow as well as a number of other diagnostics (kinetic energy, spatial spectra, etc.) in simulations helps distinguish these different regimes much better. For example, when the first non-axisymmetric 
$m=1$ mode appears, all the equatorially symmetric 
$m=1$ spherical harmonic flow coefficients can be seen to oscillate at the same frequency. An analysis of the spectral coefficients, the frequencies in a spectrogram, combined with a visualisation of the flow field is used to differentiate between the different regimes. The non-axisymmetric zonal flow fields in the three different regimes at 
$E=10^{-4}$ are illustrated in figure 4. Figure 4(a) shows the flow at 
$\Delta \varOmega /\varOmega =-1$ with the 
$m=1$ Stewartson layer instability (SI) clearly visible. Figure 4(b) shows the flow dominated by an EA (3,2) inertial mode with some small-scale features inside the TC, whereas figure 4(c) shows the flow in the turbulent regime at 
$\Delta \varOmega /\varOmega =-3$, with a lot of small-scale features near the inner boundary and a more chaotic flow field.
Spectrograms obtained from velocity time series, from simulations: (a) spectrogram from XSHELLS simulations at 
$E=1.125\times 10^{-4}$; (b) spectrogram from simulations at 
$E=3\times 10^{-5}$. Here 
$\omega$ is the angular frequency of the Fourier transform whereas 
$S(\omega )$ is the amplitude spectrum. The hydrodynamic regimes are marked and the inertial modes observed have been annotated. SI denotes Stewartson layer instability; EA denotes equatorially antisymmetric.
Isosurfaces of non-axisymmetric zonal flow from simulations at 
$E=10^{-4}$: (a) 
$\Delta \varOmega /\varOmega =-1$ is in regime (ii) with the 
$m=1$ mode; (b) 
$\Delta \varOmega /\varOmega =-2$ is in the regime with EA inertial modes; and (c) 
$\Delta \varOmega /\varOmega =-3$ is in the turbulent regime.
In the experiments of H16, the inertial modes were identified by comparing their frequencies against frequencies from theoretical works (Zhang et al. Reference Zhang, Earnshaw, Liao and Busse2001; Wicht Reference Wicht2014) as well as past experimental works (Kelley et al. Reference Kelley, Triana, Zimmerman, Tilgner and Lathrop2007, Reference Kelley, Triana, Zimmerman and Lathrop2010; Matsui et al. Reference Matsui, Adams, Kelley, Triana, Zimmerman, Buffett and Lathrop2011; Triana Reference Triana2011; Rieutord et al. Reference Rieutord, Triana, Zimmerman and Lathrop2012). Additional comparisons of morphology of modes were also made against theoretical inertial mode structures in spheres (Zhang et al. Reference Zhang, Earnshaw, Liao and Busse2001). In the case of simulations, the inertial modes can be clearly identified by a few different methods. The first is by comparing the frequencies observed in the spectrograms to the oscillation frequencies of the spectral spherical harmonic coefficients. This determines the longitudinal symmetry 
$m$ as well as the equatorial symmetry 
$(l-m)$ of the mode. The exact mode is then determined by comparing the frequency with the nearest analytical frequency of inertial modes in a sphere (Zhang et al. Reference Zhang, Earnshaw, Liao and Busse2001), as well as by spectrally filtering out the structure of the mode and comparing it with the theoretical structure.
3. Identifying transition to turbulence
In experiments as well as simulations, the temporal spectra help us determine the transition to turbulence. We examine here the spectra at individual 
$\Delta \varOmega /\varOmega$ values from the XSHELLS spectrogram presented in figure 3(a). We have selected three representative 
$\Delta \varOmega /\varOmega$ values, 
$\Delta \varOmega /\varOmega = (-0.6, -1.8, -2.7)$, which lie in regimes (iii), (iv) and (v), respectively, as shown in figure 5. At 
$\Delta \varOmega /\varOmega = -0.6$, the spectrum consists of only discrete peaks at the drift frequency of the 
$m=1$ SI and its higher multiples. In the EA inertial mode regime, at 
$\Delta \varOmega /\varOmega = -1.8$ (orange), there is a drastic change in the nature of the spectrum and it consists of a nearly flat background for 
$\omega /\varOmega \leq 2$ and a sharp decay for larger Fourier frequencies. The frequencies of the 
$m=1$ SI (around 
$\omega /\varOmega = 0.1$) and of the dominant inertial mode (
$(3,2)$ mode, around 
$\omega /\varOmega = 0.7$) are the most clearly visible peaks on top of the flat background. A flat background of energy for 
$0<\omega <2\varOmega$ and a subsequent decay demonstrates the fact that most of the kinetic energy in the flow manifests in inertial waves and is characteristic of inertial wave turbulence (e.g. Duran-Matute et al. Reference Duran-Matute, Flór, Godeferd and Jause-Labert2013; Clark di Leoni, Cobelli & Mininni Reference Clark di Leoni, Cobelli and Mininni2015). Thus, there is some amount of inertial wave turbulence already present in the EA inertial modes regime. This can be seen in the small scales visible inside the TC in figure 4(b), close to the inner boundary. However, the large-scale inertial mode still carries the dominant amount of energy in this regime.
Temporal spectra at three different values of differential rotation, each in a different hydrodynamic regime. The horizontal axis shows the Fourier frequency 
$\omega$ scaled with the outer boundary rotation rate 
$\varOmega$ whereas the vertical axis shows the amplitude spectrum 
$S(\omega )$. Vertical dotted lines mark 
$\omega /\varOmega = 0.1$, 0.7 and 2.
What we define as the ‘turbulent’ regime in this study is characterised by a further sudden increase in this flat background spectrum of inertial waves, as seen for 
$\Delta \varOmega /\varOmega = -2.7$, whereas the decay beyond 
$\omega /\varOmega = 2$ becomes less steep. Consequently, the peaks for the 
$m=1$ SI and the 
$(3,2,0.666)$ mode, despite having similar energies as for 
$\Delta \varOmega /\varOmega = -1.8$, are now less prominent with respect to the background. The small scale inertial wave turbulence is no longer limited to inside the TC, but now can be seen in the bulk as well (figure 4c) and, thus, the global large-scale inertial mode no longer carries a huge fraction of the energy. The typical decay of the spectrum beyond 
$2\varOmega$ has also been observed by H16 (figure 10 of Hoff et al. Reference Hoff, Harlander and Triana2016) and in the 3-meter experiment (figures 6.3 and 6.15 of Triana Reference Triana2011). The shallower decay of the spectrum beyond 
$2\varOmega$ in the turbulent regime shows a decrease in the influence of rotation which results in a greater content of energy for 
$\omega > 2\varOmega$. This is consistent with the fact that smaller scales and increased flow velocities in the turbulent regime lead to a dominance of advection over the effect of the Coriolis force. The change in the force balance is further discussed in § 5.
Pseudo-spectral codes provide direct information on different flow length scales and, hence, spatial spectra of kinetic energy. Figure 6 shows the change in energy spectrum in the zonal flow with 
$\Delta \varOmega /\varOmega$ at different colatitudes, very close to the inner boundary at 
$r/r_o = 0.354$. In figure 6(a), we see that for all 
$|\Delta \varOmega /\varOmega | \geq 1.5$ in the EA inertial mode regime, there already is a significant amount of energy in the smaller scales (high 
$m$) inside the TC. In figure 6(b), we find that the energy in the smaller scales are high for 
$|\Delta \varOmega /\varOmega | \geq 2.3$, indicating that the boundary layer at the inner boundary gets progressively destabilised at lower latitudes as 
$\Delta \varOmega /\varOmega$ becomes more negative. The turbulent regime sets in at 
$\Delta \varOmega /\varOmega = -2.3$ and its significance is that the boundary layer at the equator gets destabilised.
Change in the power spectrum of the zonal flow 
$u_\phi$ with 
$\Delta \varOmega /\varOmega$ at different latitudes at a fixed radius 
$r/r_o = 0.354$. The horizontal axis shows the azimuthal wavenumber 
$m$ whereas the vertical axis shows the power in a single wavenumber. Here 
$\Delta \varOmega /\varOmega =-2.3$ marking the transition to turbulence is plotted with a black line.
Figure 7 shows the total kinetic energy spectra with respect to spherical harmonic order 
$l$ at different radial levels from S3 simulations at 
$E=10^{-4}$. The onset of turbulence in the spatial spectra is characterised by an increase in energy in the small scales in general and close to the inner boundary in particular, which is the region of the highest flow speeds and thus most extreme Reynolds numbers. The system is driven by imposed differential rotation at the largest system scale. The energy then cascades to smaller scales via the different instabilities and nonlinear interactions. This cascade becomes decisively more efficient in the turbulent regime. The decrease in the influence of rotation can be seen in the spatial spectrum at large spherical harmonic degree as it gets progressively closer to a classic Kolmogorov 
$-$5/3 spectrum, as shown in figure 7(b).
Kinetic energy spectra from simulations at 
$E=10^{-4}$, shown at different radial levels: (a) the case for 
$\Delta \varOmega /\varOmega = -1.5$, before the transition to turbulence; (b) the case for 
$\Delta \varOmega /\varOmega = -3$, after the transition.
4. Flow analysis
4.1. Mean flow and angular momentum transport
The transition to turbulence is characterised by a sudden increase in the axisymmetric flow, whereas there is a drop in the non-axisymmetric kinetic energy (figure 8). The axisymmetric flow is clearly dominated by the zonal component, which is by a factor of 
$E^{-1/2}$ larger than the meridional circulation. Both components increase upon turbulence onset. The non-axisymmetric component is dominated by the EA part owing to the presence of the EA inertial modes before the transition to turbulence. This changes upon turbulence onset when EA inertial modes lose their energy. Figure 9(a,b) illustrates that the mean zonal flow not only intensifies but also starts to spread beyond the TC. The panels show the mean zonal flow, averaged in the 
$z$- and 
$\phi$-directions and in time, as a function of the cylindrical radial distance 
$s$ and the differential rotation rate for experiments E1 (a) and simulations S4 (b). Before the transition to turbulence (vertical lines), the zonal flow roughly resembles the Proudman solution for spherical Couette flow (Proudman Reference Proudman1956), staying restricted to the region inside the TC (horizontal lines). Beyond the transition, the zonal flow is significantly more vigorous and extends beyond the TC. The mean flow behaves the same way for simulations at 
$E=10^{-4}$. From table 1, we can compare the onset of turbulence for full 3-D simulations S3 and S4 and their axisymmetric counterparts, S3a and S4a. For 
$E=10^{-4}$, turbulence sets in at a 13 % lower differential rotation rate, at 
$E=3\times 10^{-5}$ the difference has reduced to 5 %. Figure 10 compares the time and azimuthally averaged zonal flow and meridional circulation in the turbulent regime at 
$E=10^{-4}$ of a 3-D simulation (a) and an axisymmetric simulation (b). Both cases show an additional pair of rolls along the inner boundary. These rolls represent an outward flow at the inner boundary which can contribute to advection towards the equator and then into the region outside the TC. In the axisymmetric case (figure 10b), the inner roll pair is more pronounced than in the 3-D case at the same parameters (figure 10a). In addition, an additional pair of rolls develops near the inner boundary equator and subsequently joins the set of rolls at higher latitudes to form a continuous radial jet. This jet-like feature is not seen in the 3-D case. The overall structure of the zonal flows is very similar in the two cases, axisymmetric turbulence is also characterised by a spreading of zonal flows beyond the TC.
Change in kinetic energy with differential rotation: (a) 
$E=10^{-4}$; (b) 
$E=3\times 10^{-5}$. As the flow transitions to turbulence (marked with vertical dotted lines), there is a sudden ‘burst’ in axisymmetric, mostly zonal, kinetic energy. The non-axisymmetric flow contributions, on the other hand, decrease in amplitude.
Plots of 
$z$, 
$\phi$ and time-averaged zonal flow velocity: (a) zonal flow from the experiments of H16 at 
$3.043\times 10^{-5}$; (b) the same from simulations at 
$E=3\times 10^{-5}$. The horizontal axis shows differential rotation whereas the vertical axis shows the cylindrical radial coordinate, scaled to the outer boundary 
$s/r_o$. The horizontal dotted line marks the TC whereas the vertical dotted line marks the critical differential rotation for the transition to turbulence.
Zonal flow and streamlines of meridional circulation from simulations at 
$E=10^{-4}$. Dashed (solid) lines represent anticlockwise (clockwise) circulation. (a) The case for a turbulent 3-D simulation at 
$\Delta \varOmega /\varOmega = -3.5$. (b) The same for an axisymmetric simulation at 
$\Delta \varOmega /\varOmega = -3.5$. Both are averaged in time and azimuth. Colours indicate zonal flow with blue being retrograde and red, if any, being prograde.
One can notice that the meridional circulation in figure 10(a) looks equatorially asymmetric compared with figure 10(b). Beyond the onset of the EA inertial modes, the 3-D simulations continue to possess an EA component of kinetic energy (figure 8), whereas the kinetic energy for the axisymmetric simulations continues to remain purely equatorially symmetric, even beyond the onset of turbulence. Thus, the symmetry breaking with respect to the equator seems unique to the presence of non-axisymmetric flows. The extension of the mean flow into the bulk along with the additional roll pair seems to push the Stewartson layer away from the TC. Whether we should still call this a Stewartson layer is unclear. As already discussed by Wicht (Reference Wicht2014), the appearance of a new pair of rolls close to the inner boundary indicates that the Coriolis force due to the outer boundary rotation ceases to be dominant. We expect this to happen when 
$\Delta \varOmega /\varOmega$ becomes smaller than 
$-1$. At 
$\Delta \varOmega /\varOmega =-1$, the inner boundary is at rest in the inertial frame. For 
$\Delta \varOmega /\varOmega \leq -1$ it rotates in the opposite direction to the outer boundary. When 
$\Delta \varOmega /\varOmega$ is negative enough, centrifugal forces drive an outward flow at the inner boundary that gives rise to the additional meridional roll pair. The transition from Coriolis-force-dominated dynamics to inertia-dominated dynamics should start close to the inner boundary where the effective rotation in the inertial frame is minimal.
Turbulence creates small-scale flow and transports angular momentum more efficiently from the inner boundary to the bulk of the flow outside the TC. This increases the viscous friction at the inner boundary so that a larger torque at the inner boundary is required to maintain the flow. Figure 11 shows the increase of the viscous torque on the inner core with 
$\Delta \varOmega /\varOmega$ for simulations at two different Ekman numbers. Before the onset of turbulence, the torque is simply proportional to 
$|\Delta \varOmega /\varOmega |$ and scales like 
$G\sim |\Delta \varOmega /\varOmega |^\alpha$ with 
$\alpha =1$, as shown with the compensated plot. In the turbulent regime, however, the torque increases more steeply with 
$\alpha \sim 2$. Zimmerman (Reference Zimmerman2010) reported approaching 
$\alpha =2$ in the 3-metre experiment in Maryland for a non-rotating outer boundary. The torque scalings for the axisymmetric simulations show a similar behaviour but the torque becomes smaller than in the 3-D simulations where non-axisymmetric instabilities provide a more efficient transport of angular momentum. In conclusion, the instability responsible for the onset of turbulence is predominantly an instability of the axisymmetric flow. The weaker non-axisymmetric flow components help stabilise the flow and yield a later onset.
Torque applied to the inner sphere and its variation with differential rotation magnitude, compensated by the linear scaling. The solid purple lines show a quadratic scaling.
4.2. Inertial modes
The large-scale EA inertial modes that get excited in regime (iv) continue to exist after the transition to turbulence. There is, however, a jump in the inertial mode frequencies. This can clearly be seen in the ‘brightest’ spectral lines in both panels of figure 3. This goes together with the sudden spreading of the background zonal flow beyond the TC causing further deformation of the inertial modes, as shown in B18. In both 3-D simulation suites that we studied, the flow is dominated by the inertial mode 
$(3,2,0.666)$ when turbulence sets in. After the transition, the mode loses at least half of its energy but still clearly dominates against the broadband turbulent background, as shown in figure 12 for experiments E1 and simulations S3. The energy estimates were determined by using a frequency filter on velocity obtained in experiments. In case of simulations, the energy in the large-scale spherical harmonic coefficients (order 
$l \leq 6$) of the equatorial and azimuthal symmetry corresponding to a mode was used to estimate the energy in a mode.
Kinetic energy of the different inertial modes: (a) data from experiments E1 at 
$E=3.043\times 10^{-5}$; (b) data from simulations S3 at 
$E=10^{-4}$. Both show the same features of the dominant inertial mode 
$(3,2)$ having a sudden drop in its kinetic energy. Vertical dotted lines mark the transition to turbulence in each case.
In both the numerical simulations S3 at 
$E=10^{-4}$ (MagIC) and S1 at 
$E=1.125\times 10^{-4}$ (XSHELLS), a new 
$m=2$ mode emerges around 
$\Delta \varOmega /\varOmega = -2.9$ with a frequency of 
$\omega /\varOmega \approx 0.4$. The mode is visualised at 
$\Delta \varOmega /\varOmega = -3$ in figure 13(a). We project snapshots of the flow velocity 
$\boldsymbol {u}$ and its non-axisymmetric part at different times onto equatorially symmetric inertial modes of a sphere 
$\boldsymbol {Q}_j \exp ({\textrm {i}(m\phi - \omega _j t)})$, similar to B18,
\begin{equation} \left.\begin{gathered} \boldsymbol{u} = \sum c_j \boldsymbol{Q}_j \exp({{\rm i}(m\phi - \omega_j t)}),\\ \boldsymbol{u} - \langle\boldsymbol{u}\rangle_\phi = \sum c_j' \boldsymbol{Q}_j \exp({{\rm i}(m\phi - \omega_j t)}). \end{gathered}\right\} \end{equation}A new 
$m=2$ mode emerging in the turbulent regime. (a) The 3-D side view of a snapshot from a MagIC simulation at 
$E=10^{-4}$, 
$\Delta \varOmega /\varOmega = -3$. The isosurfaces show non-axisymmetric zonal flow with red (blue) being positive (negative) at 
$u_\phi = \pm 500$. (b) Top view of the same. (c) The projection of the flow onto equatorially symmetric inertial modes.
The projection coefficients are normalised by 
$[\int \boldsymbol {u}\boldsymbol{\cdot }\boldsymbol {u} \,\textrm {d} V \int \boldsymbol {Q}_j\boldsymbol{\cdot }\boldsymbol {Q}_j^{{\dagger}} \,\textrm {d} V ]^{1/2}$ (or 
$[\int (\boldsymbol {u} - \langle \boldsymbol {u}\rangle _\phi ) \boldsymbol{\cdot }(\boldsymbol {u} - \langle \boldsymbol {u}\rangle _\phi ) \int \boldsymbol {Q}_j\boldsymbol{\cdot }\boldsymbol {Q}_j^{{\dagger}} \,\textrm {d} V]^{1/2}$ in the case of 
$c_j'$). The corresponding projection coefficients 
$c$ and 
$c'$, respectively, are shown in figure 13(b). It is clear that a single inertial mode cannot be used to characterise this flow structure, with dominant contributions from all modes with 
$m=2,\ l \leq 10$ that were analysed. We also could not find other modes that form triadic resonances with this mode.
5. Force balance
The transition to the turbulent regime for both S3 and S4 goes along with a sudden rise in the nonlinear term (
$\boldsymbol {u}\boldsymbol{\cdot }\boldsymbol {\nabla }\boldsymbol {u}$). As a consequence, advection rather than Coriolis becomes the dominant force. To understand the force balance at different length scales, we decompose the magnitude of each force 
$F$ into spherical harmonics,
\begin{equation} F(r) = \sum_{l=0}^{l_{max}}\sum_{m=-l}^l F_{lm}(r) Y_{lm}(\theta,\phi), \end{equation}
where 
$Y_{lm}(\theta,\phi )$ denotes a spherical harmonic of degree 
$l$ and order 
$m$. We then investigate the magnitude of forces, at different specific spherical harmonics degrees 
$l$ and radius levels, similar to Schwaiger, Gastine & Aubert (Reference Schwaiger, Gastine and Aubert2019),
\begin{equation} F_{rms}^2(l,r) = \frac{1}{V} \sum_{m=0}^{l} r^2 |F_{lm}(r)|^2, \end{equation}
where 
$V = 4/3 {\rm \pi}(r_o^3 - r_i^3)$ is the volume of the spherical shell. Figure 14 compares the respective spectra for two simulations at 
$E=10^{-4}$, one before the transition to turbulence (a,b) and one after (c,d). This is done for two different radial levels, one near the inner boundary and one in the bulk. At large scales (low 
$l$), the leading-order force balance near the inner boundary is dominated by advection and the Coriolis force whereas the dynamics in the bulk is determined by a geostrophic balance between the Coriolis force and the pressure gradient. This remains true for both before as well as after the transition to turbulence. At small scales (large 
$l$), after the transition to turbulence, there is a clear dominance of advection close to the inner boundary. This leads to the large-scale flow in the system being aligned with the rotation axis even in the turbulent regime, whereas small-scale flows dominate close to the inner boundary. This can be seen in the 3-D flow visualisation in figure 4(c) combined with the zonal flow visualised in figure 10.
Root-mean-square (RMS) spectra of different forces in the Navier–Stokes equation at two different radial levels, from simulations at 
$E=10^{-4}$: (a) 
$r/r_0=0.353$, (b) 
$r/r_0=0.834$, (c) 
$r/r_0=0.353$ and (d) 
$r/r_0=0.896$. Horizontal axes show spherical harmonic degree whereas vertical axes show RMS values. Panels (a,b) show the case at 
$\Delta \varOmega /\varOmega = -1.5$ in the laminar regime whereas panels (c,d) show the case at 
$\Delta \varOmega /\varOmega = -3$ in the turbulent regime. All plots have the same scale on the vertical axis.
We investigate the effect of the turbulent small scales on angular momentum transport using the azimuthal component of the Navier–Stokes equation. Separating the flow velocity and pressure into mean and fluctuating parts and a subsequent mean in azimuth and time gives us the Reynolds-averaged Navier–Stokes (RANS) equation for the mean zonal flow:
\begin{equation} -\frac{2}{E} \hat{\boldsymbol{\phi}}\boldsymbol{\cdot} \langle{\overline{\hat{\boldsymbol{z}}\times\boldsymbol{u}}}\rangle + \hat{\boldsymbol{\phi}}\boldsymbol{\cdot} \langle\overline{\nabla^2\boldsymbol{u}}\rangle - \hat{\boldsymbol{\phi}}\boldsymbol{\cdot} \langle\overline{\boldsymbol{\nabla} \boldsymbol{\cdot}\bar{\boldsymbol{u}}\bar{\boldsymbol{u}}}\rangle - \hat{\boldsymbol{\phi}}\boldsymbol{\cdot} \langle\overline{\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}'\boldsymbol{u}'}\rangle = 0, \end{equation}
where a bar denotes a mean in azimuth, 
$\langle \rangle$ denotes an average in time and prime denotes a non-axisymmetric part. We use about 700 snapshots of the 3-D flow at 
$\Delta \varOmega /\varOmega = -2$ in the inertial mode regime and about 1000 snapshots at 
$\Delta \varOmega /\varOmega = -3$ in the turbulent regime at 
$E=10^{-4}$, to compute the terms above, corresponding to 100 rotations of the outer boundary in each case. In addition, we use a time-averaged flow file for computing the terms for an axisymmetric case at 
$\Delta \varOmega /\varOmega = -3$, where the Reynolds stress 
$\langle \overline {\boldsymbol {\nabla }\boldsymbol{\cdot }\boldsymbol {u}'\boldsymbol {u}'}\rangle$ is absent. The results are shown in figure 15. In all cases, as expected, viscous forces are a dominant contributor to the zonal flow acceleration near the inner boundary. In the inertial mode regime (figure 15a), there is very little zonal flow generation and, hence, very little forcing outside the TC. Here the advection force 
$\langle \overline {\boldsymbol {\nabla }\boldsymbol{\cdot }\bar {\boldsymbol {u}}\bar {\boldsymbol {u}}}\rangle$ balances the viscous force near the equator and the Coriolis force away from the equator. Beyond the transition to turbulence (figure 15b), Reynolds stresses near the inner boundary balance the viscous force in this region, whereas the advective force provides the balance away from the equator. Slightly away from the boundary, the advective force balances the Coriolis force. In the axisymmetric turbulent case (figure 15c), in the absence of Reynolds stresses, the advective force balances both the Coriolis force as well as the viscous drag.
Different terms of the RANS equation (5.3) computed near the inner boundary for simulations at 
$E=10^{-4}$: (a) the case for 
$\Delta \varOmega /\varOmega = -2$, in the inertial mode regime; (b) the case for 
$\Delta \varOmega /\varOmega = -3$, in the turbulent regime; (c) a turbulent axisymmetric simulation at 
$\Delta \varOmega /\varOmega = -3$. The plots are clipped at a cylindrical radius of 
$s/r_o = 0.65$ and a vertical extent of 
$-$0.6 to 0.6.
In the case of 3-D turbulence, the small scales in the bulk of the fluid lead to Reynolds stresses that play the dominant role in forcing the zonal flow outside the TC, whereas in the axisymmetric case, the advection due to the strong radial jet plays the same role. The resultant efficient transport of angular momentum manifests itself in a change in torque scaling. Before the transition to turbulence, the zonal flow is restricted to the TC and its amplitude is linearly dependent on 
$\Delta \varOmega$, just as shown by Proudman (Reference Proudman1956). Thus, the torque on the inner sphere, 
$G = r_i \int \tau _{r\phi } \,\textrm {d} S = r_i\int {\partial }/ {\partial } r(u_\phi /r) \,\textrm {d} S$, where 
$\textrm {d} S = r_i\sin \theta \,\textrm {d}\theta \,\textrm {d}\phi$ is the differential surface area at the inner boundary, is also proportional to 
$\Delta \varOmega$. Beyond the transition to turbulence, the Reynolds stresses and the advective force contribute significantly to the zonal flow near the equator and become the major player in enhancing angular momentum transport. Their quadratic nature thus explains the quadratic scaling law in the turbulent regime.
6. Instability near the inner boundary
As noted in § 3, as 
$\Delta \varOmega /\varOmega$ becomes increasingly negative, the flow near the inner boundary first becomes unstable at high latitudes and gives rise to small-scale flows. At the transition to bulk turbulence, the flow near the inner boundary at and around the equator becomes unstable. This is illustrated in figure 16. Figure 16(a) shows radial velocity near the inner boundary before the transition to turbulence for suite S4 at 
$\Delta \varOmega /\varOmega = -1.98$ and figure 16(b) shows the same after the transition to turbulence at 
$\Delta \varOmega /\varOmega = -2$. The presence of small scales at all latitudes is markedly visible in figure 16(b). This is made even more clear if during this transition to turbulence, we track the radial velocity near the inner boundary at all latitudes and at a single longitude with respect to time. As illustrated in figure 17, when the fluid near the inner boundary spins up to 
$\Delta \varOmega /\varOmega = -2$, small-scale turbulent features start appearing near the equator. At the same time, the total and axisymmetric kinetic energies see a marked increase, as also explained in § 4.1.
Mollweide projection of radial velocity near the inner boundary at 
$r/r_o = 0.36$. Both simulations are for the suite S4 at 
$E=3\times 10^{-5}$: (a) the case for 
$\Delta \varOmega /\varOmega = -1.98$ before the transition to turbulence; (b) the case for 
$\Delta \varOmega /\varOmega = -2$, after the transition to turbulence.
Transition to turbulence through instability at the equator at 
$E=3\times 10^{-5}$, 
$\Delta \varOmega /\varOmega = -2$. Top panel shows radial velocity near the inner boundary 
$r/r_o = 0.36$ as a function of time (on the horizontal axis) and co-latitude 
$\theta$ on the vertical axis. Bottom panel shows the total, axisymmetric and non-axisymmetric kinetic energy as a function of time. The vertical dotted line marks the spin-up time based on 
$\varDelta\varOmega$.
In order to visualise the dynamics of these small scales, we also produced a movie using snapshots of the simulation suite S3 at 
$E=10^{-4}$, 
$\Delta \varOmega /\varOmega = -2.4$, in the turbulent regime. The initial condition was a solution at 
$\Delta \varOmega /\varOmega = -2.25$, without any boundary layer instability at the equator. Several snapshots at regular intervals were used to produce the movie, which is available as supplementary material (§ 7). The movie illustrates how small-scale structures of high angular momentum fluid emanate from the equatorial boundary layer and give rise to a mean flow. It also illustrates in an equatorial section how the zonal flow is close to being axisymmetric and large scale initially, destabilising soon after as it transitions into the turbulent regime. As discussed in § 4.1, a secondary pair of meridional circulation rolls are onset close to the inner boundary. However, the movie shows that their role is rather unimportant at this stage and the primary circulation is still responsible for most of the transport.
The Rayleigh stability criterion in a rotating frame is given by (Rayleigh Reference Rayleigh1917; Kloosterziel & van Heijst Reference Kloosterziel and van Heijst1991; Ghasemi et al. Reference Ghasemi, Klein, Harlander, Kurgansky, Schaller and Will2016)
\begin{equation} \varPhi = \frac{{\partial}}{{\partial} s}(u_\phi s + \varOmega s^2)^2 < 0, \end{equation}
where 
$\varOmega = 1/E$ is the outer boundary rotation rate. We use a couple of snapshots in time of the zonal flow at zero longitude to visualise the Rayleigh discriminant 
$\varPhi$ near the inner boundary. We use a simulation at 
$E=10^{-4}, \Delta \varOmega /\varOmega =-2.3$ for this purpose. Figure 18(a) shows the at the beginning of the simulation when turbulence has not yet set in, whereas figure 18(b) shows the case after 9.5 outer boundary rotations when the boundary layer is fully unstable at all latitudes. We can see that 
$\varPhi$ is strongly negative close to the inner boundary right at the start of the simulation. Though the colour bars are the same for both plots, in terms of actual values of 
$\varPhi$, at the beginning of the simulation (figure 18a), 
$\varPhi _{min} = -2.02\times 10^{10}$ and 
$\varPhi _{max} = 4.63\times 10^8$ whereas after 9.5 outer boundary rotations (figure 18b), 
$\varPhi _{min} = -4.42\times 10^{9}$ and 
$\varPhi _{max} = 7.28\times 10^8$. This implies that, in terms of extreme values of 
$\varPhi$, the fluid near the inner boundary is 4.6 times more stable and only about half as unstable in figure 18(a) compared with figure 18(b). The boundary between strongly negative and positive parts correlate quite well with the unstable regions in figure 18(b).
Visualising the Rayleigh stability criterion near the inner boundary for a couple of snapshots at 
$E=10^{-4}$ and 
$\Delta \varOmega /\varOmega = -2.3$: (a) the state at the beginning of the simulation; (b) after 9.5 rotations of the outer boundary (or 4 rotations of the inner boundary). The right panel in each case shows the zonal flow 
$u_\phi$ and the left panel shows the Rayleigh discriminant 
$\varPhi$. Here 
$N_{rot}$ denotes the number of outer boundary rotations.
We compute the thickness of the equatorial boundary layer at the inner boundary using a slope intersection method (similar to e.g. Verzicco & Camussi Reference Verzicco and Camussi1999; Gastine, Wicht & Aurnou Reference Gastine, Wicht and Aurnou2015; Barik et al. Reference Barik, Triana, Calkins, Stanley and Aurnou2023). We use time-averaged profiles of 
$\partial u_h/\partial r$, where 
$u_h = \sqrt {u_\theta ^2 + u_\phi ^2}$ is the magnitude of horizontal velocity. These profiles are obtained by averaging 
$u_h$ in azimuth and then in co-latitude with a window of 
$10^\circ$ centred at the equator. Fitting two lines to the respective profiles, one close to the inner boundary and a second one for the bulk, we assume that the boundary layer ends where both lines intersect. This is illustrated in figure 19. We then explore how the equatorial boundary layer thickness 
$\delta$ scales with differential rotation 
$|\Delta \varOmega /\varOmega |$. The thickness is compensated for by the theoretical 
$E^{2/5}$ scaling of the equatorial boundary layer thickness (Stewartson Reference Stewartson1966; Marcotte, Dormy & Soward Reference Marcotte, Dormy and Soward2016) in figure 20(a). We can see that the scaling works rather well, except for the axisymmetric suite S3a near the transition to turbulence. The equatorial boundary layer thickness increases very slowly with 
$|\Delta \varOmega /\varOmega |$ before the onset of turbulence. Close to the onset, there is an increase in the boundary layer thickness. After the transition to turbulence, the averaging in azimuth and time measures the thickness of the viscous sublayer, which is thinner than the laminar boundary layer and it decreases with Reynolds number (Landau & Lifshitz Reference Landau and Lifshitz1959; Grossmann & Lohse Reference Grossmann and Lohse2000). This is seen as a rapid decrease in 
$\delta$ beyond the transition and then a slow decrease with 
$|\Delta \varOmega /\varOmega |$. We can define a Reynolds number based on the boundary layer thickness,
\begin{equation} Re_\delta = \frac{(\varOmega_i - \varOmega_o) \delta^2}{\nu} = \frac{(\varOmega_i - \varOmega_o)}{\varOmega_o} \frac{\varOmega_o L^2}{\nu} \left(\frac{\delta}{L}\right)^2 = \frac{\Delta\varOmega}{\varOmega} \frac{1}{E} \left(\frac{\delta}{L}\right)^2. \end{equation}
If we assume that the boundary layer becomes turbulent once it exceeds a critical Reynolds number 
$Re_c$ and use the fact that 
$\delta /L = C E^{2/5}$, where 
$C$ is a constant, we find that at criticality,
\begin{equation} \left.\begin{gathered} C \left(\dfrac{\Delta\varOmega}{\varOmega}\right)_c \frac{1}{E} E^{4/5} = Re_c,\\ \Rightarrow \left(\frac{\Delta\varOmega}{\varOmega}\right)_c E^{{-}1/5} = Re_c/C. \end{gathered}\right\} \end{equation}
Figure 20(b) shows the variation of 
$Re_\delta$ with 
$|\Delta \varOmega /\varOmega |$ for all simulation suites. We find that, except for suite S3a, the rest of them peak at 
$Re_c = 42$, 
$42$ and 
$45$ for suites S3, S4 and S4a, respectively. This implies that the assumption the existence of a critical Reynolds number works fairly well. Furthermore, figure 21 shows the compensated plot of 
$|\Delta \varOmega /\varOmega |_c E^{-1/5}$ with data from both experiments of H16 as well as our simulations. We find that the spread in the compensated plot is small, especially noting the variation along the vertical axis. The higher Ekman number simulations are slightly off a flat line, with the axisymmetric suite S3a being a complete outlier.
Illustration of how the thickness of the boundary layer at the inner boundary is determined by slope intersection method, at 
$E=3\times 10^{-5}$, 
$\Delta \varOmega /\varOmega = -1.6$. The dashed lines show the slopes near the inner boundary and in the bulk, whereas the shaded region shows the boundary layer.
(a) Scaled thickness of the equatorial boundary layer 
$\delta$ as a function of 
$|\Delta \varOmega /\varOmega |$ for 3-D (solid lines, filled symbols) and axisymmetric simulations (dashed lines, open symbols). (b) Plot of 
$Re_\delta$ as a function of 
$\Delta \varOmega /\varOmega$ for 3-D simulation suites S3 and S4.
A compensated plot of 
$|\Delta \varOmega /\varOmega |_c E^{-1/5}$ vs Ekman number for the experimental data of H16 (open circles), our 3-D simulations (orange filled circles) and axisymmetric ones (green squares).
7. Conclusion
The two sets of simulations at Ekman numbers of 
$10^{-4}$ and 
$3\times 10^{-5}$ presented here yield similar results and reproduce the experimental observations of Hoff et al. (Reference Hoff, Harlander and Triana2016) in the turbulent regime. These include the generation of zonal flow in the bulk, violating the classic solution for spherical Couette flow by Proudman (Reference Proudman1956), the loss of energy in inertial modes and inertial wave turbulence. Unfortunately, the experimental data are extremely limited in spatial extent, being limited to a single plane perpendicular to the rotation axis just above the inner sphere. This makes it difficult to make more quantitative comparisons with experiments beyond what we already made in Barik et al. (Reference Barik, Triana, Hoff and Wicht2018) and in § 4.1. However, using simulations, we have been able to generate a more complete picture of the transition to turbulence. The cause of the onset of turbulence seems to be a centrifugal instability of the boundary layer at the equator of the inner boundary, giving rise to Taylor–Görtler vortices, similar to those observed by Noir et al. (Reference Noir, Hemmerlin, Wicht, Baca and Aurnou2009) and Ghasemi et al. (Reference Ghasemi, Klein, Harlander, Kurgansky, Schaller and Will2016, Reference Ghasemi, Klein, Will and Harlander2018). The hysteresis exhibited by the system (Egbers & Rath Reference Egbers and Rath1995) implies a subcritical transition. Beyond the regimes of axisymmetric flow and the first linear instability, as the differential rotation rate is made increasingly negative, we find that the boundary layer at the inner boundary first becomes unstable at high latitudes. This is seen in both in spectral space (§ 3) as well as physical space (§ 6). This instability gives place to spiral structures along the boundary layer, ejecting small-scale plume-like structures. These small-scale structures in a rotating environment are known to excite inertial waves (Davidson, Staplehurst & Dalziel Reference Davidson, Staplehurst and Dalziel2006; Davidson Reference Davidson2013), leading eventually to inertial wave turbulence as shown in § 3.
At a critical negative differential rotation, the boundary layer at the inner boundary becomes unstable at the equator and, thus, the resultant Görtler vortices can now propagate into the bulk, outside the TC. They further contribute to an increase in the energy carried by inertial waves as well as to an increase in energy in the small scales away from the inner boundary as evidenced by the temporal and spatial spectra (§ 3). A significant increase in Reynolds stresses driving zonal flow ensues, which leads to the zonal flow spreading outside the TC just as seen in experiments as well as simulations (§ 4.1). This also leads to a more efficient angular momentum transport and, thus, to an increase in the scaling exponent of the torque at the inner boundary from linear to quadratic. A second set of axisymmetric simulations at 
$E=10^{-4}$ and 
$3\times 10^{-5}$ show a very similar behaviour in terms of creation of large-scale zonal flow, torque scalings and destabilisation of the inner boundary layer near the equator. However, in this case the instability of the equatorial boundary layer at the inner boundary gives rise to an equatorial jet, which makes the subsequent evolution of the centrifugal instability markedly different than the 3-D simulations. This equatorial jet also serves to transport angular momentum in the axisymmetric cases as opposed to Reynolds stresses for the full 3-D simulations.
The Ekman layer near the inner boundary merges with the Stewartson layer into a layer that has an extent of 
$\delta _s \times \delta _z = E^{2/5}\times E^{1/5}$ (Stewartson Reference Stewartson1966; Marcotte et al. Reference Marcotte, Dormy and Soward2016) with 
$s$ and 
$z$ representing the cylindrical radius and axial direction, respectively. Using a heuristic critical Reynolds number argument for the destabilisation of the equatorial boundary layer at the inner boundary, we show that this scaling can help explain the experimental 
$E^{1/5}$ scaling for the critical differential rotation, especially at lower Ekman numbers (§ 6). The finer details of this transition are something that can still be explored and investigated, but the centrifugal instability of the equatorial boundary layer is the clear precursor. It remains to be seen whether this scaling law extends to asymptotically low Ekman numbers. The narrowing gaps between the values of 
$(\Delta \varOmega /\varOmega )_c$ for the full 3-D and the axisymmetric simulations and their similar nature of instabilities is encouraging. This could enable us to obtain an estimate of 
$(\Delta \varOmega /\varOmega )_c$ at lower Ekman numbers with cheaper axisymmetric computations.
Furthermore, our previous work (Barik et al. Reference Barik, Triana, Hoff and Wicht2018) and current study have been limited to 
$\Delta \varOmega /\varOmega < 0$. An in-depth study for 
$\Delta \varOmega /\varOmega > 0$ is still lacking. In particular, it is not clear why one obtains high-wavenumber spiral Stewartson layer instabilities for 
$\Delta \varOmega /\varOmega > 0$, but low-wavenumber instabilities trapped inside TC for 
$\Delta \varOmega /\varOmega < 0$ (Hollerbach Reference Hollerbach2003; Wicht Reference Wicht2014). More simulations and experiments are needed to establish better scaling laws pertaining to the different hydrodynamic regimes at lower Ekman numbers. This will prove helpful not only in order to extrapolate the behaviour of the spherical Couette system to real objects, but also to understand the dichotomies between 
$\Delta \varOmega /\varOmega < 0$ and 
$\Delta \varOmega /\varOmega > 0$. The theoretical foundation for spherical Couette flow is still in its infancy as compared with the more traditional Taylor–Couette system (Grossmann, Lohse & Sun Reference Grossmann, Lohse and Sun2016). Our present study shows that there is a great scope for similar studies in spherical shells as well, where the presence of spherical curvature makes the problem less tractable.
Supplementary movies
The movie for transition to turbulence can be found at doi.org/10.6084/m9.figshare.9108533. The full 4k version can be viewed as an unlisted youtube video: https://youtu.be/6vBWwYIapC8.
Acknowledgements
A.B. would like to thank A. Tilgner, J. Aurnou, N. Schaeffer and P. Wulff for insightful discussions and S. Stanley for feedback on the manuscript draft. We gratefully acknowledge A. Mazilu from the Transylvania University of Brasov (Romania) for performing most of the experiments in the frame of a Traineeship ERASMUS+ program.
Funding
A.B. would like to thank the IMPRS for Solar System Science for funding him during his time in Germany and Sabine Stanley for funding him subsequently. A.B. would also like to thank the North-German Supercomputing Alliance (HLRN), Max Planck Computing and Data Facility (MPCDF) and the Gesellschaft für wissenschaftliche Datenverarbeitung mbH Göttingen (GWDG) for letting him generously use their supercomputing facilities. S.A.T. would like to acknowledge support from the European Research Council (ERC Advanced Grant 670874 ROTANUT) and from the infrastructure program EuHIT of the European Commission. At BTU the experiment and the PhD position of M.H. was financed by a DLR grant (50 WM 0822) and by grants from the German science foundation DFG (HA 2932/7-1, HA 2937/8-2, HA 2937/10-1).
Declaration of interests
The authors report no conflict of interest.
Data availability statement
Both codes used to run the simulations listed here are openly available. MagIC is available at https://github.com/magic-sph/magic and XSHELLS at https://bitbucket.org/nschaeff/xshells. The parameters used to run the codes are available in the paper and in Barik et al. (Reference Barik, Triana, Hoff and Wicht2018).
Open access
