2.1 The spherical Couette system

The spherical Couette system is illustrated in figure 1. It consists of concentric inner and outer spheres with radii $r_{i}$ and $r_{o}$ , respectively, with a viscous fluid in between. We use a frame of reference rotating with the outer boundary. The time evolution of the flow in this system can be expressed in terms of the Navier–Stokes and continuity equations,

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}=-\frac{1}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}p-\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}-2\unicode[STIX]{x1D6FA}_{o}\hat{\boldsymbol{z}}\times \boldsymbol{u}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0. & \displaystyle\end{eqnarray}$$

Here, $\boldsymbol{u}$ , $\unicode[STIX]{x1D70C}$ and $\unicode[STIX]{x1D6FA}_{o}$ represent the velocity of the fluid, the density of the fluid and rotation rate of the outer boundary, respectively. The unit vector in the direction of rotation is denoted by $\hat{\boldsymbol{z}}$ , while $p$ denotes an effective pressure, which includes the contribution of gravity and centrifugal forces.

Using the shell width, $L=r_{o}-r_{i}$ , as a length scale and the viscous diffusion time, $L^{2}/\unicode[STIX]{x1D708}$ , as a time scale, we define the Ekman number as

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle E=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D6FA}_{o}L^{2}=1/\unicode[STIX]{x1D6FA}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FA}_{o}$ and $\unicode[STIX]{x1D6FA}$ are the dimensional and non-dimensional outer boundary rotation rates respectively, and $\unicode[STIX]{x1D708}$ is the kinematic viscosity.

The differential rotation is quantified using

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=(\unicode[STIX]{x1D6FA}_{i}-\unicode[STIX]{x1D6FA}_{o})/\unicode[STIX]{x1D6FA}_{o}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FA}_{i}$ is the dimensional inner boundary rotation rate and $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}$ is the non-dimensional differential rotation rate of the inner boundary with respect to the outer. Thus, $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}<0$ would imply that the inner sphere rotates slower or in the opposite direction with respect to the outer sphere, while $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}>0$ would mean that the inner sphere rotates faster than the outer sphere in the same direction.

The equations are completed with the impenetrable and no-slip boundary conditions:

(2.5) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \boldsymbol{u}(r_{o})=\mathbf{0},\\ \displaystyle \boldsymbol{u}(r_{i})=(u_{r},u_{\unicode[STIX]{x1D703}},u_{\unicode[STIX]{x1D719}})=(0,0,(\unicode[STIX]{x1D6FA}_{i}-\unicode[STIX]{x1D6FA}_{o})r_{i}).\end{array}\right\} & & \displaystyle\end{eqnarray}$$

2.2 Inertial modes

Inertial modes are eigenmodes of the linear inviscid Navier–Stokes equation,

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}=-\displaystyle \frac{1}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}p-2\unicode[STIX]{x1D6FA}\hat{\boldsymbol{z}}\times \boldsymbol{u}, & \displaystyle\end{eqnarray}$$

in a bounded container subject to the impenetrable boundary condition $\boldsymbol{u}\boldsymbol{\cdot }\hat{\boldsymbol{n}}=0$ , where $\unicode[STIX]{x1D6FA}$ represents the rotation rate of the container, and $\hat{\boldsymbol{n}}$ is the unit vector normal to the surface of the container.

Taking the curl of (2.6) twice and performing some algebra yields

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}}{\unicode[STIX]{x2202}t^{2}}+4\unicode[STIX]{x1D6FA}^{2}\frac{\unicode[STIX]{x2202}^{2}\boldsymbol{u}}{\unicode[STIX]{x2202}z^{2}}=0, & \displaystyle\end{eqnarray}$$

where $z$ denotes the coordinate along the axis of rotation. This equation supports wave-like modes

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=\boldsymbol{Q}(r,\unicode[STIX]{x1D703})\text{e}^{\text{i}(m\unicode[STIX]{x1D719}-\unicode[STIX]{x1D714}t)}=\boldsymbol{Q}(r,\unicode[STIX]{x1D703})\text{e}^{\text{i}m(\unicode[STIX]{x1D719}-c_{d}t)}, & \displaystyle\end{eqnarray}$$

where we have used spherical coordinates. Here, $m$ represents the azimuthal wavenumber, $\unicode[STIX]{x1D714}$ is the drift frequency and $c_{d}=\unicode[STIX]{x1D714}/m$ is the drift speed. The eigenvalues $\unicode[STIX]{x1D714}$ are real and $|\unicode[STIX]{x1D714}|\leqslant 2\unicode[STIX]{x1D6FA}$ (Greenspan Reference Greenspan1968). We use the convention of denoting retrograde (prograde) modes with positive (negative) frequencies. The inertial mode problem is mathematically ill-posed in general (Stewartson & Rickard Reference Stewartson and Rickard1969), but smooth solutions can be obtained for certain containers, like cylinders and spheres. These modes are also known to exist in slightly modified form in viscous systems when rotational effects dominate viscous effects (e.g. Kelley et al. Reference Kelley, Triana, Zimmerman, Tilgner and Lathrop2007; Triana Reference Triana2011; Wicht Reference Wicht2014).

In the following, we shall compare the observed inertial modes with those of a sphere whose modes have been studied for a long time (Bryan Reference Bryan1889) and for which explicit analytical solutions have been found by Zhang et al. (Reference Zhang, Earnshaw, Liao and Busse2001). For a spherical shell, in the absence of analytical solutions (Stewartson & Rickard Reference Stewartson and Rickard1969), one turns towards numerical eigenmode computations (e.g. Rieutord Reference Rieutord1991). For the present study we shall be using the eigenmodes computed using the code SINGE (Vidal & Schaeffer Reference Vidal and Schaeffer2015) for comparisons.

The inertial modes of a sphere can be characterised using three indices $(l,m,n)$ (indices $(n,k,m)$ of Greenspan (Reference Greenspan1968)), where $l$ and $m$ have the same meaning as the degree and order, respectively, of a spherical harmonic, and the third index denotes the $n$ th eigenmode corresponding to the combination $(l,m)$ . However, in our present study we identify a mode by using the indices $(l,m)$ and its frequency $\unicode[STIX]{x1D714}/\unicode[STIX]{x1D6FA}$ . The indices refer to the spherical harmonic structure at the surface, $r=r_{o}$ . At depth, the structure and thus the spherical harmonic contribution changes, as is illustrated in figure 2 for the mode $(5,2,0.467)$ .

Figure 2.

The $(5,2,0.467)$ mode of a full sphere. The mode has the structure of the spherical harmonic $Y_{52}$ at the surface while that of $Y_{32}$ at mid-depth. Colours represent cylindrical radial velocity $u_{s}$ , with red (blue) being away from (towards) the rotation axis. The two-dimensional maps are hammer projections of the same quantity.

Inertial modes in a full sphere can be equatorially symmetric or antisymmetric. The equatorially symmetric modes that we find in simulations are the slow inertial modes described by Zhang et al. (Reference Zhang, Earnshaw, Liao and Busse2001) which are ‘nearly geostrophic’. This implies that they satisfy the Taylor–Proudman constraint to a large extent, i.e. the pressure gradient and Coriolis forces roughly balance each other,

(2.9) $$\begin{eqnarray}\displaystyle & \displaystyle -\text{i}\unicode[STIX]{x1D714}\boldsymbol{u}=-\unicode[STIX]{x1D735}p/\unicode[STIX]{x1D70C}-2\unicode[STIX]{x1D6FA}\hat{\boldsymbol{z}}\times \boldsymbol{u}\approx 0. & \displaystyle\end{eqnarray}$$

In such a scenario, the frequencies of these modes, $\unicode[STIX]{x1D714}$ , are ‘small’ in magnitude as compared to the boundary rotation $\unicode[STIX]{x1D6FA}$ and can thus change by a large fraction even for a small mismatch between pressure gradient and the Coriolis force. This is in contrast to the equatorially antisymmetric modes for which no such slow modes are found and which completely defy the Taylor–Proudman constraint and have larger magnitudes of frequencies as compared to the slow, nearly geostrophic, modes. Owing to this clear distinction between the frequencies of these two classes of modes, the equatorially antisymmetric modes were referred to as ‘fast inertial modes’ by Wicht (Reference Wicht2014). Thus, while the fast modes are easily identified by their frequencies, the same cannot be said of the slow columnar modes whose frequencies are susceptible to the slightest of changes in the force balance.

4.1 The numerical set-up

Equations (2.1) and (2.2) with the boundary conditions (2.5) are solved using two independent codes: MagIC (Wicht Reference Wicht2002) (see https://github.com/magic-sph/magic) and XSHELLS (Figueroa et al. Reference Figueroa, Schaeffer, Nataf and Schmitt2013) (see https://bitbucket.org/nschaeff/xshells). All quantities are non-dimensionalised using the scheme mentioned in § 2. Both codes use a poloidal/toroidal decomposition of the velocity field,

(4.1) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}=\unicode[STIX]{x1D735}\times \unicode[STIX]{x1D735}\times \hat{\boldsymbol{r}}W+\unicode[STIX]{x1D735}\times \hat{\boldsymbol{r}}Z, & \displaystyle\end{eqnarray}$$

where $W$ and $Z$ are the poloidal and toroidal potentials, respectively. Spherical harmonics are used for the spectral representation of these scalar potentials and pressure in latitude and longitude; for example,

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle W(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719},t)=\mathop{\sum }_{l=0}^{l_{max}}\mathop{\sum }_{m=0}^{l}W_{lm}(r,t)Y_{lm}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719}). & \displaystyle\end{eqnarray}$$

The codes differ in the radial representation of the variables $W$ , $Z$ and $p$ . In MagIC, they are expanded in terms of Chebyshev polynomials and the radial grid consists of Gauss–Lobatto nodes. XSHELLS uses second-order finite differences in radius.

The reader is encouraged to have a look at Christensen & Wicht (Reference Christensen, Wicht and Schubert2007) and Figueroa et al. (Reference Figueroa, Schaeffer, Nataf and Schmitt2013) for further details on the numerical methods. Some of the simulations using MagIC and all the simulations using XSHELLS presented in this paper used the SHTns library (Schaeffer Reference Schaeffer2013) for spherical harmonic transforms.

The reason for presenting complementing results from two different codes using slightly different parameters is purely historical. The MagIC simulations were performed for a radius ratio describing the Earth’s core geometry while the XSHELLS simulations used the ratio appropriate for the BTU C-S experiment. The decision to publish a joint paper only came later. Both codes have been shown to give identical results for benchmark simulations (Matsui et al. Reference Matsui, Heien, Aubert, Aurnou, Avery, Brown, Buffett, Busse, Christensen and Davies2016).

4.2 Identifying inertial modes

The identified inertial modes in experiments and simulations of the spherical Couette system are fast modes. Their well-defined frequencies and structures help to identify them in experiments and make them suitable for comparisons between simulations and experiments.

One of the advantages of pseudo-spectral simulations is the direct access to the individual poloidal and toroidal coefficients, $W_{lm}$ and $Z_{lm}$ . A careful study of the behaviour of these coefficients provides an easy way of identifying the presence of fast inertial modes in the system. The equatorial symmetry of the flow field can easily be related to that of these coefficients, as shown in table 1. An equatorially antisymmetric (EA) inertial mode of azimuthal wavenumber $m_{0}$ would manifest itself as oscillations of spherical harmonic coefficients of the corresponding symmetry, i.e. $W_{l_{1}m_{0}}$ and $Z_{l_{2}m_{0}}$ , where $l_{1}+m_{0}=\text{odd}$ and $l_{2}+m_{0}=\text{even}$ . A Fourier analysis of the oscillations of these coefficients would provide the drift frequency of the inertial mode.

Figure 4.

(a) Oscillations of spherical harmonic coefficients of toroidal ( $Z(\ell ,m)$ ) and poloidal ( $W(\ell ,m)$ ) potentials; $Z(2,2)$ and $W(3,2)$ are the dominant oscillatory coefficients. (b) Fourier transform of time series of coefficients $Z(2,2)$ and $W(3,2)$ , denoted by tilde, showing a peak in frequency at the same position, $\unicode[STIX]{x1D714}/\unicode[STIX]{x1D6FA}=0.72$ . The parameters are $E=10^{-4}$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-2$ .

Table 1.

Symmetry of flow field and the corresponding conditions on degree and order of spherical harmonics of poloidal and toroidal potentials.

As an example, let us take the $(3,2,0.666)$ inertial mode identified by Wicht (Reference Wicht2014). Figure 4 shows that the toroidal coefficient, $Z_{22}$ ( $l+m=\text{even}$ ) is the dominant oscillatory coefficient, while the same is true for the poloidal coefficient $W_{32}$ ( $l+m=\text{odd}$ ). The Fourier analyses of the coefficients $Z_{22}$ and $W_{32}$ show a single strong peak at the same frequency of $\unicode[STIX]{x1D714}/\unicode[STIX]{x1D6FA}=0.72$ , where $\unicode[STIX]{x1D714}$ is the angular frequency of the Fourier transform. This is shown in figure 4(b).

Table 2.

List of the most dominant inertial modes observed and identified, and comparisons with experiments. For simulations at $E=10^{-5}$ , the modes seem to occur at similar ranges of $|\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}|$ as all the experiments spanning a huge range of Ekman numbers.

a Using SINGE (Vidal & Schaeffer Reference Vidal and Schaeffer2015).

b From Triana (Reference Triana2011).

c Frequency ranges from Kelley et al. (Reference Kelley, Triana, Zimmerman, Tilgner and Lathrop2007); $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}$ values are single representative values from S. Triana.

d From Hoff et al. (Reference Hoff, Harlander, Egbers, Triana, Egbers, Ruck, Leder and Dopheide2016a ,Reference Hoff, Harlander and Triana b ).

e Frequencies determined using spherical harmonics.

f Frequencies determined from spectrograms.

Table 3.

List of less dominant inertial modes and comparisons with experiments and theory. The modes have been identified based on frequency and structures.

a Sources of data same as in table 2.

Tables 2 and 3 give a list of the ‘fast’ inertial modes identified in our simulations and compare it with experimental data as well as data obtained from analytical predictions for a full sphere (e.g. Greenspan Reference Greenspan1968; Zhang et al. Reference Zhang, Earnshaw, Liao and Busse2001) and eigenmode computations using SINGE. The frequency ranges provided in the tables reflect the fact that the drift frequencies of the modes are not constant but change with differential rotation ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}$ ) and outer boundary rotation ( $\unicode[STIX]{x1D6FA}=1/E$ ). Differences between MagIC and XSHELLS can be traced back to the slightly different parameters. In some instances the numerical results report an earlier onset of inertial modes, which reflects the easier mode identification using spherical harmonic coefficients.

5.1 Fundamental $m=1$ instability

The first instability in the system is a Stewartson layer instability (SI) (e.g. Hollerbach Reference Hollerbach2003) with a well-defined critical $(\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA})_{c}$ and azimuthal symmetry $m_{c}$ . At $E=10^{-4}$ , the first instability has an $m=1$ and drifts in azimuth with a frequency $\unicode[STIX]{x1D714}/\unicode[STIX]{x1D6FA}=0.086$ (Wicht Reference Wicht2014). This fundamental instability is nearly geostrophic and thus drifts with a very low drift frequency (see § 2.2).

Figure 6.

Critical values for onset of first non-axisymmetric instability at different Ekman numbers: (a) the critical wavenumber $m$ ; (b) the critical differential rotation. Red circles in (b) mark the onset of the fundamental $m=1$ mode seen as a secondary instability at low Ekman numbers. A linear extrapolation shows that the 3 m experiment was always supercritical to the $m=1$ onset. Green crosses mark the critical $|\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}|$ values determined in experiments at BTU C-S.

Figure 7.

The $m=1$ fundamental mode at (a) $E=10^{-4}$ , (b)  $E=10^{-5}$ and (c) $E=3\times 10^{-6}$ . At $E=10^{-4}$ , the mode is the first non-axisymmetric instability to set in. At lower Ekman numbers, it comes in as a secondary instability. Shown here are isosurfaces of non-axisymmetric zonal velocity, $u_{\unicode[STIX]{x1D719}}$ , with red (blue) being positive (negative).

At $E=10^{-5}$ , even though the first SI has an $m_{c}=2$ and $(\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA})_{c}=-0.133$ , the same fundamental $m=1$ instability sets in as a secondary instability already at $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-0.16$ . Numerical exploration at a lower Ekman number of $3\times 10^{-6}$ and the experimental results of Triana (Reference Triana2011) confirm that this mode continues to set in as a secondary instability for smaller Ekman numbers and negative differential rotation. Figure 6(a) shows the critical wavenumber and figure 6(b) the critical differential rotation at onset of the first non-axisymmetric instability for different Ekman numbers for $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}<0$ . In addition, the onset of the secondary $m=1$ instability has also been shown. A linear fit has been used to extrapolate the onset of this instability to lower Ekman numbers. For example, at $E=3\times 10^{-6}$ , the first instability sets in for $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-0.0213$ and has a wavenumber of $m_{c}=12$ while a secondary $m=1$ instability dominates at $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-0.09$ . Figure 7 shows the non-axisymmetric solution at these parameters and clearly illustrates the dominance of the $m=1$ but with additional complexity due to the background flow. The dominance of the $m=1$ mode was confirmed by a spherical harmonic analysis (§ 4.2). The experimental parameters of the 3 m experiment are clearly located beyond the extrapolated $m=1$ onset, which may explain the prominence of this instability in their data (Triana Reference Triana2011). The switchover between $m_{c}=1$ and $m_{c}=2$ takes place between $E=3\times 10^{-5}$ and $E=10^{-5}$ . At the parameters of the experiments at BTU C-S, $m_{c}=1$ may indeed be the primary instability.

Wicht (Reference Wicht2014) and Hoff et al. (Reference Hoff, Harlander and Triana2016b ) identified the $(5,1,-0.068)$ mode as being the most probable candidate for this $m=1$ instability, based on a qualitative visual comparison of the spatial structure. However, the observed frequency has an opposite sign (opposite drift direction), which can probably be attributed to the fact that this $m=1$ mode is mostly confined to the region inside the TC where the height decreases towards the rotation axis, leading to a negative ‘ $\unicode[STIX]{x1D6FD}$ effect’ (Schaeffer & Cardin Reference Schaeffer and Cardin2005; Wicht Reference Wicht2014).

We have calculated the projection coefficients $c$ (equation (5.2)) for a wide selection of equatorially symmetric (ES) inertial modes for a more quantitative classification. For this purpose we used a single snapshot at $E=10^{-4},\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-0.4$ where the $m=1$ SI mode is simply drifting in azimuth. Since the dynamics is clearly dominated by a simple drift of the different modes, the projection analysis will not depend on the time of the snapshot. Figure 8 summarises the results and illustrates that $(5,1,-0.068)$ mode has the largest coefficient $c^{\prime }$ between 0.5 and 0.6, considering only the non-axisymmetric flow contribution. The next nearest contenders are the $(3,1,-0.177)$ and $(7,1,-0.036)$ with not much smaller coefficients between 0.3 and 0.4. In addition to the mismatch in drift direction, this shows that the identification in terms of a single inertial mode remains ambiguous. Therefore, we will refer to it as the $m=1$ SI (Stewartson layer instability).

Figure 8.

Projection on equatorially symmetric (ES) inertial modes of a velocity snapshot at $E=10^{-4},~\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-0.4$ . Here $|\cdots |$ represents absolute value and $c^{\prime }$ represents the coefficient when the axisymmetric part of the solution has been removed. Different colours and symbols denote different wavenumbers $m$ . The dominant modes have been annotated with $(l,m,\unicode[STIX]{x1D714}/\unicode[STIX]{x1D6FA})$ . The different points for the same $l$ and same colour and symbol ( $m$ ) represent different mode frequencies and structures for the same $(l,m)$ combination.

5.2 Inertial modes

As the magnitude of differential rotation is increased, at a certain critical value, ‘fast’ equatorially antisymmetric (EA) inertial modes are excited. However, depending on the Ekman number, one obtains a different set of modes and a slightly different behaviour of the kinetic energy at onset.

Analysing the MagIC results shows that, at $E=10^{-5}$ , the equatorial symmetry is broken at $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-0.35$ , but the EA kinetic energy remains too low to identify any underlying mode. A larger differential rotation leads to oscillations in the kinetic energy, attributable to nonlinear interactions between different modes. At $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-0.55$ , the first clear onset of EA inertial modes takes place. The dominant mode is the $(5,2,0.467)$ mode seen as a clear ‘bright’ line in the spectrograms of the experiment as well as the simulations (figure 5 a,b). Figure 9 compares a velocity snapshot at $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-0.6$ with eigenmodes of a sphere and a spherical shell. While the full numerical solution looks similar to the eigenmodes, it has also been evidently modified by the background flow due to the shear on and inside the TC. Figure 10(a) compares the $m=2$ equatorially antisymmetric contribution of the MagIC simulations to that of the reconstructed velocity from PIV measurements in the experiments. The PIV data have been filtered to retain only a narrow frequency band and cover only part of the planar circle (see figure 2(b) in Hoff et al. (Reference Hoff, Harlander, Egbers, Triana, Egbers, Ruck, Leder and Dopheide2016a )). The filtering process of the numerical solution and the experiments guarantees that both screen out everything else but the $(5,2,0.467)$ mode. The full circle view of the experimental results shown in figure 10 had to be reconstructed from several measurements, leaving a gap in the centre. The experimental flow reconstruction errors are estimated to be up to 15 % in the centre of the field of view and up to 25 % near corners and edges (Hoff et al. Reference Hoff, Harlander, Egbers, Triana, Egbers, Ruck, Leder and Dopheide2016a ). These limitations along with the different scaling used in the numerical simulation complicates a direct comparison. However, both agree on a qualitative level. The dominance of the $(5,2,0.467)$ mode is also testified by its relatively large projection coefficient, as illustrated in figure 12(a).

Figure 9.

Comparison of flow structures between simulation (a), analytical prediction for a full sphere (b) based on Zhang et al. (Reference Zhang, Earnshaw, Liao and Busse2001) and spherical shell eigenmode computed using SINGE (c) for the $(5,2)$ mode at $E=10^{-5}$ , $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-0.6$ . Shown here are isosurfaces of non-axisymmetric cylindrically radial velocity ( $u_{s}$ ) of a snapshot of the solution. Red is positive (outwards) while blue is negative (inwards). Note that no filtering has been performed on the solution. The whole solution is dominated by this mode at these parameters. The similarities are self-evident.

Figure 10.

Comparison of $z$ -slices of velocity magnitudes from experiments and simulations for (a) the $(5,2,0.467)$ mode and (b) the $(3,2,0.667)$ mode. The experimental data have been taken at $z/r_{o}=1/3$ , while the simulation data have been shown at $z/r_{o}=0.35$ with the centre of the spheres as the origin. The second column shows a slice from a MagIC simulation without a central part and saturated colours so as to bring out the similarities better. The third column shows the full slice from the simulation. The structures look very similar for both modes. For the $(5,2,0.467)$ mode the data have been taken at $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-0.68$ for the experiment and at $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-0.7$ for the simulation. For the $(3,2)$ mode, the corresponding numbers are $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-1.11$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-1.1$ , respectively. Colours and arrow lengths represent velocity magnitudes – horizontal velocity for experiments and total velocity for simulations. The arrows also indicate velocity direction.

Figure 11.

Same as figure 9 for the $(3,2)$ mode at $E=10^{-5}$ , $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-1$ .

Figure 12.

Projection of velocity snapshots on equatorially antisymmetric (EA) inertial modes of a sphere at different parameters. The major modes have been annotated using $(l,m,\unicode[STIX]{x1D714}/\unicode[STIX]{x1D6FA})$ . Here $c^{\prime }$ represents the projection coefficient of the snapshot with the axisymmetric part removed. Error bars show the temporal standard deviation in $|c^{\prime }|$ . We used the following number of snapshots for computing the time average: (a) five, (b) five, (c) seven and (d) six.

Table 2 compares the frequencies of the $(5,2,0.467)$ mode found in different experimental set-ups as well as the two numerical simulations presented here and the eigenvalues for a full sphere and a spherical shell. The two eigenvalues agree perfectly. Numerical and experimental values are nearly in agreement with each other but show a slightly larger frequency than the eigenvalue predictions. A possible reason is the effect of the background flow, which is completely ignored in eigenmode computations. The dependence on differential rotation will be discussed in § 6. Since the values in table 2 span four decades in Ekman number, the dependence of the frequency on Ekman number can only be weak.

Besides the two most dominant lines in the spectrogram for the $m=1$ SI and the $(5,2,0.467)$ modes, there are several other lines which are less dominant but nevertheless clearly visible against the background. Prominent among these are two ‘side branches’ around $\unicode[STIX]{x1D714}/\unicode[STIX]{x1D6FA}=0.5$ and $\unicode[STIX]{x1D714}/\unicode[STIX]{x1D6FA}=0.44$ , identified as the $(4,3,0.5)$ and an $m=1$ EA mode, respectively. The latter is difficult to identify unambiguously as a single inertial mode due to its complex structure. However, we identify it as the $(6,1,0.44)$ since the frequency is close (as seen in table 3) and the respective projection coefficient is the largest of all $m=1$ EA modes (figure 12 a). Since the time behaviour is more complex, the calculation of the projection coefficients now relies on several snapshots. Figure 12 shows mean coefficients and illustrates the standard deviation with error bars. Note that table 3 also shows a strong influence of the background flow on the frequency, as will be discussed in § 6. Furthermore, MagIC results suggest an earlier onset at around $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-0.55$ , while in the BTU C-S experiments, the mode was not detected before $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-1.54$ . A possible explanation for this discrepancy is the fact that the mode is difficult to identify in the experiment because of its complex structure and low amplitude.

However, at $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-0.8$ these modes disappear, the EA kinetic energy decreases significantly and the $m=1$ SI mode again becomes the dominant large-scale structure. This regime is followed by an onset of a different set of EA modes at $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-0.85$ consisting of the dominant $(3,2,0.667)$ mode with ‘side branches’ of $(4,1,0.612)$ and $(10,3,0.851)$ modes, as seen in the spectrograms in figures 5(a) and 5(b). The visual comparison of a numerical snapshot at $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-1$ with eigenmodes of a full sphere and a spherical shell in figure 11 clearly attests to the large $(3,2,0.667)$ contribution. Figure 10(b) shows the similarity of the respective mode in the simulation and experiment. The projection coefficients shown in figure 12(b) quantify the dominance of the $(3,2,0.667)$ mode. The projection coefficients substantiate that the $(3,2,0.667)$ contribution amounts to about 80 % of the non-axisymmetric solution. Once again the mode identification is challenging for the subdominant modes, but a combined analysis of projection coefficients and frequency matching helps us single out these modes as the best candidates. Table 2 lists the frequencies of the $(3,2,0.667)$ mode eigenmodes for a sphere and a spherical shell as well as for the solutions observed in simulations and experiments. The deviations from the eigenvalues are probably due to the effect of background flow, as will be discussed in § 6. These EA modes along with the $m=1$ SI remain the dominant modes till the largest differential rotation magnitude that could be afforded using existing computational resources. The limits were $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-1.4$ for MagIC and $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-1.17$ for XSHELLS.

However, at an Ekman number of $E=10^{-4}$ we could run simulations with larger differential rotation magnitudes and consequently venture into the turbulent regime. At this higher Ekman number, a similar onset of EA inertial modes is observed. However, unlike at $E=10^{-5}$ , these modes remain dominant right up to the onset of turbulence. After the onset of the $m=1$ SI mode at $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-0.379$ , we find kinetic energy oscillations due to nonlinear mode interactions starting at $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=1.25$ . These give way to EA inertial modes at $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-1.33$ , determined using lower steps in $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}$ . The $(3,2,0.667)$ mode is the most dominant mode accompanied by ‘side branches’ of $(4,1,0.612)$ and $(10,3,0.851)$ modes. They can be seen as three ‘bright’ lines in the spectrogram (figure 5 c) around $\unicode[STIX]{x1D714}/\unicode[STIX]{x1D6FA}=0.7,~0.6,~0.85$ , respectively. In the projections shown in figures 12(c) and 12(d), the $(3,2,0.667)$ , $(4,1,0.612)$ and $(10,3,0.851)$ modes are seen to be dominant. The modes mentioned above will hereafter be referred to using only the indices $(l,m)$ , as only a single mode for an $(l,m)$ combination is observed.

Tables 2 and 3 compare the frequencies and $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}$ ranges of the observed inertial modes in our simulations with those observed in past experiments. It is evident that the observations of these modes are very robust and similar modes are seen across several orders of magnitudes in Ekman number.

At $E=10^{-4}$ , the transition to turbulence takes place around $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-2.3$ , which is larger in magnitude as compared to the value of $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-1.73$ at $E=1.52\times 10^{-5}$ , as found by Hoff et al. (Reference Hoff, Harlander and Triana2016b ).

The frequencies of the dominant $(3,2)$ and the $m=1$ SI modes are more strongly influenced by the background flow, including a sudden jump in the frequency of the $(3,2)$ mode, similar to observations by Hoff et al. (Reference Hoff, Harlander and Triana2016b ). At around $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-2.7$ , a new equatorially symmetric mode with $m=2$ sets in but the corresponding inertial mode has not yet been determined. The turbulent regime, though very interesting, is beyond the scope of the present study and will not be discussed further.

7.1 Theoretical background

The phenomenon of triadic resonance may appear in any physical problem involving nonlinear wave interactions. Coupled harmonic oscillators, photon decay and waves in plasma (see e.g Bellan Reference Bellan2008) are just a few examples. In a geophysical context, they have been studied in relation to internal waves in the oceans (McComas & Bretherton Reference McComas and Bretherton1977) and resonant interactions of Rossby waves (Pedlosky Reference Pedlosky1987) and their instability (Zhang & Pedlosky Reference Zhang and Pedlosky2007). Triadic resonances are also observed when a flow becomes unstable to inertial modes. Examples of such instabilities are the elliptical and shear instabilities in precession (Kerswell Reference Kerswell1993).

It can be shown (see appendix A) that the nonlinear interaction of two inertial modes $\boldsymbol{Q}_{1}\text{e}^{\text{i}(m_{1}\unicode[STIX]{x1D719}-\unicode[STIX]{x1D714}_{1}t)}$ and $\boldsymbol{Q}_{2}\text{e}^{\text{i}(m_{2}\unicode[STIX]{x1D719}-\unicode[STIX]{x1D714}_{2}t)}$ can resonantly excite a third mode $\boldsymbol{Q}_{k}\text{e}^{\text{i}(m_{k}\unicode[STIX]{x1D719}-\unicode[STIX]{x1D714}_{k}t)}$ if the following conditions are met:

1. (i) The wavenumbers have to match

   (7.1) $$\begin{eqnarray}\displaystyle & \displaystyle m_{k}=m_{1}\pm m_{2}. & \displaystyle\end{eqnarray}$$

2. (ii) The frequencies have to fulfil

   (7.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}_{k}=\unicode[STIX]{x1D714}_{1}\pm \unicode[STIX]{x1D714}_{2}. & \displaystyle\end{eqnarray}$$

3. (iii) The integral over $\text{d}\unicode[STIX]{x1D70F}=r^{2}\sin \unicode[STIX]{x1D703}\,\text{d}r\,\text{d}\unicode[STIX]{x1D703}$ provides additional selection rules. An obvious one concerns the equatorial symmetry. The product of two vector fields of opposite equatorial symmetry vanishes when integrated over volume. Thus, the integral over $\text{d}\unicode[STIX]{x1D70F}$ is non-zero only when $\boldsymbol{Q}_{k}^{\dagger }$ and hence $\boldsymbol{Q}_{k}$ has the same equatorial symmetry as the nonlinear terms involving $\boldsymbol{Q}_{1}$ and $\boldsymbol{Q}_{2}$ . If $\boldsymbol{Q}_{1}$ and $\boldsymbol{Q}_{2}$ have opposite equatorial symmetries, then $\boldsymbol{Q}_{k}$ has to be equatorially antisymmetric. If both of them possess the same equatorial symmetry, $\boldsymbol{Q}_{k}$ has to be equatorially symmetric.

4. (iv) The $\text{d}\unicode[STIX]{x1D70F}$ integral may also define additional constraints for the radial and latitudinal structure that are not explored here.

7.2 Triadic resonances in the inertial mode regime

All the identified ‘fast’ inertial modes form pairwise triadic resonances with the fundamental $m=1$ SI mode. For example, at $E=10^{-4}$ , we identify the triads $[m=1\;\text{SI},(3,2),(4,1)]$ and $[m=1\;\text{SI},(10,3),(3,2)]$ . The triads obviously fulfil the wavenumber criterion and figure 17 confirms that their frequencies indeed add up to nearly zero. As for spherical harmonics, inertial modes are equatorially symmetric (antisymmetric) when the difference $l-m$ is even (odd). The equatorial symmetry criterion is thus also fulfilled.

Figure 18 shows the three-dimensional structure of the nonlinear advective acceleration of the triad $[(3,2),(4,1),m=1\;\text{SI}]$ in the left column. This is compared with the effective acceleration of the mode with the respective wavenumber, frequency and equatorial symmetry on the right. While the driving has a richer radial and latitudinal structure, only the simpler modes are resonantly driven.

Figure 17.

Triadic resonances at $E=10^{-4}$ . Here $\unicode[STIX]{x1D714}_{diff}$ refers to the lines of the differences: $\unicode[STIX]{x1D714}_{(3,2)}-\unicode[STIX]{x1D714}_{(4,1)}-\unicode[STIX]{x1D714}_{m=1\,\text{SI}}$ and $\unicode[STIX]{x1D714}_{(10,3)}-\unicode[STIX]{x1D714}_{(3,2)}-\unicode[STIX]{x1D714}_{m=1\,\text{SI}}$ ; $\unicode[STIX]{x1D714}_{diff}$ is close to zero at all times.

Figure 18.

Comparison of the three-dimensional structure of the advective acceleration due to the nonlinear interaction of two modes on the left with the structure of the effective acceleration of the third mode of the triad on the right, at $E=10^{-4}$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-2$ . Shown here are isosurfaces of cylindrically radial ( $s$ ) components. Red denotes outwards (positive), while blue denotes inwards (negative).

To check that the triads also preserve a consistent phase locking, we compute the bicoherence, defined as

(7.3) $$\begin{eqnarray}\displaystyle & \displaystyle B(\unicode[STIX]{x1D714}_{1},\unicode[STIX]{x1D714}_{2})=\frac{\displaystyle \left|\mathop{\sum }_{n}F_{n}(\unicode[STIX]{x1D714}_{1})F_{n}(\unicode[STIX]{x1D714}_{2})F_{n}^{\dagger }(\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D714}_{2})\right|^{2}}{|F_{n}(\unicode[STIX]{x1D714}_{1})F_{n}(\unicode[STIX]{x1D714}_{2})|^{2}|F_{n}(\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D714}_{2})|^{2}}, & \displaystyle\end{eqnarray}$$

where $F(\unicode[STIX]{x1D714})$ represents the Fourier transform of a signal at frequency $\unicode[STIX]{x1D714}$ and $^{\dagger }$  denotes a complex conjugate. This is a statistical measure of the quadratic phase coupling (Nikias & Raghuveer Reference Nikias and Raghuveer1987) determined by splitting a signal into several parts, computing the ‘bispectrum’ (one of the terms in the numerator) of individual parts and eventually summing them up. Wave interactions that preserve the same phase relation over the different parts of the time series result in a particularly strong $B(\unicode[STIX]{x1D714}_{1},\unicode[STIX]{x1D714}_{2})$ .

We compute the bicoherence for a time series of zonal velocity, $u_{\unicode[STIX]{x1D719}}$ , that we split into 16 segments overlapping by 50 % using the HOSA Matlab toolbox (Swami, Mendel & Nikias Reference Swami, Mendel and Nikias1998). For each segment, 4096 points are used for the fast Fourier transform. The bicoherence $B(\unicode[STIX]{x1D714}_{1},\unicode[STIX]{x1D714}_{2})$ attains a high value for all identified triadic resonances, as shown in figure 19. For example, take the peak of $(\unicode[STIX]{x1D714}_{1},\unicode[STIX]{x1D714}_{2})=(0.61,0.11)$ representing the coupling for the $(4,1)$ mode and $m=1$  SI mode seen in this regime. In this case $B(0.61,0.11)$ has a very high value, as these modes are phase coupled to the $(3,2)$ mode, having a frequency of $0.61+0.11=0.72$ . The couplings with higher frequencies denote triadic resonances between pairs of higher harmonics of the $(3,2)$ , $(4,1)$ and $(10,3)$ modes and the higher harmonics of the $m=1$  SI mode. Strong self-interactions along the main diagonal can also be seen, especially for the dominant $(3,2)$ mode.

Figure 19.

Bicoherence of $u_{\unicode[STIX]{x1D719}}$ time series taken at $r=0.7r_{o}$ , computed for $E=10^{-4}$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D6FA}=-2$ . All frequencies are non-dimensionalised by $\unicode[STIX]{x1D6FA}$ . Compare with figure 6 from Hoff et al. (Reference Hoff, Harlander and Triana2016b ).

The major modes forming triadic resonances at $E=10^{-5}$ and the degree to which they satisfy the frequency criterion are shown in figure 20. The pairs forming triads with the $m=1$  SI mode are $[(5,2),(6,1)]$ , $[(4,3),(5,2)]$ , $[(3,2),(6,1)]$ and $[(10,3),(3,2)]$ . These triads have also been found in experiments at BTU C-S (Hoff et al. Reference Hoff, Harlander, Egbers, Triana, Egbers, Ruck, Leder and Dopheide2016a ,Reference Hoff, Harlander and Triana b ).

Figure 20.

Triadic resonances at $E=10^{-5}$ . Here $\unicode[STIX]{x1D714}_{diff}$ refers to the lines of various frequency differences which form triadic resonances; $\unicode[STIX]{x1D714}_{diff}$ is close to zero at all times.