Asymptotic theory of gravity modes in rotating stars II. Impact of general differential rotation

Context. Differential rotation has a strong influence on stellar internal dynamics and evolution, notably by triggering hydrodynamical instabilities, by interacting with the magnetic field, and more generally by inducing transport of angular momentum and chemical elements. Moreover, it modifies the way waves propagate in stellar interiors and thus the frequency spectrum of these waves, the regions they probe, and the transport they generate. Aims. We investigate the impact of a general differential rotation (both in radius and latitude) on the propagation of axisymmetric gravito-inertial waves. Methods. We use a small-wavelength approximation to obtain a local dispersion relation for these waves. We then describe the propagation of waves thanks to a ray model that follows a Hamiltonian formalism. Finally, we numerically probe the properties of these gravito-inertial rays for different regimes of radial and latitudinal differential rotation. Results. We derive a local dispersion relation that includes the effect of a general differential rotation. Subsequently, considering a polytropic stellar model, we observe that differential rotation allows for a large variety of resonant cavities that can be probed by gravito-inertial waves. We identify that for some regimes of frequency and differential rotation, the properties of gravito-inertial rays are similar to those found in the uniformly rotating case. Furthermore, we also find new regimes specific to differential rotation, where the dynamics of rays is chaotic. Conclusions. As a consequence, we expect modes to follow the same trend. Some parts of oscillation spectra corresponding to regimes similar to those of the uniformly rotating case would exhibit regular patterns, while parts corresponding to the new regimes would be mostly constituted of chaotic modes with a spectrum rather characterised by a generic statistical distribution.


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.2R (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 lowfrequency 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; , 2014, 2015Triana 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 solartype 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 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, SPB, β Cephei or Be stars, however, perturbative methods fail (Reese et al. 2006;Ballot et al. 2010Ballot et al. , 2013 and a more complete treatment of rotation is needed. The eigenvalue problem is fully twodimensional (2D), in the sense that it cannot be separated into one-dimensional problems as in the non-rotating case (Rieutord 2009). In general, the problem is 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 deformation is 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 smallwavelength 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. 2011Pasek et al. , 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.

Baroclinic, rotating equilibrium model
We consider here a general background model of star that takes baroclinicity into account. The equilibrium equations are where ρ 0 , P 0 , and ψ 0 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.
(4) In the barotropic case, pressure is a function of density only, so that g 0 can be defined as a gradient, as done in Paper I; in the baroclinic case, g 0 is no longer a gradient.

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 Article number, page 2 of 17 V. Prat et al.: Asymptotic theory of gravity modes in rotating stars. II. Impact of general differential rotation where f = 2Ωe z , ρ, 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 where A is a rank-2 tensor; the double gradient of a scalar function a is defined by (∇∇a) i j = ∂ i ∂ j a in Cartesian coordinates; the symbol : denotes the double contraction of two tensors (X : Y = i j X i j Y ji ); andP is the complex amplitude of pressure fluctuations, such that 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.

Short-wavelength approximation
After some calculations explained in Appendix B, the wave equation (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 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 (∇∇P) i j −Pk 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-order term 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 lowfrequency 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 2 or mk should also be retained in the eikonal equation. The latter then reads whereω = ω − mΩ 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 2 or mk should be retained as well because they come from secondorder derivatives (including those involving ϕ). Subsequently, the eikonal equation becomes where k ϕ = m/(r sin θ). In the general case, this equation can be solved analytically forω, 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: When only axisymmetric waves are considered (k ϕ = 0, andω = ω), Eq. (14) reduces to When centrifugal deformation (which also includes effects of baroclinicity) is neglected, the eikonal equation in spherical coordinates reads From now on, we focus our study on axisymmetric waves that are described by the eikonal equation (16).

Domains of propagation
The eikonal equation (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: where 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 ω 2 that always have two real solutions, ω − 2 and ω + 2 , as proved in Appendix C. Thus, the propagation condition is equivalent to 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 with whereQ r = r∂ r Ω. 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 which depends in general on the radius. In the case of general differential rotation but in regions where N 0 is much larger than f and Q (the norm of Q), the propagation condition (20) simplifies to The left inequality means that critical co-latitudes Θ c may exist for which ω 2 = f cos Θ c ( f cos Θ 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 (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 which is possible only when there is no critical surface in latitude.
In contrast, in regions where f is much larger than Q and N 0 , Eq. (20) becomes where A = f Q z sin Θ cos Θ + Q z 2 cos 4 Θ. By definition, Q vanishes at the centre. As N 0 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 subor 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.

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 0 , k 0 ), which is set by the initial condition (x 0 , k 0 ). 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 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 2 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 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 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 r e r + k θ e θ ) is different from the natural one associated with spherical coordinates (k = k r e r + k nat θ e θ /r). These equations are made explicit in Appendix D.
Article number, page 4 of 17 V. Prat et al.: Asymptotic theory of gravity modes in rotating stars. II. Impact of general differential rotation 3.2. Phase-space visualisation: the Poincaré surface of section Rays can be studied 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 3D space 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 e 0,max /(ωR e ), where N e 0,max 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).

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 dy/dt = f (y) with the initial condition y(0) = y 0 , 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 where h is the fixed step size. The ξ j are given implicitly by The computation of the coefficients a i j and b j follows from a 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 = i a 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 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 Eqs. (36) as For clarity in identifying the component index of X, we note that s 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 at y n and ξ q j . To reduce the cost of its computation, the derivatives of X used to form the 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 K = χ(r)(H − H 0 ), where H = ω(y) and H 0 = ω(y 0 ). This implies that the symplecticity-preserving equations of motion become 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.

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.
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 Ω 0 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: where 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 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 r sin 2 θΩ(Ω + 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 is the critical angular velocity and M is the stellar mass.

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 Sect. 4.2.1 and Sect. 4.2.2, respectively.

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 subinertial 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); ; 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.

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 the envelope. 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.

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 lowmass 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 Θ( f cos Θ + Q ⊥ ) with respect to the angular frequency ω. To simplify the discussion, we neglect here the centrifugal deformation. We thus consider the quantity κ defined by κ 2 = f cos θ( f cos θ + Q θ ), which is similar to a horizontal epicyclic frequency. Figure 10 shows how κ 2 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, κ 2 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, κ 2 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, κ 2 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 12c, 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  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. 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 12b, 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 12d, 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. 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 (subinertial, at the top) until it finally propagates past the rotation axis (super-inertial, at the bottom).
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.

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 subinertial 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 in the 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 semiclassical 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 gravitoinertial 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 . 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, 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). 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(Alvan et al. , 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 equation (14) we derived for nonaxisymmetric 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. 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.   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. 11. Blue trajectories correspond to invariant tori, green ones to island chains, and red ones to chaotic zones.  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. dk r dt = − r∂ r Ak r 2 + 2(r∂ r B − B)k r k θ + [r∂ r C − 2(C − ω 2 )]k θ 2 + r∂ r Dk c 2 + (D − ω 2 )r∂ r (k c 2 ) 2rω(k 2 + k c 2 ) , (D.8) dk θ dt = − ∂ θ Ak r 2 + 2(∂ θ B + A − ω 2 )k r k θ + (∂ θ C + 2B)k θ 2 + ∂ θ Dk c 2 + (D − ω 2 )∂ θ (k c 2 ) 2rω(k 2 + k c 2 ) . (D.9)