© ESO, 2016

1. Introduction

As compared to the wealth of new information obtained from the seismology of solar-type and red-giant stars (Miglio 2015), reliable seismic diagnostics on non-evolved, intermediate-mass and massive pulsators (γ Doradus, δ Scuti, β Cephei, SPB or Be stars) are scarce and, in any case, restricted to atypical slowly rotating stars (Kurtz et al. 2014; Saio et al. 2015; Triana et al. 2015). This stems from our limited understanding of the effects of rotation on oscillation modes (Lignières 2013; Townsend 2014; Reese 2015), at least for the high rotation rates of typical non-evolved, intermediate-mass and massive stars (Royer et al. 2007; Lignières 2013). The two classical approximations to treat rotational effects have been the perturbative expansion in Ω /ω, where Ω is the rotation rate and ω is the pulsation frequency (Saio 1981), and the so-called traditional approximation of the Coriolis acceleration (Eckart 1960). The perturbative expansion happens to be valid for small rotation rates but not for those of most early-type stars (Reese et al. 2006; Ballot et al. 2010, 2013). The traditional approximation greatly simplifies the understanding of Coriolis effects on low-frequency g-modes (e.g. Townsend 2003; Bouabid et al. 2013) but the range of validity and accuracy of this approximation are not known. The centrifugal deformation in particular is not taken into account in the traditional approximation.

More recently, the development of dedicated numerical codes has made it possible to compute modes and frequencies accurately without approximating the Coriolis and centrifugal effects on the modes (Reese et al. 2006, 2009; Ouazzani et al. 2012). This is a crucial step towards a confrontation with observational data. However, in addition to accurate theoretical frequencies, mode identification of observed frequencies requires a priori knowledge of the mode properties. In the case of slowly rotating stars, our knowledge of the expected frequency patterns is important for the interpretation of observed frequency spectra. It comes from the short-wavelength asymptotic theory of oscillation modes, which in the case of spherically symmetric non-rotating stars provides analytical expressions of the asymptotic frequency patterns. When applied to rapidly rotating stars, the same short-wavelength asymptotic analysis is not as straightforward. As in geometrical optics, it first leads to a ray model that describes wave propagation. Then, the conditions to produce a mode from the positive interference between rays that need to be established. Methods and concepts to infer the properties of the modes from those of the rays have been mostly developed in the field of quantum physics. The basic idea is to study rays as trajectories of a Hamiltonian system and then to relate the modes to the phase-space structures of the ray dynamics. For acoustic modes in rotating stars, such a ray-based asymptotic theory provided a physical classification of the modes as well as quantative predictions that have been successfully confronted with numerically computed high-frequency p-modes (Lignières & Georgeot 2008, 2009; Pasek et al. 2011, 2012).

In this paper, we intend to construct a similar ray-based asymptotic theory for gravito-inertial modes. Equations governing gravito-inertial rays in rotating stars are derived and then solved to obtain a global view of the properties of axisymmetric gravito-inertial rays and their evolution with stellar rotation. The detailed confrontation with numerically computed modes will be presented in a companion paper, although a first successful comparison in the particular case of so-called rosette modes has been already presented in a proceedings paper (Ballot et al. 2012).

Ray models are usually obtained from the Wentzel-Kramers-Brillouin (WKB) approximation of the equations governing small-amplitude perturbations, and this approximation is valid when the wavelength of the perturbation is much shorter than the length scale of the background spatial variations. Gravito-inertial ray models have been first developed for the ocean and the Earth’s atmosphere. One of the main motivation in this context is to model angular momentum transport and deposition by waves (Fritts & Alexander 2003; Marks & Eckermann 1995; Broutman et al. 2004). In stellar physics, such WKB ray models have been limited to pure gravity waves (Gough 1993; Alvan et al. 2015) and have not been used to study angular momentum transport or gravito-inertial modes. Instead, this is another ray theory, the method of characteristics, which has been applied by Dintrans & Rieutord (2000) to the problem of stellar gravito-inertial waves and their relation with gravito-inertial modes. As detailed in Harlander & Maas (2006), the WKB approximation and the method of characteristics differ in general, although they can be identical in some specific cases. For what concerns stellar gravito-inertial waves, they differ in that the method of characteristics cannot take into account the compressibility effects that produce the back-refraction of the waves approaching the stellar surface. We thus expect WKB ray models to be more realistic than the method of characteristics.

The paper is organised as follows. In Sect. 2, the WKB approximation is applied to the wave equation governing small perturbations about a model of a uniformly rotating star. A general eikonal equation for gravito-inertial waves is derived and the domains of propagation of these waves inside the star are discussed. In Sect. 3, the ray model for axisymmetric rays is derived from the eikonal equation, its Hamiltonian nature is emphasised and the standard tool to investigate the phase-space structure of the Hamiltonian ray dynamics is presented. A detailed numerical study of the ray dynamics for uniformly rotating polytropic models of stars is then presented in Sect. 4. The interpretation of the dynamics in terms of gravito-inertial modes is discussed in Sect. 5 and conclusions are provided in Sect. 6.

2. Eikonal equation for gravito-inertial waves

In this section, we derive an eikonal equation for short-wavelength, adiabatic gravito-inertial waves in a uniformly rotating polytropic model of star. The stellar model is specified in Sect. 2.1, the equations for small perturbations about this model are put in a convenient form in Sect. 2.2, the eikonal equation obtained by applying the WKB approximation to the perturbation equations is presented in Sect. 2.3, and the domains of propagation of gravito-inertial waves are discussed in Sect. 2.4. Only the main equations are presented in this section because their derivation are detailed in Appendix A.

2.1. Polytropic model of rotating star

The stellar model is a self-gravitating uniformly rotating monatomic gas that verifies a polytropic relation, which is assumed to give a reasonably good approximation of the relation between pressure and density in real stars (Hansen & Kawaler 1994). The governing equations can be written as

$\begin{array}{ccc}& & {\mathit{P}}_{\mathrm{0}}\mathrm{=}\mathit{K}{\mathit{\rho }}_{\mathrm{0}}^{\mathrm{1}\mathrm{+}\mathrm{1}\mathit{/}\mathit{\mu }}\mathit{,}\\ & & \mathrm{0}\mathrm{=}\mathrm{-}{\nabla }{\mathit{P}}_{\mathrm{0}}\mathrm{-}{\mathit{\rho }}_{\mathrm{0}}{\nabla }\mathrm{(}{\mathit{\psi }}_{\mathrm{0}}\mathrm{-}{\mathrm{\Omega }}^{\mathrm{2}}{\mathit{w}}^{\mathrm{2}}\mathit{/}{\mathrm{2}}^{\mathrm{)}}\mathit{,}\\ & & \mathrm{\Delta }{\mathit{\psi }}_{\mathrm{0}}\mathrm{=}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{0}}\mathit{,}\end{array}$where P₀ is the pressure, ρ₀ the density, K the polytropic constant, μ the polytropic index, ψ₀ the gravitational potential, Ω the rotation rate, w the distance to the rotation axis, and G the gravitational constant.

The uniform rotation ensures that the flow is barotropic. A pseudo-enthalpy h₀ = ^∫dP₀/ρ₀ = (1 + μ)P₀/ρ₀ can then be introduced and quantities describing the equilibrium model, such as the sound speed ${\mathit{c}}_{\mathrm{s}}\mathrm{=}\sqrt{{\mathrm{\Gamma }}_{\mathrm{1}}{\mathit{P}}_{\mathrm{0}}\mathit{/}{\mathit{\rho }}_{\mathrm{0}}}$, where Γ₁ = 5 / 3 is the first adiabatic exponent of the gas; the effective gravity ${{g}}_{\mathrm{0}}\mathrm{=}\mathrm{-}{\nabla }{}^{\mathrm{(}}\mathit{\psi }_{\mathrm{0}}\mathrm{-}{\mathrm{\Omega }}^{\mathrm{2}}{\mathit{w}}^{\mathrm{2}}\mathit{/}{\mathrm{2}}^{\mathrm{)}}$; and the Brunt-Väisälä frequency, given by the relation ${\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{=}{{g}}_{\mathrm{0}}\mathrm{·}\left(\frac{{\nabla }{\mathit{\rho }}_{\mathrm{0}}}{{\mathit{\rho }}_{\mathrm{0}}}\mathrm{-}\frac{\mathrm{1}}{{\mathrm{\Gamma }}_{\mathrm{1}}}\frac{{\nabla }{\mathit{P}}_{\mathrm{0}}}{{\mathit{P}}_{\mathrm{0}}}\right)$ are expressed in terms of h₀: $\begin{array}{ccc}{{g}}_{\mathrm{0}}\mathrm{=}{\nabla }{\mathit{h}}_{\mathrm{0}}\mathit{,} {\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\frac{{\mathrm{\Gamma }}_{\mathrm{1}}}{\mathit{\mu }\mathrm{+}\mathrm{1}}{\mathit{h}}_{\mathrm{0}}\mathit{,} {\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{=}\left(\frac{\mathit{\mu }{\mathrm{\Gamma }}_{\mathrm{1}}}{\mathit{\mu }\mathrm{+}\mathrm{1}}\mathrm{-}\mathrm{1}\right)\frac{{\mathit{g}}_{\mathrm{0}}^{\mathrm{2}}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\mathrm{·}& & \end{array}$(4)In polytropic models of stars, the density ρ₀ vanishes at the surface. Moreover, the polytropic relation Eq. (1)implies that the sound speed c_s also vanishes, and the last of Eqs. (4)shows that the Brunt-Väisälä frequency N₀ goes to infinity. These properties are not necessarily valid in more realistic stellar models.

Specifying the mass and rotation rate of the star is not sufficient to determine the polytropic model in physical units. This requires fixing an additional parameter, for example the stellar radius (suitable parameter choices are discussed in Hansen & Kawaler 1994; for the non-rotating case, and in Christensen-Dalsgaard & Thompson 1999, for the rotating case). In the following, however, we only present dimensionless quantities that do not depend on the choice of this additional parameter. The rotation rate Ω is compared to ${\mathrm{\Omega }}_{\mathrm{K}}\mathrm{=}{{}^{\mathrm{(}}\mathit{GM}\mathit{/}{{\mathit{R}}_{\mathrm{e}}^{\mathrm{3}}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}$, the limiting rotation rate for which the centrifugal acceleration equals the gravity at the equator, where M is the stellar mass and R_e the equatorial radius. Other frequency units important for our problem are (i) the inverse of the dynamical timescale associated with hydrostatic equilibrium ${\mathit{\omega }}_{\mathrm{0}}\mathrm{=}{{}^{\mathrm{(}}\mathit{GM}\mathit{/}{{\mathit{R}}_{\mathrm{p}}^{\mathrm{3}}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}$, which provides a lower bound for acoustic wave frequencies, where R_p is the polar radius; (ii) the Brunt-Väisälä frequency N₀, an upper bound for gravity wave frequencies; and (iii) twice the rotation rate f = 2Ω, the upper bound for pure inertial wave frequencies.

Fig. 1 Top: solid black line shows the profile of the Brunt-Väisälä frequency N0 normalised by ω0 for the non-rotating stellar model. The red lines show N0 for the Ω = 0.84ΩK model along the polar (solid line) and equatorial (dashed line) radii. A thick red tick on the x-axis indicates the polar radius Rp. Top: map in a meridional plane of N0 normalised by ω0 for the Ω = 0.84ΩK model. The black line is the star surface.

The spatial distribution of the Brunt-Väisälä frequency is shown in Fig. 1. In the non-rotating spherically symmetric case (see the black curve on the upper figure), N₀ vanishes at the centre, has a local maximum, N_0,max, around r = 0.28R and then shows a local minimum, N_0,min, in the stellar envelope near r = 0.74R. A simple analytical expression of N₀, which is valid for envelope models, (see Appendix E) shows that a local minimum is present at r = 3 / 4R for all polytropic indices, suggesting that such a minimum of N₀ is a generic feature of stars with radiative envelopes. The Brunt-Väisälä frequency distribution becomes anisotropic in centrifugally deformed stars. In particular, the amplitudes and positions of the local minima are clearly different along the polar and the equatorial axes. This is apparent in Fig. 1, which shows iso-contours of N₀ (bottom) as well as the polar and equatorial radial profiles of N₀ (top) for a Ω / Ω_K = 0.84 model. In the following, the local minima along the polar and equatorial axes are denoted ${\mathit{N}}_{\mathrm{0}\mathit{,}\min }^{\mathrm{p}}$ and ${\mathit{N}}_{\mathrm{0}\mathit{,}\min }^{\mathrm{e}}$, respectively. By comparison, the variation of the local maximum of the Brunt-Väisälä frequency between its values along the pole (${\mathit{N}}_{\mathrm{0}\mathit{,}\max }^{\mathrm{p}}$) and along the equatoral radius (${\mathit{N}}_{\mathrm{0}\mathit{,}\max }^{\mathrm{e}}$) remains small. This is because the deviation from sphericity is weaker in the central layers.

2.2. Perturbation equations and boundary conditions

Time-harmonic, small-amplitude perturbations of the stellar model are studied under two approximations. The first approximation is to neglect viscous and thermal dissipations. Both processes have a very small effect on the value of oscillation frequencies, at least in the frequency range of observed modes. Non-adiabatic processes play an essential role in the oscillation amplitude, in particular, through excitation mechanisms such as the κ mechanism, but determining mode amplitudes is outside the scope of the present study. The second approximation, known as the Cowling approximation, is to neglect perturbations of the gravitational potential induced by density perturbations. It is justified for oscillations of small wavelengths and, therefore, fully compatible with the small-wavelength regime considered here. Under these two approximations, the linear equations governing the evolution of small-amplitude perturbations can be written as

$\begin{array}{ccc}& & \\ & & \\ & & \end{array}$where v, ρ, and P, are respectively the Eulerian perturbations of velocity, density, and pressure, and f = 2Ω is the rotation vector.

Before applying the short-wavelength approximation, the perturbation equations for v, ρ, and P are reduced to a single equation for , the complex amplitude of the time-harmonic pressure perturbation . According to the calculations detailed in Appendix A.1, we obtain the equation $\begin{array}{ccc}\mathrm{\Delta }\mathit{P̂}\mathrm{=}\frac{{\mathit{f}}^{\mathrm{2}}}{{\mathit{\omega }}^{\mathrm{2}}}{\mathrm{\nabla }}_{\mathit{z}}^{\mathrm{2}}\mathit{P̂}\mathrm{+}\frac{{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}}{{\mathit{\omega }}^{\mathrm{2}}}{\mathrm{\Delta }}_{\mathrm{\perp }}\mathit{P̂}\mathrm{+}\mathrm{(}{{𝒱}}^{\mathrm{\prime }}\mathrm{+}\mathit{i}{{{𝒱}}^{\mathrm{\prime }}}_{\mathrm{m}}\mathrm{)}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{+}{\mathit{M}}^{\mathrm{\prime }}\mathit{P̂}\mathit{,}& & \end{array}$(8)where the operators ${\mathrm{\nabla }}_{\mathit{z}}^{\mathrm{2}}$ and Δ_⊥ are related to the unit vectors parallel to the rotation axis, e_z = f/f, and the effective gravity, e_∥ = −g₀/g₀, by $\begin{array}{ccc}& & {\mathrm{\nabla }}_{\mathit{z}}^{\mathrm{2}}\mathrm{\equiv }{\mathrm{\nabla }}_{\mathit{z}}\mathrm{\left(}{\mathrm{\nabla }}_{\mathit{z}}\mathrm{\right)}\mathrm{with}{\mathrm{\nabla }}_{\mathit{z}}\mathrm{\equiv }{{e}}_{\mathit{z}}\mathrm{·}{\nabla }\mathit{,}\\ & & {\mathrm{\Delta }}_{\mathrm{\perp }}\mathrm{\equiv }{\nabla }\mathrm{·}\mathrm{\left(}{{\nabla }}_{\mathrm{\perp }}\mathrm{\right)}\mathrm{with}{{\nabla }}_{\mathrm{\perp }}\mathrm{\equiv }{\nabla }\mathrm{-}{{e}}_{\mathrm{\parallel }}{\mathrm{\nabla }}_{\mathrm{\parallel }}\mathrm{and}{\mathrm{\nabla }}_{\mathrm{\parallel }}\mathrm{\equiv }{{e}}_{\mathrm{\parallel }}\mathrm{·}{\nabla }\mathit{.}\end{array}$The expressions of , and M′ are given by Eqs. ((A.39)–(A.41)), respectively. As a linear combination of f and g₀, the vector is contained in a meridional plane, whereas is colinear to f ∧ g₀ and thus perpendicular to this plane. First order derivatives in the meridional plane can be suppressed from this equation by introducing a variable and by choosing the function a adequately (see Appendix A.2). We then obtain a wave equation for acoustic and gravito-inertial waves in a uniformly rotating, centrifugally deformed, barotropic star as follows: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\mathrm{\Delta }\mathrm{\Psi ̂}\mathrm{=}\hspace{0.17em}\mathrm{-}\frac{{\mathit{\omega }}^{\mathrm{2}}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\mathrm{\Psi ̂}\mathrm{+}\frac{{\mathit{f}}^{\mathrm{2}}}{{\mathit{\omega }}^{\mathrm{2}}}{\mathrm{\nabla }}_{\mathit{z}}^{\mathrm{2}}\mathrm{\Psi ̂}\mathrm{+}\frac{{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}}{{\mathit{\omega }}^{\mathrm{2}}}{\mathrm{\Delta }}_{\mathrm{\perp }}\mathrm{\Psi ̂}\mathrm{+}\mathit{i}\frac{\mathrm{\parallel }{f}\mathrm{\wedge }{{g}}_{\mathrm{0}}\mathrm{\parallel }}{\mathit{\omega }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\mathrm{𝒯}{{e}}_{\mathit{\phi }}\mathrm{·}{\nabla }\mathrm{\Psi ̂}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{+}{\mathrm{𝒞}}^{\mathrm{\prime }}\mathrm{\Psi ̂}\mathit{,}\end{array}& & \end{array}$(11)where e_φ is the unit vector in the azimuthal direction and is given by Eq. (A.53). The first term on the right-hand side (RHS) of Eq. (11) is associated with sound waves whereas the next two terms on the RHS produce gravito-inertial waves. The fourth term on the RHS is non-zero only for non-axisymmetric perturbations in rotating stars. It contains the latitudinal derivative of the Coriolis parameter fcosθ, which gives rise to Rossby and Kelvin waves. Such a simple form for the wave equation comes at the cost of very complicated expressions for a and as functions of c_s,N₀,g₀,f, and ω (see Eqs. (A.48) and (A.54)).

2.3. The short-wavelength approximation of the perturbation equations

The ray model results from the WKB approximation of the wave Eq. (11). It consists in looking for wave-like solutions of the form $\mathrm{\Psi }\mathrm{=}\mathrm{\Re }\mathrm{\left\{}\mathit{A}\mathrm{\right(}{x}\mathrm{\left)}\exp \mathrm{i}\mathrm{\right[}\mathrm{\Phi }\mathrm{\left(}{x}\mathrm{\right)}\mathrm{-}\mathit{\omega t}\mathrm{\left]}\mathrm{\right\}}\mathrm{=}\mathrm{\Re }\mathrm{\left\{}\mathrm{\Psi ̂}\mathrm{\right(}{x}\mathrm{\left)}\exp \mathrm{\right(}\mathrm{-}\mathrm{i}\mathit{\omega t}\mathrm{\left)}\mathrm{\right\}}\mathit{,}$(12)assuming that the associated wavelength λ_w ~ ∥ ∇Φ ∥ ^-1 is much shorter than the typical length scale on which the background medium varies, hereafter denoted L_b.

The phase Φ and the amplitude A of are expanded into power series of 1 / Λ, where Λ = L_b/λ_w is assumed to be large, yielding$\mathrm{\Phi }\mathrm{=}\mathrm{\Lambda }\left({\mathrm{\Phi }}_{\mathrm{0}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{\Lambda }}{\mathrm{\Phi }}_{\mathrm{1}}\mathrm{+}\mathit{...}\right) \mathrm{and} \mathit{A}\mathrm{=}{\mathit{A}}_{\mathrm{0}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{\Lambda }}{\mathit{A}}_{\mathrm{1}}\mathrm{+}\mathit{...}$(13)Introducing these expansions into the wave Eq. (11) and considering only the dominant terms, we obtain an equation for Φ₀, the so-called eikonal equation. The amplitude A₀ is determined at the next order. The eikonal equation obtained with this procedure depends on the frequency range considered. If ω is of the order of Λ, the dominant terms of the wave Eq. (11) are the left-hand side (LHS) and the first term on the RHS, leading to an eikonal equation for sound waves ω = c_sk, where k = ∇Φ₀ and k = ∥ k ∥. If , the dominant terms are those that contain second derivatives of , leading to an eikonal equation for gravito-inertial waves. The fourth term on the RHS of the wave Eq. (11)is a first-order azimuthal derivative of and is therefore negligible when . Thus, Rossby waves are not included in the eikonal equation if one considers the regime or if one restricts to axisymmetric perturbations. Before writing down the eikonal equation for gravito-inertial waves, the role of the so-called constant term must be clarified. Indeed, if a wave, acoustic or gravito-inertial, of vanishingly small wavelength approaches the stellar surface, it will cross the surface without being refracted back into the star. This regime corresponds to the limit λ_w → 0, where the constant term is negligible and the resulting eikonal equation contains no surface refraction effect.

We are interested, however, in waves that are refracted back because oscillations modes are the result of wave interferences within the star. As known from previous studies on acoustic and gravity waves (Aerts et al. 2010), this back-refraction occurs for outwards-travelling waves that encounter atmospheric layers whose pressure scale heights are smaller that the wavelength of the wave. We also know that the pressure scale height strongly decreases as one approaches the surface in such a way that if L_b is of the order of the pressure scale height in the interior and H_s is the surface pressure scale height, we obtain H_s/L_b ≪ 1. Thus, there exists a wavelength range for which we have at the same time λ_w/L_b = 1 / Λ ≪ 1 inside the star and λ_w of the order of or larger than H_s. In the wave equation, this surface refraction effect is accounted for by the constant term because is inversely proportional to the square of the pressure scale height (see Eq. (16) below). Thus, an eikonal equation containing this term is simply obtained by assuming in the expansion given by Eq. (13). This yields $\begin{array}{ccc}{\mathit{k}}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{f}}^{\mathrm{2}}}{{\mathit{\omega }}^{\mathrm{2}}}{\mathit{k}}_{\mathit{z}}^{\mathrm{2}}\mathrm{+}\frac{{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}}{{\mathit{\omega }}^{\mathrm{2}}}\mathrm{(}{\mathit{k}}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{+}{{\mathit{k}}_{\mathit{\phi }}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{-}\mathrm{𝒞} \mathrm{with}\mathrm{𝒞}\mathrm{=}{\mathrm{𝒞}}^{\mathrm{\prime }}\mathrm{-}\frac{{\mathit{\omega }}^{\mathrm{2}}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\mathit{,}& & \end{array}$(14)where k_z = e_z·k, k_φ = e_φ·k and k_⊥ = e_⊥·k, where the unit vector e_⊥ is defined as e_⊥ = e_φ ∧ e_∥. In Eq. (14), the term in ${\mathit{k}}_{\mathit{\phi }}^{\mathrm{2}}$ results from the fact that ∇_⊥ is not just the derivative along e_⊥. It also contains derivatives in φ, and so does Δ_⊥.

Although the full expression of is a very complex function of the equilibrium quantities c_s,N₀,g₀,f, and of ω (see Eq. (A.54)), we only need to estimate where it becomes of the same order of magnitude as the other terms in the dispersion relation, that is near the surface. The increase of towards the surface is due to the sharp increase of N₀ and 1 /c_s, while g₀ does not vary much in the stellar envelope. In a polytropic model of star, we also have N₀ ∝ 1 /c_s near the surface. In order to only retain the dominant term, we look for an expansion of in powers of c_s and find that can be written as the ratio between two polynomials of ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$, namely $\begin{array}{ccc}\mathrm{𝒞}\mathrm{=}\frac{{\mathit{C}}_{\mathrm{0}}\mathrm{+}{\mathit{C}}_{\mathrm{1}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{+}\mathrm{𝒪}\mathrm{\left(}{\mathit{c}}_{\mathrm{s}}^{\mathrm{4}}\mathrm{\right)}}{{\mathrm{(}\mathit{\alpha K}\mathrm{+}\mathit{G}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{4}}}\mathit{,}& & \end{array}$(15)where the terms $\mathit{\alpha }\mathrm{=}\frac{{\mathrm{\Gamma }}_{\mathrm{1}}\mathit{\mu }}{\mathit{\mu }\mathrm{+}\mathrm{1}}\mathrm{-}\mathrm{1}$, K = (f·g₀)²−ω²g₀², G = ω²(ω²−f²) and the terms C₀ and C₁ are at the surface. In the super-inertial¹ regime ω>f, K does not vanish. The dominant term of is thus $\frac{{\mathit{C}}_{\mathrm{0}}}{\mathrm{(}\mathit{\alpha K}{\mathrm{)}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{4}}}$ and we obtain $\begin{array}{ccc}\mathrm{𝒞}\mathrm{=}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{f}}^{\mathrm{2}}}{{\mathit{\omega }}^{\mathrm{2}}}{\cos }^{\mathrm{2}}\mathrm{\Theta }\right)\frac{{\mathrm{\Gamma }}_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}}\frac{\mathit{\mu }\mathrm{-}\mathrm{1}}{\mathit{\mu }\mathrm{+}\mathrm{1}}\frac{{\mathit{g}}_{\mathrm{0}}^{\mathrm{2}}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{4}}}\mathrm{+}\mathrm{𝒪}\mathrm{(}\mathrm{1}\mathit{/}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{)}\mathit{,}& & \end{array}$(16)where cosΘ = −f.g₀/fg₀. In the sub-inertial regime ω<f, both K and C₀ vanish at the critical angles Θ_c and π−Θ_c such that cos²Θ_c = ω²/f². In this regime and in the regions close to the stellar surface and the critical angles, is given by $\begin{array}{ccc}\mathrm{𝒞}\mathrm{=}\frac{{\mathit{C}}_{\mathrm{1}}\mathrm{\left(}{\mathrm{\Theta }}_{\mathrm{c}}\mathrm{\right)}}{{\mathit{G}}^{\mathrm{2}}}\frac{\mathrm{1}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{6}}\mathrm{\left(}{\mathrm{\Theta }}_{\mathrm{c}}\mathrm{\right)}}\mathrm{+}\mathrm{𝒪}\mathrm{(}\mathrm{1}\mathit{/}{\mathit{c}}_{\mathrm{s}}^{\mathrm{4}}\mathrm{)}\mathit{.}& & \end{array}$(17)The detail of the asymptotic expression of is presented in Appendices A.3 and A.4. In the following, we always use the expression (16) for . As we see in the next sections, inspection of the eikonal and ray equations near the critical angle indicates that rays tend to avoid the surface layers close to the critical angle. This is further confirmed by our computations that do not show rays reaching the surface close to the critical angle.

Using Eq. (16), the eikonal equation can be written as $\begin{array}{ccc}{\mathit{\omega }}^{\mathrm{2}}\mathrm{=}{\mathit{f}}^{\mathrm{2}}\frac{{\mathit{k}}_{\mathit{z}}^{\mathrm{2}}}{{\mathit{k}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{c}}^{\mathrm{2}}}\mathrm{+}{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}\frac{{\mathit{k}}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathit{\phi }}^{\mathrm{2}}}{{\mathit{k}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{c}}^{\mathrm{2}}}\mathrm{+}{\mathit{f}}^{\mathrm{2}}{\cos }^{\mathrm{2}}\mathrm{\Theta }\frac{{\mathit{k}}_{\mathrm{c}}^{\mathrm{2}}}{{\mathit{k}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{c}}^{\mathrm{2}}}\mathit{,}& & \end{array}$(18)with $\begin{array}{ccc}{\mathit{k}}_{\mathrm{c}}^{\mathrm{2}}\mathrm{=}\frac{{\mathrm{\Gamma }}_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}}\frac{\mathit{\mu }\mathrm{-}\mathrm{1}}{\mathit{\mu }\mathrm{+}\mathrm{1}}\frac{{\mathit{g}}_{\mathrm{0}}^{\mathrm{2}}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{4}}}\mathrm{·}& & \end{array}$(19)This equation is similar to the local dispersion relations generally used to study gravito-inertial waves in the Earth’s atmosphere (Fritts & Alexander 2003) and to construct associated ray models (Marks & Eckermann 1995). The present relation is nevertheless more general in that it takes the centrifugal deformation of the star into account and does not make the traditional approximation. In the following we shall restrict our study to axisymmetric perturbations, that is to the case k_φ = 0.

2.4. Domains of propagation

The eikonal Eq. (18) allows us to discuss the domain of propagation of gravito-inertial waves in a more general context than previous works. In particular, when compared to the work of Dintrans & Rieutord (2000), the following discussion includes the effects of the centrifugal deformation and, through the k_c term, those of the background compressibility.

According to Eq. (12), non-evanescent wave solutions are possible only when both components of the wave vector k are real numbers. The region of space where the eikonal equation admits real solutions thus determines the domain where rays can propagate. After projecting k onto an orthogonal basis, the eikonal equation takes the general form of a conic. Using (e_∥,e_⊥) as a basis and the relations ${\mathit{k}}^{\mathrm{2}}\mathrm{=}{\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{\perp }}^{\mathrm{2}}$ and k_z = k_∥cosΘ−k_⊥sinΘ, the eikonal Eq. (18) is rewritten as $\begin{array}{ccc}\mathrm{(}{{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}}^{}& \mathrm{+}& {\mathit{f}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathrm{\Theta }\mathrm{-}{\mathit{\omega }}^{\mathrm{2}}\mathrm{)}{\mathit{k}}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathit{f}}^{\mathrm{2}}\cos \mathrm{\Theta }\sin \mathrm{\Theta }{\mathit{k}}_{\mathrm{\parallel }}{\mathit{k}}_{\mathrm{\perp }}\\ & \mathrm{-}& \mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{f}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathrm{\Theta }}^{\mathrm{)}}{\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{-}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{f}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathrm{\Theta }}^{\mathrm{)}}{\mathit{k}}_{\mathrm{c}}^{\mathrm{2}}\mathrm{=}\mathrm{0.}\end{array}$(20)To discuss the constraints on the domain of propagation that can be derived from this relation, it is useful to first consider the non-rotating case. The eikonal equation then becomes $\begin{array}{ccc}\mathrm{[}{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{-}{\mathit{\omega }}^{\mathrm{2}}\mathrm{]}\frac{{\mathit{L}}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{-}{\mathit{\omega }}^{\mathrm{2}}{\mathit{k}}_{\mathit{r}}^{\mathrm{2}}\mathrm{=}{\mathit{\omega }}^{\mathrm{2}}{\mathit{k}}_{\mathrm{c}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{,}& & \end{array}$(21)with k_⊥ = k_θ = ± L/r and L is the norm of the vector x ∧ k, which is equivalent to the angular momentum. This equation corresponds to the high-sound-speed limit of eikonal equations already derived in the non-rotating case, for example by Gough (1993). The condition that ${\mathit{k}}_{\mathit{r}}^{\mathrm{2}}$ is positive translates into ${\mathit{\omega }}^{\mathrm{2}}\mathit{<}{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{r}}^{\mathrm{2}}{\mathit{k}}_{\mathrm{c}}^{\mathrm{2}}\mathit{/}{\mathit{L}}^{\mathrm{2}}\mathrm{)}$ and this inequality fully specifies the spatial limit of the resonant cavity for a wave of given ω and L. In internal regions where ${\mathit{r}}^{\mathrm{2}}{\mathit{k}}_{\mathrm{c}}^{\mathrm{2}}\mathit{/}{\mathit{L}}^{\mathrm{2}}\mathrm{\ll }\mathrm{1}$, it further simplifies into ω<N₀, a condition that no longer depends on L or k_θ. Closer to the surface, ω ≪ N₀ and the condition for propagation is ω<c₁S_L, where S_L = Lc_s/r is the Lamb frequency and c₁ a constant that can be expressed in terms of Γ₁ and μ. If this last condition is fullfilled at the surface, the gravity ray escapes the star. Otherwise, the ray is refracted back into the star and the relation ω = c₁S_L determines the outer radius of the domain of propagation. An important point is that the condition for the back-refraction in the outer layers of the star depends on L or equivalently on k_θ.

Coming back to the rotating case, the condition that k_∥ is real can be expressed as $\begin{array}{ccc}\mathit{\delta }\mathrm{=}\mathrm{\Gamma }{\mathit{k}}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{-}{\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{f}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathrm{\Theta }}^{\mathrm{)}}}^{\mathrm{2}}{\mathit{k}}_{\mathrm{c}}^{\mathrm{2}}\mathrm{\ge }\mathrm{0}\mathit{,}& & \end{array}$(22)where δ is the reduced discriminant of the eikonal equation viewed as a quadratic equation for the variable k_∥, and $\begin{array}{ccc}\mathrm{\Gamma }\mathrm{=}{\mathit{\omega }}^{\mathrm{2}}\mathrm{(}{\mathit{f}}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}^{\mathrm{2}}\mathrm{)}\mathrm{+}{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{f}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathrm{\Theta }}^{\mathrm{)}}\mathit{.}& & \end{array}$(23)The sign of Γ determines whether gravito-inertial waves can propagate in the inner region where k_c can be neglected. In contrast, the outer limit of the domain of propagation is determined by the condition δ ≥ 0 and, as expected from the condition of propagation in the non-rotating case, it depends on the value of ${\mathit{k}}_{\mathrm{\perp }}^{\mathrm{2}}$. However, contrary to the non-rotating case, ${\mathit{k}}_{\mathrm{\perp }}^{\mathrm{2}}$ is not known a priori because rotation breaks the spherical symmetry and the associated conservation of the norm of the angular momentum. The consequence is that in the general case of gravito-inertial waves in a rotating star, we are not able to predict the outer limit of the resonant cavity from the eikonal equation alone. It is thus necessary to solve the ray dynamics to find out whether a ray of a given frequency remains confined within the star and, if this is the case, the shape of the outer limit of its resonant cavity. In the ray dynamics calculation presented in this paper, all rays stay inside the stars because, at the surface of polytropic models, k_c is infinite and δ is thus necessarily negative. In real stars, k_c remains finite at the surface, and this implies that rays with sufficently high values of k_⊥ near the surface can propagate outside the stars. Waves associated with such rays are expected to have very high spatial frequencies at the surface, and are thus unlikely to be visible. However, they can play a role in the transport of angular momentum out of the stars.

In the following of the section, we thus only discuss the domain of propagation in the internal region where the ${\mathit{k}}_{\mathrm{c}}^{\mathrm{2}}$ term can be neglected, since this approximation is valid from the centre up to a certain radius that depends on the given ray. The condition of propagation Γ ≥ 0 is equivalent to $\begin{array}{ccc}{\mathit{\omega }}_{\mathrm{-}}^{\mathrm{2}}\mathit{<}{\mathit{\omega }}^{\mathrm{2}}\mathit{<}{\mathit{\omega }}_{\mathrm{+}}^{\mathrm{2}}\mathit{,}& & \end{array}$(24)with $\begin{array}{ccc}{\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{f}}^{\mathrm{2}}\mathrm{+}{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\pm }\sqrt{{\mathrm{(}{\mathit{f}}^{\mathrm{2}}\mathrm{+}{{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{2}}\mathrm{-}\mathrm{4}{\mathit{f}}^{\mathrm{2}}{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}{\cos }^{\mathrm{2}}\mathrm{\Theta }}}{\mathrm{2}}\mathrm{·}& & \end{array}$(25)If the centrifugal deformation is neglected, that is if we replace Θ by the colatitude θ, this condition of propagation is equivalent to that found by Dintrans & Rieutord (2000). Two limit cases, f ≪ N₀ and N₀ ≪ f, are relevant since the former occurs in the bulk of typical moderately rotating stellar radiative zones, whereas the latter is verified in convective zones or close to the centre of the star where N₀ vanishes. If f ≪ N₀, the conditions (24) simplify to $\begin{array}{ccc}{\mathit{f}}^{\mathrm{2}}{\cos }^{\mathrm{2}}\mathrm{\Theta }\mathit{<}{\mathit{\omega }}^{\mathrm{2}}\mathit{<}{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{+}{\mathit{f}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathrm{\Theta }\mathit{.}& & \end{array}$(26)On the other hand, if N₀ ≪ f, the conditions (24) become $\begin{array}{ccc}{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}{\cos }^{\mathrm{2}}\mathrm{\Theta }\mathit{<}{\mathit{\omega }}^{\mathrm{2}}\mathit{<}{\mathit{f}}^{\mathrm{2}}\mathrm{+}{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathrm{\Theta }\mathit{,}& & \end{array}$(27)where in the limit of vanishing N₀, the ω<f rule for pure inertial waves is recovered.

We now consider the particular case of uniformly rotating μ = 3 polytropic models. Close to the centre of the star, the approximate conditions (27) hold because N₀ vanishes there. At some distance from the centre, the Brunt-Väisälä frequency always becomes larger than f. Moreover, for moderate rotators, the approximate condition (26) is valid from a small radius up to the surface. Indeed, as long as Ω ≤ 0.38Ω_K and r> 0.1R_e, the ratio ${\mathit{f}}^{\mathrm{2}}\mathit{/}{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}$ is always smaller than 0.10.

Fig. 2 Typical domains of propagation (grey areas) given by the Γ(ω) ≥ 0 condition as a function of the gravito-inertial wave frequency ω in a Ω = 0.38ΩK polytropic model of star. In general, the outer limit of the domain, here loosely represented by the stellar surface, cannot be predicted from ω alone and is obtained by solving the ray dynamics. Panel a) corresponds to a super-inertial frequency (ω>f) and panel b) to a sub-inertial frequency (ω<f). The two bottom figures show domains with an intermediate forbidden zone when frequencies are between the local minimum and maximum of the Brunt-Väisälä frequency. Panel c) corresponds to frequencies in the ${\mathit{N}}_{\mathrm{0}\mathit{,}\min }^{\mathrm{p}}\mathit{<}\mathit{\omega }\mathit{<}{\mathit{N}}_{\mathrm{0}\mathit{,}\max }^{\mathrm{e}}$ interval and also exists in non-rotating stars, whereas panel d) corresponds to frequencies in the ${\mathit{N}}_{\mathrm{0}\mathit{,}\min }^{\mathrm{e}}\mathit{<}\mathit{\omega }\mathit{<}{\mathit{N}}_{\mathrm{0}\mathit{,}\min }^{\mathrm{p}}$ interval and is due to the centrifugal deformation of rotating stars.

Two typical domains of propagation are shown in panels (a) and (b) of Fig. 2 in the super-inertial and sub-inertial frequency ranges, respectively. The rotation of the stellar model is Ω = 0.38Ω_K. For the super-inertial case, the condition of propagation reduces to ${\mathit{\omega }}^{\mathrm{2}}\mathit{<}{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{+}{\mathit{f}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathrm{\Theta }$, which leads to the domain shown in panel (a), including the weakly prolate avoided central region. By contrast, in the sub-inertial regime the centre is never forbidden because conditions (27) are automatically fulfilled for ω<f. Away from this central region, the conditions of propagation (26) reduce to f²cos²Θ <ω² for these waves. The resulting domain of propagation is shown in panel (b). A noticeable effect of the centrifugal deformation is that the limit of the domain is no longer strictly conical away from the central region. As the centrifugal deformation of the equipotentials becomes significant towards the surface, the boundary of the domain of propagation is bent towards lower latitudes to keep Θ_c = arccos(ω/f), the angle between the direction of the effective gravity and the rotation axis, constant.

The local minimum of the Brunt-Väisälä frequency present in the stellar envelope (see Fig. 1) induces an additional forbidden region for frequencies between ${\mathit{N}}_{\mathrm{0}\mathit{,}\min }^{\mathrm{e}}$ and ${\mathit{N}}_{\mathrm{0}\mathit{,}\max }^{\mathrm{p}}$. As illustrated on panel (c) for the Ω = 0.38Ω_K model, the domains of progagation below and above the position of ${\mathit{N}}_{\mathrm{0}\mathit{,}\min }^{\mathrm{e}}$ are disjoint when N_0,min<ω<N_0,max in non-rotating stars or when ${\mathit{N}}_{\mathrm{0}\mathit{,}\min }^{\mathrm{p}}\mathit{<}\mathit{\omega }\mathit{<}{\mathit{N}}_{\mathrm{0}\mathit{,}\max }^{\mathrm{e}}$ in rotating stars. In centrifugally deformed stars, waves with frequencies in the range ${\mathit{N}}_{\mathrm{0}\mathit{,}\min }^{\mathrm{e}}\mathit{<}\mathit{\omega }\mathit{<}{\mathit{N}}_{\mathrm{0}\mathit{,}\min }^{\mathrm{p}}$ also show a forbidden region, but this time it is limited to equatorial regions so that the inner and outer domains communicate through polar regions. An example is shown in panel (d) for the Ω = 0.38Ω_K model also. We observe that the outer limit of the forbidden region is oblate whereas the inner limit is prolate. As can be inferred from the distribution of N₀ shown in Fig. 1, this property is related to the oblate and prolate shapes of the iso-contours of the Brunt-Väisälä frequency taken in the ${}^{\mathrm{[}}\mathit{N}_{\mathrm{0}\mathit{,}\min }^{\mathrm{e}}\mathit{,}{{\mathit{N}}_{\mathrm{0}\mathit{,}\min }^{\mathrm{p}}}^{\mathrm{]}}$ interval. These peculiar resonant cavities appear in a frequency range that increases with rotation as the ratio ${\mathit{N}}_{\mathrm{0}\mathit{,}\min }^{\mathrm{p}}\mathit{/}{\mathit{N}}_{\mathrm{0}\mathit{,}\min }^{\mathrm{e}}$ goes from 1 in the non-rotating case to 1.70 at Ω = 0.84Ω_K.

3. The ray model

In this section, the ray model for progressive gravito-inertial waves is derived from the eikonal equation. Its Hamiltonian structure is emphasised in Sect. 3.1 as it enables us to investigate ray paths using tools and concepts developed for Hamiltonian dynamical systems. The geometrical structure of the phase space in particular is the key to understand the nature of the dynamics. Despite its high dimensionality, the phase space can be explored numerically through cuts known as Poincaré surfaces of section (hereafter PSS). The numerical method and PSS are presented in Sects. 3.2 and 3.3, respectively.

3.1. The Hamiltonian equations describing the ray path and the wave vector evolution along the path

The eikonal Eq. (18) is a partial differential equation (PDE) for the phase function Φ₀(x). Instead of solving this equation as a PDE, the gravito-inertial ray model consists in searching for solutions along some path x(t) called the ray path. One then has to solve the coupled differential equations that govern the ray path and the evolution of k along the path.

It has been shown that this problem can be written under very general circumstances in a Hamiltonian form (Lighthill 1978; Ott 1993). In the present case, the Hamiltonian is obtained by writing the eikonal equation in the form ω = H(x,k) and the ray path is defined by the group velocity. For a given coordinate system [ x_i ], the Hamiltonian equations of the ray dynamics are (Lignières 2011) $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{x}}_{\mathit{i}}}{\mathrm{d}\mathit{t}}& \mathrm{=}& \frac{\mathit{\partial H}}{\mathit{\partial }{\mathit{k}}_{\mathit{i}}}\mathit{,}\\ \frac{\mathrm{d}{\mathit{k}}_{\mathit{i}}}{\mathrm{d}\mathit{t}}& \mathrm{=}& \end{array}$where the ${\mathit{k}}_{\mathit{i}}\mathrm{=}\frac{\mathit{\partial }{\mathrm{\Phi }}_{\mathrm{0}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}$ are the covariant components of k on the natural basis ${{e}}_{\mathit{i}}\mathrm{=}\frac{\mathit{\partial }{x}}{\mathit{\partial }{\mathit{x}}^{\mathit{i}}}$ associated with the coordinate system [ x_i ]. In Appendix B, these equations are written with spherical coordinates and the k components on the usual unit vectors e_r and e_θ.

A classical property of gravito-inertial waves is the orthogonality of the group velocity v_g = dx/ dt and the phase velocity v_p = ωk/k². Here, this property is valid in the interior but broken when k_c is not negligible. This is apparent in the following formula: ${{v}}_{\mathrm{g}}\mathrm{·}{{v}}_{\mathrm{p}}\mathrm{=}\frac{\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{f}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathrm{\Theta }}^{\mathrm{)}}{\mathit{k}}_{\mathrm{c}}^{\mathrm{2}}}{{\mathit{k}}^{\mathrm{2}}\mathrm{(}{\mathit{k}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{c}}^{\mathrm{2}}\mathrm{)}}\mathit{,}$(30)whose derivation is detailed in Appendix C.

Another important characteristic is the behaviour of the rays whenever they reach the limits of the domain of propagation. This can be discussed from the relations $\begin{array}{ccc}\frac{\mathrm{d}{x}}{\mathrm{d}\mathit{t}}\hspace{0.17em}\mathrm{·}{{e}}_{\mathrm{\parallel }}& \mathrm{=}& \mathrm{\pm }\frac{\sqrt{\mathit{\delta }}}{\mathit{\omega }\mathrm{(}{\mathit{k}}^{\mathrm{2}}\mathrm{+}{{\mathit{k}}_{\mathrm{c}}^{\mathrm{2}}}^{\mathrm{)}}}\mathit{,}\\ \frac{\mathrm{d}{x}}{\mathrm{d}\mathit{t}}\mathrm{·}{{e}}_{\mathrm{\perp }}& \mathrm{=}& \end{array}$also derived in Appendix C. Equation (31) shows that the component of the group velocity parallel to e_∥ vanishes at the limits of the domain of propagation, that is when δ = 0. If this occurs when k_c is negligible, Eq. (22) finds that Γ also vanishes there. Thus, according to Eq. (32), both components of the group velocity vanish at the limit of the domain of propagation and the turning point is an ordinary cusp point. However, if δ = 0 is reached closer to the stellar surface, where k_c cannot be neglected, Eqs. (22) and (32) imply that Γ and dx/ dt·e_⊥ do not vanish. The trajectory is then locally parallel to e_⊥, that is tangent to an equipotential, before it is redirected towards the stellar interior. Figure 3 illustrates the two types of turning points for a given trajectory in the super-inertial regime. The capacity of the numerical method to accurately compute the trajectories near cusp points has been successfully tested against an analytical solution (see Appendix D).

Fig. 3 Gravito-inertial ray showing the two types of turning points, a cusp point on the Γ = 0 boundary (in magenta) and a smooth back-refraction in the surface layers where the decreasing density scale height becomes of the order of the wavelength in that direction.

3.2. Numerical method for the ray dynamics

The gravito-inertial ray dynamics has been investigated by integrating numerically Eqs. (B.1)–(B.4) using a 5th-order Runge-Kutta method. The step size of the integration is chosen automatically to keep the local error estimate smaller than a given tolerance. The background polytropic model of star is solved numerically with an iterative scheme, which is described in Rieutord et al. (2005).

We checked that the PSS shown in Sect. 4 are not significantly modified by decreasing the maximum allowed local error of the Runge-Kutta scheme. We also checked the influence of the resolution of the background polytropic stellar model. Finally, the conservation of the Hamiltonian is used as an independent accuracy test of the computations.

3.3. Phase-space visualisation: the Poincaré surface of section

The gravito-inertial ray dynamics is governed by a Hamiltonian with two degrees of freedom, where the phase space is four-dimensional. The fact that ω = H is time independent is equivalent to energy conservation. It implies that phase-space trajectories actually stay on a three-dimensional surface. A PSS is the intersection of all phase-space trajectories with a given three-dimensional surface, defined for example by fixing one coordinate. A PSS at a given frequency is therefore a two-dimensional surface.

Different choices are possible for the PSS, but to provide a complete view of phase space, it must be intersected by most phase-space trajectories. Here we used a PSS defined by fixing the colatitude θ to π/ 2, which corresponds to half of the equatorial plane. Furthermore, to ensure that two trajectories do not intersect on the same point on the PSS, the intersection with θ = π/ 2 is taken into account only when trajectories cross it from one particular side (here from the θ<π/ 2 side). For a given frequency ω, the PSS is computed by following many different trajectories over long time intervals. Then to scan the phase space at a given rotation, the PSS is computed for frequencies spanning the ${\mathrm{[}}^{}\mathrm{0}\mathit{,}{{\mathit{N}}_{\mathrm{0}\mathit{,}\max }^{\mathrm{p}}}^{\mathrm{]}}$ interval. For non-integrable systems, the phase-space structure is complex and the features observed in a PSS are expected to depend on the resolution (in k_r, r or ω) by which the dynamics is investigated. However, modes occupy a finite phase-space volume because Fourier analysis shows that their wave-vector localisation is inversely proportional to their spatial localisation (see Lignières & Georgeot 2009). Investigating phase space with a finite resolution is therefore sufficient to infer mode properties.

4. Gravito-inertial ray dynamics in rotating stars

In this section, we study the evolution of the gravito-inertial ray dynamics with rotation. We first describe the case of pure gravity rays in a non-rotating polytropic model (Sect. 4.1) and then consider gravito-inertial rays under the traditional approximation (Sect. 4.2). These are two simple integrable systems that serve as a reference to study the general case of gravito-inertial rays in a rotating star (Sect. 4.3).

4.1. The non-rotating case Ω = 0

In the non-rotating case, the eikonal equation simplifies into Eq. (21), where L is the second invariant that makes the Hamiltonian system with two degrees of freedom integrable. The phase space of integrable systems is made of nested invariant tori. These structures are said to be invariant because any trajectory that starts on one of the structures stays on it; they are called tori because they have the topology of a two-dimensional torus. Any torus is specified by a value of the doublet (ω, L) and its intersection with the θ = π/ 2 hypersurface is the following f_L²,ω(k_r,r) = 0 curve: $\begin{array}{ccc}{\mathit{k}}_{\mathit{r}}^{\mathrm{2}}\mathrm{=}\left[\frac{{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{\omega }}^{\mathrm{2}}}\mathrm{-}\mathrm{1}\right]\frac{{\mathit{L}}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{-}{\mathit{k}}_{\mathrm{c}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{,}& & \end{array}$(33)trivially derived from the eikonal Eq. (21).

As for any integrable system, two types of invariant tori exist. Irrational tori, on which any ray covers the whole torus, or equivalently fills the f_L²,ω(k_r,r) = 0 curve on the PSS, and rational (or resonant) tori, on which any ray closes on itself before covering the torus. Rays of resonant tori are thus periodic orbits that imprint the PSS on a finite number of points.

Fig. 4 From panel a) to panel c), three PSS at ω/N0,max = [ 0.66676,0.50200,0.05522 ] in a non-rotating polytropic model of star. The fL2,ω(kr,r) = 0 curves shown are the imprint on the PSS of the invariant tori of the integrable gravity ray dynamics. They are drawn here for L = 1 / 2 + ℓ with ℓ taking integer values from 1 to 30 corresponding to degrees of spherical harmonics. The value of L of a given tori can be easily estimated from the intersection with the blue vertical dashed line, since the ordinate of the intersection point is approximately equal to 4L. The EBK quantisation enables us to specify the tori where positive interferences of rays produce modes. For example, the (n,ℓ) = (3,3) gravity mode on panel b) and the (n,ℓ) = (50,3) gravity mode on panel c) are constructed on the tori highlighted in green. The red curve on panel a) is a separatrix also called a homoclinic orbit because it connects the hyperbolic point at (kr,r/R) = (0,0.7968) to itself. As explained in Sect. 4.3.1, it is at the origin of a large chaotic phase-space region observed in the rotating case.

Figure 4 presents three PSS computed at three different frequencies. For integrable systems, modes can be related to the dynamics through the Einstein-Brillouin-Keller (EBK) quantisation procedure that has been applied to non-rotating spherical stars by Gough (1986; see also Gough 1993; Lignières & Georgeot 2009). Quantisation conditions actually select tori that support modes in the sense that rays belonging to these tori interfere positively to construct standing waves, i.e. modes. Accordingly, the degree of the spherical harmonic of the mode, ℓ, is related to the invariant L by the relation L = ℓ + 1 / 2. We therefore chose to represent in Fig. 4 tori with ℓ between 1 and 30. Two of the three PSS have been computed for frequencies corresponding to gravity modes: the (n,ℓ) = (3,3) and the (n,ℓ) = (50,3) modes. On these PSS, the torus associated with the mode is highlighted in green. These examples are useful to gain insight into the relation between the position on the PSS and the modes. In the same spirit, the k_r normalisation has been chosen in such a way that the value of L of a given k_r = f_L²,ω(r) curve can be easily retrieved. Indeed, from the eikonal Eq. (21) and provided that ω ≪ N_0,max, we obtain Rωk_r/N_0,max ≃ LR/r_max, where r_max = 0.28R is the radius of the inner maximum of N₀. Thus, the ordinate of the intersection of the k_r = f_L²,ω(r) curve with the r = 0.25R vertical line is approximatively equal to 4L.

4.2. The traditional approximation

The so-called traditional approximation (Eckart 1960; Lee & Saio 1987) consists in neglecting the contribution of the latitudinal component of the rotation vector in the Coriolis force. The solutions of the perturbation equation are then separable into the coordinates θ and r. This approximation is justified if motions are predominantly horizontal, the effect of the Coriolis force in the radial direction is negligible, and the centrifugal deformation is neglected.

Assuming a spherically symmetric stellar model (thus e_∥ = e_r and e_⊥ = e_θ) and neglecting the fsinθ terms, the eikonal Eq. (20) simplifies into $\begin{array}{ccc}\mathrm{[}{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{-}{\mathit{\omega }}^{\mathrm{2}}\mathrm{]}{\mathit{k}}_{\mathit{\theta }}^{\mathrm{2}}\mathrm{-}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{f}}^{\mathrm{2}}{\cos }^{\mathrm{2}}\mathit{\theta }\mathrm{)}\mathrm{[}{\mathit{k}}_{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{c}}^{\mathrm{2}}\mathrm{(}\mathit{r}{\mathrm{)}}^{\mathrm{]}}\mathrm{=}\mathrm{0.}& & \end{array}$(34)From this equation and the associated dynamical equations, we found a second invariant, namely $\begin{array}{ccc}\mathit{\lambda }\mathrm{=}\frac{{\mathit{\omega }}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}{\mathit{k}}_{\mathit{\theta }}^{\mathrm{2}}}{{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{f}}^{\mathrm{2}}{\cos }^{\mathrm{2}}\mathit{\theta }}\mathrm{=}\frac{{\mathit{\omega }}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}\mathrm{(}{\mathit{k}}_{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{c}}^{\mathrm{2}}\mathrm{)}}{{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}^{\mathrm{2}}}\mathit{,}& & \end{array}$(35)which implies that the associated ray dynamics is integrable. Integrability is expected since the perturbation equation in the traditional approximation is separable. Furthermore, neglecting the k_c term, the condition of propagation derived from Eq. (34) is $\begin{array}{ccc}\mathrm{(}{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}^{\mathrm{2}}\mathrm{)}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{f}}^{\mathrm{2}}{\cos }^{\mathrm{2}}\mathit{\theta }\mathrm{)}\mathit{>}\mathrm{0.}& & \end{array}$(36)For super-inertial waves, the domain of propagation is simply given by ω<N₀. This frequency range is more restricted than that given by Γ > 0, but both are qualitatively similar and even identical in the limit f/N₀ → 0 if the centrifugal deformation is neglected (see Eq. (26)). For sub-inertial waves, the domain of propagation of the traditional approximation consists of two unconnected domains above and below the radius where ω = N₀. Above this radius, the equatorial region defined by fcosθ<ω and below, the polar cone defined by ω<fcosθ. Thus, the propagation of sub-inertial waves from the N₀>f layers towards the central region N₀<f is not allowed by the traditional approximation (for a graphical account of these differences in the domain of propagation, see Fig. 11 of Gerkema et al. 2008). This is a clear limitation of the traditional approximation since, as we have seen from the discussion on the domain of propagation, sub-inertial waves are indeed allowed to go through the centre of the star (and more generally through regions where N₀<f). Near the surface, the condition is $\mathit{\omega }\mathit{<}{\mathit{c}}_{\mathrm{1}}{\mathit{S}}_{\sqrt{\mathit{\lambda }}}$ with ${\mathit{S}}_{\sqrt{\mathit{\lambda }}}\mathrm{=}\sqrt{\mathit{\lambda }}{\mathit{c}}_{\mathrm{s}}\mathit{/}\mathit{r}$, which is a condition formally similar to the non-rotating case, with λ playing the role of L².

Fig. 5 Two PSS at ω/f = [ 2.21,0.567 ] of the ray dynamics in the traditional approximation in a spherical star rotating at Ω = 0.38ΩK. The fλ,ω(kr,r) = 0 curves shown are the imprint on the PSS of the invariant tori of the integrable dynamics. They are drawn for the quantised values of λ(ℓ,ω/f), where ℓ is the degree of the spherical harmonics of the corresponding mode at Ω = 0 and is varied from 1 to 30. The increase of the gap between the fλ,ω(kr,r) = 0 curves as ω decreases is due to the increase of λ with f/ω. Tori in green correspond to gravito-inertial modes labelled by their quantum numbers at zero rotation, that is the (n,ℓ) = (3,3) mode (top) and the (n,ℓ) = (50,3) mode (bottom).

The second expression in Eq. (35) yields the imprint of the invariant tori on the PSS. We observe that it has the form f_λ,ω(k_r,r) = 0 that is the same as in the non-rotating case using λ instead of L². This relation between λ and L² is the ray-dynamics version of the known relation between the separation constant of the perturbation equation in the traditional approximation and the product ℓ(ℓ + 1) (Bouabid et al. 2013). The PSS of the traditional approximation has thus the same structure as the non-rotating PSS. The effect of rotation comes from the increase of λ with f/ω, which produces the dilation effect observed between the two PSS shown in Fig. 5. As for the non-rotating case, gravito-inertial modes computed under the traditional approximation can be associated with specific tori. This is the case in Fig. 5, where two gravito-inertial modes labelled by their quantum numbers at zero rotation, namely (n,ℓ) = (3,3) and (n,ℓ) = (50,3), are represented in green on the PSS.

4.3. Phase-space structure of rotating stars

As rotation increases the gravito-inertial ray dynamics is no longer integrable and undergoes a transition towards a mixed state with coexistence of chaotic regions and invariant tori. This is apparent in Fig. 6, which shows the phase-space structure at Ω = 0.38Ω_K and ω/f = 3.1 together with a selection of rays in the physical space. As compared to the non-rotating PSS, which only exhibits f_L²,ω(k_r,r) = 0 one-dimensional curves, we now observe three different types of structures. First of all, there are one-dimensional curves similar in shape to the non-rotating curves. The blue ray provides an example of ray associated with this type of tori. A second and new type of phase-space structures forms around stable periodic orbits and is called island chains. An orbit is said to be stable when trajectories induced by small perturbations of the orbit remain close to it. The green ray highlights one of the invariant tori formed around a central periodic orbit represented by the light green ray shown in the physical space. Finally, the third type of structure is the chaotic region, which can be easily identified on the PSS by the fact that the associated rays do no stay on a one-dimensional curve but fill a two-dimensional surface instead. The chaotic zone of Fig. 6 is highlighted by the imprint of a ray also shown in the physical space.

The apparition of these new structures in the gravito-inertial dynamics is typical of the KAM-type (in reference to the Kolmogorov-Arnold-Moser theorem) transition of Hamiltonian systems from integrability towards chaos (Ott 1993). Accordingly, in systems that experience infinitesimal departures from integrability, all the resonant tori of the integrable system are destroyed and replaced by island chains around stable periodic orbits and small chaotic regions developing around unstable periodic orbits. Meanwhile, irrational tori are continuously perturbed but not destroyed. The evolution towards larger values of the parameter controlling the departure from integrability has been studied in many cases. The common phenomenology emerging from these studies is that, when the control parameter increases, the invariant tori, either the surviving irrational tori or the island chains, are progressively destroyed while chaotic regions occupy larger phase-space volumes. As demonstrated in the simple case of the kicked rotor (Ott 1993), for each irrational torus there is a critical value of the control parameter above which the torus is destroyed and these critical values serve to quantify the progression towards chaos.

Fig. 6 PSS computed at ω/f = 3.1 for a rotation Ω / ΩK = 0.38 and four gravito-inertial rays in the physical space belonging to the three different types of phase-space structures: a surviving KAM torus (blue), a period-5 island chain (green/light green), and a chaotic zone (red). The light green ray on the middle panel is the stable periodic orbit at the centre of the period-5 island chain. The period of an island chain is defined as the number of corresponding structures on the PSS (see Sect. 4.3.2). These rays cross the PSS on points of corresponding colours. On all figures showing rays in the physical space, the limits of the domain of propagation given by Γ(ω) = 0 are drawn in magenta.

The present gravito-inertial ray dynamics clearly undergoes a KAM-type transition towards a mixed state including surviving irrational tori, island chains, and chaotic regions. Nevertheless, we find that the actual departure from integrability strongly depends on the frequency ω. Indeed, the surface occupied by the island chains and chaotic regions strongly decreases towards low frequencies. This is obvious when comparing the phase-space structure of a Ω = 0.38Ω_K star at three frequencies, ω/f = { 3.1,2.4,0.8 }, as shown in Figs. 6–8, respectively. This observation does not mean that island chains and chaotic regions are absent below some frequency but rather that they are too small to be visible when the PSS are shown on large k_r and r scales. The same phenomenology also holds at higher rotation as shown by the PSS of a Ω = 0.84Ω_K star computed for three decreasing frequencies, ω/f = { 2,1.5,0.5 } and shown in Figs. 9–11, respectively. We therefore conclude that the departure from integrability as measured by the presence of large island chains and chaotic regions strongly diminishes in the low part of the ${\mathrm{[}}^{}\mathrm{0}\mathit{,}{{\mathit{N}}_{\mathrm{0}\mathit{,}\max }^{\mathrm{p}}}^{\mathrm{]}}$ interval and is weakly dependent on the centrifugal deformation.

Overall we may say that the tori that results from the continuous evolution of the non-rotating tori still shape the large-scale organisation of the phase space even at large rotation rates. This is in contrast with the axisymmetric acoustic ray dynamics where, as rotation increases, the phase space is progressively dominated by large islands and a large chaotic region, while irrational tori only survive for large values of k_θ/ω corresponding to whispering gallery modes with low disc-integrated visibilities.

Fig. 7 PSS computed at ω/f = 2.4 for a rotation Ω / ΩK = 0.38 and two gravito-inertial rays belonging to a period-1 island chain including the central stable periodic orbit (light blue).

Fig. 8 PSS computed at ω/f = 0.8 for a rotation Ω / ΩK = 0.38 and two gravito-inertial rays belonging to a period-11 island chain (blue) and to a surviving KAM torus (green), respectively.

In the following, the properties of each of the three types of phase-space structures are described in more detail. Considering them separately is relevant for the study of gravito-inertial modes because, according to Percival’s conjecture (Percival 1973; Berry & Robnik 1984) and following the results obtained for acoustic modes (Lignières & Georgeot 2009), we expect that each type of phase-space structure can be asymptotically related to families of modes with specific properties.

4.3.1. Chaotic regions

PSS computed at different rotation rates show that a large chaotic region is systematically present in the upper part of the ${\mathrm{[}}^{}\mathrm{0}\mathit{,}{{\mathit{N}}_{\mathrm{0}\mathit{,}\max }^{\mathrm{p}}}^{\mathrm{]}}$ frequency interval. In the following we argue that this chaotic region is related to the existence of unstable periodic orbits, which are themselves a direct consequence of the presence of a local minimum of the Brunt-Väisälä frequency in the stellar envelope. This interpretation is supported by our ability to predict the frequency range for which large chaotic regions are seen on the PSS.

In smoothly perturbed integrable Hamiltonian systems, chaotic motions appear near separatrices that connect hyperbolic points, that is unstable periodic orbits or unstable equilibrium points. In the framework of the KAM transition theory, these hyperbolic points naturally originate from the destruction of the resonant tori and chaos first appears around these hyperbolic points and along the structures that connect them (Ott 1993). Chaotic motions, however, can also originate from hyperbolic points already present in the unperturbed system. In our case, the unperturbed system does possess hyperbolic points as can be inferred from the PSS at Ω = 0 shown in panel (a) of Fig. 4. The vicinity of the point (k_r,r) = (0,0.8) is indeed typical of the phase portrait around a hyperbolic point with nearby orbits forming branches of hyperbolas. The separatrix (the red curve in panel (a) of Fig. 4) is also called a homoclinic orbit because it connects the hyperbolic point to itself. The apex of the motion of a pendulum provides a classical example of a hyperbolic point in a system with one degree of freedom, which is in this case a hyperbolic fixed point or equivalently an unstable equilibrium point. In our system with two degrees of freedom, the hyperbolic point is rather an unstable periodic orbit because it corresponds to a circular trajectory at a fixed radius (r ~ 0.8 on panel (a) of Fig. 4). Since the ray dynamics in a spherical star is separable in the radial and latitudinal coordinates, the radius of the hyperbolic point is actually a fixed point of the radial ray dynamics. As in the case of the pendulum, the radial Hamiltonian takes the usual form ${\mathit{H}}_{\mathit{r}}\mathrm{=}\mathrm{1}\mathit{/}\mathrm{2}{\mathit{k}}_{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{V}}_{\mathit{r}}$ and the unstable fixed point corresponds to the maximum of the potential V_r.

In Appendix E, analytical solutions of this problem are found via the envelope approximation, which involves the assumption that the mass enclosed in a sphere of given r does not depend on r, and we show that these solutions provide a very good approximation to the numerical solutions obtained without the envelope approximation. Accordingly, the gravity ray dynamics possesses a hyperbolic point for any frequency ω within the [ 0.93N_0,min,N_0,min ] interval. The radius of the hyperbolic point varies monotonically from r = 0.84R for ω = 0.93N_0,min to r = 0.75R for ω = N_0,min, while the value of L characterising the separatrix varies also monotonically from L = 16.92 for ω = 0.93N_0,min to L = + ∞ for ω = N_0,min.

To confirm that chaos first occurs around the separatrix, we computed gravito-inertial rays close to this orbit for a slowly rotating model Ω / Ω_K = 0.02. As illustrated by Fig. 12, a small chaotic region is found in the expected frequency range and at the expected location. We further make the hypothesis that the large chaotic regions observed at Ω / Ω_K = 0.38 and Ω / Ω_K = 0.84 are due to the same mechanism. To test this hypothesis, PSS have been computed for frequencies inside and outside $\left[\mathrm{0.93}\sqrt{\mathrm{(}{\mathit{N}}_{\mathrm{0}\mathit{,}\min }^{\mathrm{e}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{f}}^{\mathrm{2}}}\mathit{,}{\mathit{N}}_{\mathrm{0}\mathit{,}\min }^{\mathrm{p}}\right]$. To define the lower and upper bounds of this interval, rotational effects have been taken into account on phenomenological grounds. First, rotation changes the frequency range for propagation, and following the conditions (26), a f²sin²Θ term must be added to ${\mathit{N}}_{\mathrm{0}\mathit{,}\min }^{\mathrm{2}}$. Second, we included the ${}^{\mathrm{[}}\mathit{N}_{\mathrm{0}\mathit{,}\min }^{\mathrm{e}}\mathit{,}{{\mathit{N}}_{\mathrm{0}\mathit{,}\min }^{\mathrm{p}}}^{\mathrm{]}}$ interval opened by the centrifugal force. Indeed, we expect that unstable periodic orbits exist for these frequencies because there is always a range of latitudes where the condition $\mathit{\omega }\mathit{/}\sqrt{{\mathit{N}}_{\mathrm{0}\mathit{,}\min }^{\mathrm{2}}\mathrm{+}{\mathit{f}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathrm{\Theta }}\mathrm{\in }\mathrm{\left[}\mathrm{0.93}\mathit{,}\mathrm{1}\mathrm{\right]}$ is fulfilled. Considering frequencies 5% apart from the interval limits, we found that the large chaotic region is present inside and absent outside the $\left[\mathrm{0.93}\sqrt{\mathrm{(}{\mathit{N}}_{\mathrm{0}\mathit{,}\min }^{\mathrm{e}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{f}}^{\mathrm{2}}}\mathit{,}{\mathit{N}}_{\mathrm{0}\mathit{,}\min }^{\mathrm{p}}\right]$ interval.

Fig. 9 PSS computed at ω/f = 2.0 for a rotation Ω / ΩK = 0.84 and two gravito-inertial rays belonging to a period-9 island chain (blue) and to a chaotic zone (red), respectively.

Fig. 10 PSS computed at ω/f = 1.5 for a rotation Ω / ΩK = 0.84 and two gravito-inertial rays belonging to a period-5 island chain including the central stable periodic orbit (light blue).

Fig. 11 PSS computed at ω/f = 0.5 for a rotation Ω / ΩK = 0.84 and a gravito-inertial ray belonging to an island chain of large (>26) period.

Fig. 12 PSS computed at ω/f = 30.5 for a rotation Ω / ΩK = 0.02 showing the development of a small chaotic zone in the vicinity of the separatrix present in the non-rotating case (red curve on panel a) of Fig. 4).

At a given frequency, the extent of the chaotic surface seen on the PSS is not simple to predict quantitatively. Comparing Figs. 12, 6, and 9 we nevertheless observe that, as expected, the surface of the chaotic region increases with rotation.

4.3.2. Island chains

As they originate from the destruction of given rational tori, island chains can be associated with the rational number p/q of their parent torus. Then, the central periodic orbit of the p/q island imprints q times the PSS (Ott 1993); in this study we call the number q the period. The number p can be found by constructing a different PSS (for instance by fixing the radius r instead of the colatitude). For super-inertial rays, p corresponds to the number of loops visible in real space. On the six PSS displayed in Figs. 6–11, we highlighted island chains of various periods: one period-1 at ω/f = 2.4, two of period-5 (but with different p values) at ω/f = 3.1 and ω/f = 1.5, respectively, one period-9 at ω/f = 2, one period-11 at ω/f = 0.8, and one with a period larger than 26 at ω/f = 0.5. As expected from theory, the width of the island decreases with their period. Moreover, there is a clear tendency to find only small and large-period island chains in the sub-inertial regime, although there is no obvious change of behaviour at ω/f = 1.

In a typical KAM transition, island chains form immediately for a non-zero pertubation, grow with the amplitude of the perturbation, and are finally destroyed. We followed the evolution of some low-period island chains with rotation and found that their width slightly increase but no destruction of island chains above a certain rotation rate has been observed. For example, the period-1 island chain is still clearly visible at 99% of the critical velocity.

In the physical space, rays belonging to island chains do not fill their resonant cavity as they remain concentrated around the stable periodic trajectory. From Ballot et al. (2012), we already know that numerically computed modes can be associated with island chains. A detailed confrontation of the island properties with numerically computed modes will be performed in a future paper.

4.3.3. Surviving tori

The irrational tori of the non-rotating case are continuously deformed by the effect of rotation but may nevertheless have survived destruction even at high rotation rates. This is the case, in particular, for high-L tori, since for most frequencies no large island chain or chaotic region is observed in the PSS in this domain, which corresponds to high values of the normalised k_r. Some irrational tori with lower L also survive the perturbation induced by rotation. This is obvious at low frequencies where phase space is very much structured by one-dimensional curves broadly similar to the non-rotating curves (see Figs. 8 and 11 and panel (c) of Fig. 4).

This nearly integrable behaviour of the low-frequency, gravito-inertial, ray dynamics might be attributed to the proximity of a hidden integrable system. To test whether the ray dynamics in the traditional approximation could be this integrable system, we computed the “traditional” invariant λ along the path of gravito-inertial rays. Considering a sub-inertial ray in a Ω / Ω_K = 0.38 star, we found that λ varies by several orders of magnitude along a given ray. As expected, the largest variations occur in regions of space that are outside the “traditional” domain of propagation, but even within the resonant cavity of the traditional approximation, λ still varies by a factor of two. In the super-inertial regime, the variations of λ are of the order of 1.5, which is of the same order of magnitude as the variations of L; the parameter L is the invariant of the non-rotating case, which is not supposed to provide a valid approximation. We conclude that the ray dynamics computed in the traditional approximation does not provide an accurate approximation even in the low-frequency regime, where the full gravito-inertial ray dynamics is close to integrable.

5. From rays to modes

In this section, we summarise the properties of gravito-inertial modes that can be inferred from the present study of the ray dynamics. Methods to construct modes from the superposition of positively interfering rays have been first proposed for integrable systems. The EBK procedure has been established in the field of quantum physics as a generalisation of the Bohr-Sommerfeld method to quantise systems from their classical trajectories. This method only works in the presence of invariant tori and leads to a regular frequency spectrum in the sense that the spectrum can be described by a function of N integers, where N is the number of degrees of freedom of the system. For chaotic dynamical systems, the EBK quantisation does not work and must be replaced by the Gutzwiller trace formula, which is difficult to apply in practice (Gutzwiller 1990). Nevertheless, generic statistical properties of spectra associated with chaotic dynamical systems, such as the distribution of spacings between consecutive frequencies, have been discovered in quantum systems (Bohigas et al. 1984). For mixed dynamical systems, it has been conjectured that chaotic regions and phase-space regions structured by invariant tori, such as island chains or regions of surviving irrational tori, can be quantised separately. Consequently, the spectrum of mixed systems is expected to be a superposition of frequency subsets, each of which are associated with different phase-space structures. Frequency subsets associated with regions structured by invariant tori should be regular while those associated with chaotic regions should possess generic statistical properties (Percival 1973; Berry & Robnik 1984). A more detailed description of the concepts presented above together with the full references is available in Lignières & Georgeot (2009). Percival’s predictions have been applied to acoustic rays in rotating stars and successfully confronted with numerically computed high-frequency acoustic modes (Lignières & Georgeot 2009).

In our case, Percival’s conjecture implies that the frequency spectrum of short-wavelength gravito-inertial modes is a superposition of frequency subsets associated with the phase-space structures of the gravito-inertial ray dynamics, notably a regular sub-spectrum associated with the surviving irrational tori, regular sub-spectra associated with island chains, and a sub-spectrum with generic statistical properties associated with the large chaotic region. In principle, quantitative predictions on the regular sub-spectra can be obtained by conducting the EBK quantisation of the invariant tori. This has been carried out in Pasek et al. (2012) for the period-2 island chain of the acoustic ray dynamics and could be also attempted for the gravito-inertial island chains and surviving tori.

Previous numerical explorations of gravito-inertial modes already provide some insights into the relation between the asymptotic predictions of the gravito-inertial dynamics and the numerically computed modes. Using μ = 3 polytropic models, Ballot et al. (2010) followed low-degree (ℓ ≤ 3) modes in both low- and high-order ranges from Ω = 0 to Ω = 0.7Ω_K. The vast majority of these modes undergo a smooth evolution of their spatial distribution with rotation, in which the main change is the equatorial mode confinement in the sub-inertial regime. These modes most probably correspond to the phase-space region structured by the surviving irrational tori. On the other hand, some modes experience a striking modification of their spatial distribution, leading to a concentration of their energy. The connection between these modes called rosette modes and the central periodic orbits of island chains has been clearly established in Ballot et al. (2012). Since then, Takata & Saio (2013) proposed a different approach in which rosette modes result from the rotationally induced coupling of modes, which were nearly degenerate gravity modes at Ω = 0 (see also Saio & Takata 2014; Takata 2014). The link between these two apparently very different approaches to model rosette modes remains to be understood. In our study, the number K introduced by Takata & Saio (2013) to characterise families of rosette modes corresponds to the ratio 2q/p.² The authors found only integer values of K, which they justified by the fact that the only effect of rotation they considered was the Coriolis force. We also present island chains with non-integer values. This is probably because we took the centrifugal deformation into account. Modes that could be related to the large chaotic region that we describe have not been identified yet. By providing an estimate of both the frequencies and wave vectors of chaotic rays, ray dynamics will help to find these modes. More generally, a dedicated numerical study will be necessary to test the predictions of the ray-based asymptotic theory. In this context, tools allowing us to construct phase-space representations for numerically computed modes, such as Husimi distributions, are available. As illustrated by Lignières & Georgeot (2009), those can be used to show the imprint of modes on the PSS.

6. Discussion and conclusions

The main results presented in the previous sections are (i) the derivation of an eikonal equation for gravito-inertial rays in polytropic uniformly rotating stars; (ii) the analysis of the properties of this equation, including constraints on the domains of propagation for k_φ = 0 rays; (iii) the numerical computation of the k_φ = 0 ray dynamics for μ = 3 polytropic models with increasing rotation, and (iv) the interpretation of the numerical results using tools and concepts of Hamiltonian dynamics. In particular, the exploration of the phase space provides us with a global view of the properties of gravito-inertial waves in rotating stars. Qualitative predictions on the frequency spectrum of gravito-inertial modes have been made and should be a powerful tool to interpret the properties of numerically computed gravito-inertial modes. A quantitative result already provided by the ray theory is the domain of propagation. Other quantitative predictions are expected from the EBK quantisation of the surviving irrational tori and of the island chains.

A first restriction of the present work is that we did not consider k_φ ≠ 0 rays. This would be necessary for an asymptotic study of non-axisymmetric gravito-inertial modes. Another restriction is the use of μ = 3 polytropic models although the formalism can be extended to more realistic rotating stellar models. In particular, the effects of a convective core and/or a convective envelope would be interesting to consider.

The present ray dynamics formalism can be improved by taking the effect of background mean motions and, in particular, differential rotation into account. For example, self-consistent models including baroclinic effects computed by the ESTER code (Rieutord & Espinosa Lara 2013) could be used. According to atmospheric or oceanic studies, the eikonal equation including mean horizontal motions is identical to (18), except that the frequency is changed into the Doppler-shifted frequency, which is the frequency that would be observed in a frame of reference moving with the mean horizontal motions (Fritts & Alexander 2003). Analagous with the effect of wind shear, differential rotation should cause some back-refraction and thus will modify the resonant cavity. Another improvement to the eikonal equation would be to include non-axisymmetric Rossby waves.

Ray models might also be used to study the transport of angular momentum by waves. This would require us to extend the ray formalism to treat the evolution of the wave amplitude A along the ray path (the next order of the WKB approximation). Such amplitude equations are already used to investigate the driving of the atmospheric circulation by waves (Marks & Eckermann 1995). To conduct similar investigations in stars, it will be necessary to model energy deposition processes such as thermal dissipation and critical layers (Fritts et al. 1998; Alvan et al. 2013) and to take phase changes at caustics into account (Broutman et al. 2004).

Angular momentum transport by gravito-inertial waves that occasionally escape the star has been invoked as a possible explanation of the Be phenomenon (Ishimatsu & Shibahashi 2013). Specific properties of chaotic gravito-inertial waves such as those observed in our calculations could be interesting in this context. At any location and in particular at the surface, the horizontal wavelength k_⊥ of chaotic waves can indeed take many different values. Thus a chaotic wave can be back-refracted many times at the surface before it escapes the star because k_⊥ has reached a high enough value.

When compared to the ray theory of Dintrans & Rieutord (2000), the present work includes the effects of centrifugal deformation as well as a realistic back-refraction of rays in the stellar envelope. This last point has a crucial effect on the comparison between the two approaches. In Dintrans & Rieutord (2000), the rays of characteristic are confined into the star by a rigid and stress-free surface and this boundary condition destroys the Hamiltonian character of the dynamics of characteristics. As a consequence, the volume of phase-space elements is not conserved, which in turns enables rays to focus on attractors. This focusing effect is at the origin of the singular gravito-inertial modes that are effectively observed in the presence of rigid boundaries (Maas et al. 1997; Dintrans et al. 1999). However, there are no such rigid boundaries in real stars and singular gravito-inertial modes focused on attractors of characteristics by unrealistic rigid boundaries should not exist in stars. As we have seen in the present paper, the ray dynamics ought to be Hamiltonian in stars as long as dissipative processes are neglected.

Acknowledgments

The authors thank Michel Rieutord for his careful reading of the manuscript, and an anonymous referee for contributing to the improvement of this paper. This work was supported by the Centre National de la Recherche Scientifique (CNRS) through the Programme National de Physique Stellaire (PNPS).

References

Appendix A: Derivation of the perturbation Eq. (11)

In this appendix, our goal is to transform the perturbation Eqs. (5)–(7) into Eq. (11) for . In Sect. A.1, density and velocity perturbations appearing in Eqs. (5)–(7) are eliminated to obtain an equation for pressure perturbations only. Then in Sect. A.2, the new variable is introduced to eliminate the first-order derivatives in the meridional plane from the perturbation equation on . Finally, in Sects. A.3 and A.4, the constant term is expanded in powers of ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$ to determine the dominant $\mathrm{𝒪}\mathrm{(}\mathrm{1}\mathit{/}{\mathit{H}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{)}$ approximation of .

Appendix A.1: Derivation of an equation for pressure perturbations only

We first rewrite the perturbation Eqs. (5)–(7), using u = ρ₀v and the definition of the Brunt-Väisälä frequency ${\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{=}{{g}}_{\mathrm{0}}\mathrm{·}\left(\frac{{\nabla }{\mathit{\rho }}_{\mathrm{0}}}{{\mathit{\rho }}_{\mathrm{0}}}\mathrm{-}\frac{\mathrm{1}}{{\mathrm{\Gamma }}_{\mathrm{1}}}\frac{{\nabla }{\mathit{P}}_{\mathrm{0}}}{{\mathit{P}}_{\mathrm{0}}}\right)$ as follows: $\begin{array}{ccc}{\mathit{\partial }}_{\mathit{t}}\mathit{\rho }& \mathrm{+}& \\ {\mathit{\partial }}_{\mathit{t}}{u}& \mathrm{+}& \\ {\mathit{\partial }}_{\mathit{t}}\mathit{P}& \mathrm{=}& \end{array}$In order to eliminate u from these equations, we first calculate f ∧() as ${f}\mathrm{\wedge }{\mathit{\partial }}_{\mathit{t}}{u}\mathrm{+}\mathrm{(}{f}\mathrm{·}{u}\mathrm{)}{f}\mathrm{-}{\mathit{f}}^{\mathrm{2}}{u}\mathrm{=}\mathrm{-}{f}\mathrm{\wedge }{\nabla }\mathit{P}\mathrm{+}\mathit{\rho }{f}\mathrm{\wedge }{{g}}_{\mathrm{0}}\mathit{.}$(A.4)The first term on the LHS is then developed using ∂_t() to give ${\mathit{\partial }}_{\mathit{tt}}^{\mathrm{2}}{u}\mathrm{+}{\mathit{f}}^{\mathrm{2}}{u}\mathrm{-}\mathrm{(}{f}\mathrm{·}{u}\mathrm{)}{f}\mathrm{=}\mathrm{-}{\nabla }{\mathit{\partial }}_{\mathit{t}}\mathit{P}\mathrm{+}{\mathit{\partial }}_{\mathit{t}}\mathit{\rho }{{g}}_{\mathrm{0}}\mathrm{+}{f}\mathrm{\wedge }{\nabla }\mathit{P}\mathrm{-}\mathit{\rho }{f}\mathrm{\wedge }{{g}}_{\mathrm{0}}\mathit{.}$(A.5)We use a similar method to eliminate the scalar product f·u. We compute f·()${f}\mathrm{·}{\mathit{\partial }}_{\mathit{t}}{u}\mathrm{=}\mathrm{-}{f}\mathrm{·}{\nabla }\mathit{P}\mathrm{+}\mathit{\rho }{f}\mathrm{·}{{g}}_{\mathrm{0}}\mathit{,}$(A.6)and introduce it into ∂_t() to get ${\mathit{\partial }}_{\mathit{ttt}}^{\mathrm{3}}{u}\mathrm{+}{\mathit{f}}^{\mathrm{2}}{\mathit{\partial }}_{\mathit{t}}{u}\mathrm{=}\mathrm{-}\mathrm{(}{f}\mathrm{·}{\nabla }\mathit{P}\mathrm{)}{f}\mathrm{+}\mathit{\rho }\mathrm{(}{f}\mathrm{·}{{g}}_{\mathrm{0}}\mathrm{)}{f}\mathrm{-}{\nabla }{\mathit{\partial }}_{\mathit{tt}}^{\mathrm{2}}\mathit{P}\mathrm{+}{\mathit{\partial }}_{\mathit{tt}}^{\mathrm{2}}\mathit{\rho }{{g}}_{\mathrm{0}}\mathrm{+}{f}\mathrm{\wedge }{\nabla }{\mathit{\partial }}_{\mathit{t}}\mathit{P}\mathrm{-}{\mathit{\partial }}_{\mathit{t}}\mathit{\rho }{f}\mathrm{\wedge }{{g}}_{\mathrm{0}}\mathit{.}$(A.7)Defining the operator $\mathrm{ℒ}\mathrm{\equiv }{\mathit{\partial }}_{\mathit{ttt}}^{\mathrm{3}}\mathrm{+}{\mathit{f}}^{\mathrm{2}}\mathit{\partial }t$(A.8)enables us to rewrite Eqs. (A.1) and (A.3) as $\begin{array}{ccc}{\mathit{\partial }}_{\mathit{t}}\mathrm{ℒ}\mathrm{\left(}\mathit{\rho }\mathrm{\right)}& \mathrm{=}& \\ {\mathit{\partial }}_{\mathit{t}}\mathrm{ℒ}\mathrm{\left(}\mathit{P}\mathrm{\right)}& \mathrm{=}& \end{array}$where ℒ(u) is given by Eq. (A.7) and α, defined by $\mathit{\alpha }\mathrm{=}\frac{{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{2}}{{\mathit{N}}_{\mathrm{0}}}^{\mathrm{2}}}{{{\mathit{g}}_{\mathrm{0}}}^{\mathrm{2}}}\mathrm{=}\frac{\mathit{\mu }{\mathrm{\Gamma }}_{\mathrm{1}}}{\mathit{\mu }\mathrm{+}\mathrm{1}}\mathrm{-}\mathrm{1}\mathit{,}$(A.11)is uniform in a polytropic model. Replacing ℒ(u) in Eq. (A.10) by the RHS of Eq. (A.7), we get ${\mathit{\partial }}_{\mathit{t}}\mathrm{ℒ}\mathrm{\left(}\mathit{P}\mathrm{\right)}\mathrm{+}\mathit{\alpha }\mathrm{\left[}\mathrm{\right(}{f}\mathrm{·}{{g}}_{\mathrm{0}}\mathrm{\left)}\mathrm{\right(}{f}\mathrm{·}{\nabla }\mathit{P}\mathrm{)}\mathrm{+}\mathrm{(}{f}\mathrm{\wedge }{{g}}_{\mathrm{0}}\mathrm{)}\mathrm{·}{\nabla }{\mathit{\partial }}_{\mathit{t}}\mathit{P}\mathrm{+}{{g}}_{\mathrm{0}}\mathrm{·}{\nabla }{\mathit{\partial }}_{\mathit{tt}}^{\mathrm{2}}{\mathit{P}}^{\mathrm{]}}\mathrm{=}{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{2}}{\mathit{\partial }}_{\mathit{t}}\mathrm{ℒ}\mathrm{\left(}\mathit{\rho }\mathrm{\right)}\mathrm{+}\mathit{\alpha }\mathrm{\left[}\mathrm{\right(}{f}\mathrm{·}{{g}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}\mathit{\rho }\mathrm{+}{{\mathit{g}}_{\mathrm{0}}}^{\mathrm{2}}{\mathit{\partial }}_{\mathit{tt}}^{\mathrm{2}}{\mathit{\rho }}^{\mathrm{]}}\mathrm{=}{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{2}}\mathrm{𝒟}\mathrm{\left(}\mathit{\rho }\mathrm{\right)}\mathit{,}$(A.12)where $\mathrm{𝒟}\mathrm{\equiv }{\mathit{\partial }}_{\mathit{tttt}}^{\mathrm{4}}\mathrm{+}\mathrm{(}{\mathit{f}}^{\mathrm{2}}\mathrm{+}{{\mathit{N}}_{\mathrm{0}}}^{\mathrm{2}}\mathrm{)}{\mathit{\partial }}_{\mathit{tt}}^{\mathrm{2}}\mathrm{+}\frac{{{\mathit{N}}_{\mathrm{0}}}^{\mathrm{2}}\mathrm{(}{f}\mathrm{·}{{g}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{{{\mathit{g}}_{\mathrm{0}}}^{\mathrm{2}}}\mathrm{·}$(A.13)To simplify calculations, we assume that P and ρ are harmonic in time, i.e. and . Then the operator is simply a multiplication by $\mathit{D}\mathrm{=}{\mathit{\omega }}^{\mathrm{4}}\mathrm{-}{\mathit{\omega }}^{\mathrm{2}}\mathrm{(}{\mathit{f}}^{\mathrm{2}}\mathrm{+}{{\mathit{N}}_{\mathrm{0}}}^{\mathrm{2}}\mathrm{)}\mathrm{+}\frac{{{\mathit{N}}_{\mathrm{0}}}^{\mathrm{2}}\mathrm{(}{f}\mathrm{·}{{g}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{{{\mathit{g}}_{\mathrm{0}}}^{\mathrm{2}}}\mathit{,}$(A.14)and ∂_tℒ a multiplication by $\mathit{G}\mathrm{=}{\mathit{\omega }}^{\mathrm{2}}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{f}}^{\mathrm{2}}\mathrm{)}\mathit{.}$(A.15)From Eq. (A.12), is expressed as a function of as follows: $\mathit{\rho ̂}\mathrm{=}\frac{\mathit{G}\mathit{P̂}\mathrm{+}{H}\mathrm{·}{\nabla }\mathit{P̂}}{{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{2}}\mathit{D}}\mathit{,}$(A.16)with ${H}\mathrm{=}\mathit{\alpha }\mathrm{\left[}\mathrm{\right(}{f}\mathrm{·}{{g}}_{\mathrm{0}}\mathrm{)}{f}\mathrm{-}\mathrm{i}\mathit{\omega }{f}\mathrm{\wedge }{{g}}_{\mathrm{0}}\mathrm{-}{\mathit{\omega }}^{\mathrm{2}}{{g}}_{\mathrm{0}}\mathrm{]}\mathit{.}$(A.17)We then compute ∇·ℒ(u) from Eq. (A.7), using the general identity ∇·(ba) = a·∇b + b∇·a and the fact that f is uniform ${\nabla }\mathrm{·}\mathrm{ℒ}\mathrm{\left(}{u}\mathrm{\right)}\mathrm{=}\mathrm{-}{f}\mathrm{·}{\nabla }\mathrm{(}{f}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}\mathrm{-}\mathrm{i}\mathit{\omega }{\nabla }\mathrm{·}\mathrm{(}{f}\mathrm{\wedge }{\nabla }\mathit{P̂}\mathrm{)}\mathrm{+}{\mathit{\omega }}^{\mathrm{2}}\mathrm{\Delta }\mathit{P̂}\mathrm{+}{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho ̂}{F}\mathrm{\right)}\mathit{,}$(A.18)with ${F}\mathrm{=}\mathrm{(}{f}\mathrm{·}{{g}}_{\mathrm{0}}\mathrm{)}{f}\mathrm{+}\mathrm{i}\mathit{\omega }{f}\mathrm{\wedge }{{g}}_{\mathrm{0}}\mathrm{-}{\mathit{\omega }}^{\mathrm{2}}{{g}}_{\mathrm{0}}\mathit{.}$(A.19)The expression (A.18) can be simplified using ${\nabla }\mathrm{·}\mathrm{(}{f}\mathrm{\wedge }{\nabla }\mathit{P̂}\mathrm{)}\mathrm{=}{\nabla }\mathit{P̂}\mathrm{·}\mathrm{(}{\nabla }\mathrm{\wedge }{f}\mathrm{)}\mathrm{-}{f}\mathrm{·}\mathrm{(}{\nabla }\mathrm{\wedge }{\nabla }\mathit{P̂}\mathrm{)}\mathrm{=}\mathrm{0}\mathit{,}$(A.20)which finally yields ${\nabla }\mathrm{·}\mathrm{ℒ}\mathrm{\left(}{u}\mathrm{\right)}\mathrm{=}{\mathit{\omega }}^{\mathrm{2}}\mathrm{\Delta }\mathit{P̂}\mathrm{-}{f}\mathrm{·}{\nabla }\mathrm{(}{f}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}\mathrm{+}{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho ̂}{F}\mathrm{\right)}\mathit{.}$(A.21)Replacing the expression (A.21) in Eq. (A.9) and using Eq. (A.16) to express , we obtain the following equation for : $\frac{\mathit{G}\mathrm{(}\mathit{G}\mathit{P̂}\mathrm{+}{H}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}}{{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{2}}\mathit{D}}\mathrm{+}{\mathit{\omega }}^{\mathrm{2}}\mathrm{\Delta }\mathit{P̂}\mathrm{-}{f}\mathrm{·}{\nabla }\mathrm{(}{f}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}\mathrm{+}{\nabla }\mathrm{·}\left[\frac{\mathrm{(}\mathit{G}\mathit{P̂}\mathrm{+}{H}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}{F}}{{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{2}}\mathit{D}}\right]\mathrm{=}\mathrm{0.}$(A.22)We use the identity ∇·(a ∧ b) = b·(∇ ∧ a)−a·(∇ ∧ b), which implies that ∇·(f ∧ g₀) = 0, as f is uniform and g₀ derives from a potential. Defining the vector F₀ = (f·g₀)f−ω²g₀, the vectors F and H are written as F = F₀ + iωf ∧ g₀ and H = α(F₀−iωf ∧ g₀). Thus, Eq. (A.22) becomes $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\frac{\mathit{G}\mathrm{(}\mathit{G}\mathit{P̂}\mathrm{+}\mathit{\alpha }{{F}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}}{{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{2}}\mathit{D}}\mathrm{+}{\mathit{\omega }}^{\mathrm{2}}\mathrm{\Delta }\mathit{P̂}\mathrm{-}{f}\mathrm{·}{\nabla }\mathrm{(}{f}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}\mathrm{+}{\nabla }\mathrm{·}\left[\frac{\mathrm{(}\mathit{G}\mathit{P̂}\mathrm{+}\mathit{\alpha }{{F}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}{{F}}_{\mathrm{0}}}{{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{2}}\mathit{D}}\right]\mathrm{+}\mathrm{i}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }\mathrm{)}\frac{\mathit{\omega G}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{D}}\mathrm{(}{f}\mathrm{\wedge }{{g}}_{\mathrm{0}}\mathrm{)}\mathrm{·}{\nabla }\mathit{P̂}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array} \mathrm{-}\mathrm{i}\mathit{\alpha \omega }{\nabla }\mathrm{·}\left[\mathrm{(}{f}\mathrm{\wedge }{{g}}_{\mathrm{0}}\mathrm{)}\mathrm{·}{\nabla }\mathit{P̂}\frac{{F}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{D}}\right]\mathrm{+}\mathrm{i}\mathit{\omega }\frac{\mathit{\alpha }}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{D}}\mathrm{(}{f}\mathrm{\wedge }{{g}}_{\mathrm{0}}\mathrm{)}\mathrm{·}{\nabla }\mathrm{(}{{F}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}\mathrm{=}\mathrm{0}\mathit{,}\end{array}$(A.23)where we have used the fact that α is uniform and f ∧ g₀ is parallel to e_φ. We then use the relation ${{F}}_{\mathrm{0}}\mathrm{·}{\nabla }\left(\frac{{{F}}_{\mathrm{0}}}{{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{2}}\mathit{D}}\mathrm{·}{\nabla }\mathit{P̂}\right)\mathrm{=}\frac{\mathrm{1}}{{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{2}}\mathit{D}}{{F}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathrm{(}{{F}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}\mathrm{+}\left[{{F}}_{\mathrm{0}}\mathrm{·}{\nabla }\left(\frac{\mathrm{1}}{{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{2}}\mathit{D}}\right)\right]{{F}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathit{P̂}$(A.24)to get $\begin{array}{ccc}\mathit{D}{\mathit{\omega }}^{\mathrm{2}}\mathrm{\Delta }\mathit{P̂}\mathrm{-}\mathit{D}{f}\mathrm{·}{\nabla }\mathrm{(}{f}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}& \mathrm{+}& \frac{{{\mathit{N}}_{\mathrm{0}}}^{\mathrm{2}}}{{{\mathit{g}}_{\mathrm{0}}}^{\mathrm{2}}}{{F}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathrm{(}{{F}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}\mathrm{+}\frac{\mathit{N}}{{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{2}}}{{F}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{+}\frac{\mathit{M}}{{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{2}}}\mathit{P̂}\mathrm{+}\mathrm{i}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }\mathrm{)}\frac{\mathit{\omega G}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\mathrm{(}{f}\mathrm{\wedge }{{g}}_{\mathrm{0}}\mathrm{)}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{-}\mathrm{i}\mathit{\alpha \omega D}{\nabla }\mathrm{·}\left[\mathrm{(}{f}\mathrm{\wedge }{{g}}_{\mathrm{0}}\mathrm{)}\mathrm{·}{\nabla }\mathit{P̂}\frac{{F}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{D}}\right]\\ & \mathrm{+}& \mathrm{i}\mathit{\omega }\frac{\mathit{\alpha }}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\mathrm{(}{f}\mathrm{\wedge }{{g}}_{\mathrm{0}}\mathrm{)}\mathrm{·}{\nabla }\mathrm{(}{{F}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}\mathrm{=}\mathrm{0}\mathit{,}\end{array}$(A.25)where $\begin{array}{ccc}\mathit{M}& \mathrm{=}& \mathit{G}\left[\mathit{G}\mathrm{+}{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{2}}\mathit{D}{\nabla }\mathrm{·}\left(\frac{{{F}}_{\mathrm{0}}}{{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{2}}\mathit{D}}\right)\right]\mathit{,}\\ \mathit{N}& \mathrm{=}& \mathrm{(}\mathrm{1}\mathrm{+}\mathit{\alpha }\mathrm{)}\mathit{G}\mathrm{+}{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{2}}\mathit{D}{\nabla }\left(\frac{\mathit{\alpha }{{F}}_{\mathrm{0}}}{{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{2}}\mathit{D}}\right)\mathrm{·}\end{array}$From the definition of F₀, we have $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{{F}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathrm{(}{{F}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}\mathrm{=}\mathrm{(}{f}\mathrm{·}{{g}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}{f}\mathrm{·}{\nabla }\mathrm{(}{f}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}\mathrm{+}{\mathit{\omega }}^{\mathrm{4}}{{g}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathrm{(}{{g}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}\mathrm{-}{\mathit{\omega }}^{\mathrm{2}}\mathrm{(}{f}\mathrm{·}{{g}}_{\mathrm{0}}\mathrm{)}\left[{f}\mathrm{·}{\nabla }\mathrm{(}{{g}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}\mathrm{+}{{g}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathrm{(}{f}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}\right]\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array} \mathrm{+}\mathrm{[}\mathrm{(}{f}\mathrm{·}{{g}}_{\mathrm{0}}\mathrm{)}{f}\mathrm{·}{\nabla }\mathrm{(}{f}\mathrm{·}{{g}}_{\mathrm{0}}\mathrm{)}\mathrm{-}{\mathit{\omega }}^{\mathrm{2}}{{g}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathrm{(}{f}\mathrm{·}{{g}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{]}}{f}\mathrm{·}{\nabla }\mathit{P̂}\mathit{,}\end{array}$(A.28)where, thanks to the fact that f is uniform, the last term can be also written $\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{f}\mathrm{·}{\nabla }\mathit{K}\mathrm{)}{f}\mathrm{·}{\nabla }\mathit{P̂}$ with $\mathit{K}\mathrm{=}\mathrm{(}{f}\mathrm{·}{{g}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}^{\mathrm{2}}{{\mathit{g}}_{\mathrm{0}}}^{\mathrm{2}}\mathrm{=}{{F}}_{\mathrm{0}}\mathrm{·}{{g}}_{\mathrm{0}}\mathit{.}$(A.29)In the end, we get $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\mathit{D}\mathrm{\Delta }\mathit{P̂}\mathrm{-}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{f}}^{\mathrm{2}}\mathrm{-}{{\mathit{N}}_{\mathrm{0}}}^{\mathrm{2}}\mathrm{)}{f}\mathrm{·}{\nabla }\mathrm{(}{f}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}\mathrm{+}\frac{{{\mathit{N}}_{\mathrm{0}}}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}}{{{\mathit{g}}_{\mathrm{0}}}^{\mathrm{2}}}{{g}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathrm{(}{{g}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}\mathrm{-}\frac{{{\mathit{N}}_{\mathrm{0}}}^{\mathrm{2}}}{{{\mathit{g}}_{\mathrm{0}}}^{\mathrm{2}}}\mathrm{(}{f}\mathrm{·}{{g}}_{\mathrm{0}}\mathrm{)}\mathrm{[}{f}\mathrm{·}{\nabla }\mathrm{(}{{g}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}\mathrm{+}{{g}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathrm{(}{f}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{\left)}\mathrm{\right]}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array} \mathrm{+}\frac{\mathit{N}{f}\mathrm{·}{{g}}_{\mathrm{0}}\mathrm{+}\frac{\mathit{\alpha }}{\mathrm{2}}{f}\mathrm{·}{\nabla }\mathit{K}}{{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}}{f}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{-}\frac{\mathit{N}}{{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{2}}}{{g}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{+}\frac{\mathit{M}}{{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}}\mathit{P̂}\mathrm{+}\mathrm{i}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }\mathrm{)}\frac{\mathit{G}}{\mathit{\omega }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\mathrm{(}{f}\mathrm{\wedge }{{g}}_{\mathrm{0}}\mathrm{)}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{-}\mathrm{i}\mathit{\alpha }\frac{\mathit{D}}{\mathit{\omega }}{\nabla }\mathrm{·}\left[\mathrm{(}{f}\mathrm{\wedge }{{g}}_{\mathrm{0}}\mathrm{)}\mathrm{·}{\nabla }\mathit{P̂}\frac{{F}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{D}}\right]\\ \begin{array}{c}\end{array} \mathrm{+}\mathrm{i}\frac{\mathit{\alpha }}{\mathit{\omega }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\mathrm{(}{f}\mathrm{\wedge }{{g}}_{\mathrm{0}}\mathrm{)}\mathrm{·}{\nabla }\mathrm{(}{{F}}_{\mathrm{0}}\mathrm{·}{\nabla }\mathit{P̂}\mathrm{)}\mathrm{=}\mathrm{0}\mathit{,}\end{array}$(A.30)which is an equation for a single unknown, as we were looking for.