1. Introduction
The interiors of density-stratified or rotating fluids are anisotropic media (Maas Reference Maas2022). These allow for stable hydrostatic and cyclostrophic equilibria, respectively. In the former equilibrium, the fluid density increases monotonically in the direction of gravity. In the latter equilibrium, the angular momentum increases radially outwards from the rotation axis. The single direction in which gravity or the rotation axis points (vertical, say) is distinct from the two orthogonal horizontal directions. In a gravitational field, it is energetically costly to displace fluid parcels vertically. In rotating fluids, this applies to radial displacements. The disturbance is restored by gravity or centrifugal/centripetal forces, as well as by a pressure gradient force imposed by neighbouring fluid.
In the fluid interior, anisotropy manifests itself in disturbances that propagate as internal waves under a fixed angle relative to gravity or the rotation axis (Görtler Reference Thorpe1943; Oser Reference Thorpe1958). The wave’s frequency, restricted by an upper threshold imposed by stratification and/or rotation (Gerkema & Zimmerman Reference Thorpe2008), determines this inclination. Upon reflection from a solid boundary, these waves preserve their frequency and hence their inclination, thus obeying a non-specular type of reflection in the vertical (Bretherton Reference Bretherton1964). This implies that, in an enclosed fluid domain possessing sloping boundaries, reflecting waves exhibit wave focusing and defocusing. In a two-dimensional enclosure, focusing dominates over defocusing. Generically (i.e. for nearly any wave frequency below the threshold) internal waves will approach a limit cycle (Bretherton Reference Bretherton1964; Stewartson Reference Thorpe1971, Reference Thorpe1972) – a wave attractor (Maas & Lam Reference Thorpe1995). In ideal fluids, this attractor represents a spatial singularity, where velocities and pressure gradients increase without bound. This evokes a nonlinear adjustment in which a new pair of waves is generated via a triadic resonant instability (Dauxois et al. Reference Thorpe2017). In real, viscous fluids, this process is accompanied or superseded by the development of free (internal) viscous boundary layers surrounding the wave attractor (Hazewinkel et al. Reference Thorpe2008; Rieutord & Valdettaro Reference Thorpe2010; Le Dizes Reference Thorpe2015; He et al. Reference Thorpe2022).
When the shape of the two-dimensional fluid domain has residual symmetries (such as, e.g. a trapezoidal basin), for certain distinct wave frequencies regular wave modes exist (Maas & Lam Reference Thorpe1995). In such ‘globally resonant modes’, wave focusing is exactly balanced by subsequent defocusing. These modes resemble the classical, complete orthogonal set of eigenmodes, encountered in elliptic problems where the spatial structure is governed by a Helmholtz equation (Riley, Hobson & Bence Reference Thorpe2006). It is therefore tempting to search for such ‘eigenmodes’ in density-stratified and rotating three-dimensional fluids, regardless of the shape of the fluid domain. This search is fuelled by the discovery of such eigenmodes in rotating spheres and then ellipsoids (Bryan Reference Bryan1889; Zhang et al. Reference Zhang, Earnshaw, Liao and Busse2001; Colin de Verdìère & Vidal Reference Thorpe2025). It has been proven that in the full sphere the inertial modes form a complete basis (Ivers, Jackson & Winch Reference Thorpe2015). The completeness of the set of inertial modes in a cylinder rotating along its axis has been discussed in Friedlander & Siegmann (Reference Thorpe1982).
Still, could there exist wave modes in the three-dimensional fluid domain that avoid trapping onto wave attractors? Due to limitations of the separation-of-variables method employed to construct globally resonant wave solutions, recourse is here taken to exploring the short-wave-ray dynamics, comparable to a geometric optics approach. This assumes that internal waves are short compared with basin scales. In the fluid interior, these waves follow straight lines along the energy propagation direction imposed by the group velocity. Upon meeting a boundary, obliquely incident waves reflect and instantaneously refract (Phillips Reference Thorpe1963), i.e. they change their horizontal propagation direction. This is because, in an ideal fluid, the combined fluid particle motion along incident and reflected rays needs to vanish in a direction normal to the boundary. Hence, decomposing the fluid motions into components in the vertical plane normal to the boundary and a component parallel to the boundary, the latter is unchanged during the reflection process, while the former alters due to (de)focusing reflections. Upon iteration of this process, the internal wave ray bounces around the fluid domain. Many rays then reach a limit cycle, a wave attractor. In three-dimensional domains, because of basin symmetries (such as in a paraboloid), this can be a two-dimensional attracting manifold (Maas et al. Reference Thorpe1997; Maas Reference Thorpe2005). However, when this symmetry is absent, the waves can be further focused onto a one-dimensional manifold, a ‘super-attractor’ (Pillet, Maas & Dauxois Reference Thorpe2019; Favier & Le Dizès Reference Thorpe2024).
Interestingly, there exist rays that upon iteration avoid trapping and that resemble whispering gallery modes (WGMs), initially also termed ‘edge waves’ (Manders & Maas Reference Thorpe2004; Drijfhout & Maas Reference Thorpe2007; Rabitti & Maas Reference Pillet, Maas and Dauxois2013, Reference Rabitti and Maas2014; Pillet et al. Reference Thorpe2019). These modes were first discovered in acoustics (Rayleigh Reference Thorpe1878), where they represent sound waves that are guided along an outer wall. These spread much less than sound emitted into the interior of an acoustic chamber due to the lack of geometric spreading when propagating on a surface of lower dimension (the gallery’s wall counter to its interior). The WGMs are also found in other media, such as electromagnetic waves (Mie Reference Thorpe1908; Debye Reference Thorpe1909), matter waves for neutrons (Nesvizhevsky et al. Reference Thorpe2010) and more. For internal waves, some WGMs are instead guided by the critical line, the line on which the boundary slope is equal to the internal wave-ray inclination, that we call localised WGMs. We distinguish these from non-localised WGMs that occur in basins of trapezoidal shape, for which the internal waves first need to traverse the basin back-and-forth upon reflecting from a wall at the opposite side. The mechanism enabling attenuation-less propagation in internal wave WGMs differs from that in acoustic ones. Here, we study WGMs in several basin shapes. This is done by geometric construction of periodic orbits and by examining deviations from said periodic orbits. By perturbing these orbits, we find cases that represent isolated periodic orbits, as well as cases having some meta-stability, that might be used to collectively construct a spatially extended beam of internal waves. Moreover, in the case of localised WGMs, these deviations lead us to discover a new kind of attractor not yet predicted, which lies parallel to the channel’s direction instead of perpendicular to it. We use this new attractor to explain energy accumulation near critical slopes in submarine settings. In Appendix A, we discuss the case of WGMs appearing in enclosed basins, relevant for lakes.
2. Three-dimensional internal wave-ray reflection law
Our system consists of a three-dimensional (3-D) inviscid Boussinesq fluid of linear stratification, with a squared Brunt–Väisälä frequency
$N^2 = -({g}/{\rho _0})({\text{d} \bar {\rho }(z)}/{\text{d}z})$
, where
$\rho _0$
and
$\bar {\rho }$
denote the spatio-temporally uniform and static
$z$
-dependent parts of the density field. Alongside the incompressibility condition,
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0$
, the fluid can be described by momentum and mass conservation equations, respectively,
\begin{equation} \rho _0\frac {{\textrm{D}}\boldsymbol{u}}{{\textrm{D}}t} = -\boldsymbol{\nabla }\!p - \rho {g} \hat {z},\,\,\,\,\,\,\frac {{\textrm{D}}\rho }{{\textrm{D}}t} =0, \end{equation}
where
${\textrm{D}}/{\textrm{D}}t$
denotes the total derivative. Here,
$g$
denotes the acceleration due to gravity, while the non-uniform part of the density,
$\rho=\bar {\rho }+\rho ^{\prime}$
, is split into a static and a dynamic part,
$\bar {\rho }$
and
$\rho ^{\prime}$
, respectively. In Appendix B, we discuss the case in which the fluid is also rotating. We assume that all fields are periodic in time and scale as
$e^{-i\omega t}$
, where
$\omega$
denotes the frequency. After linearising to the first order with respect to a small
$\boldsymbol{u}$
and performing some algebraic manipulations, using subscripts to denote derivatives, we arrive at the linearised Poincaré equation for the spatial perturbation pressure (Maas Reference Thorpe2005)
\begin{equation} p_{\textit{xx}}+p_{\textit{yy}}-\gamma ^2(\omega )p_{\textit{zz}}=0. \end{equation}
Its corresponding dispersion relation reads
\begin{equation} \gamma ^2(\omega ) \equiv \frac {\omega ^2}{N^2-\omega ^2}=\frac {k_x^2+k_y^2}{k_z^2}=\tan ^2{\beta }, \end{equation}
in terms of wave vector
\begin{equation} \boldsymbol{k} = (k_x,k_y,k_z) =\sqrt {k_x^2+k_y^2}\big(\cos {\phi _k},\sin {\phi _k},\pm \gamma ^{-1}\big), \end{equation}
where
$\beta$
denotes the wave vector’s angle from the vertical and the
$\pm$
sign depends on whether the wave vector is pointing up or down, and N is defined above. Ray inclination with respect to the horizontal,
$\gamma$
, depends on the frequency of the wave,
$\omega$
. By scaling horizontal coordinates
$x^{\prime}=x/L,y^{\prime}=y/L$
by horizontal domain scale
$L$
and by stretching the vertical coordinate,
$z^{\prime}= z/(\gamma L)$
, rays of any frequency
$\omega \lt N$
follow characteristics with inclination 1. When examining internal wave rays in a closed basin, this scaling and stretching come at the price of changing its characteristic height
$H$
to the stretched height
$\tau =H/(\gamma L)$
. A full derivation of equations of motion and stretched height from the Boussinesq approximation can be found in Maas & Lam (Reference Thorpe1995). Here,
$\phi _k$
is the horizontal angle of the wave relative to the positive
$x$
direction. The vertical directions of group velocity and phase velocity are opposite to each other. A wave with a wave vector component that points upward will propagate energy with a downward component, and vice versa. When hitting a locally linear wall with slope
$s = \tan {\mu }$
, where
$\mu$
is the wall's slope with the horizontal, a ray characterised by wave vector
$\boldsymbol{k}$
reflects to a different direction
$\boldsymbol{m}$
. In terms of a scaled wall slope
\begin{equation} \eta \equiv \frac {s}{\gamma } \end{equation}
the reflection law determining the horizontal exit angle
$\phi _{\textit{out}}$
can be presented in 3 different, equivalent ways
\begin{align} \tan {\phi _{\textit{out}}} &= \frac {(1 - \eta ^2)\sin {\phi _{\textit{in}}}}{(1 + \eta ^2)\cos {\phi _{\textit{in}}} + 2\textit{sign}(k_z)\eta }, \end{align}
\begin{align} \sin {\phi _{\textit{out}}} &= \textit{sign}(1-\eta )\frac {(1-\eta ^2)\sin {\phi _{\textit{in}}}}{(1 + \eta ^2) + 2\textit{sign}(k_z)\eta \cos {\phi _{\textit{in}}}}, \end{align}
\begin{align} \cos {\phi _{\textit{out}}} &= \textit{sign}(1-\eta )\frac {(1+\eta ^2)\cos {\phi _{\textit{in}}} + 2\textit{sign}(k_z)\eta }{2\textit{sign}(k_z)\eta \cos {\phi _{\textit{in}}} +(1+\eta ^2)}. \end{align}
The bottom slope is called subcritical when the bottom slope is less than the ray slope,
$\eta \lt 1$
, and supercritical if it is steeper,
$\eta \gt 1$
. This is illustrated in panels (a) and (b) of figure 1, where
$\phi _k = \phi _{\textit{in}}$
and
$\phi _m = \phi _{\textit{out}}$
represent the horizontal angles relative to the up-slope
$x$
-direction of the rays that are incident and refracted/reflected from a sloping wall, respectively. Note the reciprocity of these angles when reversing their propagation direction, indicated by arrows. When working in stretched coordinates, in which
$\gamma ^2=1$
, the scaled wall slope,
$\eta$
, becomes equal to
$s$
.
Figure 1. (a) Sub-critical internal wave-ray reflection. (b) Super-critical internal wave-ray reflection. Red arrows indicate wave energy propagation direction. Adapted from Thorpe (Reference Thorpe1997). The angles between arrow tails and the down-slope direction are equivalent to the angles between arrow heads and the up-slope direction, which are used in the reflection law. In both figures
$\beta$
is the ray’s inclination angle and
$\mu$
the wall’s angle with the horizontal plane.
3. Internal wave rays in parabolic channels
In this paper, we look at internal wave rays in infinite channels. The channels we are interested in are symmetric along the along-channel direction, and have a piecewise continuous transverse bottom profile and a flat surface at which we employ the rigid lid approximation. In these channels of stretched depth
$\tau$
, we examine rays of inclination
$\gamma =1$
. From here on, without loss of generality, all channels lie along the
$\hat {y}$
direction.
As the system is homogeneous in the
$\hat {y}$
direction, it is natural to look at the 2-D projection on the
$x{-}z$
plane (see figure 2). In that plane, we can define the ray’s effective inclination
\begin{equation} \gamma^{\prime} \equiv \frac {-1}{\cos {\phi }}. \end{equation}
We note that a ray propagating parallel to the
$y$
-direction (i.e.
$\phi =\pm \pi /2$
) has an effective inclination
$\gamma^{\prime} =\infty$
, and a ray propagating parallel to the
$x$
-direction (i.e. perpendicular to the sloping bottom,
$\phi =0 \mbox{ or } \pi$
) has
$\gamma^{\prime} =\pm 1$
.
In a sense, this projection transforms the 3-D internal wave-ray billiard system to a 2-D billiard with a complex reflection law. A periodic orbit of the 2-D billiard in this
$x{-}z$
projection plane, where rays constituting the periodic orbit have values of
$\gamma^{\prime} \ne 1$
, proves the existence of an internal wave ray in three dimensions that forever propagates along the channel without getting trapped onto a wave attractor: an internal wave WGM; see figure 2.
3.1. Whispering gallery mode existence
When examining cross-sections of a parabolic channel,
$z=\tau x^2$
, it is natural to use the parametrisation used in Maas & Lam (Reference Thorpe1995). The parabolic cross-section is characterised by two parameters, the half-width of the channel, with a scaling equal to 1, and its stretched depth
$\tau$
. As in parabolic channels with
$\tau \lt 0.5$
, all reflections are subcritical, meaning that all rays get attracted to the beaches, we focus on parabolic channels with
$\tau \gt 0.5$
. This also ensures the existence of a critical slope in the channel: the line where the local bottom slope equals the ray slope, which we scaled to 1.
Figure 2. (a) The WGMs in the right half of a parabolic channel, obtained by launching rays from the red dots into the
$y$
-direction. These rays reflect (and, at the bottom, refract) from its intersections with the parabolic boundary and rigid surface. The critical line is shown in magenta. (b) Two-dimensional projection of (a) onto the
$x{-}z$
plane. Vertical lines correspond to rays propagating parallel to the channel.
Although different WGMs may also be constructed, inspired by the ones proposed by Drijfhout & Maas (Reference Thorpe2007), we examine rays launched down channel from the rigid lid surface, i.e. parallel to the
$y$
-direction. A ray launched from
$x_+\gt x_{crit}\equiv ({1}/{2\tau })$
hits a wall with slope
$s=2\tau x_+$
, which focuses it to the angle
$\varphi$
. By plugging
$\phi _{\textit{in}}={\pi }/{2}$
into the reflection law, and remembering that for a ray propagating downwards
$\textit{sign}(k_z)=+1$
,
$\varphi$
and
$x_+$
are obtained from
\begin{equation} \cos {\varphi }=\frac {-2s}{1+s^2}=\frac {-4\tau x_+}{1+(2\tau x_+)^2}\Longleftrightarrow x_+=x_{\textit{crit}}\left(\frac {-\cos {\varphi }}{1+\sin {\varphi }}\right)\!. \end{equation}
After propagating in the 2-D projection billiard with effective inclination
$\gamma^{\prime} \equiv {-1}/{\cos {\varphi }}$
, the ray intersects the channel at
\begin{equation} x_-=x_{\textit{crit}}\left(\frac {-\cos {\varphi }}{1-\sin {\varphi }}\right)=\frac {x_{\textit{crit}}}{s}. \end{equation}
This makes sense as
$s\gt 1$
. At that point, the slope
$2\tau x_-$
reflects the ray to be parallel to the channel due to the reflection law being symmetric; a ray backtracking will perfectly retrace its path. After the ray hits
$(x_-,\tau x_-^2)$
and a subsequent surface reflection, it reverses its own path back to
$x_+$
and parallel to the
$y$
-direction. Therefore, we have found a continuum of periodic orbits in the 2-D billiard, meaning we have found a continuum of WGMs in parabolic channels. These WGMs can be seen in figure 2, both in the 3-D system and in its 2-D projection.
These WGM pairs of
\begin{equation} x_{\pm }=x_{\textit{crit}}\left(\frac {-\cos {\varphi }}{1\pm \sin {\varphi }}\right) \end{equation}
are capped by the width of the channel
$x_+\lt 1$
that in turn dictates
$x_-\gt ({1}/{4\tau ^2})$
. Rays launched outside this bound get trapped onto wave attractors. As these WGMs are centred around the critical line, they are localised WGMs.
We define the WGM step size as the distance a ray travels down the channel during a single period in the 2-D projection billiard. As can be seen in figure 2, all such WGMs share a step size
\begin{equation} \varLambda =4\tau \left(1 - \frac {1}{4\tau ^{2}}\right)\!. \end{equation}
As
$\varLambda$
is independent of launching location,
$x_\pm$
, it remains consistent up until the boundary separating WGMs and cross-isobath attractors.
3.2. Stability analysis of ray trajectories
3.2.1. Spatial perturbations
As we have established the existence of launching conditions leading to periodic orbits, we shall now examine the response to small perturbations from these launching conditions. As the system is invariant under translations in the
$y$
-direction, the WGMs are neutrally stable with respect to perturbations of the launching point’s
$y$
-coordinate. As a ray launched in the
$\hat {y}$
direction from
$(x_0,y_0,\tau - \delta )$
is identical to a ray launched from
$(x_0,y_0- \delta ,\tau )$
, WGMs are also neutrally stable with respect to perturbations of the launching point’s
$z$
coordinate. This indicates the existence of a band of WGMs for all values of
$x$
between
${1}/({4\tau ^2})$
and
$1$
. A ray launched in the
$\hat {y}$
direction outside of this band will get trapped onto a wave attractor parallel to the
$x-z$
plane. This proves that WGMs are also neutrally stable with respect to perturbations of the launching point’s
$x$
coordinate.
If, instead of perturbing the launching parameters, the perturbation is in the wave frequency, one may see it as a perturbation of the basin itself. Specifically, this will augment its effective height
$\tau$
, after the re-stretching of the vertical coordinate to achieve
$\gamma =1$
. In that case, the ray’s trajectory is still a WGM; however, its step size has changed due to the change in
$\tau$
.
3.2.2. Angular perturbations
We now go on to numerically examine rays with a perturbed launching angle
$\varphi = {\pi }/{2} + \epsilon$
with respect to the
$\hat {x}$
direction.
Figure 3. (a) Top view of internal wave beam in a parabolic channel of scaled depth
$\tau =0.7$
launched inwards with respect to the WGM with
$\epsilon \gt 0$
. Red points indicate launching locations, and numbered black points represent successive reflections from the channel’s surface. The magenta critical line is dashed and horizontal angle
$\varphi$
with respect to the depth gradient at surface reflection 1 is denoted in green. (b) Side view of two rays focusing after being launched inwards with respect to the WGM. In both panels, particle velocities increase, as implied by the focusing of the rays. Colour qualitatively highlights diminishing distance between rays in a beam.
We observe that a ray launched ‘outwards’ at the surface, from
$x\gt x_{\textit{crit}}$
, and with
$\epsilon \lt 0$
, upon reflection from the bottom, will refract towards the rim. Upon subsequent reflections, the ray will increase moving into the cross-slope direction until it gets trapped by an across-isobath wave attractor. This is unsurprising as wave attractors in channels are well documented (Manders & Maas Reference Thorpe2004; Drijfhout & Maas Reference Thorpe2007; Pillet et al. Reference Thorpe2019; Bratspiess, Maas & Heifetz Reference Bratspiess, Maas and Heifetz2025). We now observe a ray launched ‘inwards’ from the same location, i.e. having
$\epsilon \gt 0$
. This time the ray will drift inwards, as depicted in figure 3, in the direction of the critical slope. As long as
$\epsilon \lt \epsilon _{\textit{crit}}$
, determined below, the ray will be trapped at and above the critical slope, on a 2-D ‘critical-slope wave attractor’ that, unlike those known up until now, lies parallel to the channel’s direction instead of perpendicular to it. A ray launched with
$\epsilon \gt \epsilon _{\textit{crit}}$
will focus onto an across-isobath wave attractor.
As long as a ray launched from the surface intersects the channel’s bottom on the same side of the critical line as its launching point, it is in the critical-slope attractor’s basin of attraction. In order to find
$\epsilon _{\textit{crit}}$
for a ray launched from
$(x_0,1)$
we calculate the critical effective 2-D inclination for which the ray hits exactly the critical point
$({1}/{2\tau },{1}/{4\tau })$
\begin{equation} \gamma^{\prime} _{\textit{crit}}=\frac {4\tau -1}{4\tau x_0-2}\Rightarrow \frac {\pi }{2}+\epsilon _{\textit{crit}}=\arccos {\left(-\frac {4\tau x_0 - 2}{4\tau - 1}\right)}. \end{equation}
We conclude that the border between trapping at the critical-slope attractor or onto a classical cross-channel attractor is given by rays launched exactly down channel, which constitute the WGMs.
When a beam of parallel rays is launched inwards and drifts towards the critical location, the rays move closer to each other and become denser. From conservation of momentum, it follows that, upon focusing of the internal wave beam, its amplitude increases , as shown by Dauxois & Young (Reference Thorpe1999), and its energy diverges, as happens in the vicinity of the well-known cross-isobath attractors. From this, we can expect increased velocities and energy accumulation around the critical slopes.
At the critical slope itself, the ray dynamics breaks down. The asymptotic analysis in Bratspiess et al. (Reference Bratspiess, Maas and Heifetz2025) predicts that at the critical level, a ray coming from above will be focused perfectly along the cross-channel
$x$
-direction, towards the centre of the channel. However, this trajectory is unavailable due to the channel being convex, rendering the trajectories no longer physically possible. Once the rays approach the vicinity of the critical slope, velocities diverge, and the local length scale goes to zero, indicating that the original scaling analysis no longer holds and nonlinear and dissipative effects are no longer negligible. These new effects will smooth out the consequences of the dynamics’ breakdown, similar to the asymptotic theory in He et al. (Reference Thorpe2022).
The velocity profile of a ray is not uniform along the WGM cycle. From ray densities in figure 3, we can deduce that particle velocities are lowest in the outer vertical section of the orbit, where rays are most sparse, and highest in the diagonal section, where rays are the most dense. In the same manner, the focusing and velocity amplification are non-monotonic during the whole focusing depicted in figure 3, but only when comparing the same sections in different periods. Wave group and phase velocities go to zero as rays approach the critical slope, as can be deduced from taking the limits
\begin{equation} \lim _{|\boldsymbol{k}|\to \infty } \boldsymbol{v}_\phi =\frac {\omega }{|\boldsymbol{k}^2|}\boldsymbol{k}, \qquad \boldsymbol{v}_g=\boldsymbol{\nabla} _k\:\omega, \end{equation}
using the given dispersion relation.
4. Internal wave rays in trapezoidal channels
4.1. Whispering gallery mode existence
The cross-section of a trapezoidal channel is defined by three parameters, stretched height
$\tau$
, base length
$a$
and wall slope
$s = \tan {\mu }$
. A related trapezoidal geometry, but with a sloping bottom instead of a sloping wall, was considered in Pillet et al. (Reference Thorpe2019). By scaling coordinates, one can choose
$\tau =1$
, without loss of generality. This rescaling leaves the two other parameters as sole definers of channel geometry. To avoid rays getting trapped in corners, we focus on channels with
$s\gt 1$
.
Similar to the WGMs in the previous subsection, we once again examine rays launched down channel from the rigid lid surface. Rays launched from above the flat base are uninteresting, so we focus on those launched above the sloping wall. As vertical and horizontal walls reflect internal wave rays specularly, we use the corners in which they intersect as mirrors to achieve the desired periodic orbits.
A ray launched above the sloping wall parallel to the
$y$
-direction will be reflected at a horizontal angle related to the effective inclination
\begin{equation} \gamma^{\prime} =\frac {1+s^2}{2s}. \end{equation}
To derive the geometrical condition restricting the existence of WGMs we examine its potential central ray, a ray launched from a corner with the 2-D effective inclination from (3.6). A WGM exists if this ray hits the sloping wall from below, effectively backtracking its path, as described in the previous paragraph. A simple geometric calculation in Appendix C. yields the condition
\begin{equation} a\frac {s^2+1}{2s} - \left \lfloor { a\frac {s^2+1}{2s}}\right \rfloor \lt \frac {s^2-1}{2s^2}, \end{equation}
that, when satisfied, indicates that a WGM indeed exists. Here, the floor of a given real number
$q$
,
$\lfloor q \rfloor$
, indicates the greatest integer less than or equal to
$q$
. As these WGMs are not centred around the critical line they are non-localised WGMs. As in the previous section, using geometry one can prove that all WGMs in a trapezoidal channel share a step size
\begin{equation} \varLambda = 1+M\frac {s^2+1}{s^2-1}-as\frac {2}{s^2-1}, \end{equation}
Figure 4. (a) The WGMs in a trapezoidal channel, obtained by launching rays down channel from the red dots into the
$y$
direction. All WGMs share
$M=1$
. (b) Two-dimensional projection of the first half of the channel depicted in (a) onto the
$x{-}z$
plane. Vertical lines correspond to rays propagating parallel to the channel.
where
$M$
is the number of flat bottom and surface reflections before the ray reflects specularly from the corner. Figure 4 shows an example, see also Manders & Maas (Reference Thorpe2004). Pillet et al. (Reference Thorpe2019) presented a brief analysis dealing with channels with a different kind of trapezoidal profile, having
$a=0$
, and a subcritical bottom intersecting a vertical wall. In their analysis, they constructed the horizontal angle at which a subcritically reflecting ray undergoes a specular reflection from the sloping bottom by solving
\begin{equation} \begin{aligned} \phi _{\textit{in}}& =\pi -\phi _{\textit{out}}\Rightarrow \cos {\phi _{\textit{in}}+\cos {\phi _{\textit{out}}}}=0\Rightarrow \cos {\phi _{\textit{in}}} \\ & \quad +\, \frac {(1+\eta ^2)\cos {\phi _{\textit{in}}} + 2\eta }{2\eta \cos {\phi _{\textit{in}}} +(1+\eta ^2)}=0\Rightarrow \cos {\phi _{\textit{in}}}=-\eta, \end{aligned} \end{equation}
where the second transition was derived by substituting the reflection law and remembering that for a ray travelling downwards
$\textit{sign}({k_z})=1$
. Here, we extend their analysis by addressing stability with respect to spatial perturbations to launching parameters and by observing that all WGMs in a channel share a step size.
4.2. Stability analysis of ray trajectories
Similar to the behaviour seen in a parabolic channel, WGMs in a trapezoidal channel are neutrally stable with respect to spatial perturbations and unstable with respect to angular perturbations at launch.
Unlike parabolic channels, where horizontal angle instability was demonstrated by numerical simulations, in the trapezoidal channel, instability can be understood analytically. A ray launched in the channel having a single inclined wall has only one slope,
$s$
, to focus and defocus from. This means that reflections from above and below are reciprocal in focusing and defocusing terms. In that case, a ray launched at a small angle
$\epsilon$
with respect to the along-channel
$\hat {y}$
direction preserves its original horizontal angle after a full WGM cycle. This means that, even with a small angle, the ray will consistently drift away from the periodic WGM, eventually becoming trapped at a wave attractor in a plane orthogonal to the sloping wall.
This means that, unlike parabolic channels where WGMs are at the border between two basins of attraction of two different types of wave attractors, in trapezoid channels, where only a single basin of attraction exists, WGMs can be seen as its boundary, having infinite focusing time.
4.3. Whispering gallery beams
Since we can find continuous groups of rays that propagate together while remaining coherent, both in parabolic as well as trapezoidal channels, in both cases, we can define whispering gallery beams (WGBs). Similar to a beam of light, these beams can transport energy and information coherently across long distances. The difference between standard beams and WGBs is that the second type’s existence depends on the existence of a critical slope or corner acting as a critical slope.
5. Relating observations with ray tracing simulations
Our analysis of internal wave WGMs lets us intuitively understand several phenomena in submarine canyons and open ocean basins, specifically along-canyon tidal energy fluxes, similar to ones measured in Aslam, Hall & Dye (Reference Aslam, Hall and Dye2018), and energy and turbulence accumulation, similar to what was observed in Horn & Meincke (Reference Thorpe1976) and Nash et al. (Reference Thorpe2004).
5.1. Internal tide energy fluxes in submarine canyons
Similar to the way in which electromagnetic waves propagate through optical fibres (Senior & Jamro Reference Thorpe2009) and whale sounds propagate across the ocean through the SOund Fixing And Ranging channel (Payne & Webb Reference Thorpe1971; Braun et al. Reference Braun, Arostegui, Thorrold, Papastamatiou, Gaube, Fontes and Afonso2022), so can internal waves propagate along submarine canyons via WGMs, albeit with different mechanisms enabling their propagation. Although our theoretical analysis deals with idealised parabolic channels, the dynamics of real ocean currents in submarine canyons can be approximated near their critical levels. There, channel slopes match ray inclinations, and the channel can locally be seen as straight in along-channel and parabolic in cross-channel directions.
As WGMs enable attenuation-less propagation, one expects that even localised forcing will lead to an energy distribution all along the channel dispersed by said WGMs. Numerical simulations by Drijfhout & Maas (Reference Thorpe2007) and ocean measurements by Aslam et al. (Reference Aslam, Hall and Dye2018) further strengthen this hypothesis.
Figure 5. Adaptations of figure 6 (left) and figure 7 (right) of Aslam et al. (Reference Aslam, Hall and Dye2018). (a) Depth-integrated along-canyon baroclinic M2 (semi-diurnal) energy flux (blue), horizontal kinetic energy (HKE, green), and available potential energy (APE, red) with distance along the thalweg, which is the curve of the lowest elevation in the channel. (b) Along-canyon and (c) across-canyon baroclinic M2 energy flux with distance along the thalweg. Positive along-canyon values are toward the head of the canyon limb. Positive across-canyon values are to the left when looking up-canyon. (d) Along-thalweg slope criticality to the M2 internal tide (blue) and smoothed, using a 5 km running mean (black). Near-critical values (
$0.8 \lt \eta \lt 1.3$
; McPhee-Shaw & Kunze (Reference Thorpe2002)) are indicated in grey. The dashed grey lines indicate the criticality range over which bottom intensification occurs.
Aslam et al. (Reference Aslam, Hall and Dye2018) measured internal tides in the Bay of Biscay canyons adjacent to the Irish Sea and investigated their along- and across-channel energy fluxes. They observed that along-canyon energy fluxes coincide with channel wall slopes being near criticality with respect to the semidiurnal tidal frequency, as is now predicted, which could not be explained by across-isobath wave attractors. Given that our proposed model demands critical slopes in order to transport energy along the canyon, we expect to see no energy flux when the critical slopes are absent, as can be seen in figure 5.
An alternative explanation for the observed energy fluxes may be rectified tides that induce a stable one-way current (Huthnance Reference Thorpe1973). In the context of rectified barotropic tides, we do not find this alternative explanation convincing, as it is known that rectification by barotropic currents is strong only in shallow shelf seas where the barotropic tide is strong (Maas & Zimmerman Reference Thorpe1989). This leaves only the option of rectified internal tides, at locations where the internal tidal beams reflect. As this paper deals with ray tracing within the linearised regime only, we leave a proper treatment of the nonlinear phenomenon of rectified internal tides for future work.
5.2. Energy accumulation near critical slopes
The new mechanism of focusing internal waves towards the critical-slope attractor rather than towards the well-known cross-slope wave attractors can explain the accumulation of energy, turbulence and nonlinear interactions near critical slopes. The diverging ray density near the critical point seen in figure 3, indicating diverging velocities and pressure, is reminiscent of the accumulations observed in nature (Horn & Meincke Reference Thorpe1976; Nash et al. Reference Thorpe2004) and in numerical simulations (Drijfhout & Maas Reference Thorpe2007). Even though critical-slope energy accumulation may be explainable by across-isobath wave attractors, we find it less likely in the settings of continental slopes adjacent to oceans, such as in Horn & Meincke (Reference Thorpe1976), as the lack of horizontal closure allows only along-isobath critical-slope wave attractors.
An alternative explanation for critical-slope intensifications is single reflection focusing. In this proposed scenario, energy accumulates after a single focusing downward-directed reflection from the supercritical side of a slope, which guides the ray towards the critical line. Once focused near the critical line, the wave breaks without defocusing or approaching the attractor, leaving the energy accumulated where it is observed. Although this alternative explanation might be valid, we believe it presupposes an unlikely coincidence, namely that a ray, intersecting a sloping curved boundary, follows a path that precisely hits this boundary at the point where its slope is identical to that of the ray. While it is true that internal tidal beams possess a finite width in cross-beam direction, these widths are likely still very small compared with the depth of the ocean and the radius of curvature of the curved boundary, so that the previous counter-argument still holds.
6. Conclusions
To conclude, via direct geometric construction, we have analytically proved the conjectured existence of internal wave WGMs, i.e. periodic orbits in the 3-D internal wave-ray billiard that do not get absorbed by any wave attractor, in parabolic and trapezoidal channels. Utilising similar techniques, we have investigated a rotationally invariant truncated cone basin (Appendix A).
In channels, we have shown that WGMs are neutrally stable with respect to perturbations in launching location. From this, we deduce the existence of WGBs, continuous groups of rays propagating coherently through the channels. These WGBs can be used to transport energy and information, similarly to electromagnetic or acoustic waves in different waveguides.
By examining perturbations in launching angles in localised WGMs, we have discovered a new kind of wave attractor in parabolic channels. The new wave attractor discovered lies along the critical line and above it, counter to the channel attractors known until now, which lie perpendicular to the channel’s direction. It is worth noting that although the channel basins analysed are invariant along one of their axes, this phenomenon is purely three-dimensional and cannot be reproduced in 2-D basins, which are cross-sections of the full 3-D channels.
We suggest that these WGBs and critical-slope attractors may explain several phenomena observed in nature. Specifically, as a mechanism explaining accumulations and fluxes of energy and nonlinearities in the vicinity of critical slopes, WGBs may account for along-channel energy fluxes observed by Aslam et al. (Reference Aslam, Hall and Dye2018) and critical-slope intensifications.
Regarding the question of whether we can relate the ray approach to a modal approach (as in quantum chaos studies by means of the Weyl–Wigner correspondence) we need to remark that each singular attractor field is described by a double infinite set of modes. Moreover, in a 2-D study of a uniformly stratified fluid contained in a tilted square, an initial regular disturbance projects onto a continuous band of wave attractor frequencies of a type associated with the disturbance’s initial structure (Bajars, Frank & Maas Reference Bajars, Frank and Maas2013). In this context one could expect internal wave systems with localised WGMs to have a discrete branch between two continuum sets of singular eigenmodes, one for each type of wave attractor, and a divergence in the density of states within these continua.
To further connect analytical results from the ray tracing perspective, a future research direction might be a ray tracing investigation performed in natural basins. One may numerically launch rays in the Bay of Biscay submarine canyons using the bathymetry from Aslam et al. (Reference Aslam, Hall and Dye2018). If the proposed connection stands true we will expect to see rays traversing and accumulating, similarly to the kinetic and potential energies measured.
A different way to connect the theoretical work with natural observations will be to perform Computational Fluid Dynamics simulations in both simplified and natural geometries, If our theories hold we should expect to be able to reproduce the phenomena in both linear simulations, equivalent to the ray tracing approximation, and in nonlinear simulations.
One of the other future directions the research may take is answering the unresolved question regarding the fate of non-converging rays in a paraboloidal basin introduced by Maas (Reference Thorpe2005). A potential solution to this question may be that some rays approach a circular critical-slope attractor, bound by circular WGMs, similar to those in the parabolic channel. However, further research is required.
Another potential future direction is a revision of previous works in which 3-D basins are approximated by 2-D ones. Re-examinations may reveal phenomena unexpectedly hidden by the simplification, even though intuitively the 2-D and 3-D systems seem equivalent. A specific basin that may be worthwhile re-investigating is the spherical shell, in which an examination of meridional cross-sections and across-isobath wave attractors might miss WGMs and critical-slope attractors hidden by the 2-D projection that does not preserve the full dynamics of the system, unlike the projections presented in this work.
Acknowledgements
N. Bratspiess wishes to thank Ilias Sibgatullin for insightful discussions, which helped re-imagine the analysis. N. Bratspiess also wishes to thank Laurette Tuckerman and Jacob Sonnenschein for helpful comments given during the judgement of this work in the context of an MSc thesis. We thank the reviewers for perceptive comments on a previous version of the paper.
Declaration of interests
The authors report no conflicts of interest.
Appendix A
Now that we have concluded the analysis of straight channels, basins invariant with respect to a single Cartesian coordinate, we move on to basins invariant with respect to the angular coordinate. Rotationally invariant basins are common in the literature on internal waves (Maas Reference Thorpe2005; Rabitti & Maas Reference Pillet, Maas and Dauxois2013, Reference Rabitti and Maas2014; Sibgatullin et al. Reference Thorpe2017; Pacary et al. Reference Thorpe2023), and are, among other reasons, relevant in the context of stellar dynamics in astrophysics (Rieutord & Valdettaro Reference Thorpe2010; Le Bars et al. Reference Thorpe2010; Ogilvie Reference Thorpe2013) and lakes in the context of geophysical fluid dynamics (Boegman et al. Reference Boegman, Imberger, Ivey and Antenucci2003; Palshin et al. Reference Thorpe2018), such as Lake Geneva (Reiss, Lemmin & Barry Reference Thorpe2023).
A.1. Two-dimensional projection
In a cylindrical coordinate system,
$(\rho ,\theta ,z)$
, bodies of rotation such as cones, spheres, paraboloids and others can be represented as a single variable function
$\rho (z)$
or
$z(\rho )$
. Intuitively, this persuades us to look for the effective 2-D billiard of the 3-D systems, an example of which is depicted in figure 6. We shall do so by examining a differential step of a ray propagating at inclination
$\pm 1$
. In the horizontal plane, a straight path satisfies
\begin{equation} \rho =\frac {\alpha }{\cos {(\theta -\theta _0)}}, \end{equation}
where
$\alpha$
and
$\theta _0$
are determined by the ray’s launching location and direction. The horizontal differential step satisfies
\begin{equation} \text{d}\boldsymbol{l} =(\text{d}\rho ,\rho \text{d}\theta ) \Rightarrow \gamma ^{- 1} \text{d}z=|\text{d}\boldsymbol{l}|=\sqrt {1+\rho ^2\Big(\frac {\text{d}\theta }{\text{d}\rho }\Big)^2}\text{d}\rho =\frac {\rho }{\sqrt {\rho ^2-\alpha ^2}}\text{d}\rho. \end{equation}
From this we calculate the curve which the ray follows in
$(\rho ,z)$
-space
Figure 6. (a) Side view of a wave attractor (yellow) and a non-localised WGM (blue) in a frustum basin. Surface reflections are indicated by black dots. (b) Top view of (a).
\begin{equation} \rho ^2=\alpha ^2+\left ({\gamma ^{- 1}(z-z_0)+\sqrt {\rho _0^2-\alpha ^2}}\right )^2. \end{equation}
The intersection of this hyperbola with the boundary is the next collision point, and the new angle
$\theta$
can be calculated from (4.1). Although this projection is much less elegant due to the projected characteristics being hyperbolas and not straight lines, it is still a simplification of sorts which may be found useful in the future.
A.2. Beyond the 2-D projection
Due to the awkwardness of the 2-D projection explored in the previous subsection, it is sometimes more convenient to approach the 3-D billiard directly. An example for one such case is the frustum, a truncated cone which can be defined by floor radius
$R_{\textit{in}}$
, ceiling radius
$R_{\textit{out}}$
and stretched height
$\tau$
. The wall’s slope is once again
$s = {\tau }/({R_{\textit{out}}-R_{\textit{in}}})$
. In the frustum, one can find periodic or quasi-periodic WGMs by imposing a condition on the sum of incoming and outgoing angles, and then finding the reflection radius on the frustum that results in the desired trajectory. One such WGM is depicted in figure 6. As the frustum basin presented in figure 6 lacks a critical line and the presented WGM passes through the bulk of the basin, one can deduce that it is a non-localised WGM.
Following the same approach, in a paraboloidal basin, results in a bowtie WGM similar to the ones in Rabitti & Maas (Reference Thorpe2014). In rotationally invariant basins that do feature a critical line, like the paraboloidal basin, we expect to see along-isobath critical-slope attractors, however, confirming this statement requires further investigation.
This approach can also be used in the case of basins with discrete symmetries, such as the truncated elliptic cone in Favier & Le Dizès (Reference Favier and Le Dizès2024) or the 3-D stadium in Bratspiess et al. (Reference Bratspiess, Maas and Heifetz2025). In these cases, no simple projection to two dimensions can exist, but the existing symmetries may still enable us to find WGMs.
Appendix B
In the case in which the fluid is subject to rotation at a rate
$\Omega$
, having its axis antiparallel to the direction of gravity an additional term, representing the Coriolis force, must be added to the Euler equation
\begin{equation} \rho _0\left (\frac {{\textrm{D}}\boldsymbol{u}}{{\textrm{D}}t} + \boldsymbol{f}\times \boldsymbol{u} \right ) = -\boldsymbol{\nabla }\!p - g \rho \hat {z},\,\,\,\,\frac {{\textrm{D}}\rho }{{\textrm{D}}t}=0. \end{equation}
Here,
$\boldsymbol{f}=f\hat {z}$
, with
$f= 2 \Omega$
the Coriolis parameter. As before, we assume that all fields are periodic in time and scale as
$e^{-i\omega t}$
. We once again linearise to the first order with respect to a small
$\boldsymbol{u}$
and perform some algebraic manipulations, using subscripts to denote derivatives, and arrive at the linearised Poincaré equation for the spatial perturbation pressure (Maas Reference Thorpe2005)
\begin{equation} p_{\textit{xx}}+p_{\textit{yy}}-\gamma ^2(\omega )p_{\textit{zz}}=0. \end{equation}
However, now that rotation is taken into account, the dispersion relation reads
\begin{equation} \gamma ^2(\omega ) = \frac {\omega ^2-f^2}{N^2-\omega ^2}=\frac {k_x^2+k_y^2}{k_z^2}=\tan ^2{\beta } \end{equation}
in terms of wave vector
\begin{equation} \boldsymbol{k} = (k_x,k_y,k_z) =\sqrt {k_x^2+k_y^2}\big(\cos {\phi _k},\sin {\phi _k},\pm \gamma ^{-1}\big) \end{equation}
and where
$\beta$
once again denotes the wave vector’s angle from the vertical.
As the dispersion relation is of the same form as in cases where rotation is negligible, all analyses presented in the main body of the paper are still valid, with the only changes being in the allowed frequencies of internal waves, which are now also limited by the Coriolis parameter
\begin{equation} f\lt \omega \lt N \qquad \text{or} \qquad N\lt \omega \lt f, \end{equation}
depending whether
$N\gt f$
or vice versa. The group velocity and phase velocity are still orthogonal to each other. When
$N\gt f$
, a wave with a wave vector component that points upward will propagate energy with a downward component, and vice versa. For
$N\lt f$
, phase and energy propagate in the same vertical direction but opposite horizontal direction, meaning that the reflection law must be adjusted accordingly, specifically
$\textit{sign}(k_z) \rightarrow -\textit{sign}(k_z)$
.
It is also worth mentioning that although the fluid rotates either clockwise or counterclockwise in the horizontal plane, there is no asymmetry in horizontal ray propagation between the two options. This can be understood from the fact that the rotation element of the dispersion relation is
$f^2$
, which is invariant under
$f\rightarrow -f$
.
Figure 7. Two-dimensional projection of a trapezoidal billiard with the constructed central WGM (red solid line). As there is a single reflection from the flat bottom before the ray arrives at the corner,
$M=1$
. The dashed red line and bottom trapezoid indicate the geometric relation used for the derivation.
Figure 8. Diagram illustrating different areas in
$(s,a)$
space where WGMs are supported. Respective
$M$
values are annotated in each area. White areas in
$(s,a)$
space represent trapezoidal basins which do not support WGMs.
Appendix C
The corners of a trapezoidal domain act as the reflection points and centres of the WGB. Therefore, the analytical proof of the existence of a WGM in a trapezoidal channel begins at the corners opposite to the sloping wall. We again examine the 2-D projection of a 3-D trapezoidal channel billiard with base length
$a$
, stretched height
$\tau =1$
and wall slope
$s$
. In this projection plane, a ray launched from a corner at an inclination
$\gamma^{\prime} = ({1+s^2}/{2s})$
, that hits the sloping wall, reflects with the new inclination
$\gamma^{\prime \prime}=+\infty$
, which is equivalent to down-channel propagation in the full 3-D system. The vertical ray hits the rigid lid from which it reflects specularly back to the same point on the sloping wall, and from there back to the corner from which it was launched, closing a periodic orbit.
By defining the horizontal distance from the base’s end to the reflection point on the sloping wall
$\delta$
, as illustrated in figure 7, one can derive the connection
\begin{equation} (a+\delta )\gamma^{\prime} =M+\delta s\Rightarrow (a+\delta )(1+s^2)=2Ms+2\delta s^2\Rightarrow \delta =a\frac {s^2+1}{s^2-1}-M\frac {2s}{s^2-1}, \end{equation}
where
$M$
is a non-negative integer representing the number of times the ray reflects specularly from the flat bottom and the rigid lid before hitting the corner. By imposing the geometric constraint
$0\lt \delta \lt {1}/{s}$
, one can derive the condition trapezoidal channels must satisfy in order to support WGMs
\begin{equation} a\frac {s^2+1}{2s}\gt M\gt a\frac {s^2+1}{2s}-\frac {s^2-1}{2s^2}, \end{equation}
which can be shown to be equivalent to
\begin{equation} a\frac {s^2+1}{2s} - \left \lfloor { a\frac {s^2+1}{2s}}\right \rfloor \lt \frac {s^2-1}{2s^2}. \end{equation}
Given a pair of parameters
$(a,s)$
, satisfying the condition above, one can calculate
$M$
using
\begin{equation} M = \left \lfloor { a\frac {s^2+1}{2s}}\right \rfloor \!. \end{equation}
Figure 8 shows areas of parameter space where trapezoidal channels support WGMs and their corresponding values of
$M$
.
Open access
