1. Introduction
In rotating flows, inertial waves with frequency smaller than twice the rotation rate propagate at a fixed angle with respect to the rotation axis (Greenspan Reference Greenspan1968). The frequency and the angle are preserved when inertial waves reflect on a boundary. However, an inertial wave beam may contract or expand as it reflects. This linear contraction effect is responsible for inviscid singularities in the inertial wave field (Rieutord, Georgeot & Valdettaro Reference Rieutord, Georgeot and Valdettaro2001; Ogilvie Reference Ogilvie2020).
There are two types of inviscid singularities concerned in the present work. One is at the critical latitude of a sphere where the rays are tangent to the boundary and where Ekman pumping blows up (Roberts & Stewartson Reference Roberts and Stewartson1963). This singularity propagates within the fluid along the tangent critical line at the critical latitude (Kerswell Reference Kerswell1995). When regularised by viscosity, it forms concentrated internal shear layers around the critical line. The viscous self-similar solution of Moore & Saffman (Reference Moore and Saffman1969) and Thomas & Stevenson (Reference Thomas and Stevenson1972) is expected to describe the viscous structure of these thin layers for small Ekman numbers. For a librating spheroid, Le Dizès & Le Bars (Reference Le Dizès and Le Bars2017) derived the singularity strength and the amplitude of the self-similar solution by matching asymptotically the shear layer solution with the Ekman layer solution. The self-similar solution was found to be in agreement with direct numerical simulation. The same self-similar solution (with the singularity strength and the amplitudes derived in an open domain) was also used to describe the solution on a periodic orbit in a spherical shell geometry (He et al. Reference He, Favier, Rieutord and Le Dizès2022, hereafter HFRL22). In that case, the solution was obtained by considering its propagation along the periodic orbit for an infinite number of cycles. It was found to agree very well with the numerical solutions obtained for low Ekman numbers. In particular, both the internal shear layer structure and its amplitude scaling in 
$E^{1/12}$ were recovered by the numerical results using Ekman numbers as low as 
$10^{-10}$.
The singularity obtained from the critical latitude on the outer sphere gives rise to different internal shear layers. These layers are weaker and thicker, and do not possess a self-similar structure (Kerswell Reference Kerswell1995; Lin & Noir Reference Lin and Noir2021). Kida (Reference Kida2011) obtained their asymptotic structure for a precessing sphere.
Besides libration and precession, which drive the flows through viscosity, non-viscous forcing associated with translating or deforming bodies has also been analysed. Many studies have been performed in the context of stratified fluids for applications to tidal flows. Analytic results were obtained for the cylinder and the sphere in an unbounded geometry (Hurley Reference Hurley1997; Hurley & Keady Reference Hurley and Keady1997; Voisin Reference Voisin2003) and validated experimentally in both two (Sutherland & Linden Reference Sutherland and Linden2002; Zhang, King & Swinney Reference Zhang, King and Swinney2007) and three (Flynn, Onu & Sutherland Reference Flynn, Onu and Sutherland2003; Voisin, Ermanyuk & Flór Reference Voisin, Ermanyuk and Flór2011; Ghaemsaidi & Peacock Reference Ghaemsaidi and Peacock2013) dimensions. Hurley & Keady (Reference Hurley and Keady1997) and Voisin (Reference Voisin2003) also showed that in the far field, the solution takes the self-similar form predicted by Moore & Saffman (Reference Moore and Saffman1969). The singularity strength, however, varies with respect to the nature of the forcing. Machicoane et al. (Reference Machicoane, Cortet, Voisin and Moisy2015) discussed this effect for pulsating and oscillating spheres.
The other inviscid singularity is the attractor in a closed container onto which inertial waves tend to focus (Rieutord & Valdettaro Reference Rieutord and Valdettaro1997). The presence of such singularities is related to the hyperbolic character of the Poincaré equation describing the wave structure: it leads to an ill-posed Cauchy problem, except for a few geometries such as the cylinder or the ellipsoid (Rieutord, Georgeot & Valdettaro Reference Rieutord, Georgeot and Valdettaro2000). Contrary to the singularity attached to the boundary at the critical latitude, the attractor is a limit cycle that is not directly dependent on the nature of the forcing. Attractors also generate intense internal shear layers, as first observed in a trapezoidal tank for a stably stratified fluid (Maas et al. Reference Maas, Benielli, Sommeria and Lam1997). The asymptotic structure of these layers was analysed in a forced regime in two dimensions by Ogilvie (Reference Ogilvie2005) (hereafter O05). Under a few technical hypotheses, he was able to derive the functional equation describing the inviscid streamfunction and to provide the viscous asymptotic expression of the streamfunction close to the attractor. In particular, O05 showed that for his quadrilateral geometry possessing a unique attractor, the main contribution to the solution is associated with the logarithmic singularity of the inviscid streamfunction. We will use and adapt his results to our geometry. His results were confirmed by a numerical study of an inclined rotating square in Jouve & Ogilvie (Reference Jouve and Ogilvie2014).
In a spherical shell, there may exist both critical-latitude and attractor singularities at the same time. In HFRL22, we have considered a case where no attractor was present. We have assumed that the fluid was forced by librating the inner core at a frequency such that inertial waves propagated in a direction oriented at 
$45^{\circ }$ with respect to the rotation axis. All the ray trajectories were periodic in that case, and the (critical) path issued from the critical latitude on the inner core was just a rectangle in the upper left meridional cut of the shell. For other frequencies, the rays issued from the critical latitude are expected to perform a more complex pattern and possibly converge to an attractor (Tilgner Reference Tilgner1999; Ogilvie & Lin Reference Ogilvie and Lin2004; Ogilvie Reference Ogilvie2009). It is this situation that we want to address in the present work. We consider the same framework as in HFRL22, where local asymptotic solutions propagated in the volume are compared with global numerical results, but for a frequency for which an attractor is now present.
The paper is organised as follows. The framework is introduced in § 2. In § 2.1, we describe the three-dimensional (3-D) and two-dimensional (2-D) configurations that we have considered, and provide the governing equations. In § 2.2, the numerical method used to integrate the equations for each configuration is explained. In § 3, we first analyse the wave beams emitted from the critical latitude on the inner core. The asymptotic solution built by propagating the self-similar solution is compared to the numerical solution. Discrepancies are observed close to the attractors for some of the cases. In § 4, we then focus on the solution close to the attractors. We construct an asymptotic solution based on the theory of O05 for an attractor without phase shift in § 4.1, and provide a numerical validation in § 4.2. A brief conclusion is provided finally in § 5.
2. Framework
2.1. Configurations
In this paper, we consider the flow of an incompressible fluid of constant kinematic viscosity 
$\nu ^*$ rotating around the axis 
$\boldsymbol {e}_z$ with a uniform rotation rate 
$\varOmega ^*$. We consider two different configurations. The first is the axisymmetric flow filling a 3-D spherical shell, as in HFRL22. The other configuration is the 2-D flow, but with three velocity components, between two co-axial cylinders whose axis is horizontal, as in Rieutord, Valdettaro & Georgeot (Reference Rieutord, Valdettaro and Georgeot2002) and Rieutord & Valdettaro (Reference Rieutord and Valdettaro2010). In the following, geometries, governing equations and forcings are described separately for the two configurations.
2.1.1. Three-dimensional configuration
The geometry of the 3-D spherical shell is shown in figure 1(a), whose meridional plane can be found in figure 2 of HFRL22. The radii of the outer and inner spheres are 
$\rho ^*$ and 
$\eta \rho ^*$ (with 
$0<\eta <1$ the aspect ratio), respectively. Lengths are non-dimensionalised by the outer radius 
$\rho ^*$ such that the inner and outer dimensionless radii are 
$\eta$ and 
$1$, respectively. Time is non-dimensionalised by the angular period 
$1/\varOmega ^*$. The imposed harmonic forcing is the libration of one of the two boundaries, with amplitude 
$\varepsilon =\varepsilon ^*/\varOmega ^*$ (
$\varepsilon \ll 1$) and frequency 
$\omega =\omega ^*/\varOmega ^*$. The Ekman number is defined as
\begin{equation} E=\frac{\nu^*}{\varOmega^*\rho^{*2}} , \end{equation}
with 
$\nu ^*$ being the kinematic viscosity.
Figure 1. Configurations: (a) a 3-D spherical shell subject to the longitudinal libration on the inner core; (b) a 2-D cylindrical annulus subject to the symmetric forcing on the inner core; (c) a 2-D cylindrical annulus subject to the antisymmetric forcing on the inner core. The red arrows show the magnitudes and directions of the forcings at one instant.
As in HFRL22, we care about the linear harmonic response when the Ekman number is extremely small. We look for solutions that are harmonic in time:
\begin{equation} \varepsilon(\boldsymbol{v},p)\,{\rm e}^{-\mathrm{i}\omega t}+{\rm c.c.}, \end{equation}
with 
${\rm c.c.}$ denoting complex conjugation. The velocity 
$\boldsymbol {v}$ and pressure 
$p$ satisfy the linearised incompressible Navier–Stokes equations in the rotating frame:
\begin{gather} -\mathrm{i}\omega\boldsymbol{v}+2\boldsymbol{e}_z\times\boldsymbol{v}={-}\boldsymbol{\nabla} p+E\,\nabla^2\boldsymbol{v} , \end{gather}
\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}=0 . \end{gather}
In terms of the velocity components and pressure, the governing equations in the cylindrical coordinate system 
$(r,z,\phi )$ become
\begin{gather} -\mathrm{i}\omega v_r-2v_\phi+\frac{\partial p}{\partial r}-E\left(\nabla^2-\frac{1}{r^2}\right)v_r=0 , \end{gather}
\begin{gather}-\mathrm{i}\omega v_z+\frac{\partial p}{\partial z}-E\,\nabla^2v_z=0 , \end{gather}
\begin{gather}-\mathrm{i}\omega v_\phi+2v_r-E\left(\nabla^2-\frac{1}{r^2}\right)v_\phi=0\ , \end{gather}
\begin{gather}\frac{\partial v_r}{\partial r}+\frac{v_r}{r}+\frac{\partial v_z}{\partial z}=0 , \end{gather}with the Laplacian operator
\begin{equation} \nabla^2=\frac{\partial^2}{\partial r^2}+\frac{1}{r}\,\frac{\partial}{\partial r}+\frac{\partial^2}{\partial z^2} . \end{equation}The fluid is subject to no-slip boundary conditions on all boundaries. One of the two boundaries is subject to the longitudinal libration, as shown by the red arrows in figure 1(a), which corresponds to the oscillating solid body rotation of the boundary according to
\begin{equation} \boldsymbol{v}(\rho)=r\boldsymbol{e}_\phi \quad \text{at} \ \rho=\eta\ \mbox{or}\ 1, \end{equation}while the other boundary is not moving,
\begin{equation} \boldsymbol{v}(\rho)=\boldsymbol{0} \quad \mbox{at} \ \rho=1\ \mbox{or}\ \eta, \end{equation}
where 
$r$ is the distance to the rotation axis of the cylindrical coordinate system 
$(r,z,\phi )$, while 
$\rho =\sqrt {r^2+z^2}$ is the distance to the centre in the spherical coordinate system.
2.1.2. Two-dimensional configuration
We also consider a 2-D simplification of the 3-D axisymmetric configuration discussed above. The geometry can be viewed as a slender cored torus with the principal radius tending to infinity (Rieutord et al. Reference Rieutord, Valdettaro and Georgeot2002; Rieutord & Valdettaro Reference Rieutord and Valdettaro2010), which is effectively equivalent to two co-axial cylinders whose principal axis is horizontal, as shown in figures 1(b,c). The flow between the two cylinders satisfies governing equations similar to (2.3), while the curvature terms in the differential operators, such as 
$1/r$, 
$(1/r)(\partial /\partial r)$ and 
$1/r^2$, are omitted. Explicitly, in terms of the velocity components and pressure, the governing equations are
\begin{gather} -\mathrm{i}\omega v_x-2v_y+\frac{\partial p}{\partial x}-E\,\nabla^2 v_x=0, \end{gather}
\begin{gather}-\mathrm{i}\omega v_z+\frac{\partial p}{\partial z}-E\,\nabla^2 v_z=0, \end{gather}
\begin{gather}-\mathrm{i}\omega v_y+2v_x-E\,\nabla^2 v_y=0, \end{gather}
\begin{gather}\frac{\partial v_x}{\partial x}+\frac{\partial v_z}{\partial z}=0, \end{gather}with the Laplacian operator
\begin{equation} \nabla^2=\partial^2/\partial x^2+\partial^2/\partial z^2. \end{equation}
We use 
$(x,y,z)$ to denote the Cartesian coordinates, where 
$Ox$ and 
$Oz$ are the horizontal and vertical axes, respectively, and 
$Oy$ is along the direction perpendicular to the 
$Oxz$ plane, as shown in figures 1(b,c). Note that although we use the same symbol for the Laplacian operators in two and three dimensions, there is no ambiguity since the 2-D and 3-D operators are used independently in the corresponding dimensions.
Similar to libration in the 3-D configuration, the imposed forcing should be on the boundary. We also require that it drives the fluid in the bulk through viscous coupling only. The direction of the forcing is thus aligned with that of 
$\boldsymbol {e}_y$ perpendicular to the 
$Oxz$ plane. We consider two options for the amplitude of the forcing. One option is that the amplitude is a constant, which is
\begin{equation} \boldsymbol{v}(\varrho)=\boldsymbol{e}_y \quad \mbox{at} \ \varrho=\eta\ \mbox{or}\ 1, \end{equation}
where 
$\varrho =\sqrt {x^2+z^2}$. The cylinder subject to this forcing is expected to oscillate uniformly along the direction 
$\boldsymbol {e}_y$, as shown by the red arrows in figure 1(b). The other option is that the amplitude of the forcing depends linearly on the horizontal coordinate 
$x$, which is
\begin{equation} \boldsymbol{v}(\varrho)=x\boldsymbol{e}_y \quad \mbox{at} \ \varrho=\eta\ \mbox{or}\ 1. \end{equation}
The cylinder subject to this forcing oscillates non-uniformly, inducing shear at the inner boundary, as shown by the red arrows in figure 1(c). While unrealistic from an experimental point of view, it is a mathematically well-posed boundary condition and provides another symmetry, as discussed later. While the formula for the 2-D antisymmetric forcing (2.11) is similar to the 3-D libration case (2.6), they differ in that the horizontal coordinate 
$x$ in the 2-D configuration can be negative.
Both forcings are symmetric about the horizontal axis 
$Ox$. However, the former forcing (2.10) is symmetric about the vertical axis 
$Oz$, while the latter (2.11) is antisymmetric about 
$Oz$; see the red arrows in figures 1(b) and 1(c), respectively. These two forcings are thus referred to as symmetric and antisymmetric forcings, respectively, according to their symmetries about the 
$Oz$ axis. They are also imposed on one of the two boundaries, while the other boundary condition is no-slip.
In summary, we consider three different forcings, which are referred to as the 3-D libration (2.6), 2-D symmetric (2.10) and 2-D antisymmetric (2.11) forcings. The first is defined in the 3-D spherical shell, while the latter two correspond to the 2-D cylindrical annulus. Note that we restrict our study to purely axisymmetric situations, so that we ignore zonally propagating waves that require azimuthal inhomogeneities as discussed by Rabitti & Maas (Reference Rabitti and Maas2013).
2.2. Numerical methods
The governing equations (2.3) are solved numerically by spectral methods for both the 3-D and 2-D configurations. We actually solve the vorticity equation, which is the curl of the momentum equations (2.3a):
\begin{equation} -\mathrm{i}\omega\boldsymbol{\nabla}\times\boldsymbol{v}+2\boldsymbol{\nabla}\times (\boldsymbol{e}_z\times\boldsymbol{v})=E\,\boldsymbol{\nabla}\times(\nabla^2\boldsymbol{v}). \end{equation}
In the 2-D configuration, the curl is taken only in the 
$Oxz$ plane. The numerical methods are different for the two configurations. Therefore, they are presented separately.
2.2.1. Three-dimensional configuration
In the 3-D configuration, the numerical method is similar to that in our former work (HFRL22). The governing equations are solved in the spherical coordinates 
$(\rho,\theta,\phi )$, with 
$\rho$ the distance to the centre, 
$\theta$ the colatitude, and 
$\phi$ the azimuthal angle. The velocity is expanded onto the vector spherical harmonics in the angular directions:
\begin{equation} \boldsymbol{v}=\sum_{l=0}^{+\infty}\sum_{m={-}l}^{{+}l}u^l_m(\rho)\,\boldsymbol{R}_l^m+v^l_m(\rho)\,\boldsymbol{S}_l^m+w^l_m(\rho)\,\boldsymbol{T}_l^m, \end{equation}with
\begin{equation} \boldsymbol{R}_l^m=Y_l^m(\theta,\phi)\,\boldsymbol{e}_\rho, \quad \boldsymbol{S}_l^m=\boldsymbol{\nabla} Y_l^m, \quad T_l^m=\boldsymbol{\nabla}\times\boldsymbol{R}_l^m. \end{equation}
The gradients are taken on the unit sphere. The vorticity equation (2.12) is projected onto the basis, and 
$u^l$ and 
$w^l$ satisfy a set of ordinary differential equations
\begin{gather} E\varDelta_lw^l+\mathrm{i}\omega w^l={-}2A_l\rho^{l-1}\,\frac{\partial}{\partial \rho}\left(\frac{u^{l-1}}{\rho^{l-2}}\right)-2A_{l+1}\rho^{{-}l-2}\,\frac{\partial}{\partial \rho}\left(\rho^{l+3}u^{l+1}\right), \end{gather}
\begin{gather}E\varDelta_l\varDelta_l(\rho u^l)+\mathrm{i}\omega\varDelta_l(\rho u^l)=2B_l\rho^{l-1}\,\frac{\partial}{\partial \rho}\left(\frac{w^{l-1}}{\rho^{l-1}}\right)+2B_{l+1}\rho^{{-}l-2}\,\frac{\partial}{\partial \rho}\left(\rho^{l+2}w^{l+1}\right), \end{gather}with
\begin{equation} A_l=\frac{1}{l^2\sqrt{4l^2-1}}, \quad B_l=l^2(l^2-1)A_l, \quad \varDelta_l=\frac{\mathrm{d}^2}{\mathrm{d}\rho^2}+\frac{2}{\rho}\,\frac{\mathrm{d}}{\mathrm{d}\rho}-\frac{l(l+1)}{\rho^2} \end{equation}
(e.g. Rieutord Reference Rieutord1991). Axisymmetry (
$m=0$) is employed. Here, 
$v^l$ is related to 
$u^l$ through the continuity equation
\begin{equation} v^l=\frac{1}{\rho l(l+1)}\,\frac{\mathrm{d}\rho^2u^l}{\mathrm{d}\rho}. \end{equation}One of the two boundaries is subject to the no-slip boundary condition
\begin{equation} w^l=u^l=\frac{\mathrm{d}u^l}{\mathrm{d}\rho}=0 \quad \mbox{at} \ \rho=1\ \mbox{or}\ \eta. \end{equation}The other boundary is subject to the libration (2.6), whose projection onto the spherical harmonics yields the inhomogeneous boundary condition
\begin{equation} w^l=2\,\sqrt{\frac{\rm \pi}{3}}\,\rho\delta_{1,l}, \quad u^l=\frac{\mathrm{d}u^l}{\mathrm{d}\rho}=0 \quad \mbox{at} \ \rho=\eta \ \mbox{or}\ 1, \end{equation}
where 
$\delta _{1,l}$ is the Kronecker symbol. Note that the libration is imposed on the spherical harmonic degree 
$l=1$.
Equations (2.15)–(2.19) are truncated to the spherical harmonic degree 
$L$. The derivatives with respect to the radial coordinate 
$\rho$ are replaced by the Chebyshev differentiation matrices at 
$N+1$ collocation points of the Gauss–Lobatto grid. Then a block tridiagonal system is obtained as
\begin{equation} \begin{bmatrix} \boldsymbol{D}_1 & \boldsymbol{C}_1 & & & \\ \boldsymbol{B}_1 & \boldsymbol{D}_2 & \boldsymbol{C}_2 & & \\ & {\ddots} & {\ddots} & {\ddots} & \\ & & \boldsymbol{B}_{L-1} & \boldsymbol{D}_{L-1} & \boldsymbol{C}_{L-1} \\ & & & \boldsymbol{B}_L & \boldsymbol{D}_L \end{bmatrix} \begin{bmatrix} \boldsymbol{w}^1\\ \boldsymbol{\rho u}^2\\ {\vdots}\\ \boldsymbol{w}^{L-1}\\ \boldsymbol{\rho u}^{L} \end{bmatrix}= \begin{bmatrix} \boldsymbol{b}_1\\ \boldsymbol{b}_2\\ {\vdots}\\ \boldsymbol{b}_{L-1}\\ \boldsymbol{b}_{L} \end{bmatrix}. \end{equation}
The blocks within the coefficient matrix and the vectors are 
$(N+1)\times (N+1)$ and 
$(N+1)\times 1$, respectively. The order of the coefficient matrix is 
$(N+1)L$, and the number of non-zero elements is 
$(N+1)^2(3L-2)$. This block tridiagonal system is usually solved by an LU solver (Rieutord & Valdettaro Reference Rieutord and Valdettaro1997), by which the coefficient matrix is stored in the banded matrix format, and the number of elements in memory is 
$(N+1)^2(4L-4)-(N+1)(L-2)$. On the other hand, the block tridiagonal system can be solved by the block version of the standard tridiagonal algorithm (also called the Thomas algorithm), which is Gaussian elimination on a block tridiagonal system. This method has been utilised by Ogilvie & Lin (Reference Ogilvie and Lin2004). The algorithm can be found in Engeln-Mèullges & Uhlig (Reference Engeln-Mèullges and Uhlig1996, p. 121). The elimination is advanced forwards from the lowest spherical harmonic degree to the highest, and the block tridiagonal matrix is reduced to a block upper bidiagonal one, then the solution is obtained by backward substitution. During the forward elimination, the updated diagonal block 
$\boldsymbol {D}_l$ is factorised by the LU solver. A partial pivoting of the block is employed in order to improve the numerical stability.
The three blocks 
$\boldsymbol {B}_l$, 
$\boldsymbol {D}_l$ and 
$\boldsymbol {C}_l$, and the inhomogeneous term 
$\boldsymbol {b}_l$ at the spherical harmonic degree 
$l$, are needed only when they take part in the forward elimination. Hence the storage of the whole coefficient matrix is unnecessary. However, all the updated super-diagonal blocks 
$C_l$ should be reserved in memory for the backward substitution. Their size is 
$(N+1)^2(L-1)$, which is almost one-third of that of non-zero elements in the original coefficient matrix, and one-quarter of that in the banded matrix format required by the global LU solver. Therefore, the memory usage of the block tridiagonal algorithm is much less than that of the global LU solver, especially when 
$L$ and 
$N$ are very large, as required for very low Ekman numbers. We develop a code based on the block tridiagonal algorithm using the efficient dynamic programming language Julia (Bezanson et al. Reference Bezanson, Edelman, Karpinski and Shah2017). For now, we can reach 
$E=10^{-11}$ by using 
$8000$ spherical harmonics and 
$2500$ Chebyshev polynomials using double-precision floating-point format. The memory footprint is around 750 GB.
2.2.2. Two-dimensional configuration
In the 2-D configuration, we take the numerical method similar to that adopted by Rieutord et al. (Reference Rieutord, Valdettaro and Georgeot2002) and Rieutord & Valdettaro (Reference Rieutord and Valdettaro2010). The vorticity equation (2.12) is solved in the polar coordinates 
$(\varrho,\vartheta )$, with 
$\varrho$ the distance to the centre, and 
$\vartheta$ the angle measured from the horizontal axis 
$Ox$. In terms of the streamfunction 
$\psi$ and the associated variable 
$\chi$,
\begin{equation} v_\varrho={-}\frac{1}{\varrho}\,\frac{\partial\psi}{\partial\vartheta}, \quad v_\vartheta=\frac{\partial\psi}{\partial\varrho}, \quad v_y=\chi, \end{equation}and the vorticity equation is recast to
\begin{gather} -\mathrm{i}\omega\,\nabla^2\psi+2\left(\sin{\vartheta}\,\frac{\partial\chi}{\partial \varrho}+\frac{\cos{\vartheta}}{\varrho}\,\frac{\partial\chi}{\partial \vartheta}\right)-E\,\nabla^4\psi=0, \end{gather}
\begin{gather}-\mathrm{i}\omega\chi-2\left(\sin{\vartheta}\,\frac{\partial\psi}{\partial \varrho}+\frac{\cos{\vartheta}}{\varrho}\,\frac{\partial\psi}{\partial \vartheta}\right)-E\,\nabla^2\chi=0, \end{gather}with the operator
\begin{equation} \nabla^2=\frac{\partial^2}{\partial \varrho^2}+\frac{1}{\varrho}\,\frac{\partial}{\partial \varrho}+\frac{1}{\varrho^2}\,\frac{\partial^2}{\partial \vartheta^2}. \end{equation}
The streamfunction 
$\psi$ and the associated variable 
$\chi$ are expanded by Fourier series in the angular direction as
\begin{equation} \psi=\sum_{l={-}\infty}^{+\infty}\psi_l(\varrho)\,{\rm e}^{\mathrm{i}l\vartheta}, \quad \chi={-}\mathrm{i}\sum_{l={-}\infty}^{+\infty}\chi_l(\varrho)\,{\rm e}^{\mathrm{i}l\vartheta}. \end{equation}The projection of the governing equations (2.22) onto this basis is
\begin{gather} \mathrm{i}\omega\,\nabla_l^2\psi_l+(\chi_{l-1}^\prime-\chi_{l+1}^\prime)-\frac{1}{\varrho}\left[(l-1)\chi_{l-1}+(l+1)\chi_{l+1}\right]+E\,\nabla_l^4\psi_l=0, \end{gather}
\begin{gather}\mathrm{i}\omega\chi_l+(\psi_{l-1}^\prime-\psi_{l+1}^\prime)-\frac{1}{\varrho}\left[(l-1)\psi_{l-1}+(l+1)\psi_{l+1}\right]+E\,\nabla_l^2\chi_l=0, \end{gather}with
\begin{equation} \nabla_l^2=\frac{\mathrm{d}^2}{\mathrm{d}\varrho^2}+\frac{1}{\varrho}\,\frac{\mathrm{d}}{\mathrm{d}\varrho}-\frac{l^2}{\varrho^2}. \end{equation}The unforced boundary is subject to the no-slip boundary condition
\begin{equation} \psi_l=\frac{\mathrm{d}\psi_l}{\mathrm{d}\varrho}=\chi_l=0 \quad \mbox{at} \ \varrho=1 \ \mbox{or}\ \eta. \end{equation}
The other boundary is subject to the viscous boundary forcings (2.10) and (2.11). Both forcings are symmetric about the horizontal axis 
$Ox$, which leads to
\begin{equation} \psi_{{-}l}={-}\psi_l, \quad \chi_{{-}l}=\chi_l. \end{equation}Only the non-negative Fourier components are necessary to be computed. The symmetric forcing (2.10) imposes the boundary condition
\begin{equation} \psi_l=\frac{\mathrm{d}\psi_l}{\mathrm{d}\varrho}=0, \quad \chi_l=\mathrm{i}\delta_{0,l} \quad \mbox{at} \ \varrho=\eta\ \mbox{or}\ 1. \end{equation}
Note that the forcing is imposed at 
$l=0$. Therefore, the following Fourier components are excited:
\begin{equation} \chi_0,\psi_1,\chi_2,\psi_3,\ldots. \end{equation}On the other hand, the antisymmetric forcing (2.11) imposes the boundary condition
\begin{equation} \psi_l=\frac{\mathrm{d}\psi_l}{\mathrm{d}\varrho}=0, \quad \chi_l=\mathrm{i}\,\frac{\varrho}{2}\,\delta_{1,l} \quad \mbox{at} \ \varrho=\eta\ \mbox{or}\ 1. \end{equation}
Note that the forcing is imposed at 
$l=1$. Therefore, the following Fourier components are excited:
\begin{equation} \chi_1,\psi_2,\chi_3,\psi_4,\ldots. \end{equation} As in the 3-D configuration, the equations are truncated at the Fourier component 
$L$, and the derivatives to 
$\varrho$ are replaced by the Chebyshev differentiation matrices with order 
$N+1$. The resulting block tridiagonal system is solved by the same block tridiagonal algorithm as before.
The verification of the two spectral codes used in this paper can be found in Appendix A.
3. Wave beams from the critical latitude on the inner core
The aforementioned forcings are imposed on the inner core. The forcing frequency 
$\omega$ is chosen in the inertial range such that inertial waves propagate at an inclined angle 
$\theta _c=\arccos {\omega /2}$ relative to the horizontal plane. As in HFRL22, two concentrated wave beams are expected to be generated from the critical latitude localised at 
$(r,z)=(\eta \sqrt {1-\omega ^2/4},\eta \omega /2)$ on the inner core. These wave beams travel along the tangential line at the critical latitude in two opposite directions (northwards and southwards), reflect on the boundaries, and form a ray pattern in the spherical shell geometry. In general, for a fixed inclined angle 
$\theta _c$, any ray pattern is composed of the four rays with opposite propagation directions, which are referred to as the northward, outward, southward and inward, as shown in figure 2. In HFRL22, we considered the case where the ray pattern is a simple periodic pattern. Here, we consider a more general situation where the wave beams converge towards an attractor. Our first objective is to analyse whether an asymptotic solution can be constructed by propagating the self-similar solution describing the concentrated wave beam emitted from the critical latitude, as it was done in HFRL22.
Figure 2. Four propagation directions of the rays in a closed domain. The local vectors attached to each ray are the orientations of the local frames 
$(x_\parallel,x_\perp )$.
In § 3.1, the asymptotic theory is presented. The properties of the self-similar solution and of the reflection laws are first recalled and adapted to the 2-D configurations that we also consider before analysing the propagation towards the attractor. The asymptotic solution is then compared to numerical results in § 3.2.
3.1. Asymptotic theory
3.1.1. Viscous self-similar solution and scaling
The concentrated ray beams emitted from the critical latitude are associated with an inviscid singularity along the critical ray (Le Dizès Reference Le Dizès2023). It is the viscous smoothing of this singularity that gives rise in the limit of small Ekman numbers to a self-similar expression for the dominant wave beam velocity components (Moore & Saffman Reference Moore and Saffman1969).
The natural way to describe this self-similar solution is to introduce the local coordinates 
$(x_\parallel,x_\perp )$ on the critical ray path, with 
$x_\parallel$ measuring the travelled distance from the source along the critical ray, and 
$x_\perp$ measuring the displacement relative to the critical ray (
$x_\perp =0$ is the critical ray equation). The orientation of 
$x_\perp$ is chosen as indicated in figure 2. It is assumed not to change during the beam propagation.
The wave beam is centred on the critical ray and has width of order 
$E^{1/3}$. In the 
$(r,z)$ plane, its main velocity component is oriented along 
${\boldsymbol {e}_\parallel }$ and can be written at leading order in 
$E^{1/3}$ in the 3-D axisymmetric geometry as (see details in Le Dizès & Le Bars Reference Le Dizès and Le Bars2017)
\begin{equation} v_\parallel = \frac{1}{\sqrt{r}}\,C_0\,H_m(x_\parallel,x_\perp)=\frac{1}{\sqrt{r}}\,C_0\left(\frac{x_\parallel}{2\sin\theta_c}\right)^{{-}m/3}h_m(\zeta), \end{equation}with the similarity variable
\begin{equation} \zeta=x_\perp E^{{-}1/3}\left(\frac{2\sin\theta_c}{x_\parallel}\right)^{1/3}, \end{equation}and the special function introduced by Moore & Saffman (Reference Moore and Saffman1969),
\begin{equation} h_m(\zeta)=\frac{{\rm e}^{-\mathrm{i}m{\rm \pi}/2}}{(m-1)!}\int^{+\infty}_0{\rm e}^{\mathrm{i}p\zeta-p^3}p^{m-1}\,{\rm d}p. \end{equation}
The parameters 
$C_0$ and 
$m$ denote the amplitude and singularity strength, respectively, which will be specified below. Note that this meaning of 
$m$ should not be confused with the order of the spherical harmonics in the spectral expansion (2.13), which is not used here since we focus on purely axisymmetric solutions. The velocity across the critical rays 
$v_\perp$ and the pressure 
$p$ are 
$O(E^{1/3})$ smaller. However, the wave beam has a velocity component normal to the 
$(r,z)$ plane of the same order, which is given by (see Rieutord et al. Reference Rieutord, Georgeot and Valdettaro2001; Le Dizès & Le Bars Reference Le Dizès and Le Bars2017)
\begin{equation} v_\phi={\pm}\mathrm{i} v_\parallel. \end{equation}
The sign corresponds to the sign of the projection of the local unit vector 
$\boldsymbol {e}_\parallel$ onto the global unit vector 
$\boldsymbol {e}_r$. For the northward and inward rays, the sign is 
$-$; for the southward and outward rays, the sign is 
$+$ (see figure 2).
The inviscid singularity that gives rise to the self-similar viscous solution is recovered by taking the limit 
$\zeta \rightarrow \infty$ in (3.1):
\begin{equation} v_\parallel\rightarrow \frac{1}{\sqrt{r}}\,C_0\,x_\perp^{{-}m}E^{m/3} \quad \mbox{as} \ \zeta \rightarrow +\infty. \end{equation} As we will see, it is also useful to introduce the streamfunction 
$\psi$ that can be defined for axisymmetric flows by
\begin{equation} v_\parallel = \epsilon\,\frac{1}{r}\,\frac{\partial \psi}{\partial x_\perp}, \quad v_\perp = {-}\epsilon\,\frac{1}{r}\,\frac{\partial \psi}{\partial x_\parallel}, \end{equation}
where 
$\epsilon =1$ for the rays propagating northwards and southwards, and 
$\epsilon =-1$ for the rays propagating inwards and outwards (see figure 2). Equation (3.6a) can be integrated to give at leading order
\begin{equation} \psi=\epsilon\sqrt{r}\,\frac{C_0E^{1/3}}{m-1}\,H_{m-1}(x_\parallel,x_\perp). \end{equation}
Note that the streamfunction 
$\psi$ is 
$E^{1/3}$ smaller than the parallel velocity 
$v_\parallel$.
The above expressions are valid for 3-D axisymmetric geometries. For 2-D configurations, the term 
$\sqrt {r}$ is not present in the velocity and streamfunction expressions. We get
\begin{gather} v_\parallel^{(2\text{-}D)}=C_0\,H_m(x_\parallel,x_\perp), \end{gather}
\begin{gather}\psi^{(2\text{-}D)}=\epsilon\,\frac{C_0E^{1/3}}{m-1}\,H_{m-1}(x_\parallel,x_\perp). \end{gather} The velocity component 
$v_y$ perpendicular to the 
$(x,z)$ plane differs from 
$v_\parallel ^{(2\text {-}D)}$ by a 
$\pm {\rm \pi}/2$ phase factor, as the relation (3.4) between 
$v_\phi$ and 
$v_\parallel$ in three dimensions.
In the self-similar solution (3.1), the free parameters, the singularity strength 
$m$ and the amplitude 
$C_0$ depend on the nature of the forcing. For a viscous forcing – that is, a forcing induced by Ekman pumping – these parameters can be obtained in closed form for the northward and southward beams generated from the critical latitude (Le Dizès & Le Bars Reference Le Dizès and Le Bars2017; Le Dizès Reference Le Dizès2023). For a librating sphere, they are given by (Le Dizès & Le Bars Reference Le Dizès and Le Bars2017)
\begin{equation} m=5/4, \end{equation}and
\begin{gather} C_0=\frac{E^{1/12}}{8(2\sin\theta_c)^{1/2}(2{/\eta})^{1/4}}\,{\rm e}^{\mathrm{i}{\rm \pi}/2} \quad \mbox{for the northward beam}, \end{gather}
\begin{gather}C_0=\frac{E^{1/12}}{8(2\sin\theta_c)^{1/2}(2{/\eta})^{1/4}}\,{\rm e}^{\mathrm{i}3{\rm \pi}/4} \quad \mbox{for the southward beam}. \end{gather}
These expressions can be applied to our geometry for the three forcings (2.6), (2.10) and (2.11) imposed on the inner core. Considering the different non-dimensionalisation of lengths adopted by Le Dizès & Le Bars (Reference Le Dizès and Le Bars2017) and this work (the radial distance of the critical latitude to the rotation axis versus the outer radius), the absolute value of the complex amplitude 
$C_0$ should be adapted as indicated in table 1 for the three forcings. The factor 
$\eta \sin {\theta _c}$ is the distance of the critical latitude to the axis 
$Oz$.
Table 1. Absolute value of the complex amplitude 
$C_0$ for different forcings.
Note that the amplitude 
$C_0$ of the parallel velocity scales as 
$E^{1/12}$. This scaling has been validated by HFRL22 for Ekman numbers down to 
$10^{-10}$. The amplitude of the streamfunction is weaker and of order 
$E^{5/12}$.
3.1.2. Reflections on the boundaries and on the axis
The reflection of a self-similar wave beam on a boundary has been discussed in Le Dizès (Reference Le Dizès2020) and HFRL22. Le Dizès (Reference Le Dizès2020) showed that the wave beam keeps its self-similar form when it reflects on a boundary. More precisely if the incident beam is written as 
$v_\parallel ^{(i)}=C_0^{(i)}\,H_m(x_\parallel ^{(i)},x_\perp ^{(i)})$, then the reflected beam can also be written as 
$v_\parallel ^{(r)}=C_0^{(r)}\,H_m(x_\parallel ^{(r)},x_\perp ^{(r)})$, with
\begin{equation} \frac{x_{{\parallel} _b}^{(r)}}{x_{{\parallel} b}^{(i)}}=\alpha^3, \quad \frac{C_0^{(r)}}{C_0^{(i)}}=\alpha^{m-1}, \end{equation}
where the subscript 
$b$ indicates values taken at the reflection point. The reflection factor 
$\alpha$ at the reflection point is given by
\begin{equation} \alpha = \frac{\sin \theta^{(r)}}{\sin\theta^{(i)}}, \end{equation}
where 
$\theta ^{(r)}$ and 
$\theta ^{(i)}$ are the angles of the reflected and incident beams with respect to the boundary (see figure 2). This factor is smaller than 1 (resp. larger than 1) when there is a contraction (resp. expansion) of the beam. A reflection on a boundary then just modifies the travelled distance from the source and the amplitude of the beam. In particular, it has no effect on its phase.
Note, however, that this reflection law assumes implicitly that the beam is not forced at the boundary where it reflects. In particular, this implies a simple relation on the streamfunction of the incident and reflected beams at the boundary that can be written as
\begin{equation} \psi^{(r)} (x_{{\parallel} b}^{(r)},x_{{\perp} b}^{(r)}) +\psi ^{(i)}(x_{\parallel_b}^{(i)},x_{{\perp} b}^{(i)}) =0. \end{equation}We will see below that this relation is no longer valid when we get very close to an attractor.
The crossing of the wave beam with the rotation axis is of different nature. In the 3-D axisymmetric geometry, the self-similar solution diverges on the axis, but it can nevertheless be continued as if there were a reflection. The relation between the incident and reflected beams is obtained by a matching condition with the solution obtained close to the axis (see Le Dizès & Le Bars Reference Le Dizès and Le Bars2017). In that case, we obtain a phase shift 
${\rm \pi} /2$ between the reflected and incident beams:
\begin{equation} x_{{\parallel} _b}^{(r)} = x_{{\parallel} b}^{(i)}, \quad C_0^{(r)}={\rm e}^{\mathrm{i}\,\varphi}C_0^{(i)}, \end{equation}
with 
$\varphi ={\rm \pi} /2$.
In the 2-D configurations, the condition of reflection to apply on the axis 
$Oz$ is related directly to the property of symmetry of the forcing. On the axis 
$Oz$, the projections of propagation directions of the incident and reflected rays onto the global unit vector 
$\boldsymbol {e}_x$ are of opposite sign. According to the formula (3.4), we then have the relations
\begin{equation} v_y^{(r)}={\pm}\mathrm{i} v_\parallel^{(r)}, \quad v_y^{(i)}={\mp}\mathrm{i} v_\parallel^{(i)}. \end{equation}
For the 2-D symmetric forcing (2.10) where 
$v_y$ is forced in a symmetric way about the axis 
$Oz$, we have 
$v_y^{(r)}=v_y^{(i)}$. Therefore, the parallel velocities are of opposite sign, which means that
\begin{equation} \varphi={\rm \pi} \end{equation}
in (3.14). For the 2-D antisymmetric forcing (2.11) with 
$v_y^{(r)}=-v_y^{(i)}$, the parallel velocity is unchanged, which means that
\begin{equation} \varphi= 0. \end{equation} In order to consider a quarter of the domain in the 
$(r,z)$ or 
$(x,z)$ plane, the horizontal axis 
$Or$ (or 
$Ox$) also has to be considered as a place of reflection. Applying the same approach, we can show easily that no phase shift is created between reflected and incident beams on this axis for all the three forcings.
3.1.3. Propagation of critical-latitude beams
Having provided the structure of the wave beam and how it reflects on the boundaries and the axis, we are now in a position to analyse its propagation in a closed geometry. As explained above, we consider a frequency such that the rays emitted from the critical latitude on the inner core end up on an attractor. Our objective is to obtain the property of the self-similar beam centred on the critical ray as it moves towards the attractor. An example of a critical ray is shown in figure 3, where the ray (blue lines) propagates northwards from the critical latitude and spirals into one side of the attractor (red lines). In the following, we use this figure for explanation purposes, but the methodology is applicable for any type of wave pattern.
Figure 3. Schematic of propagation of a critical ray towards an attractor for 
$\eta =0.35$ and 
$\omega =0.8317$. The symbol 
$\ast$ denotes an arbitrary point, from which the local coordinates of each segment 
$(0,n)$ are measured.
The reflection positions on the axes and the boundaries during every loop are indicated as 
$P_{j,n}$, where 
$j$ denotes the reflection position and ranges from 
$0$ to 
$J-1$ (
$J=8$ in figure 3). The index 
$n$ denotes the number of the cycle and ranges from 
$1$ to 
$\infty$. For example, the reflection points on the rotation axis are 
$P_{1,1}, P_{1,2}, \ldots$ and 
$P_{1,\infty }$. To simplify the formula, we assume that the initial point of a cycle is the position 
$J$ of the former cycle, that is, 
$P_{0,n+1}= P_{J,n}$. The critical latitude corresponds to 
$P_{0,1}$. The critical ray follows the following path during propagation:
\begin{equation} P_{0,1}P_{1,1}\cdots P_{J-1,1} \rightarrow P_{0,2}P_{1,2}\cdots P_{J-1,2} \rightarrow \cdots \rightarrow P_{0,\infty}P_{1,\infty}\cdots P_{J-1,\infty} . \end{equation}
The critical ray ends up on the attractor denoted by 
$P_{0,\infty }P_{1,\infty }\cdots P_{J-1,\infty }$ after an infinite number of cycles.
The solution obtained by propagating the self-similar beam along the critical ray is expected to be composed of as many contributions as the number of segments between two reflection points. We use the subscript 
$(j,n)$ to denote the parameters associated with the segment 
$P_{j,n}P_{j+1,n}$ (with 
$j$ between 0 and 
$J-1$). Finding the parameters characterising this contribution requires tracking the variation of the travelled distance and of the amplitude during all the previous reflections. For this purpose, it is useful to write the travelled distance 
$x_{\parallel (\,j,n)}$ as
\begin{equation} x_{{\parallel} (j,n)} = L_{j,n}^{(s)} + x_{{\parallel} (j,n)}', \end{equation}
where 
$x_{\parallel (j,n)}'$ is the distance from 
$P_{j,n}$, and 
$L_{j,n}^{(s)}$ is the distance of the ‘virtual’ source 
$P_{j,n}^{(s)}$ from 
$P_{j,n}$. The condition of reflection (3.11a) applied in 
$P_{j+1,n}$ implies that
\begin{equation} L_{j+1,n}^{(s)} = (L_{j,n}^{(s)} + L_{j,n})\alpha_{j+1,n}^3, \end{equation}
where 
$L_{j,n}$ is the length of the segment 
$(j,n)$, and 
$\alpha _{j+1,n}$ is the reflection factor at 
$P_{j+1,n}$. Concerning the amplitude 
$C_{j,n}$ of the self-similar solution, we obtain from (3.11b) with (3.9) that
\begin{equation} C_{j+1,n} = C_{j,n} \alpha_{j+1,n}^{1/4}\,{\rm e}^{\mathrm{i} \varphi_{j+1}}, \end{equation}
where 
$\varphi _j$ is the phase shift obtained at the reflection at 
$P_{j,n}$. For the critical ray shown in figure 3, this phase shift is null except for 
$j=1$ (because the reflection is on the axis), for which it can be 
${\rm \pi} /2$ (3-D case), 
${\rm \pi}$ (2-D symmetric case) or 0 (2-D antisymmetric case).
In the following, we will consider the solution in a section perpendicular to the segments 
$(0,n)$. It is therefore useful to consider the evolution of the beam after each cycle for this particular segment as a function of 
$n$. Using (3.20), we can write
\begin{equation} L_{0,n+1}^{(s)} \equiv L_{J,n}^{(s)} = (L_{0,n}^{(s)} + \varLambda _{n} )\alpha_n^3, \end{equation}with
\begin{equation} \varLambda_n=L_{0,n}+\frac{L_{1,n}}{\alpha_{1,n}^3}+\frac{L_{2,n}}{\alpha_{1,n}^3\alpha_{2,n}^3}+\cdots+\frac{L_{J-1,n}}{\alpha_{1,n}^3\alpha_{2,n}^3\cdots\alpha_{J-1,n}^3} \end{equation}and
\begin{equation} \alpha_{n}=\alpha_{1,n}\alpha_{2,n}\cdots\alpha_{J,n} . \end{equation}Similarly, we obtain
\begin{equation} C_{0,n+1} \equiv C_{J,n} = C_{0,n} \alpha _n^{1/4}\,{\rm e}^{\mathrm{i} \varphi}, \end{equation}with
\begin{equation} \varphi=\varphi_1+\varphi_2+\cdots+\varphi_J. \end{equation}
Note that 
$\alpha _{J,n}=\alpha _{0,n+1}$ and 
$\varphi _J=\varphi _0$. For the first segment of the first cycle, the source is at 
$P_{0,1}$, so 
$L_{0,1}=0$ and the amplitude 
$C_{0,1}$ is given by the expression (3.10) of 
$C_0$.
Although a given parameter 
$\alpha _{j,n}$ can be larger than 1, the product (3.24) that defines 
$\alpha _n$ is necessarily smaller than 1 (for 
$n$ sufficiently large) because the critical ray converges towards an attractor. Its limit value 
$\alpha _{\infty }$ corresponds to the contraction factor of the attractor. The amplitude of the beam therefore goes rapidly to zero as one gets close to the attractor. This guarantees that although the various contributions superimpose on each other close to the attractor, the sum will remain finite on the attractor. The expression obtained by summing all the contributions coming from the segments 
$(0,n)$ with 
$n$ ranging from 1 to 
$\infty$ is then well defined. It can be written as
\begin{equation} v_\parallel \sim \sum_{n=1}^{\infty}v_{{\parallel} (0,n)}, \quad \psi\sim\sum_{n=1}^{\infty}\psi_{0,n} \end{equation}
for the parallel velocity and the streamfunction, respectively. These expressions are expected to provide an asymptotic solution close to segments 
$(0,n)$. In the following, they will be referred to as the critical-latitude solution. In the next section, they are plotted and compared to numerical solutions.
3.2. Results
The numerical solutions are obtained for Ekman numbers as low as 
$10^{-11}$, for which the scale separation between the wave beams and the domain size is clear. For simplicity, the velocity components 
$v_\phi$ in three dimensions and 
$v_y$ in two dimensions are used for comparison. Other velocity components follow a similar trend.
We consider the wave pattern with two coexisting attractors in a spherical shell as discussed by Tilgner (Reference Tilgner1999) and Rieutord et al. (Reference Rieutord, Georgeot and Valdettaro2001). The aspect ratio and the frequency for this case are 
$\eta =0.35$ and 
$\omega =0.8102$, respectively. The forcing is imposed on the inner core. The numerical results of the 3-D libration at 
$E=10^{-9}$ and 
$E=10^{-11}$ illustrated by the amplitude of 
$v_\phi$ are shown in figure 4(a). The wave beams at the lower Ekman number are more separated. The wave pattern is consistent with the ray paths from the critical latitude on the forced inner core (see figure 4b). The ray propagating northwards from the critical latitude (in blue) converges onto the polar attractor 
$P_{0,\infty }^{(P)}\cdots P_{7,\infty }^{(P)}$, while that propagating southwards from the critical latitude (in green) converges onto the equatorial attractor 
$P_{0,\infty }^{(E)}\cdots P_{5,\infty }^{(E)}$. The corresponding 2-D results are not shown because the ray paths are identical and the wave pattern is qualitatively similar for the same aspect ratio and frequency. However, one should note that the phase shift 
$\varphi$ varies for different attractors and forcings. For the polar attractor with one vertex on the axis 
$Oz$, the phase shifts are 
${\rm \pi} /2$, 
${\rm \pi}$ and 
$0$ for the 3-D libration, 2-D symmetric and 2-D antisymmetric forcings, respectively. For the equatorial attractor, there is no phase shift for any of the forcing as this attractor does not touch the axis 
$Oz$.
Figure 4. Results corresponding to the 3-D libration of the inner core for 
$\eta =0.35$ and 
$\omega =0.8102$. (a) Numerical results of the amplitude of 
$v_\phi$ at 
$E=10^{-9}$ and 
$E=10^{-11}$. (b) Ray paths from the critical latitude on the inner core.
Two cuts crossing the two attractors are chosen in order to validate the critical-latitude asymptotic solution given by (3.27a,b). Figures 5 and 6 compare the velocity profiles between the asymptotic solutions and the numerical solutions at 
$E=10^{-11}$ on the cuts 
$S_1$ and 
$S_2$, respectively (see figure 4). The cut 
$S_1$ on the polar attractor is crossed only by the ray propagating northwards from the critical latitude (blue lines in figure 4b), while the cut 
$S_2$ on the equatorial attractor is crossed only by the ray propagating southwards (green lines in figure 4b). In figures 5 and 6, the vertical lines show the positions of the northward and southward critical rays when they cross 
$S_1$ and 
$S_2$, respectively. These critical positions correspond to different successive loops. From the rightmost critical position (
$r_1$) to the leftmost one (
$r_\infty$), the critical ray propagates from the first loop (
$n=1$) to the final loop (
$n=\infty$), and from the critical latitude to the final attractor. The critical-latitude solutions on these two cuts 
$S_1$ and 
$S_2$ are built by propagating the self-similar solutions from the critical latitude northwards and southwards respectively, by using the infinite sum of self-similar solutions (3.27a,b); see the dashed lines in figures 5 and 6, respectively. Since the amplitude decreases exponentially, the summation is conducted over a large enough number of loops in order to ensure convergence (around 150 loops in practice). The amplitudes are rescaled according to table 1, in order to make sure that the wave beams from the critical latitude possess the same amplitudes for all three forcings. Note that the radial dependence of the 3-D configuration is removed by multiplying the velocity with 
$\sqrt {r}$.
Figure 5. Comparison of velocity profiles between the critical-latitude asymptotic solutions and the numerical solutions on the cut 
$S_1$ of the polar attractor shown in figure 4 at 
$E=10^{-11}$ for three forcings: (a,b) 3-D libration (phase shift 
$\varphi ={\rm \pi} /2$); (c,d) 2-D symmetric forcing (phase shift 
$\varphi ={\rm \pi}$); (e,f) 2-D antisymmetric forcing (no phase shift). Here, (a,c,e) are the real parts, and (b,d,f) are the imaginary parts. Amplitudes are rescaled by the expressions for 
$|C_0|$ in table 1.
Figure 6. Same as for figure 5 but on the cut 
$S_2$ of the equatorial attractor. There is no phase shift for all three forcings.
Around the first critical position 
$r_1$, the wave beam from the critical latitude is within the first loop and has not experienced any contraction or expansion on the boundaries. It takes the same shape for all the three forcings, and the self-similar solution agrees with the numerical solution very well for both cuts and all three forcings, as shown by the profiles around 
$r_1$ in figures 5 and 6. After one loop, the wave beam moves on to the next critical position 
$r_2$. The amplitude decreases as expected from (3.25) for 
$\alpha _n<1$. The agreement between the self-similar solution and the numerical solution is still good; see the profiles around 
$r_2$ in figures 5 and 6. However, its shape is now dependent on the phase shift that it has experienced during the first loop. For the polar attractor on 
$S_1$, the profile around 
$r_2$ changes compared to that around 
$r_1$ for the 3-D libration and 2-D symmetric forcing since there is a non-zero phase shift for both these cases (see figures 5a–d). The profiles between 
$r_1$ and 
$r_2$ remain similar for the 2-D antisymmetric forcing since there is no phase shift (figures 5e,f). For the equatorial attractor on 
$S_2$, there is no phase shift and the profiles remain similar from 
$r_1$ to 
$r_2$ as shown in figure 6 for all three forcings. A similar behaviour can be observed from the critical position 
$r_2$ to the next position 
$r_3$. Interestingly, as shown in figure 6, the equatorial attractor profiles on 
$S_2$, which do not cross the rotation axis and therefore do not alter the phase, are almost the same for all forcing types, after rescaling the amplitudes according to table 1. This is another confirmation that these cases should be describable by a unique theory.
When the wave beam moves on to the position of the attractor (
$r_\infty$), successive critical positions become very close to each other and the profiles from different loops are not well separated. Finally, the wave beam just propagates on the attractor, and the summation of the self-similar solutions is conducted there. As shown by the profiles around the positions of the attractors 
$r_\infty$ in figures 5 and 6, the attractor with phase shift (figures 5a–d) is weaker than that without phase shift (figures 5e,f and 6). This phenomenon can be explained by the summation of the self-similar solutions on the attractor. When there is a phase shift on the attractor path, the self-similar solutions of successive loops with different phases cancel out, which makes the local solution in the vicinity of the attractor negligible after summation. Otherwise, the self-similar solutions with the same phase accumulate on the attractor, which makes the solution much stronger. For the polar attractor with phase shift (
$r_\infty$ in figures 5a–d), the critical-latitude solution (3.27a,b) is consistent with the numerical solution. However, for the polar attractor without phase shift (
$r_\infty$ in figures 5e,f), the asymptotic solution does not perform as well as for the positions far from the attractor. This deviation is much more obvious for the equatorial attractor without phase shift (
$r_\infty$ in figure 6), where the amplitudes are largely overestimated and the critical-latitude solution (3.27a,b) deviates from the numerical solutions gradually as the ray converges towards the attractor.
In order to investigate what is happening around the attractor, the velocity amplitude scalings with Ekman number of both the critical-latitude asymptotic solution (3.27a,b) and the numerical solution at the critical positions on the cuts 
$S_1$ and 
$S_2$ are shown in figure 7 for the 3-D libration. Similar behaviour can be observed for the other two forcings and are not shown here. The critical-latitude solution (3.27a,b) and the numerical solution at the first critical position 
$r_1$ follow the expected scaling 
$E^{1/12}$, which validates the Ekman number scaling of the wave beams from the critical latitude for a libration frequency different from that in HFRL22. Then the amplitude on 
$S_1$ (polar attractor with phase shift) decreases to a weaker level as the ray gets closer to the position of the attractor 
$r_\infty$, and the scaling is eventually closer to 
$E^{1/6}$. The scaling 
$E^{1/6}$ around 
$r_\infty$ is more obvious for the other cut 
$S_2$ of the equatorial attractor without phase shift. However, the corresponding prefactor is overpredicted by the critical-latitude solution.
Figure 7. Velocity amplitude scalings of (a,b) the critical-latitude solution and (c,d) the numerical solution, at the critical positions on the cuts (a,c) 
$S_1$ and (b,d) 
$S_2$ for the 3-D libration, at the Ekman numbers 
$[10^{-11}, 10^{-6}]$. Here, 
$r_1, r_2, \ldots$ and 
$r_\infty$ are the critical positions shown in figures 5 and 6.
The change of scaling of the critical-latitude solution could have been anticipated from (3.25) for its amplitude. We have seen that because 
$\alpha _n <1$, the amplitude of the self-similar beam decreases as it gets closer to the attractor. But the beam has also a finite width of order 
$E^{1/3}$, so the contributions obtained from each cycle superimpose on each other when the critical ray gets at a distance of this order from the attractor. This stops the decreasing of the amplitude after a number 
$n_s$ of cycles that can be estimated approximately by 
$(\log E^{1/3})/(\log \alpha _\infty )$, which corresponds to the number of contractions needed to go from 1 to 
$E^{1/3}$ with the contraction factor 
$\alpha _\infty$. The amplitude 
$C_{0,n_s}$ has then decreased from its initial value 
$C_0$ by a factor 
$\alpha _1^{1/4} \alpha _2^{1/4} \cdots \alpha _{n_s}^{1/4}$, which is close to 
$\alpha _{\infty }^{n_s/4}\approx E^{1/12}$. The velocity amplitude of the critical-latitude solution is therefore expected to become 
$O(E^{1/12})$ smaller close to the attractor and therefore of order 
$E^{1/6}$, as observed.
The amplitude of the streamfunction of the critical-latitude solution also decreases from 
$O(E^{5/12})$ to 
$O(E^{1/2})$ as we get close to the attractor. Then it becomes of the order of the Ekman pumping (see Appendix B). This means that the hypothesis of negligible Ekman pumping that has been used to obtain the reflection laws of the beam in § 3.1.2 breaks down. In particular, (3.13) should not be valid close to the attractor. We suspect that the discrepancies observed close to the attractor between the critical-latitude solution and the numerical solution are due to this effect.
In the next section, we develop a new asymptotic theory to describe the solution close to the attractor. This theory, which takes into account the Ekman pumping close to the attractor, is based on ideas developed originally by O05.
4. Solution close to the attractors
4.1. Asymptotic theory
In order to apply the results of O05, we first consider 2-D configurations without phase shift or reflection on the axis. In this framework, the inviscid problem can be solved using the 2-D streamfunction only, and a global solution for the streamfunction is obtained as a sum of two functions that are constant along the characteristics of the problem (e.g. Maas & Lam Reference Maas and Lam1995).
In the present work, we are looking for a solution localised near an attractor, say 
$P_{0,\infty }\cdots P_{J-1,\infty }$ as illustrated in figure 3, which can be written as a sum of local solutions 
$\psi _{j,\infty }^{(2\text {-}D)}(x_{\perp (j,\infty )})$ valid close to the segment 
$(j,\infty )$ only. In the 2-D inviscid framework, these local solutions just mean that the solution is transported from 
$P_{j,\infty }$ to 
$P_{j+1, \infty }$ along the lines 
$x_{(\perp j,\infty )} = \text {constant}$ without modification. Whereas the critical-latitude solution was reflected on boundary without modification, the solution constructed by O05 is forced directly by the boundary condition in the neighbourhood of the attractor. Close to a point 
$P_{j, \infty }$, the solution is expected to be composed of the incident solution 
$\psi _{j,\infty }^{(2\text {-}D)}(x_{\perp ( j,\infty )})$ and the reflected solution 
$\psi _{j+1,\infty }^{(2\text {-}D)}(x_{\perp (j+1,\infty )})$, and satisfies the boundary condition on the surface
\begin{equation} \psi_{j, \infty}^{(2\text{-}D)}(x_{{\perp} (j, \infty)}) + \psi_{j+1,\infty}^{(2\text{-}D)}(x_{{\perp} (j+1,\infty)}) = \psi_{j, \infty}^{(EP)} , \end{equation}
where the condition of being on the surface close to 
$P_{j,\infty }$ means that
\begin{equation} x_{{\perp} (j+1,\infty)} = \alpha _{j,\infty} x_{{\perp} (j, \infty)}. \end{equation}
In (4.1), 
$\psi _{j,\infty }^{(EP)}$ is the value of the streamfunction prescribed at 
$P_{j,\infty }$. In our case, this prescribed value is given by the Ekman pumping (EP) at the surface (or is zero if there is no Ekman pumping).
If we apply these conditions at each point 
$P_{0,\infty }, P_{1,\infty }, \ldots\!, P_{J-1,\infty}$ of the attractor, we end up, after a complete cycle, with an equation for each 
$\psi _{j,\infty }$ that can be written as
\begin{equation} \psi_{j,\infty}^{(2\text{-}D)}(\alpha_\infty x_{{\perp} (j, \infty)}) - \psi_{j,\infty}^{(2\text{-}D)}(x_{{\perp} (j, \infty)}) =\epsilon_{j,\infty} \delta, \end{equation}with
\begin{equation} \delta = \sum_{k=1}^{J/2}\left( \psi_{2k, \infty}^{(EP)} -\psi_{2k+1, \infty}^{(EP)} \right) \end{equation}and
\begin{equation} \alpha_\infty = \alpha_{1,\infty} \alpha_{2,\infty} \cdots \alpha_{J,\infty}, \end{equation}
where we have used the fact that the point 
$P_{j+J,\infty }$ corresponds to the point 
$P_{j,\infty }$ (implying that 
$\psi _{j,\infty }=\psi _{j+J,\infty }$ and 
$\psi _{j+J, \infty }^{(EP)} = \psi _{j, \infty }^{(EP)}$ ). The parameter 
$\epsilon _{j,\infty }$ is the sign in the streamfunction definition (3.6). The parameter 
$\alpha _\infty$ is the contraction factor of the attractor. Equation (4.3) is a functional constraint on the inviscid solution close to a 2-D attractor. It is identical to equation (3.17) of O05.
The general solution of this equation can be written as
\begin{equation} \psi_{j,\infty}^{(2\text{-}D)}(x_{{\perp} (j, \infty)}) = \frac{\epsilon_{j,\infty} \delta}{\ln\alpha_\infty}\ln |x_{{\perp} (j, \infty)}| +\sum_{n={-}\infty}^{+\infty}h_n^{{\pm}} |x_{{\perp} (j, \infty)}|^{{2n{\rm \pi}\mathrm{i}}/{\ln \alpha_\infty}} , \end{equation}
where the 
$\pm$ sign is for positive or negative 
$x_{\perp (j, \infty )}$. Interestingly, the dominant logarithmic part of this solution has a simple expression that depends only on the contraction factor 
$\alpha _\infty$ and the forcing term 
$\delta$. This part corresponds to a particular solution of (4.3), while the sum is a general homogeneous solution determined by the global ray mapping. Contrary to O05, we will not try to determine this homogeneous solution, since the global ray mapping for the attractors in a spherical shell is not expected to be simple. We will keep only the dominant logarithmic term to describe the inviscid solution close to the attractor. This hypothesis is not justified from an asymptotical point of view, but O05 showed that the correction associated with the homogeneous part was very small for his case.
If we keep only the particular solution, then we get a simple inviscid expression for the parallel velocity:
\begin{equation} v_{{\parallel} (j, \infty)}^{(2\text{-}D)} \sim \frac{\delta}{\ln\alpha_\infty}\,x_{{\perp} (j, \infty)}^{{-}1} . \end{equation}As explained already above, this singular behaviour can be smoothed by viscosity by introducing the self-similar solution of Moore & Saffman (Reference Moore and Saffman1969). The viscous solution that matches with the singular behaviour (4.7) is
\begin{equation} v_{{\parallel} (j, \infty)}^{(2\text{-}D)} \sim C_0^{(A)}\,H_1(x_{{\parallel} (j, \infty)}, x_{{\perp} (j, \infty)}) , \end{equation}with
\begin{equation} C_0^{(A)} = \frac{\delta}{\ln\alpha_\infty}\,E^{{-}1/3} , \end{equation}
where 
$H_m(x_\parallel,x_\perp )$ has been defined in (3.1).
In the above expression, the virtual source of the beam – that is, the position where 
$x_{\parallel (j, \infty )} =0$ – is, however, not known. This position can be obtained by using the argument developed in § 3.1.3 for the critical-latitude solution. In particular, if as above, 
$x_{\parallel (j, \infty )}$ is written as
\begin{equation} x_{{\parallel} (j, \infty)} = L_{j,\infty}^{(s)} + x_{{\parallel} (j, \infty)}', \end{equation}
with 
$x_{\parallel (j, \infty )}'=0$ at 
$P_{j,\infty }$, then the distance 
$L_{j,\infty }^{(s)}$ satisfies (3.20) and (3.22) with 
$n\rightarrow \infty$. For the first segment between 
$P_{0, \infty }$ and 
$P_{1, \infty }$, we then get, using (3.22),
\begin{equation} L_{0,\infty} ^{(s)} = (L _{0,\infty}^{(s)} + \varLambda_\infty )\alpha_\infty^3 , \end{equation}that is,
\begin{equation} L_{0,\infty} ^{(s)} =\varLambda_\infty\,\frac{\alpha_\infty^3}{1-\alpha_\infty^3}, \end{equation}
where 
$\varLambda _\infty$ and 
$\alpha _\infty$ are given by (3.23) and (3.24) with 
$n \rightarrow \infty$.
It is worth noting that the amplitude of the attractor solution does not change along the cycle. This is due to the particular value of the index of 
$m$ of the self-similar solution – that is, 
$m=1$ for the attractor solution – which guarantees that the amplitude does not change when the beam reflects on the boundary, as prescribed by the reflection law (3.11b).
For the 3-D configurations or when there is a phase shift during a cycle, the above considerations have to be modified. First note that in 3-D axisymmetric geometries, the streamfunction is not propagated identically along characteristics as it is in two dimensions. It evolves according to a propagator defined by the Riemann function, which is a Legendre function of index 
$-1/2$ for axisymmetric solutions considered here (but see § 2.3.2 of Rieutord et al. (Reference Rieutord, Georgeot and Valdettaro2001, for more details)). In other words, there is no simple global expression of the inviscid solution for the streamfunction in three dimensions. However, if the solution varies on small scales compared to the distance to the axis, then the propagation is almost as in two dimensions: in that case, an approximate local solution can be obtained far from the axis in the form
\begin{equation} \psi=\sqrt{r}\,\tilde{\psi}(x_\perp) ,\quad v_\parallel = \frac{\tilde{v}_\parallel (x_\perp)}{\sqrt{r}}, \end{equation}
where the 
$\sqrt {r}$ factor guarantees that these approximations are valid up to second-order corrections. The same analysis as above can then be performed for 
$\tilde {\psi }$ as long as we are far from the axis. This leads to a 3-D expression for the local solution near an attractor without phase shift that is obtained directly from (4.8) as
\begin{equation} v_{{\parallel} (j, \infty)} \sim \frac{C_0^{(A)}}{\sqrt{r}}\,H_1(x_{{\parallel} (j, \infty)}, x_{{\perp} (j, \infty)}) , \end{equation}
with 
$C_0^{(A)}$ given by (4.9) but with a slightly different expression for 
$\delta$:
\begin{equation} \delta = \sum_{k=1}^{J/2}\left( \frac{\psi_{2k, \infty}^{(EP)}}{\sqrt{r_{2k,\infty}}} -\frac{\psi_{2k+1, \infty}^{(EP)}}{\sqrt{r_{2k+1,\infty}}}\right) . \end{equation}
For the solution along the first segment 
$P_{0, \infty }P_{1, \infty }$, 
$x_{\parallel (j, \infty )}$ is still defined by (4.10) and (4.12).
For the attractor without reflection on the axis, it is the expressions (4.8) for 2-D configurations and (4.14) for 3-D configurations that we will use and compare to our numerical data.
We now want to consider the case of an attractor touching the axis. For the critical-latitude solution, we have seen that a phase shift could be generated as the beam reflects on the axis. A similar phenomenon is expected for the attractor solution. Let us first consider the 2-D configurations. We have seen that in that case, the presence of a phase shift depends on the symmetry of the forcing with respect to the 
$Oz$ axis. A phase shift 
$\varphi ={\rm \pi}$ is expected for a symmetric forcing, while no phase shift is present for an antisymmetric forcing. This is easily understood as changing 
$x$ into 
$-x$ changes 
$x_{\perp (0,\infty )}$ into 
$x_{\perp (1,\infty )}$. The local solutions 
$\psi _{0,\infty }^{(2\text {-}D)}(x_{\perp (0,\infty )})$ and 
$\psi _{1,\infty }^{(2\text {-}D)} (x_{\perp (1,\infty )})$, valid close to the two lines 
$x_{\perp (0,\infty )}=0$ and 
$x_{\perp (1,\infty )}=0$, respectively, should therefore satisfy the same symmetry as the forcing, that is,
\begin{equation} \psi_{1,\infty}^{(2\text{-}D)} (x_{{\perp} (1,\infty)}) = \psi_{0,\infty}^{(2\text{-}D)}(x_{{\perp} (0,\infty)}) \end{equation}for the symmetric forcing, and
\begin{equation} \psi_{1,\infty}^{(2\text{-}D)} (x_{{\perp} (1,\infty)}) ={-} \psi_{0,\infty}^{(2\text{-}D)}(x_{{\perp} (0,\infty)}) \end{equation}
for the antisymmetric forcing, when 
$x_{\perp (0,\infty )}=x_{\perp (1,\infty )}$. For the antisymmetric case, the axis is therefore as a vertical boundary without Ekman pumping (compare (4.17) to (4.1)). The same solution as for a 2-D attractor not reflecting on the axis can therefore be used. For the symmetric case, this is no longer the case. Owing to (4.16), (4.3) now becomes
\begin{equation} \psi_{j,\infty}^{(2\text{-}D)}(\alpha_\infty x_{{\perp} (j, \infty)}) + \psi_{j,\infty}^{(2\text{-}D)}(x_{{\perp} (j, \infty)}) =\epsilon_{j,\infty} \delta \end{equation}
if there is an even number of reflections on the axis. A particular solution to this equation is just 
$\psi _{j,\infty }^{(2\text {-}D)} ( x_{\perp (j, \infty )}) = \epsilon _{j,\infty }\delta /2$, so there is no longer a logarithmic singularity in the solution. Thus the inviscid expression (4.7) is not obtained, and neither is its viscous counterpart (4.8).
In three dimensions, as a phase shift of 
${\rm \pi} /2$ appears when the ray reflects on the axis, a similar phenomenon is expected if the number of reflections on the axis is not a multiple of 4. In that case, no logarithmic singularity should be present in the function 
$\tilde {\psi }_{j,\infty }$, and the analysis performed above should also break down. A weaker attractor solution is probably obtained in that case, which could explain why no significant attractor contribution was observed in the numerical solution when there is a phase shift. Finding the correct asymptotic form of the attractor solution in the presence of a phase shift is not an easy task. We leave it for future studies, probably in a simpler geometry.
4.2. Results
We now try to assess the performance of the attractor solution discussed above, and see whether the rather poor performance of the critical-latitude solution (3.27a,b) for the cases without phase shift (the polar attractor forced by the 2-D antisymmetric forcing in figures 5(e,f) and the equatorial attractors in figure 6) is improved. Since all three forcings are imposed on the inner core, only one vertex for each attractor is forced, namely 
$P^{(P)}_{0,\infty }$ for the polar attractor, and 
$P^{(E)}_{0,\infty }$ for the equatorial attractor (see figure 4b). The forcing term 
$\delta$ (see (4.4) and (4.15)) can be simplified to
\begin{equation} \delta = \left\{ \begin{array}{@{}ll} \psi^{(EP)}_{0,\infty}, & \mbox{2-D}, \\[3pt] \psi^{(EP)}_{0,\infty}/\sqrt{r_{0,\infty}}, & \mbox{3-D}. \end{array} \right. \end{equation}
The values of the Ekman pumping at the positions 
$P^{(P)}_{0,\infty }$ and 
$P^{(E)}_{0,\infty }$ correspond to the formulae of the inner core in table 2 of Appendix B. As shown in figure 8, the attractor solution (in red) performs much better than the previous critical-latitude solution. It demonstrates the necessity of including the Ekman pumping into the asymptotic solution as one gets close to the attractor, especially for the attractors without phase shift. However, it is unnecessary for the attractors with phase shift, as shown by figures 5(a–d), since there is no logarithmic singularity.
Table 2. Expressions of Ekman pumping 
$\psi ^{(EP)}$ (divided by 
$(({1}/{\lambda _+})-({1}/{\lambda _-}))\sqrt {E}$) generated by different forcings.
Figure 8. Comparison of velocity profiles between the numerical solution, the critical-latitude solution and the attractor solution at 
$E=10^{-11}$: (a,b) the cut 
$S_1$ of the polar attractor forced by the 2-D antisymmetric forcing; (c,d) the cut 
$S_2$ of the equatorial attractor forced by the 3-D libration. (The other 2-D forcings for the equatorial attractor show the same results as in (c,d).)
Note that we now face the problem of defining a transition between the critical-latitude solution valid during the first cycles and the attractor solution eventually valid close to the attractor. While we do not discuss this aspect of the problem in this paper, further studies would be required to clarify this transition.
The forcings considered up to now were imposed on the inner core, which excites the wave beams from the critical latitude on the inner core. These wave beams propagate towards the attractors and coexist with them. In order to further validate the attractor solution, it is helpful to consider a configuration where the attractor is not affected by the propagation of the wave beams from the forced critical latitude. Fortunately, such a configuration exists for the same aspect ratio and frequency but with the forcing imposed on the outer boundary. The wave pattern of the 3-D libration on the outer boundary is shown in figure 9(a), which can be compared to the ray paths in figure 9(b). Since only the critical latitude on the outer boundary is forced, the only option for the initial propagation direction is pointing into the bulk. As shown in figure 9(b), this ray (in cyan) propagates onto the polar attractor. The corresponding wave beam from the critical latitude on the outer boundary should possess 
$E^{1/5}$ width and 
$E^{1/5}$ amplitude (Roberts & Stewartson Reference Roberts and Stewartson1963; Noir, Jault & Cardin Reference Noir, Jault and Cardin2001; Lin & Noir Reference Lin and Noir2021), but it will not be our concern here. More importantly, figure 9(a) shows that the equatorial attractor is still present, although it is not connected to the ray emerging from the critical latitude on the outer boundary. It should thus be forced directly by the Ekman pumping at the positions of the attractor on the outer boundary. The attractor solution can be built for this attractor since there is no phase shift associated with it. Because the vertices 
$P_{2,\infty }^{(E)}$, 
$P_{3,\infty }^{(E)}$ and 
$P_{5,\infty }^{(E)}$ of the equatorial attractor are forced, the forcing term 
$\delta$ (see (4.4) and (4.15)) can be simplified to
\begin{equation} \delta = \left\{ \begin{array}{@{}ll} \psi^{(EP)}_{2, \infty}-\psi^{(EP)}_{3, \infty}-\psi^{(EP)}_{5, \infty}, & \mbox{2-D}, \\[3pt] \psi^{(EP)}_{2, \infty}/\sqrt{r_{2,\infty}}-\psi^{(EP)}_{3, \infty}/\sqrt{r_{3,\infty}}-\psi^{(EP)}_{5, \infty}/\sqrt{r_{5,\infty}}, & \mbox{3-D}. \end{array} \right. \end{equation}
The values of the Ekman pumping at the positions 
$P_{2, \infty }^{(E)}$, 
$P_{3, \infty }^{(E)}$ and 
$P_{5, \infty }^{(E)}$ correspond to the formulae of the outer boundary in table 2 of Appendix B. Figure 10 shows the comparison between the attractor solutions and the numerical solutions for the three forcings at three Ekman numbers. The amplitudes of the three forcings are rescaled. Good performance of the attractor solution is observed. As the Ekman number decreases, the agreement between the two solutions becomes better. The small ripples on the negative side of the similarity variable at the low Ekman numbers are wave beams from the critical latitude on the unforced inner core. They are much weaker, and the accumulation of them on the attractor remains negligible compared to the attractor beam. Figure 11 shows the Ekman number scalings of the attractor beams with beam width and velocity amplitude in 
$E^{1/3}$ and 
$E^{1/6}$, respectively, as expected.
Figure 9. Results corresponding to the 3-D libration of the outer boundary for 
$\eta =0.35$ and 
$\omega =0.8102$. (a) Numerical result of the amplitude of 
$v_\phi$ at 
$E=10^{-11}$. (b) Ray paths from the critical latitude on the outer boundary.
Figure 10. Comparison of velocity profiles between the attractor solutions and the numerical solutions on the cut 
$S_2$ of the equatorial attractor excited by the three forcings imposed on the outer boundary at three Ekman numbers.
Figure 11. Ekman number scalings of the equatorial attractor excited by the three forcings imposed on the outer boundary: (a) beam width measured by the distance of the peaks of the profiles in figure 10; (b) velocity amplitude taken at the critical position.
To summarise, we have seen that the solution close to an attractor without phase shift is described well by our asymptotic solution obtained by keeping only the logarithmic singularity contribution of the inviscid expression of the streamfunction. This has been observed for all types of forcing, in two and three dimensions, and for configurations where the attractor is connected to the critical latitude or not.
5. Conclusion
Using asymptotic analysis and numerical integration, we have studied the linear harmonic solution obtained in a rotating spherical shell by librating the inner or outer boundary for very small Ekman numbers. We have considered a shell aspect ratio and a forcing frequency such that the ray beams converge towards either a polar attractor touching the rotation axis, or an equatorial attractor not touching the rotation axis. Both 3-D axisymmetric and 2-D configurations with different types of forcing have been considered to analyse the effect of the geometric singularity on the axis (obtained in three dimensions) and the influence of a phase shift (present in the polar attractor in three dimensions, and in two dimensions with a symmetric forcing). We have focused our interest on the concentrated internal shear layers that appear along the ray emitted from the critical latitude on the inner core, and close to the attractors.
We have first shown that when the forcing is performed on the inner core, the dominant part of the solution is associated with a critical-latitude beam. We have shown that the characteristics of this beam are obtained by propagating the self-similar solution issued from the critical latitude on the inner core, as in an unbounded geometry (Le Dizès & Le Bars Reference Le Dizès and Le Bars2017) or for periodic ray paths (HFRL22). This self-similar solution has a width in 
$E^{1/3}$, a well-defined velocity amplitude in 
$E^{1/12}$, and a velocity structure corresponding to the singularity index 
$m=5/4$. As it propagates and reflects on boundaries (and possibly on the axis), its amplitude decreases down to 
$E^{1/6}$ until it reaches one of the two attractors.
We have then observed that the numerical solution departs from the asymptotic critical-latitude solution when we get close to the attractor, for some of the attractors. We have seen that the departure was present when the rays do not exhibit a phase shift along the attractor, that is, for the equatorial attractor and for the polar attractor in two dimensions with an antisymmetric forcing. We have then constructed a new asymptotic solution to describe the solution close to such an attractor, using results from O05. The main idea is based on the derivation of an inviscid functional equation for the streamfunction obtained by propagating the solution on a complete cycle on the attractor, taking into account contraction/expansion as well as Ekman pumping from the boundaries. The equation that is obtained when there is no phase shift is the equation obtained by O05. We have solved this equation by keeping only the logarithmic singular part. When smoothed by viscosity, this singular behaviour leads to a self-similar expression for the velocity with singularity index 
$m=1$ and amplitude in 
$E^{1/6}$. Contrary to the critical-latitude solution, the amplitude of this attractor solution depends on the Ekman pumping at the locations where the attractor touches the boundaries. We have shown that it describes correctly the numerical solution close to the attractor for all the attractors without phase shift.
From an asymptotic point of view, it would now be useful to obtain a solution that describes both the critical solution and the attractor solution in order to understand how the index characterising the self-similar solution changes from 
$m=5/4$ to 
$m=1$.
When the attractor exhibits a phase shift, the analysis of O05 cannot be applied completely. We have seen that a different functional equation is obtained for the streamfunction, which does not possess any logarithmically singular solution. We suspect that the amplitude of the solution is weaker in that case, which could explain why its contribution is not visible in the numerical solution close to the attractor. Obtaining an asymptotic expression describing the solution in that case still constitutes one of the important remaining issues.
Funding
J.H. acknowledges the China Scholarship Council for financial support (CSC 202008440260). The Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high-performance computing resources.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Verification of the spectral codes
The spectral codes are verified against the open-source spectral-element software Nek5000 (Fischer Reference Fischer1997; Nek5000 Version 19.0, Argonne National Laboratory, Illinois, available at https://nek5000.mcs.anl.gov). This code has already been used in the context of inertial wave propagation (Favier et al. Reference Favier, Barker, Baruteau and Ogilvie2014). Linear temporal simulations with the time-harmonic forcing are implemented by Nek5000. After a large enough number of periods, the time-harmonic steady state is reached, and the results at different instants are extracted to compare with the real and imaginary parts of the spectral results. Since it is almost impossible to reach the very low Ekman number 
$10^{-11}$ when solving the initial-value problem with Nek5000, the comparison is done at the relatively high Ekman number 
$10^{-6}$. The simulations are run for all three forcings considered in this work. In the 3-D configuration, the simulation is run in the upper right quarter of an annulus, with the axisymmetric and symmetric boundary conditions set on the two straight boundaries. In the 2-D configuration, the simulations are run in the upper half of an annulus, with symmetric boundary conditions set on the two straight boundaries. One of the curved boundaries is subject to the harmonic forcing, while the other is subject to the no-slip boundary condition. The aspect ratio and the forcing frequency are chosen to be 
$0.35$ and 
$\sqrt {2}$, respectively, so that the wave pattern is a simple periodic orbit as in HFRL22. The comparisons are shown in figure 12. The results of Nek5000 are shown on the left-hand side, while those of the spectral codes are shown on the right-hand side. They agree with each other very well.
Figure 12. Comparison between the direct numerical results given by Nek5000 and the spectral codes: (a) 3-D libration; (b) 2-D symmetric forcing; (c) 2-D antisymmetric forcing. The combination of the aspect ratio and forcing frequency is 
$(\eta,\omega )=(0.35,\sqrt {2})$.
On the other hand, the convergence of the spectral codes is verified by the spectra of the spherical harmonic (or Fourier) components and the Chebyshev coefficients, as in Rieutord & Valdettaro (Reference Rieutord and Valdettaro1997). Figure 13 shows the spectra for the 3-D libration imposed on the inner core with aspect ratio 
$0.35$ and forcing frequency 
$0.8102$ at the lowest Ekman number 
$10^{-11}$. The 2-D results are similar and omitted.
Figure 13. Spectra of (a) the spherical harmonic components and (b) the Chebyshev coefficients at 
$E=10^{-11}$, with resolution 
$(N,L)=(2500,8000)$. For each 
$l$ or 
$n$, the maximum value over the other spectral component is taken. The forcing corresponds to the 3-D libration imposed on the inner core. The aspect ratio is 
$\eta =0.35$, and the forcing frequency is 
$\omega =0.8102$.
Appendix B. Ekman pumping
The viscous forcing generates an Ekman layer adjacent to the boundary. The Ekman pumping plays a role in the generation of wave beams in the bulk. In order to derive the formula for it, it is convenient to use the streamfunction expression in spherical or polar coordinates.
B.1. Three-dimensional configuration
We first consider the libration imposed on the inner core. In the spherical coordinates 
$(\rho,\theta,\phi )$, the streamfunction 
$\psi$ and the associated variable 
$\chi$ are defined as
\begin{equation} v_\rho=\frac{1}{\rho^2\sin{\theta}}\,\frac{\partial\psi}{\partial\theta}, \quad v_\theta={-}\frac{1}{\rho\sin{\theta}}\,\frac{\partial\psi}{\partial\rho}, \quad v_\phi=\frac{\chi}{\rho\sin{\theta}}. \end{equation}The governing equations (2.3) are recast to
\begin{gather} -\mathrm{i}\omega\, D^2\psi+2\left(\cos{\theta}\,\frac{\partial\chi}{\partial\rho}-\frac{\sin{\theta}}{\rho}\,\frac{\partial\chi}{\partial\theta}\right)-E\,D^4\psi=0, \end{gather}
\begin{gather}-\mathrm{i}\omega\chi-2\left(\cos{\theta}\,\frac{\partial\psi}{\partial \rho}-\frac{\sin{\theta}}{\rho}\,\frac{\partial\psi}{\partial\theta}\right)-E\,D^2\chi=0, \end{gather}with the operator
\begin{equation} D^2=\frac{\partial^2}{\partial \rho^2}-\frac{1}{\rho^2\tan\theta}\,\frac{\partial}{\partial\theta}+\frac{1}{\rho^2}\,\frac{\partial^2}{\partial\theta^2}, \end{equation}and the boundary conditions
\begin{gather} \psi=\partial\psi/\partial\rho=0, \quad \chi=\eta^2\sin^2\theta \quad \mbox{at} \ \rho=\eta, \end{gather}
\begin{gather}\psi=\partial\psi/\partial\rho=\chi=0 \quad \mbox{at} \ \rho=1. \end{gather} The length scale of the Ekman layer is 
$\sqrt {E}$. The radial distance to the centre is rescaled as
\begin{equation} \hat{\rho}=(\rho-\eta)/\sqrt{E}. \end{equation}
The streamfunction 
$\psi$ and the associated 
$\chi$ are expanded as to the leading order:
\begin{equation} \psi=\sqrt{E}\,\hat{\psi}^{(1)}(\hat{\rho},\theta), \quad \chi=\hat{\chi}^{(0)}(\hat{\rho},\theta). \end{equation}In the leading order, the governing equations (B2) become
\begin{gather} -\mathrm{i}\omega\,\frac{\partial^2\hat{\psi}^{(1)}}{\partial\hat{\rho}^2}+2\cos{\theta}\,\frac{\partial\hat{\chi}^{(0)}}{\partial\hat{\rho}}-\frac{\partial^4\hat{\psi}^{(1)}}{\partial\hat{\rho}^4}=0, \end{gather}
\begin{gather}-\mathrm{i}\omega\hat{\chi}^{(0)}-2\cos{\theta}\,\frac{\partial\hat{\psi}^{(1)}}{\partial\hat{\rho}}-\frac{\partial^2\hat{\chi}^{(0)}}{\partial\hat{\rho}^2}=0, \end{gather}with the boundary conditions
\begin{gather} \psi=\partial\psi/\partial\rho=0, \quad \chi=\eta^2\sin^2\theta \quad \mbox{at} \ \hat{\rho}=0, \end{gather}
\begin{gather}\partial\psi/\partial\rho\rightarrow0, \quad \chi\rightarrow0 \quad \mbox{as} \ \hat{\rho}\rightarrow \infty. \end{gather}The solution of the streamfunction is obtained as
\begin{equation} \psi=\sqrt{E}\,\mathrm{i}\,\frac{\eta^2\sin^2\theta}{2}\left(-\frac{{\rm e}^{-\lambda_+\hat{\rho}}}{\lambda_+}+\frac{1}{\lambda_+}+\frac{{\rm e}^{-\lambda_-\hat{\rho}}}{\lambda_-}-\frac{1}{\lambda_-}\right), \end{equation}
with 
$\lambda _\pm$ defined as
\begin{equation} \lambda_+ = (1-i)\sqrt{\omega/2+\cos{\theta}} \end{equation}and
\begin{equation} \lambda_- = \left\{ \begin{array}{@{}ll} (1-i)\sqrt{\omega/2-\cos{\theta}}, & \omega/2>\cos{\theta}, \\[2pt] (1+i)\sqrt{\cos{\theta}-\omega/2}, & \omega/2<\cos{\theta}. \end{array} \right. \end{equation}
When 
$\hat {\rho }$ goes to 
$+\infty$, the Ekman pumping is obtained as
\begin{equation} \psi^{(EP)}=\frac{\mathrm{i}\eta^2\sin^2\theta}{2}\left(\frac{1}{\lambda_+}-\frac{1}{\lambda_-}\right)\sqrt{E}. \end{equation}
The Ekman pumping blows up at the critical colatitude 
$\theta _c=\arccos {\omega /2}$.
When the libration is imposed on the outer boundary, the boundary conditions become different, and the corresponding Ekman pumping can be obtained similarly:
\begin{equation} \psi^{(EP)}=\frac{-\mathrm{i}\sin^2\theta}{2}\left(\frac{1}{\lambda_+}-\frac{1}{\lambda_-}\right)\sqrt{E}. \end{equation}B.2. Two-dimensional configuration
In the 2-D configuration, the governing equations (2.22) of the streamfunction and the associated variable 
$\chi$ in polar coordinates are solved asymptotically using the same scaling 
$E^{1/2}$ as in the 3-D configuration. The expressions of the Ekman pumping generated by different forcings are given in table 2. Note that 
$\vartheta$ in the 2-D configuration has been replaced by 
${\rm \pi} /2-\theta$ in order to keep expressions similar to the 3-D counterpart.
Note that the Ekman pumping is 
$O(E^{1/2})$, except at the critical latitude where the Ekman pumping blows up.
