2. Governing equations

A cuboid with square cross-section of side length $L$ and height ${A{\kern-4pt}R} L$, where ${A{\kern-4pt}R}$ is the height-to-base aspect ratio, is completely filled with an incompressible fluid of kinematic viscosity $\nu$. The container is rotating around an axis passing through its centre and the midpoints of opposite horizontal edges at a mean rate $\varOmega$ that is modulated harmonically at a frequency $\sigma$ with amplitude $\delta \varOmega$. 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]^2\times [-0.5{A{\kern-4pt}R},0.5{A{\kern-4pt}R} ]$ that is fixed in the container, with the origin at the centre. The corresponding non-dimensional velocity field is $\boldsymbol {u}=(u,v,w)$. In these coordinates, the non-dimensional angular velocity is

(2.1a,b)\begin{equation} \boldsymbol{\varOmega}(t)=[1+{Ro} \cos(2\omega t)] {\boldsymbol{\varOmega}_0}, \quad {\boldsymbol{\varOmega}_0}=(0,1,{A{\kern-4pt}R})/\sqrt{1+{A{\kern-4pt}R}^2}, \end{equation}

where $2\omega =\sigma /\varOmega >0$ is the non-dimensional libration frequency, and the relative amplitude ${Ro}=\delta \varOmega /\varOmega$ is the Rossby number. Figure 1 shows a schematic of the system for the three aspect ratios ${A{\kern-4pt}R} =1/2$, 1 and 2 considered in the DNS. The two walls of the container at $x=\pm 0.5$ are parallel to the axis of rotation. The four remaining walls are inclined at angles $\alpha =\text {arccot}({A{\kern-4pt}R} ^{\pm 1})$ relative to the rotation axis $\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 other two are the tropical edges. The remaining eight edges are inclined at angles ${\rm \pi} /2-\alpha$ relative to $\boldsymbol {\varOmega }_0$.

Figure 1.

Schematic of the cuboid librating around its axis of rotation $\boldsymbol {\varOmega }(t)$, for aspect ratios ${A{\kern-4pt}R}$ as indicated; the meridional plane $x=0$ is shown in blue.

The non-inertial frame of reference attached to the librating cuboid introduces Coriolis and Euler body forces into the (non-dimensional) governing equations:

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

where ${E}=\nu /(\varOmega L^2)$ is the Ekman number, and $p$ is the reduced pressure that incorporates the centrifugal force. In this frame of reference, the no-slip boundary conditions are $\boldsymbol {u}=\boldsymbol {0}$ on all six walls of the container. For small but non-zero ${Ro}$, it is convenient to introduce $\boldsymbol {v}=\boldsymbol {u}/{Ro}$ and $p_r=p/{Ro}$, as in Lopez et al. (Reference Lopez, Shen, Welfert and Wu2022). Using (2.1a,b), the governing equations then become

(2.3)\begin{align} \frac{\partial \boldsymbol{v}}{\partial t} +{Ro}\,(\boldsymbol{v}\boldsymbol{\cdot}{\boldsymbol{\nabla}})\boldsymbol{v} +2\boldsymbol{\varOmega} \times \boldsymbol{v} - 2\omega\sin(2\omega t)\,\boldsymbol{\varOmega}_0 \times \boldsymbol{x}={-}{\boldsymbol{\nabla}}p_r +{E}\,{\nabla}^2\boldsymbol{v}, \quad {\boldsymbol{\nabla}}\boldsymbol{\cdot}\boldsymbol{v}=0. \end{align}

For $0<{E}={Ro}\ll 1$, solutions $\boldsymbol {v}$ and $p_r$ of (2.3) are typically synchronous limit cycles that respect the centrosymmetry of the governing equations and boundary conditions,

(2.4)\begin{equation} \mathcal{C}: [\boldsymbol{v},p_r](\boldsymbol{x},t) \mapsto[-\boldsymbol{v},p_r](-\boldsymbol{x},t), \end{equation}

corresponding to a reflection through the origin. In the linear inviscid setting (${E}={Ro}=0$), the system (2.3) is formally reduced to a non-homogeneous linear system whose forcing is cognisant of the libration of the container, i.e. the Euler force persists.

The numerical scheme and code used here is essentially the same as that used in closely related problems (Lopez et al. Reference Lopez, Shen, Welfert and Wu2022; Wu et al. Reference Wu, Welfert and Lopez2022b). The governing equations (2.3) are discretised using a third-order linearly implicit scheme in time, and a spectral scheme in space. At each time step, the velocity components are updated via a third-order backwards difference formula scheme, with the viscous terms as well as the Coriolis and Euler forces treated implicitly, while the pressure gradient and nonlinear terms are evaluated via third-order extrapolation. This leads to a coupled system of elliptic equations with variable coefficients. The pressure is then updated via the solution of a Poisson problem posed in the cuboid, including the boundary, leading to a constant-coefficient Poisson equation. The number of time steps $\delta t$ used per librational forcing period is $200$ for all Ekman and Rossby numbers considered here.

The elliptic and Poisson problems are discretised in each of the three spatial directions using Legendre polynomials of degree $M$ for the velocity components, and degree $M-2$ for the pressure, with $M$ ranging from $M=100$ for ${E}=10^{-5}$ to $M=650$ for ${E}=10^{-8}$. Specifically, the velocity and pressure at times $t=n \delta t$ are expanded as

(2.5a,b)\begin{equation} \boldsymbol{v}^n(x,y,z) = \sum_{i,j,k=2}^{M} a^n_{ijk}\,\psi_{ijk}(x,y,z), \quad p_r^n (x,y,z) = \sum_{i,j,k=0}^{M-2} b^n_{ijk}\,\psi_{ijk}(x,y,z), \end{equation}

in terms of basis functions

(2.6)\begin{equation} \psi_{ijk}(x,y,z) = \psi_i(2x)\,\psi_j(2y)\,\psi_k(2z/{A{\kern-4pt}R}), \end{equation}

with

(2.7) \begin{equation} \psi_i(\xi)= r_i[\mathfrak{L}_{i-2}(\xi)-\mathfrak{L}_i(\xi)], \quad i\ge0, \ {-1}\le\xi\le1, \end{equation}

where $\mathfrak {L}_k$ is the Legendre polynomial of degree $k$, with $\mathfrak {L}_k(\xi )=0$ for $k<0$, and normalising constants $r_i$. The property $\psi _i(\pm 1)=0$ for $i\ge 2$ guarantees $\boldsymbol {v}^n=\boldsymbol {0}$ on the walls of the cuboid, while the lower degree of $p_r^n$ compared to $\boldsymbol {v}^n$ guarantees the well-posedness of the resulting discrete problems. The resulting banded linear systems for the velocity and pressure updates are then diagonalised and solved. An optional scalar auxiliary variable procedure developed in Wu, Huang & Shen (Reference Wu, Huang and Shen2022a) introduces a third-order correction of the velocity that guarantees unconditional stability and improves the robustness of the scheme, in particular during transient dynamics.

3. Overview of ray tracing analysis

Ray tracing is normally performed on the unforced homogeneous system, neglecting the Euler force, but the forcing determines the locations from which beams are emitted. Viscous interactions between the oscillatory boundary layers driven by the small-amplitude librational forcing lead to inertial wave beams being emitted into the cuboid from vertices and/or edges. Which vertices and/or edges emit depends on ${A{\kern-4pt}R}$ and $\omega$. In general, a vertex may emit along a double cone of directions with its apex at the emission point forming a conical sheet. However, a point on an edge may emit in only at most four directions because of continuity requirements with beams emitted from neighbouring points and tangentiality with conical sheets emitted from its endpoints, if any. The inviscid theory presented in Welfert et al. (Reference Welfert, Lopez and Wu2023) shows that a beam emitted from a vertex or an edge tends to be attracted, after many reflections, towards a plane with constant $x$, parallel to the axis of rotation. Ultimately, the beam focuses either onto a corner point in this planar cross-section, which is a point on an edge orthogonal to the axis of rotation, or onto a closed-loop interior attractor in this $x$-plane. In the inviscid setting, this $x$-plane depends non-smoothly on the location of the point or direction of emission because of boundary singularities (edges and vertices) creating jumps in the wall normal direction. However, in practice, these are regularised by viscous effects.

A closed-loop attractor is characterised by the number of reflections $m$ at the wall $z=0.5{A{\kern-4pt}R}$ and $n$ at the wall $y=0.5$, with an equal number of reflections at walls $z=-0.5{A{\kern-4pt}R}$ and $y=-0.5$ by the centrosymmetry. Such an attractor is denoted $m:n$. Figure 2 shows the inviscid $(\omega ^2,{A{\kern-4pt}R} )$-regime diagram for wave beams contained completely in the meridional plane $x=0$. The region of existence of interior attractors is delimited by the curves $\omega ^2 = 1/(1+{A{\kern-4pt}R} ^{\pm 2})$, along which either the walls $y=\pm 0.5$ or the walls $z=\pm 0.5{A{\kern-4pt}R}$ have critical reflection slope. Inside this region, all reflections are supercritical, whereas outside the region, reflections on the walls at $y=\pm 0.5$ and $z=\pm 0.5{A{\kern-4pt}R}$ are subcritical. For ${A{\kern-4pt}R} =1$, reflections are not supercritical for any $\omega ^2$, and they are critical only for $\omega ^2=1/2$. Outside these regions, but for $\omega ^2<1$, inertial wave beams focus onto edges and/or vertices. For ${A{\kern-4pt}R} =1$, there are no interior attractors, and as $\omega ^2$ is varied across $1/2$, there is a switch in the edges and vertices to which wave beams focus. Note that in the limits ${A{\kern-4pt}R} ^{\pm 1}\to \infty$, the walls of the container tend to be either parallel or orthogonal to the rotation axis, yet the range of existence of interior attractors extends to $0<\omega ^2<1$.

Figure 2.

Inviscid $(\omega ^2,{A{\kern-4pt}R} )$-regime diagram for wave beams in the plane $x=0$: the bold black curves $\omega ^2=1/(1+{A{\kern-4pt}R} ^{\pm 2})$ delimit regions of differing criticality of wave beam reflections on the walls at $y=\pm 0.5$ or $z=\pm 0.5{A{\kern-4pt}R}$ and focusing of the beams onto $m:n$ or $n:m$ attractors. The red stars correspond to cases discussed in this study. The figure has been adapted from Welfert et al. (Reference Welfert, Lopez and Wu2023).

Although the attractors in a given $x$-plane in the DNS and VEBA associated with ${A{\kern-4pt}R} =a$ and ${A{\kern-4pt}R} =1/a$ can be mapped to each other via a $90^\circ$ rotation and rescaling, other details differ away from the attractors due to conical shears originating from the vertices of the cuboids. These differences are due to the spanwise widths of the two configurations. A perfect match exists between the $1\times 1\times a$ and $1/a \times 1 \times 1/a$, rather than $1\times a\times 1/a$, configurations, albeit with a different effective Ekman number; see Appendix A for details.

4. Direct numerical simulations of simple attractors

We consider simple attractors, either a point attractor at edges for ${A{\kern-4pt}R} =1$, or a 1 : 1 attractor for ${A{\kern-4pt}R} =1/2$ and ${A{\kern-4pt}R} =2$, all at the same libration frequency corresponding to $\omega =0.55$ ($\omega ^2=0.3025$); these cases are indicated by the red stars in figure 2. They are examined via DNS over the range ${E}={Ro}\in [10^{-8}, 10^{-5}]$. The spatial and temporal resolutions used for these DNS are listed in table 1. For all cases, ${E}$ and ${Ro}$ are sufficiently small that the responses to librational forcing are $\mathcal {C}$-invariant synchronous limit cycles. It is important to reduce both ${E}$ and ${Ro}$ in lock-step; fixing $0<{Ro}\ll 1$ and reducing only ${E}$ will lead to the fixed libration amplitude ${Ro}$ being sufficiently large for a sufficiently small ${E}$ to trigger flow instability. On the other hand, fixing $0<{E}\ll 1$ and reducing only ${Ro}$ ultimately results in a synchronous forced response flow whose magnitude $\|\boldsymbol {u}\|$ scales with ${Ro}$, so that $\|\boldsymbol {v}\|$ remains ${O}(1)$ and tends towards being temporally harmonic, corresponding to the solution of (2.3) with ${Ro}=0$. The interest in studying the response as both ${Ro}$ and ${E}$ are reduced in problems involving geometric focusing onto attractors is that even as ${Ro}$ is reduced, the gradients in the flow velocity grow as ${E}$ is reduced and the combined result is not a priori obvious.

Table 1.

Degree $M$ of the Legendre polynomials in each of the three spatial directions used in the DNS with $\omega =0.55$ and ${E}={Ro}$, and ${A{\kern-4pt}R}$ as indicated.

It is informative to put into perspective the ranges of ${E}$ and ${Ro}$ considered here compared to current state-of-the-art laboratory experiments investigating inertial wave attractors (e.g. Brunet, Dauxois & Cortet Reference Brunet, Dauxois and Cortet2019; Boury et al. Reference Boury, Sibgatullin, Ermanyuk, Shmakova, Odier, Joubaud, Maas and Dauxois2021). The smallest Ekman numbers achieved experimentally are ${E}\gtrsim 2 \times 10^{-6}$, corresponding to a background rotation $\varOmega \approx 2\ {\rm rad}\ {\rm s}^{-1}$ in a container of length scale $L\approx 0.5$ m filled with water of kinematic viscosity $\nu \approx 10^{-6}\ {\rm m}^2\ {\rm s}^{-1}$. To achieve ${E}= 10^{-8}$ with water at the same background rotation rate, the container length scale needs to be $L\approx 7$ m. As for the Rossby number ${Ro}$ characterising the forcing amplitude in the experiments, this is no smaller than approximately $7 \times 10^{-3}$ due to signal-to-noise ratio issues; this is small enough for the flows to be strongly affected by the Coriolis force, but is nevertheless large enough to correspond to a developed nonlinear regime.

The flows are characterised primarily by using their enstrophy density $\boldsymbol {\omega }^2=|\boldsymbol {\nabla }\times \boldsymbol {v}|^2$, its mean

(4.1)\begin{equation} \overline{\boldsymbol{\omega}^2}=\frac{1}{\tau} \int_{t^*}^{t^*+\tau} \boldsymbol{\omega}^2\, {\rm d} t, \end{equation}

and its standard deviation

(4.2)\begin{equation} \boldsymbol{\omega}^2_{SD} = \sqrt{\frac{1}{\tau} \int_{t^*}^{t^*+\tau} (\boldsymbol{\omega}^2-\overline{\boldsymbol{\omega}^2})^2\, {\rm d} t}, \end{equation}

where $\tau ={\rm \pi} /\omega$ is the libration period, and $t^*$ is a time by which the response flow is $\tau$-periodic (typically of the order of $10^{3}\tau$).

At aspect ratios ${A{\kern-4pt}R} =1/2$ and 2, in the limit with both ${E}\to 0$ and ${Ro}\to 0$, $1:1$ attractors exist in the range $\omega ^2\in (0.2,0.64)$ (Welfert et al. Reference Welfert, Lopez and Wu2023). Figure 3 illustrates these cases for $\omega =0.55$ and ${E}={Ro}\in [10^{-8}, 10^{-5}]$, along with the ${A{\kern-4pt}R} =1$ case for which focusing is to point attractors on the tropical edges. The figure shows $\overline {\boldsymbol {\omega }^2}$ and $\boldsymbol {\omega }^2_{SD}$ in the meridional plane $x=0$, seen from the positive $x$ direction, which summarise the synchronous response flow. Supplementary movie 1 animates $\boldsymbol {\omega }^2$ for the ${E}={Ro}= 10^{-8}$ cases over one libration period $\tau ={\rm \pi} /\omega$. For ${A{\kern-4pt}R} =1/2$ and 2, the traces of wave beams emanating from the four corners (the north and south poles and the midpoints of the tropical edges) focus onto an interior attractor that becomes thinner with deceasing ${E}$, while for ${A{\kern-4pt}R} =1$, the focusing is into the tropical corners. For all cases, $\overline {\boldsymbol {\omega }^2}$ and $\boldsymbol {\omega }^2_{SD}$ on the attractor seem to grow slightly faster than ${E}^{-1}$, while the response in the interior away from the attractor remains essentially independent of ${E}$ as ${E}$ is decreased. In the $x=0$ plane, VEBA of wave beams emitted from the four corners into the plane is also performed and compared to the DNS results. The VEBA in the last column of figure 3(a) captures well how the beams (depicted in blue) focus towards the attractor (depicted in red). Traces of the primary conic beams emitted from the vertices of the cuboid and of planar beams emitted from edges at $x=\pm 0.5$ intersecting the plane $x=0$ are also shown in light blue in the VEBA; these are also evident in the DNS. Secondary traces due to reflections of these beams, which are also visible in the DNS, are not included in VEBA (see Appendix B for further details concerning these beams).

Figure 3.

(a) Mean enstrophy density $\overline {\boldsymbol {\omega }^2}$ and (b) standard deviation of the enstrophy density $\sqrt {2}\boldsymbol {\omega }^2_{SD}$ in the meridional plane $x=0$ (where the rotation axis is vertical through the top and bottom corners), for ${A{\kern-4pt}R}$ and ${E}={Ro}$ as indicated. The last column in (a) has the VEBA results showing beams emitted from the four corners in blue and the attractor in red; traces of the primary conic beams emitted from the vertices of the cuboid and planar beams emitted from edges at $x=\pm 0.5$ are in light blue. Secondary traces due to reflections of these beams, which are visible in the DNS, are not included in VEBA. The black lines through the origin in panels in the first column of (a), $(0,y,0)$ and $(0,0,z)$, are used to plot profiles of $\overline {\boldsymbol {\omega }^2}$ and $\boldsymbol {\omega }^2_{SD}$ in figure 4. Supplementary movie 1 (available at https://doi.org/10.1017/jfm.2023.772) animates $\boldsymbol {\omega }^2$ for the ${E}={Ro}= 10^{-8}$ cases over one libration period $\tau ={\rm \pi} / \omega$.

Away from the attractor region, the mean and standard deviation of the enstrophy density are related by $\overline {\boldsymbol {\omega }^2}\approx \sqrt {2}\boldsymbol {\omega }^2_{SD}$. This is a strong indication that in regions where this relationship holds, the periodic flow is very close to being temporally harmonic (see Appendix C for details). On the other hand, along the attractor for decreasing ${E}$ and ${Ro}$, $\overline {\boldsymbol {\omega }^2}$ becomes increasingly larger than $\sqrt {2}\boldsymbol {\omega }^2_{SD}$, indicating that the oscillations are increasingly non-harmonic and nonlinear. To quantify this further, figure 4 shows profiles of $\overline {\boldsymbol {\omega }^2}$ and $\sqrt {2}\boldsymbol {\omega }^2_{SD}$ along the line segments from $(x,y,z)=(0,-0.5,0)$ to $(0,0,0)$ to $(0,0,0.5{A{\kern-4pt}R} )$, which traverse the two branches of the attractor region. These profiles provide some insight into how the shape and strength of the response across the viscous shear layers associated with the attractor regions vary with decreasing ${E}$ and ${Ro}$. At larger ${E}$, the wave beams emanating from the corners of the $x=0$ meridional plane can hardly be differentiated from their long-term location, with a broad diffuse transverse profile of the shear. With decreasing ${E}$, the emitted shears become thinner, making it possible to distinguish them clearly from their ultimate location (the attractor) following multiple reflections on the walls of the cuboid. As a result, the transverse profiles become more spatially oscillatory, with each local maximum representing the crossing of successive reflections of the shears. Figure 4(b) shows zoom-ins of the attractor regions ($y\sim -0.4$ and $z\sim 0.8$) for the ${A{\kern-4pt}R} =2$ cases, with additional intermediate values of ${E}={Ro}$. These clearly show how, for the largest ${E}= 10^{-5}$, the shear layer is a fusion of the beams emitted from the corners and all of their reflections. By ${E}\sim 10^{-6}$, the layer separates into two with one more intense than the other. Further decreasing ${E}$, the more intense layer further splits at ${E}\sim 10^{-7}$, again with one of them becoming more intense. By ${E}\sim 10^{-8}$, the weaker of the two layers that split at ${E}\sim 10^{-6}$ is about to split. The various split shear layers that appear with decreasing ${E}$ follow the VEBA trajectories of the beams emitted from the corners as they reflect off the walls and focus towards the attractor; the most intense of the split layers converges towards the VEBA attractor, which in the linear inviscid setting of VEBA is a Delta-like distribution.

Figure 4.

(a) Profiles of $\overline {\boldsymbol {\omega }^2}$ and $\sqrt {2}\boldsymbol {\omega }^2_{SD}$ along the line segments shown in figure 3, $(x,y,z)=(0,y,0)$ and $(x,y,z)=(0,0,z)$, for ${A{\kern-4pt}R}$ and ${E}={Ro}$ as indicated. (b) Zoom-ins for the ${A{\kern-4pt}R} =2$ cases in the attractor regions, with additional intermediate values of ${E}$ and ${Ro}$ as indicated. The vertical lines correspond to the locations as determined by VEBA of the attractor (in red) and primary (before any reflection, farther away from the red line) or secondary (after one wall reflection, closer to the red line) vortex sheets originating from tropical or polar edges (in blue); compare with the VEBA panel from figure 3(a). Note the more rapid growth of $\overline {\boldsymbol {\omega }^2}$ at the attractor.

Returning to figure 4(a), the secondary peaks due to the conical shears originating from the vertices of the cuboid, located at $z\approx 0.35$ for ${A{\kern-4pt}R} =2$ and $y\approx -0.25$ for ${A{\kern-4pt}R} =1/2$, grow, without splitting, at a much smaller rate with decreasing ${E}$ than the shears in the attractor region. These secondary peaks correspond to the primary intersection of the conical shears with the $x=0$ plane, and show a self-similar convergence with decreasing ${E}$ towards a Delta-like distribution.

Figure 4(a) also shows that while the $\boldsymbol {\omega }^2_{SD}$ profiles have essentially the same shape as the $\overline {\boldsymbol {\omega }^2}$ profiles, the ratio $\overline {\boldsymbol {\omega }^2}/\boldsymbol {\omega }^2_{SD}$ becomes increasingly larger with decreasing ${E}$ in the attractor region. However, this ratio remains close to $\sqrt {2}$ away from the attractor region. This is a clear indication of the oscillations in the shear layers in the attractor region becoming increasingly nonlinear and non-harmonic with decreasing ${E}$ and ${Ro}$.

Figure 5 shows $\overline {\boldsymbol {\omega }^2}$ and $\boldsymbol {\omega }^2_{SD}$ along the attractor in the $x=0$ plane, using the $(y,z)$ location of the attractor predicted from VEBA, which is parametrised by the arc length along the attractor $s$. The mean enstrophy along the attractor scales as $\overline {\boldsymbol {\omega }^2}\sim {E}^{-4/3}$, whereas the standard deviation scales as $\boldsymbol {\omega }^2_{SD}\sim {E}^{-3/4}$. While both grow with decreasing ${E}$, the mean grows much faster, an indication that the shear flow in the attractor region becomes more nonlinear as ${E}$ and ${Ro}$ are reduced. These scalings, however, do not apply in the localised regions where the attractor reflects at the walls. In these regions, $s\approx 0.94$ and $s\approx 2.48$, the spikes in both $\overline {\boldsymbol {\omega }^2}$ and $\boldsymbol {\omega }^2_{SD}$ are several orders of magnitude larger than in other locations along the attractor. Figure 6 zooms in on these locations, and shows that $\overline {\boldsymbol {\omega }^2}\sim \boldsymbol {\omega }^2_{SD}\sim {E}^{-5/3}$ at the peaks. Furthermore, the arc length distance over which this very different scaling holds scales as ${E}^{1/2}$, suggesting that the ${E}^{-5/3}$ scaling may be due to interactions between the shear in the attractor and the oscillatory Ekman layers at the walls.

Figure 5.

Scaled $\overline {\boldsymbol {\omega }^2}$ and $\boldsymbol {\omega }^2_{SD}$ as functions of arc length $s$ along the attractor (localised using VEBA), for ${A{\kern-4pt}R} =2$ and ${E}={Ro}$ as indicated.

Figure 6.

Rescaled profiles of $\overline {\boldsymbol {\omega }^2}$ and $\boldsymbol {\omega }^2_{SD}$ from figure 5, showing (a,b) zoom-ins around $s=s_1=(\sqrt {1-\omega ^2}-0.75\omega )\sqrt {5} \approx 0.9451$ together with the same zoom-ins in terms of a scaled arc length centred at $s_1$, and (c,d) zoom-ins around $s=s_2=(\sqrt {1-\omega ^2}+0.5\omega )\sqrt {5}\approx 2.4824$ together with the same zoom-ins in terms of a scaled arc length centred at $s_2$, showing self-similarity of the peaks associated with the reflections on the short and long walls, respectively. The total length of the attractor path is $L=2s_2\approx 4.9648$ (Welfert et al. Reference Welfert, Lopez and Wu2023, (2.34)).

The $\overline {\boldsymbol {\omega }^2}$ and $\boldsymbol {\omega }^2_{SD}$ boundary layers at the walls $y=-0.5$ and $z=0.5{A{\kern-4pt}R}$ are very thin and difficult to discern in figure 4. Instead, boundary layer profiles are plotted in figure 7, which illustrates how their thickness and intensity vary with decreasing ${E}={Ro}$ along the line segments shown in figure 3, $(x,y,z)=(0,y,0)$ and $(x,y,z)=(0,0,z)$, which are away from locations where the attractor reflects on the walls. In the oscillatory Ekman boundary layers, away from where any interior shear layers reflect, we find that the relationship $\overline {\boldsymbol {\omega }^2}\approx \sqrt {2}\boldsymbol {\omega }^2_{SD}$ holds, so that the oscillations are very nearly temporally harmonic, both $\overline {\boldsymbol {\omega }^2}$ and $\boldsymbol {\omega }^2_{SD}$ scale with ${E}^{-1}$, and the boundary layer thickness for both scales with ${E}^{1/2}$.

Figure 7.

Profiles of scaled $\overline {\boldsymbol {\omega }^2}$ and $\boldsymbol {\omega }^2_{SD}$ in the boundary layer along the line segments shown in figure 3, $(x,y,z)=(0,y,0)$ and $(x,y,z)=(0,0,z)$, for ${A{\kern-4pt}R}$ and ${E}={Ro}$ as indicated.

Figure 8 shows $\overline {\boldsymbol {\omega }^2}$ along the cuboid walls, with distance parametrised by the arc length $s\in [0,6]$ (only $s\in [0,3]$ is shown as $s\in [3,6]$ is the same by $\mathcal {C}$ symmetry), in the meridional plane $x=0$ for ${A{\kern-4pt}R} =2$ and several values of ${E}={Ro}$. Almost everywhere on the wall, $\overline {\boldsymbol {\omega }^2}\sim {E}^{-1}$. At the corners $s=0$ and 1, $\overline {\boldsymbol {\omega }^2}=0$, and in the regions where the shear layers in the attractor regions reflect at the walls at $s\approx 0.2$ and $s\approx 1.5$, the scaling for the peak changes to $\overline {\boldsymbol {\omega }^2}\sim {E}^{-5/3}$, as noted earlier. It is of interest to consider $\overline {\boldsymbol {\omega }^2}$ and $\boldsymbol {\omega }^2_{SD}$ on the whole cuboid boundary. These are shown in figure 9 for the three aspect ratios ${A{\kern-4pt}R} =1/2$, 1 and 2. Both 3-D perspective views and unfolded views of the walls are illustrated (only three walls are shown; the other three walls are the same by symmetry). The $x=0.5$ wall, which is parallel to the rotation axis (indicated by a black arrow), shows the footprints of the attractor shear layers depicted in the $x=0$ planes in figure 3. On the walls oblique to the rotation axis, footprints of the conical beams emitted from the vertices are evident. On these walls for ${A{\kern-4pt}R} =1/2$ and 2, $\overline {\boldsymbol {\omega }^2}$ and $\boldsymbol {\omega }^2_{SD}$ have strong peaks in the regions where the shear layers in the attractor region reflect at the walls, whereas for ${A{\kern-4pt}R} =1$, the peaks are at the tropical edges $(y,z)=(\pm 0.5, \pm 0.5{A{\kern-4pt}R} )$ to which beams focus. There is considerable variation in the $x$ direction for all cases. Also, away from localised regions where reflections occur, $\overline {\boldsymbol {\omega }^2}\sim \sqrt {2}\boldsymbol {\omega }^2_{SD}$, indicating that the oscillatory Ekman layers are nearly harmonic. Supplementary movie 2 animates $\boldsymbol {\omega }^2$ on the surface over one libration period $\tau ={\rm \pi} /\omega$, and shows that the nearly harmonic oscillations are dominated by progressive waves that sweep through the regions, with more intense $\overline {\boldsymbol {\omega }^2}$ corresponding to the shear layer footprints on $x=0.5$ and the regions where the shear layers reflect on the $y=0.5$ and $z=0.5{A{\kern-4pt}R}$ walls.

Figure 8.

Scaled $\overline {\boldsymbol {\omega }^2}$ along the boundary at $x=0$, for ${A{\kern-4pt}R} =2$ and ${E}={Ro}$ as indicated.

Figure 9.

Perspective and unfolded views of (a) $\overline {\boldsymbol {\omega }^2}$ and (b) $\sqrt {2}\boldsymbol {\omega }^2_{SD}$ on the cuboid surfaces for ${E}={Ro}= 10^{-8}$ and ${A{\kern-4pt}R}$ as indicated. Supplementary movie 2 animates $\boldsymbol {\omega }^2$ on the surface over a libration period $\tau ={\rm \pi} /\omega$.

Now consider the structure of the attractor shear layer. Figure 10 shows $\overline {\boldsymbol {\omega }^2}$ and $\sqrt {2}\boldsymbol {\omega }^2_{SD}$ on two branches of the attractor (the distributions on the other two branches are the same with $x\to -x$ by the $\mathcal {C}$ symmetry), whose location is predicted from VEBA, for ${A{\kern-4pt}R} =2$ and ${E}={Ro}\in [10^{-8}, 10^{-5}]$. The general trends on these quantities along the attractor intersection with the meridional plane $x=0$ reported in figures 5 and 6 are borne out, but there is considerable variation in the $x$ direction, which becomes larger with decreasing ${E}$ and ${Ro}$. Also evident are signatures of where the conical beams emitted from vertices, and a number of their reflections, intersect the attractor shear layer. For the larger ${E}$, these are evident faintly in the mean enstrophy, but as ${E}$ is reduced, the amplification of $\overline {\boldsymbol {\omega }^2}$ swamps these completely. However, as the enstrophy in the attractor shear layer tends towards a steady state with decreasing ${E}$, the conical beam intersections become much more evident in $\sqrt {2}\boldsymbol {\omega }^2_{SD}$, as they are locally the dominant contribution to the flow oscillations. The last row in figure 10 shows the intersections with the attractor of the wave beams emitted from edges and vertices, and their subsequent reflections, as determined via VEBA. This captures the fine details in the DNS $\sqrt {2}\boldsymbol {\omega }^2_{SD}$ on the attractor, especially for the smaller ${E}$ and ${Ro}$. This is clear evidence that in the DNS, the wave beams are oscillatory, but as they eventually focus onto the attractor, they fuse into a mean shear flow that locally overwhelms the inertial oscillations.

Figure 10.

Illustrations of $\overline {\boldsymbol {\omega }^2}$ and $\sqrt {2}\boldsymbol {\omega }^2_{SD}$ on two branches of the attractor (the other two are the $\mathcal {C}$ conjugates), for ${A{\kern-4pt}R} =2$ and ${E}={Ro}$ as indicated, as well as a 3-D perspective of $\overline {\boldsymbol {\omega }^2}$ on the attractor for the case ${E}= 10^{-8}$. The bottom row presents intersections with the attractor of the wave beams emitted from edges and vertices, and their subsequent reflections, as determined via VEBA. Supplementary movie 3 includes animations of $\boldsymbol {\omega }^2$ over one libration period.

Figures 5, 6, 9 and 10 show that there is a considerable increase in $\overline {\boldsymbol {\omega }^2}$ as the attractor shear reflects on the walls of the cuboid, and this gain appears to increase as ${E}$ and ${Ro}$ are reduced. In the linear inviscid regime, VEBA gives what the gains are due solely to geometric focusing at the reflections. For ${A{\kern-4pt}R} =2$, the gains $g_1$ at reflections on the long walls at $y=\pm 0.5$, and $g_2$ on the short walls at $z=\pm 1$, are (Welfert et al. Reference Welfert, Lopez and Wu2023, (2.5) and (2.7))

(4.3)\begin{equation} g_1 = \frac{\sin(\theta+\alpha)}{\sin(\theta-\alpha)} =\frac{{A{\kern-4pt}R}\omega+\sqrt{1-\omega^2}}{{A{\kern-4pt}R}\omega-\sqrt{1-\omega^2}}\approx7.31 \end{equation}

and

(4.4)\begin{equation} g_2 = \frac{\cos(\theta+\alpha)}{\cos(\theta-\alpha)} =\frac{\omega+{A{\kern-4pt}R}\sqrt{1-\omega^2}}{\omega-{A{\kern-4pt}R}\sqrt{1-\omega^2}}\approx1.98. \end{equation}

The enstrophy density scales with the fourth power of the gain, and hence increases by approximate factors $2851$ on the long walls and $15.43$ on the short walls. Figure 11 shows how the gains in $\overline {\boldsymbol {\omega }^2}$ at these wall reflections vary with decreasing ${E}$ and ${Ro}$ in the DNS. These gains were determined simply on the $x=0$ plane, and clearly figures 9 and 10 show that there is considerable variation in $\overline {\boldsymbol {\omega }^2}$ with $x$ at the reflections. Nevertheless, the point to be made is that the $\overline {\boldsymbol {\omega }^2}$ gains in the DNS with decreasing ${E}$ and ${Ro}$ surpass the geometric gains from VEBA, indicating that nonlinear processes are playing an increasing role.

Figure 11.

Gains in $\overline {\boldsymbol {\omega }^2}$ at reflections at $x=0$ on the long and short walls as functions of ${E}$ and ${Ro}$ for ${A{\kern-4pt}R} =2$.

The primary nonlinear processes contributing to the enhanced gains in $\overline {\boldsymbol {\omega }^2}$ at the reflections are most likely to be vortex tilting and stretching. These nonlinear mechanisms are inherently 3-D and play a prominent role in hydrodynamics in regimes where the inertial time scale is many orders of magnitude smaller than the viscous time scale, in our case ${E}\ll 1$ (Frisch Reference Frisch1996; Majda & Bertozzi Reference Majda and Bertozzi2002; Davidson Reference Davidson2013). Figure 12 shows the magnitudes of the mean and standard deviation of the vortex stretching/tilting term, $|\overline {(\boldsymbol {\omega }\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {v}}|$ and $|[(\boldsymbol {\omega }\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {v}]_{SD}|$, on a long and short attractor branch; supplementary movie 3 includes animations of $|{(\boldsymbol {\omega }\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {v}}|$ over one libration period, and figure 13 shows their profiles at the meridional plane $x=0$. As with the enstrophy, $|{(\boldsymbol {\omega }\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {v}}|$ on the attractor tends towards a steady state with decreasing ${E}$ and ${Ro}$, unsteadiness becoming localised to regions where the attractor reflects on the walls and to a lesser extent where beams emitted from edges and vertices intersect the attractor shear layer. As was the case with $\sqrt {2}\boldsymbol {\omega }^2_{SD}$ on the attractor, $|[(\boldsymbol {\omega }\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {v}]_{SD}|$ also matches the VEBA results shown in figure 10.

Figure 12.

Illustrations of $|\overline {(\boldsymbol {\omega }\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {v}}|$ and $|[(\boldsymbol {\omega }\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {v}]_{SD}|$ on two branches of the attractor (the other two are the $\mathcal {C}$ conjugates), for ${A{\kern-4pt}R} =2$ and ${E}={Ro}$ as indicated. Supplementary movie 3 includes animations of $|(\boldsymbol {\omega }\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {v}|$ over one libration period.

Figure 13.

Scaled $|\overline {(\boldsymbol {\omega }\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {v}}|$ and $|[(\boldsymbol {\omega }\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {v}]_{SD}|$ as functions of arc length $s$ along the attractor (localised using VEBA), for ${A{\kern-4pt}R} =2$ and ${E}={Ro}$ as indicated.

With nonlinearities not vanishing in the attractor with decreasing ${E}$ and ${Ro}$, it is interesting to determine to what extent $\boldsymbol {v}$ and $\boldsymbol {\omega }$ are aligned on the attractor, i.e. to determine $\alpha =\arccos (\hat {\boldsymbol {\omega }}\boldsymbol {\cdot }\hat {\boldsymbol {v}})$. In the linear inviscid theory of a single circularly polarised inertial wave, either $\alpha =0^\circ$, $\boldsymbol {\omega }\boldsymbol {\cdot } \boldsymbol {v}>0$ and the group velocity is from south to north, or $\alpha =180^\circ$, $\boldsymbol {\omega }\boldsymbol {\cdot } \boldsymbol {v}<0$ and the group velocity is from north to south. In open space, such a Beltrami flow has $(\boldsymbol {v}\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {v}=0$. Appendix D shows that there is in general no expectation of alignment between $\boldsymbol {v}$ and $\boldsymbol {\omega }$ for superpositions of circularly polarised waves, such as on an attractor. Figure 14 shows the mean and standard deviation of $\alpha$ along the attractor for ${A{\kern-4pt}R} =2$ and various ${E}={Ro}$. As ${E}$ and ${Ro}$ are reduced from $10^{-5}$, the flow on the attractor tends to become more Beltrami-like, but as these are reduced below $10^{-7}$, this trend is reversed. This departure from Beltrami-like flow for increasingly smaller ${E}$ and ${Ro}$ is consistent with the flow in the attractor shear layer region becoming increasingly more nonlinear, and the thickness of the shear layers becoming increasingly thinner as they focus onto Delta-like distributions. Figure 15 shows $\bar {\alpha }$ and $\alpha _{SD}$ on the attractor; the 3-D perspective view for the ${E}={Ro}= 10^{-8}$ case illustrates how $\bar {\alpha }\approx 180^\circ$ almost everywhere on the branches with the group velocity going north, and $\bar {\alpha }\approx 0^\circ$ almost everywhere on the branches with the group velocity going south.

Figure 14.

Variations of the mean and standard deviation of $\alpha =\arccos (\hat {\boldsymbol {v}}\boldsymbol {\cdot }\hat {\boldsymbol {\omega }})$ with arc length $s$ along the attractor (localised using VEBA) at $x=0$, for ${A{\kern-4pt}R} =2$ and ${E}={Ro}$ as indicated. Only two branches of the attractor are shown; on the other two branches, the mean angle is $180^\circ -\bar {\alpha }$ and the standard deviation is the same.

Figure 15.

Mean and standard deviation of the angle between velocity and vorticity, $\alpha =\arccos (\hat {\boldsymbol {\omega }}\boldsymbol {\cdot }\hat {\boldsymbol {v}})$, on two branches of the attractor, for ${A{\kern-4pt}R} =2$ and ${E}={Ro}$ as indicated, as well as a 3-D perspective of $\bar {\alpha }$ on the attractor for the case ${E}= 10^{-8}$. Supplementary movie 3 includes animations of $\alpha$ over one libration period.

5. Zeroth law: ${E}\mathcal {E}\ne 0$ as ${E}\to 0$

In this section, the global energy balance is considered for decreasing ${E}$ and ${Ro}$. Using $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {v}=0$ and the identities

(5.1a,b)\begin{equation} (\boldsymbol{v}\boldsymbol{\cdot}{\boldsymbol{\nabla}})\boldsymbol{v}= {\boldsymbol{\nabla}}\left(\frac{\boldsymbol{v}^2}{2}\right)+ \boldsymbol{\omega}\times\boldsymbol{v} \quad \textrm{and}\quad {\nabla}^2\boldsymbol{v}= {\boldsymbol{\nabla}}({\boldsymbol{\nabla}}\boldsymbol{\cdot}\boldsymbol{v})-{\boldsymbol{\nabla}}\times\boldsymbol{\omega}={-}{\boldsymbol{\nabla}}\times\boldsymbol{\omega}, \end{equation}

the Navier–Stokes equation (2.3) can be rewritten as

(5.2)\begin{equation} \frac{\partial \boldsymbol{v}}{\partial t} +{\boldsymbol{\nabla}}\left(p_r+{Ro}\,\frac{\boldsymbol{v}^2}{2}\right) +({Ro}\,\boldsymbol{\omega}+2\boldsymbol{\varOmega}) \times \boldsymbol{v} +{E}\,{\boldsymbol{\nabla}}\times\boldsymbol{\omega} = 2\omega\sin(2\omega t)\,\boldsymbol{\varOmega}_0\times\boldsymbol{x}. \end{equation}

Multiplying (5.2) by $\boldsymbol {v}$, using $\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {v}=0$ and the identity

(5.3)\begin{equation} \boldsymbol{v}\boldsymbol{\cdot}({\boldsymbol{\nabla}}\times\boldsymbol{\omega})= \boldsymbol{\omega}^2+{\boldsymbol{\nabla}}\boldsymbol{\cdot}(\boldsymbol{\omega}\times\boldsymbol{v}), \end{equation}

yields the energy equation (Wu et al. Reference Wu, Welfert and Lopez2022b)

(5.4)\begin{align} \frac{\partial}{\partial t}\left(\frac{\boldsymbol{v}^2}{2}\right) +{E} \boldsymbol{\omega}^2 - 2\omega\sin(2\omega t)\, \boldsymbol{v}\boldsymbol{\cdot}(\boldsymbol{\varOmega}_0\times\boldsymbol{x}) ={-}{\boldsymbol{\nabla}}\boldsymbol{\cdot}\left[ \left(p_r +{Ro}\,\frac{\boldsymbol{v}^2}{2}\right)\boldsymbol{v} +E (\boldsymbol{\omega}\times\boldsymbol{v})\right]. \end{align}

Integrating (5.4) over the volume of the cuboid $\mathcal {V}$ leads to

(5.5)\begin{equation} {{\rm d} \mathcal{K}}/{{\rm d} t} + {E}\mathcal{E} + \mathcal{W} = \mathcal{D}(\boldsymbol{v}), \end{equation}

where

(5.6)\begin{equation} \mathcal{K} =\int_{\mathcal{V}} \frac{\boldsymbol{v}^2}{2}\,{\rm d}\boldsymbol{x} \end{equation}

is the kinetic energy,

(5.7)\begin{equation} \mathcal{E} =\int_{\mathcal{V}} \boldsymbol{\omega}^2 \, {\rm d} \boldsymbol{x} \end{equation}

is the enstrophy,

(5.8)\begin{equation} \mathcal{W}={-}2\omega\sin(2\omega t) \int_{\mathcal{V}}\boldsymbol{v}\boldsymbol{\cdot}(\boldsymbol{\varOmega}_0 \times\boldsymbol{x})\,{\rm d} \boldsymbol{x} \end{equation}

is the work done by the librational (Euler) force, with $\boldsymbol {v}\boldsymbol {\cdot }(\boldsymbol {\varOmega }_0\times \boldsymbol {x})= \boldsymbol {\varOmega }_0\boldsymbol {\cdot }(\boldsymbol {x}\times \boldsymbol {v})$ the (relative) axial angular momentum, and

(5.9)\begin{equation} \mathcal{D}(\boldsymbol{v})={-}\int_{\mathcal{V}}{\boldsymbol{\nabla}}\boldsymbol{\cdot}\left[ \left(p_r +{Ro}\,\frac{\boldsymbol{v}^2}{2}\right)\boldsymbol{v} +E (\boldsymbol{\omega}\times\boldsymbol{v})\right]\, {\rm d} \boldsymbol{x} \end{equation}

is a residual term. When $\boldsymbol {v}$ has sufficient regularity, Gauss's divergence theorem together with the no-slip boundary conditions results in $\mathcal {D}(\boldsymbol {v})=0$. This is the case for numerical simulations with sufficient spatial resolution obtained at small but non-zero ${E}$.

Figure 16 plots the three terms on the left-hand side of (5.5) over one libration period $\tau ={\rm \pi} /\omega$ for ${A{\kern-4pt}R} =2$ and ${E}={Ro}\in [10^{-8}, 10^{-5}]$, with ${\rm d} \mathcal {K}/ {\rm d} t$ determined from the time series of $\mathcal {K}$ via Fourier spectral differentiation. The figure shows that ${{\rm d} \mathcal {K}}/{{\rm d} t}$ and $\mathcal {W}$ almost balance each other. The dissipation term ${E}\mathcal {E}$ is much smaller but appears to remain bounded away from zero with decreasing ${E}$. This property is known as the zeroth law or the dissipative anomaly in the context of turbulent flows (Onsager Reference Onsager1949; Dubrulle Reference Dubrulle2019; Drivas Reference Drivas2022). The same property holds here as a result of large gradients that develop in the direction transverse to the increasingly thinner interior viscous shear layers that ultimately lead to the low-regularity properties of the attractor as ${E}\to 0$. The finite size of ${E}\mathcal {E}$ is consistent with the scaling $\overline {\boldsymbol {\omega }^2}\sim {E}^{-4/3}$ obtained in figure 5 along the interior viscous shear layer, over a shear layer thickness of order ${E}^{1/3}$ at small non-zero ${E}$. Note that in the regions where the beams forming the shear layer reflect on the walls, they are more intense, with a scaling $\overline {\boldsymbol {\omega }^2}\sim {E}^{-5/3}$ (see figure 6). This occurs in localised regions with an extent scaling with ${E}^{(1/2+1/2)}$, corresponding to the combination of viscous boundary layer depth (${E}^{1/2}$) and the along-wall region of the reflections (${E}^{1/2}$), resulting in a vanishing contribution to ${E}\mathcal {E}$ as ${E}\to 0$.

Figure 16.

Variations over one libration period of the three terms on the left-hand side of (5.5), ${{\rm d} \mathcal {K}}/{{\rm d} t}$, ${E}\mathcal {E}$ and $\mathcal {W}(\boldsymbol {v})$, and their sum, for ${A{\kern-4pt}R} =2$ and ${E}={Ro}$ as indicated.

The sum of the terms on the left-hand side of (5.5), shown in figure 16, although at least four orders of magnitude smaller than both ${\rm d} \mathcal {K}/ {\rm d} t$ and $\mathcal {W}$, and an order of magnitude smaller than the mean dissipation ${E}\mathcal {E}$, does not seem to vanish as ${E}\to 0$. In this limit, (5.5) has to be understood in a distributional sense, and $\mathcal {D}(\boldsymbol {v})$ may not vanish, resulting in anomalous dissipation.

Figure 17 plots the three terms on the left-hand side of (5.5) and their sum for the ${E}= 10^{-8}$ case computed with a range of spatial resolution (using Legendre polynomials of degree ranging from $M=500$ to $M=650$). The results are virtually indistinguishable, indicating that the numerical solutions are well converged, and the finite dissipation and anomaly observed in figure 16 are not numerical artefacts.

Figure 17.

Variations over one libration period of the three terms on the left-hand side of (5.5), ${{\rm d} \mathcal {K}}/{{\rm d} t}$, ${E}\mathcal {E}$ and $\mathcal {W}(\boldsymbol {v})$, and their sum, for ${A{\kern-4pt}R} =2$ and ${E}={Ro}= 10^{-8}$, and spatial resolution $M$ as indicated.

6. Summary and conclusions

A detailed numerical investigation of the structure of inertial wave attractors in rectilinear cuboids is conducted at very small Ekman and Rossby numbers in order to assess which aspects are captured by the linear inviscid theory, implemented in the vertex and edge beam analysis (VEBA) presented in Welfert et al. (Reference Welfert, Lopez and Wu2023), and which are not. VEBA captures most of the geometric features of the intricate network of shear layers emanating from edges and vertices, and their ultimate focusing, at the forcing frequency considered, onto an attractor. This attractor has a quasi-invariant structure in the $x$ direction as a result of the container walls at $x=\pm 0.5$ being parallel to the rotation axis. VEBA also predicts correctly the blow-up of enstrophy density at the attractor due to the overall geometric focusing of wave beams at wall reflections.

The simulations of the full incompressible Navier–Stokes system enable a quantification of this blow-up as the Ekman and Rossby numbers are decreased, with a resulting global viscous dissipation that does not vanish, but instead approaches a static, positive mean value. The development of a non-trivial mean flow in attractor shear layers has also been reported in forced rotating spherical shells (Tilgner Reference Tilgner2007). Our simulations also show that vortex stretching and tilting persist within the attractor, as well as where beams of energy interact with each other, with decreasing Ekman and Rossby numbers, so that the flows are fundamentally three-dimensional and nonlinear. Even for extremely small forcing amplitudes, nonlinear effects have important dynamic consequences for the shear layers forming on attractors. Such a dichotomy between predictions of attractors based on linear inviscid theory and their intrinsically nonlinear and dissipative nature has long been recognised (Maas Reference Maas2005). The configuration studied here, in which the flow remains a perfectly symmetric synchronous limit cycle, not only offers a simple set-up for studying the roles of dissipative and nonlinear effects on inertial waves, but is also a promising candidate for assessing numerically the existence of anomalous dissipation and/or work in the limits of vanishing Ekman and Rossby numbers.