© ESO 2018

1 Introduction

Differential rotation plays a key role in stellar dynamics, magnetism, and evolution. It triggers instabilities (such as Rayleigh-Taylor, Goldreich-Schubert-Friecke, or shear instabilities) and interacts with large-scale flows (such as meridional circulation) that induce transport of chemicals and angular momentum in stellar interiors (e.g. Zahn 1992; Maeder & Zahn 1998; Mathis & Zahn 2004). This transport greatly modifies the structural, rotational, chemical and magnetic evolution of stars (e.g. Maeder 2009; Mathis et al. 2013; Rieutord 2006, and references therein). It is therefore crucial to constrain the amount of differential rotation present in stars to better understand these transport processes. Further, differential rotation can drastically impact the propagation and the frequency spectrum of gravity waves, especially in the case of strong gradients of angular velocity (Ando 1985; Lee & Saio 1993; Mathis 2009; Mirouh et al. 2016; Guenel et al. 2016). This is of great importance for stars that exhibit gravity waves, since these waves allow for probing stellar interiors and rotation thanks to seismic diagnoses based on the asymptotic properties of the waves (e.g. Bouabid et al. 2013; Prat et al. 2017). Moreover, these waves are able to transport angular momentum. In particular, they provide a mechanism to explain the weak differential rotation revealed by the CoRoT and Kepler space missions in solar-type (Schatzman 1993; Zahn et al. 1997; Talon & Charbonnel 2005), evolved (Talon & Charbonnel 2008; Fuller et al. 2014; Pinçon et al. 2017) and massive stars (Lee et al. 2014; Rogers et al. 2013; Rogers 2015).

Seismology of slow rotators, in which rotation can be considered as a perturbation of the non-rotating case, provides us with many constraints on their internal structure and rotation. Within this approximation, the main effect of rotation on oscillations is to lift the degeneracy between modes of the same radial order and degree but different azimuthal orders (Saio 1981). These rotational splittings depend directly on the internal rotation profile of the stars (Aerts et al. 2010, and references therein).

In the case of the Sun, splittings of solar p modes have highlighted a latitudinal differential rotation in the convective envelope and a relatively flat rotation profile in the radiative zone down to 0.2 R_⊙ (Thompson et al. 1996; Couvidat et al. 2003). Splittings of g-mode candidates also suggested that the solar core rotates around four times as fast as the bulk of the radiative zone (García et al. 2007). This has recently been confirmed by Fossat et al. (2017) using frequency modulations of p modes by low-frequency g modes.

For distant stars, which we cannot spatially resolve, progress has also been made in this direction. Many subgiant and red giant stars exhibit rotational splittings of mixed modes that allow us to estimate the contrast in rotation between their core and their surface (Beck et al. 2012; Mosser et al. 2012; Deheuvels et al. 2012, 2014, 2015; Triana et al. 2017). It is found that the core of such stars typically rotates between 2 and 20 times as fast as their envelope. Constraints on differential rotation inside solar-type stars can be put by combining information from rotational splittings with estimates of the surface rotation rate obtained using other methods (Benomar et al. 2015; Nielsen et al. 2017). In some more massive stars, rotational splittings of g modes provide constraints on their internal rotation (Triana et al. 2015; Murphy et al. 2016). When p-mode splittings are also observed, the contrast in rotation between the core and the envelope can be estimated (Kurtz et al. 2014; Saio et al. 2015). Some of these stars exhibit a core that rotates slower than the envelope, down to a rotation ratio of 30%.

For fast rotators, such as γ Doradus, δ Scuti, slowly pulsating B-type (SPB), β Cephei or Be stars, however, perturbative methods fail (Reese et al. 2006; Ballot et al. 2010, 2013) and a more complete treatment of rotation is needed. The eigenvalue problem is fully two-dimensional (2D), in the sense that it cannot beseparated into one-dimensional problems as in the non-rotating case (Rieutord 2009). In general, the problemis not solvable analytically and is computationally more expensive than separable problems. A first attempt to propose new diagnoses for g modes based on 2D computations of modes has been made by Ouazzani et al. (2017).

To simplify the problem of computing eigenmodes for moderate rotators, the so-called traditional approximation (Eckart 1960) can be used. It consists in neglecting the horizontal component of the rotation vector, which makes the problem separable again if the star is assumed to rotate uniformly and if the centrifugal deformationis not taken into account. Within this approximation, the effect of the rotation on low-frequency g modes is greatly simplified (Lee & Saio 1997; Townsend 2003; Bouabid et al. 2013) and the approximation has been used to interpret seismic data of γ Doradus (Van Reeth et al. 2016) and SPB stars (Pápics et al. 2017).

The effect of differential rotation on the waves and on the transport they generate has been studied by various authors (Ando 1985; Lee & Saio 1993; Dzhalilov & Staude 2004; Mathis 2009) but most of them focus on a given type of differential rotation profile (e.g. cylindrical or radial), or use simplifying assumptions that are not completely justified, such as the traditional approximation.

To progress in the understanding of the effect of a general differential rotation on stellar oscillations, it is also possible to build asymptotic theories, which make approximations to gain physical insight. Studies of gravito-inertial waves based on characteristics have been done for uniform (Friedlander & Siegmann 1982; Dintrans & Rieutord 2000) and differential Mirouh et al. (2016) rotation. In the case of purely inertial waves (without stratification), radial and cylindrical differential rotation profiles (Baruteau & Rieutord 2013) and conical ones (Guenel et al. 2016) have been considered.

Contrary to the method of characteristics, the small-wavelength Jeffrey-Wentzel-Kramers-Brillouin (JWKB) asymptotic theory can take into account the compressibility effects that produce the back-refraction of the waves approaching the stellar surface. This approach has been followed by Lignières & Georgeot (2009) to study acoustic modes in rapidly rotating stars. Using semi-classical quantisation concepts to link acoustic rays to pressure modes, they predicted spectral patterns that have been successfully confronted to bi-dimensional numerical computations of modes (see also Pasek et al. 2011, 2012).

More recently, a JWKB asymptotic theory for g modes in a uniformly rotating star was built (Prat et al. 2016, hereafter referred to as Paper I). For the first time, this theory took into account the full effect of rotation (both Coriolis and centrifugal accelerations) and a realistic back-refraction of gravito-inertial waves near the stellar surface. The main prediction of this paper is that modes can be classified as either (i) regular modes, which are similar to modes found in the non-rotating case; (ii) island modes, where the energy is localised around a so-called periodic orbit; or (iii) chaotic modes, which are expected to be characterised by an irregular spatial pattern and a generic statistical distribution of the frequency spacings between nearest modes. Using results of Paper I, Prat et al. (2017) derived theoretical period spacings in the low-frequency regime, where regular modes are expected to dominate.

The aim of the current paper is to investigate the impact of general differential rotation on axisymmetric gravito-inertial waves. First, we derive an eikonal equation (i.e. a local dispersion relation) for axisymmetric gravito-inertial waves in a differentially rotating star (Sect. 2). Second, we describe the ray model associated with the eikonal equation and the tools used to numerically investigate the nature of the ray dynamics (Sect. 3). Third, we explore the ray dynamics for various types of differential rotation profiles in Sect. 4. Finally, we conclude in Sect. 5.

Fig. 1 Coordinate systems.

2 Eikonal equation for gravito-inertial waves

2.1 Baroclinic, rotating equilibrium model

We consider here a general background model of star that takes baroclinicity into account. The equilibrium equations are $$\begin{array}{l}-{\rho }_{0}s{\Omega }^2{e}_s=-\nabla P_0-{\rho }_0\nabla {\psi }_0,\\ {\nabla }^2{\psi }_0=4\pi G{\rho }_0,\end{array}$$

where ρ₀, P₀, and ψ₀ are the equilibrium density, pressure and gravitational potential, respectively; s is the cylindrical radial coordinate, e_s the unit vector associated with it, Ω is the rotation rate, and G the gravitational constant. Equation of state and an entropy equation have to be added to these equations to close the system. The effects of an induced meridional circulation are neglected. In practice, Eqs. (1) and (2) are sufficient to specify the background terms that are involved in the dispersion relation of waves.

Equation (1) can be written in the form $$0=-\nabla P_0+{\rho }_0{g}_0,$$(3)

where g₀ is the effective gravity defined by $$g_0=-\nabla {\psi }_0+s{\Omega }^2{e}_s.$$(4)

In the barotropic case, pressure is a function of density only, so that g₀ can be defined as a gradient, as done in Paper I; in the baroclinic case, g₀ is no longer a gradient.

2.2 Perturbation equations

We use here two major approximations: the adiabatic approximation, which neglects dissipative processes so that the evolution of the fluid is isentropic, and the Cowling approximation, which neglects the variations of the gravitational potential induced by waves (Cowling 1941). The equations for small perturbations around the equilibrium state in the inertial frame are $$\begin{array}{l}({\partial }_t+\Omega {\partial }_{\phi })\rho +\nabla \cdot ({\rho }_{0}v)=0,\\ ({\partial }_t+\Omega {\partial }_{\phi })v_i+{(f\wedge v)}_i+(v\cdot Q){\delta }_{\phi i}=-\frac{{\nabla }_{i}P}{{\rho }_0}+\frac{\rho }{{\rho }_0}{g}_{0i},\\ c_{\text{s}}{}^2\left[({\partial }_t+\Omega {\partial }_{\phi })\rho +v\cdot \nabla {\rho }_0\right]=({\partial }_t+\Omega {\partial }_{\phi })P+v\cdot \nabla P_0,\end{array}$$

where $$\begin{array}{ll}f\hfill & =2\Omega e_z,\hfill \\ Q\hfill & =s\nabla \Omega ,\hfill \end{array}$$

ρ, v, and P are density, velocity, and pressure fluctuations, respectively; c_s is the sound speed, and e_z and e_φ are the unit vectors aligned with the rotation axis and in the azimuthal direction (associated with the azimuthal angle φ), respectively.As shown in Appendix A, it is possible to obtain a single equation for pressure fluctuations only of the form $$\mathcal{A}:\nabla \nabla \widehat{P}+B\cdot \nabla \widehat{P}+C\widehat{P}=0,$$(10)

where $\mathcal{A}$ is a rank-2 tensor; the double gradient of a scalar function a is defined by ${(\nabla \nabla a)}_{ij}={\partial }_i{\partial }_{j}a$ in Cartesian coordinates; the symbol: denotes the double contraction of two tensors ($\mathcal{X}:\mathcal{Y}={\sum }_{ij}{\mathcal{X}}^{ij}{\mathcal{Y}}_{ji}$); and $\widehat{P}$ is the complex amplitude of pressure fluctuations, such that $$P(x,t)=\Re \left[\widehat{P}(r,\theta )e^{i(m\phi -\omega t)}\right],$$(11)

where m is the azimuthal order of the waves and ω their angular frequency. Although the problem is initially three-dimensional (3D), for a given m, it becomes 2D due to the quantisation in the azimuthal direction.

2.3 Short-wavelength approximation

After some calculations explained in Appendix B, the wave Eq. (10) is transformed into a local dispersion relation, also called an eikonal equation. This is done by using the JWKB approximation, which consists in searching for wave-like solutions of the form $$\widehat{P}(x)=\Re \left[A(x)e^{i\Phi (x)}\right],$$(12)

where the wavelength associated with the wavevector k = ∇Φ (Φ is the phase) is much smaller than the typical length of variations of the amplitude A associated with the variations of the background model. As discussed in Paper I, this would normally be equivalent to retaining only second-order derivatives, with ${(\nabla \nabla \widehat{P})}_{ij}\simeq -\widehat{P}{k}_i{k}_j$. However, to ensure that waves are back-refracted near the stellar surface (which is needed to study modes) we also retain a zeroth-orderterm that becomes large near the surface. This means that the JWKB approximation is valid only in an intermediate regime of wavelength, where the latter has to be larger than the inverse of the density scaleheight. Moreover, since we focus here on the asymptotic regime of low-frequency gravity waves, we can neglect as a first step the acoustic part of the wave equation.

There is a priori no reason to assume that m is large enough so that terms in m² or mk should also be retained in the eikonal equation. The latter then reads $$\begin{array}{ll}(k^2+k_{\text{c}}{}^2){\widehat{\omega }}^2=& f(f+Q_s)k_z{}^2+N_0{}^2{k}_{\perp }{}^2-f{Q}_z(k_s{k}_z+k_{\parallel }{k}_{\perp })\\ & +f\cos \Theta (f\cos \Theta +Q_{\perp })k_{\text{c}}{}^2,\end{array}$$(13)

where $\widehat{\omega }=\omega -m\Omega$ is the Doppler-shifted angular frequency, ∥ and ⊥ refer to components along e_∥ and e_⊥, unit vectors parallel and orthogonal to the effective gravity, and Θ is the angle between the rotation axis and e_∥, as illustrated in Fig. 1; and k_c, defined in Eq. (B.15), is a term accounting for the back-refraction of waves near the stellar surface. The only dependence on m is included in the Doppler-shifted frequency.

It is possible, however, to obtain a more complete eikonal equation by considering that the aforementioned terms in m² or mk should be retained as well because they come from second-order derivatives (including those involving φ). Subsequently, the eikonal equation becomes $$\begin{array}{l}i{\widehat{\omega }}^3(k^2+k_{\text{c}}{}^2)+{\widehat{\omega }}^2{k}_{\phi }Q\cdot k\\ -i\widehat{\omega }[f(f+Q_s)k_z{}^2+N_0{}^2(k_{\perp }{}^2+k_{\phi }{}^2)\\ -f{Q}_z(k_s{k}_z+k_{\parallel }{k}_{\perp })+f\cos \Theta (f\cos \Theta +Q_{\perp })k_{\text{c}}{}^2]\\ +k_{\phi }\left[(f{Q}_z{Q}_{\parallel }-N_0{}^2{Q}_{\perp })k_{\perp }-f^2{Q}_z{k}_z\right]=0,\end{array}$$(14)

where k_φ = m∕(rsinθ). In the general case, this equation canbe solved analytically for $\widehat{\omega }$, but the solutions are non-trivial functions of structure, rotation, and k.

When there is no differential rotation, we recover the eikonal equation obtained in Paper I: $${\widehat{\omega }}^2=\frac{f^2{k}_z{}^2+N_0{}^2(k_{\perp }{}^2+k_{\phi }{}^2)+f^2{\cos }^2\Theta k_{\text{c}}{}^2}{k^2+k_{\text{c}}{}^2}.$$(15)

When only axisymmetric waves are considered (k_φ = 0, and $\widehat{\omega }=\omega$), Eq. (14) reduces to $$\begin{array}{ll}(k^2+k_{\text{c}}{}^2){\omega }^2=& f(f+Q_s)k_z{}^2+N_0{}^2{k}_{\perp }{}^2-f{Q}_z(k_s{k}_z+k_{\parallel }{k}_{\perp })\\ & +f\cos \Theta (f\cos \Theta +Q_{\perp })k_{\text{c}}{}^2.\end{array}$$(16)

When centrifugal deformation (which also includes effects of baroclinicity) is neglected, the eikonal equation in spherical coordinates reads $$\begin{array}{ll}(k^2+k_{\text{c}}{}^2){\omega }^2=\hfill & (k_r{}^2+k_{\text{c}}{}^2)f\cos \theta (f\cos \theta +Q_{\theta })\hfill \\ \hfill & +k_{\theta }{}^2\left[N_0{}^2+f\sin \theta (f\sin \theta +Q_r)\right]\hfill \\ \hfill & -k_r{k}_{\theta }[f\cos \theta (f\sin \theta +Q_r)\hfill \\ \hfill & +f\sin \theta (f\cos \theta +Q_{\theta })].\hfill \end{array}$$(17)

From now on, we focus our study on axisymmetric waves that are described by the eikonal Eq. (16).

2.4 Domains of propagation

The eikonal Eq. (16) can be seen as a quadratic equation in k_∥ or k_⊥. It is thus possible to derive a condition for them to be real: $$\Gamma k_{\perp }{}^2-{\left[{\omega }^2-f\cos \Theta (f\cos \Theta +Q_{\perp })\right]}^2{k}_{\text{c}}{}^2\ge 0,$$(18)

where $$\begin{array}{cc}\Gamma \hfill & =f\cos \Theta \left[f{Q}_z(f\sin \Theta +Q_{\parallel })-N_0{}^2(f\cos \Theta +Q_{\perp })\right]\hfill \\ \hfill & +{\omega }^2\left[N_0{}^2+f(f+Q_s)\right]-{\omega }^4.\hfill \end{array}$$(19)

In the bulk of the star, k_c is negligible, so the propagation condition reduces to Γ ≥ 0. In the following, we focus on the region where this is the case. Γ can be seen as a quadratic function of ω² that always have two real solutions, ${\omega }_{-}{}^2$ and ${\omega }_{+}{}^2$, as proved inAppendix C. Thus, the propagation condition is equivalent to $${\omega }_{-}{}^2\le {\omega }^2\le {\omega }_{+}{}^2.$$(20)

In the particular case of a purely radial differential rotation and when the background model is assumed to be spherically symmetric, one can write the propagation condition (based on the eikonal equation given in Eq. (17), which neglects all effects of centrifugal deformation, including baroclinic effects) in the form $$A(r){\cos }^4\theta +B(r){\cos }^2\theta +C(r)\ge 0,$$(21)

with $$\begin{array}{ll}A(r)\hfill & =-{\Omega }^2{\widehat{Q}}_r{}^2,\hfill \\ B(r)\hfill & ={\Omega }^2{\widehat{Q}}_r{}^2-f^2{N}_0{}^2-{\omega }^{2}f{\widehat{Q}}_r,\hfill \\ C(r)\hfill & ={\omega }^2\left[N_0{}^2+f(f+{\widehat{Q}}_r)-{\omega }^2\right],\hfill \end{array}$$

where ${\widehat{Q}}_r=r{\partial }_r\Omega$. Roots of Eq. (21) correspond to critical co-latitudes, which delimit regions where waves can propagate. Thus, when a critical co-latitude θ_c exists, it satisfies $${\theta }_{\text{c}}=\arccos \sqrt{\frac{-B\pm \sqrt{B^2-4AC}}{2A}},$$(25)

which depends in general on the radius.

In the case of general differential rotation but in regions where N₀ is much larger than f and Q (the norm of Q), the propagation condition Eq. (20) simplifies to $$f\cos \Theta (f\cos \Theta +Q_{\perp })\le {\omega }^2\le N_0{}^2+f\sin \Theta (f\sin \Theta +Q_{\parallel }).$$(26)

The left inequality means that critical co-latitudes Θ_c may exist for which ω² = fcosΘ_c(fcosΘ_c + Q_⊥). In solid-body rotation, they exist if and only if ω < f. Here, since f and Q_⊥ vary with space, the picture is much more complex, and many different situations may occur, with either none, one, or several critical surfaces. By definition, Q_⊥ vanishes along the rotation axis. Further, if the rotation is symmetric with respect to the equator, Q_⊥ also vanishes at the equator. This means that the lower part of the condition of propagation Eq. (26) becomes 0 ≤ ω at the equator (always propagation) and f ≤ ω at the poles. As a consequence, if ω < f(θ = 0), waves cannot propagate near the rotation axis, and there is at least one critical surface in latitude between θ = 0 and θ = π∕2. When critical pseudo co-latitudes exist, they verify $$\cos {\Theta }_{\text{c}}=\frac{-Q_{\perp }\pm \sqrt{Q_{\perp }{}^2+4{\omega }^2}}{2f},$$(27)

which is possible only when ω² ≤ f(f + Q_⊥) if Θ ≤ π∕2, or ω² ≤f(f − Q_⊥) if Θ ≥ π∕2. Thus, if $\omega >{\max }_{\theta }[f(f+\left|Q_{\perp }\right|)]$, there is no criticalsurface in latitude.

In contrast, in regions where f is much larger than Q and N₀, Eq. (20) becomes $$N_0{}^2{\cos }^2\Theta -A\le {\omega }^2\le f(f+Q_s)+N_0{}^2{\sin }^2\Theta +A,$$(28)

where $A=f{Q}_z\sin \Theta \cos \Theta +Q_z{}^2{\cos }^4\Theta$. By definition, Q vanishes at the centre. As N₀ also vanishes, the right-hand inequality means that the waves cannot propagate to the centre when ω >f.

To summarise the above discussion on the propagation domains, we find that waves that are sub-inertial throughout (ω < f) avoid a region around the rotation axis while waves that are super-inertial throughout (ω > f) avoid the stellar centre. This property is verified in differentially rotating stars as in uniformly rotating ones, but new features appear when waves are not sub- or super-inertial throughout. For instance, propagation domains with multiple avoided regions in latitude are only possible in differentially rotating stars. Examples of these new types of domains are illustrated in Sect. 4.

Mirouh et al. (2016) have also studied propagation domains of gravito-inertial waves. However, they considered a Boussinesq fluid between two concentric rigid spheres for one particular radial differential rotation, which makes the comparison with the present case not very instructive.

3 Ray model for axisymmetric waves

3.1 Hamiltonian formalism

The angular frequency of waves ω is constant. Thus, the scalar function ω(x, k) must remain constant when waves propagate and be equal to ω(x₀, k₀), which is set by the initial condition (x₀, k₀). This implies that the propagation of waves can be described with the Hamiltonian formalism.

We define a ray as a trajectory tangent to the group velocity v_g = ∇_kω, where ∇_k is the gradient with respectto k. It reads $$\frac{\text{d}x}{\text{d}t}={\nabla }_k\omega .$$(29)

This choice is motivated by the fact that the group velocity characterises the transport of energy by the waves, whereas the phase velocity v_p = ωk∕k² is the velocity of the wave front. The constancy of ω then requires the evolution of the wavevector along the ray path to be governed by $$\frac{\text{d}k}{\text{d}t}=-{\nabla }_x\omega ,$$(30)

where ∇_x is the spatial gradient. Equations (29) and (30) thus have a Hamiltonian form, where the Hamiltonian is ω.

In Paper I, it was shown that the ray dynamics can be written in spherical coordinates (r, θ) as $$\begin{array}{ll}\frac{\text{d}r}{\text{d}t}\hfill & =\frac{\partial \omega }{\partial k_r},\hfill \\ \frac{\text{d}\theta }{\text{d}t}\hfill & =\frac{1}{r}\frac{\partial \omega }{\partial k_{\theta }},\hfill \\ \frac{\text{d}{k}_r}{\text{d}t}\hfill & =-\frac{\partial \omega }{\partial r}+\frac{k_{\theta }}{r}\frac{\partial \omega }{\partial k_{\theta }},\hfill \\ \frac{\text{d}{k}_{\theta }}{\text{d}t}\hfill & =-\frac{1}{r}\frac{\partial \omega }{\partial \theta }-\frac{k_{\theta }}{r}\frac{\partial \omega }{\partial k_r}.\hfill \end{array}$$

The presence of the last term in each of the last two equations comes from the fact that the basis used for the wavevector (k = k_re_r + k_θe_θ) is different from the natural one associated with spherical coordinates ($k=k_r{e}_r+k_{\theta }^{\text{nat}}{e}_{\theta }/r$). These equations are made explicit in Appendix D.

Fig. 2 Rotation profile with ΩC∕ΩR = 2, rC = 0.2, and αC= 25 (blue solid line). The black dashed lines are the three frequencies for which we computed PSS: from bottom to top ω∕(2ΩR) = 0.8, 1.6, and 2.4. They correspond to the sub-, trans-, and super-inertial regimes, respectively.

3.2 Phase-space visualisation: the Poincaré surface of section

Rays can bestudied as dynamical systems with two degrees of freedom, so the phase space is four-dimensional (r, θ, k_r, k_θ). At a given frequency, all trajectories stay on a 3Dspace because ω is constant in time. To visualise the structure of the phase space, it is convenient to consider the intersection of the latter with a given hyper-surface, usually defined by fixing one phase-space coordinate. This is called a Poincaré section or surface of section (PSS). To have a complete view of the structure of the phase space, the intersecting hyper-surface must be chosen so that most trajectories intersect it several times. For our problem, the equatorial plane (θ = π∕2) appears as a logical choice, because most of the ray trajectories cross it, including low-frequency waves, which are trapped near the equatorial plane. It is important to retain points for the PSS only when the intersecting surface is crossed from a given side (in our case, from the region where θ < π∕2). This ensures that two trajectories never cross on a PSS.

In the following, PSS are represented in the plane (r, k_r). For consistency with Paper I, the radius is normalised by the equatorial radius of the star R_e, while the radial wavevector is normalised by the quantity $N_{0,\text{max}}^{\text{e}}/(\omega R_{\text{e}})$, where $N_{0,\text{max}}^{\text{e}}$ is the value of the inner maximum of the radial profile of the Brunt-Väisälä frequency in the equatorial plane (see Paper I for more detail, in particular the simple analytical description of the PSS in the non-rotating case).

3.3 Numerical method

To properly characterise the dynamics of gravito-inertial rays of a system with differential rotation and a non-trivial background state such as those considered here, it is necessary to follow a sufficient number of rays of varying frequency. In particular, to identify the different regions of the phase space, invariant tori and chaotic regions, each of those rays needs to be well-sampled both in time and space.

The spatial accuracy of the ray dynamics is limited by the method employed to compute the background medium, that is the model of rapidly rotating star. For the simulations carried out here, the background thermodynamic state is computed assuming a polytropic equation of state and uniform rotation, thus taking into account the rotationally modified gravitational potential. This state is expressed as the coefficients of a truncated spectral expansion in radius as Chebyshev polynomials and in latitude as Legendre polynomials, thereby limiting the accuracy of its spatial reconstruction. Changes in the values of the background state due to the spectral truncation error are very small in absence of discontinuity in the background state.

Generating a statistically significant sample of intersections with a given PSS requires a numerical integrator that is both stable and accurate enough to follow the ray for a sufficiently long time. This integrator must also be robust in regions where the equations become stiff, such as near the coordinate origin, as well as near turning surfaces like the stellar surface. Fortuitously, the ray dynamics equations employed here have a symplectic or Hamiltonian character, as described in Sect. 3.1. Therefore, we may appeal to an extant class of implicit symplectic integrators that have the property that the simplectic structure is preserved under the discrete map of the numerical method. This means that the volume of the phase space for each ray is also preserved.

One particular set of implicit symplectic integrators are the Gauss-Legendre-Runge-Kutta (GLRK) methods. To make this explicit for ray dynamics, one can formulate the problem as a coupled set of ordinary differential equations d y∕dt = f(y) with the initial condition y(0) = y₀, which also selects the frequency of the ray, where y is the tuple (x, k) and f(y) = (∇_kω, −∇_xω). In a discretised form, the s-stage implicit GLRK method yields $$y_{n+1}=y_n+h{\sum }_{j=1}^s{b}_{j}f({\xi }_j),$$(35)

where h is the fixed step size. The ξ_j are given implicitly by $${\xi }_j-y_n-h{\sum }_{i=1}^s{a}_{ji}f({\xi }_i)=0.$$(36)

The computation of the coefficients a_ij and b_j follows froma set of algebraic equations that themselves are derived from the coefficients of the Taylor expansion of ξ_j and f (Butcher 1963). Even so, there are unconstrained parameters that allow for the construction of a variety of RK methods with various error bounds. Here the c_j = ∑_ia_ji are those free parameters and are chosen to be the zeros of the shifted Legendre polynomial of degree s. This choice leads to the symplectic nature of the method in a non-trivial way (e.g. Sanz-Serna 1988). We have implemented this method with a choice of stage between s = 1 and 5, as the Butcher tableau for these methods are widely available, where the symmetry of the coefficients implies that these methods have a truncation error of $\mathcal{O}(h^{2s})$.

Since the ξ_j are implicit and f is a non-linear function, we have implemented a Newton-Raphson solver to compute them. To do this, we construct a new vector function X(y_n, ξ_j) that concatenates Eq. (36) as $$\begin{array}{l}{X}_{s(j-1)+k}^q={\xi }_{jk}^q-y_{nk}-h{\sum }_{i=1}^s{a}_{ji}{f}_k({\xi }_i^q),\\ J^q\left(y_n,{\xi }_j^q\right)\delta X^q=-X^q,\\ X^{q+1}=X^q+\delta X^q.\end{array}$$

For clarity in identifying the component index of X, we note thats is the total number of stages, j is the index of the current stage, and k is the index of the component of the tuple y. Furthermore, the index q denotes the current step of the iterative solver and J^q is the Jacobian of the function X evaluated aty_n and ${\xi }_j^q$. To reduce the cost of its computation, the derivatives of X used to formthe Jacobian are finite differences, rather than analytical derivatives. The iterative implicit solver necessarily places greater restrictions on the step size h that follow from the properties of f and from the requirement of uniqueness of the solutions ξ_j. Hence, these restrictions will also vary depending upon the initial conditions of each ray, which means that choosing an appropriate h is a nontrivial exercise that often has to be done heuristically. Moreover, enforcing global uniqueness typically requires h to be chosen to be smaller than one would otherwise expect.

We have further implemented an adaptive time stepping method to partially circumvent the time step restrictions. To do this, we solve a modified Hamiltonian problem following the method described in Hairer (1997), but one that retains the symplectic nature of the original equations. The idea behind this is that one has a step size h = χ(y)Δt, which corresponds to a remapping of the time variable as τ = χ(y)t. Since the bulk of the time stepping issues occur near the coordinate origin, we have chosen χ(y) = χ(r) ∝ 1 − e^−λr with r being the radius and λ a constant. A modified Hamiltonian is constructed that depends upon χ as $\mathcal{K}=\chi (r)(\mathcal{H}-{\mathcal{H}}_0)$, where $\mathcal{H}=\omega (y)$ and ${\mathcal{H}}_0=\omega (y_0)$. This implies that the symplecticity-preserving equations of motion become $$\begin{array}{ll}\frac{\text{d}x}{\text{dt}}\hfill & =\chi \left(r\right){\nabla }_k\mathcal{H},\hfill \\ \frac{\text{d}k}{\text{dt}}\hfill & =-\chi \left(r\right){\nabla }_x\mathcal{H}-\frac{\partial \chi \left(r\right)}{\partial r}\left(\mathcal{H}-{\mathcal{H}}_0\right).\hfill \end{array}$$

Therefore, the only notable change in the method to accommodate these time step changes is in the right-hand-side vector function f of the GLRK method.

This new method has been compared to the fourth-order explicit Runge-Kutta method used in Paper I and is more accurate. However, it comes at the expense of more function evaluations, and thus a higher computational cost.

Fig. 3 PSS (left) and examples of ray trajectories (right) for the three regimes identified in Fig. 2: sub-, trans-, and super-inertial, respectively. Blue trajectories correspond to invariant tori (except for the second case, which corresponds to an island chain), green ones to island chains, and the red one to a chaotic zone. Lighter blue and green trajectories correspond to periodic orbits at the centre of the islands. Magenta lines are the limits of the domain of propagation. The imprints of the rays are shown on the PSS with colours corresponding to the rays.

4 Gravito-inertial ray dynamics in differentially rotating stars

For simplicity, the background stellar structure we used for our numerical computations is a centrifugally deformed polytropic (so barotropic) model of a uniformly rotating star, and the rotation profile is given by a prescribed analytical formula presented hereafter. This is not fully self-consistent because (i) the centrifugal deformation is slightly incorrect, and (ii) the rotation profiles considered here would lead to a baroclinic structure. This may for example affect the profile of the Brunt-Väisälä frequency. However, this simplification allows us to easily investigate the various types of general differential rotation.

In the following we present the general law of differential rotation we use (Sect. 4.1) and then study the ray dynamics, considering the two cases of radial (Sect. 4.2) and latitudinal (Sect. 4.3) differential rotation separately.

4.1 Rotation profile

For the computations shown in the present paper, we considered the following three-zones rotation profile: $${\Omega }_0(r,\theta )=\frac{{\Omega }_{\text{C}}{e}^{-{\alpha }_{\text{C}}(r-r_{\text{C}})}+{\Omega }_{\text{R}}+[{\Omega }_{\text{R}}+{\Omega }_{\text{D}}\cos (2\theta )]e^{{\alpha }_{\text{E}}(r-r_{\text{E}})}}{e^{-{\alpha }_{\text{C}}(r-r_{\text{C}})}+1+e^{{\alpha }_{\text{E}}(r-r_{\text{E}})}}.$$(42)

The core is characterised by an almost uniform rotation rate Ω_C up to a radius r_C, and α_C is the steepness of the transition with the bulk of the radiative zone. The latter rotates at a uniform rate Ω_R. The envelope starts at a radius r_E and has a mean rotation rate Ω_R and a latitudinal differential rotation characterised by Ω_D. The steepness of the transition between the envelope and the bulk of the radiative zone is given by α_E. In solar-type stars, the envelope is convective, but here, since we use a polytropic model with a single polytropic index, the full star is radiative. Ω_D < 0 corresponds to a solar differential rotation, with the equator faster than the poles, whereas Ω_D > 0 corresponds to an anti-solar differential rotation. Both configurations are supported by 3D numerical simulations (see Brun & Toomre 2002; Brown et al. 2008; Matt et al. 2011; Gastine et al. 2014; Brun et al. 2017).

At the centre of the star, this profile may not be well defined. First, if Ω_D ≠0, it gives different values of Ω₀ for different values of θ. Second, the presence of differential rotation causes the radial gradient of rotation to be non-zero at the centre. To eliminate these pathological features, we regularise the rotation profile, adding the following corrections to the original profile: $$\Omega (r,\theta )={\Omega }_0(r,\theta )-\gamma e^{-{\alpha }_{\text{E}}{r}^2/2+\delta r}\cos (2\theta )+\epsilon \frac{e^{-{\alpha }_{\text{C}}r}}{{\alpha }_{\text{C}}},$$(43)

where $$\begin{array}{ll}\gamma \hfill & ={\Omega }_{\text{D}}\frac{e^{-{\alpha }_{\text{E}}{r}_{\text{E}}}}{e^{{\alpha }_{\text{C}}{r}_{\text{C}}}+1+e^{-{\alpha }_{\text{E}}{r}_{\text{E}}}},\hfill \\ \delta \hfill & =\frac{{\alpha }_{\text{E}}+({\alpha }_{\text{C}}+{\alpha }_{\text{E}})e^{{\alpha }_{\text{C}}{r}_{\text{C}}}}{e^{{\alpha }_{\text{C}}{r}_{\text{C}}}+1+e^{-{\alpha }_{\text{E}}{r}_{\text{E}}}},\hfill \\ \epsilon \hfill & =({\Omega }_{\text{R}}-{\Omega }_{\text{C}})\frac{{\alpha }_{\text{C}}{e}^{{\alpha }_{\text{C}}{r}_{\text{C}}}+({\alpha }_{\text{C}}+{\alpha }_{\text{E}})e^{{\alpha }_{\text{C}}{r}_{\text{C}}-{\alpha }_{\text{E}}{r}_{\text{E}}}}{{(e^{{\alpha }_{\text{C}}{r}_{\text{C}}}+1+e^{-{\alpha }_{\text{E}}{r}_{\text{E}}})}^2}.\hfill \end{array}$$

We need to consider the stability of the rotation profile with respect to hydrodynamical instabilities (e.g. Zahn 1983; Maeder 2009). In particular, we made sure that the profiles used in the present paper are Rayleigh-Taylor stable, that is, the quantity $$N_{\Omega }{}^2=f(f+Q_s)$$(47)

is positive in the whole star. We also verified that the rotation profile is never centrifugally unstable. In a spherical model, this means that the quantity rsin²θΩ(Ω + r∂_rΩ) must always be smaller than gravity. At the equator, the radial differential rotation is negligible, and the previous stability criterion reduces to Ω < Ω_K, where $${\Omega }_{\text{K}}=\sqrt{\frac{GM}{R_{\text{eq}}{}^3}},$$(48)

is the critical angular velocity and M is the stellarmass.

4.2 Radial differential rotation near the core

In this section we investigate the influence of radial differential rotation near the core on the structure of the phase space. This is motivated by a significant number of observations showing a contrast in the rotation rate between the core and the envelope of stars (see Sect. 1).

We takeΩ_D = 0 and Ω_R∕Ω_K = 0.38. The two cases (fast or slow core) are studied in Sects. 4.2.1 and 4.2.2, respectively.

4.2.1 Fast core

First, we consider the case Ω_C∕Ω_R = 2. We thus have Ω_R < Ω < Ω_C, which defines three different frequency regimes. To understand these regimes, we computed PSS at the three different frequencies shown in Fig. 2.

First, when ω < 2Ω_R, ω < 2Ω throughout the star, so according to Sect. 2.4, rays avoid a region around the rotation axis and can propagate near the centre, as for sub-inertial rays in the uniformly rotating case. We therefore call this the purely sub-inertial regime. As illustrated in Fig. 3a, the structure of the phase space is dominated by invariant tori, as it is in the corresponding regime with uniform rotation. Invariant tori are the only kind of structures present in integrable systems.

Second, when ω > 2Ω_C, ω > 2Ω everywhere in the star. It follows that rays can propagate everywhere but near the centre, as in the uniformly rotating super-inertial case. For that reason, we call it the purely super-inertial regime. Figure 3c shows an example of PSS computed in this regime featuring the same kind of structure that was found in Paper I: mainly invariant tori and island chains, such as those associated with rosette modes, which have been identified by Ballot et al. (2012) and further described by Takata & Saio (2013); Saio & Takata (2014); Takata (2014).

Third, when 2Ω_R < ω < 2Ω_C, ω < 2Ω in most of the core, and ω > 2Ω in most of the radiative zone. As a consequence, rays can propagate near the centre, and can propagate at any latitude in the radiative zone. This new regime, that we call trans-inertial, is dominated by chaotic regions and a few island chains, as can be seen in Fig. 3b. The presence of chaos may be explained by the fact that near the transition between the sub- and super-inertial regions, a small difference in position or momentum results in two completely different behaviours when approaching the rotation axis: either propagation if the ray is locally super-inertial, or reflection on the critical surface if the ray is locally sub-inertial. This is illustrated in Fig. 4, which shows the detail of the red trajectory of Fig. 3b.

To determine if chaos persists in this regime even at low differential rotation, we computed a PSS with Ω_C ∕Ω_R = 1.1 and ω∕(2Ω_R) = 1.05. As can be seen on Fig. 5, chaos is still dominant, with some island chains.

However, it seems to be the result of the juxtaposition of small chaotic zones, whereas we found a large chaotic zone when Ω_C ∕Ω_R = 2. We also computed a PSS with Ω_C∕Ω_R = 4 and ω∕(2Ω_R) = 2.4, which is shown in Fig. 6.

One can see that the dynamics is dominated by chaos, except for a large island chain.

Fig. 4 Detail of a chaotic trans-inertial trajectory showing the transition between the sub-inertial and super-inertial behaviours. The ray bounces several times on the turning surface (sub-inertial, at the top) until it finally propagates past the rotation axis (super-inertial, at the bottom).

4.2.2 Slow core

We now consider the case Ω_C∕Ω_R = 0.5. Again, there are three different regimes, which are investigated through PSS computed at three frequencies shown in Fig. 7.

In the sub-inertial regime, Fig. 8a shows that the dynamics is dominated by invariant tori, as in the case of a fast core. This suggests that the near integrability at low frequencies that was observed with uniform rotation is still valid when considering a general radial differential rotation. If this is indeed the case, it means that new seismic diagnoses similar to those obtained by Prat et al. (2017) in the uniformly rotating case could be derived.

The super-inertial regime is also very similar to what we find for the case of the fast core, as can be seen in Fig. 8c. In both sub- and super-inertial regimes, the main difference with the case of the fast core is the shape of theenvelope. As in the case of a fast core, the trans-inertial regime is largely dominated by chaos, with some stability islands. However, the propagation domain, with avoided regions both at the centre (although it is rather small) and around the rotation axis, is different from the one obtained with a fast core. This is illustrated in Fig. 8b.

4.3 Latitudinal differential rotation in the envelope

In this section we investigate the effect of latitudinal differential rotation in the envelope on the structure of the phase space. Such differential rotation takes place in convective envelopes of low-mass stars, as mentioned earlier, but also possibly in radiative envelopes of massive stars (see e.g. Rieutord & Espinosa Lara 2013). Besides, at the interface between the envelope and the bulk of the radiative zone the latitudinal differential rotation of the envelope generates locally a potentially strong radial differential rotation, as it is in the solar tachocline.

For simplicity, we assume a zero radial differential rotation near the core (Ω_C = Ω_R). An example of rotation profile used in this section is given in Fig. 9.

As mentioned in Sect. 2.4, the existence and the number of critical surfaces in latitude depends on the value of f cosΘ(fcosΘ + Q_⊥) with respect to the angular frequency ω. To simplify the discussion, we neglect here the centrifugal deformation. We thus consider the quantity κ defined by κ² = fcosθ(fcosθ + Q_θ), which is similar to a horizontal epicyclic frequency. Figure 10 shows how κ² depends on the latitude and on the degree of differential rotation η = Ω_D∕Ω_R.

First, κ is always zero at the equator.

When − 1∕7 < η < 1∕3, κ² is positive and monotonic, κ is real, and its maximum value κ_max is reached along the rotation axis. When ω < κ_max, there is only one latitude at which ω = κ, so there is only one critical latitude (sub-inertial regime). The smaller the frequency, the closer the critical surface to the equatorial plane, as in the uniformly rotating case. When ω >κ_max, there is no critical surface (super-inertial regime).

When η > 1∕3, κ² is negative near the equatorial plane and becomes positive and monotonic closer to the rotation axis, where it has a maximum value. As in the previous case, we can define sub- and super-inertial regimes. The main difference is that when the frequency tends to zero, the critical surface does not tend to the equatorial plane, but to another limit surface.

When η < −1∕7, κ² is positive, but not monotonic. κ reaches a maximum value at a certain latitude, then decrease and reaches a positive value κ_p along the rotation axis. When ω < κ_p, there is only one critical surface. When κ_p < ω < κ_max, there are two critical surfaces, and waves can propagate near the equatorial plane and the rotation axis. We name this regime mid-inertial. When ω > κ_max, there is no critical surface.

In addition to the number of critical surfaces in the envelope, waves are either sub- or super-inertial in the inner region. There are thus many different regimes, as illustrated in Fig. 11.

In the following, unless mentioned otherwise, we used r_E = 0.7 and α_E = 40. For Figs. 12a and c, which correspond to regimes with a strong anti-solar differential rotation, this would lead to the existence of regions that are unstable with respect to the Rayleigh-Taylor instability. Since for real stars this instability would very rapidly change the background rotation towards a stable configuration, we chose to use α_E = 22, which is a smoother transition, for the concerned calculations to avoid the existence of unstable regions. However, the nature of the dynamics is not affected by this difference.

Regimes 1 and 2 are purely sub-inertial. The main difference between the two regimes is that at low frequency, waves are trapped very close to the equatorial plane in regime 2, and not in regime 1. As shown in Figs. 12a and b, these two regimes are largely dominated by invariant tori and tiny island chains, as expected for purely sub-inertial waves.

Regimes 3 and 4 are both sub-inertial in the envelope and super-inertial in the core, thus trans-inertial. In Figs. 12c and d, one can see that in contrast with regimes 1 and 2, there is not much difference in the domain of propagation, because the frequency is not small compared to the rotation frequency. Both regimes are dominated by invariant tori when the maximum value of k_r on the PSS for a given trajectory is low and by chaotic structures at higher values. This can be explained by the fact that rays with low enough maximum values of k_r do not propagate into the differentially rotating region, and thus behave as in the uniformly rotating case. This was also the case in Sect. 4.2, but it was less visible due to the limited size of the core. The main difference between regimes 3 and 4 is that several island chains of significant size are clearly visible inside the chaotic region in regime 3, but not in regime 4.

Regime 5 is sub-inertial in the core and mid-inertial in the envelope (two critical surfaces). Because of the critical surface in the core, the two regions of the envelope where waves can propagate are not connected to each other. As a consequence, the domain of propagation of waves that are trapped near the equatorial plane is similar to the one found in the purely sub-inertial regimes. However, Fig. 12e shows that in regime 5, invariant tori that dominate the phase space (as in the purely sub-inertial case) coexist with island chains and relatively small chaotic zones. Moreover, the invariant tori seem to be less smooth than in the purely sub-inertial regime. Although waves can also be trapped in the polar regions, such waves are not shown because the PSS based on the equatorial plane cannot capture their dynamics.

Regime 6 is super-inertial in the core and mid-inertial in the envelope. In contrast to regime 5, here the two regions of the envelope where waves can propagate are connected to each other. As illustrated in Fig. 12f, the structure of the phase space is very different depending on whether rays propagate in the envelope or not, as in regimes 3 and 4.

Regime 7 is sub-inertial in the core and super-inertial in the envelope. Although the domain of propagation in this regime is very different for those in regimes 3, 4, and 6, the structure of the phase space is very similar, with different behaviours for rays that propagate into the envelope (and are chaotic) and for those that do not (and are nearly integrable). This can be seen in Fig. 12g.

Regime 8 is purely super-inertial. Figure 12h shows that in this regime, the phase space is dominated by invariant tori and island chains, similarly to the purely super-inertial regime of Sect. 4.2.1

Fig. 5 PSS at ΩC∕ΩR = 1.1 and ω∕(2ΩR) = 1.05 (left) and two trans-inertial rays: one belonging to an island chain (blue) and one belonging to a chaotic zone (red).

Fig. 6 PSS at ΩC∕ΩR = 4 and ω∕(2ΩR) = 2.4 (left) and one trans-inertial ray belonging to an island chain (blue).

Fig. 7 Rotation profile with ΩC∕ΩR = 0.5, rC = 0.2, and αC= 25 (blue solid line). The black dashed lines are the three frequencies for which we computed PSS: from bottom to top ω∕(2ΩR) = 0.4, 0.8, and 1.6.They correspond to the sub-, trans- and super-inertial regimes, respectively.

Fig. 8 PSS (left) and examples of ray trajectories (right) for the three regimes shown in Fig. 7: sub-, trans-, and super-inertial, respectively. Blue trajectories correspond to invariant tori (except for the second case, which corresponds to an island chain), green ones to island chains, and the red one to a chaotic zone. The lighter green trajectory corresponds to a periodic orbits at the centre of the island.

5 Conclusions

In this paper we generalised the work done in Paper I by introducing the effect of a general differential rotation on the ray dynamics for gravito-inertial waves. In contrast to previous studies with differential rotation, we considered here the full Coriolis acceleration (i.e. without the traditional approximation) and the full centrifugal acceleration. Focusing on axisymmetric waves as a first step, we wrote the equations governing the rays dynamics and implemented them in a ray-tracing code. We then numerically investigated the domain of propagation of rays and the nature of their dynamics in various regimes of differential rotation.

We find that differential rotation can generate a large variety of domains of propagation. However, one can distinguish between three main regimes for the nature of the ray dynamics. At low frequency, we observe a regime similar to the sub-inertial regime in the case of uniform rotation, where waves are trapped near the equatorial plane. The dynamics is dominated by invariant tori, even though these structures are deformed with respect to the uniformly rotating case. At high frequency, rays behave as inthe super-inertial regime of the uniformly rotating case, with mostly invariant tori and island chains, and sometimes chaotic zones in a narrow range of frequency (in the vicinity of the inner maximum of the Brunt-Väisälä frequency, see Paper I for more details). Between these two regimes, we find a new regime, called trans-inertial. This regime is characterised by a dynamics dominated by chaotic zones (sometimes with some island chains) for rays that propagate into differentially rotating regions and by invariant tori for rays that stay in regions with negligible differential rotation.

The properties of the modes, which result from the superposition of positively interfering rays, can be inferred from semi-classical quantisation methods and concepts mostly developed in the domain of quantum physics (see for example Sect. 5 of Paper I). Accordingly, the spectrum of axisymmetric gravito-inertial modes in fast rotators should have the following properties: (i) at low frequency, mostly regular modes; (ii) at high frequency, mostly regular modes and island modes, and a few chaotic modes; and (iii) in between, mostly chaotic modes, in a range that depends on the amount of differential rotation. The complete picture should nevertheless be more complicated, since non-axisymmetric modes may behave differently.

The large variety of domains of propagation generated by differential rotation means that gravito-inertial waves can probe many different resonant cavities in differentially rotating stars. The properties of these cavities have a key influence on the visibility of modes and on the transport of angular momentum by waves, which need to extract angular momentum in excitation regions and deposit it in damping regions (see e.g. Mathis 2009).

As already mentioned, the purely super-inertial regime is dominated by invariant tori. This suggests the existence of a nearby integrable system that could provide us with seismic diagnoses for low-frequency gravito-inertial modes, as done in the case of uniform rotation in Prat et al. (2017).

In the present work we did not use self-consistent rotating stellar models, as we neglected the effect of differential rotation on stellar structure. Determining to what extent our results still apply for more realistic background models will require studying the ray dynamics in baroclinic and centrifugally deformed stellar models such as ESTER models (Rieutord et al. 2016). Another possible extension of the present work would be to consider the effect of differential rotation on acoustic rays. As the Coriolis force has a negligible effect on high-frequency acoustic waves, we expect the main effect to come from the way differential rotation modifies the stellar background model. Finally, our predictions on the properties of gravito-inertial modes based on ray dynamics should also be compared with numerically computed modes, using bi-dimensional codes such as TOP (Reese et al. 2006) or ACOR (Ouazzani et al. 2012), knowing that for acoustic modes (Lignières & Georgeot 2009; Pasek et al. 2011, 2012) and some gravito-inertial modes (Ballot et al. 2012), such comparisons have so far been successful.

The ray dynamics can also be used to interpret the huge amount of data produced by global 3D time-dependent simulations of the excitation, propagation, and dissipation of non-linear waves (Alvan et al. 2014, 2015).

In the near future, we plan to generalise the present work to non-axisymmetric waves. This will require a significant amount of work, since the eikonal Eq. (14) we derived for non-axisymmetric waves is more complex than the one for axisymmetric waves, and the Doppler-shifting of waves needs to be taken into account. Since the transport of angular momentum by waves comes from the difference in the excitation and the damping between prograde and retrograde waves, it is crucial to understand the physics of non-axisymmetric waves. Eventually, it will also be necessary to include the effect of dissipative processes in the present formalism to describe the transport generated by the waves.

Another physical ingredient that has been neglected in this study is the magnetic field. In addition to modifying the background structure of the stars, it also affects the propagation of gravito-inertial waves, which become magneto-gravito-inertial waves (Mathis & de Brye 2011), and the transport of angular momentum due to these waves (Mathis & de Brye 2012).

Fig. 9 Rotation profile with ΩD∕ΩR = 0.2, αE = 40, rE = 0.7.

Fig. 10 Value of κ2 = fcosθ(fcosθ + Qθ) as a function of η = ΩD∕ΩR and cos2θ in the envelope. The solid line corresponds to the contour κ2 = 0. The dashed lines denote the position of the extremum of κ2 as a funtion of θ at fixed η. The dotted lines delimits the region between η = −1∕7 and η = 1∕3 where such an extremum does not exist.

Fig. 11 Different regimes as a function of the frequency and the degree of differential rotation in the envelope. See main text for the description of these regimes. The markers correspond to the sets of parameters for which PSS have been computed.

Fig. 12 PSS (left) and examples of ray trajectories (right) for the four first regimes shown in Fig. 11. Blue trajectories correspond to invariant tori, green ones to island chains, and red ones to chaotic zones.

Fig. 12 continued. PSS (left) and examples of ray trajectories (right) for the four last regimes shown in Fig. 11. Blue trajectories correspond to invariant tori (except for the last case, where it corresponds to an island chain, and where the light blue corresponds to the periodic orbit at the centre of the island), green ones to island chains, and red ones to chaotic zones.

Acknowledgements

V.P., S.M., and K.A. acknowledge support from the European Research Council through ERC grant SPIRE 647383. V.P., F.L., and J.B. acknowledge the International Space Science Institute (ISSI) for supporting the SoFAR international team¹. The authorsacknowledge funding by SpaceInn, PNPS (CNRS/INSU), and CNES CoRoT/Kepler and PLATO grants at DAp and IRAP. The authors also thank the referee for giving useful comments that have allowed us to improve the article.

Appendix A Derivation of the wave equation

In this appendix, we want to derive a wave equation for one variable only from the perturbation Eqs. (5)–(7).

First, we assume time-harmonicity, and look for solutions of the form A(x) = ℜ[Â(r, θ)e^i(mφ−ωt)]. Density fluctuations can be expressed in terms of pressure and velocity fluctuations using Eq. (7): $$\widehat{\rho }=\frac{\widehat{P}}{c_{\text{s}}{}^2}-i\frac{\widehat{u}}{\widehat{\omega }}\cdot \left(\frac{\nabla {\rho }_0}{{\rho }_0}-\frac{\nabla P_0}{c_{\text{s}}{}^2{\rho }_0}\right),$$(A.1)

where $\widehat{u}={\rho }_0\widehat{v}$ and $\widehat{\omega }=\omega -m\Omega$. It is possible to link the right-hand side of the dot product to the Brunt-Väisälä frequency defined by $$N_0{}^2=g_0\cdot \left(\frac{\nabla {\rho }_0}{{\rho }_0}-\frac{\nabla P_0}{{\Gamma }_1{P}_0}\right),$$(A.2)

where Γ₁ is the first adiabatic exponent of the fluid. To do so, we use the relation $$\frac{\nabla {\rho }_0}{{\rho }_0}=\left(\frac{\nabla {\rho }_0}{{\rho }_0}\cdot g_0\right)\frac{g_0}{g_0{}^2}+\frac{(\nabla \wedge g_0)\wedge g_0}{g_0{}^2},$$(A.3)

which derives from Eq. (3). Using the definition of the effective gravity given in Eq. (4), Eq. (A.1) becomes $$\widehat{\rho }=\frac{\widehat{P}}{c_{\text{s}}{}^2}+i\frac{\widehat{u}}{\widehat{\omega }{g}_0}\cdot (N_0{}^2{e}_{\parallel }+f{Q}_z{e}_{\perp }),$$(A.4)

where e_∥ = −g₀∕g₀ and e_⊥ is the unit vector orthogonal to e_∥ such that (e_∥, e_⊥, e_φ) is a direct orthonormal basis (see Fig. 1).

The previous equation can now be used to eliminate ρ in Eqs. (5) and (6): $$\begin{array}{l}-i\widehat{\omega }\frac{\widehat{P}}{c_{\text{s}}{}^2}+\frac{\widehat{u}}{g_0}\cdot (N_0{}^2{e}_{\parallel }+f{Q}_z{e}_{\perp })+\nabla \cdot \widehat{u}+\frac{im{u}_{\phi }}{r\sin \theta }=0,\\ \mathcal{M}\cdot \widehat{u}=-\nabla \widehat{P}-\frac{im\widehat{P}{e}_{\phi }}{r\sin \theta }+\frac{\widehat{P}{g}_0}{c_{\text{s}}{}^2},\end{array}$$

where the tensor $\mathcal{M}$ is described in the basis (e_∥, e_⊥, e_φ) by the matrix $$\left[\begin{array}{ccc}-\frac{i}{\widehat{\omega }}({\widehat{\omega }}^2-N_0{}^2)& \frac{i}{\widehat{\omega }}f{Q}_z& -f\sin \Theta \\ 0& -i\widehat{\omega }& -f\cos \Theta \\ f\sin \Theta +Q_{\parallel }& f\cos \Theta +Q_{\perp }& -i\widehat{\omega }\end{array}\right],$$(A.7)

and Θ is the angle between e_z and e_∥, which is equal to the co-latitude θ when the centrifugal deformation is neglected. When $\mathcal{M}$ is invertible, Eq. (A.6) can be used to express û as a function of $\widehat{P}$ and its gradient: $$\widehat{u}={\mathcal{M}}^{-1}\cdot \left(-\nabla \widehat{P}-\frac{im\widehat{P}{e}_{\phi }}{r\sin \theta }+\frac{\widehat{P}{g}_0}{c_{\text{s}}{}^2}\right),$$(A.8)

where $${\mathcal{M}}^{-1}=-\frac{1}{\Gamma }\left[\begin{array}{ccc}i\widehat{\omega }[{\widehat{\omega }}^2-f\cos \Theta (f\cos \Theta +Q_{\perp })]& i\widehat{\omega }f\cos \Theta (f\sin \Theta +Q_{\parallel })& -f({\widehat{\omega }}^2\sin \Theta +f{Q}_z\cos \Theta )\\ i\widehat{\omega }f\cos \Theta (f\sin \Theta +Q_{\parallel })& i\widehat{\omega }[{\widehat{\omega }}^2-N_0{}^2-f\sin \Theta (f\sin \Theta +Q_{\parallel })]& -f\cos \Theta ({\widehat{\omega }}^2-N_0{}^2)\\ {\widehat{\omega }}^2(f\sin \Theta +Q_{\parallel })& ({\widehat{\omega }}^2-N_0{}^2)(f\cos \Theta +Q_{\perp })+f{Q}_z(f\sin \Theta +Q_{\parallel })& i\widehat{\omega }({\widehat{\omega }}^2-N_0{}^2)\end{array}\right].$$(A.9)

and Γ is defined in Eq. (19). Finally, the combination of Eqs. (A.5) and (A.8) leads to a single equation for $\widehat{P}$ only: $$-i\widehat{\omega }\frac{\widehat{P}}{c_{\text{s}}{}^2}+\left(\frac{N_0{}^2{e}_{\parallel }+f{Q}_z{e}_{\perp }}{g_0}+\frac{im{e}_{\phi }}{r\sin \theta }\right)\cdot \left[{\mathcal{M}}^{-1}\cdot \left(-\nabla \widehat{P}-\frac{im\widehat{P}{e}_{\phi }}{r\sin \theta }+\frac{\widehat{P}{g}_0}{c_{\text{s}}{}^2}\right)\right]+\nabla \cdot \left[{\mathcal{M}}^{-1}\cdot \left(-\nabla \widehat{P}-\frac{im\widehat{P}{e}_{\phi }}{r\sin \theta }+\frac{\widehat{P}{g}_0}{c_{\text{s}}{}^2}\right)\right]=0.$$(A.10)

Equivalently, using the tensor identity $\nabla \cdot (\mathcal{A}\cdot b)=(\nabla \cdot \mathcal{A})\cdot b+{\mathcal{A}}^T:\nabla b$, where $\mathcal{A}$ is a tensor, yields the Poincaré equation $$\mathcal{A}:\nabla \nabla \widehat{P}+B\cdot \nabla \widehat{P}+C\widehat{P}=0,$$(A.11)

where $$\begin{array}{ll}\mathcal{A}\hfill & =-{({\mathcal{M}}^{-1})}^T,\hfill \\ B\hfill & =-\left(\frac{N_0{}^2{e}_{\parallel }+f{Q}_z{e}_{\perp }}{g_0}+\frac{im{e}_{\phi }}{r\sin \theta }\right)\cdot {\mathcal{M}}^{-1}-\nabla \cdot ({\mathcal{M}}^{-1})+\left(\frac{g_0}{c_{\text{s}}{}^2}-\frac{im{e}_{\phi }}{r\sin \theta }\right)\cdot {({\mathcal{M}}^{-1})}^T,\hfill \\ C\hfill & =-\frac{i\widehat{\omega }}{c_{\text{s}}{}^2}+\left(\frac{N_0{}^2{e}_{\parallel }+f{Q}_z{e}_{\perp }}{g_0}+\frac{im{e}_{\phi }}{r\sin \theta }\right)\cdot {\mathcal{M}}^{-1}\cdot \left(\frac{g_0}{c_{\text{s}}{}^2}-\frac{im{e}_{\phi }}{r\sin \theta }\right)+\nabla \cdot \left[{\mathcal{M}}^{-1}\cdot \left(\frac{g_0}{c_{\text{s}}{}^2}-\frac{im{e}_{\phi }}{r\sin \theta }\right)\right].\hfill \end{array}$$

Appendix B Determination of the surface term

Close to the surface, c_s becomes very small. Let us consider here the case where N₀ remains finite near the surface. The dominant term of B is $$B_0=\frac{g_0}{c_{\text{s}}{}^2}\cdot {({\mathcal{M}}^{-1})}^T,$$(B.1)

and the one of C is $$C_0=g_0\cdot {({\mathcal{M}}^{-1})}^T\cdot \nabla \left(\frac{1}{c_{\text{s}}{}^2}\right).$$(B.2)

Using the relation $c_{\text{s}}{}^2={\Gamma }_1{P}_0/{\rho }_0$, along with Eqs. (A.2) and (A.3), one can write $$\nabla (c_{\text{s}}{}^2)=\beta g_0+c_{\text{s}}{}^2\left[\frac{\nabla {\Gamma }_1}{{\Gamma }_1}-\frac{(\nabla \wedge g_0)\wedge g_0}{g_0{}^2}\right],$$(B.3)

where β = Γ₁ − 1 − α and $\alpha =c_{\text{s}}{}^2{N}_0{}^2/g_0{}^2$, which tendsto zero at the surface for convective envelopes and to a finite value for polytropic radiative envelopes. Assuming that the term in brackets in the equation above remains finite and that c_s vanishes at the surface implies that near the surface $\nabla (c_{\text{s}}{}^2)\simeq \beta g_0$. Thus, $$C_0\simeq -\beta \frac{g_0}{c_{\text{s}}{}^4}\cdot {({\mathcal{M}}^{-1})}^T\cdot g_0.$$(B.4)

Now, the wave equation Eq. (A.11) is put in the so-called normal form to get rid of first order terms. This is done by writing $\widehat{P}=a\psi$ with a well-chosen function a. If one retains only dominant terms of B and C, and using the fact that the double gradient is symmetric, Eq. (A.11) becomes $$-a\mathcal{N}:\nabla \nabla \psi +(-2\nabla a\cdot \mathcal{N}+a{B}_0)\cdot \nabla \psi +(-\mathcal{N}:\nabla \nabla a+B_0\cdot \nabla a+C_{0}a)\psi =0,$$(B.5)

where $$\mathcal{N}=-\frac{1}{\Gamma }\left[\begin{array}{ccc}i\widehat{\omega }[{\widehat{\omega }}^2-f\cos \Theta (f\cos \Theta +Q_{\perp })]& i\widehat{\omega }f\cos \Theta (f\sin \Theta +Q_{\parallel })& \frac{1}{2}({\widehat{\omega }}^2{Q}_{\parallel }-f^2{Q}_z\cos \Theta )\\ i\widehat{\omega }f\cos \Theta (f\sin \Theta +Q_{\parallel })& i\widehat{\omega }[{\widehat{\omega }}^2-N_0{}^2-f\sin \Theta (f\sin \Theta +Q_{\parallel })]& \frac{1}{2}[({\widehat{\omega }}^2-N_0{}^2)Q_{\perp }+f{Q}_z(f\sin \Theta +Q_{\parallel })]\\ \frac{1}{2}({\widehat{\omega }}^2{Q}_{\parallel }-f^2{Q}_z\cos \Theta )& \frac{1}{2}[({\widehat{\omega }}^2-N_0{}^2)Q_{\perp }+f{Q}_z(f\sin \Theta +Q_{\parallel })]& i\widehat{\omega }({\widehat{\omega }}^2-N_0{}^2)\end{array}\right].$$(B.6)

is the symmetric part of ${\mathcal{M}}^{-1}$. To make all first-order terms vanish, one would need $$2\nabla a\cdot \mathcal{N}=a{B}_0.$$(B.7)

Supposing that this is possible, this would lead to $$\nabla \nabla a\simeq \frac{a}{4c_{\text{s}}{}^4}\left[g_0\cdot {({\mathcal{M}}^{-1})}^T{\mathcal{N}}^{-1}-2\beta g_0\right]\otimes \left[g_0\cdot {({\mathcal{M}}^{-1})}^T{\mathcal{N}}^{-1}\right],$$(B.8)

where ⊗ denotes the outer product of two vectors, which is a tensor defined by ${(a\otimes b)}_{ij}=a_i{b}_j$. The right-hand side of Eq. (B.8) is not symmetric, and thus cannot be the double gradient of a function. This proves that no function a verifies Eq. (B.7). Instead we choose a such that $$2\nabla a\cdot \mathcal{N}=a{B}_{\text{s}},$$(B.9)

where $$B_{\text{s}}=\frac{g_0}{c_{\text{s}}{}^2}\cdot \mathcal{N}.$$(B.10)

This yields $$\nabla a=a\frac{g_0}{2c_{\text{s}}{}^2},$$(B.11)

and Eq. (B.5) becomes $$-\mathcal{N}:\nabla \nabla \psi +B_{\text{a}}\cdot \nabla \psi +k_{\text{c}}{}^2(e_{\parallel }\cdot \mathcal{N}\cdot e_{\parallel })\psi =0,$$(B.12)

where $$B_{\text{a}}=B_0-B_{\text{s}}=\frac{g_0}{c_{\text{s}}{}^2\Gamma }\left[{\widehat{\omega }}^2\left(f\sin \Theta +\frac{Q_{\parallel }}{2}\right)+\frac{f^2{Q}_z\cos \Theta }{2}\right]e_{\phi },$$(B.13)

and $$k_{\text{c}}{}^2=\frac{(1-2\beta )g_0{}^2}{4c_{\text{s}}{}^4}=\frac{(3+2\alpha -2{\Gamma }_1)g_0{}^2}{4c_{\text{s}}{}^4}.$$(B.14)

We note that the term in B_a in Eq. (B.12) will be neglected when applying the JWKB approximation.

The expression for $k_{\text{c}}{}^2$ given in Eq. (B.14) is negative when α = 0 and Γ₁ > 3∕2, which is the case in convective envelopes. However, stars with a convective envelope always have a thin stably stratified surface layer. In the case of a radiative envelope where N₀ becomes very large near the surface and scales as $1/c_{\text{s}}{}^2$, computing the constant term is very lengthy, so we assume here that Eq. (19) of Paper I is still valid in the presence of differential rotation, namely $$k_{\text{c}}{}^2=\frac{[{(1+\alpha )}^2-{\beta }^2]g_0{}^2}{4c_{\text{s}}{}^4}.$$(B.15)

Appendix C Number of solutions of the equation Γ = 0

In the case of uniform rotation, Γ = 0 reduces to $${\omega }^4-{\omega }^2(N_0{}^2+f^2)+N_0{}^2{f}^2{\cos }^2\Theta =0,$$(C.1)

which always have two real solutions for ω².

In presence of differential rotation, some calculations are required. First, the discriminant of Γ = 0 seen as a quadratic equation for ω² reads $$\Delta ={\left[N_0{}^2+f(f+Q_s)\right]}^2+4f\cos \Theta \left[f{Q}_z(f\sin \Theta +Q_{\parallel })-N_0{}^2(f\cos \Theta +Q_{\perp })\right].$$(C.2)

It can be rewritten in terms of Q_∥ and Q_⊥: $$\begin{array}{ll}\Delta & =N_0{}^4+2N_0{}^2\left[f^2(1-2{\cos }^2\Theta )+f(\sin \Theta Q_{\parallel }-\cos \Theta Q_{\perp })\right]+f^4+2f^3\left[\sin \Theta Q_{\parallel }(1+2{\cos }^2\Theta )+\cos \Theta Q_{\perp }(1-2{\sin }^2\Theta )\right]\\ & +f^2\left[Q_{\parallel }{}^2({\sin }^2\Theta +4{\cos }^2\Theta )-2\sin \Theta \cos \Theta Q_{\parallel }{Q}_{\perp }+{\cos }^2\Theta Q_{\perp }{}^2\right].\end{array}$$(C.3)

To go further, we need to consider the reduced discrimant of the equation Δ = 0 seen as a quadratic equation for $N_0{}^2$: $$\begin{array}{ll}{\Delta }^{\prime }& ={\left[f^2(1-2{\cos }^2\Theta )+f(\sin \Theta Q_{\parallel }-\cos \Theta Q_{\perp })\right]}^2-f^4-2f^3\left[\sin \Theta Q_{\parallel }(1+2{\cos }^2\Theta )+\cos \Theta Q_{\perp }(1-2{\sin }^2\Theta )\right]\\ & -f^2\left[Q_{\parallel }{}^2({\sin }^2\Theta +4{\cos }^2\Theta )-2\sin \Theta \cos \Theta Q_{\parallel }{Q}_{\perp }+{\cos }^2\Theta Q_{\perp }{}^2\right].\end{array}$$(C.4)

After simplification, one finally obtains $${\Delta }^{\prime }=-4f^2{\cos }^2\Theta {(f\sin \Theta +Q_{\parallel })}^2,$$(C.5)

which is always negative. As a consequence, the discriminant Δ never changes sign, and remains always positive. Thus, the quadratic equation Γ = 0 always has two real solutions.

Appendix D Ray dynamics equations in spherical coordinates

The first step is to express the general eikonal Eq. (16) for axisymmetric waves in spherical coordinates: $${\omega }^2=\frac{A{k}_r{}^2+2B{k}_r{k}_{\theta }+C{k}_{\theta }{}^2+D{k}_{\text{c}}{}^2}{k_r{}^2+k_{\theta }{}^2+k_{\text{c}}{}^2},$$(D.1)

where $$\begin{array}{ll}A\hfill & =f^2{\cos }^2\theta +f\{-\cos \theta \sin (\theta -\Theta )\cos (\theta -\Theta )Q_r+[\cos \Theta \cos (\theta -\Theta )+\cos \theta {\sin }^2(\theta -\Theta )]Q_{\theta }\}+N_0{}^2{\sin }^2(\theta -\Theta ),\hfill \\ B\hfill & =-f^2\sin \theta \cos \theta -f[\cos \theta {\cos }^2(\theta -\Theta )Q_r+\sin \theta {\sin }^2(\theta -\Theta )Q_{\theta }]+N_0{}^2\sin (\theta -\Theta )\cos (\theta -\Theta ),\hfill \\ C\hfill & =f^2{\sin }^2\theta +f\{[\cos \Theta \sin (\theta -\Theta )+\sin \theta {\cos }^2(\theta -\Theta )]Q_r-\sin \theta \sin (\theta -\Theta )\cos (\theta -\Theta )Q_{\theta }\}+N_0{}^2{\cos }^2(\theta -\Theta ),\hfill \\ D\hfill & =f^2{\cos }^2\Theta +f\cos \Theta [\sin (\theta -\Theta )Q_r+\cos (\theta -\Theta )Q_{\theta }].\hfill \end{array}$$

Equations (31)–(34) then give $$\begin{array}{ll}\frac{\text{d}r}{\text{d}t}\hfill & =\frac{(A-{\omega }^2)k_r+B{k}_{\theta }}{\omega (k^2+k_{\text{c}}{}^2)},\hfill \\ \frac{\text{d}\theta }{\text{d}t}\hfill & =\frac{B{k}_r+(C-{\omega }^2)k_{\theta }}{r\omega (k^2+k_{\text{c}}{}^2)},\hfill \\ \frac{\text{d}{k}_r}{\text{d}t}\hfill & =-\frac{r{\partial }_{r}A{k}_r{}^2+2(r{\partial }_{r}B-B)k_r{k}_{\theta }+[r{\partial }_{r}C-2(C-{\omega }^2)]k_{\theta }{}^2+r{\partial }_{r}D{k}_{\text{c}}{}^2+(D-{\omega }^2)r{\partial }_r(k_{\text{c}}{}^2)}{2r\omega (k^2+k_{\text{c}}{}^2)},\hfill \\ \frac{\text{d}{k}_{\theta }}{\text{d}t}\hfill & =-\frac{{\partial }_{\theta }A{k}_r{}^2+2({\partial }_{\theta }B+A-{\omega }^2)k_r{k}_{\theta }+({\partial }_{\theta }C+2B)k_{\theta }{}^2+{\partial }_{\theta }D{k}_{\text{c}}{}^2+(D-{\omega }^2){\partial }_{\theta }(k_{\text{c}}{}^2)}{2r\omega (k^2+k_{\text{c}}{}^2)}.\hfill \end{array}$$

References