# Inertial wave super-attractor
in a truncated elliptical cone

Benjamin Favier1 and Stéphane Le Dizès1<EMAIL_ADDRESS>1 Aix
Marseille Univ, CNRS, Centrale Marseille, IRPHE, Marseille, France

###### Abstract

We consider inertial waves propagating in a fluid contained in a non-
axisymmetric three-dimensional rotating cavity. We focus on the particular
case of a fluid enclosed inside a truncated cone or frustum, which is the
volume that lies between two horizontal parallel planes cutting an upright
cone. While this geometry has been studied in the past, we generalise it by
breaking its axisymmetry and consider the case of a truncated elliptical cone
for which the horizontal sections are elliptical instead of circular. The
problem is first tackled using ray tracing where local wave packets are
geometrically propagated and reflected within the closed volume without
attenuation. We complement these results with a local asymptotic analysis and
numerical simulations of the original linear viscous problem. We show that the
attractors, well-known in two dimensional or axisymmetric domains, can be
trapped in a particular plane in three-dimension provided that the axisymmetry
of the domain is broken. Contrary to previous examples of attractors in three-
dimensional domains, all rays converge towards the same limit cycle regardless
of initial conditions, and it is localised in the bulk of the fluid.

##  1 Introduction

Inertial waves and internal gravity waves are waves propagating in rotating
and stratified fluids, respectively. In three dimensions, when forced locally
in a uniformly rotating (or stratified) fluid, these waves have the
particularity to propagate along a double cone that makes a constant angle
with respect to the axis of rotation (or to the direction of stratification).
This property allows to show that, when propagating within a closed container,
wave packets may converge after multiple reflections on solid boundaries to a
particular surface called attractor in 2D (Maas & Lam, 1995) or axisymmetric
geometries such as the spherical shell (Rieutord & Valdettaro, 1997). One of
the simplest geometries giving rise to attractors is a rectangular container
with a sloping boundary. This trapezoidal geometry has been the subject of
numerous works since the first study in a stratified fluid by Maas et al.
(1997). The theoretical attractors obtained by ray tracing have been observed
both experimentally (Hazewinkel et al., 2008, 2010; Scolan et al., 2013) and
numerically (Drijfhout & Maas, 2007; Grisouard et al., 2008). Similar results
were also obtained in rotating fluids in the same geometry (Maas, 2001;
Manders & Maas, 2003, 2004) or in its axisymmetric version (Klein et al.,
2014; Sibgatullin et al., 2019; Boury et al., 2021; Pacary et al., 2023).
Attractors have also been found to be generic features of inertial waves in
spherical shells (Rieutord & Valdettaro, 1997; Rieutord et al., 2000). Early
studies concerned attractors localised close to the equator (Stern, 1963;
Bretherton, 1964; Stewartson, 1972) which possess similar properties as 2D
attractors (Rieutord et al., 2002). But the picture seems more complex in a
spherical shell owing to the presence of the rotation axis in the domain and
of a critical latitude singularity issued from the inner sphere (Rieutord &
Valdettaro, 2010, 2018; He et al., 2022, 2023). The robustness of attractors
has been analysed with respect to wall friction (Beckebanze et al., 2018),
nonlinearity (Grisouard et al., 2008; Jouve & Ogilvie, 2014; Favier et al.,
2014; Beckebanze et al., 2021; Ryazanov et al., 2021) and instabilities
(Brouzet et al., 2016b, 2017; Dauxois et al., 2018). The 2D framework has also
been used to obtain most of the available mathematical results (Maas & Lam,
1995; Manders et al., 2003; Maas, 2009; Bajars et al., 2013; Beckebanze &
Keady, 2016; Colin de Verdière & Saint-Raymond, 2020; Dyatlov et al., 2022;
Makridin et al., 2023).

While fundamental in nature, studies about attractors are motivated by their
ability to focus energy at small length scales, which could potentially impact
the dissipative properties of many geophysical systems. Internal gravity waves
and attractors can break and lead to turbulence (Staquet & Sommeria, 2002;
Brouzet et al., 2016a), and as such are known to play an important role in the
energy budget of the ocean where they are often excited by tides (Wunsch,
1975). Attractors can for example be excited by tidal waves in a paraboloidal
basin (Maas, 2005) or between two ridges (Echeverri et al., 2011). Their
occurrence has been analysed in the configurations of the Mozambique Channel
(Manders et al., 2004) and of the Luzon strait (Tang & Peacock, 2010; Wang et
al., 2015). In the astrophysical context, inertial waves and attractors could
be important in the synchronisation processes of rapidly-rotating
astrophysical objects as they provide a way to rapidly dissipate energy (Zahn,
1975; Ogilvie & Lin, 2004).

Most of the works mentioned above have considered a 2D or a 3D framework with
symmetry (axisymmetric or invariant along one direction). In that case, rays
were assumed to remain confined within a particular 2D plane since the system
is effectively invariant along the third direction. As soon as the rays are
assumed to propagate out of this 2D plane, it becomes important to consider
the 3D reflection law of localised wave beams (Manders & Maas, 2004; Maas,
2005). This was done in an axisymmetric geometry by Maas (2005) and Rabitti &
Maas (2013, 2014), and in a 3D rectangular domain with a sloping wall but
invariant along one direction by Manders & Maas (2004) and Pillet et al.
(2018). For both the axisymmetric spherical shell and the 3D trapezoidal basin
that has a uniform shape in a transverse direction, the authors showed that
the attractor is possibly trapped in a specific plane which depends on the
initial conditions of the ray tracing protocol. To our knowledge, the only
study who considered a completely 3D geometry is the recent work of Pillet et
al. (2019). They considered a trapezoidal geometry, but for which the sloping
plane is also inclined in the transverse direction, thus breaking the symmetry
on which previous studies were implicitly constructed. In that case, the ray
beams tend to drift along the twice-inclined boundary towards the vertical
boundary eventually closing the domain. The rays accumulate there around a
trapped attractor close to the vertical boundary (Pillet et al., 2019), which
is not due to a global focusing but instead due to the arrest by the vertical
side boundary of the continuous shift of the wave beams.

In the present work, focusing on the purely rotating case with inertial waves
only, we consider a 3D non-axisymmetric geometry and show that we obtain a
global focusing, regardless of initial conditions, of the wave beams along a
local curve (as opposed to a surface in the axisymmetric case). Contrary to
the work of Pillet et al. (2019), our attractor does not rely on the
confinement induced by a vertical boundary but is trapped in the bulk of the
fluid domain. We thereby provide the first evidence of a super-attractor, that
is a 1D curve on which all rays tend to focus after multiple reflections on
solid boundaries, regardless of their initial positions and orientations
within the fluid volume. In §2, we present the framework, the 3D reflection
law and our geometry (a truncated elliptic cone). In §3, we demonstrate that,
in our geometry, ray beams corresponding to waves of a given frequency do
converge for almost all initial conditions to a unique limit cycle. This limit
cycle only depends on the geometry and the frequency of the waves. An analytic
formula for the contracting factors obtained in the two focusing directions is
derived in this section and compared to numerical ray tracing results. In §4,
numerical simulations of the linear viscous problem are presented. We show
that the viscous response to a global forcing tends to be focused onto the
super-attractor. A preliminary scaling with respect to the Ekman number is
also proposed for the velocity amplitude, before we finally present our
conclusions.

##  2 Formulation of the problem and methods

In this section, we describe the equations of motion, the reflection law of a
localised wave packet and the particular geometry considered for the rest of
the study.

###  2.1 General problem and equations

We consider the incompressible flow of a fluid of constant kinematic viscosity
$\nu$ contained inside a closed container rotating at a constant rate
$\bm{\Omega}=\Omega\bm{e}_{z}$. While we expect our results to remain valid
when the fluid is stratified in density, due to the similarity in the
dispersion relations and propagation properties, we focus on the purely
rotating case here for simplicity. Using $1/(2\Omega)$ as the time scale and
the characteristic length scale $a$ of the container as the reference length
scale, and focusing on infinitesimal perturbations to the solid body rotation
flow, the linearised dimensionless equations of motion in the rotating frame
are

$\displaystyle\frac{\partial\bm{u}}{\partial t}+\bm{e}_{z}\times\bm{u}$
$\displaystyle=-\nabla P+E\nabla^{2}\bm{u}$ (1)
$\displaystyle\nabla\cdot\bm{u}$ $\displaystyle=0$ (2)

where $\bm{u}$ is the perturbation velocity, $P$ is the pressure incorporating
the centrifugal acceleration and $E=\nu/(2\Omega a^{2})$ is the Ekman number.
We shall focus on the response to an harmonic boundary forcing of
dimensionless frequency $\omega$ and for small Ekman numbers. The frequency of
inertial waves is bounded by twice the rotation rate and we therefore focus on
dimensionless frequencies comprised between 0 and 1.

Information on the solution can be obtained in an inviscid framework by
monitoring the propagation of localised wave beams, as detailed in the next
section.

###  2.2 Reflection law of a localised beam

Our description mainly follows the approach of Manders & Maas (2004), Maas
(2005), Rabitti & Maas (2014) and Pillet et al. (2018, 2019). Yet, we slightly
modify their approach to obtain a simpler reflection law.

As in Rabitti & Maas (2014), we consider a wave beam of frequency $\omega$
that is sufficiently localised such that it travels along a geometrical ray
path. Each ray is characterised by its angle of propagation $\varphi$ with
respect to the vertical rotation axis $\bm{e}_{z}$ and by its azimuthal angle
$\phi$ with respect to the axis $\bm{e}_{x}$. We assume $0\leq\varphi\leq\pi$
and $0\leq\phi<2\pi$. The angle $\varphi$ is given by $\pi/2\pm\theta$ where
the angle $\theta$ (between 0 and $\pi/2$) is fixed by the frequency of the
inertial oscillations according to the dispersion relation of inertial waves,
which in our dimensionless formulation is simply

$\omega=\cos\theta\ .$ (3)

This condition means that the rays propagate in an axial cone that makes an
angle $\theta$ with respect to the horizontal plane. This property is
maintained during the ray propagation including when it reflects on solid
boundaries (Phillips, 1966). In two-dimensional or axisymmetric domains, this
non-specular reflection can lead to limit cycles towards which all rays
converge. These cycles are called attractors.

The horizontal (or azimuthal) angle $\phi$ is irrelevant in two dimensions
since the ray always propagates inside the same predefined plane. In three
dimensions however, the ray is free to move through the whole volume. Contrary
to $\varphi$, $\phi$ is not conserved during reflections on solid boundaries.
The reflection law is not specular but corresponds to a tendency for the ray
to converge towards the vertical plane containing the steepest descent
direction. This has been extensively discussed in Manders & Maas (2004),
Rabitti & Maas (2013, 2014) and Pillet et al. (2018, 2019).

In the present study, we shall use the following relation between the incident
and reflected angles $\phi_{i}$ and $\phi_{r}$ when a ray reflects on a plane
surface inclined by an angle $\alpha$ (between 0 and $\pi/2$) with respect to
the horizontal plane (see figure 1)

$\tan(\phi_{r}-\phi_{\bf
n})=\frac{(\tan^{2}\theta-\tan^{2}\alpha)\sin(\phi_{i}-\phi_{\bf
n})}{(\tan^{2}\theta+\tan^{2}\alpha)\cos(\phi_{i}-\phi_{\bf
n})+2\zeta\tan\alpha\tan\theta}\ .$ (4)

In this expression, $\phi_{\bf n}$ is the azimuthal angle of the normal vector
${\bf n}$ of the surface oriented towards the fluid and $\zeta={\rm
sgn}(n_{z}V_{iz})$, where $n_{z}$ is the vertical component of the normal
vector and $V_{iz}$ is the vertical component of the incident velocity. A
similar expression was derived in Manders & Maas (2004) for wavevector angles
in the context of plane wave reflection. In term of wave propagation
directions, the reflection law has never been written in this form in the
literature. We provide a derivation in Appendix A.

(a) (b)

Figure 1: Reflection of an inertial beam on an inclined surface. The incident
and reflected beams are defined by their velocity vectors ${\bf V_{i}}$ and
${\bf V_{r}}$, the surface by its normal vector ${\bf n}$ oriented toward the
fluid. (a) 3D view; (b) Projected view on the horizontal plane.

###  2.3 Geometry of the fluid domain

We consider the volume contained within a truncated elliptic cone defined by

$x^{2}+\frac{y^{2}}{b^{2}}\leq\frac{z^{2}}{\tan^{2}\alpha}~{},$ (5)

with $\tan\alpha\leq z\leq\tan\alpha+H$. The angle $\alpha$ is the angle made
by the conical surface with respect to the horizontal in the $x$-direction.
The base of the cone always has a unit radius along the $x$-axis while it
reaches $b$ along the $y$-axis. Without loss of generality, we assume $b\geq
1$ so that the short axis is always along the $x$-axis. As it will become
apparent later, the attractor will localise in that case onto the particular
plane $(Oxz)$. The case $b=1$ corresponds to the classical axisymmetric
truncated cone geometry. This type of axisymmetric geometries involving a
truncated conical surface, sometimes called a frustum, has already been
studied previously (Beardsley, 1970; Henderson & Aldridge, 1992; Borcia &
Harlander, 2013; Klein et al., 2014; Sibgatullin et al., 2019; Pacary et al.,
2023). The novelty of our study is to extend the linear wave dynamics to the
case of non-axisymmetric domains which correspond in our case to $b>1$.
Although theoretical results will be derived for any values of $H$ and
$\alpha$, the numerical investigations will focus on the particular parameters
$H=1$ and $\alpha=\pi/4$ and consider super-critical slopes only for which
$\theta<\alpha$.

##  3 Ray tracing and local analysis

In this section, we discuss the properties of the three-dimensional ray paths
satisfying the reflection laws discussed above. We start by discussing the
axisymmetric case ($b=1$) before considering ray paths in the non-axisymmetric
geometry ($b>1$).

###  3.1 Axisymmetric truncated cone

In this section, we consider the case $b=1$ in equation (5) which corresponds
to an axisymmetric truncated cone. This geometry has been recently considered
in Pacary et al. (2023). It corresponds to the axisymmetric version of the 2D
trapezoidal geometry that has been studied in numerous works, as described in
the introduction. Owing to the axisymmetry, the normal vector of the conical
surface is always oriented toward the axis of symmetry $Oz$. This implies that
if a ray is oriented towards this axis, it reaches the conical surface with a
meridional angle $\phi_{i}=\phi_{\bf n}+\pi$, and reflects with an angle
$\phi_{r}=\phi_{\bf n}$ as prescribed by (4). It therefore continues to be
oriented towards the axis. This means that if a ray initially lies within a
meridional plane, it remains confined to this plane forever as in a 2D
geometry.

(a)(b)

Figure 2: (a) 2D ray paths for two particular trajectories starting from the
orange empty circles and for $\omega=0.8$. The red path corresponds to the
last half of the many reflections and shows the final attractor. (b) Position
of each reflection along the $x$ axis of the last few reflections as a
function of $\omega$. The two vertical lines correspond to $\omega=1/\sqrt{2}$
and $\omega=2/\sqrt{5}$.

The ray paths confined within a particular meridional 2D plane of the
axisymmetric cone, which is equivalent to a 2D trapezoid, have been analysed
in many papers (Maas & Lam, 1995; Maas et al., 1997; Grisouard et al., 2008;
Brunet et al., 2019; Pacary et al., 2023). They converge towards attractors
owing to the focusing effect resulting from the reflection along the conical
slope. This is illustrated in figure 2(a) where we show two ray paths
converging towards an attractor for the particular case $\omega=0.8$, $H=1$
and $\tan\alpha=1$. This particular attractor (shown in red in figure 2(a)) is
composed of two symmetrical quadrangles and is the simplest attractor (with
the lowest number of reflections) that can be obtained in this geometry. It
exists for $\alpha>\theta$ (i.e. when the reflection on the sloping boundary
is subcritical) and in a finite range of parameters that can be obtained by
finding the coordinates of the reflection points. For instance, the points
$(\pm x_{a},z_{a})$ where the attractor reflects on the inclined boundary (see
figure 6 below) can be expressed in terms of the angles $\alpha$ and $\theta$
and the height $H$ as

$x_{a}=\frac{H}{\tan\theta}~{},~{}~{}z_{a}=\frac{H\tan\alpha}{\tan\theta}\ .$
(6)

Writing the condition that $z_{a}$ is between $\tan\alpha$ and $\tan\alpha+H$,
implies the following condition

$\tan\theta<H<\frac{\tan\alpha\tan\theta}{\tan\alpha-\tan\theta}\ .$ (7)

This condition can be written in term of the frequency $\omega=\cos\theta$ as

$\frac{1}{\sqrt{1+H^{2}}}<\omega<\frac{\tan\alpha+H}{\sqrt{(\tan\alpha+H)^{2}+H^{2}\tan^{2}\alpha}}$
(8)

which defines the frequency range for which this type of attractor exists. For
the parameters used in figure 2 ($H=1$ and $\tan\alpha=1$), it corresponds to
the interval $1/\sqrt{2}<\omega<2\sqrt{5}$. This interval can be seen in
figure 2(b) where we show the positions on the $x$-axis of each reflection
(after a large number of reflections in order to focus on the attractor path)
for many random initial conditions uniformly distributed across the surface
and as a function of frequency. There are no attractors for frequencies below
$\omega<1/\sqrt{2}$. The relatively empty regions for these low frequencies
correspond to rays trapped in the upper corners. For $\omega>2/\sqrt{5}$,
other attractors are still observed but they are now characterised by a more
complex path involving multiple reflections on each of the boundaries. In the
following, we will focus on the simpler attractor observed for frequencies
satisfying the conditions (8).

If we now authorise the ray beams to deviate from a particular meridional
plane, their path becomes more complex. One has to monitor the horizontal
angle $\phi$ of the direction of the ray and the meridional angle $\psi$ of
the position where it reflects on the boundaries (see figure 6(b) below). This
problem has been recently studied in Pacary et al. (2023) for this particular
geometry (see also section 7 in Maas (2005)). They showed that the 2D
attractor is still obtained but its location along the azimuthal direction is
now dependent on the choice of initial conditions. This phenomenon is
illustrated in figure 3 where we show that two rays emitted from the same
location with different initial horizontal angles of propagation $\phi$ end on
the same 2D attractor but in different meridional planes. However, the
geometry being azimuthally invariant, all the meridional planes are possible
and should be obtained with the same probability. A fully axisymmetric
attractor would then be obtained (with a non-uniform distribution of trapping
planes though, see Maas (2005)) if all the possible initial conditions were
simultaneously considered (with the exception of whispering gallery modes
(Pillet et al., 2019) which we did not observe here). More details can be
found in Pacary (2023); Pacary et al. (2023), where such an axisymmetric
system is explored using both ray tracing and experiments.

(a) (b)

Figure 3: Ray path for two particular trajectories initialised at the same
point indicated by the empty circle. The two trajectories differ only by the
initial horizontal angle of propagation. The thin lines correspond to the
transient propagation while the thick lines correspond to the final limit
cycle after many reflections. The conical surface is coloured in light grey to
distinguish it from the horizontal planes. The height of the cone is $H=1$,
its opening angle is $\tan\alpha=1$ and the frequency is $\omega=0.8$. (a)
Side view. (b) View from the top. A movie (Movie1) showing the propagation of
many randomly initialised rays can be found in Supplementary materials.

###  3.2 Elliptic cone

In this section, we consider the unexplored case $b>1$ in equation (5) which
corresponds to a truncated elliptic cone. As the axisymmetry is now broken, we
do not expect the presence of any axisymmetric attractor. In the following, we
focus on the case where the reflection on the conical surface remains
subcritical which leads to the following upper bound on $b$

$b<b_{c}=\frac{\tan\alpha}{\tan\theta}\ .$ (9)

The first effect of the elliptic deformation is to modify the orientation of
the normal vector of the conical surface. Except in the vertical planes $x=0$
and $y=0$ corresponding to the directions of the principal axes of the
elliptical cone, the normal vector to the cone is no longer oriented towards
the vertical rotation axis. This means that contrary to the axisymmetric case,
no ray can stay trapped in a particular meridional plane apart from the planes
$x=0$ and $y=0$. When a ray is initialised inside any other meridional plane,
it deviates from it after its first reflection on the conical surface. Whether
the planes $x=0$ and $y=0$ are stable or unstable equilibrium will be
discussed below.

(a) (b)
(c) (d)

Figure 4: Same as figure 3 but for a non-axisymmetric domain with $b=1.2$ in
equation (5). The other parameters are $\omega=0.8$, $H=1$ and $\tan\alpha=1$.
(a) Side view in the $(y,z)$ plane. (b) Side view in the $(x,z)$ plane. (c) 3D
view. (d) Top view in the $(x,y)$ plane. A movie (Movie2) showing the
propagation of many randomly initialised rays can be found in Supplementary
materials.

(a)(b)

Figure 5: Positions of the reflection points along the $x$ axis in red and
along the $y$ axis in blue for the final limit cycle obtained from many random
initial conditions. (a) Axisymmetric case $b=1$ and (b) elliptic case $b=1.2$.
The cone height is $H=1$ and its opening angle is $\tan\alpha=1$. The graph
focuses on the frequency range defined by equation (8). The thin lines
correspond to the analytical position of the six reflection points of the 2D
attractor.

We first repeat the ray tracing experiment done previously for $b=1$ but in a
non-axisymmetric geometry with $b=1.2$. The same frequency $\omega=0.8$ is
considered. An example is shown in figures 4, where two ray paths are
generated from the same initial position but with different horizontal angles.
We observe that the two ray paths now converge towards the same unique
attractor localised in the plane $y=0$ containing the semi-minor axis of the
elliptical cone. One can repeat the experiment for many initial positions and
horizontal angles $\phi$ showing that all rays eventually end up in this
particular plane $y=0$ and along the same limit cycle. For this particular
case, it thus seems that the plane $y=0$ is a stable equilibrium while the
plane $x=0$ is unstable. This is further confirmed in figure 5 where we show
the positions $x$ and $y$ of the reflection points of the limit cycle
(obtained by considering the last ten iterations of a total of five hundred)
for many frequencies and many random initial conditions. We compare the
axisymmetric case $b=1$ in figure 5(a) with the elliptic case $b=1.2$ in
figure 5(b). The thin lines indicate the $x$ coordinates of the six reflection
points of the 2D attractor path, as already shown in Figure 2(b), in the
particular plane $y=0$ which is the same for all ellipticities $b$ in our
case. For the axisymmetric case, each independent ray path converges towards a
similar 2D attractor in a particular meridional plane, depending on its
initial conditions, leading to a dense pattern of reflection points when
plotting their $(x,y)$ coordinates. Note that the projection from the
axisymmetric domain to the arbitrary coordinates $x$ and $y$ leads to a non-
uniform density of reflection points which tend to accumulate close to the
expected positions of the 2D attractors (shown by the thin lines). For the
non-axisymmetric case $b=1.2$ however, all the rays converge towards the
particular plane $y=0$, regardless of their initial conditions in terms of
position and initial horizontal angle, while we recover the same attractor
structure as in 2D when considering the $x$ coordinates. Note that we recover
this peculiar property for all frequencies within the attractor range given by
equation (8).

The attractor in the particular plane $y=0$ seems to attract all the rays,
regardless of their initial positions within the volume, at the exception of
rays initialised within the $x=0$ plane which is also an equilibrium, albeit
unstable. Contrary to the axisymmetric case, there is a second azimuthal
convergence (related to the varying curvature of the elliptical cone along the
azimuthal direction) in addition to the meridional convergence at the origin
of the classical 2D attractor (related to the inclination of the conical
surface). For this reason, we call this final limit cycle a “super-attractor”,
to differentiate it from the classical attractor surface observed in the
axisymmetric case. The super-attractor results from a convergence of rays in
both axial and azimuthal directions while there is no azimuthal convergence
for an attractor in an axisymmetric geometry, as its azimuthal position
depends on the initial position and horizontal angle of propagation.

####  3.2.1 Local analysis of the super-attractor

In this section, we analyse the local properties of the super-attractor. Our
objective is to confirm the attracting character of the super-attractor both
towards a particular limit cycle in the meridional plane and towards a
particular meridional plane regardless of initial conditions. While the first
property is shared with regular attractors, we show that the second is
specific to super-attractors. We aim to obtain an analytic expression for the
attraction rate (i.e. Lyapunov coefficient) as a function of the geometrical
parameters.

We focus on the super-attractor which is made of two symmetrical quadrangles
as illustrated in figure 2(a). This 2D attractor has been studied in section
§3.1. It exists in the plane $x=0$ for frequencies in the interval defined in
(8). The points $(\pm x_{a},0,z_{a})$ where reflection occurs has been given
in (6).

(a)(b)

Figure 6: (a) Perspective view of the elliptical cone in thick solid line. The
attractor is indicated in red line. The cone is twice symmetrically duplicated
to avoid considering the reflections on the horizontal planes. (b) View from
the top of a ray emitted close to the attractor. The red dashed lines indicate
the section of the cone at the altitude $z_{a}$ of the attractor. The ray
emitted from the point $(x_{0},y_{0},z_{0})$ with a horizontal angle
$\phi_{0}$ is reflected at $(x_{1},y_{1},z_{1})$ with the angle $\phi_{1}$.

We now want to analyse the behaviour of a ray close to the attractor. We
consider a ray emitted from a point $(x_{0},y_{0},z_{0})$ on the inclined
surface close to $(x_{a},0,z_{a})$. This ray is emitted downward with a
vertical angle $\varphi_{0}=\pi/2+\theta$ and an azimuthal angle $\phi_{0}$ as
illustrated in figure 6. If this ray is not too far from the attractor, it
first reflects on the lower plane surface, then reflects on the upper plane
surface before reaching the inclined surface on the other side in a point
$(x_{1},y_{1},z_{1})$ close to $(-x_{a},0,z_{a})$. Both points
$(x_{0},y_{0},z_{0})$ and $(x_{1},y_{1},z_{1})$ are on the cone defined by (5)
so their horizontal coordinates can be written as & x_0 = z_0 cosψ0tanα , y_0
= b z_0 sinψ0tanα ,
x_1 = z_1 cosψ1tanα , y_1 = b z_1 sinψ1tanα , where $\psi_{0}$ and $\psi_{1}$
are their azimuthal angles (see figure 6(b)).

As the reflections on the plane surfaces do not modify the azimuthal angle of
the ray (see equation (4) for the particular case $\alpha=0$), the coordinates
of the point $(x_{1},y_{1},z_{1})$ can be obtained by continuing the ray in
horizontal mirror images of the cones. Such a ray reaches the boundary of the
second cone image in a point which has just been shifted vertically by twice
the height of the cone, that is in $(x_{1},y_{1},z_{1}-2H)$ (see figure 6).
This property means that there exists $\lambda$ such that &x_1 = x_0 -
λcosθcosϕ_0 ,
y_1 = y_0 - λcosθsinϕ_0 ,
z_1 -2H = z_0 + λsinθ .

When the ray is close to the attractor, $(x_{0},y_{0},z_{0})$ and
$(x_{1},y_{1},z_{1})$ are assumed to be close to $(x_{a},0,z_{a})$ and
$(-x_{a},0,z_{a})$, respectively. Moreover, $\psi_{0}$ is assumed to be small,
and $\phi_{0}$ and $\psi_{0}$ are close to $\pi$. At first order, equations
(6) give for the relative distances & δx_0 ∼δz0tanα , δy_0 ∼b z_a ψ0tanα ,
δx_1 ∼- δz1tanα , δy_1 ∼- b z_a (ψ1\- π) tanα . From (6c), one gets $\lambda$
as

$\lambda\sim\frac{-2H+\delta z_{1}-\delta z_{0}}{\sin\theta}\ .$ (10)

Inserting this expression in (6a), one obtains using (6a,b)

$\delta x_{1}\sim-K\delta x_{0}~{},~{}~{}\delta z_{1}\sim K\delta z_{0}\ ,$
(11)

with

$K=\frac{\tan\alpha-\tan\theta}{\tan\alpha+\tan\theta}=\frac{\sin(\alpha-\theta)}{\sin(\alpha+\theta)}\
,$ (12)

which corresponds to the 2D contraction factor of the attractor.

In (6b), we obtain

$\delta y_{1}\sim\delta y_{0}-\frac{2H(\phi_{0}-\pi)}{\tan\theta},$ (13)

which gives a first relation between the angles

$\psi_{1}-\pi=-\psi_{0}+\frac{2}{b}(\phi_{0}-\pi)~{}.$ (14)

A second relation that expresses the angle $\phi_{1}$ of the reflected ray in
terms of $\phi_{0}$ and $\psi_{0}$ is obtained by applying the condition of
reflection (4) at $(x_{1},y_{1},z_{1})$. Close to the attractor, the azimuthal
angle $\phi_{\bf n}$ of the normal vector ${\bf n}$ is given at leading order
by

$\phi_{\bf n}\sim\frac{(\psi_{1}-\pi)}{b}~{}.$ (15)

The normal vector ${\bf n}$ and the incident rays are oppositely oriented with
respect to the vertical, so $\zeta={\rm sgn}(n_{z}V_{iz})=-1$. Expression (4)
then gives for small angles

$\phi_{1}-\phi_{\bf
n}\sim-\frac{\tan\alpha-\tan\theta}{\tan\alpha+\tan\theta}(\phi_{0}-\pi-\phi_{\bf
n})=-K(\phi_{0}-\pi-\phi_{\bf n}).$ (16)

Using (14) and (15), this expression finally reduces to

$\phi_{1}\sim-\frac{1+K}{b}\psi_{0}+\left(\frac{2(1+K)}{b^{2}}-K\right)(\phi_{0}-\pi)~{}.$
(17)

Expressions (14) and (17) can be written as

$\bm{\Psi}_{1}=\bm{\mathcal{M}}\bm{\Psi}_{0}$ (18)

for the vectors $\bm{\Psi}_{1}=(\psi_{1}-\pi,\phi_{1})^{\top}$ and
$\bm{\Psi}_{0}=(\psi_{0},\phi_{0}-\pi)^{\top}$ where the matrix
$\bm{\mathcal{M}}$ is

$\bm{\mathcal{M}}=\left(\begin{array}[]{cc}-1&2/b\\\\[2.84544pt]
-(1+K)/b&2(1+K)/b^{2}-K\end{array}\right).$ (19)

The operation can be repeated and after $2N$ reflections on the cone boundary,
that is $N$ cycles, we get

$\bm{\Psi}_{2N}=\bm{\mathcal{M}}^{2N}\bm{\Psi}_{0}$ (20)

where $\bm{\Psi}_{2N}=(\psi_{2N},\phi_{2N}-\pi)^{\top}$.

Introducing the eigenvalues $\lambda_{\pm}$ and associated eigenvectors
$\bm{\Psi}_{\pm}$ of the matrix $\bm{\mathcal{M}}$ which are defined by & λ_±
= (1+K)(2/b2-1) ±((1+K)2(2/b2-1)2\- 4K)1/22 ,
Ψ_±= (2/b, 1 + λ_±)^⊤ , $\bm{\Psi}_{2N}$ can be written as

$\bm{\Psi}_{2N}=C_{0}(\lambda_{+})^{2N}\bm{\Psi}_{+}+D_{0}(\lambda_{-})^{2N}\bm{\Psi}_{-}\
,$ (21)

where $C_{0}$ and $D_{0}$ are constants depending on the initial condition
only & C_0= 2(ϕ0-π) - b(1+λ-)ψ02(λ+-λ-) ,
D_0= 2(ϕ0-π) - b(1+λ+)ψ02(λ--λ+) .

The functions $\lambda_{+}$ and $\lambda_{-}$ characterise the behaviour of
the angles $\psi_{2N}$ and $\phi_{2N}-\pi$ as a function of the cycle number
$N$. These angles go to zero if and only if $|\lambda_{\pm}|<1$. This
condition is here equivalent to $b>1$ (since $0<K<1$). The functions
$\lambda_{\pm}$ depend on the value of $b$ with respect to the particular
values

$b_{c\pm}=\frac{\sqrt{2(1+K)}}{1\pm\sqrt{K}}\ .$ (22)

The two $\lambda_{\pm}$ are real positive for $1<b<b_{c+}$, and real negative
for $b>b_{c-}$. For $b_{c+}<b<b_{c-}$, they are complex conjugates with a
constant absolute values equal to $\sqrt{K}$. When $b=b_{c\pm}$, the solution
evolves differently as shown in appendix B.

####  3.2.2 Comparison between local analysis and numerical ray tracing

This section compares the prediction of the local analysis with that of the
global ray tracing approach described above. To do so, we randomise the
initial position inside the volume, the horizontal angle of propagation and
the sign of the vertical velocity component. We then track the ray for
$10^{4}$ reflections on boundaries which is enough to reach the attractor in
most cases. The horizontal angle is computed at each reflection on the conical
surface and for $x_{i}>0$ according to
$\psi_{i}=\textrm{arg}(x_{i}+\textrm{i}y_{i})=\textrm{atan2}(y_{i},x_{i})$
where $(x_{i},y_{i})$ are the coordinates of the $i^{\textrm{th}}$ reflection
point on the positive half $x>0$ of the conical surface. The horizontal angle
difference $\Delta\psi_{i}=|\psi_{i}-\psi_{i-1}|$ is then tracked as a
function of the number of cycles. Note that we track the convergence of the
horizontal angle after a complete cycle around the attractor for comparison
with the prediction (6) of the local analysis presented in section 3.2.1. In
that case, the analysis is performed on the number of cycles around the
attractor and not on the number of reflections on the boundaries (there are 6
reflections on boundaries and 2 reflections on the conical surface per cycle
for the particular attractor considered here).

(a) (b)

Figure 7: (a) Angle difference $\Delta\psi_{N}=|\psi_{N}-\psi_{N-1}|$ in
logarithmic scale as a function of the number of cycles $N$. We show 100
independent realisations with different initial conditions for each value of
$b$. (b) Decay rate $\sigma$ defined by $\Delta\psi_{N}\sim\exp(\sigma N)$ as
a function of the long axis $b$. Symbols correspond to the ray tracing
approach while lines correspond to the linear local analysis (6). The upper
branches correspond to $\sigma=\ln|\lambda_{+}|$ while the lower branches
correspond to $\sigma=\ln|\lambda_{-}|$. Note that the data for $\omega=0.76$
stop because of the upper bound $b_{c}\approx 1.17$ as given by (9).

The evolution of the horizontal angle increment is typically characterised by
a short transient followed by an exponential decay until machine accuracy is
eventually reached and the ray becomes trapped in a particular plane. Some
examples for $\omega=0.8$, $H=1$, $\tan\alpha=1$ and three values of $b$ are
shown in Figure 7(a) to illustrate the convergence process. We observe that
all rays converge towards the super-attractor at the same exponential rate,
irrespective of the initial conditions. The duration of the transient before
the ray effectively converges towards the attractor depends on the initial
conditions and corresponds to the early reflections far from the final limit
cycle. Note that for this particular case, the particular value of $b$ defined
by equation (22) is $b_{c+}=1.097$. Consistent with the theoretical
prediction, we observe real values for the decay rate when $b<b_{c+}$ (see the
two cases $b=1.02$ and $b=1.05$ in figure 7(a)) and complex values when
$b>b_{c+}$ (see the case $b=1.2$).

A best fit $\Delta\psi_{N}\sim\exp(\sigma N)$, where $N$ is the number of
cycles, is computed in the range $10^{-14}<\Delta\psi_{N}<10^{-3}$ and
averaged over $10^{3}$ independent rays initialised randomly within the whole
volume. Results are shown in Figure 7(b). An excellent agreement between the
local theoretical prediction and the ray tracing approach is observed for
various values of $b$ and three different frequencies within the attractor
range. Interestingly, the decay rate is $\sigma=\ln|\lambda_{-}|=\ln|K|$ when
$b=1$ while it is smaller when $b>1$. There is therefore a drastic difference
between the case $b=1$, for which all rays converge rapidly toward a plane
different for each ray, and the case $b>1$ for which the convergence rate is
smaller and actually tends to 0 when $b\rightarrow 1$. Although the breaking
of the axisymmetry does create a globally attracting plane as already
discussed previously, the rate of convergence towards this unique plane
increases with $b$ and is actually maximum once $b>b_{c+}$. In that case, it
is actually half of that observed when $b=1$.

##  4 Numerical simulations of the linear viscous problem

Up to now, we have only discussed the properties of the ray paths which are
only valid in the limit of vanishing Ekman numbers. The link between the
properties of the ray paths and the actual viscous solution of the original
linear Navier-Stokes equations (1) is not obvious. In this context, it is
desirable to check whether the globally-attracting solutions discussed in
previous sections have any counterpart when considering the direct solution of
equations (1).

To that end, we solve the linear viscous equations (1) using the spectral
element solver Nek5000 (Fischer, 1997; Deville et al., 2002). The domain is
discretised using a number $\mathcal{E}$ of hexahedral elements. Elements have
been refined close to all boundaries to properly resolve viscous Ekman layers.
The velocity is discretised within each element using Lagrange polynomial
interpolants based on tensor-product arrays of Gauss-Lobatto-Legendre
quadrature points. The polynomial order $l_{d}$ of the expansion basis on each
element is fixed to $11$ in this study, while the number of elements goes up
to $\mathcal{E}=29952$ for $E=10^{-7}$. Convergence has been tested by
gradually increasing the polynomial order for a fixed number of elements. The
Coriolis term is treated explicitly by a third-order extrapolation scheme
whereas the viscous terms are treated implicitly by a third-order backward
differentiation scheme. Similarly to the ray tracing approach discussed
previously, we focus on the particular case $H=1$, $\alpha=\pi/4$ and
$\omega=0.8$.

(a)(b)

Figure 8: Azimuthal profile of the velocity amplitude averaged over time and
over the 2D attractor path on each meridional plane. (a) Variable long axis
$b$ at constant Ekman number $E=10^{-7}$. The amplitude is normalised by its
azimuthal average. (b) Variable Ekman number at constant long axis $b=1.2$.
The forcing frequency is $\omega=0.8$ and the forcing pattern is defined by
equation (23) in all cases.

Figure 9: Vertical velocity component in the planes (a) $y=0$, (b) $x=0$ and
(c) $z=1.9$. The same color scale is used for the three plots. Parameters are
$E=10^{-7}$, $\omega=0.8$ and $b=1.2$. The inclined line in (a) indicates the
profile of the plane shown in Figure 10. The dotted lines in (a) and (b)
indicate the plane $z=1.9$ shown in (c) while the dash lines show the 2D
attractor path in each meridional plane.

Since we want to compare the axisymmetric conical geometry with its elliptical
counterpart, care must be taken when choosing the forcing. In particular,
latitudinal libration is not well suited since the forcing would be of viscous
nature in the axisymmetric case (since the velocity at the boundaries would be
purely tangential) while there would be a non-zero normal velocity and hence
pressure coupling in the non-axisymmetric elliptical case. Another constraint
comes from the corners which will inevitably contribute to the viscous
dynamics by emitting their own singular shear layers. Taking into account
these considerations, we opted for a rotating vertical forcing at the bottom
plane of the cone, sometimes called negative nutation (Sibgatullin et al.,
2017), defined by

$u_{z}(z=1)=\begin{cases}\big{(}x\cos(\omega t)-y\sin(\omega
t)\big{)}f(r)&\text{if $r<1$}\\\ 0&\text{if $r>1$}\end{cases}\ ,$ (23)

where $f(r)=2r^{3}-3r^{2}+1$ is a smoothing function to ensure that the
forcing vanishes close to the corners and $r=\sqrt{x^{2}+y^{2}}$ is the
cylindrical radius. The other two velocity components are zero and the other
boundaries are all no-slip. A similar forcing has already been used both
experimentally (Pacary et al., 2023) and numerically (Sibgatullin et al.,
2019). Note that we obtained qualitatively similar results when considering an
axisymmetric forcing similar to that used by Boury et al. (2021) for example.

From a fluid at rest, we run the simulation until a periodic response is
obtained. In order to quantify the inhomogeneities between different attractor
planes, we first define the 2D attractor path for each plane obtained from the
intersection between a vertical plane containing the origin and the frustum.
Each plane is parametrised with its azimuthal angle with respect to the $x$
direction. The plane crossing the short elliptical axis thus corresponds to
$\phi=0$ while the plane crossing the long elliptical axis corresponds to
$\phi=\pm\pi/2$. Note that the 2D attractor path on each plane depends on the
ellipticity $b$. While the same attractor path is expected on each individual
plane when $b=1$, different attractors (but the same rectangular topology) are
expected when $b>1$. We then compute the averaged velocity amplitude along
each path by averaging over time once the periodic response is obtained and
along the attractor path. This averaged amplitude is a constant in the
axisymmetric case $b=1$ but depends on the orientation of the plane once
$b\neq 1$. As we vary the ellipticity of the domain, the amplitude of the
response also varies. In order to focus on the azimuthal inhomogeneities
induced by the wave attractor, we further normalise the amplitude by its
average over all azimuthal angles $\phi$.

This ratio is plotted in figure 8(a) for a fixed $E=10^{-7}$ and varying
$b\geq 1$. As expected, it is unity for the axisymmetric case $b=1$. As $b$
increases, a clear focusing of the energy along the short axis corresponding
to $\phi=0$ is observed. Note that a second local maximum is also observed
along the long axis $\phi=\pm\pi/2$. The same simulations are now run at a
fixed $b=1.2$ and varying Ekman number from $E=10^{-4}$ down to $E=10^{-7}$.
Results are displayed in figure 8(b). While focusing is observed for all Ekman
numbers, it is more and more pronounced as the Ekman number decreases. Note
again that a residual localisation also persists around the long axis
$\phi=\pm\pi/2$. This secondary localisation of the energy is not expected
from the ray tracing approach only since it corresponds to an unstable
equilibrium as discussed previously in section 3.2.1. We have observed such
residual localisation strongly depends on the particular choice of forcing and
seems to be partially dependent on contributions from the bottom corner (see
Figure 9(b)), which goes well beyond our current understanding mostly based on
local ray tracing.

(a) (b) (c)
(d) (e) (f)

Figure 10: Maps of averaged velocity amplitude normalised by its maximum value
on the local plane shown in Figure 9(a). The Ekman number is decreasing from
left to right. The top row (a-c) corresponds to the non-axisymmetric case
$b=1.2$ while the bottom row (d-f) corresponds to the axisymmetric case $b=1$.
The line shows the iso-contour corresponding to 90% of the local maximum value
around the theoretical attractor location which is the origin of the local
frame used here.

A visual inspection of the wave field qualitatively confirms this conclusion.
Figure 9 shows the vertical velocity component at a particular time once a
periodic response has been obtained. One can see that the response along the
attractor path is much stronger in the $y=0$ plane (see Figure 9(a)) than it
is in the $x=0$ plane (see Figure 9(b)). Note of course that the 2D attractor
path (shown as dash lines in Figures 9(a) and (b)) is not the same in both
planes since the domain is non-axisymmetric. This explains why we observe
apparently curved wave beams. They correspond to the gradual modification of
the attractor path as the geometry changes along the azimuthal direction due
to the weak elliptical deformation. Similar structures were for example
observed by Bühler & Muller (2007). We also see in Figure 9(b) that the
attractor path in the unstable plane $x=0$ is very close to the bottom corner,
which perhaps explains why we have observed a secondary peak of energy in
Figure 8. The bottom corner forces a local shear layer which might overlap
with the attractor beam in a non-trivial way. Nevertheless, Figure 9(c)
clearly shows that the amplitude of the stable attractor localised around the
plane $y=0$ is significantly larger than that of any other planes, including
the unstable attractor localised in the $x=0$ plane.

In order to show the local structure of the super-attractor, we define a plane
perpendicular to the attractor path along the short axis, as shown in Figure
9(a). The velocity amplitude is averaged over time and normalised by its local
maximum closest to the attractor path to obtain the maps displayed in Figure
10. We compare the axisymmetric case $b=0$ on the bottom row with the non-
axisymmetric case $b=1.2$ on the top row. For the axisymmetric case, one
observes the gradual focusing of the axisymmetric wave beam around the
attractor position as the Ekman number is decreased. The observed curvature is
due to intersection between the flat plane and the curved axisymmetric wave
beam. This axisymmetric pattern is clearly broken for $b=1.2$. At $E=10^{-5}$,
while we observe a local maximum close to the theoretical position of the
super-attractor (which corresponds to the origin with our choice of
coordinates), no local beam is observed. At lower Ekman numbers however, a
localised beam is observed with a complex anisotropic structure. It is more
elongated along the transverse direction $y$ than along the in-plane
coordinate $s$. This confirms that the energy injected by the global large-
scale forcing is eventually focused preferentially onto the super-attractor
path localised around the $y=0$ plane. While it is too early to convincingly
discuss possible scalings with the Ekman number, we nevertheless report the
amplitude scaling observed in our simulations in Figure 11. We consider three
different measures of amplitudes. The first is obtained by averaging over all
2D attractor paths of each meridional section, which we refer to as global.
The second corresponds to the average on the short axis attractor only while
the last corresponds to the local amplitude maximum obtained from the maps
displayed in Figure 10. We observe that for all three measures, the amplitude
increases as the Ekman number decreases. The local measure might follow the
scaling $E^{-1/6}$, which is expected from classical 2D attractors (He et al.,
2023) forced by inviscid forcing, although a much larger range of Ekman
numbers should be explored to convincingly conclude on this matter. This is
left for future works.

Figure 11: Amplitude scaling with Ekman for three different measures: the
average is performed over all attractor paths irrespective of their azimuthal
orientation (Global), only along the ray path corresponding to the short axis
(Short axis) or focusing on the local maximum observed on the maps shown in
Figure 10 (Local). Note that the local value at $E=10^{-4}$ is absent since we
cannot unambiguously measure a local maximum in that case.

##  5 Discussion

We have shown how a 2D attractor can be localised in 3D into a particular
plane by breaking the axisymmetry of the fluid domain. All ray paths have been
observed to converge to a particular 1D curve contained within a specific
plane. A local analysis has allowed us to confirm the exponential decay of the
horizontal angle and has been successfully compared with ray tracing
computations. While numerical simulations have shown a localisation of the
energy onto the attracting plane, the harmonic response also contains
contributions from the corners which might be at the origin of a secondary
peak of energy along the long-axis which is unstable according to our ray
tracing analysis. By construction of our particular geometry, the secondary
unstable attractor is more and more affected by the corner contributions as
$b$ increases, which might explain why we observe a secondary energy increase
there. However, the observed increase of kinetic energy along the stable
super-attractor can only be explained by a secondary focusing since its
topology within the particular plane $y=0$ remains the same regardless of the
elliptical deformation $b$.

While our example has convincingly shown the existence of such super-
attractors, generalising this result to other geometries would be valuable. We
have focused on super-critical reflections only since $\alpha>\theta$ in our
case while it is known that the focusing properties of the 3D reflection is
non-trivially depending on the criticality of the slope (Pillet et al., 2019).
We suspect that similar super-attractors exist in other conical geometries,
including those involving sub-critical slopes, such as those used by Klein et
al. (2014), Sibgatullin et al. (2017) and Boury et al. (2021), provided that
the axisymmetry used up to now is broken. Other geometries known to support 2D
or axisymmetric attractors, such as the spherical shell or the
paraboloidal/parabolic-shaped stadium, can probably generate super-attractors
once elliptically deformed. In that respect, much remains to be done to bridge
the gap between the local analysis close to the attracting plane and the
global properties for which the whole geometry, and not just its asymptotic
behaviour close to the attracting plane, matters.

We have considered the particular case of an elliptic deformation in order to
break the axisymmetry. While it has the advantage of being a well-defined
perturbation of the axisymmetric reference case, it is certainly not the only
way to proceed. It would be interesting to consider other types of non-
axisymmetric geometries like cuboids (Wu et al., 2023). Possible connections
with quantum chaos in stadium billiards resulting in scarred patterns could
also be explored (Kudrolli et al., 2001). One could also further investigate
the required geometrical properties of the fluid container necessary for the
emergence of super-attractors. Finally, it is important to realise that, while
the final super-attractor is indeed one-dimensional in the sense that all rays
converge towards a one-dimensional parametric curve, this particular attractor
is contained inside a plane since it corresponds to the localisation of the
otherwise 2D attractor. An interesting question would be to know whether
general 1D curves, not necessarily contained within the same plane, can be
super-attractors and if so in which geometries. This is a tremendous
geometrical problem and we hope this preliminary work will motivate future
studies in that direction.

## Acknowledgments

The authors would like to acknowledge Leo Maas for fruitful discussions and
relevant comments on an earlier version of the manuscript. Centre de Calcul
Intensif d’Aix-Marseille is acknowledged for granting access to its high-
performance computing resources. This work was granted access to the HPC
resources of IDRIS under the allocation 2023-A0140407543 made by GENCI.

## Declaration of Interests

The authors report no conflict of interest.

##  Appendix A 3D reflection law

In this section, the 3D reflection law is obtained. The relation between the
incident and reflected angles $\phi_{i}$ and $\phi_{r}$ can be obtained by
writing down the conditions on the velocity vectors ${\bf V_{i}}$ and ${\bf
V_{r}}$ during the reflection process. We consider a plane surface inclined by
an angle $\alpha$ (between 0 and $\pi/2$) with respect to the horizontal plane
and defined it using its normal vector ${\bf n}$ oriented towards the fluid.
The first condition is the condition of non-penetrability which reads

$({\bf V_{r}}+{\bf V_{i}})\cdot{\bf n}=0.$ (24)

The second condition is the conservation of the horizontal velocity
perpendicular to the normal which can be written as

$({\bf V_{r}}\times{\bf e_{z}})\cdot{\bf n}=({\bf V_{i}}\times{\bf
e_{z}})\cdot{\bf n}.$ (25)

They are two other conditions that prescribe the direction of propagation of
the incident and reflected beams which should be on the axial cone of angle
$\theta$

$\tan^{2}\theta||{\bf V_{r}}\times{\bf e_{z}}||^{2}=|{\bf V_{r}}\cdot{\bf
e_{z}}|^{2},$ (26) $\tan^{2}\theta||{\bf V_{i}}\times{\bf e_{z}}||^{2}=|{\bf
V_{i}}\cdot{\bf e_{z}}|^{2}.$ (27)

To obtain a general formula, the idea is to express the vectors ${\bf n}$,
$\bf V_{r}$ and $\bf V_{i}$ in the frame $(x^{\prime},y^{\prime},z)$ oriented
along the steepest descent direction (see figure 1).

The above equations can then be written as &(V_rx’ + V_ix’) tanφ_n \+ (V_rz +
V_iz) = 0 ,
V_ry’ = V_iy’ ,
tan^2θ(V_rx’^2 + V_ry’^2) = V_rz^2 ,
tan^2θ(V_ix’^2 + V_iy’^2) = V_iz^2 , where $\varphi_{\bf n}$ is the angle
(between 0 and $\pi$) of the vector ${\bf n}$ with respect to the vertical
axis. Using (Ab) to eliminate $V_{ry^{\prime}}$ and $V_{iy^{\prime}}$ in the
equation obtained by subtracting (Ad) to (Ac) gives

$\tan^{2}\theta(V_{rx^{\prime}}^{2}-V_{ix^{\prime}}^{2})=V_{rz}^{2}-V_{iz}^{2}~{},$
(28)

which, after eliminating $v_{ix^{\prime}}$ using (Aa), reduces to

$(\tan^{2}\theta-\tan^{2}\varphi_{\bf n})V_{rx^{\prime}}=2\tan\varphi_{\bf
n}V_{iz}+(\tan^{2}\theta+\tan^{2}\varphi_{\bf n})V_{ix^{\prime}}~{}.$ (29)

Using

$\frac{V_{iz}}{V_{iy^{\prime}}}=\frac{{\rm
cotan}\varphi_{i}}{\sin(\phi_{i}-\phi_{\bf
n})}~{},~{}~{}\frac{V_{ix^{\prime}}}{V_{iy^{\prime}}}={\rm
cotan}(\phi_{i}-\phi_{\bf
n})~{},~{}~{}\frac{V_{rx^{\prime}}}{V_{ry^{\prime}}}={\rm
cotan}(\phi_{r}-\phi_{\bf n})$ (30)

we get an expression that gives the angle $\phi_{r}-\phi_{\bf n}$:

$\tan(\phi_{r}-\phi_{\bf n})=\frac{(\tan^{2}\theta-\tan^{2}\varphi_{\bf
n})\sin(\phi_{i}-\phi_{\bf n})}{(\tan^{2}\theta+\tan^{2}\varphi_{\bf
n})\cos(\phi_{i}-\phi_{\bf n})+2\tan\varphi_{\bf n}{\rm cotan}\varphi_{i}}.$
(31)

This formula applies to all configurations. Only the term $\tan\varphi_{n}{\rm
cotan}\varphi_{i}$ depends on the orientation of the normal ${\bf n}$ and of
the incident beam ${\bf V_{i}}$ with respect to the vertical. As
$\varphi_{i}=\pi/2-\theta$ and $\varphi_{\bf n}=\alpha$ (resp.
$\varphi_{i}=\pi/2+\theta$ and $\varphi_{\bf n}=\pi-\alpha$) when these
vectors are oriented upwards (resp. downwards), we get

$\tan(\phi_{r}-\phi_{\bf
n})=\frac{(\tan^{2}\theta-\tan^{2}\alpha)\sin(\phi_{i}-\phi_{\bf
n})}{(\tan^{2}\theta+\tan^{2}\alpha)\cos(\phi_{i}-\phi_{\bf n})\pm
2\tan\alpha\tan\theta},$ (32)

where the sign $-$ is taken when the two vectors ${\bf n}$ and ${\bf V_{i}}$
are oriented oppositely with respect to the vertical (one of the two vectors
is oriented upward and the other downward). This formula works for all types
of reflections (subcritical or supercritical), and for horizontal and vertical
surfaces as well. For this reason, we think that it is more convenient than
the formulas that have been used in Rabitti & Maas (2014) and Pillet et al.
(2018).

##  Appendix B Local behavior for $b=b_{c\pm}$

For the values $b=b_{c\pm}$, the two eigenvalues are equal (to
$\lambda_{c\pm}=\pm\sqrt{K}$) and equations (21a,b) giving $C_{0}$ and $D_{0}$
break down. In that case, the two eigenvectors given by (6b) are also equal
(to $\bm{\Psi}_{c\pm}$). One should introduce another vector, say
$\bm{E}_{1}=(1,0)^{\top}$ to decompose the vector $\bm{\Psi}_{2p}$. Using the
fact that when $b=b_{c\pm}$

$\bm{\mathcal{M}}^{2p}\bm{E}_{1}=K^{p}\bm{E}_{1}\mp
pK^{p-1/2}\sqrt{2(1+K)}\bm{\Psi}_{c\pm},$ (33)

one can easily show that

$\bm{\Psi}_{2p}=\left(C_{c\pm}K^{p}\mp
D_{c}\,pK^{p-1/2}\sqrt{2(1+K)}\right)\bm{\Psi}_{c\pm}+D_{c}K^{p}\bm{E}_{1}$
(34)

where

$C_{c\pm}=\frac{\phi_{0}-\pi}{1\pm\sqrt{K}}~{},~{}~{}D_{c}=\psi_{0}-\sqrt{\frac{2}{1+K}}(\phi_{0}-\pi).$
(35)

## References

  * Bajars et al. (2013) Bajars, J., Frank, J. & Maas, L.R.M. 2013 On the appearance of internal wave attractors due to an initial or parametrically excited disturbance. J. Fluid Mech. 714, 283–311.
  * Beardsley (1970) Beardsley, R.C. 1970 An experimental study of inertial waves in a closed cone. Studies in Applied Mathematics 49 (2), 187–196.
  * Beckebanze et al. (2018) Beckebanze, F., Brouzet, C., Sibgatullin, I. N. & Maas, L. R. M. 2018 Damping of quasi-two-dimensional internal wave attractors by rigid-wall friction. J. Fluid Mech. 841, 615–635.
  * Beckebanze et al. (2021) Beckebanze, F., Grayson, K. M., Maas, L. R. M. & Dalziel, S. B. 2021 Experimental evidence of internal wave attractor signatures hidden in large-amplitude multi-frequency wave fields. J. Fluid Mech. 915, A41.
  * Beckebanze & Keady (2016) Beckebanze, F. & Keady, G. 2016 On functional equations leading to exact solutions for standing internal waves. Wave Motion 60, 181–195.
  * Borcia & Harlander (2013) Borcia, I.D. & Harlander, U. 2013 Inertial waves in a rotating annulus with inclined inner cylinder: comparing the spectrum of wave attractor frequency bands and the eigenspectrum in the limit of zero inclination. Theoretical and Computational Fluid Dynamics 27, 397–413.
  * Boury et al. (2021) Boury, S., Sibgatullin, I., Ermanyuk, E., Shmakova, N., Odier, P., Joubaud, S., Maas, L. R. M. & Dauxois, T. 2021 Vortex cluster arising from an axisymmetric inertial wave attractor. J. Fluid Mech. 926, A12.
  * Bretherton (1964) Bretherton, F. P. 1964 Low frequency oscillations trapped near the equator. Tellus 16, 181–185.
  * Brouzet et al. (2017) Brouzet, C., Ermanyuk, E., Joubaud, S., Pillet, G. & Dauxois, T. 2017 Internal wave attractors: different scenarios of instability. J. Fluid Mech. 811, 544–568.
  * Brouzet et al. (2016a) Brouzet, C., Ermanyuk, E. V., Joubaud, S., Sibgatullin, I. & Dauxois, T. 2016a Energy cascade in internal-wave attractors. Europhysics Letters 113 (4), 44001\.
  * Brouzet et al. (2016b) Brouzet, C., Sibgatullin, I. N., Scolan, H., Ermanyuk, E. V. & Dauxois, T. 2016b Internal waves attractors examined using laboratory experiments and 3D simulations. J. Fluid Mech. 793, 109–131.
  * Brunet et al. (2019) Brunet, M., Dauxois, T. & Cortet, P.-P. 2019 Linear and nonlinear regimes of an inertial wave attractor. Phys. Rev. Fluids 4 (3), 034801.
  * Bühler & Muller (2007) Bühler, O. & Muller, C.J. 2007 Instability and focusing of internal tides in the deep ocean. J. Fluid Mech. 588, 1–28.
  * Dauxois et al. (2018) Dauxois, T., Joubaud, S., Odier, P. & Venaille, A. 2018 Instabilities of internal gravity wave beams. Annu. Rev. Fluid Mech. 50, 1–28.
  * Deville et al. (2002) Deville, M. O., Fischer, P. F. & Mund, E. H. 2002 High-Order Methods for Incompressible Fluid Flow. Cambridge University Press.
  * Drijfhout & Maas (2007) Drijfhout, S. & Maas, L. R. M. 2007 Impact of channel geometry and rotation on the trapping of internal tides. J. Phys. Ocean. 37, 2740–63.
  * Dyatlov et al. (2022) Dyatlov, S., Wang, J. & Zworski, M. 2022 Mathematics of internal waves in a 2d aquarium arXiv:2112.10191.
  * Echeverri et al. (2011) Echeverri, P., T., Yokossi, Balmforth, N. J. & Peacock, T. 2011 Tidally generated internal-wave attractors between double ridges. J. Fluid Mech. 669, 354–374.
  * Favier et al. (2014) Favier, B., Barker, A. J., Baruteau, C. & Ogilvie, G. I. 2014 Non-linear evolution of tidally forced inertial waves in rotating fluid bodies. Mon. Not. R. Astron. Soc. 439, 845–860.
  * Fischer (1997) Fischer, P. F. 1997 An overlapping Schwarz method for spectral element solution of the incompressible Navier–Stokes equations. Journal of Computational Physics 133, 84–101.
  * Grisouard et al. (2008) Grisouard, N., Staquet, C. & Pairaud, I. 2008 Numerical simulation of a two-dimensional internal wave atractor. J. Fluid Mech. 614, 1–14.
  * Hazewinkel et al. (2010) Hazewinkel, J., Tsimitri, C., Maas, L. R. M. & Dalziel, S. B. 2010 Observations on the robustness of internal wave attractors to perturbations. Phys. Fluids 22, 107102\.
  * Hazewinkel et al. (2008) Hazewinkel, J., Van Breevoort, P., Dalziel, S. B. & Maas, L. R. M. 2008 Observations of the wavenumber spectrum and evolution of an internal wave attractor. J. Fluid Mech. 598, 373–382.
  * He et al. (2022) He, J., Favier, B., Rieutord, M. & Le Dizès, S. 2022 Internal shear layers in librating spherical shells: the case of periodic characteristic paths. J. Fluid Mech. 939, A3.
  * He et al. (2023) He, J., Favier, B., Rieutord, M. & Le Dizès, S. 2023 Internal shear layers in librating spherical shells: the case of attractors. J. Fluid Mech. (submitted).
  * Henderson & Aldridge (1992) Henderson, G.A. & Aldridge, K.D. 1992 A finite-element method for inertial waves in a frustum. Journal of Fluid Mechanics 234, 317–327.
  * Jouve & Ogilvie (2014) Jouve, L. & Ogilvie, G.I. 2014 Direct numerical simulations of an inertial wave attractor in linear and nonlinear regimes. J. Fluid Mech. 745, 223–250.
  * Klein et al. (2014) Klein, M., Seelig, T., Kurgansky, M., Ghasemi, A. V. & others 2014 Inertial wave excitation and focusing in a liquid bounded by a frustrum and a cylinder. J. Fluid Mech. 751, 255–297.
  * Kudrolli et al. (2001) Kudrolli, A., Abraham, M. C. & Gollub, J. P. 2001 Scarred patterns in surface waves. Physical Review E 63 (2), 026208.
  * Maas (2001) Maas, L.R.M. 2001 Wave focusing and ensuing mean flow due to symmetry breaking in rotating fluids. J. Fluid Mech. 437, 13–28.
  * Maas (2005) Maas, L.R.M. 2005 Wave attractors: linear yet nonlinear. International Journal of Bifurcation and Chaos 15 (09), 2757–2782.
  * Maas (2009) Maas, L. R. M. 2009 Exact analytic self-similar solution of a wave attractor field. Physica D 238, 502–505.
  * Maas et al. (1997) Maas, L. R. M., Benielli, D., Sommeria, J. & Lam, F.-P. A. 1997 Observation of an internal wave attractor in a confined, stably stratified fluid. Nature 388, 557–561.
  * Maas & Lam (1995) Maas, L. R. M. & Lam, F.-P. 1995 Geometric focusing of internal waves. J. Fluid Mech. 300, 1–41.
  * Makridin et al. (2023) Makridin, Z. V., Khe, A. K., Sibgatullin, I. N. & Ermanyuk, E. V. 2023 Forced internal wave attractors: Linear inviscid theory. Phys. Rev. Fluids (in press).
  * Manders et al. (2003) Manders, A.M., Duistermaat, J.J. & Maas, L.R.M. 2003 Wave attractors in a smooth convex enclosed geometry. Physica D: Nonlinear Phenomena 186 (3-4), 109–132.
  * Manders & Maas (2003) Manders, A. M. & Maas, L. R. M. 2003 Observations of intertial waves in a rectangular basin with one sloping boundary. J. Fluid Mech. 493, 59–88.
  * Manders & Maas (2004) Manders, A. M. & Maas, L. R. M. 2004 On the three-dimansonal structure of the inertial wave filed in a rectangular basin with one sloping boundary. Fluid Dyn. Res. 35, 1–21.
  * Manders et al. (2004) Manders, A. M., Maas, L. R. M. & Gerkema, T. 2004 Observation of internal tides in the Mozambique Channel. J. Geophys. Res. 109, C12034.
  * Ogilvie & Lin (2004) Ogilvie, G. I. & Lin, D. N. C. 2004 Tidal dissipation in rotating giant planets. Astroph. J. 610, 477–509.
  * Pacary (2023) Pacary, C. 2023 Dynamique des milieux stratifiés ou tournants: Attracteurs d’ondes axisymétriques, interaction non-linéaire et panaches en rotation. PhD thesis, Ecole Normale Supérieur de Lyon.
  * Pacary et al. (2023) Pacary, C., Dauxois, T., Ermanyuk, E., Metz, P., Moulin, M. & Joubaud, S. 2023 Observation of inertia-gravity wave attractors in an axisymmetric enclosed basin. J. Fluid Mech. (submitted).
  * Phillips (1966) Phillips, O. M. 1966 The Dynamics of the Upper Ocean. Cambridge University Press.
  * Pillet et al. (2018) Pillet, G., Ermanyuk, E.V., Maas, L.R.M., Sibgatullin, I.N. & Dauxois, T. 2018 Internal wave attractors in three-dimensional geometries: trapping by oblique reflection. J. Fluid Mech. 845, 203–225.
  * Pillet et al. (2019) Pillet, G., Maas, L.R.M. & Dauxois, T. 2019 Internal wave attractors in 3d geometries: a dynamical systems approach. Eur. J. Mech. B/Fluids 77, 1–16.
  * Rabitti & Maas (2013) Rabitti, A. & Maas, L.R.M. 2013 Meridional trraping and zonal propagation of inertial waves in a rotating shell. J. Fluid Mech. 729, 445–470.
  * Rabitti & Maas (2014) Rabitti, A. & Maas, L.R.M. 2014 Inertial wave rays in rotating spherical fluid domains. J. Fluid Mech. 758, 621–654.
  * Rieutord et al. (2000) Rieutord, M., Georgeot, B. & Valdettaro, L. 2000 Wave Attractors in Rotating Fluids: A Paradigm for Ill-Posed Cauchy Problems. Phys. Rev. Lett. 85, 4277–4280.
  * Rieutord & Valdettaro (1997) Rieutord, M. & Valdettaro, L. 1997 Inertial waves in a rotating spherical shell. J. Fluid Mech. 341, 77–99.
  * Rieutord & Valdettaro (2010) Rieutord, M. & Valdettaro, L. 2010 Viscous dissipation by tidally forced inertial modes in a rotating spherical shell. J. Fluid Mech. 643, 363–394.
  * Rieutord & Valdettaro (2018) Rieutord, M. & Valdettaro, L. 2018 Axisymmetric inertial modes in a spherical shell at low Ekman numbers. J. Fluid Mech. 844, 597–634.
  * Rieutord et al. (2002) Rieutord, M., Valdettaro, L. & Georgeot, B. 2002 Analysis of singular inertial modes in a spherical shell: the slender toroidal shell model. J. Fluid Mech. 463, 345–360.
  * Ryazanov et al. (2021) Ryazanov, D. A., Providukhina, M. I., Sibgatullin, I. N. & Ermanyuk, E. V. 2021 Biharmonic attractors of internal gravity waves. Fluid Dyn. 56, 403–412.
  * Scolan et al. (2013) Scolan, H., Ermanyuk, E. & Dauxois, T. 2013 Nonlinear fate of internal waves attractors. Phys. Rev. Lett. 110, 234501.
  * Sibgatullin et al. (2017) Sibgatullin, I., Ermanyuk, E., Maas, L.R.M., Xu, X. & Dauxois, T. 2017 Direct numerical simulation of three-dimensional inertial wave attractors. In 2017 Ivannikov ISPRAS Open Conference (ISPRAS), pp. 137–143. IEEE.
  * Sibgatullin et al. (2019) Sibgatullin, I., Xu, X., Tretyakov, A. & Ermanyuk, E. 2019 Influence of geometry on energy flow and instability in inertial wave attractors for rotating annular frustum. In AIP Conference Proceedings, , vol. 2116. AIP Publishing.
  * Staquet & Sommeria (2002) Staquet, C. & Sommeria, J. 2002 Internal gravity waves: from instabilities to turbulence. Annual Review of Fluid Mechanics 34 (1), 559–593.
  * Stern (1963) Stern, M. E. 1963 Trapping of low frequency oscillations in an equatorial “boundary layer”. Tellus 15, 246–250.
  * Stewartson (1972) Stewartson, K. 1972 On trapped oscillations in a slightly viscous rotating fluid. J. Fluid Mech. 54, 749–761.
  * Tang & Peacock (2010) Tang, W. & Peacock, T. 2010 Lagrangian coherent structures and internal wave attractors. Chaos 20, 017508\.
  * Colin de Verdière & Saint-Raymond (2020) Colin de Verdière, Y. & Saint-Raymond, L. 2020 Attractors for two-dimensional waves with homogeneous Hamiltonians of degree 0. Comm. Pure Appl. Math. 73, 421–462.
  * Wang et al. (2015) Wang, G., Zheng, Q., M., Lin., Dai, D. & Qiao, F. 2015 Three dimensional simulation of internal wave attractors in the Luzon Strait. Acta Oceanol. Sin 34, 14–21.
  * Wu et al. (2023) Wu, K., Welfert, B.D. & Lopez, J.M. 2023 Inertial wave attractors in librating cuboids. J. Fluid Mech. 973, A20.
  * Wunsch (1975) Wunsch, C. 1975 Internal tides in the ocean. Rev. Geophys. Space Phys. 13, 167–182.
  * Zahn (1975) Zahn, J.-P. 1975 Tidal friction in close binary systems. A& A 57, 383–394.