Asymptotic theory of gravity modes in rotating stars
II. Impact of general differential rotation
^{1}
IRFU, CEA, Université ParisSaclay,
91191
GifsurYvette,
France
email: vincent.prat@cea.fr
^{2}
Université Paris Diderot, AIM, Sorbonne Paris Cité, CEA, CNRS,
91191
GifsurYvette,
France
^{3}
Université de Toulouse, UPSOMP, IRAP,
Toulouse,
France
^{4}
CNRS, IRAP,
14 avenue Édouard Belin,
31400
Toulouse,
France
Received:
2
January
2018
Accepted:
11
March
2018
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 gravitoinertial waves.
Methods. We use a smallwavelength 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 gravitoinertial 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 gravitoinertial waves. We identify that for some regimes of frequency and differential rotation, the properties of gravitoinertial 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.
Key words: asteroseismology / waves / chaos / stars: oscillations / stars: rotation
© ESO 2018
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0;), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Differential rotation plays a key role in stellar dynamics, magnetism, and evolution. It triggers instabilities (such as RayleighTaylor, GoldreichSchubertFriecke, or shear instabilities) and interacts with largescale 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 solartype (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 nonrotating 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 gmode 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; 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 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 pmode 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 Btype (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 twodimensional (2D), in the sense that it cannot beseparated into onedimensional problems as in the nonrotating 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 socalled 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 lowfrequency 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 gravitoinertial 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 JeffreyWentzelKramersBrillouin (JWKB) asymptotic theory can take into account the compressibility effects that produce the backrefraction 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 semiclassical quantisation concepts to link acoustic rays to pressure modes, they predicted spectral patterns that have been successfully confronted to bidimensional 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 backrefraction of gravitoinertial 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 nonrotating case; (ii) island modes, where the energy is localised around a socalled 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 lowfrequency 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 gravitoinertial waves. First, we derive an eikonal equation (i.e. a local dispersion relation) for axisymmetric gravitoinertial 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. 

Open with DEXTER 
2 Eikonal equation for gravitoinertial 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
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.
Equation (1) can be written in the form (3)
where g_{0} is the effective gravity defined by (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.
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
ρ, 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 (10)
where is a rank2 tensor; the double gradient of a scalar function a is defined by in Cartesian coordinates; the symbol: denotes the double contraction of two tensors (); and is the complex amplitude of pressure fluctuations, such that (11)
where m is the azimuthal order of the waves and ω their angular frequency. Although the problem is initially threedimensional (3D), for a given m, it becomes 2D due to the quantisation in the azimuthal direction.
2.3 Shortwavelength 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 wavelike solutions of the form (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 secondorder derivatives, with . However, to ensure that waves are backrefracted near the stellar surface (which is needed to study modes) we also retain a zerothorderterm 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 (13)
where is the Dopplershifted 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 backrefraction of waves near the stellar surface. The only dependence on m is included in the Dopplershifted 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 (14)
where k_{φ} = m∕(rsinθ). In the general case, this equation canbe solved analytically for , but the solutions are nontrivial functions of structure, rotation, and k.
When there is no differential rotation, we recover the eikonal equation obtained in Paper I: (15)
When only axisymmetric waves are considered (k_{φ} = 0, and ), Eq. (14) reduces to (16)
When centrifugal deformation (which also includes effects of baroclinicity) is neglected, the eikonal equation in spherical coordinates reads (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: (18)
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, and , as proved inAppendix C. Thus, the propagation condition is equivalent to (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 (21)
where . Roots of Eq. (21) correspond to critical colatitudes, which delimit regions where waves can propagate. Thus, when a critical colatitude θ_{c} exists, it satisfies (25)
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 Eq. (20) simplifies to (26)
The left inequality means that critical colatitudes Θ_{c} may exist for which ω^{2} = fcosΘ_{c}(fcosΘ_{c} + Q_{⊥}). In solidbody 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 colatitudes exist, they verify (27)
which is possible only when ω^{2} ≤ f(f + Q_{⊥}) if Θ ≤ π∕2, or ω^{2} ≤f(f − Q_{⊥}) if Θ ≥ π∕2. Thus, if , there is no criticalsurface in latitude.
In contrast, in regions where f is much larger than Q and N_{0}, Eq. (20) becomes (28)
where . By definition, Q vanishes at the centre. As N_{0} also vanishes, the righthand 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 subinertial throughout (ω < f) avoid a region around the rotation axis while waves that are superinertial 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 superinertial 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 gravitoinertial 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_{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 v_{g} = ∇_{k}ω, where ∇_{k} is the gradient with respectto k. It reads (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^{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 (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
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 (). These equations are made explicit in Appendix D.
Fig. 2 Rotation profile with Ω_{C}∕Ω_{R} = 2, r_{C} = 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 superinertial regimes, respectively. 

Open with DEXTER 
3.2 Phasespace visualisation: the Poincaré surface of section
Rays can bestudied as dynamical systems with two degrees of freedom, so the phase space is fourdimensional (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 hypersurface, usually defined by fixing one phasespace 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 hypersurface 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 lowfrequency 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 , where is the value of the inner maximum of the radial profile of the BruntVäisälä frequency in the equatorial plane (see Paper I for more detail, in particular the simple analytical description of the PSS in the nonrotating case).
3.3 Numerical method
To properly characterise the dynamics of gravitoinertial rays of a system with differential rotation and a nontrivial 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 wellsampled 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 GaussLegendreRungeKutta (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_{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 sstage implicit GLRK method yields (35)
where h is the fixed step size. The ξ_{j} are given implicitly by (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} = ∑_{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 nontrivial way (e.g. SanzSerna 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 .
Since the ξ_{j} are implicit and f is a nonlinear function, we have implemented a NewtonRaphson solver to compute them. To do this, we construct a new vector function X(y_{n}, ξ_{j}) that concatenates Eq. (36) as
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 . 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 , where and . This implies that the symplecticitypreserving equations of motion become
Therefore, the only notable change in the method to accommodate these time step changes is in the righthandside vector function f of the GLRK method.
This new method has been compared to the fourthorder explicit RungeKutta 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 superinertial, 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. 

Open with DEXTER 
4 Gravitoinertial 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 selfconsistent 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 BruntVä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 threezones rotation profile: (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 solartype 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 antisolar 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 nonzero at the centre. To eliminate these pathological features, we regularise the rotation profile, adding the following corrections to the original profile: (43)
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 RayleighTaylor stable, that is, the quantity (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^{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 (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 subinertial 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 superinertial case. For that reason, we call it the purely superinertial 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 transinertial, 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 superinertial 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 superinertial, or reflection on the critical surface if the ray is locally subinertial. 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 transinertial trajectory showing the transition between the subinertial and superinertial behaviours. The ray bounces several times on the turning surface (subinertial, at the top) until it finally propagates past the rotation axis (superinertial, at the bottom). 

Open with DEXTER 
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 subinertial 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 superinertial 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 superinertial 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 transinertial 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 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Θ(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 κ^{2} = fcosθ(fcosθ + 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 (subinertial 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 (superinertial 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 superinertial 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 midinertial. When ω > κ_{max}, there is no critical surface.
In addition to the number of critical surfaces in the envelope, waves are either sub or superinertial 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 antisolar differential rotation, this would lead to the existence of regions that are unstable with respect to the RayleighTaylor 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 subinertial. 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 subinertial waves.
Regimes 3 and 4 are both subinertial in the envelope and superinertial in the core, thus transinertial. 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 subinertial in the core and midinertial 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 subinertial regimes. However, Fig. 12e shows that in regime 5, invariant tori that dominate the phase space (as in the purely subinertial case) coexist with island chains and relatively small chaotic zones. Moreover, the invariant tori seem to be less smooth than in the purely subinertial 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 superinertial in the core and midinertial 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 subinertial in the core and superinertial 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 superinertial. Figure 12h shows that in this regime, the phase space is dominated by invariant tori and island chains, similarly to the purely superinertial regime of Sect. 4.2.1
Fig. 5 PSS at Ω_{C}∕Ω_{R} = 1.1 and ω∕(2Ω_{R}) = 1.05 (left) and two transinertial rays: one belonging to an island chain (blue) and one belonging to a chaotic zone (red). 

Open with DEXTER 
Fig. 6 PSS at Ω_{C}∕Ω_{R} = 4 and ω∕(2Ω_{R}) = 2.4 (left) and one transinertial ray belonging to an island chain (blue). 

Open with DEXTER 
Fig. 7 Rotation profile with Ω_{C}∕Ω_{R} = 0.5, r_{C} = 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 superinertial regimes, respectively. 

Open with DEXTER 
Fig. 8 PSS (left) and examples of ray trajectories (right) for the three regimes shown in Fig. 7: sub, trans, and superinertial, 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. 

Open with DEXTER 
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 gravitoinertial 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 raytracing 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 inthe superinertial 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 BruntVäisälä frequency, see Paper I for more details). Between these two regimes, we find a new regime, called transinertial. 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 nonaxisymmetric modes may behave differently.
The large variety of domains of propagation generated by differential rotation means that gravitoinertial 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 superinertial regime is dominated by invariant tori. This suggests the existence of a nearby integrable system that could provide us with seismic diagnoses for lowfrequency gravitoinertial modes, as done in the case of uniform rotation in Prat et al. (2017).
In the present work we did not use selfconsistent 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 highfrequency 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 gravitoinertial modes based on ray dynamics should also be compared with numerically computed modes, using bidimensional 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 gravitoinertial 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 timedependent simulations of the excitation, propagation, and dissipation of nonlinear waves (Alvan et al. 2014, 2015).
In the near future, we plan to generalise the present work to nonaxisymmetric waves. This will require a significant amount of work, since the eikonal Eq. (14) we derived for nonaxisymmetric waves is more complex than the one for axisymmetric waves, and the Dopplershifting 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 nonaxisymmetric 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 gravitoinertial waves, which become magnetogravitoinertial 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, r_{E} = 0.7. 

Open with DEXTER 
Fig. 10 Value of κ^{2} = fcosθ(fcosθ + Q_{θ}) as a function of η = Ω_{D}∕Ω_{R} and cos^{2}θ 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. 

Open with DEXTER 
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. 

Open with DEXTER 
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. 

Open with DEXTER 
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. 

Open with DEXTER 
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^{1}. 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 timeharmonicity, 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): (A.1)
where and . It is possible to link the righthand side of the dot product to the BruntVäisälä frequency defined by (A.2)
where Γ_{1} is the first adiabatic exponent of the fluid. To do so, we use the relation (A.3)
which derives from Eq. (3). Using the definition of the effective gravity given in Eq. (4), Eq. (A.1) becomes (A.4)
where e_{∥} = −g_{0}∕g_{0} 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):
where the tensor is described in the basis (e_{∥}, e_{⊥}, e_{φ}) by the matrix (A.7)
and Θ is the angle between e_{z} and e_{∥}, which is equal to the colatitude θ when the centrifugal deformation is neglected. When is invertible, Eq. (A.6) can be used to express û as a function of and its gradient: (A.8)
and Γ is defined in Eq. (19). Finally, the combination of Eqs. (A.5) and (A.8) leads to a single equation for only: (A.10)
Equivalently, using the tensor identity , where is a tensor, yields the Poincaré equation (A.11)
Appendix B Determination of the surface term
Close to the surface, c_{s} becomes very small. Let us consider here the case where N_{0} remains finite near the surface. The dominant term of B is (B.1)
Using the relation , along with Eqs. (A.2) and (A.3), one can write (B.3)
where β = Γ_{1} − 1 − α and , 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 . Thus, (B.4)
Now, the wave equation Eq. (A.11) is put in the socalled normal form to get rid of first order terms. This is done by writing with a wellchosen 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 (B.5)
is the symmetric part of . To make all firstorder terms vanish, one would need (B.7)
Supposing that this is possible, this would lead to (B.8)
where ⊗ denotes the outer product of two vectors, which is a tensor defined by . The righthand 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 (B.9)
and Eq. (B.5) becomes (B.12)
We note that the term in B_{a} in Eq. (B.12) will be neglected when applying the JWKB approximation.
The expression for given in Eq. (B.14) is negative when α = 0 and Γ_{1} > 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_{0} becomes very large near the surface and scales as , 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 (B.15)
Appendix C Number of solutions of the equation Γ = 0
In the case of uniform rotation, Γ = 0 reduces to (C.1)
which always have two real solutions for ω^{2}.
In presence of differential rotation, some calculations are required. First, the discriminant of Γ = 0 seen as a quadratic equation for ω^{2} reads (C.2)
It can be rewritten in terms of Q_{∥} and Q_{⊥}: (C.3)
To go further, we need to consider the reduced discrimant of the equation Δ = 0 seen as a quadratic equation for : (C.4)
After simplification, one finally obtains (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: (D.1)
References
 Aerts, C., ChristensenDalsgaard, J., & Kurtz, D. W. 2010, Asteroseismology, (Dordrecht: Springer Science+Business Media B.V.) [CrossRef] [Google Scholar]
 Alvan, L., Brun, A. S., & Mathis, S. 2014, A&A, 565, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Alvan, L., Strugarek, A., Brun, A. S., Mathis, S., & Garcia, R. A. 2015, A&A, 581, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ando, H. 1985, PASJ, 37, 47 [NASA ADS] [Google Scholar]
 Ballot, J., Lignières, F., Reese, D. R., & Rieutord, M. 2010, A&A, 518, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ballot, J., Lignières, F., Prat, V., Reese, D. R., & Rieutord, M. 2012, in Progress in Solar/Stellar Physics with Helio and Asteroseismology, eds. H. Shibahashi, M. Takata, & A. E. LynasGray, ASP Conf. Ser., 462, 389 [NASA ADS] [Google Scholar]
 Ballot, J., Lignières, F., & Reese, D. R. 2013, in Numerical Exploration of Oscillation Modes in Rapidly Rotating Stars, Lect. Notes Phys., (Springer), 865, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Baruteau, C., & Rieutord, M. 2013, J. Fluid Mech., 719, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Beck, P. G., Montalban, J., Kallinger, T., et al. 2012, Nature, 481, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Benomar, O., Takata, M., Shibahashi, H., Ceillier, T., & García, R. A. 2015, MNRAS, 452, 2654 [NASA ADS] [CrossRef] [Google Scholar]
 Bouabid, M.P., Dupret, M.A., Salmon, S., et al. 2013, MNRAS, 429, 2500 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, B. P., Browning, M. K., Brun, A. S., Miesch, M. S., & Toomre, J. 2008, ApJ, 689, 1354 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., & Toomre, J. 2002, ApJ, 570, 865 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., Strugarek, A., Varela, J., et al. 2017, ApJ, 836, 192 [NASA ADS] [CrossRef] [Google Scholar]
 Butcher, J. C. 1963, J. Aust. Math. Soc., 3, 185 [CrossRef] [Google Scholar]
 Couvidat, S., García, R. A., TurckChièze, S., et al. 2003, ApJ, 597, L77 [NASA ADS] [CrossRef] [Google Scholar]
 Cowling, T. G. 1941, MNRAS, 101, 367 [NASA ADS] [Google Scholar]
 Deheuvels, S., García, R. A., Chaplin, W. J., et al. 2012, ApJ, 756, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Deheuvels, S., Doğan, G., Goupil, M. J., et al. 2014, A&A, 564, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Deheuvels, S., Ballot, J., Beck, P. G., et al. 2015, A&A, 580, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dintrans, B., & Rieutord, M. 2000, A&A, 354, 86 [NASA ADS] [Google Scholar]
 Dzhalilov, N. S., & Staude, J. 2004, A&A, 421, 305 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eckart, C. 1960, Hydrodynamics of oceans and atmospheres (New York: Pergamon Press) [Google Scholar]
 Fossat, E., Boumier, P., Corbard, T., et al. 2017, A&A, 604, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Friedlander, S., & Siegmann, W. L. 1982, Geophys. Astro. Fluid, 19, 267 [NASA ADS] [CrossRef] [Google Scholar]
 Fuller, J., Lecoanet, D., Cantiello, M., & Brown, B. 2014, ApJ, 796, 17 [NASA ADS] [CrossRef] [Google Scholar]
 García, R. A., TurckChièze, S., JiménezReyes, S. J., et al. 2007, Science, 316, 1591 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Gastine, T., Yadav, R. K., Morin, J., Reiners, A., & Wicht, J. 2014, MNRAS, 438, L76 [NASA ADS] [CrossRef] [Google Scholar]
 Guenel, M., Baruteau, C., Mathis, S., & Rieutord, M. 2016, A&A, 589, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hairer, E. 1997, Appl. Numer. Math., 25, 219 [CrossRef] [Google Scholar]
 Kurtz, D. W., Saio, H., Takata, M., et al. 2014, MNRAS, 444, 102 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, U., Neiner, C., & Mathis, S. 2014, MNRAS, 443, 1515 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, U., & Saio, H. 1993, MNRAS, 261, 415 [Google Scholar]
 Lee, U., & Saio, H. 1997, ApJ, 491, 839 [NASA ADS] [CrossRef] [Google Scholar]
 Lignières, F.,& Georgeot, B. 2009, A&A, 500, 1173 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maeder, A. 2009, Physics, Formation and Evolution of Rotating Stars, (Berlin Heidelberg: SpringerVerlag) [Google Scholar]
 Maeder, A., & Zahn, J.P. 1998, A&A, 334, 1000 [NASA ADS] [Google Scholar]
 Mathis, S. 2009, A&A, 506, 811 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., & de Brye N. 2011, A&A, 526, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., & de Brye N. 2012, A&A, 540, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., & Zahn, J.P. 2004, A&A, 425, 229 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., Decressin, T., Eggenberger, P., & Charbonnel, C. 2013, A&A, 558, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Matt, S. P., Do Cao, O., Brown, B. P., & Brun, A. S. 2011, Astron. Nachr., 332, 897 [NASA ADS] [CrossRef] [Google Scholar]
 Mirouh, G. M., Baruteau, C., Rieutord, M., & Ballot, J. 2016, J. Fluid Mech., 800, 213 [NASA ADS] [CrossRef] [Google Scholar]
 Mosser, B., Goupil, M. J., Belkacem, K., et al. 2012, A&A, 548, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Murphy, S. J., Fossati, L., Bedding, T. R., et al. 2016, MNRAS, 459, 1201 [NASA ADS] [CrossRef] [Google Scholar]
 Nielsen, M. B., Schunker, H., Gizon, L., Schou, J., & Ball, W. H. 2017, A&A, 603, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ouazzani, R.M., Dupret, M.A., & Reese, D. R. 2012, A&A, 547, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ouazzani, R.M., Salmon, S. J. A. J., Antoci, V., et al. 2017, MNRAS, 465, 2294 [NASA ADS] [CrossRef] [Google Scholar]
 Pápics, P. I., Tkachenko, A., Van Reeth, T., et al. 2017, A&A, 598, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pasek, M., Georgeot, B., Lignières, F., & Reese, D. R. 2011, Phys. Rev. Lett., 107, 11101P [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Pasek, M., Lignières, F., Georgeot, B., & Reese, D. R. 2012, A&A, 546, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pinçon, C., Belkacem, K., Goupil, M. J., & Marques, J. P. 2017, A&A, 605, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Prat, V., Lignières, F., & Ballot, J. 2016, A&A, 587, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Prat, V., Mathis, S., Lignières, F., Ballot, J., & Culpin, P.M. 2017, A&A, 598, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reese, D., Lignières, F., & Rieutord, M. 2006, A&A, 455, 621 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rieutord, M. 2006, in EAS Pub. Ser., eds. M. Rieutord & B. Dubrulle, 21, 275 [CrossRef] [EDP Sciences] [Google Scholar]
 Rieutord, M. 2009, in The Rotation of Sun and Stars (Berlin: Springer Verlag), Lect. Notes Phys., 765, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Rieutord, M., & Espinosa Lara F. 2013, in Lect. Notes Phys., eds. M. Goupil, K. Belkacem, C. Neiner, F. Lignières, & J. J. Green, (Berlin: Springer Verlag), 865 49 [NASA ADS] [CrossRef] [Google Scholar]
 Rieutord, M., Espinosa Lara, F., & Putigny, B. 2016, J. Comput. Phys., 318, 277 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, T. M. 2015, ApJ, 815, L30 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, T. M., Lin, D. N. C., McElwaine, J. N., & Lau, H. H. B. 2013, ApJ, 772, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Saio, H. 1981, ApJ, 244, 299 [NASA ADS] [CrossRef] [Google Scholar]
 Saio, H., & Takata, M. 2014, PASJ, 66, 58 [NASA ADS] [Google Scholar]
 Saio, H., Kurtz, D. W., Takata, M., et al. 2015, MNRAS, 447, 3264 [NASA ADS] [CrossRef] [Google Scholar]
 SanzSerna, J. M. 1988, BIT Numer. Math., 28, 877 [CrossRef] [Google Scholar]
 Schatzman, E. 1993, A&A, 279, 431 [NASA ADS] [Google Scholar]
 Takata, M. 2014, PASJ, 66, 80 [NASA ADS] [Google Scholar]
 Takata, M., & Saio, H. 2013, PASJ, 65, 68 [NASA ADS] [Google Scholar]
 Talon, S., & Charbonnel, C. 2005, A&A, 440, 981 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Talon, S., & Charbonnel, C. 2008, A&A, 482, 597 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Thompson, M. J., Toomre, J., Anderson, E. R., et al. 1996, Science, 272, 1300 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Townsend, R. H. D. 2003, MNRAS, 340, 1020 [NASA ADS] [CrossRef] [Google Scholar]
 Triana, S. A., Moravveji, E., Pápics, P. I., et al. 2015, ApJ, 810, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Triana, S. A., Corsaro, E., De Ridder, J., et al. 2017, A&A, 602, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Van Reeth, T., Tkachenko, A., & Aerts, C. 2016, A&A, 593, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zahn, J.P. 1983, in SaasFee Advanced Course 13: Astrophysical Processes in Upper Main Sequence Stars, eds. A. N. Cox, S. Vauclair, & J. P. Zahn, 253 [Google Scholar]
 Zahn, J.P. 1992, A&A, 265, 115 [NASA ADS] [Google Scholar]
 Zahn, J.P., Talon, S., & Matias, J. 1997, A&A, 322, 320 [Google Scholar]
All Figures
Fig. 1 Coordinate systems. 

Open with DEXTER  
In the text 
Fig. 2 Rotation profile with Ω_{C}∕Ω_{R} = 2, r_{C} = 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 superinertial regimes, respectively. 

Open with DEXTER  
In the text 
Fig. 3 PSS (left) and examples of ray trajectories (right) for the three regimes identified in Fig. 2: sub, trans, and superinertial, 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. 

Open with DEXTER  
In the text 
Fig. 4 Detail of a chaotic transinertial trajectory showing the transition between the subinertial and superinertial behaviours. The ray bounces several times on the turning surface (subinertial, at the top) until it finally propagates past the rotation axis (superinertial, at the bottom). 

Open with DEXTER  
In the text 
Fig. 5 PSS at Ω_{C}∕Ω_{R} = 1.1 and ω∕(2Ω_{R}) = 1.05 (left) and two transinertial rays: one belonging to an island chain (blue) and one belonging to a chaotic zone (red). 

Open with DEXTER  
In the text 
Fig. 6 PSS at Ω_{C}∕Ω_{R} = 4 and ω∕(2Ω_{R}) = 2.4 (left) and one transinertial ray belonging to an island chain (blue). 

Open with DEXTER  
In the text 
Fig. 7 Rotation profile with Ω_{C}∕Ω_{R} = 0.5, r_{C} = 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 superinertial regimes, respectively. 

Open with DEXTER  
In the text 
Fig. 8 PSS (left) and examples of ray trajectories (right) for the three regimes shown in Fig. 7: sub, trans, and superinertial, 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. 

Open with DEXTER  
In the text 
Fig. 9 Rotation profile with Ω_{D}∕Ω_{R} = 0.2, α_{E} = 40, r_{E} = 0.7. 

Open with DEXTER  
In the text 
Fig. 10 Value of κ^{2} = fcosθ(fcosθ + Q_{θ}) as a function of η = Ω_{D}∕Ω_{R} and cos^{2}θ 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. 

Open with DEXTER  
In the text 
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. 

Open with DEXTER  
In the text 
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. 

Open with DEXTER  
In the text 
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. 

Open with DEXTER  
In the text 