2. Governing equations

Figure 1. Schematic of the cube librating around the rotation axis $\boldsymbol {\varOmega }(t)$ , with the meridional plane $x=0$ shown in blue and the equatorial plane $y=-z$ shown in green.

A cube of side length $L$ is completely filled with an incompressible fluid of kinematic viscosity $\nu$ . It rotates around an axis passing through the midpoints of opposite edges at a mean rate $\varOmega$ that is modulated harmonically at a frequency $\sigma$ with amplitude $\delta$ . The system is non-dimensionalised using $L$ as the length scale and $1/\varOmega$ as the time scale, and is described in terms of a non-dimensional Cartesian coordinate system $\boldsymbol {x}=(x,y,z)\in [-0.5,0.5]^3$ that is fixed in the cube, with the origin at the centre. The non-dimensional angular velocity is

(2.1) \begin{equation} \boldsymbol {\varOmega }(t)=[1+ Ro \cos (2\omega t)]{\boldsymbol {\varOmega }_0}\,, \quad {\boldsymbol {\varOmega }_0}=(0,1,1)/\sqrt {2}, \end{equation}

where the non-dimensional libration frequency is

(2.2) \begin{equation} 2\omega ={\sigma }/{\varOmega } \end{equation}

and the relative amplitude is the Rossby number

(2.3) \begin{equation} Ro ={\delta }/{\varOmega }. \end{equation}

Figure 1 shows a schematic of the system. The two walls at $x=\pm 0.5$ are parallel to the rotation axis $\boldsymbol {\varOmega }_0$ and the four remaining walls are inclined at angles $\pm {45}^{\circ }$ relative to $\boldsymbol {\varOmega }_0$ . Four of the edges are orthogonal to $\boldsymbol {\varOmega }_0$ . The two bisected by the rotation axis are termed the north and south polar edges and the vertices at their ends are polar vertices. The other two edges are the equatorial edges with equatorial vertices at their ends. The two vertices at $(x,y,z)=(\pm 0.5,\mp 0.5,\pm 0.5)$ are the eastern equatorial vertices and the two vertices at $(\mp 0.5,\mp 0.5,\pm 0.5)$ are the western equatorial vertices. The remaining eight edges are inclined at angles $\pm {45}^{\circ }$ relative to $\boldsymbol {\varOmega }_0$ ; these are tropical edges, each with a polar and an equatorial vertex at their ends. The depth is the distance parallel to the rotation axis between the top and bottom inclined walls; this varies from $\sqrt {2}$ between the polar edges to 0 at the equatorial edges. Also shown in the schematic are the meridional plane $x=0$ in blue and the equatorial plane $y=-z$ in green.

In the frame of reference fixed to the cube, the non-dimensional Navier–Stokes equations for the relative velocity $\boldsymbol u$ are

(2.4) \begin{equation} \frac {\partial \boldsymbol {u}}{\partial t} + (\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {u} + \boldsymbol {\nabla }p -{E}\nabla ^2\boldsymbol {u} +2\boldsymbol {\varOmega }\! \times \! \boldsymbol {u} = -\frac {{\rm d}\boldsymbol {\varOmega }}{{\rm d}t}\! \times \! \boldsymbol {x}, \quad \boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}=0, \end{equation}

where $p$ is the reduced pressure incorporating the centrifugal force and

(2.5) \begin{equation} {E}=\nu /(\varOmega \, L^2) \end{equation}

is the Ekman number. The last term on the left-hand side of (2.4) is the Coriolis force that includes a time-harmonic component from the librational forcing. The right-hand-side term is the time-harmonic Euler force also corresponding to the librational forcing. The boundary conditions are no-slip, ${\boldsymbol u}={\boldsymbol 0}$ , on all walls of the cube.

In the small-forcing regime, $0\lt Ro \ll 1$ , the magnitude of the response flow is of the order of $ Ro$ , and it is convenient to rescale the velocity and pressure by introducing

(2.6) \begin{equation} {\boldsymbol v}={\boldsymbol u}/ Ro\quad \text {and}\quad q=p/ Ro , \end{equation}

with corresponding vorticity ${\boldsymbol \nabla }\times {\boldsymbol v}$ and enstrophy density $({\boldsymbol \nabla }\times {\boldsymbol v})^2$ . The system (2.4) then becomes

(2.7) \begin{equation} \dfrac {\partial \boldsymbol {v}}{\partial t} + Ro (\boldsymbol {v}\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {v} +\boldsymbol {\nabla }q -{E}\boldsymbol {\nabla }^2\boldsymbol {v} +2[1+ Ro \cos (2\omega t)]\boldsymbol {\varOmega }_0 \times \boldsymbol {v} = \boldsymbol {f},\quad \boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {v}=0, \end{equation}

where the scaled Euler force is

(2.8) \begin{equation} \boldsymbol {f}=2\omega \sin (2\omega t)\,\boldsymbol {\varOmega }_0\times \boldsymbol {x}. \end{equation}

The systems (2.4) and (2.7) are only equivalent for $ Ro \ne 0$ . Indeed, for $ Ro =0$ (2.4) remains nonlinear, but unforced, with trivial solution $(\boldsymbol {u},p)=(\boldsymbol {0},0)$ , while (2.7) reduces to the linear forced system

(2.9) \begin{equation} \frac {\partial \boldsymbol {v}_0}{\partial t} +\boldsymbol {\nabla }q_0 - {E}\boldsymbol {\nabla }^2\boldsymbol {v}_0 +2\boldsymbol {\varOmega }_0\times \boldsymbol {v}_0 = \boldsymbol {f}, \quad \boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {v}_0=0, \end{equation}

with non-trivial solutions $(\boldsymbol {v}_0,q_0)=\left .(\boldsymbol {v},q)\right |_{Ro =0}$ . The system (2.9) provides a convenient way to explore the small- $ Ro$ limit in an asymptotic sense, namely $(\boldsymbol {u},p)\approx Ro (\boldsymbol {v}_0,q_0)$ for sufficiently small $ Ro$ .

3. Numerics

The numerical scheme for solving the governing equations (2.7) is based on a consistent splitting scheme for the Navier–Stokes equations with no-slip boundary conditions, initially proposed in Guermond & Shen (Reference Guermond and Shen2003). The current scheme incorporates the higher-order temporal accuracy and improved conditional stability for both velocity and pressure introduced in Wu et al. (Reference Wu, Huang and Shen2022a ). Briefly, the momentum equation is discretised in time using a third-order backwards difference formula, BDF3, with the Euler force and the viscous terms treated implicitly and the nonlinear, Coriolis force and pressure terms evaluated explicitly via third-order extrapolation. An updated pressure is obtained by integrating against a gradient field, with the viscous term evaluated using the alternate formulation ${\boldsymbol \nabla }^2{\boldsymbol v}=-{\boldsymbol \nabla }\times ({\boldsymbol \nabla }\times {\boldsymbol v})$ , which improves the accuracy of the computed pressure. The treatment of each stage leads to a decoupled system of Helmholtz-type equations with constant coefficients for the velocity and a Poisson equation with constant coefficients for the pressure. These are then spatially discretised using Legendre polynomials of degree $\ell$ in each of the three spatial directions for the velocity components and of degree $\ell -2$ for the pressure. The nonlinear advection term is evaluated in physical space via projection from spectral space, and a zero-mean pressure condition is imposed by setting the constant basis pressure coefficient to $0$ . This results in a symmetric positive-definite banded system for the pressure, which is solved efficiently via diagonalisation.

In this study, we only present results for small $\omega =0.05$ for which the foliation patterned topographic Rossby waves are prominent (Wu et al. Reference Wu, Welfert and Lopez2022b ). The spatial and temporal resolutions needed for numerical stability and accuracy are much larger for small $\omega$ than for $\omega \sim 1$ at the same $E$ and $ Ro$ . Furthermore, these requirements increase for both decreasing $E$ and increasing $ Ro$ . Table 1 lists the Legendre polynomial degree, $\ell$ , and number of time steps per libration period, $n_t$ , needed in the simulations for $\omega =0.05$ and $ Ro \lesssim 10{E}$ for $E\in [10^{-9},10^{-6}]$ .

Table 1. Legendre polynomial degree, $\ell$ , and number of time steps per libration period, $n_t$ , used for $\omega =0.05$ , $ Ro \lesssim 10{E}$ and $E$ as indicated.

The same numerical approach is used to solve the linear system (2.9), but without computing the nonlinear or time-dependent Coriolis terms. Not having to compute those terms represents significant computational savings.

4. The linear forced viscous response

For $E\gt 0$ , solutions to the homogeneous linear system, (2.9) with $\boldsymbol {f}=\boldsymbol {0}$ , viscously decay in time, $\boldsymbol {v}_0\to \boldsymbol {0}$ (Greenspan Reference Greenspan1968), whereas the librationally forced linear system, (2.9) with $\boldsymbol {f}\ne \boldsymbol {0}$ , drives a synchronous temporally harmonic response. This response may include cavity modes as well as shear layers that for small $E$ may be described as beams comprised of circularly polarised waves. The fate of the wavebeams in the container at $E=0$ can be estimated via ray tracing using the linear inviscid unforced equations together with the reflection laws for circularly polarised waves. A problem with this approach, however, is that with this system being unforced, it says nothing about where the wavebeams originate from in the container, and being linear it says nothing about their initial energy. Considerations from small- $E$ simulations of the librationally forced linear system help inform where the beams are emitted from. This is detailed in Wu et al. (Reference Wu, Welfert and Lopez2022b ) and Welfert et al. (Reference Welfert, Lopez and Wu2023). Briefly, for the cube orientation under consideration and small $\omega$ , conic sheets of beams are emitted from polar vertices and planar sheets of beams are emitted from polar and tropical edges. These reflect on walls and eventually focus onto the equatorial edges and vertices. Figure 2 shows the traces of beams on the surface of the cube, in the equatorial plane $y=-z$ and in the meridional plane $x=0$ that are respectively emitted from polar vertices, polar edges and tropical edges, as well as their superposition. The colours correspond to the enstrophy density level of the circularly polarised wave relative to its initial enstrophy density at emission, which is arbitrary. In the absence of viscous effects, $E=0$ , the enstrophy density of a beam remains the same between consecutive reflections, but changes upon reflections at the oblique walls due to focusing.

Figure 2. Traces of wavebeams, originating from various locations as indicated, on the cube surface and their intersections with the equatorial plane $y=-z$ and the meridional plane $x=0$ , obtained via ray tracing at $\omega =0.05$ . The colours correspond to the enstrophy density of a wavebeam relative to its initial enstrophy density at emission.

Figure 3. Snapshots of the linear response flows at $\omega =0.05$ and $E$ as indicated, showing $E(\boldsymbol {\nabla }\times {\boldsymbol v}_0)^2$ and $q_0$ on the surface and in the equatorial and meridional planes. The cessation of focusing in a vicinity $d$ of equatorial edges is indicated by a black line in the $(\boldsymbol {\nabla }\times {\boldsymbol v}_0)^2$ equatorial planes. Supplementary movie 1 $(a,b)$ animates the surface and equatorial plane views for $E=10^{-7}$ and $10^{-9}$ over one period.

We now consider the librationally forced linear viscous system (2.9) at $\omega =0.05$ . Figure 3 shows snapshots of the enstrophy density $({\boldsymbol \nabla }\times {\boldsymbol v}_0)^2$ , scaled by $E$ , and the pressure $q_0$ on the cube surface and the equatorial and meridional planes for a range of $E$ . The enstrophy density on the equatorial and meridional planes for $E=10^{-9}$ has much in common with the linear inviscid ray tracing results shown in figure 2. The linear inviscid wavebeams correspond to the thin viscous vortex sheets that are emitted from the same edges and vertices. One difference is that for the linear inviscid ray tracing there is focusing along the entire equatorial edges, whereas the non-zero- $E$ flows have intense peaks in $(\boldsymbol {\nabla }\times {\boldsymbol v}_0)^2$ localised near the eastern equatorial vertices; $(\boldsymbol {\nabla }\times {\boldsymbol v}_0)^2$ vanishes at the edges due to the no-slip boundary conditions. This difference is particularly evident on the surface plots. Moreover, the viscous shear layers do not focus all the way to the equatorial edges, particularly for increasing $E$ : their focusing is viscously suppressed. Comparing the viscous results shown in the equatorial plane in figure 3 with the inviscid wavebeam analysis shown in figure 2, a loss of left–right and up–down symmetry is apparent. This is a consequence of the viscous term breaking the time-reversal symmetry of the inviscid problem (Wu et al. Reference Wu, Welfert and Lopez2022b , equation (4.7)).

The distance $d$ from the equatorial edge to where focusing appears to cease is indicated by a thin line in the plot of $(\boldsymbol {\nabla }\times {\boldsymbol v}_0)^2$ in the equatorial plane (second row of figure 3). The distinct difference in flow dynamics either side of this line is better appreciated by viewing the animations in supplementary movie 1. This distance decreases with $E$ , with $d=500{E}^{1/2}$ ( $d=0.5$ for ${E}=10^{-6}$ ). For decreasing $E$ , the wavebeams associated with the thinner viscous shear layers have larger wavenumbers and their focusing reaches closer towards the equatorial edges. As they focus, the successive reflections on the oblique walls occur closer to each other, eventually resulting in viscous cross-annihilation (Stanaway, Shariff & Hussain Reference Stanaway, Shariff and Hussain1988; Terrington, Hourigan & Thompson Reference Terrington, Hourigan and Thompson2021) between incident and reflected vortex sheets.

Figure 4. Zoom-ins ( $\times 4$ ) around the equatorial vertex $(0.5,-0.5,0.5)$ in the equatorial plane shown in figure 3 of $E(\boldsymbol {\nabla }\times {\boldsymbol v}_0)^2$ (top) and $q_0$ (bottom).

In the regions near the equatorial edges devoid of the focusing shear layers, a very different flow behaviour is evident, corresponding to the foliated progressive waves mentioned in the Introduction. The focusing shear layers are related to the circularly polarised waves analysed in the linear inviscid ray tracing and hence their $(\boldsymbol {\nabla }\times {\boldsymbol v}_0)^2$ traces remain fixed in space for all time. In contrast, the arms of the foliated waves emerge from the equatorial edges, with one end more-or-less fixed at the eastern equatorial vertices. They sweep out an arc in a retrograde fashion relative to the rotation axis and terminate at the walls $x=\pm 0.5$ . This is illustrated in supplementary movie 1 $(a,b)$ where $(\boldsymbol {\nabla }\times {\boldsymbol v}_0)^2$ for ${E}=10^{-7}$ and $10^{-9}$ is animated over one libration period. The stationary sharp structures seen in the equatorial plane, particularly evident in the scaled enstrophy density at $E=10^{-9}$ (top-right panel of figure 4), coincide very well with the linear inviscid ray tracings in figure 2. Superimposed with these are the retrograde foliated waves. For increasing $E$ , the stationary peaks corresponding to the beams are absent in the regions near the equatorial edges and the foliated progressive waves remain the only contribution to the flow within the distance $d$ of the equatorial edges.

To further explore these foliated waves, figure 3 also shows the pressure $q_0$ on the cube surface and in the equatorial and meridional planes. The pressure plots are particularly useful as the pressure contribution from the foliated waves completely dominates the contribution from the focusing shear layers. The meridional plane views at $x=0$ as well as at the wall $x=0.5$ show that the foliation pattern is quasi-independent of the axial direction. This is true for both $(\boldsymbol {\nabla }\times {\boldsymbol v}_0)^2$ and $q_0$ . Each alternating arm of the foliation has opposite signed $q_0$ , and each half libration period an arm is mapped to the location of its neighbour in the retrograde direction, as shown in supplementary movie 1 $(a,b)$ .

Figure 4 presents plots of $(\boldsymbol {\nabla }\times {\boldsymbol v}_0)^2$ and $q_0$ in the equatorial plane shown in figure 3 zoomed in by a factor of 4 near the eastern equatorial vertex at $(x,y,z)=(0.5,-0.5,0.5)$ . For the largest ${E}=10^{-6}$ , the inertial wavebeams are absent and for diminishing $E$ their presence becomes more evident in $(\boldsymbol {\nabla }\times {\boldsymbol v}_0)^2$ . The zoom-ins make it clearer that the spatial shape of $(\boldsymbol {\nabla }\times {\boldsymbol v}_0)^2$ associated with the foliated waves remains invariant with $E$ , while the wavebeams at smaller $E$ focus closer to the equatorial edge and their superpositions with the foliated waves lead to an interference pattern. On the other hand, $q_0$ is almost entirely due to the foliated waves and it is even clearer that their structure is independent of $E$ , at least for $E\leqslant 10^{-6}$ . All these features point to the foliated progressive waves being forced topographic Rossby waves.

Figure 5. Variation with $E$ of $\|(\boldsymbol {\nabla }\times {\boldsymbol v}_0)^2\|_\infty$ and $\|q_0\|_\infty$ .

From figures 3 and 4, it is evident that $(\boldsymbol {\nabla }\times {\boldsymbol v}_0)^2$ scales with $E^{-1}$ almost everywhere throughout the cube, except in the vicinity of the eastern equatorial vertices where it grows much larger with decreasing $E$ . In order to quantify this localised growth, we use an $\infty$ -norm, which is the maximum of a quantity $Q$ throughout the cube over one libration period:

(4.1) \begin{equation} \|{Q}\|_\infty = \max \{ |{{Q}({\boldsymbol x},\,t)|: {\boldsymbol x} \in [-0.5,0.5]^{3},\,\, t \in [t_a,\,t_a+\pi /\omega ]}\}. \end{equation}

The time $t_a$ is a time after the flow has settled to a synchronous periodic state, typically corresponding to several dozen libration periods, with the number of periods increasing with decreasing $E$ .

Figure 5 shows how the $\infty$ -norms of $(\boldsymbol {\nabla }\times {\boldsymbol v}_0)^2$ and $q_0$ vary with $E$ . The pressure scales weakly with $E$ , with $\|q_0\|_\infty \propto {E}^{-1/6}$ . In contrast, $\|(\boldsymbol {\nabla }\times {\boldsymbol v}_0)^2\|_\infty \propto {E}^{-8/3}$ grows much faster with decreasing $E$ . The scaling of $(\boldsymbol {\nabla }\times {\boldsymbol v}_0)^2$ with $E^{-1}$ almost everywhere throughout the cube, evident in figure 3, is to be expected from the development of thinner and more intense boundary layers and shear layers with decreasing $E$ , whereas the much faster $E^{-8/3}$ scaling of $\|(\boldsymbol {\nabla }\times {\boldsymbol v}_0)^2\|_\infty$ is indicative of developing singularities near the two eastern equatorial vertices about which the foliation patterns pivot.

5. Nonlinear responses at small $ Ro$

Figure 6. Snapshots of the symmetric synchronous response flows at $\omega =0.05$ and $E$ and $ Ro$ as indicated, showing $E(\boldsymbol {\nabla }\times {\boldsymbol v})^2$ and $q$ on the surface, equatorial and meridional planes. Supplementary movie 1 $(c,d)$ animates the surface and equatorial plane views for $E=10^{-7}$ and $10^{-9}$ over one period.

The development towards singularities in the enstrophy density strongly localised at the two eastern equatorial vertices with decreasing $E$ raises the question of how the nonlinearities in the Navier–Stokes equations, so far neglected, affect their $E^{-8/3}$ scaling. Figure 6 shows snapshots of the enstrophy density $(\boldsymbol {\nabla }\times {\boldsymbol v})^2$ and pressure $q$ on the cube surface and on the equatorial and meridional planes for the same $E$ as in figure 3, but for $ Ro = Ro _{\textit{max}}({E})$ , the largest values of $ Ro$ for each $E$ for which we were able to compute a stable synchronous response flow. The spatial and temporal resolutions needed at $ Ro = Ro _{\textit{max}}({E})$ are listed in table 2. For the largest $E=10^{-6}$ , further increasing $ Ro$ leads to an instability introducing quasi-periodic behaviour. For smaller $E$ , increasing $ Ro$ leads to temporal and spatial resolution requirements that made simulations of quasi-periodic responses prohibitive. This is due, as demonstrated below, to the developing singularities as $E$ is decreased having large nonlinear effects even for very small non-zero $ Ro$ .

Table 2. Legendre polynomial degree, $\ell$ , and number of time steps per libration period, $n_t$ , used for $\omega =0.05$ , $ Ro = Ro _{\textit{max}}({E})$ and $E$ as indicated.

Comparing figures 3 and 6, as well as supplementary movies 1 $(a,b)$ and 1 $(c,d)$ , there is virtually no difference between the enstrophy density and pressure at $ Ro =0$ and the larger $ Ro$ , except in the vicinity of the developing singularities near the two eastern equatorial vertices. To quantify these differences the synchronous flow $({\boldsymbol v}, q)$ , computed from (2.7) at a given $E$ and $0\lt Ro \ll 1$ , is decomposed as

(5.1) \begin{equation} {\boldsymbol v}={\boldsymbol v}_0 +\widetilde {\boldsymbol v}, \quad q=q_0+\widetilde {q}, \end{equation}

where $({\boldsymbol v}_0, q_0)$ is the linear temporally harmonic flow (i.e. it has zero mean and no higher harmonics) computed from (2.9) at the same $E$ and $(\widetilde {\boldsymbol v},\widetilde {q})$ is time periodic containing all the nonlinear contributions from $({\boldsymbol v}, q)$ , i.e. the mean $(\overline {\boldsymbol v},\overline {q})$ and all higher harmonics.

In order to minimise numerical errors in evaluating $(\widetilde {\boldsymbol v},\widetilde {q})$ , $({\boldsymbol v}_0, q_0)$ is recomputed here using the same higher resolution used for $({\boldsymbol v},q)$ , listed in table 2, compared with that from table 1, which was sufficient to obtain well-resolved $({\boldsymbol v}_0,q_0)$ .

Figure 7. (a–f) Variation with $ Ro$ of the enstrophy and pressure $\infty$ -norms of the various components in (5.1), as well as for the mean flow, for $E$ as indicated.

Figure 8. Snapshots of (top rows) the enstrophy density associated with the mean flow $\overline {\boldsymbol v}$ , nonlinear contribution $\widetilde {\boldsymbol v}$ , and (bottom rows) the mean pressure $\overline {q}$ and the pressure nonlinear contribution $\widetilde {q}$ , in the equatorial plane for $E$ and $ Ro$ as indicated. The $d\times d$ square corner insets, with $d$ as in figure 3, delimit the zoom-in regions shown beneath each plot. Supplementary movie 2 animates the zoom-ins of $({\boldsymbol \nabla }\times \widetilde {\boldsymbol v})^2$ and $\widetilde {q}$ over one libration period.

The $\infty$ -norms of the various components in (5.1) are shown in figure 7 $(a,b)$ . The limits as $ Ro \to 0$ of $\|(\boldsymbol {\nabla }\times {\boldsymbol v})^2\|_\infty$ and $\|q\|_\infty$ are the values of $\|(\boldsymbol {\nabla }\times {\boldsymbol v}_0)^2\|_\infty$ and $\|q\|_\infty$ shown in figure 5; $\|(\boldsymbol {\nabla }\times {\boldsymbol v})^2\|_\infty$ and $\|q\|_\infty$ do not vary very much with $ Ro$ until $ Ro \approx Ro _{\textit{max}}$ . On the other hand, $\|({\boldsymbol \nabla }\times \widetilde {\boldsymbol v})^2\|_\infty \propto Ro ^2$ and $\|\widetilde {q}\|_\infty \propto Ro$ , as expected from the quadratic nonlinearity in the Navier–Stokes equations. Figure 7 $(c,d)$ rescales $\|({\boldsymbol \nabla }\times \widetilde {\boldsymbol v})^2\|_\infty$ and $\|\widetilde {q}\|_\infty$ showing that the nonlinear contributions in the enstrophy density scale with $E^{-16/3} Ro ^2$ and the nonlinear contributions in the pressure scale with $E^{-3/2} Ro$ . These scalings are used in the colourmap scalings of $({\boldsymbol \nabla }\times \widetilde {\boldsymbol v})^2$ and $\widetilde {q}$ in figure 8, as well as the enstrophy density of the mean flow $({\boldsymbol \nabla }\times \overline {\boldsymbol v})^2$ and the mean pressure $\overline {q}$ , which have the same scalings but slightly lower magnitudes, confirming that, as $E$ decreases, the weakly nonlinear response flow is captured everywhere in the container by the linear component $({\boldsymbol v}_0,q_0)$ , except in regions around the two eastern equatorial vertices of size scaling with $d\sim {E}^{1/2}$ . This is not surprising since $ Ro$ is quite small and the enstrophy density is dominated almost everywhere (i.e. everywhere except within a distance $d$ of the equatorial edges) by the focusing shear layers which at these small values of $E$ and $ Ro$ are essentially circularly polarised wavebeams, with $(\boldsymbol {v\cdot \nabla })\boldsymbol {v}=0$ for individual beams.

Zoom-ins of the enstrophy densities associated with the nonlinear contributions $(\widetilde {\boldsymbol v},\widetilde {q})$ in a $d\times d$ square region of the eastern equatorial vertices, indicated in the plots of $({\boldsymbol \nabla }\times \widetilde {\boldsymbol v})^2$ and $\widetilde {q}$ in figure 8, are also shown in that figure and animated over one libration period in supplementary movie 2. In order to avoid pixelation effects resulting from repeated zoom-ins, the spectral solution was projected onto finer physical grids for each panel. The similar structure of the last three panels for each quantity confirms both the spatial extent, $\mathcal {O}({E}^{1/2})$ , and magnitudes, $\mathcal {O}({E}^{-16/3} Ro ^2)$ , for the enstrophy density and $\mathcal {O}({E}^{-3/2} Ro )$ for the pressure, corresponding to the nonlinear contributions.

The above $( Ro ,{E})$ -dependent scalings raise the question of how small $ Ro$ should be compared with $E$ in order for the linear approximation $(\boldsymbol {v},q)\simeq (\boldsymbol {v}_0,q_0)$ to be globally valid. From the scalings inferred from figures 5 and 7, we have $\|({\boldsymbol \nabla }\times {\boldsymbol v}_0)^2\|_\infty \sim {E}^{-8/3}$ and $\|({\boldsymbol \nabla }\times \widetilde {\boldsymbol v})^2\|_\infty \sim {E}^{-16/3} Ro ^2$ , so that $\|({\boldsymbol \nabla }\times \widetilde {\boldsymbol v})^2\|_\infty /\|({\boldsymbol \nabla }\times {\boldsymbol v}_0)^2\|_\infty \sim {E}^{-8/3} Ro ^2$ . For this to be smaller than some acceptable tolerance $\mathcal {T}^2$ , we must have

(5.2) \begin{equation} Ro \lt {E}^{4/3}\mathcal {T}. \end{equation}

The scalings from figures 5 and 7 also show $\|q_0\|_\infty \sim {E}^{-1/6}$ and $\|\widetilde {q}\|_\infty \sim {E}^{-3/2} Ro$ , so that $\|\widetilde {q}\|_\infty /\|q_0\|_\infty \sim {E}^{-4/3} Ro$ , which is smaller than the tolerance $\mathcal {T}$ under the same condition (5.2). So, for the small-forcing-amplitude flows to be globally well approximated by the linear flows, keeping $ Ro$ approximately smaller than $E^{4/3}$ is necessary.