3.1. The uniform vorticity approach

The study of Zhang et al. (Reference Zhang, Chan and Liao2012) suggests that the laminar solution in a spheroidal container remains essentially of uniform vorticity, provided that $E\ll 1$ . The same was already observed in precessing ellipsoids (Poincaré Reference Poincaré1910; Stewartson & Roberts Reference Stewartson and Roberts1963; Busse Reference Busse1968; Noir & Cébron Reference Noir and Cébron2013). Along the same line of thought, we propose to seek for solutions of uniform vorticity ${\bf\omega}$ . In the following, however, it will be more convenient to use the rotation rate $\boldsymbol{q}={\bf\omega}/2$ as primary variable. As such, we may write

with ${\it\zeta}$ a gauge function that allows the velocity field $\boldsymbol{u}$ to be solenoidal and satisfy the impermeability condition (2.12). One can show that such a flow is an exact solution of the equations of motion (2.1) and (2.2) in the bulk of the cavity. This solution satisfies the non-penetration condition (2.12), but not necessarily the no-slip boundary condition (2.11). As a consequence we expect a boundary layer to develop in the vicinity of the walls, which in turn drives a secondary flow of order $E^{1/2}$ in the interior, the so-called Ekman pumping. In the limit of small Ekman number considered in this study, a classical approach would be to neglect the boundary layer and its associated Ekman pumping. Under this assumption, the vorticity equation (2.13) takes the inviscid form

Hence, the final state of the system always depends on the initial conditions. This can be avoided by reintroducing the viscosity, which requires an exact description of the Ekman layer in a triaxial ellipsoid. Such a calculation is a formidable task, and leads to cumbersome calculations. We thus choose to adopt the simpler but successful approach pioneered by Noir & Cébron (Reference Noir and Cébron2013) in the case of precessing triaxial ellipsoids. It consists of parametrizing the effect of viscous Ekman torques on the bulk flow by adding a dissipative term in the vorticity equation (3.2). These torques act so as to inhibit differential rotation between the motion of the fluid and the container. This term is proportional to $E^{1/2}$ and to the differential rotation between the container and the fluid, i.e.  $\boldsymbol{q}-\boldsymbol{q}_{M}$ , where $\boldsymbol{q}_{M}$ is the container rotation rate in the considered frame of reference. In the frame attached to the container, $\boldsymbol{q}_{M}=\mathbf{0}$ , and (3.2) simply becomes

where ${\it\lambda}>0$ parametrizes the rate of dissipation. In § 3.3, we will determine ${\it\lambda}$ by fitting the outcome of the uniform vorticity model to three-dimensional (3D) numerical simulations. As noted by Noir & Cébron (Reference Noir and Cébron2013), the reduced model does not take into account secondary viscous effects (e.g. internal shear layers).

Following Hough (Reference Hough1895), Roberts & Wu (Reference Roberts and Wu2011) and Noir & Cébron (Reference Noir and Cébron2013), the uniform vorticity velocity field satisfying the non-penetration condition (2.12) is given by

Substituting (3.4) into (3.3) yields

In contrast to (3.3), this formulation of the equations of motion takes the form of a set of ordinary differential equations (ODEs), and is computationally thus far less demanding.

3.2. Analysis in the limit of weak forcing

In this section, we seek an analytical solution of (3.5)–(3.7) in the limit of small libration amplitude, i.e.  ${\it\varepsilon}\ll 1$ . We assume an asymptotic development of the form

We now substitute this ansatz into (3.5)–(3.7), and successively analyse for increasing orders of ${\it\varepsilon}$ . Starting at order 0, and invoking that $\cos {\it\theta}=1-O({\it\varepsilon}^{2})$ , $\sin {\it\theta}=O({\it\varepsilon})$ and $\dot{{\it\theta}}={\it\varepsilon}\cos ({\it\omega}_{L}t)=O({\it\varepsilon})$ , $\ddot{{\it\theta}}=-{\it\omega}_{L}{\it\varepsilon}\sin ({\it\omega}_{L}t)=O({\it\varepsilon})$ , we find

To analyse this system, we introduce the Lyapunov function $F$ (Manneville Reference Manneville2010),

Using (3.9)–(3.11) and after some algebra, we obtain

which is strictly negative definite. From this, it follows immediately that the solution of (3.9)–(3.11) tends to $\boldsymbol{q}^{(0)}=\mathbf{0}$ for any initial condition.

We now proceed with our analysis to first order in ${\it\varepsilon}$ . Using the solution $\boldsymbol{q}^{(0)}=\mathbf{0}$ , we find

As a first step, we consider the inviscid case $E=0$ , for which (3.14)–(3.16) has the particular integral with The most remarkable feature of this solution is the singularity that occurs for ${\it\omega}_{L}=f$ , where $f$ is the eigenfrequency of the so-called spin-over mode (Vantieghem Reference Vantieghem2014). This suggests that there is a close connection between the present problem and the theory of inertial modes. In order to substantiate this, we first recast (3.14)–(3.16) with $E=0$ in matrix–vector form:

This notation emphasizes that the dynamical system that we study is simply a harmonic oscillator driven by the Poincaré force. Resonance, i.e. unbounded growth of the flow, takes place if the driving frequency ${\it\omega}_{L}$ matches one of the natural frequencies ${\it\mu}$ of the unforced system. These are solutions of the eigenvalue problem

where $\text{i}$ is the imaginary unit. This yields ${\it\mu}_{1,2}=\pm f$ and ${\it\mu}_{3}=0$ . Following Vantieghem (Reference Vantieghem2014) and using (3.4), (3.23) is the projection of the inertial mode equation in vorticity formulation,

on the subspace of solenoidal uniform vorticity flows that satisfy the impermeability condition (2.12). The eigenmodes associated with ${\it\mu}_{1,2}=\pm f$ are associated with the so-called spin-over mode, i.e. an inertial mode whose vorticity is uniform and located in the equatorial plane. Resonant coupling is not possible for the eigenvalue ${\it\mu}_{3}=0$ , since its eigenvector $(0,0,1)$ is orthogonal to the right-hand side of (3.22). Summarizing, the solution (3.17)–(3.19) embodies that latitudinal libration drives a spin-over mode, and that resonance takes place if the libration frequency matches the spin-over frequency. The linearized uniform vorticity theory thus confirms the inviscid theory of Chan et al. (Reference Chan, Liao and Zhang2011) and Zhang et al. (Reference Zhang, Chan and Liao2012), and extends it to triaxial ellipsoids.

For the biaxial ellipsoid with $b=c$ , i.e.  ${\it\beta}_{bc}=0$ , the solution (3.17)–(3.19) reduces after some algebra to

for any value of ${\it\omega}_{L}$ and does not exhibit resonant behaviour. In the frame of mean rotation (denoted by subscript $R$ ), this yields

This shows that the motion of the container cannot be transmitted to the fluid in the absence of container topography in planes perpendicular to the libration axis.

A further interesting choice is ${\it\omega}_{L}=f_{0}$ , which cancels $q_{x}^{(1)}$ . This implies that the streamlines of this flow are ellipses located in planes $y=\mathit{cst}$ . For $a=c$ , i.e.  ${\it\beta}_{ac}=0$ , this situation occurs if ${\it\omega}_{L}=1$ , and is characterized by circular streamlines.

To fix the non-physical behaviour at resonance, we return to the viscous case $(E\neq 0)$ . It is still possible to work out a particular solution for (3.14) and (3.15), although it is more intricate:

We can now distinguish between two cases, depending on whether the first term dominates the second one in the denominator of (3.29) and (3.30) or vice versa. Assuming that ${\it\omega}_{L},f,{\it\lambda}=O(1)$ , the former occurs if $|{\it\omega}_{L}-f|\gg \sqrt{E}$ , i.e. for libration frequencies far away from resonance. The first term in each of (3.29) and (3.30) is then of order of magnitude 1, whereas the second ones are $O(\sqrt{E})$ . For $E\rightarrow 0$ , we recover the inviscid solution (3.17) and (3.18).

For $|{\it\omega}_{L}-f|\ll \sqrt{E}$ , the first term in (3.29) and (3.30) is of order of magnitude 1, as the leading-order contributions in their numerator and denominator both scale as $O(E)$ . The second term in (3.29) and (3.30) now scales as $E^{-1/2}$ , and thus dominates the first term. We find thus that for $E\ll 1$ , the total solution scales as $E^{-1/2}$ and is established within a frequency window of width $O(E^{1/2})$ . These scalings are consistent with Zhang et al. (Reference Zhang, Chan and Liao2012). Readers whose main interest is in the dynamical stability of the solutions derived above, can now skip the remainder of this section without great loss.

Anticipating the numerical validation of our analysis, we introduce $e_{k}$ and $\overline{e_{k}}$ , the mean kinetic energy density and its time average, which can be readily computed from (3.4), and read

where $V=4abc{\rm\pi}/3$ is the volume of the ellipsoid, and $T=2{\rm\pi}/{\it\omega}_{L}$ the libration period. Upon substitution of (3.29) and (3.30) into this expression, we obtain the following result for resonant flows:

The second term in $O(E)$ may dominate if ${\it\beta}_{bc}\ll \sqrt{E}$ . In particular, if ${\it\beta}_{bc}=0$ , the expression for $\overline{e_{k}}$ for resonant and non-resonant flows is

For an observer in the frame of mean rotation, $\overline{e_{k}}$ for ${\it\beta}_{bc}=0$ is then given by

This indicates that, in the absence of topography (i.e. for ${\it\beta}_{bc}=0$ ), viscosity can transmit the librating motion of the mantle to the fluid. However, the flow becomes vanishingly weak for non-resonant frequencies as $E\rightarrow 0$ . At resonance, the amplitude of the flow is independent of $E$ , and thus $O(E^{-1/2})$ times stronger than the non-resonant flow. Motivated by the experimental study of Aldridge & Toomre (Reference Aldridge and Toomre1969), Zhang et al. (Reference Zhang, Chan, Liao and Aurnou2013) have argued that similar scalings hold for the case of longitudinal libration of a spherical container.

We recall that the above analytical solutions have been obtained under the a priori assumption that the non-linear terms in the equation of motion remain negligible with respect to the linear ones. A posteriori, we find that this is indeed established off-resonance, given that the flow scales as ${\it\varepsilon}$ . The amplitude of the resonant flow, however, scales as ${\it\varepsilon}{\it\beta}_{bc}E^{-1/2}$ , and the non-linear terms will therefore be negligible if ${\it\varepsilon}{\it\beta}_{bc}E^{-1/2}\ll 1$ .

In figure 4 we compare the solution (3.29)–(3.31) of the linearized ODEs to numerical solutions of the non-linear system (3.5)–(3.7) for $(a,b,c)=(1,\sqrt{5/4},\sqrt{3/4})$ . The libration frequency ${\it\omega}=\sqrt{10/7}$ is chosen such that resonance occurs; the value of ${\it\lambda}$ is set to 2.7 and the Ekman number to $5\times 10^{-5}$ . This choice of parameters is representative of the numerical results that will be discussed in § 3.3. The non-linear solutions for $q_{x}$ , $q_{y}$ and $e_{k}$ virtually collapse with their linear counterparts at ${\it\varepsilon}=0.005$ . For ${\it\varepsilon}=0.05$ , the amplitudes of $q_{x}$ and $q_{y}$ are approximately 10–20 % lower than predicted by the linear theory. We see furthermore that a non-zero retrograde solution for the vorticity component $q_{z}$ emerges when ${\it\varepsilon}$ increases; this feature is generated by non-linear interactions between the components of the linear solutions $\boldsymbol{q}^{(1)}$ . Indeed, in this case the quantity ${\it\varepsilon}{\it\beta}_{bc}E^{-1/2}\approx 1.77$ and thus non-linearities are no longer negligible, as argued in the previous paragraph. The $q_{z}$ -component consists of the superposition of a non-zero mean, referred to as a zonal flow (Busse Reference Busse2010; Sauret et al. Reference Sauret, Cébron, Morize and Le2010) and a harmonic modulation of frequency $2{\it\omega}_{L}$ . However, a complete investigation of the zonal flow requires us to take into account higher-order effects, such as non-linear interactions in the viscous Ekman layer. Such a detailed study would involve an exact description of the Ekman boundary layer and will not be pursued here.

Figure 4. Comparison between the exact solution (3.29) and (3.30) of the linearized ODEs (solid line) and direct integration of the non-linear ODEs (3.5)–(3.7) for ${\it\varepsilon}=0.005$ (circles) and 0.05 (squares) for $(a,b,c)=(1,\sqrt{5/4},\sqrt{3/4})$ and ${\it\omega}_{L}=f=\sqrt{10/7}$ . The distance between the dashed vertical lines spans one libration cycle $T=2{\rm\pi}/{\it\omega}_{L}$ .

3.3. Numerical validation and discussion

We now wish to validate our model of uniform vorticity flow. To this end, we will compare the results presented above against 3D non-linear numerical simulations obtained by means of a non-structured finite-volume code (Vantieghem Reference Vantieghem2011). It uses a collocated arrangement of the variables, and a second-order centred-finite-difference-like discretization stencil for the spatial differential operators. The time advancement algorithm is based on a canonical fractional-step method (Kim & Moin Reference Kim and Moin1985). More specifically, the procedure to compute the velocity and reduced pressure $\boldsymbol{u}^{j+1},p^{\,j+1}$ at time step $t^{j+1}=t^{j}+{\rm\Delta}t$ , given the respective variables at time step $t^{j}$ , is as follows.

1. (i) We first solve the intermediate velocity $\boldsymbol{u}^{\star }$ from the equation

   (3.40) $$\begin{eqnarray}\displaystyle \frac{\boldsymbol{u}^{\star }-\boldsymbol{u}^{j}}{{\rm\Delta}t} & = & \displaystyle -\boldsymbol{u}_{AB}^{j+1/2}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{CN}^{j+1/2}-\boldsymbol{{\rm\nabla}}p^{j}\nonumber\\ \displaystyle & & \displaystyle -\,2{\it\bf\Omega}^{j+1/2}\times \boldsymbol{u}^{j+1/2}+E{\rm\nabla}^{2}\boldsymbol{u}_{CN}^{j+1/2}+\dot{{\it\bf\Omega}}^{j+1/2}\times \boldsymbol{r},\end{eqnarray}$$
   with no-slip boundary condition $\boldsymbol{u}_{bnd}^{\star }=\mathbf{0}$ . In this expression, $\boldsymbol{u}_{AB}^{j+1/2}$ and $\boldsymbol{u}_{CN}^{j+1/2}$ denote the velocity at time step $j+(1/2)$ computed using a second-order Adams–Bashforth, respectively Crank–Nicolson, approach, i.e. 
   (3.41) $$\begin{eqnarray}\displaystyle & \boldsymbol{u}_{AB}^{j+1/2}={\textstyle \frac{3}{2}}\boldsymbol{u}^{j}-{\textstyle \frac{1}{2}}\boldsymbol{u}^{j-1}, & \displaystyle\end{eqnarray}$$

   (3.42) $$\begin{eqnarray}\displaystyle & \boldsymbol{u}_{CN}^{j+1/2}={\textstyle \frac{1}{2}}(\boldsymbol{u}^{j}+\boldsymbol{u}^{\star }). & \displaystyle\end{eqnarray}$$
   The mixed Adams–Bashforth/Crank–Nicolson formulation for the advective term has the advantage of being kinetic-energy conserving and time-stable for any ${\rm\Delta}t$ , and it does not require the solution of a non-linear system for the unknown $\boldsymbol{u}^{\star }$ .

2. (ii) The new velocity $\boldsymbol{u}^{j+1}$ is then related to $\boldsymbol{u}^{\star }$ by

   (3.43) $$\begin{eqnarray}\boldsymbol{u}^{j+1}=\boldsymbol{u}^{\star }-{\rm\Delta}t\left({\rm\Delta}p^{j+1}\right)\!,\end{eqnarray}$$
   with ${\rm\Delta}p^{j+1}=p^{j+1}-p^{j}$ . Imposing the incompressibility constraint on $\boldsymbol{u}^{j+1}$ leads to a Poisson equation for ${\rm\Delta}p^{j+1}$ ,
   (3.44) $$\begin{eqnarray}{\rm\nabla}^{2}\left({\rm\Delta}p^{j+1}\right)=({\rm\Delta}t)^{-1}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}^{\star },\end{eqnarray}$$
   with boundary condition $\hat{\boldsymbol{n}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}({\rm\Delta}p^{j+1})=0$ . This Poisson equation is solved with the algebraic multigrid method BoomerAMG (Henson & Yang Reference Henson and Yang2002). Discussions regarding the performance of this code can be found in Cébron, Vantieghem & Herreman (Reference Cébron, Vantieghem and Herreman2014) and Marti et al. (Reference Marti, Schaeffer, Hollerbach, Cébron, Nore, Luddens, Guermond, Aubert, Takehiro, Sasaki, Hayashi, Simitev, Busse, Vantieghem and Jackson2014).

The computational grid consists of a mixture of tetrahedral, prismatic and hexahedral elements. The total number of control volumes (CVs) at the lowest Ekman number, $E=5\times 10^{-5}$ , is approximately 2.3 million. This corresponds to a typical CV size ${\rm\Delta}x$ of about 0.012. In the bulk of the mesh, the CVs are quasi-isotropic. Near the boundaries, however, they become more anisotropic as we reduce the wall-normal spacing to ${\rm\Delta}x\approx 0.0012\lesssim 0.2E^{1/2}$ in order to resolve the viscous Ekman layers.

In the following, we will discuss numerical results for four different geometries, which we label with Roman numerals I–IV, and whose properties are summarized in table 1. Case I is exactly the same geometry as studied by Zhang et al. (Reference Zhang, Chan and Liao2012) so as to allow direct comparison between both theories. In case II, we have a biaxial ellipsoid with $b=c$ for which no resonance can occur. Case III is also a biaxial ellipsoid, but now with $a=c$ . Finally, case IV is a triaxial ellipsoid with three different semi-major axes, and will allow us to demonstrate the generality of the present theoretical and numerical approach.

In figure 5, we compare the total kinetic energy density (3.32) retrieved from our numerical simulations with the linear result (3.29)–(3.31), (3.34). We emphasize that the numerically obtained energy density also contains non-uniform vorticity contributions, as the integral over the ellipsoidal volume also takes into account the viscous boundary layers, for example. Nevertheless, there is excellent agreement between the two approaches. The coefficient ${\it\lambda}_{num}$ used to trace the theoretical curve is chosen to minimize the discrepancy with the numerical data points, and is given in the rightmost column of table 1. This value can be compared with the theoretical asymptotic decay rate ${\it\lambda}_{asym}$ of the spin-over mode, given in the penultimate column of table 1. We see that there is a striking quantitative resemblance for the three geometries with topographic coupling (i.e. cases I, III and IV). The slight systematic higher value of the dissipation rate is due to the simplified expression of the dissipation used in our reduced model, and the moderate non-asymptotic values of the Ekman number. Therefore, at first order, one can interpret the dissipation rate in our model as the decay rate of the spin-over mode. We note finally that we have not inverted ${\it\lambda}$ for the non-resonant case II. According to (3.36) and (3.37), the energy density is independent of ${\it\lambda}$ to leading order in $E$ for this geometry. Far away from resonance, we see that $\overline{e_{k}}\approx 0.1$ ; near the spin-over frequency, however, $\overline{e_{k}}$ tends to 0.05. This agrees well with (3.36) and (3.37).

Figure 5. Mean kinetic energy density $\overline{e_{k}}$ for $E=5\times 10^{-5}$ and ${\it\varepsilon}=0.005$ at resonance. Comparison between the result of the linear theory (3.29)–(3.31), (3.33) and (3.34) and the numerical simulations (symbols). Results for geometries I–IV. Note the difference in scale between the different subfigures.

In figure 6, we compare the time series of $e_{k}$ retrieved from the numerical simulations with those from the linear theory. In order to trace the theoretical curve, the coefficient ${\it\lambda}$ was set to the value extracted from the frequency scan, i.e. the one provided in the rightmost column of table 1. Both curves virtually collapse, and this further validates our initial ansatz of uniform vorticity flow and the reduced viscosity model.

Figure 6. Time series of the mean kinetic energy density $e_{k}$ for $E=5\times 10^{-5}$ , ${\it\varepsilon}=0.005$ , and geometry IV. Results for ${\it\omega}_{L}=\sqrt{10/7}$ (resonant, a) and ${\it\omega}_{L}=1.5$ (non-resonant, b). Comparison between the result of the linear theory (3.33) (solid line) and the numerical simulations (dashed line).

We now investigate how well the analytical expression (3.33) is established at different values of the Ekman number; our results are summarized in figure 7. For the spheroidal geometry I, we can also compare our results to the asymptotic expression of Zhang et al. (Reference Zhang, Chan and Liao2012). Generally, we find that the agreement between the theory and the numerics is excellent. Noticeable differences are only observed for $E\gtrsim 4\times 10^{-4}$ for the triaxial case IV.

Figure 7. Mean kinetic energy density $\overline{e_{k}}$ for $E=5\times 10^{-5}$ and ${\it\varepsilon}=0.005$ . Comparison between the result of the linear theory (3.29)–(3.31), (3.33) and (3.34) (solid line) and the numerical simulations (symbols) for the cases with topographic coupling I (diamonds), III (squares) and IV (circles). Stars correspond to the asymptotic theory of Zhang et al. (Reference Zhang, Chan and Liao2012) for a spheroidal geometry.

The mean kinetic energy density is a global measure of the amplitude of the driven flow, but we also wish to assess how well the uniform vorticity solution is established locally. To this end, we show snapshots of $u_{y}$ in the plane $x=0$ and $u_{z}$ in the plane $z=0$ in figure 8. According to the uniform vorticity assumption (3.4), we expect these isolines to be straight curves, with $u_{y}$ being independent of $y$ in the plane $x=0$ . For non-resonant frequencies, we see that the uniform vorticity solution is reasonably well established in the bulk of the flow, whereas there are noticeable differences for the resonant flow. Finally, in both cases, the theoretical solution breaks down near the container walls, i.e. in the viscous Ekman layers. Similar behaviour was observed by Zhang et al. (Reference Zhang, Chan and Liao2012). The resonant case is further investigated in figure 9, where we have removed the uniform vorticity component from the flow and show a snapshot of $u_{y}$ in the plane $x=0$ . We observe a conical-like structure that is reminiscent of internal shear layers (Kerswell Reference Kerswell1995). As the thickness and amplitude of these layers scale as $E^{1/5}$ and $E^{3/10}$ respectively, this pattern is diffuse and smeared out over a broad area.

Figure 8. Isocontours of $u_{y}$ in the plane $x=0$ (a,c) and $u_{z}$ in the plane $z=0$ (b,d) for geometry IV and libration frequencies ${\it\omega}=1$ (a,b) and ${\it\omega}=f=\sqrt{10/7}$ (c,d).

Figure 9. Instantaneous isocontours of the numerically obtained magnitude of the non-uniform vorticity contribution to the flow in the meridional plane $x=0$ for geometry IV, $E=5\times 10^{-5}$ and the resonant driving frequency ${\it\omega}_{L}=f=\sqrt{10/7}$ .

Finally, we also assess the stability of the numerically obtained solutions. To this end, we adopt the following approach. For a given solution, we consider a perturbation that consists of a random velocity field with zero mean and a root-mean-square amplitude exceeding that of the original solution by a factor of approximately two. Starting from this ‘initial condition’, we integrate the Navier–Stokes equation in time. Systematically, we observe that the system quickly returns to its original state, which gives evidence for the stability of the numerically obtained solutions. This is illustrated in figure 10, where we show time series of the kinetic energy density $e_{k}$ before and after the perturbation for the triaxial geometry IV, ${\it\omega}_{L}=f$ and $E=5\times 10^{-5}$ .

Figure 10. Kinetic energy density $e_{k}{\it\varepsilon}^{-2}$ before and after perturbation with a random velocity field (at $t=500$ ) for the triaxial geometry IV, $E=5\times 10^{-5}$ and the resonant driving frequency ${\it\omega}_{L}=f=\sqrt{10/7}$ .

4.1. Energy considerations

Prior to seeking direct solutions for (4.7), we will first investigate it through the prism of the power statement

which is obtained by taking the inner product of (4.7), and introduces the strain rate tensor $\unicode[STIX]{x1D643}$ ,

Here we have introduced the notation $\langle \,\boldsymbol{f},\boldsymbol{g}\rangle =\int _{V}\boldsymbol{f}\boldsymbol{\cdot }\boldsymbol{g}^{\dagger }\text{d}V$ . This teaches us that amplification of $\boldsymbol{v}$ , and thus the presence of instability, requires that the right-hand side of (4.8) does not vanish. Using (3.4), (3.17) and (3.18), $\unicode[STIX]{x1D643}$ reads

with

By virtue of the Cauchy–Schwarz inequality, we infer that

This allows us to bound the growth rate ${\it\sigma}$ of a perturbation $\boldsymbol{v}\propto \exp ({\it\sigma}t)$ , since

where $|h_{1}({\it\omega}_{L})|,|h_{2}({\it\omega}_{L})|\leqslant 1$ , and where the bar notation denotes an adequate time average that allows us to account for an additional harmonic time-dependence of $\boldsymbol{v}$ . Given that $h_{1}$ and $h_{2}$ only depend on ${\it\omega}_{L}$ , any dependence of the instability on the geometry and ${\it\varepsilon}$ is captured by the parameters ${\it\chi}_{1}$ and ${\it\chi}_{2}$ . We immediately observe stability for ${\it\beta}_{bc}=0$ . Note that for ${\it\omega}_{L}=f_{0}^{2}=1+{\it\beta}_{ac}$ , ${\it\chi}_{1}={\it\chi}_{2}$ , and for ${\it\omega}_{L}=1$ and arbitrary ${\it\beta}_{ac}$ , ${\it\chi}_{1}=-{\it\chi}_{2}$ . In these cases, there is thus only a single control parameter. In particular, ${\it\chi}_{1}={\it\chi}_{2}=0$ for ${\it\omega}_{L}=1$ and ${\it\beta}_{ac}=0$ . No instability can grow for such a configuration regardless of the value of the topographic coupling ${\it\beta}_{bc}$ . Intuitively, this can be understood from inspection of the flow, (3.4), (3.17) and (3.18), which has only unstrained, circular streamlines in planes $y=\mathit{cst}$ .

4.2. Local method: short-wavelength Lagrangian stability analysis

The approach we follow here is based on the short-wavelength Lagrangian theory, used by Bayly (Reference Bayly1986) and Craik & Criminale (Reference Craik and Criminale1986), and then generalized in Friedlander & Vishik (Reference Friedlander and Vishik1991), Lifschitz & Hameiri (Reference Lifschitz and Hameiri1991, Reference Lifschitz and Hameiri1993) and Lifschitz (Reference Lifschitz1994), where the whole theory is thoroughly explained. This theory is now rather classical in stability studies of flows (e.g. Bayly, Holm & Lifschitz Reference Bayly, Holm and Lifschitz1996; Lebovitz & Lifschitz Reference Lebovitz and Lifschitz1996; Leblanc & Cambon Reference Leblanc and Cambon1997), and we thus only recall below some basic elements of the stability analysis, following the approach of Le Dizès & Eloy (Reference Le Dizès and Eloy1999). We found it simplest to work in the inertial frame of reference. The perturbation velocity $\boldsymbol{v}$ is written in the geometrical optics, or WKB (Wentzel–Kramers–Brillouin) form:

Here, the amplitude $\boldsymbol{a}(\boldsymbol{r},t)$ and phase ${\it\psi}(\boldsymbol{r},t)$ are real functions dependent on space $\boldsymbol{r}$ and time $t$ . The characteristic wavelength ${\it\vartheta}\ll 1$ is the small parameter used for the asymptotic (WKB) expansion. In the inviscid limit, the evolution of (4.15) is governed by the linearized Euler equations. Along the pathlines of the flow $\boldsymbol{U}$ in the inertial frame of reference, the leading-order problem can then be written in Lagrangian form as a system of ODEs (Lifschitz Reference Lifschitz1994):

with constraint

Here $\text{d}/\text{d}t=\partial _{t}+\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}$ are Lagrangian derivatives, $\unicode[STIX]{x1D644}$ is the identity matrix, $\boldsymbol{K}=\boldsymbol{{\rm\nabla}}{\it\chi}$ is the (local) wavevector along the Lagrangian trajectory $\boldsymbol{X}$ . The incompressibility condition (4.19) is always fulfilled if the initial condition $(\boldsymbol{X}_{0},\boldsymbol{K}_{\mathbf{0}},\boldsymbol{a}_{0})$ satisfies $\boldsymbol{K}_{\mathbf{0}}\boldsymbol{\cdot }\boldsymbol{a}_{0}=0$ (Le Dizès Reference Le Dizès2000). As shown by Lifschitz & Hameiri (Reference Lifschitz and Hameiri1991), the existence of an unbounded solution for $\boldsymbol{a}$ provides a sufficient condition of instability. Assuming closed pathlines, stability is naturally analysed over one turnover period $T$ along the pathline. Note that this system of equations can be seen as an extension of rapid distortion theory (RDT) to non-homogeneous flows (Cambon, Teissedre & Jeandel Reference Cambon, Teissedre and Jeandel1985; Cambon et al. Reference Cambon, Benoit, Shao and Jacquin1994; Sipp & Jacquin Reference Sipp and Jacquin1998).

In practice, (4.16) has to be solved as a first step to know the trajectory $\boldsymbol{X}$ emerging out of initial position $\boldsymbol{X}_{0}$ . Knowing $\boldsymbol{X}$ , one can solve the wavevector equation (4.17) for an initial vector $\boldsymbol{K}_{0}$ . As the magnitude of $\boldsymbol{K}_{0}$ does not influence the growth of $\boldsymbol{a}$ , we only have to consider the spherical surface $\Vert \boldsymbol{K}_{\mathbf{0}}\Vert =1$ in wave space to find the maximum growth rate. Knowledge of $\boldsymbol{X}$ and $\boldsymbol{K}$ finally allows us to solve (4.18) for the amplitudes $\boldsymbol{a}$ and to look for growing solutions.

In figure 11, we present typical results of the local stability analysis of the basic flow (3.17) and (3.18) for different geometries. We vary the libration frequency ${\it\omega}_{L}$ between 0 and 4, but exclude a small frequency band around the spin-over frequency $f$ where the inviscid base flow diverges. We choose to present our results as a function of $({\it\omega}_{L}-f)/(4-f)$ so that this frequency band collapses for all geometries. More precisely, we display ${\it\sigma}/{\it\chi}$ , where ${\it\chi}=|{\it\chi}_{1}|+|{\it\chi}_{2}|$ . Both figures 11(a) and 11(b) clearly illustrate that ${\it\sigma}/{\it\chi}<1$ , which is consistent with the upper bound (4.14). We consider the simple case of a rotating spheroid $a=b$ in figure 11(a). Since two control parameters ( ${\it\chi}_{1}$ , ${\it\chi}_{2}$ ) are available in this case, it is not possible to collapse all the curves on a single one. For the case $a=c$ , however, there is only a single control parameter ${\it\chi}_{2}$ , and thus ${\it\chi}=|{\it\chi}_{2}|$ . Figure 11(b) clearly shows that all curves then collapse onto a single master curve, which actually depends only on the function $h_{2}({\it\omega}_{L})$ introduced in § 4.1. One can notice that the local stability results indicate that, for $a=c$ , the function $h_{2}({\it\omega}_{L})$ is bounded by $1/2$ when we are not on the resonance peak.

Figure 11. Results of the inviscid local stability analysis, carried out at fixed ${\it\varepsilon}=0.04$ . The thick vertical black line indicates the resonance area. (a) ${\it\sigma}/{\it\chi}$ for a spheroid ( $a=b=1$ ). (b) ${\it\sigma}/{\it\chi}$ for a spheroid ( $a=c=1$ ). In this case, there is only one control parameter, i.e. ${\it\chi}=|{\it\chi}_{2}|$ .

4.3. Global method: Gledzer–Ponomarev (GP) polynomial perturbation analysis

Anticipating the results of the GP method, we first put forward that the base flow $\boldsymbol{u}_{0}$ is prone to inertial instabilities, similar to those that can grow upon uniform vorticity flows driven by precession, tides and longitudinal libration (Kerswell Reference Kerswell1993, Reference Kerswell1994; Lacaze, Le Gal & Le Dizès Reference Lacaze, Le Gal and Le Dizès2004; Le Bars et al. Reference Le Bars, Lacaze, Le Dizès, Le Gal and Rieutord2010; Cébron et al. Reference Cébron, Le Bars, Noir and Aurnou2012b ; Wu & Roberts Reference Wu and Roberts2013). Underlying these instabilities is a parametric resonance mechanism, whereby the strain imposed by the geometry couples two inertial modes of the unforced system, provided certain resonance conditions are met.

We follow roughly an approach set out by Tilgner (Reference Tilgner2007), and start our analysis by recasting (4.7) in a more condensed form,

where $\mathscr{L}$ is a time-dependent linear operator, which acts on $\boldsymbol{v}$ , and is characterized by a single frequency ${\it\omega}_{L}$ . It will now be instructive to consider perturbations $\boldsymbol{v}$ that are superpositions of two inertial modes of the cavity, i.e. 

Here $({\it\mu}_{A},\boldsymbol{V}_{A})$ and $({\it\mu}_{B},\boldsymbol{V}_{B})$ are eigenvalue–eigenmode pairs of the inviscid inertial mode problem

previously encountered in its ‘vorticity’ form (3.24). The modes $\boldsymbol{V}_{A,B}$ satisfy (2.12), and are mutually orthogonal, i.e.  $\langle \boldsymbol{V}_{B},\boldsymbol{V}_{A}\rangle ={\it\delta}_{AB}$ , with ${\it\delta}_{AB}$ the Kronecker delta. Substituting (4.21) into (4.20) and projection on $\boldsymbol{V}_{A}$ , $\boldsymbol{V}_{B}$ , respectively, leads to We note that all terms on the right-hand side of (4.24) and (4.25) are of order of magnitude ${\it\varepsilon}$ . This suggests that its solution will be characterized by the time scale ${\it\varepsilon}^{-1}$ . Thus, the rationale for the ansatz (4.21) is that it permits a time scale separation between the amplitudes $A$ and $B$ , characterized by the slow time scale ${\it\varepsilon}^{-1}$ , and the fast rotation scale of order of magnitude 1, captured by the inertial mode factors $\exp (\text{i}{\it\mu}_{A,B}t)$ . The fact that $A$ and $B$ only depend on the slow time scale ${\it\varepsilon}^{-1}$ suggests the following asymptotic expansion:

Substituting this into (4.24) and (4.25) yields, to order 0 in ${\it\varepsilon}$ , to ${\dot{A}}_{0}={\dot{B}}_{0}=0$ . To the next order in ${\it\varepsilon}$ , we obtain, after projection on $\boldsymbol{V}_{A}$ , $\boldsymbol{V}_{B}$ respectively:

The time-dependence of the diagonal terms in the above system is $\exp (\pm \text{i}{\it\omega}_{L}t)$ ; for the off-diagonal terms, it is $\exp (\pm \text{i}({\it\mu}_{A}-{\it\mu}_{B}\pm {\it\omega}_{L})t)$ . Resonant coupling between the two modes, tantamount to the occurrence of instability, requires that $A_{k},B_{k}$ can grow unbounded as $t\rightarrow \infty$ . This demands that

Since $-2\leqslant {\it\mu}_{A,B}\leqslant 2$ , couplings of the form (4.29) are limited to the libration frequency range $|{\it\omega}_{L}|\leqslant 4$ .

A further resonance condition can be found in the case of spheroidal geometries, for which the inertial modes $\boldsymbol{V}_{A,B}(\boldsymbol{r})$ are of the form $\boldsymbol{W}_{A,B}(r,z)\exp (\text{i}m_{A,B}{\it\varphi})$ , with $m_{A,B}\in \mathbb{Z}$ , and $(r,{\it\varphi},z)$ conventional cylindrical coordinates (Zhang, Liao & Earnshaw Reference Zhang, Liao and Earnshaw2004). After substituting this into (4.21), we find that the requirement that the right-hand side of (4.8) does not cancel, leads to the resonance condition

The above theoretical considerations are closely related to the GP method, which solves (4.7) in ellipsoidal enclosures. Devised originally by Gledzer & Ponomarev (Reference Gledzer and Ponomarev1977) and Lebovitz (Reference Lebovitz1989), it has since been employed by, amongst others, Kerswell (Reference Kerswell1993), Roberts & Wu (Reference Roberts and Wu2011) and Wu & Roberts (Reference Wu and Roberts2011, Reference Wu and Roberts2013). More precisely, this method restricts $\boldsymbol{v}$ to the vector space of perturbations that (i) are divergence-free, (ii) satisfy the impermeability condition, and (iii) have a vorticity field whose Cartesian components are polynomials that have degree less than $n$ in the Cartesian coordinates. For a given polynomial degree $n$ , this defines a vector space of finite dimension $N$ , and thus one can expand the perturbation in a finite series

In this expression, the vectors $\boldsymbol{V}_{j}^{\star }$ denote the basis vectors of the GP subspace of degree $n$ , and are Cartesian polynomials of degree $n$ or less. As a consequence of the particular nature of the Coriolis operator and of the uniform vorticity character of the base flow $\boldsymbol{u}_{0}$ , the expansion (4.31) reduces (4.7) to a closed system of $N$ linear ODEs for the functions $a_{1}(t),a_{2}(t),\dots ,a_{N}(t)$ . Since the coefficients of this system are time-periodic with period $T=2{\rm\pi}/{\it\omega}_{L}$ , one can solve it by means of Floquet analysis, or directly by brute numerical force; the last method was applied in this work. First, we integrate the ODEs in time from $t=0$ to $t=3000$ . Then, for $k=1,2,\dots ,N$ , we fit $\log |a_{k}(t)|$ with $\log |a_{k}(t)|=A_{k}+{\it\sigma}_{k}t$ , which gives the fit parameters $A_{k}$ and ${\it\sigma}_{k}$ , from which we determine the growth rate ${\it\sigma}=\max _{k}{\it\sigma}_{k}$ . Since this procedure is the source of small uncertainties in the measurement of ${\it\sigma}$ , we filter away values of ${\it\sigma}<0.02{\it\varepsilon}$ .

Recently, Vantieghem (Reference Vantieghem2014) established that one can always construct an orthogonal basis of inertial modes for the (finite-dimensional) GP subspaces of polynomial vector fields. Hence, the GP method describes the dynamics of the perturbation $\boldsymbol{v}_{2}$ in terms of a superposition of inertial modes, and thus allows identification of possible resonant couplings. This sets out the basic principles for a global analysis that is valid for ellipsoids of arbitrary shape.

To illustrate this, we solve (4.7) using quadratic ( $n=2$ ) and cubic ( $n=3$ ) GP subspaces for two configuration: (i) a spheroid with $a=b=1,c=0.5$ , (ii) a fully triaxial geometry with $a=1,b=1.5,c=0.5$ . We consider the libration frequency band between 0 and 4.5; the value of ${\it\varepsilon}$ is fixed at 0.1. Because the base flow $\boldsymbol{u}_{0}$ diverges as ${\it\omega}_{L}$ approaches the spin-over frequency $f$ , we exclude a small frequency window of width ${\rm\Delta}{\it\omega}=0.15$ around $f$ . In figures 12 and 13, we plot ${\it\sigma}/{\it\chi}$ for the spheroidal and triaxial geometry, respectively. These figures also provide the inertial mode spectrum of quadratic and cubic polynomial inertial modes. For the spheroidal geometry, we tabulate the azimuthal wave number $m$ of the mode as well.

Figure 12. GP analysis for the spheroidal geometry with $(a,b,c)=(1,1,0.5)$ using quadratic (red dashed) and cubic (blue) polynomial perturbations. Inviscid normalized growth rate ${\it\sigma}/{\it\chi}$ as a function of the libration frequency ${\it\omega}_{L}$ . A–F indicate the main resonance peaks and are given in the right-hand table. The black rectangle indicates the frequency band centred around the direct resonance at ${\it\omega}_{L}=1.6$ where the inviscid base flow $\boldsymbol{U}_{0}$ diverges.

Figure 13. GP analysis for the triaxial geometry with $(a,b,c)=(1,1.5,0.5)$ using quadratic (red dashed) and cubic (blue) polynomial perturbations. Inviscid normalized growth rate ${\it\sigma}/{\it\chi}$ as a function of the libration frequency ${\it\omega}_{L}$ . A–K indicate the main resonance peaks and are given in the right-hand table. The black rectangle indicates the frequency band centred around the direct resonance at ${\it\omega}_{L}=1.68$ where the inviscid base flow $\boldsymbol{U}_{0}$ diverges.

For both configurations, we observe a spectrum characterized by sharp peaks, which reflects the parametric-resonance-like nature of the instability. We see that the upper bound (4.14), i.e.  ${\it\sigma}\leqslant {\it\chi}=|{\it\chi}_{1}|+|{\it\chi}_{2}|$ , indeed holds. Further inspection of the main peaks (denoted by capital letters in figures 12 and 13) shows that these resonant frequencies are indeed associated with the criterion (4.29), where ${\it\mu}_{A}$ and ${\it\mu}_{B}$ correspond to values that are listed in the left-hand side of the respective figures. For the spheroidal geometry (figure 12), we see that every resonance is associated with an inertial mode pair that also satisfies the condition (4.30). The inertial mode spectrum of the spheroid is such that one resonance peak can sometimes be related to multiple pairs of modes. The triaxial configuration discussed in figure 13 exhibits more resonance peaks. In terms of mode interactions, one can argue intuitively that the equatorial deformation induces additional strain components that give birth to an extended set of mode couplings.

The results of the GP analysis thus give evidence that libration in latitude indeed drives parametric resonances that arise from inertial mode couplings. Having focused on quadratic and cubic inertial modes, we have (only) found a handful of resonances. More resonances would be expected, however, if we were to include higher-order polynomial GP bases. Indeed, since the inertial mode spectrum is dense (Greenspan Reference Greenspan1968), it is always possible to find two inertial mode frequencies that satisfy the criterion (4.29) for a given frequency ${\it\omega}_{L}$ . As such, we expect that instability can take place for any libration frequency $|{\it\omega}_{L}|\leqslant 4$ .

4.4. Direct numerical simulations

In this subsection, we compare the results of the LH and GP analysis against direct numerical simulations. In contrast to the preceding theoretical approaches, it is not feasible to carry out numerical simulations in the inviscid regime. Therefore, we have to reintroduce a viscous term into the perturbation equation (4.7). However, we replace the no-slip boundary condition with the stress-free condition

The rationale for this choice is that the no-slip condition may give rise to centrifugal-type instabilities if ${\it\varepsilon}$ is too high and/or $E$ too low (Noir et al. Reference Noir, Hemmerlin, Wicht, Baca and Aurnou2009; Sauret et al. Reference Sauret, Cébron, Le Bars and Le Dizès2012). Stress-free boundary conditions, on the other hand, allow isolation of inertial instabilities, such as the parametric resonances discussed above. A further advantage of this choice is that we avoid the computationally expensive task of resolving thin Ekman boundary layers. We keep the value of the Ekman number fixed at $E=10^{-4}$ , and we will focus on one particular geometry, characterized by the length of the semi-axes $(a,b,c)=(1,1.5,0.5)$ , i.e. the same triaxial geometry as studied in § 4.3. The deformation chosen is thus rather large. The motivation for this choice resides in the fact that our theoretical methods are not restricted to asymptotically small deformations, in contrast to other theoretical analyses (e.g. Kerswell Reference Kerswell1994; Le Bars et al. Reference Le Bars, Lacaze, Le Dizès, Le Gal and Rieutord2010).

In a first series of numerical simulations, we seek solutions of the linear perturbation equation (4.7). As initial condition we impose a random velocity field with a root-mean-square amplitude of $5\times 10^{-3}$ . We use the same numerical method as described in § 3.3. The computational mesh consists of uniformly distributed tetrahedral elements. In order to verify that our spatial resolution is adequate, two different mesh sizes have been considered, containing 1.3–10.5 million CVs respectively. The results reported hereafter are those obtained with the highest resolution; the quantities of interest do not differ more than ${\sim}1\,\%$ from those at lower resolution. This is shown in figure 14(a), where we compare the evolution of perturbation kinetic energy $E_{k}$ for ${\it\varepsilon}=0.3,{\it\omega}_{L}=2.1$ for the two different resolutions. Figure 14(b) illustrates the spatial structure of the growing instability for this parameter set.

Figure 14. Direct numerical solution of the (viscous) perturbation equation for ${\it\varepsilon}=0.3$ and ${\it\omega}_{L}=2.1$ . (a) Time series of the kinetic energy $E_{k}$ for a simulation using 1.3 million (dashed line) and 10.5 million (solid line) mesh elements. Also shown is an exponential fit to determine the growth rate ${\it\sigma}$ (dot-dashed line). (b) Snapshot of the velocity field $v_{y}$ in the planes $y=0$ and $z=0$ , illustrating the spatial structure of the growing instability.

We are now concerned with comparing the growth rates between the LH, GP and numerical approaches. To extract growth rates from our simulations, we proceed as follows. For a given value of ${\it\varepsilon}$ , we fit a curve of the form $\log E_{k}=A+2{\it\sigma}t$ to the numerical time series, as can be seen from figure 14(a). The fitting window is chosen large with respect to ${\it\omega}_{L}$ . The slope of this line gives a certain value for ${\it\sigma}({\it\varepsilon},{\it\omega}_{L},E=10^{-4})$ . Because we use stress-free boundary conditions and ${\it\varepsilon}$ is small compared to one, we can argue that (i) ${\it\sigma}$ is linear in ${\it\varepsilon}$ , and (ii) the viscously modified growth rate ${\it\sigma}_{visc}$ is related to the inviscid growth rate ${\it\sigma}$ by

with ${\it\alpha}$ independent of ${\it\varepsilon}$ . To eliminate the unknown coefficient ${\it\kappa}$ , we perform simulations for different values of ${\it\varepsilon}$ and report ${\it\alpha}={\it\sigma}/{\it\varepsilon}={\rm\Delta}{\it\sigma}_{visc}/{\rm\Delta}{\it\varepsilon}$ in figure 15.

Figure 15. Stability analysis in a triaxial configuration with $(a,b,c)=(1,1.5.0.5)$ . Comparison of growth rates resulting from the local LH analysis (solid line), the global GP analysis for quadratic (dashed line) and cubic (dot-dashed line) polynomial perturbations and direct numerical simulation (circles).

The values of the growth rate retrieved from the numerical simulations are located in the non-shaded area in figure 15, i.e. they are contained between the GP lower and LH upper bound. The numerics and the local LH analysis exhibit the same trend for the dependence of ${\it\sigma}$ on ${\it\omega}_{L}$ . However, the simulations yield growth rates that are about 35 % smaller than the ones found in the local LH analysis. This is consistent with the interpretation of the LH growth rates as upper bounds. Indeed, we recall that the LH approach does not constrain the perturbations to be bounded, and therefore provides a supremum rather than a maximum for the growth rate. It is worth noting that discrepancies of the same order of magnitude have been observed between the local and global stability analysis for flows driven by precession, libration and tides (Kerswell Reference Kerswell1993, Reference Kerswell1994; Cébron et al. Reference Cébron, Le Bars, Moutou and Le Gal2012a ).

According to figure 15 (or equivalently figure 13), the GP analysis does not predict instability for ${\it\omega}_{L}$ between 2 and roughly 2.9. Our simulations, on the other hand, show that instability occurs for almost the whole frequency window ${\it\omega}_{L}\in [2,4]$ . This implies that the numerically observed instabilities for ${\it\omega}_{L}=2.1,2.4,2.7$ should be associated with inertial modes of polynomial degree $n>3$ . This is indeed the case, as is illustrated in figure 16, where we show a snapshot of $u_{z}$ in the equatorial plane for ${\it\omega}_{L}=2.7$ and 3.3. The velocity profile for ${\it\omega}_{L}=2.7$ (figure 16 a) is too complex to be well captured by a quadratic or cubic polynomial. The corresponding profile for ${\it\omega}_{L}=2.7$ , on the other hand, has a spatially simpler structure.

Figure 16. Instantaneous profile of $u_{z}$ in the $z=0$ -plane for (a)  ${\it\varepsilon}=0.3,{\it\omega}_{L}=2.7$ and (b)  ${\it\varepsilon}=0.3,{\it\omega}_{L}=3.3$ .

In the earlier discussion we have shown that the instability arises by virtue of a parametric resonance mechanism coupling two inertial modes. We now investigate whether this can also be achieved in our numerical simulations. To this end, we consider time series of the velocity field at a single point $\boldsymbol{r}_{0}=(0.25,0.6,0.15)$ . In figure 17, we trace the frequency content $|S(\mathfrak{f})|$ of $\boldsymbol{s}(t)=\boldsymbol{v}(\boldsymbol{r}_{0},t)\exp (-{\it\sigma}t)$ . The factor $\exp (-{\it\sigma}t)$ is a compensation for the exponentially growing amplitude of $\boldsymbol{v}$ . $|S(\mathfrak{f})|$ is defined by $|S(\mathfrak{f})|^{2}=|S_{x}(\mathfrak{f})|^{2}+|S_{y}(\mathfrak{f})|^{2}+|S_{z}(\mathfrak{f})|^{2}$ , where $S_{x,y,z}(\mathfrak{f})$ denote the discrete Fourier transforms of the respective components of $\boldsymbol{s}(t)$ . We find that the spectrum is clearly dominated by two peaks, that moreover satisfy the resonance condition $\mathfrak{f}_{1}+\mathfrak{f}_{2}\approx {\it\omega}_{L}$ . This confirms that a parametric resonance mechanism is operating. In figure 17(a), we also observe smaller peaks around ${\it\omega}\approx 3.5,3.7$ . These are subharmonic effects that are due to the term $\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{U}_{0}+\boldsymbol{U}_{0}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{v}$ in (4.7) that couples the perturbation $\boldsymbol{v}$ with frequency content $f=1.104,1.289$ and the base flow oscillating at ${\it\omega}_{L}=2.4$ . A similar phenomenon is also present in figure 17(b).

Figure 17. Discrete Fourier transform of $\boldsymbol{s}(t)=\boldsymbol{v}(\boldsymbol{r}_{0},t)\exp (-{\it\sigma}t)$ at location $\boldsymbol{r}_{0}=(0.25,0.6,0.15)$ for (a)  ${\it\varepsilon}=0.2,{\it\omega}_{L}=2.4$ and (b)  ${\it\varepsilon}=0.2,{\it\omega}_{L}=3.3$ . Dashed lines indicate the driving frequency ${\it\omega}_{L}$ .

Finally, we also consider direct numerical simulations of the non-linear equations (2.1) and (2.2) with ${\it\varepsilon}=0.3$ and ${\it\omega}_{L}=2.7$ , supplemented with the stress-free boundary condition (4.32). In figure 18(a), we display the evolution in time of the kinetic energy, as well as the energy related to the non-uniform vorticity component of the flow. This quantity is obtained as the difference between $\boldsymbol{u}$ and the projection of $\boldsymbol{u}$ on the subspace of uniform vorticity flows. Starting from the solution (3.4), (3.17) and (3.18), we see that a weak non-uniform vorticity flow immediately emerges, which is related to a weak boundary layer flow required to match the stress-free boundary condition. However, until $t\approx 50$ , the flow remains essentially of uniform vorticity, as can be seen from figure 19(a), which shows an instantaneous velocity profile of $u_{y}$ in the planes $x=0$ and $y=0$ at $t=35$ . For $t\gtrsim 50$ , the non-uniform vorticity component undergoes exponential growth over several orders of magnitude. This gives evidence of the presence of instability. Furthermore, we see that the growth rate of the the non-uniform vorticity component in the non-linear simulation is in good agreement with the results from the corresponding simulation of the linear equation (4.7). In figure 18(b), we depict an instantaneous profile of the non-uniform vorticity component of $u_{z}$ in the plane $z=0$ at $t=125$ ; we find that its spatial structure is similar to its counterpart from the linear stability analysis, displayed in figure 16(a). The growth of the non-uniform vorticity component also affects the structure of the total velocity field. This can be seen from figure 19(b), where we see slightly rippled isolines of $u_{y}$ at $t=125$ . Eventually, the instability saturates when the non-uniform vorticity kinetic energy is approximately of the same order of magnitude as the total kinetic energy. The flow has now completely lost its uniform vorticity character, as can be seen from figure 18(c,d).

Figure 18. (a) Time series of non-linear simulations for ${\it\varepsilon}=0.3$ and ${\it\omega}_{L}=2.7$ , showing the total kinetic energy (solid thick line) and the kinetic energy of the non-uniform vorticity flow (solid thin line). The dashed line shows the evolution of perturbation kinetic energy in a corresponding linear simulation. (b) Snapshot of the $z$ -component of the non-uniform vorticity flow in the plane $z=0$ at $t=125$ .

Figure 19. Non-linear simulations of latitudinally libration driven flow for ${\it\varepsilon}=0.2$ and ${\it\omega}_{L}=2.7$ . Instantaneous profiles of $u_{y}$ in the planes $x=0$ and $y=0$ at (a)  $t=35$ , (b)  $t=125$ , (c)  $t=215$ and (d)  $t=275$ .