Asymptotic theory of gravity modes in rotating stars
I. Ray dynamics
^{1} MaxPlanck Institut für Astrophysik, KarlSchwarzschildStr. 1, 85748 Garching bei München, Germany
email: vprat@mpagarching.mpg.de
^{2} Université de Toulouse, UPSOMP, IRAP, 31028 Toulouse, France
^{3} CNRS, IRAP, 14 avenue Édouard Belin, 31400 Toulouse, France
Received: 13 November 2015
Accepted: 22 December 2015
Context. The seismology of earlytype stars is limited by our incomplete understanding of gravitoinertial modes.
Aims. We develop a shortwavelength asymptotic analysis for gravitoinertial modes in rotating stars.
Methods. The WentzelKramersBrillouin approximation was applied to the equations governing adiabatic small perturbations about a model of a uniformly rotating barotropic star.
Results. A general eikonal equation, including the effect of the centrifugal deformation, is derived. The dynamics of axisymmetric gravitoinertial rays is solved numerically for polytropic stellar models of increasing rotation and analysed by describing the structure of the phase space. Three different types of phasespace structures are distinguished. The first type results from the continuous evolution of structures of the nonrotating integrable phase space. It is predominant in the lowfrequency region of the phase space. The second type of structures are island chains associated with stable periodic rays. The third type of structures are large chaotic regions that can be related to the envelope minimum of the BruntVäisälä frequency.
Conclusions. Gravitoinertial modes are expected to follow this classification, in which the frequency spectrum is a superposition of subspectra associated with these different types of phasespace structures. The detailed confrontation between the predictions of this raybased asymptotic theory and numerically computed modes will be presented in a companion paper.
Key words: asteroseismology / chaos / waves / stars: oscillations / stars: rotation
© ESO, 2016
1. Introduction
As compared to the wealth of new information obtained from the seismology of solartype and redgiant stars (Miglio 2015), reliable seismic diagnostics on nonevolved, intermediatemass and massive pulsators (γ Doradus, δ Scuti, β Cephei, SPB or Be stars) are scarce and, in any case, restricted to atypical slowly rotating stars (Kurtz et al. 2014; Saio et al. 2015; Triana et al. 2015). This stems from our limited understanding of the effects of rotation on oscillation modes (Lignières 2013; Townsend 2014; Reese 2015), at least for the high rotation rates of typical nonevolved, intermediatemass and massive stars (Royer et al. 2007; Lignières 2013). The two classical approximations to treat rotational effects have been the perturbative expansion in Ω /ω, where Ω is the rotation rate and ω is the pulsation frequency (Saio 1981), and the socalled traditional approximation of the Coriolis acceleration (Eckart 1960). The perturbative expansion happens to be valid for small rotation rates but not for those of most earlytype stars (Reese et al. 2006; Ballot et al. 2010, 2013). The traditional approximation greatly simplifies the understanding of Coriolis effects on lowfrequency gmodes (e.g. Townsend 2003; Bouabid et al. 2013) but the range of validity and accuracy of this approximation are not known. The centrifugal deformation in particular is not taken into account in the traditional approximation.
More recently, the development of dedicated numerical codes has made it possible to compute modes and frequencies accurately without approximating the Coriolis and centrifugal effects on the modes (Reese et al. 2006, 2009; Ouazzani et al. 2012). This is a crucial step towards a confrontation with observational data. However, in addition to accurate theoretical frequencies, mode identification of observed frequencies requires a priori knowledge of the mode properties. In the case of slowly rotating stars, our knowledge of the expected frequency patterns is important for the interpretation of observed frequency spectra. It comes from the shortwavelength asymptotic theory of oscillation modes, which in the case of spherically symmetric nonrotating stars provides analytical expressions of the asymptotic frequency patterns. When applied to rapidly rotating stars, the same shortwavelength asymptotic analysis is not as straightforward. As in geometrical optics, it first leads to a ray model that describes wave propagation. Then, the conditions to produce a mode from the positive interference between rays that need to be established. Methods and concepts to infer the properties of the modes from those of the rays have been mostly developed in the field of quantum physics. The basic idea is to study rays as trajectories of a Hamiltonian system and then to relate the modes to the phasespace structures of the ray dynamics. For acoustic modes in rotating stars, such a raybased asymptotic theory provided a physical classification of the modes as well as quantative predictions that have been successfully confronted with numerically computed highfrequency pmodes (Lignières & Georgeot 2008, 2009; Pasek et al. 2011, 2012).
In this paper, we intend to construct a similar raybased asymptotic theory for gravitoinertial modes. Equations governing gravitoinertial rays in rotating stars are derived and then solved to obtain a global view of the properties of axisymmetric gravitoinertial rays and their evolution with stellar rotation. The detailed confrontation with numerically computed modes will be presented in a companion paper, although a first successful comparison in the particular case of socalled rosette modes has been already presented in a proceedings paper (Ballot et al. 2012).
Ray models are usually obtained from the WentzelKramersBrillouin (WKB) approximation of the equations governing smallamplitude perturbations, and this approximation is valid when the wavelength of the perturbation is much shorter than the length scale of the background spatial variations. Gravitoinertial ray models have been first developed for the ocean and the Earth’s atmosphere. One of the main motivation in this context is to model angular momentum transport and deposition by waves (Fritts & Alexander 2003; Marks & Eckermann 1995; Broutman et al. 2004). In stellar physics, such WKB ray models have been limited to pure gravity waves (Gough 1993; Alvan et al. 2015) and have not been used to study angular momentum transport or gravitoinertial modes. Instead, this is another ray theory, the method of characteristics, which has been applied by Dintrans & Rieutord (2000) to the problem of stellar gravitoinertial waves and their relation with gravitoinertial modes. As detailed in Harlander & Maas (2006), the WKB approximation and the method of characteristics differ in general, although they can be identical in some specific cases. For what concerns stellar gravitoinertial waves, they differ in that the method of characteristics cannot take into account the compressibility effects that produce the backrefraction of the waves approaching the stellar surface. We thus expect WKB ray models to be more realistic than the method of characteristics.
The paper is organised as follows. In Sect. 2, the WKB approximation is applied to the wave equation governing small perturbations about a model of a uniformly rotating star. A general eikonal equation for gravitoinertial waves is derived and the domains of propagation of these waves inside the star are discussed. In Sect. 3, the ray model for axisymmetric rays is derived from the eikonal equation, its Hamiltonian nature is emphasised and the standard tool to investigate the phasespace structure of the Hamiltonian ray dynamics is presented. A detailed numerical study of the ray dynamics for uniformly rotating polytropic models of stars is then presented in Sect. 4. The interpretation of the dynamics in terms of gravitoinertial modes is discussed in Sect. 5 and conclusions are provided in Sect. 6.
2. Eikonal equation for gravitoinertial waves
In this section, we derive an eikonal equation for shortwavelength, adiabatic gravitoinertial waves in a uniformly rotating polytropic model of star. The stellar model is specified in Sect. 2.1, the equations for small perturbations about this model are put in a convenient form in Sect. 2.2, the eikonal equation obtained by applying the WKB approximation to the perturbation equations is presented in Sect. 2.3, and the domains of propagation of gravitoinertial waves are discussed in Sect. 2.4. Only the main equations are presented in this section because their derivation are detailed in Appendix A.
2.1. Polytropic model of rotating star
The stellar model is a selfgravitating uniformly rotating monatomic gas that verifies a polytropic relation, which is assumed to give a reasonably good approximation of the relation between pressure and density in real stars (Hansen & Kawaler 1994). The governing equations can be written as
where P_{0} is the pressure, ρ_{0} the density, K the polytropic constant, μ the polytropic index, ψ_{0} the gravitational potential, Ω the rotation rate, w the distance to the rotation axis, and G the gravitational constant.
The uniform rotation ensures that the flow is barotropic. A pseudoenthalpy h_{0} = ^{∫}dP_{0}/ρ_{0} = (1 + μ)P_{0}/ρ_{0} can then be introduced and quantities describing the equilibrium model, such as the sound speed , where Γ_{1} = 5 / 3 is the first adiabatic exponent of the gas; the effective gravity ; and the BruntVäisälä frequency, given by the relation are expressed in terms of h_{0}: (4)In polytropic models of stars, the density ρ_{0} vanishes at the surface. Moreover, the polytropic relation Eq. (1)implies that the sound speed c_{s} also vanishes, and the last of Eqs. (4)shows that the BruntVäisälä frequency N_{0} goes to infinity. These properties are not necessarily valid in more realistic stellar models.
Specifying the mass and rotation rate of the star is not sufficient to determine the polytropic model in physical units. This requires fixing an additional parameter, for example the stellar radius (suitable parameter choices are discussed in Hansen & Kawaler 1994; for the nonrotating case, and in ChristensenDalsgaard & Thompson 1999, for the rotating case). In the following, however, we only present dimensionless quantities that do not depend on the choice of this additional parameter. The rotation rate Ω is compared to , the limiting rotation rate for which the centrifugal acceleration equals the gravity at the equator, where M is the stellar mass and R_{e} the equatorial radius. Other frequency units important for our problem are (i) the inverse of the dynamical timescale associated with hydrostatic equilibrium , which provides a lower bound for acoustic wave frequencies, where R_{p} is the polar radius; (ii) the BruntVäisälä frequency N_{0}, an upper bound for gravity wave frequencies; and (iii) twice the rotation rate f = 2Ω, the upper bound for pure inertial wave frequencies.
Fig. 1 Top: solid black line shows the profile of the BruntVäisälä frequency N_{0} normalised by ω_{0} for the nonrotating stellar model. The red lines show N_{0} for the Ω = 0.84Ω_{K} model along the polar (solid line) and equatorial (dashed line) radii. A thick red tick on the xaxis indicates the polar radius R_{p}. Top: map in a meridional plane of N_{0} normalised by ω_{0} for the Ω = 0.84Ω_{K} model. The black line is the star surface. 

Open with DEXTER 
The spatial distribution of the BruntVäisälä frequency is shown in Fig. 1. In the nonrotating spherically symmetric case (see the black curve on the upper figure), N_{0} vanishes at the centre, has a local maximum, N_{0,max}, around r = 0.28R and then shows a local minimum, N_{0,min}, in the stellar envelope near r = 0.74R. A simple analytical expression of N_{0}, which is valid for envelope models, (see Appendix E) shows that a local minimum is present at r = 3 / 4R for all polytropic indices, suggesting that such a minimum of N_{0} is a generic feature of stars with radiative envelopes. The BruntVäisälä frequency distribution becomes anisotropic in centrifugally deformed stars. In particular, the amplitudes and positions of the local minima are clearly different along the polar and the equatorial axes. This is apparent in Fig. 1, which shows isocontours of N_{0} (bottom) as well as the polar and equatorial radial profiles of N_{0} (top) for a Ω / Ω_{K} = 0.84 model. In the following, the local minima along the polar and equatorial axes are denoted and , respectively. By comparison, the variation of the local maximum of the BruntVäisälä frequency between its values along the pole () and along the equatoral radius () remains small. This is because the deviation from sphericity is weaker in the central layers.
2.2. Perturbation equations and boundary conditions
Timeharmonic, smallamplitude perturbations of the stellar model are studied under two approximations. The first approximation is to neglect viscous and thermal dissipations. Both processes have a very small effect on the value of oscillation frequencies, at least in the frequency range of observed modes. Nonadiabatic processes play an essential role in the oscillation amplitude, in particular, through excitation mechanisms such as the κ mechanism, but determining mode amplitudes is outside the scope of the present study. The second approximation, known as the Cowling approximation, is to neglect perturbations of the gravitational potential induced by density perturbations. It is justified for oscillations of small wavelengths and, therefore, fully compatible with the smallwavelength regime considered here. Under these two approximations, the linear equations governing the evolution of smallamplitude perturbations can be written as
where v, ρ, and P, are respectively the Eulerian perturbations of velocity, density, and pressure, and f = 2Ω is the rotation vector.
Before applying the shortwavelength approximation, the perturbation equations for v, ρ, and P are reduced to a single equation for , the complex amplitude of the timeharmonic pressure perturbation . According to the calculations detailed in Appendix A.1, we obtain the equation (8)where the operators and Δ_{⊥} are related to the unit vectors parallel to the rotation axis, e_{z} = f/f, and the effective gravity, e_{∥} = −g_{0}/g_{0}, by The expressions of , and M′ are given by Eqs. ((A.39)–(A.41)), respectively. As a linear combination of f and g_{0}, the vector is contained in a meridional plane, whereas is colinear to f ∧ g_{0} and thus perpendicular to this plane. First order derivatives in the meridional plane can be suppressed from this equation by introducing a variable and by choosing the function a adequately (see Appendix A.2). We then obtain a wave equation for acoustic and gravitoinertial waves in a uniformly rotating, centrifugally deformed, barotropic star as follows: (11)where e_{φ} is the unit vector in the azimuthal direction and is given by Eq. (A.53). The first term on the righthand side (RHS) of Eq. (11) is associated with sound waves whereas the next two terms on the RHS produce gravitoinertial waves. The fourth term on the RHS is nonzero only for nonaxisymmetric perturbations in rotating stars. It contains the latitudinal derivative of the Coriolis parameter fcosθ, which gives rise to Rossby and Kelvin waves. Such a simple form for the wave equation comes at the cost of very complicated expressions for a and as functions of c_{s},N_{0},g_{0},f, and ω (see Eqs. (A.48) and (A.54)).
2.3. The shortwavelength approximation of the perturbation equations
The ray model results from the WKB approximation of the wave Eq. (11). It consists in looking for wavelike solutions of the form (12)assuming that the associated wavelength λ_{w} ~ ∥ ∇Φ ∥ ^{1} is much shorter than the typical length scale on which the background medium varies, hereafter denoted L_{b}.
The phase Φ and the amplitude A of are expanded into power series of 1 / Λ, where Λ = L_{b}/λ_{w} is assumed to be large, yielding(13)Introducing these expansions into the wave Eq. (11) and considering only the dominant terms, we obtain an equation for Φ_{0}, the socalled eikonal equation. The amplitude A_{0} is determined at the next order. The eikonal equation obtained with this procedure depends on the frequency range considered. If ω is of the order of Λ, the dominant terms of the wave Eq. (11) are the lefthand side (LHS) and the first term on the RHS, leading to an eikonal equation for sound waves ω = c_{s}k, where k = ∇Φ_{0} and k = ∥ k ∥. If , the dominant terms are those that contain second derivatives of , leading to an eikonal equation for gravitoinertial waves. The fourth term on the RHS of the wave Eq. (11)is a firstorder azimuthal derivative of and is therefore negligible when . Thus, Rossby waves are not included in the eikonal equation if one considers the regime or if one restricts to axisymmetric perturbations. Before writing down the eikonal equation for gravitoinertial waves, the role of the socalled constant term must be clarified. Indeed, if a wave, acoustic or gravitoinertial, of vanishingly small wavelength approaches the stellar surface, it will cross the surface without being refracted back into the star. This regime corresponds to the limit λ_{w} → 0, where the constant term is negligible and the resulting eikonal equation contains no surface refraction effect.
We are interested, however, in waves that are refracted back because oscillations modes are the result of wave interferences within the star. As known from previous studies on acoustic and gravity waves (Aerts et al. 2010), this backrefraction occurs for outwardstravelling waves that encounter atmospheric layers whose pressure scale heights are smaller that the wavelength of the wave. We also know that the pressure scale height strongly decreases as one approaches the surface in such a way that if L_{b} is of the order of the pressure scale height in the interior and H_{s} is the surface pressure scale height, we obtain H_{s}/L_{b} ≪ 1. Thus, there exists a wavelength range for which we have at the same time λ_{w}/L_{b} = 1 / Λ ≪ 1 inside the star and λ_{w} of the order of or larger than H_{s}. In the wave equation, this surface refraction effect is accounted for by the constant term because is inversely proportional to the square of the pressure scale height (see Eq. (16) below). Thus, an eikonal equation containing this term is simply obtained by assuming in the expansion given by Eq. (13). This yields (14)where k_{z} = e_{z}·k, k_{φ} = e_{φ}·k and k_{⊥} = e_{⊥}·k, where the unit vector e_{⊥} is defined as e_{⊥} = e_{φ} ∧ e_{∥}. In Eq. (14), the term in results from the fact that ∇_{⊥} is not just the derivative along e_{⊥}. It also contains derivatives in φ, and so does Δ_{⊥}.
Although the full expression of is a very complex function of the equilibrium quantities c_{s},N_{0},g_{0},f, and of ω (see Eq. (A.54)), we only need to estimate where it becomes of the same order of magnitude as the other terms in the dispersion relation, that is near the surface. The increase of towards the surface is due to the sharp increase of N_{0} and 1 /c_{s}, while g_{0} does not vary much in the stellar envelope. In a polytropic model of star, we also have N_{0} ∝ 1 /c_{s} near the surface. In order to only retain the dominant term, we look for an expansion of in powers of c_{s} and find that can be written as the ratio between two polynomials of , namely (15)where the terms , K = (f·g_{0})^{2}−ω^{2}g_{0}^{2}, G = ω^{2}(ω^{2}−f^{2}) and the terms C_{0} and C_{1} are at the surface. In the superinertial^{1} regime ω>f, K does not vanish. The dominant term of is thus and we obtain (16)where cosΘ = −f.g_{0}/fg_{0}. In the subinertial regime ω<f, both K and C_{0} vanish at the critical angles Θ_{c} and π−Θ_{c} such that cos^{2}Θ_{c} = ω^{2}/f^{2}. In this regime and in the regions close to the stellar surface and the critical angles, is given by (17)The detail of the asymptotic expression of is presented in Appendices A.3 and A.4. In the following, we always use the expression (16) for . As we see in the next sections, inspection of the eikonal and ray equations near the critical angle indicates that rays tend to avoid the surface layers close to the critical angle. This is further confirmed by our computations that do not show rays reaching the surface close to the critical angle.
Using Eq. (16), the eikonal equation can be written as (18)with (19)This equation is similar to the local dispersion relations generally used to study gravitoinertial waves in the Earth’s atmosphere (Fritts & Alexander 2003) and to construct associated ray models (Marks & Eckermann 1995). The present relation is nevertheless more general in that it takes the centrifugal deformation of the star into account and does not make the traditional approximation. In the following we shall restrict our study to axisymmetric perturbations, that is to the case k_{φ} = 0.
2.4. Domains of propagation
The eikonal Eq. (18) allows us to discuss the domain of propagation of gravitoinertial waves in a more general context than previous works. In particular, when compared to the work of Dintrans & Rieutord (2000), the following discussion includes the effects of the centrifugal deformation and, through the k_{c} term, those of the background compressibility.
According to Eq. (12), nonevanescent wave solutions are possible only when both components of the wave vector k are real numbers. The region of space where the eikonal equation admits real solutions thus determines the domain where rays can propagate. After projecting k onto an orthogonal basis, the eikonal equation takes the general form of a conic. Using (e_{∥},e_{⊥}) as a basis and the relations and k_{z} = k_{∥}cosΘ−k_{⊥}sinΘ, the eikonal Eq. (18) is rewritten as (20)To discuss the constraints on the domain of propagation that can be derived from this relation, it is useful to first consider the nonrotating case. The eikonal equation then becomes (21)with k_{⊥} = k_{θ} = ± L/r and L is the norm of the vector x ∧ k, which is equivalent to the angular momentum. This equation corresponds to the highsoundspeed limit of eikonal equations already derived in the nonrotating case, for example by Gough (1993). The condition that is positive translates into and this inequality fully specifies the spatial limit of the resonant cavity for a wave of given ω and L. In internal regions where , it further simplifies into ω<N_{0}, a condition that no longer depends on L or k_{θ}. Closer to the surface, ω ≪ N_{0} and the condition for propagation is ω<c_{1}S_{L}, where S_{L} = Lc_{s}/r is the Lamb frequency and c_{1} a constant that can be expressed in terms of Γ_{1} and μ. If this last condition is fullfilled at the surface, the gravity ray escapes the star. Otherwise, the ray is refracted back into the star and the relation ω = c_{1}S_{L} determines the outer radius of the domain of propagation. An important point is that the condition for the backrefraction in the outer layers of the star depends on L or equivalently on k_{θ}.
Coming back to the rotating case, the condition that k_{∥} is real can be expressed as (22)where δ is the reduced discriminant of the eikonal equation viewed as a quadratic equation for the variable k_{∥}, and (23)The sign of Γ determines whether gravitoinertial waves can propagate in the inner region where k_{c} can be neglected. In contrast, the outer limit of the domain of propagation is determined by the condition δ ≥ 0 and, as expected from the condition of propagation in the nonrotating case, it depends on the value of . However, contrary to the nonrotating case, is not known a priori because rotation breaks the spherical symmetry and the associated conservation of the norm of the angular momentum. The consequence is that in the general case of gravitoinertial waves in a rotating star, we are not able to predict the outer limit of the resonant cavity from the eikonal equation alone. It is thus necessary to solve the ray dynamics to find out whether a ray of a given frequency remains confined within the star and, if this is the case, the shape of the outer limit of its resonant cavity. In the ray dynamics calculation presented in this paper, all rays stay inside the stars because, at the surface of polytropic models, k_{c} is infinite and δ is thus necessarily negative. In real stars, k_{c} remains finite at the surface, and this implies that rays with sufficently high values of k_{⊥} near the surface can propagate outside the stars. Waves associated with such rays are expected to have very high spatial frequencies at the surface, and are thus unlikely to be visible. However, they can play a role in the transport of angular momentum out of the stars.
In the following of the section, we thus only discuss the domain of propagation in the internal region where the term can be neglected, since this approximation is valid from the centre up to a certain radius that depends on the given ray. The condition of propagation Γ ≥ 0 is equivalent to (24)with (25)If the centrifugal deformation is neglected, that is if we replace Θ by the colatitude θ, this condition of propagation is equivalent to that found by Dintrans & Rieutord (2000). Two limit cases, f ≪ N_{0} and N_{0} ≪ f, are relevant since the former occurs in the bulk of typical moderately rotating stellar radiative zones, whereas the latter is verified in convective zones or close to the centre of the star where N_{0} vanishes. If f ≪ N_{0}, the conditions (24) simplify to (26)On the other hand, if N_{0} ≪ f, the conditions (24) become (27)where in the limit of vanishing N_{0}, the ω<f rule for pure inertial waves is recovered.
We now consider the particular case of uniformly rotating μ = 3 polytropic models. Close to the centre of the star, the approximate conditions (27) hold because N_{0} vanishes there. At some distance from the centre, the BruntVäisälä frequency always becomes larger than f. Moreover, for moderate rotators, the approximate condition (26) is valid from a small radius up to the surface. Indeed, as long as Ω ≤ 0.38Ω_{K} and r> 0.1R_{e}, the ratio is always smaller than 0.10.
Fig. 2 Typical domains of propagation (grey areas) given by the Γ(ω) ≥ 0 condition as a function of the gravitoinertial wave frequency ω in a Ω = 0.38Ω_{K} polytropic model of star. In general, the outer limit of the domain, here loosely represented by the stellar surface, cannot be predicted from ω alone and is obtained by solving the ray dynamics. Panel a) corresponds to a superinertial frequency (ω>f) and panel b) to a subinertial frequency (ω<f). The two bottom figures show domains with an intermediate forbidden zone when frequencies are between the local minimum and maximum of the BruntVäisälä frequency. Panel c) corresponds to frequencies in the interval and also exists in nonrotating stars, whereas panel d) corresponds to frequencies in the interval and is due to the centrifugal deformation of rotating stars. 

Open with DEXTER 
Two typical domains of propagation are shown in panels (a) and (b) of Fig. 2 in the superinertial and subinertial frequency ranges, respectively. The rotation of the stellar model is Ω = 0.38Ω_{K}. For the superinertial case, the condition of propagation reduces to , which leads to the domain shown in panel (a), including the weakly prolate avoided central region. By contrast, in the subinertial regime the centre is never forbidden because conditions (27) are automatically fulfilled for ω<f. Away from this central region, the conditions of propagation (26) reduce to f^{2}cos^{2}Θ <ω^{2} for these waves. The resulting domain of propagation is shown in panel (b). A noticeable effect of the centrifugal deformation is that the limit of the domain is no longer strictly conical away from the central region. As the centrifugal deformation of the equipotentials becomes significant towards the surface, the boundary of the domain of propagation is bent towards lower latitudes to keep Θ_{c} = arccos(ω/f), the angle between the direction of the effective gravity and the rotation axis, constant.
The local minimum of the BruntVäisälä frequency present in the stellar envelope (see Fig. 1) induces an additional forbidden region for frequencies between and . As illustrated on panel (c) for the Ω = 0.38Ω_{K} model, the domains of progagation below and above the position of are disjoint when N_{0,min}<ω<N_{0,max} in nonrotating stars or when in rotating stars. In centrifugally deformed stars, waves with frequencies in the range also show a forbidden region, but this time it is limited to equatorial regions so that the inner and outer domains communicate through polar regions. An example is shown in panel (d) for the Ω = 0.38Ω_{K} model also. We observe that the outer limit of the forbidden region is oblate whereas the inner limit is prolate. As can be inferred from the distribution of N_{0} shown in Fig. 1, this property is related to the oblate and prolate shapes of the isocontours of the BruntVäisälä frequency taken in the interval. These peculiar resonant cavities appear in a frequency range that increases with rotation as the ratio goes from 1 in the nonrotating case to 1.70 at Ω = 0.84Ω_{K}.
3. The ray model
In this section, the ray model for progressive gravitoinertial waves is derived from the eikonal equation. Its Hamiltonian structure is emphasised in Sect. 3.1 as it enables us to investigate ray paths using tools and concepts developed for Hamiltonian dynamical systems. The geometrical structure of the phase space in particular is the key to understand the nature of the dynamics. Despite its high dimensionality, the phase space can be explored numerically through cuts known as Poincaré surfaces of section (hereafter PSS). The numerical method and PSS are presented in Sects. 3.2 and 3.3, respectively.
3.1. The Hamiltonian equations describing the ray path and the wave vector evolution along the path
The eikonal Eq. (18) is a partial differential equation (PDE) for the phase function Φ_{0}(x). Instead of solving this equation as a PDE, the gravitoinertial ray model consists in searching for solutions along some path x(t) called the ray path. One then has to solve the coupled differential equations that govern the ray path and the evolution of k along the path.
It has been shown that this problem can be written under very general circumstances in a Hamiltonian form (Lighthill 1978; Ott 1993). In the present case, the Hamiltonian is obtained by writing the eikonal equation in the form ω = H(x,k) and the ray path is defined by the group velocity. For a given coordinate system [ x_{i} ], the Hamiltonian equations of the ray dynamics are (Lignières 2011) where the are the covariant components of k on the natural basis associated with the coordinate system [ x_{i} ]. In Appendix B, these equations are written with spherical coordinates and the k components on the usual unit vectors e_{r} and e_{θ}.
A classical property of gravitoinertial waves is the orthogonality of the group velocity v_{g} = dx/ dt and the phase velocity v_{p} = ωk/k^{2}. Here, this property is valid in the interior but broken when k_{c} is not negligible. This is apparent in the following formula: (30)whose derivation is detailed in Appendix C.
Another important characteristic is the behaviour of the rays whenever they reach the limits of the domain of propagation. This can be discussed from the relations also derived in Appendix C. Equation (31) shows that the component of the group velocity parallel to e_{∥} vanishes at the limits of the domain of propagation, that is when δ = 0. If this occurs when k_{c} is negligible, Eq. (22) finds that Γ also vanishes there. Thus, according to Eq. (32), both components of the group velocity vanish at the limit of the domain of propagation and the turning point is an ordinary cusp point. However, if δ = 0 is reached closer to the stellar surface, where k_{c} cannot be neglected, Eqs. (22) and (32) imply that Γ and dx/ dt·e_{⊥} do not vanish. The trajectory is then locally parallel to e_{⊥}, that is tangent to an equipotential, before it is redirected towards the stellar interior. Figure 3 illustrates the two types of turning points for a given trajectory in the superinertial regime. The capacity of the numerical method to accurately compute the trajectories near cusp points has been successfully tested against an analytical solution (see Appendix D).
Fig. 3 Gravitoinertial ray showing the two types of turning points, a cusp point on the Γ = 0 boundary (in magenta) and a smooth backrefraction in the surface layers where the decreasing density scale height becomes of the order of the wavelength in that direction. 

Open with DEXTER 
3.2. Numerical method for the ray dynamics
The gravitoinertial ray dynamics has been investigated by integrating numerically Eqs. (B.1)–(B.4) using a 5thorder RungeKutta method. The step size of the integration is chosen automatically to keep the local error estimate smaller than a given tolerance. The background polytropic model of star is solved numerically with an iterative scheme, which is described in Rieutord et al. (2005).
We checked that the PSS shown in Sect. 4 are not significantly modified by decreasing the maximum allowed local error of the RungeKutta scheme. We also checked the influence of the resolution of the background polytropic stellar model. Finally, the conservation of the Hamiltonian is used as an independent accuracy test of the computations.
3.3. Phasespace visualisation: the Poincaré surface of section
The gravitoinertial ray dynamics is governed by a Hamiltonian with two degrees of freedom, where the phase space is fourdimensional. The fact that ω = H is time independent is equivalent to energy conservation. It implies that phasespace trajectories actually stay on a threedimensional surface. A PSS is the intersection of all phasespace trajectories with a given threedimensional surface, defined for example by fixing one coordinate. A PSS at a given frequency is therefore a twodimensional surface.
Different choices are possible for the PSS, but to provide a complete view of phase space, it must be intersected by most phasespace trajectories. Here we used a PSS defined by fixing the colatitude θ to π/ 2, which corresponds to half of the equatorial plane. Furthermore, to ensure that two trajectories do not intersect on the same point on the PSS, the intersection with θ = π/ 2 is taken into account only when trajectories cross it from one particular side (here from the θ<π/ 2 side). For a given frequency ω, the PSS is computed by following many different trajectories over long time intervals. Then to scan the phase space at a given rotation, the PSS is computed for frequencies spanning the interval. For nonintegrable systems, the phasespace structure is complex and the features observed in a PSS are expected to depend on the resolution (in k_{r}, r or ω) by which the dynamics is investigated. However, modes occupy a finite phasespace volume because Fourier analysis shows that their wavevector localisation is inversely proportional to their spatial localisation (see Lignières & Georgeot 2009). Investigating phase space with a finite resolution is therefore sufficient to infer mode properties.
4. Gravitoinertial ray dynamics in rotating stars
In this section, we study the evolution of the gravitoinertial ray dynamics with rotation. We first describe the case of pure gravity rays in a nonrotating polytropic model (Sect. 4.1) and then consider gravitoinertial rays under the traditional approximation (Sect. 4.2). These are two simple integrable systems that serve as a reference to study the general case of gravitoinertial rays in a rotating star (Sect. 4.3).
4.1. The nonrotating case Ω = 0
In the nonrotating case, the eikonal equation simplifies into Eq. (21), where L is the second invariant that makes the Hamiltonian system with two degrees of freedom integrable. The phase space of integrable systems is made of nested invariant tori. These structures are said to be invariant because any trajectory that starts on one of the structures stays on it; they are called tori because they have the topology of a twodimensional torus. Any torus is specified by a value of the doublet (ω, L) and its intersection with the θ = π/ 2 hypersurface is the following f_{L2,ω}(k_{r},r) = 0 curve: (33)trivially derived from the eikonal Eq. (21).
As for any integrable system, two types of invariant tori exist. Irrational tori, on which any ray covers the whole torus, or equivalently fills the f_{L2,ω}(k_{r},r) = 0 curve on the PSS, and rational (or resonant) tori, on which any ray closes on itself before covering the torus. Rays of resonant tori are thus periodic orbits that imprint the PSS on a finite number of points.
Fig. 4 From panel a) to panel c), three PSS at ω/N_{0,max} = [ 0.66676,0.50200,0.05522 ] in a nonrotating polytropic model of star. The f_{L2,ω}(k_{r},r) = 0 curves shown are the imprint on the PSS of the invariant tori of the integrable gravity ray dynamics. They are drawn here for L = 1 / 2 + ℓ with ℓ taking integer values from 1 to 30 corresponding to degrees of spherical harmonics. The value of L of a given tori can be easily estimated from the intersection with the blue vertical dashed line, since the ordinate of the intersection point is approximately equal to 4L. The EBK quantisation enables us to specify the tori where positive interferences of rays produce modes. For example, the (n,ℓ) = (3,3) gravity mode on panel b) and the (n,ℓ) = (50,3) gravity mode on panel c) are constructed on the tori highlighted in green. The red curve on panel a) is a separatrix also called a homoclinic orbit because it connects the hyperbolic point at (k_{r},r/R) = (0,0.7968) to itself. As explained in Sect. 4.3.1, it is at the origin of a large chaotic phasespace region observed in the rotating case. 

Open with DEXTER 
Figure 4 presents three PSS computed at three different frequencies. For integrable systems, modes can be related to the dynamics through the EinsteinBrillouinKeller (EBK) quantisation procedure that has been applied to nonrotating spherical stars by Gough (1986; see also Gough 1993; Lignières & Georgeot 2009). Quantisation conditions actually select tori that support modes in the sense that rays belonging to these tori interfere positively to construct standing waves, i.e. modes. Accordingly, the degree of the spherical harmonic of the mode, ℓ, is related to the invariant L by the relation L = ℓ + 1 / 2. We therefore chose to represent in Fig. 4 tori with ℓ between 1 and 30. Two of the three PSS have been computed for frequencies corresponding to gravity modes: the (n,ℓ) = (3,3) and the (n,ℓ) = (50,3) modes. On these PSS, the torus associated with the mode is highlighted in green. These examples are useful to gain insight into the relation between the position on the PSS and the modes. In the same spirit, the k_{r} normalisation has been chosen in such a way that the value of L of a given k_{r} = f_{L2,ω}(r) curve can be easily retrieved. Indeed, from the eikonal Eq. (21) and provided that ω ≪ N_{0,max}, we obtain Rωk_{r}/N_{0,max} ≃ LR/r_{max}, where r_{max} = 0.28R is the radius of the inner maximum of N_{0}. Thus, the ordinate of the intersection of the k_{r} = f_{L2,ω}(r) curve with the r = 0.25R vertical line is approximatively equal to 4L.
4.2. The traditional approximation
The socalled traditional approximation (Eckart 1960; Lee & Saio 1987) consists in neglecting the contribution of the latitudinal component of the rotation vector in the Coriolis force. The solutions of the perturbation equation are then separable into the coordinates θ and r. This approximation is justified if motions are predominantly horizontal, the effect of the Coriolis force in the radial direction is negligible, and the centrifugal deformation is neglected.
Assuming a spherically symmetric stellar model (thus e_{∥} = e_{r} and e_{⊥} = e_{θ}) and neglecting the fsinθ terms, the eikonal Eq. (20) simplifies into (34)From this equation and the associated dynamical equations, we found a second invariant, namely (35)which implies that the associated ray dynamics is integrable. Integrability is expected since the perturbation equation in the traditional approximation is separable. Furthermore, neglecting the k_{c} term, the condition of propagation derived from Eq. (34) is (36)For superinertial waves, the domain of propagation is simply given by ω<N_{0}. This frequency range is more restricted than that given by Γ > 0, but both are qualitatively similar and even identical in the limit f/N_{0} → 0 if the centrifugal deformation is neglected (see Eq. (26)). For subinertial waves, the domain of propagation of the traditional approximation consists of two unconnected domains above and below the radius where ω = N_{0}. Above this radius, the equatorial region defined by fcosθ<ω and below, the polar cone defined by ω<fcosθ. Thus, the propagation of subinertial waves from the N_{0}>f layers towards the central region N_{0}<f is not allowed by the traditional approximation (for a graphical account of these differences in the domain of propagation, see Fig. 11 of Gerkema et al. 2008). This is a clear limitation of the traditional approximation since, as we have seen from the discussion on the domain of propagation, subinertial waves are indeed allowed to go through the centre of the star (and more generally through regions where N_{0}<f). Near the surface, the condition is with , which is a condition formally similar to the nonrotating case, with λ playing the role of L^{2}.
Fig. 5 Two PSS at ω/f = [ 2.21,0.567 ] of the ray dynamics in the traditional approximation in a spherical star rotating at Ω = 0.38Ω_{K}. The f_{λ,ω}(k_{r},r) = 0 curves shown are the imprint on the PSS of the invariant tori of the integrable dynamics. They are drawn for the quantised values of λ(ℓ,ω/f), where ℓ is the degree of the spherical harmonics of the corresponding mode at Ω = 0 and is varied from 1 to 30. The increase of the gap between the f_{λ,ω}(k_{r},r) = 0 curves as ω decreases is due to the increase of λ with f/ω. Tori in green correspond to gravitoinertial modes labelled by their quantum numbers at zero rotation, that is the (n,ℓ) = (3,3) mode (top) and the (n,ℓ) = (50,3) mode (bottom). 

Open with DEXTER 
The second expression in Eq. (35) yields the imprint of the invariant tori on the PSS. We observe that it has the form f_{λ,ω}(k_{r},r) = 0 that is the same as in the nonrotating case using λ instead of L^{2}. This relation between λ and L^{2} is the raydynamics version of the known relation between the separation constant of the perturbation equation in the traditional approximation and the product ℓ(ℓ + 1) (Bouabid et al. 2013). The PSS of the traditional approximation has thus the same structure as the nonrotating PSS. The effect of rotation comes from the increase of λ with f/ω, which produces the dilation effect observed between the two PSS shown in Fig. 5. As for the nonrotating case, gravitoinertial modes computed under the traditional approximation can be associated with specific tori. This is the case in Fig. 5, where two gravitoinertial modes labelled by their quantum numbers at zero rotation, namely (n,ℓ) = (3,3) and (n,ℓ) = (50,3), are represented in green on the PSS.
4.3. Phasespace structure of rotating stars
As rotation increases the gravitoinertial ray dynamics is no longer integrable and undergoes a transition towards a mixed state with coexistence of chaotic regions and invariant tori. This is apparent in Fig. 6, which shows the phasespace structure at Ω = 0.38Ω_{K} and ω/f = 3.1 together with a selection of rays in the physical space. As compared to the nonrotating PSS, which only exhibits f_{L2,ω}(k_{r},r) = 0 onedimensional curves, we now observe three different types of structures. First of all, there are onedimensional curves similar in shape to the nonrotating curves. The blue ray provides an example of ray associated with this type of tori. A second and new type of phasespace structures forms around stable periodic orbits and is called island chains. An orbit is said to be stable when trajectories induced by small perturbations of the orbit remain close to it. The green ray highlights one of the invariant tori formed around a central periodic orbit represented by the light green ray shown in the physical space. Finally, the third type of structure is the chaotic region, which can be easily identified on the PSS by the fact that the associated rays do no stay on a onedimensional curve but fill a twodimensional surface instead. The chaotic zone of Fig. 6 is highlighted by the imprint of a ray also shown in the physical space.
The apparition of these new structures in the gravitoinertial dynamics is typical of the KAMtype (in reference to the KolmogorovArnoldMoser theorem) transition of Hamiltonian systems from integrability towards chaos (Ott 1993). Accordingly, in systems that experience infinitesimal departures from integrability, all the resonant tori of the integrable system are destroyed and replaced by island chains around stable periodic orbits and small chaotic regions developing around unstable periodic orbits. Meanwhile, irrational tori are continuously perturbed but not destroyed. The evolution towards larger values of the parameter controlling the departure from integrability has been studied in many cases. The common phenomenology emerging from these studies is that, when the control parameter increases, the invariant tori, either the surviving irrational tori or the island chains, are progressively destroyed while chaotic regions occupy larger phasespace volumes. As demonstrated in the simple case of the kicked rotor (Ott 1993), for each irrational torus there is a critical value of the control parameter above which the torus is destroyed and these critical values serve to quantify the progression towards chaos.
Fig. 6 PSS computed at ω/f = 3.1 for a rotation Ω / Ω_{K} = 0.38 and four gravitoinertial rays in the physical space belonging to the three different types of phasespace structures: a surviving KAM torus (blue), a period5 island chain (green/light green), and a chaotic zone (red). The light green ray on the middle panel is the stable periodic orbit at the centre of the period5 island chain. The period of an island chain is defined as the number of corresponding structures on the PSS (see Sect. 4.3.2). These rays cross the PSS on points of corresponding colours. On all figures showing rays in the physical space, the limits of the domain of propagation given by Γ(ω) = 0 are drawn in magenta. 

Open with DEXTER 
The present gravitoinertial ray dynamics clearly undergoes a KAMtype transition towards a mixed state including surviving irrational tori, island chains, and chaotic regions. Nevertheless, we find that the actual departure from integrability strongly depends on the frequency ω. Indeed, the surface occupied by the island chains and chaotic regions strongly decreases towards low frequencies. This is obvious when comparing the phasespace structure of a Ω = 0.38Ω_{K} star at three frequencies, ω/f = { 3.1,2.4,0.8 }, as shown in Figs. 6–8, respectively. This observation does not mean that island chains and chaotic regions are absent below some frequency but rather that they are too small to be visible when the PSS are shown on large k_{r} and r scales. The same phenomenology also holds at higher rotation as shown by the PSS of a Ω = 0.84Ω_{K} star computed for three decreasing frequencies, ω/f = { 2,1.5,0.5 } and shown in Figs. 9–11, respectively. We therefore conclude that the departure from integrability as measured by the presence of large island chains and chaotic regions strongly diminishes in the low part of the interval and is weakly dependent on the centrifugal deformation.
Overall we may say that the tori that results from the continuous evolution of the nonrotating tori still shape the largescale organisation of the phase space even at large rotation rates. This is in contrast with the axisymmetric acoustic ray dynamics where, as rotation increases, the phase space is progressively dominated by large islands and a large chaotic region, while irrational tori only survive for large values of k_{θ}/ω corresponding to whispering gallery modes with low discintegrated visibilities.
Fig. 7 PSS computed at ω/f = 2.4 for a rotation Ω / Ω_{K} = 0.38 and two gravitoinertial rays belonging to a period1 island chain including the central stable periodic orbit (light blue). 

Open with DEXTER 
Fig. 8 PSS computed at ω/f = 0.8 for a rotation Ω / Ω_{K} = 0.38 and two gravitoinertial rays belonging to a period11 island chain (blue) and to a surviving KAM torus (green), respectively. 

Open with DEXTER 
In the following, the properties of each of the three types of phasespace structures are described in more detail. Considering them separately is relevant for the study of gravitoinertial modes because, according to Percival’s conjecture (Percival 1973; Berry & Robnik 1984) and following the results obtained for acoustic modes (Lignières & Georgeot 2009), we expect that each type of phasespace structure can be asymptotically related to families of modes with specific properties.
4.3.1. Chaotic regions
PSS computed at different rotation rates show that a large chaotic region is systematically present in the upper part of the frequency interval. In the following we argue that this chaotic region is related to the existence of unstable periodic orbits, which are themselves a direct consequence of the presence of a local minimum of the BruntVäisälä frequency in the stellar envelope. This interpretation is supported by our ability to predict the frequency range for which large chaotic regions are seen on the PSS.
In smoothly perturbed integrable Hamiltonian systems, chaotic motions appear near separatrices that connect hyperbolic points, that is unstable periodic orbits or unstable equilibrium points. In the framework of the KAM transition theory, these hyperbolic points naturally originate from the destruction of the resonant tori and chaos first appears around these hyperbolic points and along the structures that connect them (Ott 1993). Chaotic motions, however, can also originate from hyperbolic points already present in the unperturbed system. In our case, the unperturbed system does possess hyperbolic points as can be inferred from the PSS at Ω = 0 shown in panel (a) of Fig. 4. The vicinity of the point (k_{r},r) = (0,0.8) is indeed typical of the phase portrait around a hyperbolic point with nearby orbits forming branches of hyperbolas. The separatrix (the red curve in panel (a) of Fig. 4) is also called a homoclinic orbit because it connects the hyperbolic point to itself. The apex of the motion of a pendulum provides a classical example of a hyperbolic point in a system with one degree of freedom, which is in this case a hyperbolic fixed point or equivalently an unstable equilibrium point. In our system with two degrees of freedom, the hyperbolic point is rather an unstable periodic orbit because it corresponds to a circular trajectory at a fixed radius (r ~ 0.8 on panel (a) of Fig. 4). Since the ray dynamics in a spherical star is separable in the radial and latitudinal coordinates, the radius of the hyperbolic point is actually a fixed point of the radial ray dynamics. As in the case of the pendulum, the radial Hamiltonian takes the usual form and the unstable fixed point corresponds to the maximum of the potential V_{r}.
In Appendix E, analytical solutions of this problem are found via the envelope approximation, which involves the assumption that the mass enclosed in a sphere of given r does not depend on r, and we show that these solutions provide a very good approximation to the numerical solutions obtained without the envelope approximation. Accordingly, the gravity ray dynamics possesses a hyperbolic point for any frequency ω within the [ 0.93N_{0,min},N_{0,min} ] interval. The radius of the hyperbolic point varies monotonically from r = 0.84R for ω = 0.93N_{0,min} to r = 0.75R for ω = N_{0,min}, while the value of L characterising the separatrix varies also monotonically from L = 16.92 for ω = 0.93N_{0,min} to L = + ∞ for ω = N_{0,min}.
To confirm that chaos first occurs around the separatrix, we computed gravitoinertial rays close to this orbit for a slowly rotating model Ω / Ω_{K} = 0.02. As illustrated by Fig. 12, a small chaotic region is found in the expected frequency range and at the expected location. We further make the hypothesis that the large chaotic regions observed at Ω / Ω_{K} = 0.38 and Ω / Ω_{K} = 0.84 are due to the same mechanism. To test this hypothesis, PSS have been computed for frequencies inside and outside . To define the lower and upper bounds of this interval, rotational effects have been taken into account on phenomenological grounds. First, rotation changes the frequency range for propagation, and following the conditions (26), a f^{2}sin^{2}Θ term must be added to . Second, we included the interval opened by the centrifugal force. Indeed, we expect that unstable periodic orbits exist for these frequencies because there is always a range of latitudes where the condition is fulfilled. Considering frequencies 5% apart from the interval limits, we found that the large chaotic region is present inside and absent outside the interval.
Fig. 9 PSS computed at ω/f = 2.0 for a rotation Ω / Ω_{K} = 0.84 and two gravitoinertial rays belonging to a period9 island chain (blue) and to a chaotic zone (red), respectively. 

Open with DEXTER 
Fig. 10 PSS computed at ω/f = 1.5 for a rotation Ω / Ω_{K} = 0.84 and two gravitoinertial rays belonging to a period5 island chain including the central stable periodic orbit (light blue). 

Open with DEXTER 
Fig. 11 PSS computed at ω/f = 0.5 for a rotation Ω / Ω_{K} = 0.84 and a gravitoinertial ray belonging to an island chain of large (>26) period. 

Open with DEXTER 
Fig. 12 PSS computed at ω/f = 30.5 for a rotation Ω / Ω_{K} = 0.02 showing the development of a small chaotic zone in the vicinity of the separatrix present in the nonrotating case (red curve on panel a) of Fig. 4). 

Open with DEXTER 
At a given frequency, the extent of the chaotic surface seen on the PSS is not simple to predict quantitatively. Comparing Figs. 12, 6, and 9 we nevertheless observe that, as expected, the surface of the chaotic region increases with rotation.
4.3.2. Island chains
As they originate from the destruction of given rational tori, island chains can be associated with the rational number p/q of their parent torus. Then, the central periodic orbit of the p/q island imprints q times the PSS (Ott 1993); in this study we call the number q the period. The number p can be found by constructing a different PSS (for instance by fixing the radius r instead of the colatitude). For superinertial rays, p corresponds to the number of loops visible in real space. On the six PSS displayed in Figs. 6–11, we highlighted island chains of various periods: one period1 at ω/f = 2.4, two of period5 (but with different p values) at ω/f = 3.1 and ω/f = 1.5, respectively, one period9 at ω/f = 2, one period11 at ω/f = 0.8, and one with a period larger than 26 at ω/f = 0.5. As expected from theory, the width of the island decreases with their period. Moreover, there is a clear tendency to find only small and largeperiod island chains in the subinertial regime, although there is no obvious change of behaviour at ω/f = 1.
In a typical KAM transition, island chains form immediately for a nonzero pertubation, grow with the amplitude of the perturbation, and are finally destroyed. We followed the evolution of some lowperiod island chains with rotation and found that their width slightly increase but no destruction of island chains above a certain rotation rate has been observed. For example, the period1 island chain is still clearly visible at 99% of the critical velocity.
In the physical space, rays belonging to island chains do not fill their resonant cavity as they remain concentrated around the stable periodic trajectory. From Ballot et al. (2012), we already know that numerically computed modes can be associated with island chains. A detailed confrontation of the island properties with numerically computed modes will be performed in a future paper.
4.3.3. Surviving tori
The irrational tori of the nonrotating case are continuously deformed by the effect of rotation but may nevertheless have survived destruction even at high rotation rates. This is the case, in particular, for highL tori, since for most frequencies no large island chain or chaotic region is observed in the PSS in this domain, which corresponds to high values of the normalised k_{r}. Some irrational tori with lower L also survive the perturbation induced by rotation. This is obvious at low frequencies where phase space is very much structured by onedimensional curves broadly similar to the nonrotating curves (see Figs. 8 and 11 and panel (c) of Fig. 4).
This nearly integrable behaviour of the lowfrequency, gravitoinertial, ray dynamics might be attributed to the proximity of a hidden integrable system. To test whether the ray dynamics in the traditional approximation could be this integrable system, we computed the “traditional” invariant λ along the path of gravitoinertial rays. Considering a subinertial ray in a Ω / Ω_{K} = 0.38 star, we found that λ varies by several orders of magnitude along a given ray. As expected, the largest variations occur in regions of space that are outside the “traditional” domain of propagation, but even within the resonant cavity of the traditional approximation, λ still varies by a factor of two. In the superinertial regime, the variations of λ are of the order of 1.5, which is of the same order of magnitude as the variations of L; the parameter L is the invariant of the nonrotating case, which is not supposed to provide a valid approximation. We conclude that the ray dynamics computed in the traditional approximation does not provide an accurate approximation even in the lowfrequency regime, where the full gravitoinertial ray dynamics is close to integrable.
5. From rays to modes
In this section, we summarise the properties of gravitoinertial modes that can be inferred from the present study of the ray dynamics. Methods to construct modes from the superposition of positively interfering rays have been first proposed for integrable systems. The EBK procedure has been established in the field of quantum physics as a generalisation of the BohrSommerfeld method to quantise systems from their classical trajectories. This method only works in the presence of invariant tori and leads to a regular frequency spectrum in the sense that the spectrum can be described by a function of N integers, where N is the number of degrees of freedom of the system. For chaotic dynamical systems, the EBK quantisation does not work and must be replaced by the Gutzwiller trace formula, which is difficult to apply in practice (Gutzwiller 1990). Nevertheless, generic statistical properties of spectra associated with chaotic dynamical systems, such as the distribution of spacings between consecutive frequencies, have been discovered in quantum systems (Bohigas et al. 1984). For mixed dynamical systems, it has been conjectured that chaotic regions and phasespace regions structured by invariant tori, such as island chains or regions of surviving irrational tori, can be quantised separately. Consequently, the spectrum of mixed systems is expected to be a superposition of frequency subsets, each of which are associated with different phasespace structures. Frequency subsets associated with regions structured by invariant tori should be regular while those associated with chaotic regions should possess generic statistical properties (Percival 1973; Berry & Robnik 1984). A more detailed description of the concepts presented above together with the full references is available in Lignières & Georgeot (2009). Percival’s predictions have been applied to acoustic rays in rotating stars and successfully confronted with numerically computed highfrequency acoustic modes (Lignières & Georgeot 2009).
In our case, Percival’s conjecture implies that the frequency spectrum of shortwavelength gravitoinertial modes is a superposition of frequency subsets associated with the phasespace structures of the gravitoinertial ray dynamics, notably a regular subspectrum associated with the surviving irrational tori, regular subspectra associated with island chains, and a subspectrum with generic statistical properties associated with the large chaotic region. In principle, quantitative predictions on the regular subspectra can be obtained by conducting the EBK quantisation of the invariant tori. This has been carried out in Pasek et al. (2012) for the period2 island chain of the acoustic ray dynamics and could be also attempted for the gravitoinertial island chains and surviving tori.
Previous numerical explorations of gravitoinertial modes already provide some insights into the relation between the asymptotic predictions of the gravitoinertial dynamics and the numerically computed modes. Using μ = 3 polytropic models, Ballot et al. (2010) followed lowdegree (ℓ ≤ 3) modes in both low and highorder ranges from Ω = 0 to Ω = 0.7Ω_{K}. The vast majority of these modes undergo a smooth evolution of their spatial distribution with rotation, in which the main change is the equatorial mode confinement in the subinertial regime. These modes most probably correspond to the phasespace region structured by the surviving irrational tori. On the other hand, some modes experience a striking modification of their spatial distribution, leading to a concentration of their energy. The connection between these modes called rosette modes and the central periodic orbits of island chains has been clearly established in Ballot et al. (2012). Since then, Takata & Saio (2013) proposed a different approach in which rosette modes result from the rotationally induced coupling of modes, which were nearly degenerate gravity modes at Ω = 0 (see also Saio & Takata 2014; Takata 2014). The link between these two apparently very different approaches to model rosette modes remains to be understood. In our study, the number K introduced by Takata & Saio (2013) to characterise families of rosette modes corresponds to the ratio 2q/p.^{2} The authors found only integer values of K, which they justified by the fact that the only effect of rotation they considered was the Coriolis force. We also present island chains with noninteger values. This is probably because we took the centrifugal deformation into account. Modes that could be related to the large chaotic region that we describe have not been identified yet. By providing an estimate of both the frequencies and wave vectors of chaotic rays, ray dynamics will help to find these modes. More generally, a dedicated numerical study will be necessary to test the predictions of the raybased asymptotic theory. In this context, tools allowing us to construct phasespace representations for numerically computed modes, such as Husimi distributions, are available. As illustrated by Lignières & Georgeot (2009), those can be used to show the imprint of modes on the PSS.
6. Discussion and conclusions
The main results presented in the previous sections are (i) the derivation of an eikonal equation for gravitoinertial rays in polytropic uniformly rotating stars; (ii) the analysis of the properties of this equation, including constraints on the domains of propagation for k_{φ} = 0 rays; (iii) the numerical computation of the k_{φ} = 0 ray dynamics for μ = 3 polytropic models with increasing rotation, and (iv) the interpretation of the numerical results using tools and concepts of Hamiltonian dynamics. In particular, the exploration of the phase space provides us with a global view of the properties of gravitoinertial waves in rotating stars. Qualitative predictions on the frequency spectrum of gravitoinertial modes have been made and should be a powerful tool to interpret the properties of numerically computed gravitoinertial modes. A quantitative result already provided by the ray theory is the domain of propagation. Other quantitative predictions are expected from the EBK quantisation of the surviving irrational tori and of the island chains.
A first restriction of the present work is that we did not consider k_{φ} ≠ 0 rays. This would be necessary for an asymptotic study of nonaxisymmetric gravitoinertial modes. Another restriction is the use of μ = 3 polytropic models although the formalism can be extended to more realistic rotating stellar models. In particular, the effects of a convective core and/or a convective envelope would be interesting to consider.
The present ray dynamics formalism can be improved by taking the effect of background mean motions and, in particular, differential rotation into account. For example, selfconsistent models including baroclinic effects computed by the ESTER code (Rieutord & Espinosa Lara 2013) could be used. According to atmospheric or oceanic studies, the eikonal equation including mean horizontal motions is identical to (18), except that the frequency is changed into the Dopplershifted frequency, which is the frequency that would be observed in a frame of reference moving with the mean horizontal motions (Fritts & Alexander 2003). Analagous with the effect of wind shear, differential rotation should cause some backrefraction and thus will modify the resonant cavity. Another improvement to the eikonal equation would be to include nonaxisymmetric Rossby waves.
Ray models might also be used to study the transport of angular momentum by waves. This would require us to extend the ray formalism to treat the evolution of the wave amplitude A along the ray path (the next order of the WKB approximation). Such amplitude equations are already used to investigate the driving of the atmospheric circulation by waves (Marks & Eckermann 1995). To conduct similar investigations in stars, it will be necessary to model energy deposition processes such as thermal dissipation and critical layers (Fritts et al. 1998; Alvan et al. 2013) and to take phase changes at caustics into account (Broutman et al. 2004).
Angular momentum transport by gravitoinertial waves that occasionally escape the star has been invoked as a possible explanation of the Be phenomenon (Ishimatsu & Shibahashi 2013). Specific properties of chaotic gravitoinertial waves such as those observed in our calculations could be interesting in this context. At any location and in particular at the surface, the horizontal wavelength k_{⊥} of chaotic waves can indeed take many different values. Thus a chaotic wave can be backrefracted many times at the surface before it escapes the star because k_{⊥} has reached a high enough value.
When compared to the ray theory of Dintrans & Rieutord (2000), the present work includes the effects of centrifugal deformation as well as a realistic backrefraction of rays in the stellar envelope. This last point has a crucial effect on the comparison between the two approaches. In Dintrans & Rieutord (2000), the rays of characteristic are confined into the star by a rigid and stressfree surface and this boundary condition destroys the Hamiltonian character of the dynamics of characteristics. As a consequence, the volume of phasespace elements is not conserved, which in turns enables rays to focus on attractors. This focusing effect is at the origin of the singular gravitoinertial modes that are effectively observed in the presence of rigid boundaries (Maas et al. 1997; Dintrans et al. 1999). However, there are no such rigid boundaries in real stars and singular gravitoinertial modes focused on attractors of characteristics by unrealistic rigid boundaries should not exist in stars. As we have seen in the present paper, the ray dynamics ought to be Hamiltonian in stars as long as dissipative processes are neglected.
Acknowledgments
The authors thank Michel Rieutord for his careful reading of the manuscript, and an anonymous referee for contributing to the improvement of this paper. This work was supported by the Centre National de la Recherche Scientifique (CNRS) through the Programme National de Physique Stellaire (PNPS).
References
 Aerts, C., Kurtz, D., & ChristensenDalsgaard, J. 2010, Asteroseismology (Springer) [Google Scholar]
 Alvan, L., Mathis, S., & Decressin, T. 2013, A&A, 553, A86 [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]
 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 [Google Scholar]
 Ballot, J., Lignières, F., & Reese, D. R. 2013, in Numerical Exploration of Oscillation Modes in Rapidly Rotating Stars (Springer), Lect. Notes Phys., 865, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Berry, M. V., & Robnik, M. 1984, J. Phys. A Math. Gen., 17, 2413 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Bohigas, O., Giannoni, M. J., & Schmit, C. 1984, Phys. Rev. Lett., 52, 1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Bouabid, M.P., Dupret, M.A., Salmon, S., et al. 2013, MNRAS, 429, 2500 [NASA ADS] [CrossRef] [Google Scholar]
 Broutman, D., Rottman, J. W., & Eckermann, S. D. 2004, Ann. Rev. Fluid Mech., 36, 233 [NASA ADS] [CrossRef] [Google Scholar]
 ChristensenDalsgaard, J., & Thompson, M. J. 1999, A&A, 350, 852 [NASA ADS] [Google Scholar]
 Cox, J. P. 1968, Principles of stellar structure (New York: Gordon and Breach) [Google Scholar]
 Dintrans, B., & Rieutord, M. 2000, A&A, 354, 86 [NASA ADS] [Google Scholar]
 Dintrans, B., Rieutord, M., & Valdettaro, L. 1999, J. Fluid Mech., 398, 271 [NASA ADS] [CrossRef] [Google Scholar]
 Eckart, C. 1960, Hydrodynamics of oceans and atmospheres (New York: Pergamon Press) [Google Scholar]
 Fritts, D. C., & Alexander, M. J. 2003, Rev. Geophys., 41, 1003 [NASA ADS] [CrossRef] [Google Scholar]
 Fritts, D. C., Vadas, S. L., & Andreassen, O. 1998, A&A, 333, 343 [NASA ADS] [Google Scholar]
 Gerkema, T., Zimmerman, J. T. F., Maas, L. R. M., & van Haren, H. 2008, Rev. Geophys., 46, 2004 [NASA ADS] [CrossRef] [Google Scholar]
 Gough, D. 1993, in Astrophysical Fluid Dynamics – Les Houches 1987, 399 [Google Scholar]
 Gough, D. O. 1986, in Hydrodynamic and Magnetodynamic Problems in the Sun and Stars, ed. Y. Osaki, 117 [Google Scholar]
 Gutzwiller, M. C. 1990, Chaos in classical and quantum mechanics (Springer) [Google Scholar]
 Hansen, C. J., & Kawaler, S. D. 1994, Stellar Interiors. Physical Principles, Structure, and Evolution (Berlin: SpringerVerlag) [Google Scholar]
 Harlander, U., & Maas, L. R. M. 2006, Meteorol. Z., 15, 439 [CrossRef] [Google Scholar]
 Ishimatsu, H., & Shibahashi, H. 2013, in Progress in Physics of the Sun and Stars: A New Era in Helio and Asteroseismology, eds. H. Shibahashi, & A. E. LynasGray, ASP Conf. Ser., 479, 325 [Google Scholar]
 Kurtz, D. W., Saio, H., Takata, M., et al. 2014, MNRAS, 444, 102 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, U., & Saio, H. 1987, MNRAS, 224, 513 [NASA ADS] [CrossRef] [Google Scholar]
 Lighthill, J. 1978, Waves In Fluids (Cambridge University Press) [Google Scholar]
 Lignières, F. 2011, in The Pulsations of the Sun and the Stars (Springer), Lect. Notes Phys., 832, 259 [NASA ADS] [CrossRef] [Google Scholar]
 Lignières, F. 2013, in Stellar Pulsations: Impact of New Instrumentation and New Insights, eds. J. C. Suárez, R. Garrido, L. A. Balona, & J. ChristensenDalsgaard, Astrophys. Space Sci. P., 31, 43 [Google Scholar]
 Lignières, F., & Georgeot, B. 2008, Phys. Rev. E, 78, 016215 [NASA ADS] [CrossRef] [Google Scholar]
 Lignières, F., & Georgeot, B. 2009, A&A, 500, 1173 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maas, L. R. M., Benielli, D., Sommeria, J., & Lam, F.P. A. 1997, Nature, 388, 557 [NASA ADS] [CrossRef] [Google Scholar]
 Marks, C. J., & Eckermann, S. D. 1995, J. Atm. Sci., 52, 1959 [NASA ADS] [CrossRef] [Google Scholar]
 Miglio, A. 2015, IAU General Assembly, 22, 51618 [NASA ADS] [Google Scholar]
 Ott, E. 1993, Chaos in dynamical systems (Cambridge University Press) [Google Scholar]
 Ouazzani, R.M., Dupret, M.A., & Reese, D. R. 2012, A&A, 547, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pasek, M., Georgeot, B., Lignières, F., & Reese, D. R. 2011, Phys. Rev. Lett., 107, 121101 [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]
 Percival, I. C. 1973, J. Phys. B At. Mol., 6, L229 [NASA ADS] [CrossRef] [Google Scholar]
 Reese, D. R. 2015, in Eur. Phys. J. Web Conf., 101, 5007 [Google Scholar]
 Reese, D., Lignières, F., & Rieutord, M. 2006, A&A, 455, 621 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reese, D. R., MacGregor, K. B., Jackson, S., Skumanich, A., & Metcalfe, T. S. 2009, A&A, 506, 189 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rieutord, M., & Espinosa Lara, F. 2013, in Lect. Not. Phys. 865 (Berlin: Springer Verlag), eds. M. Goupil, K. Belkacem, C. Neiner, F. Lignières, & J. J. Green, 49 [Google Scholar]
 Rieutord, M., Corbard, T., Pichon, B., Dintrans, B., & Lignières, F. 2005, in SF2A2005: Semaine de l’Astrophysique Francaise, eds. F. Casoli, T. Contini, J. M. Hameury, & L. Pagani, 759 [Google Scholar]
 Royer, F., Zorec, J., & Gómez, A. E. 2007, A&A, 463, 671 [NASA ADS] [CrossRef] [EDP Sciences] [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]
 Takata, M. 2014, PASJ, 66, 80 [NASA ADS] [Google Scholar]
 Takata, M., & Saio, H. 2013, PASJ, 65, 68 [NASA ADS] [Google Scholar]
 Townsend, R. 2014, in IAU Symp. 301, eds. J. A. Guzik, W. J. Chaplin, G. Handler, & A. Pigulski, 153 [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]
Appendix A: Derivation of the perturbation Eq. (11)
In this appendix, our goal is to transform the perturbation Eqs. (5)–(7) into Eq. (11) for . In Sect. A.1, density and velocity perturbations appearing in Eqs. (5)–(7) are eliminated to obtain an equation for pressure perturbations only. Then in Sect. A.2, the new variable is introduced to eliminate the firstorder derivatives in the meridional plane from the perturbation equation on . Finally, in Sects. A.3 and A.4, the constant term is expanded in powers of to determine the dominant approximation of .
Appendix A.1: Derivation of an equation for pressure perturbations only
We first rewrite the perturbation Eqs. (5)–(7), using u = ρ_{0}v and the definition of the BruntVäisälä frequency as follows: In order to eliminate u from these equations, we first calculate f ∧() as (A.4)The first term on the LHS is then developed using ∂_{t}() to give (A.5)We use a similar method to eliminate the scalar product f·u. We compute f·()(A.6)and introduce it into ∂_{t}() to get (A.7)Defining the operator (A.8)enables us to rewrite Eqs. (A.1) and (A.3) as where ℒ(u) is given by Eq. (A.7) and α, defined by (A.11)is uniform in a polytropic model. Replacing ℒ(u) in Eq. (A.10) by the RHS of Eq. (A.7), we get (A.12)where (A.13)To simplify calculations, we assume that P and ρ are harmonic in time, i.e. and . Then the operator is simply a multiplication by (A.14)and ∂_{t}ℒ a multiplication by (A.15)From Eq. (A.12), is expressed as a function of as follows: (A.16)with (A.17)We then compute ∇·ℒ(u) from Eq. (A.7), using the general identity ∇·(ba) = a·∇b + b∇·a and the fact that f is uniform (A.18)with (A.19)The expression (A.18) can be simplified using (A.20)which finally yields (A.21)Replacing the expression (A.21) in Eq. (A.9) and using Eq. (A.16) to express , we obtain the following equation for : (A.22)We use the identity ∇·(a ∧ b) = b·(∇ ∧ a)−a·(∇ ∧ b), which implies that ∇·(f ∧ g_{0}) = 0, as f is uniform and g_{0} derives from a potential. Defining the vector F_{0} = (f·g_{0})f−ω^{2}g_{0}, the vectors F and H are written as F = F_{0} + iωf ∧ g_{0} and H = α(F_{0}−iωf ∧ g_{0}). Thus, Eq. (A.22) becomes (A.23)where we have used the fact that α is uniform and f ∧ g_{0} is parallel to e_{φ}. We then use the relation (A.24)to get (A.25)where From the definition of F_{0}, we have (A.28)where, thanks to the fact that f is uniform, the last term can be also written with (A.29)In the end, we get (A.30)which is an equation for a single unknown, as we were looking for.
Using the relation , from the definition of D, Eq. (A.30) is rewritten as (A.31)where (A.32)To greatly simplify the terms involving secondorder derivatives, one can express in the (nonorthonormal) basis defined by f, g_{0}, and e_{φ} and take the divergence of the resulting equation. It gives (A.33)where we define the following operators and expressions: with p = f·g_{0}, , n = f^{2} and . If in addition, the operators ∇_{z} and Δ_{⊥} are defined as in Eqs. (9) and (10), the wave Eq. (A.31) becomes (A.38)where, using the fact that f ∧ g_{0} = −fg_{0}sinΘe_{φ},
Appendix A.2: Normal form
Here, Eq. (A.41) for pressure perturbations is written in a normal form by eliminating firstorder derivatives in the meridional plane. This can be done by introducing a new variable defined by , where the axisymmetric function a is chosen to eliminate the firstorder term . To carry out this substitution, we use the following relations: where the operators ∇_{∥} and ∇_{z} are defined by Eqs. (9) and (10). Equation (A.38) becomes (A.45)where The condition translates into a firstorder linear differential equation for a. The projection of this equation on e_{∥} and e_{z} leads to two equations for ∇_{∥}a/a and ∇_{z}a/a. These equations can be solved as a linear system for these two unknowns provided that the determinant of the system, which is equal to D, does not vanish. We then use these expressions of ∇_{∥}a/a and ∇_{z}a/a to write (A.48)where where and are defined by . With this choice for the function a, the wave equation Eq. (A.45) simplifies into (A.52)where the term is now expressed as a function of the dimensionless quantity , defined by (A.53)The term equals M′′/a and can be expressed in terms of to be written as (A.54)From the definitions of M′ and M, it is easy to see that can be written as to recover Eq. (11).
Appendix A.3: The constant term
To take into account the backrefraction of outgoing waves in the shortwavelength approximation, is expanded in powers of H = ℛT/g_{0} or equivalently and we only retain the dominant term. Among the quantities involved in , the inverse of the sound speed and the BruntVäisälä frequency tend to become very large towards the stellar surface, whereas f,g_{0}, and Γ_{1} remain finite. To derive the dominant term of approaching the stellar surface, we use the following expression for : (A.55)where (A.56)and w_{0} = w_{0}e_{∥} and w_{1} are both with respect to the expansion.
The dominant terms in the expression (A.54) of are
whereas is negligible. Summing these terms we obtain, in the K ≠ 0 case, (A.62)As shown below in Sect. A.4, the calculation of w_{0} leads to (A.63)thus to (A.64)In the subinertial regime (ω<f), the dominant term of vanishes at the critical angles Θ_{c} and π−Θ_{c} such that cosΘ_{c} = ω/f. Coming back to the full expression of , Eq. (A.54), it can be shown that (A.65)and that C_{1} does not vanish when K = 0. Thus, as one approaches the surface along the critical angles, the dominant term in the expression of is given by (A.66)In practice, we found that rays do not reach the surface layers near the critical angle. This can be understood as the dispersion relation and the ray equations show that near the critical angle, k_{⊥} vanishes, and this implies that dr(θ) / dθ vanishes as well.
Appendix A.4: The w_{0} term
According to Eqs. (A.55) and (A.49), w_{0} is the dominant term of . From the definitions of and , we obtain (A.67)This expression is then developed using Thus, as vanishes, (A.75)The expression of w_{0} becomes where (A.78)Developing the last expression, we show that (A.79)Thus,
Appendix B: The ray dynamics equations
In this section, the ray dynamics equations for axisymmetric rays are written using the spherical coordinates [ r,θ ] and the wavevector components [ k_{r},k_{θ} ] on the usual orthonormal basis e_{r},e_{θ}. To derive them, it is convenient to start from the general Hamiltonian Eqs. (28) and (29) for spherical coordinates [ r,θ ]. These equations are expressed in terms of the covariant components of k on the natural basis , and a change of variables from to [ r,θ,k_{r},k_{θ} ] is necessary. From the relations and , Eqs. (28) and (29) take the following form: where the Hamiltonian H(r,θ,k_{r},k_{θ}) is given by the eikonal equation written in the form ω = H(r,θ,k_{r},k_{θ}).
To develop the previous equations, one possible way is to start from the original eikonal Eq. (18) in the case k_{φ} = 0 written as (B.5)and to calculate the partial derivatives of this expression with respect to the coordinates [ r,θ,k_{r},k_{θ} ]. To do so, one needs to express k_{z} and k_{⊥} as functions of [ r,θ,k_{r},k_{θ} ], namely where the last relation has been obtained from and e_{s} is the classical cylindrical unit vector pointing away from the zaxis.
The dynamical Eqs. (B.1)–(B.4) in spherical coordinates are These equations have been used for the numerical computation of the gravitoinertial dynamics in barotropic models of uniformly rotating stars.
Appendix C: Derivation of relations (30), (31), and (32)
It is also useful to derive other relations to understand ray dynamics. Below, the relation (30) between the group velocity and the phase velocity is derived. Using the first two equations of the Hamiltonian system (B.1) and (B.2), we have (C.1)Then, taking the derivatives of the eikonal Eq. (B.5) with respect to k_{r} and k_{θ}, we obtain It follows that (C.4)Then, from Eqs. (B.6) and (B.7), we have
and from there we deduce that which corresponds to relation (30).
Similar calculations enable us to derive the expression (31). Indeed, using the first two equations of the Hamiltonian system (B.1) and (B.2) and the expression of e_{∥} in Eq. (B.8), we have (C.11)Then, from Eqs. (C.2), (C.3), and (C.5)–(C.8), we find (C.12)Solving the eikonal equation (20) as a quadratic equation for the variable k_{∥}, we find (C.13)where δ has been already defined by Eq. (22). Then, the relation (31) is obtained replacing k_{z} by k_{z} = k_{∥}cosΘ−k_{⊥}sinΘ in Eq. (C.12). From the difference of Eq. (30) and Eq. (31), an expression for dx/ dt·e_{⊥} is obtained. Using the above expression of and Eq. (22), we derive Eq. (32).
Appendix D: Test of the numerical method on a regular cusp point
We considered the simple case of gravity rays in a planeparallel atmosphere characterised by an exponential variation of the BruntVäisälä frequency in the z direction. The ray dynamics equations derived from the dispersion relation are where k_{x} is constant, k_{z} = k_{z0}−γωt, and γ = dlnN_{0}/ dz. An analytical solution of these equations going through a regular cusp point of coordinates (x_{r},z_{r}) is (D.3)where the time variable t has been replaced by u = k_{z}/k_{x}. The cusp point is reached at u = 0. The code used in the present paper was able to reproduce this analytical solution with a controlled level of accuracy.
Appendix E: Unstable fixed point of the nonrotating radial ray dynamics
In nonrotating stars, the ray dynamics equations can be solved separately in the radial and latitudinal coordinates. The radial dynamics is governed by Eq. (21). For some values of the frequency, the phase portrait of the radial dynamics is typical of doublewellpotential systems with two elliptic regions around the two fixed points at the minima of the potential separated from highenergy motions by a separatrix that goes through the unstable hyperbolic fixed point at the maximum of the potential. This is illustrated by panel (a) of Fig. 4. Equation (21) can be seen as a classical Hamiltonian with H_{r} = 0 and V_{r} a potential that depends on two parameters, ω and L. In our case, the separatrix separates lowL from highL motions. When such a system is moved away from integrability, chaos appears first near the unstable hyperbolic point. Here, we show that using an envelope model we can analytically determine the frequency range where the nonrotating radial dynamics shows a doublewellpotential behaviour as well as the separatrix and the position of the unstable point.
In envelope models of stars, the radial dependence of the mass inside a given radius is neglected. Within this approximation, the hydrostatic equation can be integrated to provide simple expressions for the profile of the thermodynamic quantities (Cox 1968). In particular, the BruntVäisälä frequency can be written as (E.1)and can be shown to possess a local minimum N_{0,min} at r/R = 0.75. This local minimum is expected to be a generic feature of radiative stellar envelopes. It is at the origin of the existence of an unstable fixed point of the ray dynamics.
Using Eq. (E.1), Eq. (21) governing the radial dynamics becomes (E.2)where we have introduced the dimensionless quantities , x = r/R and the notations a = α(μ + 1) / Γ_{1} and b = (μ + 1)(μ−1) / 4. In our onedimensional system, k_{r} and the radial derivative of V_{r} must vanish at fixed points. As H_{r} = 0, the first condition implies that V_{r} also vanishes at the fixed point. The equilibrium will be unstable if V_{r} is maximum there. Asking that both V_{r} and its radial derivative vanish is equivalent to solve the following system of equations for the variable x excluding the x = 0 case: (E.3)Using two new variables u = 1−x and v = x^{3}, the system above is equivalent to solving a cubic equation for u, a second equation giving the relation between and L as follows: (E.4)The cubic equation has a unique real solution if L is smaller than a minimum value L_{min}, and this solution is not physical as it gives a point outside the star. Above L_{min}, there are in addition two real solutions among which the smaller solution corresponds to the unstable point and the larger solution to the outer stable point. Simple analytical solutions are obtained for the two extreme cases L = L_{min} and L = + ∞.
In conclusion, an unstable hyperbolic point is present in the frequency range 0.9348N_{0,min} ≤ ω ≤ N_{0,min}. As ω increases in this interval, the position of the unstable point decreases from r = 0.8406R to r = 0.75R, while the L value of the separatrix increases from to L = + ∞. It is interesting that some properties of the unstable hyperbolic point, namely the radius and the frequency interval normalised by N_{0,min}, do not depend on the polytropic index. For a μ = 3 polytropic model, L_{min} = 16.92.
We compared these analytical solutions to numerical determinations of the hyperbolic point obtained for the full μ = 3 polytropic model, that is without the envelope approximation. The agreement is fairly good since the numerical solution gives c_{2}N_{0,min} ≤ ω ≤ N_{0,min} with 0.9371 <c_{2}< 0.9372 and 16.80 <L_{min}< 16.81.
All Figures
Fig. 1 Top: solid black line shows the profile of the BruntVäisälä frequency N_{0} normalised by ω_{0} for the nonrotating stellar model. The red lines show N_{0} for the Ω = 0.84Ω_{K} model along the polar (solid line) and equatorial (dashed line) radii. A thick red tick on the xaxis indicates the polar radius R_{p}. Top: map in a meridional plane of N_{0} normalised by ω_{0} for the Ω = 0.84Ω_{K} model. The black line is the star surface. 

Open with DEXTER  
In the text 
Fig. 2 Typical domains of propagation (grey areas) given by the Γ(ω) ≥ 0 condition as a function of the gravitoinertial wave frequency ω in a Ω = 0.38Ω_{K} polytropic model of star. In general, the outer limit of the domain, here loosely represented by the stellar surface, cannot be predicted from ω alone and is obtained by solving the ray dynamics. Panel a) corresponds to a superinertial frequency (ω>f) and panel b) to a subinertial frequency (ω<f). The two bottom figures show domains with an intermediate forbidden zone when frequencies are between the local minimum and maximum of the BruntVäisälä frequency. Panel c) corresponds to frequencies in the interval and also exists in nonrotating stars, whereas panel d) corresponds to frequencies in the interval and is due to the centrifugal deformation of rotating stars. 

Open with DEXTER  
In the text 
Fig. 3 Gravitoinertial ray showing the two types of turning points, a cusp point on the Γ = 0 boundary (in magenta) and a smooth backrefraction in the surface layers where the decreasing density scale height becomes of the order of the wavelength in that direction. 

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

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

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

Open with DEXTER  
In the text 
Fig. 7 PSS computed at ω/f = 2.4 for a rotation Ω / Ω_{K} = 0.38 and two gravitoinertial rays belonging to a period1 island chain including the central stable periodic orbit (light blue). 

Open with DEXTER  
In the text 
Fig. 8 PSS computed at ω/f = 0.8 for a rotation Ω / Ω_{K} = 0.38 and two gravitoinertial rays belonging to a period11 island chain (blue) and to a surviving KAM torus (green), respectively. 

Open with DEXTER  
In the text 
Fig. 9 PSS computed at ω/f = 2.0 for a rotation Ω / Ω_{K} = 0.84 and two gravitoinertial rays belonging to a period9 island chain (blue) and to a chaotic zone (red), respectively. 

Open with DEXTER  
In the text 
Fig. 10 PSS computed at ω/f = 1.5 for a rotation Ω / Ω_{K} = 0.84 and two gravitoinertial rays belonging to a period5 island chain including the central stable periodic orbit (light blue). 

Open with DEXTER  
In the text 
Fig. 11 PSS computed at ω/f = 0.5 for a rotation Ω / Ω_{K} = 0.84 and a gravitoinertial ray belonging to an island chain of large (>26) period. 

Open with DEXTER  
In the text 
Fig. 12 PSS computed at ω/f = 30.5 for a rotation Ω / Ω_{K} = 0.02 showing the development of a small chaotic zone in the vicinity of the separatrix present in the nonrotating case (red curve on panel a) of Fig. 4). 

Open with DEXTER  
In the text 