Issue 
A&A
Volume 631, November 2019



Article Number  A140  
Number of page(s)  17  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201936459  
Published online  07 November 2019 
Regularities in the spectrum of chaotic pmodes in rapidly rotating stars
^{1}
Institut de Recherche en Astrophysique et Planétologie, Université de Toulouse, CNRS, CNES, UPS, France
email: francois.lignieres@irap.omp.eu
^{2}
Laboratoire de Physique Théorique, IRSAMC, Université de Toulouse, CNRS, UPS, France
Received:
5
August
2019
Accepted:
20
September
2019
Context. Interpreting the oscillations of massive and intermediate mass stars remains a challenging task. In fast rotators, the oscillation spectrum of pmodes is a superposition of subspectra which corresponds to different types of modes, among which island modes and chaotic modes are expected to be the most visible. This paper is focused on chaotic modes, which have not been thoroughly studied before.
Aims. We study the properties of high frequency chaotic pmodes in a polytropic model. Unexpected peaks appear in the frequency autocorrelations of the spectra. Our goal is to find a physical interpretation for these peaks and also to provide an overview of the mode properties.
Methods. We used the 2D oscillation code “TOP” to produce the modes and acoustic ray simulations to explore the wave properties in the asymptotic regime. Using the tools developed in the field of quantum chaos (or wave chaos), we derived an expression for the frequency autocorrelation involving the travel time of acoustic rays.
Results. Chaotic mode spectra were previously thought to be irregular, that is, described only through their statistical properties. Our analysis shows the existence, in chaotic mode spectra, of a pseudo large separation. This means that chaotic modes are organized in series, such that the modes in each series follow a nearly regular frequency spacing. The pseudo large separation of chaotic modes is very close to the large separation of island modes. Its value is related to the sound speed averaged over the meridional plane of the star. In addition to the pseudo large separation, other correlations appear in the numerically calculated spectra. We explain their origin by the trapping of acoustic rays near the stable islands.
Key words: asteroseismology / waves / chaos / stars: oscillations / stars: rotation
© B. Evano et al. 2019
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
Despite the many successful advancements made in stellar seismology, we are still unable to unlock most of the information contained in the acoustic oscillations of typical massive and intermediate mass stars. The observed pulsational behavior of these stars in the range of acoustic frequencies is diverse and has so far resisted a meaningful empirical classification. In addition to evolutionary effects, both the nonlinear processes that set the mode amplitudes and the sensitivity of the spectrum organization to the unknown rotation rate contribute to the observed diversity (Bowman & Kurtz 2018).
Progress in modeling the rotational effects together with the flow of high quality data from ultra precise space photometry missions, nevertheless, revive the interest for the seismology of massive and intermediate mass stars. An important step was the detection of regular frequency patterns analogous to the solarlike large separation in δ Scuti stars (García Hernández et al. 2009, 2013, 2015; Paparó et al. 2016; Michel et al. 2017). These patterns were predicted from the first oscillation model, that took realistic centrifugal distortion into account (Lignières et al. 2006), up to the more realistic models to date (Reese et al. 2017). Further progress is expected from TESS and PLATO missions that will include bright stars that can be better constrained through spectroscopy.
To model the acoustic spectrum of rapidly rotating stars, 2D oscillation codes that take full account of the effect of rotation on the oscillations were developed (Reese et al. 2006, 2009; Ouazzani et al. 2012). They can be run for different models of centrifugally deformed stars, from polytropic models (Reese et al. 2006) to more sophisticated ones (Reese et al. 2009; Ouazzani et al. 2015). Exploring the acoustic spectrum as a function of the star rotation is not a straightforward process, in particular, because the density of frequencies increases with the numerical resolution and the mode classification can be tedious (Ballot et al. 2013). First attempts were restricted to low sphericalharmonic degree modes, which were followed by progressively increasing the rotation rate (Lignières et al. 2006; Reese et al. 2006, 2008; Pasek et al. 2012). More complete spectra have then been obtained at fixed rotation rates (Lignières & Georgeot 2009; Reese et al. 2009, 2017; Ouazzani et al. 2015). An automatic classification of the computed mode using neural network methods has recently been proposed (Mirouh et al. 2019).
Understanding the spectrum organization is key to constructing seismic tools, the detection of regular patterns in δ Scuti stars being a good example. The structure of the spectrum is expected to best reveal itself in the asymptotic regime, that is, at highfrequencies for acoustic modes. This motivates highfrequency computations even if massive and intermediate mass stars oscillate at lower acoustic frequencies. The asymptotic regime is also amenable to theoretical descriptions. In the short wavelength “WKB” limit, acoustic waves are described by rays whose propagation obeys Hamiltonian dynamics. In Lignières & Georgeot (2009), the acoustic ray dynamics of rapidly rotating stars was studied and semiclassical concepts and methods developed in quantum physics were used to infer the asymptotic properties of the associated acoustic spectrum. Reflecting the phase space structure, different families of modes were identified. Above some rotation, two families are most likely to be observed: the 2period island modes and the chaotic modes. The first family was studied in detail through numerical computations (Lignières et al. 2006; Reese et al. 2008, 2009; Ouazzani et al. 2015) and semianalytical models (Pasek et al. 2011, 2012). These modes show regular frequency spacings and should be the most visible in the observed spectra. They have been cited to explain the frequency patterns observed in some δ Scuti stars (García Hernández et al. 2015).
Chaotic modes are called chaotic because they result from the constructive interference of waves which, in the short wavelength limit, have chaotic trajectories. They were first studied in quantum physics, showing distinctive features, such as ergodicity or the universality of the nearest neighbor statistics of their spectrum. Lignières & Georgeot (2009) identify ∼200 axisymmetric chaotic modes at a given rotation rate and verify that they follow the expected statistics.
In this paper the asymptotic properties of chaotic acoustic modes in rapidly rotating stars is investigated in detail. We construct a large set of highfrequency chaotic modes, computed at various rotation rates and azimuthal numbers m, analyze the mode properties, and interpret them using semiclassical methods.
Among these properties, the presence of peaks in the autocorrelation of the chaotic mode spectrum had not been reported in the experimental or modeled wave systems considered in the fields of quantum chaos. This justified a publication dedicated to this particular point in a physics journal (Evano et al. 2019). The present paper complements and adds to this publication, putting emphasis on the asteroseismic applications.
The paper is organized as follows: we first introduce the formalism and numerical methods, then present the properties of high frequency chaotic modes (Sect. 3), and develop the semiclassical theory which explains the features observed (Sect. 4). The appendix presents detailed derivations needed for Sect. 4.
2. Formalism and numerical methods
In this section, we first introduce the equations and tools used to study propagating pressure waves in stars. Then we present the range of numerically computed highfrequency pressure modes and finally the method used to classify the modes and obtain a set of highfrequency chaotic modes for different rotation rates.
2.1. Pressure waves and their ray limit in rapidly rotating stars
We consider adiabatic pressure perturbations in a selfgravitating gas. We are focusing on pmodes in the highfrequency regime, and as is usual in this case we apply the Cowling approximation, neglecting the perturbations of the gravitational potential. This approximation is known to be valid for highfrequency perturbations in non rotating stars (Aerts et al. 2010). We also neglect the Coriolis force since its influence on the pulsation frequencies is known to be weak in the highfrequency regime, see Reese et al. (2006). Finally, we can discard the gravity waves since the BruntVäisälä frequency is very small compared to the high pmode frequencies. With all these approximations taken into account, the adiabatic pressure perturbations obey a Helmholtzlike equation derived in Lignières & Georgeot (2009):
where , with the complex amplitude associated with the timeharmonic pressure perturbation , f a function of the background model, ω_{c} is a cutoff frequency of the model and c_{s} is the inhomogeneous sound velocity.
The shortwavelength approximation of Eq. (1) leads to the eikonal equation
This equation can be put in the form of a Hamiltonian system describing the propagation of acoustic rays (Ott 1993), by setting . The motion takes place in a meridional plane rotating with the ray at an angular velocity , where is a constant of motion and d is the distance to the axis of rotation. To compute the ray paths we use an alternative Hamiltonian form derived in Lignières & Georgeot (2009) as:
where is the frequencyscaled wave vector projected onto the meridional plane and W is the potential. The corresponding dynamical system is
where , is a timelike variable, s being the curvilinear coordinate along the ray. This system is then integrated using a fifth order RungeKutta method. To simplify the notation, we now refer to as simply .
The motion of acoustic rays inside the meridional plane corresponds to a dynamical system with N_{d} = 2 degrees of freedom. Because the wave frequency is a conserved quantity, the phase space trajectories evolve on a 2N_{d} − 1 surface. To reveal the properties of a dynamical system, it is customary to define a surface of lower dimension, called a Poincaré Surface of section (PSS), by fixing an additional parameter (Gutzwiller 1990; Ott 1993). There are multiple possible choices of PSS, the recommendations for a faithful representation are given in (Ott 1993). We have chosen a surface at constant radial distance z from the star surface r_{s}(θ): r_{p}(θ) = r_{s}(θ)−z, where r_{p}(θ) is the radial coordinate of the PSS. An example of such a PSS is shown in Fig. 1. Each dot on this figure corresponds to the crossing of the PSS with a ray approaching the surface. Chaotic zones correspond to trajectories filling densely an area, whereas regular trajectories are constrained to remain on a closed curve. The chaotic dynamics implies ergodicity in a certain phase space zone, and exponential separation of nearby trajectories, whose rate is measured by the Lyapunov exponents (Ott 1993).
Fig. 1. PSS at rotation rate Ω/Ω_{k} = 0.589 with , where θ is the colatitude and k_{θ} the projection of the wave vector on the line tangent to the r_{p}(θ) = r_{s}(θ)−z curve. Two stable islands are embedded in a chaotic region, which is itself surrounded by whispering gallery rays. The different types of rays are color coded, blue: whispering gallery, green: 6period island, yellow: 2period island and red: ergodic. 
The system undergoes a KAMlike transition (named from a famous theorem from Kolmogorov, Arnold and Moser describing the typical transition to chaos, see for example Ott 1993), with the rotation rate Ω playing the role of the perturbation. It is known that in this case the phase space contains chaotic zones, island orbits around stable periodic orbits and KAM tori reminiscent of the Ω = 0 system. This is represented in Fig. 1. We did not draw the small island structures in the domain where k_{θ}/ω is high. Ergodic trajectories appear as soon as the spherical symmetry is broken. At low rotation rates the chaotic zone is very small. However, it expands considerably as the star is flattened by rotation. The rotation rate of the model of star will be given in terms of the Keplerian breakup rotation rate , where G is Newton’s constant, M is the mass of the star and R_{eq} the equatorial radius. Near Ω/Ω_{k} ≃ 0.32 the chaotic zone is comparable in size to the main stable island and continues to grow afterward (see Fig. 2). The ratio 𝒱_{I}/𝒱_{II}, where 𝒱_{I} is the volume occupied by the main stable island and 𝒱_{II} is the volume of the main chaotic zone, is a relevant quantity as it indicates (as we shall see) the proportion of island modes to chaotic modes. Measuring the areas 𝒜_{I} and 𝒜_{II} on the PSS, corresponding respectively to the main island and chaotic zones, we found out that the function 𝒜_{I}(Ω)/𝒜_{II}(Ω) is not monotonic. In fact, the ratio of the 2period island over the chaotic zone shrinks beyond Ω/Ω_{k} = 0.658, reaches a minimum around Ω/Ω_{k} = 0.682, then increases again from this point onward. This nonmonotonic behavior was not detected in Lignières & Georgeot (2009) as rotation rates around Ω/Ω_{k} = 0.682 were not considered.
Fig. 2. Chaotic zone of the PSS at increasing values of the rotation rate. The chaotic zone grows monotonically, which is not the case for the main island zone. 
The non axisymmetric ray dynamics , shown in Fig. 3, is qualitatively similar, but the expansion of the chaotic zone is retarded. Moreover, the domain of propagation is reduced. Indeed, the ray has azimuthal velocity and, as a result, cannot come arbitrarily close to the rotation axis. In the meridional plane, this phenomenon translates into a condition on the colatitude (Lignières & Georgeot 2009) , being the frequencyscaled angular momentum norm.
Fig. 3. PSS at Ω/Ω_{k} = 0.589 with , where θ is the colatitude and k_{θ} the projection of the wave vector on the line tangent to the r_{p}(θ) = r_{s}(θ)−z curve. The phase space structures are similar to those presented in Fig. 1. 
The semiclassical theory of such a mixed phase space with both regular and chaotic zones, due to Berry and Robnik (Berry & Robnik 1984), indicates that regular and chaotic modes can be associated to the different phase space areas and form independent subspectra. This prediction was verified in Lignières & Georgeot (2009) by computing highfrequency modes and classifying them according to the phase space structure.
2.2. The set of numerically computed highfrequency pmodes
As for the ray dynamics, the star is modeled by a uniformly rotating, selfgravitating monoatomic perfect gaz of adiabatic exponent Γ = 5/3. We impose that the pressure and density satisfy a polytropic relation P ∝ ρ^{1 + 1/μ}, with μ = 3 (Hansen et al. 2004). As we are interested in highfrequency, and thus smallwavelength, pmodes, attention has to be paid to numerical resolution. The stellar model is calculated using spectral methods with Chebyshev polynomials in a pseudo radial direction and Legendre polynomial in latitude, corresponding to N_{j} = 96, 128, 140 radial points and N_{t} ≤ 185 latitudinal points. To compute its oscillation modes, we use the code TOP described in (Reese et al. 2006, 2009). Modes are computed using N_{j} Chebyshev polynomials T_{j} and N_{l} spherical harmonics , through the decomposition:
where the degrees l are either odd or even integers (Reese et al. 2006). The needed resolution is determined by the smallest scale on which the mode amplitude varies. Physically, it is directly related to the highest values taken by the components of the wave vector k_{θ} and k_{r}. Using the PSS of a given stellar model, we can determine the upper limit of the wave vector components or associated with a particular phase space structure. The chaotic zone growth seen in Fig. 2 indicates that the required angular resolution to compute chaotic modes increases rapidly with rotation. Therefore, computing highfrequency chaotic modes is more and more demanding as the model approaches Ω_{k}. The central chaotic zone is bounded by whispering gallery rays in the full phase space. Thus we know that, in a given frequency domain, the required resolution to compute all chaotic modes is attained if we are able to produce a few well resolved whispering gallery modes. With these constraints, we produced spectra at six different rotations: Ω/Ω_{k} = 0.481, 0.545, 0.589, 0.658, 0.706 and 0.809 with frequencies ranging from ∼23 ω_{p} to ∼47 ω_{p}, with , where R_{p} is the polar radius. A typical example of each kind of mode is shown in Fig. 4: a whispering gallery mode corresponding to a KAM torus, a 6period island mode corresponding to the 6period island chain, a 2period island mode corresponding to the 2period island chain and a chaotic mode corresponding to the chaotic zone. To construct a large set of chaotic modes requires to identify them among all the computed modes.
Fig. 4. Four odd axisymmetric modes at rotation Ω/Ω_{k} = 0.589: (a) chaotic mode, (b) whispering gallery mode, (c) 2period island mode (ℓ = 0), with a black line indicating the central periodic orbit, and (d) 6period island mode. The figure shows the scaled pressure amplitude , with d the distance to the rotation axis and ρ_{0} the equilibrium density. 
2.3. Mode identification
In this section we explain our methodology to isolate chaotic modes. To achieve this goal we proceed by elimination. The basic idea is to identify all regular modes and remove them from the dataset until only chaotic modes remain.
To begin with, we identify 2period island modes using the fact that their frequency spectrum is regular and of the form ω_{n, l, m} = nδ_{n}(m)+ℓδ_{ℓ}(m)+α, where n (resp. ℓ) is the number of nodes along (resp. transverse to) the central orbit γ of the island, shown in Fig. 4, and α is a constant (Pasek et al. 2012). Thus the frequency spectrum is completely determined by the regular spacings δ_{n} and δ_{ℓ}, with the half large separation ^{1}. However, the regularity of the 2period island modes may be altered by two aspects: first the deviations from asymptotic theory, expected at finite frequency, and secondly the occurrence of avoided crossings between an island mode and a chaotic mode. In practice, the island mode spacing is rigid enough to identify them. In the case of an avoided crossing, we choose arbitrarily one of the two modes as being an island mode and discard the other one. The impact of this choice on the spectrum is weak since the frequencies of two modes in the process of an avoided crossing are very close. In the same way, 6period island modes are identified on the basis of their regularity (Lignières & Georgeot 2009).
Then we remove the whispering gallery modes automatically. The PSS shows that the of any whispering gallery ray is higher than the of the chaotic zone. This means that whispering gallery modes vary on smaller latitudinal scales than chaotic modes, meaning that they have more angular nodes. Therefore, we expect their spherical harmonic expansion, Eq. (6), to be dominated by the high degree components.
For a given m, we sum over j to get the averaged coefficients . Thus the dominant degree l_{max} corresponds to the value of l where is the greatest. In the bottom panel of Fig. 5, the distribution of whispering gallery modes with respect to l_{max} is shown for a given rotation. There is indeed no whispering gallery mode below l_{max} = 10. The top panel displays the same distribution for chaotic modes, showing that they are instead dominated by their low degree components. In fact, there is a critical value l_{c} of l_{max} such that l_{max} ≤ l_{c} for every chaotic mode and l_{max} > l_{c} for any whispering gallery mode. Therefore, it is sufficient to find l_{c} to remove all whispering gallery modes from the dataset.
Fig. 5. Histograms showing the distribution of three types of odd modes as a function of their dominant degree l_{max}, at Ω/Ω_{k} = 0.589 in the frequency domain 25.6 ω_{p} to 33.5 ω_{p}. The existence of a gap between the highest value of l_{max} for chaotic modes, marked by a dashed line, and the lowest value of l_{max} for whispering gallery modes is used to automatically separate these two types of modes. 
At finite frequency, here ω ≤ 47ω_{p}, the BerryRobnik regime is not perfectly satisfied (Vidmar et al. 2007). Therefore we expect that some modes cannot be classified in a unique category. We find such modes and call them interface modes. There are two types of interface modes: whispering gallerylike and islandlike. An example of islandlike interface mode is presented in the central panel of Fig. 6. By nature, interface modes have properties of both chaotic modes and either island or whispering gallery modes. To distinguish between “well defined” modes and interface modes, we can look either at their amplitude pattern or at their Husimi distribution, which is a projection of the eigenmode in phase space (Chang & Shi 1986; Lignières & Georgeot 2009). The distinction between chaotic modes and interface modes is not clearcut, whether we use the amplitude of the modes or the Husimi distributions as a guide. Any time we encountered an ambivalent mode, we applied the following rules to categorize it. If the amplitude is high near an internal caustic, then the mode is a whispering gallerylike interface mode. The difference between chaotic, islandlike and island modes is illustrated in Fig. 6. We consider that a mode is not chaotic if its amplitude is concentrated around the central orbit of the main island. In parallel, we examine the placement of the mode in the spectrum, and look for regularities. Even at the highest frequencies considered, ℓ = 4 island modes cannot be distinguished from interface modes without ambiguity from the amplitudes.
Fig. 6. Three modes at Ω/Ω_{k} = 0.589 belonging to series 1. From left to right: we see a chaotic mode, an interface mode and an island mode with ℓ = 4. The thick white line is the main island central orbit. The modes intensity Ψ^{2} is represented, where Ψ is the scaled pressure amplitude as in Fig. 4. 
Our dataset of chaotic frequencies is described in Table 1. All axisymmetric modes were identified using the methodology described above. For non axisymmetric modes the last step, i.e. the identification of interface modes through the amplitude of the mode or the Husimi distribution, was not performed. The modes listed in Table 1 are separated into symmetry classes associated with the quantum number m and the parity with respect to the equator (m^{−} odd, m^{+} even). A distinctive property of chaotic mode spectra is to have a universal statistical behavior (Gutzwiller 1990). In particular, their nearest neighbor distribution , with , where ⟨ω_{n + 1} − ω_{n}⟩ corresponds to the mean frequency spacing, is expected to follow the Wigner surmise of Random Matrix Theory (Gutzwiller 1990). This is checked in Fig. 7 and it confirms that our method for selecting chaotic modes works properly. We aggregated the eight independent spectra to obtain better statistics (left panel). Using this standard aggregating procedure enables us to obtain much better statistics than in previous studies (Lignières & Georgeot 2008). The integrated distribution , displayed in the right panel, shows that the agreement is good for individual spectra as well. We did not take into account the specific case of Ω/Ω_{k} = 0.706 which corresponds to an anomalous statistics, intermediate between the Wigner and Poisson distributions, with no level repulsion (Evano et al. 2019). Such a distribution is characteristic of a superposition of independent spectra, and will be explained later (see Sect. 4.2).
Characteristics of chaotic frequency spectra.
Fig. 7. Left: nearest neighbors spacing distribution , with 1344 frequency levels obtained from eight independent spectra: Ω/Ω_{k} = 0.481 (206 odd levels), Ω/Ω_{k} = 0.545 (223 odd levels, 105 even levels), Ω/Ω_{k} = 0.589 (217 odd levels, 96 even levels), Ω/Ω_{k} = 0.658 (207 odd levels, 120 even levels) and Ω/Ω_{k} = 0.809 (170 odd levels). Right: integrated distribution N(s) for all eight independent spectra. In both panels the dashed line the Wigner surmise and the dotted line is the prediction for Poissonian spectra. 
3. Properties of highfrequency chaotic modes
3.1. Spatial distribution of the mode amplitudes
We start with a qualitative description of the main spatial features of chaotic modes. Representative examples of their amplitude distribution in a meridional plane are displayed in Fig. 8. A first observation is that their nodal pattern is complex in the sense that, contrary to regular modes, there is no simple way to count the number of nodes. Nevertheless we also notice that while the inner part of the chaotic modes looks random, the outer part is much more structured. Indeed it is possible, near the surface, to count radial and angular nodes. The nodes appears regularly spaced in the radial direction but are unevenly distributed in the angular direction. The amplitude distribution of chaotic modes is therefore mainly characterized by its irregular nature, but with radial regularities near the surface. In the next subsection, we will show that the chaotic spectra present some regularities as well. Another notable property of the chaotic modes presented in Fig. 8 is that they spread out in the entire stellar interior. This distinguishes chaotic modes from regular modes that are confined in a narrow part of the star: whispering gallery modes stay close to the star surface and island modes are trapped in the vicinity of a central periodic orbit (see Fig. 4). Chaotic modes are thus the only class of pmodes able to probe the star center at high frequency.
Fig. 8. Scaled pressure amplitude Ψ shown for twelve odd axisymmetric chaotic modes at Ω/Ω_{k} = 0.589 with quantum number m = 0. 
3.2. Regularities in the spectra
We computed the autocorrelation R_{2}(ξ) = ∫d(ω − 1/2ξ)d(ω + 1/2ξ) dω of a variety of chaotic spectra, where the density d(ω) was obtained by convoluting the spectra with a Gaussian function of small standard deviation compared to the mean frequency spacing and of height unity. The results are shown in Fig. 9 at five different rotations for the case of axisymmetric modes with odd equatorial parity. The figure clearly shows peaks emerging from the noise level, which are not predicted by Random Matrix Theory nor seen in generic chaotic spectra. The peak that appears at every rotation, and which is usually the most visible, is referred to as the “main peak”. Its position, denoted Δ_{c}, slightly decreases with increasing rotation. Moreover, additional peaks of significant amplitude appear in the autocorrelations, their relative amplitude being large especially at Ω/Ω_{k} = 0.706. In the following we characterize the spectrum structure behind these peaks (the next two subsections) and then extend our analysis to the nonaxisymmetric modes.
Fig. 9. Autocorrelations R_{2}(ξ), where ξ is a displacement in frequency, of chaotic spectra with odd parity: (a) 206 levels from 28.35 ω_{p} to 46.89 ω_{p}, (b) 223 levels from 28.15 ω_{p} to 44.09 ω_{p}, (c) 217 levels from 26.02 ω_{p} to 40.29 ω_{p}, (d) 283 levels from 23.57 ω_{p} to 36.22 ω_{p} and (e) 170 levels from 24.02 ω_{p} to 30.01 ω_{p}. The autocorrelations have been normalized such that their value at the origin is unity. The main peak’s position is labeled Δ_{c} and marked with a solid line. Secondary peaks are marked in panels c, d and e with a dashed line. 
3.2.1. Main peak
The systematic presence of the main peak in the autocorrelations is a hint that the spectra are structured by the characteristic frequency spacing Δ_{c}. To verify this idea, we make use of socalled échelle diagrams. On such diagrams, two modes distant in frequency of Δ will be represented by two points sharing the same abscissa and separated by one unit on the y axis. We consider a portion of spectrum at Ω/Ω_{k} = 0.589, and analyze its structure in a detailed way. In Fig. 10, the frequencies of all odd parity axisymmetric chaotic modes between 25.07 ω_{p} and 33.78 ω_{p} are represented on an échelle diagram, using a folding value equal to the main peak Δ_{c} = 1.0899 ω_{p} of the autocorrelation. The frequencies are grouped in approximate vertical lines, some of which are very well aligned, whereas others seem to form diagonal or wavy tracks (Changing slightly the folding value breaks the best alignments but other tracks line up instead). The échelle diagram shows that the chaotic spectrum can be split into series of modes separated in frequency by approximately Δ_{c}.
Fig. 10. Échelle diagram of chaotic modes at Ω/Ω_{k} = 0.589 in the range 25.60 ω_{p} to 33.54 ω_{p}, with odd parity. 
This property is wellknown for lowdegree modes in nonrotating stars as well as island modes in rotating ones (see Fig. 11). The frequency spacing between two consecutive modes is the socalled large separation and their amplitude distributions only differs by the number of nodes along a particular direction (radial for modes in a non rotating star and along the central periodic orbit for island modes). In such regular spectra, this structuring of the frequencies is directly related to the classical dynamics through socalled EinsteinBrillouinKeller theory (Pasek et al. 2011, 2012). Its appearance in a chaotic spectrum is more mysterious, and will be explained in Sect. 4.
Fig. 11. Échelle diagram of island modes ℓ = 0, 1, 2, 3, 4 at Ω/Ω_{k} = 0.589, with only odd modes. The lowest point of the ℓ = 4 track corresponds to an islandlike interface mode. 
In the autocorrelations of even chaotic mode spectra, both the main peak and secondary peaks usually have a smaller amplitude. Nonetheless, even mode frequencies produce ridges in the échelle diagram, just like odd mode frequencies. In Lignières et al. (2006), it was pointed out that even modes are more strongly impacted by avoided crossings. This fact could explain why even modes are less regular than odd modes.
Comparing the amplitude patterns of the chaotic modes that belong to the same track on the échelle diagram, we find out that consecutive modes are often very similar. This is illustrated by two examples of consecutive modes in Fig. 12. From a systematical search for couples of modes with a similar amplitude distribution that are separated in frequency by approximately Δ_{c} ^{2}, we were able to label the vertical series of chaotic mode échelle diagram, endingup with the sixteen series shown in Fig. 13. The very few modes left without a label are modes undergoing an avoided crossing.
Fig. 12. Mode intensity Ψ^{2} at Ω/Ω_{k} = 0.589, where Ψ is the scaled pressure amplitude, showing the similarity between consecutive modes. Top: two consecutive modes that belong to series 8. Bottom: two consecutive modes that belong to series 3. 
Fig. 13. Échelle diagram showing all chaotic frequencies in the range 25.60 ω_{p} to 33.54 ω_{p} modulo the large separation Δ_{c}. Series of modes are labeled by numbers 1 to 16. Black dots correspond to modes that do not fit into a series. 
This analysis indicates that the spectrum organization in series of modes separated by a fixed frequency spacing and showing similarities in their spatial distribution is also relevant for chaotic modes. In this context, Δ_{c} can be interpreted as a large separation for chaotic modes. Its value turns out to be close to the island mode large separation in the same frequency range Δ_{i} = 1.0907 ω_{p}.
There are nevertheless some important differences with nonrotating and island modes. First the frequency spacing is much more regular for nonrotating or island modes than for chaotic modes. This is obvious from the comparison of the chaotic mode échelle diagram with the échelle diagram of the 2period island modes shown in Fig. 11. The similarity of the amplitude distributions along a vertical series is also much stronger for island modes than for chaotic modes. In agreement with the asymptotic theory of island modes (Pasek et al. 2012), island modes of the same series have the same number of nodes, denoted ℓ, in the direction perpendicular to the periodic orbit. By contrast, the comparison between two chaotic modes of the same series but separated by a few Δ_{c} is not as clear, since the patterns slowly evolves from one mode to the next. Another important difference comes from the characteristics of the series. As expected from the theory and as observed in Fig. 11, island mode series should not stop toward high frequencies as modes with the same ℓ and higher n remain aligned in the échelle diagram. This is not the case for chaotic modes as some series like series number 1, 2, and 13 on Fig. 13 come to an end in the frequency range considered. In parallel, some series starts above a given frequency, for example the series number 16 appearing above ω = 27Δ_{c}. Typically, interface modes appear at the start of a series as whispering gallerylike modes or at the end of a series as islandlike modes. In our data, the end of a series of chaotic modes is also the start of a series of island modes of given ℓ; for instance, series 1 is followed by ℓ = 4 island modes and series 2 by ℓ = 5 island modes. This transition is visible by comparing Fig. 13, where series 1 ends just below 27Δ_{c}, and Fig. 11 where the ℓ = 4 island modes start just above 27Δ_{c}.
3.2.2. Secondary peaks
In addition to the main peak discussed above, we see other peaks that we call secondary peaks. At Ω/Ω_{k} = 0.589, Ω/Ω_{k} = 0.706 and Ω/Ω_{k} = 0.809, one can clearly distinguish these peaks from the background noise. They are marked by a dashed line in Fig. 9. The peaks occur at αδ + βΔ_{c}, where δ is the position of the first secondary peak and α and β are integers. For instance, the peaks of Fig. 9, panel c, occur at (α = +1, β = 0) and (α = +1, β = +1). At Ω/Ω_{k} = 0.706, secondary peaks are numerous and close in amplitude to the main peak. At this rotation, we have already noticed the nearest neighbors statistics is anomalous (Evano et al. 2019), intermediate between the Wigner and Poisson distributions with no level repulsion. This anomalous statistics is a sign that independent subspectra coexist at Ω/Ω_{k} = 0.706. We will propose an explanation for these secondary peaks in Sect. 4.
3.3. Nonaxisymmetric modes
Because the model is cylindrically symmetric, modes are quantized in the azimuthal direction, the quantum number being denoted m. In the ray dynamics, it corresponds to the quantization of the invariant associated with this symmetry, that is the projection L_{z} of the angular momentum on the axis of rotation. For a given rotation and frequency, increasing has the effect of reducing the size of the chaotic zone. This can be seen by comparing the PSS of Fig. 1 with m = 0 and the PSS of Fig. 3 with m = 4, both with ω = 24.41 ω_{p}. From the mode numerical computations, we find for m = 1 the main peak is still clearly visible in the autocorrelation, and almost exactly at the same position as for the axisymmetric case. However for m = 4 the autocorrelation shows a forest of peaks, the main peak being slightly shifted toward low values (see Fig. 14).
Fig. 14. Top panel: autocorrelation at Ω/ω_{k} = 0.589 for quantum number m = 1, in the frequency domain 30.51 ω_{p} to 38.48 ω_{p}. Bottom panel: autocorrelation at Ω/ω_{k} = 0.589 for quantum number m = 4, in the frequency domain 30.53 ω_{p} to 38.51 ω_{p}. In both cases, the dashed line is the position Δ_{c} of the main peak for axisymmetric modes at the same rotation rate. 
4. Semiclassical interpretation
In Sect. 3, we gave a description of the most important properties of chaotic modes in our model of rotating stars. Our goal in the present section is to provide a theoretical understanding of these results based on asymptotic methods. The statistical properties of chaotic spectra are well described by Random Matrix Theory (Bohigas 1991; Mehta 2004), which studies the distribution of the eigenvalues of matrices filled with Gaussianly distributed random numbers. Real symmetric random matrices are of particular interest. Indeed, they form the Gaussian Orthogonal Ensemble (GOE) that models the spectra of timereversal symmetric systems. Comparing an autocorrelation of the stellar model with the GOE autocorrelation (see Fig. 15) makes it clear that the main peak and secondary peaks are not generic features of wave chaos. If Random Matrix Theory is efficient at modeling the generic properties of chaotic systems, it is not able to grasp specific behaviors that may arise in particular systems (see e.g., Bogomolny et al. 1992, 1997). However, semiclassical methods based on the propagation of rays are well suited for this task.
Fig. 15. Comparison of the “stellar” autocorrelation R_{2}(ξ) at Ω/Ω_{k} = 0.589, represented in solid line, with the GOE autocorrelation represented by a dashed line. To compare the “stellar” and GOE autocorrelations, two changes have been made. First, the autocorrelation of our model has been rescaled so that the value of the plateau, i.e. the line around which there are fluctuations, is approximately unity. The value of the plateau otherwise depends on the mean density and the Gaussian smoothing of the spectrum. Secondly, the GOE distribution has been rescaled, in the direction of the x axis, by the mean level spacing of the Ω/Ω_{k} = 0.589 spectrum in the frequency domain considered. 
Using the semiclassical periodic orbit theory, which relates the mode properties to the acoustic ray dynamics, we show in Sect. 4.1 that the Δ_{c} regularity is caused by the strong decrease of the sound speed at the surface, and present a theory which predicts the occurrence of the peak and its characteristics from the ray dynamics. We then explain the presence of secondary peaks by the transport properties of the phase space in Sect. 4.2. In Sect. 4.3 we show, using a simplified model, how the behavior of acoustic rays near the surface may induce structure in the nodal pattern. In Sect. 4.4 we discuss how the families of chaotic modes evolve when the frequency domain changes. At last in Sect. 4.5, we show that the proximity of Δ_{c} and Δ_{i} is not accidental, and we propose a way to differentiate chaotic modes and island modes using the symmetries of the system.
4.1. Main autocorrelation peak
4.1.1. Periodic orbit theory
Before diving into the subject of periodic orbit theory, let us introduce a few quantities that will be later used to characterize the statistical properties of the spectra. The staircase function N(ω) counts the number of modes below a given frequency ω. From the staircase function one can define the spectral density d(ω)=dN/dω. Both the staircase function and the spectral density are often written as the sum of two contributions called the smooth part (or ) and the fluctuating (or oscillating) part N^{osc}(ω) (or d^{osc}(ω)), where the smooth part is obtained by locally averaging the function in the neighborhood of a target frequency.
Periodic orbit theory is an asymptotic semiclassical theory developed in the limit of high frequency or short wavelength. It is based on the trace formula (Gutzwiller 1990), which relates the spectral density to a sum over the periodic trajectories of the Hamiltonian system. Though it was originally derived in the framework of quantum mechanics, the trace formula can be adapted to any wave system with a ray approximation. Indeed trace formulas have been built and tested, for instance in optics using microwave cavities (Kudrolli et al. 1994) or for plate vibrations (Bogomolny & Hugues 1998). In the same spirit, we rederive the trace formula for the system considered here in Appendix A. Using j as a label for the periodic orbits γ_{j}, the formula reads:
where is the action and is the acoustic travel time of γ_{j}. The amplitude term is expressed as
where I denotes the identity matrix and M_{j} is the monodromy matrix which describes the linearized motion around the periodic orbit and whose eigenvalues give the stability of the orbit. The periodic orbits are either socalled primitive orbits or repetitions of them. In many systems, the density ρ(T) of orbits with a time period T grows exponentially as
where λ is the average lyapunov exponent of the system, which describes the rate at which nearby trajectories diverge. Long orbits are less stable and their amplitude drops down as
The exponential growth of the number of long orbits makes it numerically difficult to attain good precision on the spectral density using the trace formula. Moreover in many systems, including the stellar model studied here, there is no automatic procedure to find the periodic orbits. To circumvent these issues, Berry showed (Berry 1985) that some statistical quantities such as the variance, or more importantly the autocorrelation of the spectrum, can be approximated using the trace formula in systems where individual periodic orbits are not known, provided one knows their distribution.
The autocorrelation R_{2}(ξ) = ⟨d(ω−1/2 ξ) d(ω+1/2 ξ)⟩, with the average ⟨f⟩ = ∫fdω, can be rewritten as , where ϕ(ξ) is a smooth function. The theory aims at deriving an approximate expression of the form factor:
which is the Fourier transform of the autocorrelation
Getting rid of the real parts through the fact that with and then , gives
Inserting this expression in Eq. (11) and expanding the action as S_{j}(ω ± 1/2ξ)≈S_{j}(ω)±1/2 ξ (∂S_{j}/∂ω) = S_{j}(ω)±1/2 ξT_{j}(ω), where T_{j} is the acoustic travel time introduced before, leads to a new expression of the form factor:
Under the frequency average ⟨.⟩, the contribution of the offdiagonal terms i ≠ j in the double sum can be neglected. This diagonal approximation is valid for “short times” below the Ehrenfest time T_{E} ≈ (1/λ)ln(ω) (Bogomolny et al. 1997). Beyond T_{E}, there are pairs of orbits with very close actions S_{i} ≈ S_{j} and higher order terms need to be computed (Bogomolny & Keating 1996; Sieber & Richter 2001). In the regime where the diagonal approximation is valid, the expression of the form factor is reduced to
In this expression K(T), as a Fourier transform of a correlation function in frequency, is a function of time T. Equation (16) means that K(T) is related to the distribution of travel times of all the periodic orbits labeled by j (and not of their lengths). This will be crucial in the subsequent analysis. For long orbits Eq. (16) can be written as
where ρ(T) dT is the number of periodic orbit with a time period between T and T + dT. A prescription for the amplitude A(T) and the density ρ(T) is given in Eqs. (10) and (9), and can be generalized (Hannay & Ozorio De Almeida 1984), leading to the conclusion that for generic chaotic systems (Berry 1985) the form factor is a linear function of T i.e. K(T)∝T (up to the Ehrenfest time), which is in accordance with the predictions of GOE (Bohigas 1991). Such a result does not predict the occurrence of peaks in the autocorrelation. We will show that the same theory but using the specific density of periodic orbits of our system leads to a different result predicting autocorrelation peaks.
4.1.2. Distribution of acoustic travel times
To this aim, we need to model the density ρ(T) of periodic orbit travel times in our system. First, we define a chord as a portion of orbit between two consecutive rebounds at the surface. From this definition, a trajectory that bounces n times at the surface will be called a nchord. For any integer n, periodic trajectories with n rebounds are a subset of nchord trajectories. As already mentioned, we do not know any systematic way of finding all the periodic orbits in the system. Nonetheless, we may infer some of their properties by studying large samples of nchord trajectories.
At a given rotation rate, we compute the acoustic travel times of thousands of nchord trajectories. The 1chord distribution is shown in the right panel of Fig. 16, at Ω/Ω_{k} = 0.481. The distribution is a narrow packet: its standard deviation σ_{0} is small compared to its mean value T_{0}. We interpret this effect as a consequence of the strong decrease of the sound speed near the surface^{3}. Indeed, the trajectories travel rapidly through the core of the star and the acoustic time is dominated by the surface behavior, so the actual length of the trajectory has very little impact on its travel time. We expect the characteristic time T_{0} to increase with rotation, because of the increasing volume of the star. This is indeed confirmed, with T_{0} = 5.19/ω_{p} at Ω/Ω_{k} = 0.481 and T_{0} = 6.30/ω_{p} at Ω/Ω_{k} = 0.706.
Fig. 16. Left panels: number of nchord trajectories, with n = 1, …, 20, vs their travel time T at rotations Ω/Ω_{k} = 0.481 (top) and Ω/Ω_{k} = 0.809 (bottom), with 300 bins for the total distribution. The nchord samples contain ∼4200 chords each. Right panel: Ω/Ω_{k} = 0.481 1chord distribution in more details, with ∼84 000 chords and 100 bins. The mean value T_{0} of the distribution is marked with a dashed line and the standard deviation σ_{0} is shown. 
Moreover, due to the centrifugal deformation, the region where the sound speed is very small widens at the equator. Hence the dispersion of travel times must increase as well with rotation, as can be seen by comparing, for instance, the chord joining the two poles with the chord joining opposite points on the equator. The total distribution of nchords with n ≤ 20, with the same number of orbits for all values of n, is represented in Fig. 16, panels a and b at rotations Ω/Ω_{k} = 0.481 and Ω/Ω_{k} = 0.809. The packets are evenly spaced out (T_{n} ≈ nT_{0}) but get thicker as n increases (), since each packet can be seen as the sum of n independent variables. Thus this packet structure will disappear when n becomes large, i.e. for long travel times. An important difference between the two distributions shown in Fig. 16 is the rate at which adjacent packets overlap, leading to the disappearance of the packet structure. This is quantified by the ratio σ_{0}/T_{0} which is 0.049 at Ω/Ω_{k} = 0.481 and 0.097 at Ω/Ω_{k} = 0.809.
The constraints imposed by the nchord acoustic time distribution could be strong enough to impose a kind of periodic oscillation in the distribution of periodic orbits, and thus in the form factor K(T). This oscillation could then produce a peak in the autocorrelation in virtue of the Fourier relation^{4}. We will now show that it is indeed the case
To this aim, we model the distribution of travel times as a sum of Gaussian functions: P_{Ω}(T) = ∑_{n}P_{n, Ω}(T), with P_{n, Ω}(T) the probability distribution of nchords travel times:
where T_{0}, in the prefactor, has been added for normalization purposes. It is a crude approximation but it encapsulates in a simple form all the relevant properties of the distributions, namely the presence of regularly spaced packets and their widening with increasing n. One has to keep in mind that the dependency on Ω is not explicit but hidden in the values of T_{0} and σ_{0}. Finally, we need to take into account the fact that the number of trajectories of acoustic time T grows with T. Thus, we write ρ(T) = (1/T)e^{λT} × P_{Ω}(T) and the form factor becomes K(T)∝A^{2}(T) ρ(T) = T P_{Ω}(T).
The Fourier transform of this quantity
is shown in the middle panel of Fig. 17 for six rotations corresponding to the simulated chaotic spectra. This quantity is a semiclassical approximation of C(ξ) (see Eq. (12)) and is closely related to the autocorrelation R_{2}(ξ) (see Fig. 9). The results show that the periodic orbit theory based on the ray model indeed predicts a peak in the autocorrelation. In the middle panel of Fig. 17, its theoretical position moves with rotation, as in the numerical mode computations, from 1.1458 ω_{p} at Ω/Ω_{k} = 0.481 to 0.9254 ω_{p} at Ω/Ω_{k} = 0.809. In the top panel, is compared to the peak position of the numerical modes Δ_{c}, with good agreement. The position, height and width of the peaks are completely determined by T_{0} and σ_{0}. The position is found straightforwardly as . The height and width are controlled by the ratio σ_{0}/T_{0}. Indeed, the condition σ_{0} ≪ T_{0} is necessary to clearly distinguish the packets in the distribution of travel times, as in Fig. 16. Due to the increase of σ_{0}/T_{0} with rotation, the peak gets less visible at high rotation rates.
Fig. 17. Upper panel: theoretical large separation (upward triangles) compared to the numerical peak’s position (downward triangles) calculated for axisymmetric modes at six rotation rates. Middle panel: theoretical autocorrelations with quantum number m = 0, from right to left the rotation rate increases. Bottom panel: theoretical autocorrelations at Ω/Ω_{k} = 0.589 with m = 0, 2, 4, 6 and frequency ω = 24.41 ω_{p}. 
The theoretical large separation is not only a function of the rotation rate Ω, but also of the projected angular momentum . In the bottom panel of Fig. 17 we show the predicted peaks for and 0.2458/ω_{p} at Ω/Ω_{k} = 0.589. The shift in position is governed by the value of T_{0} and the change in amplitude by the value of σ_{0} as increases. The main difference with the axisymmetric case is a reduction of the domain of propagation, resulting in a smaller standard deviation σ_{0}. As the mean travel time T_{0} between two rebounds depends on the impact of on the ray paths, it changes slightly in a non monotonic fashion.
We have seen that the sound speed profile in the star imposes restrictions on the travel times of acoustic rays. In spite of the chaotic nature of the dynamics, the small dispersion of travel times produces order in the spectra, in the form of a peak in the autocorrelation at a value that corresponds to the mean travel time between two points at the surface. It is fundamentally a radial phenomenon, since the strong variation of the sound speed occurs in the radial direction. Introducing the radial acoustic time at a given colatitude , where r is the radial coordinate, the mean acoustic time between two rebounds can be estimated without using ray dynamics, as
where the integral is performed over a quarter of the star because of the axial and equatorial symmetries. In the m ≠ 0 case, the domain of integration must be adapted, since the acoustic ray cavity is reduced in size. We find that 2π/T^{av} is a good estimation of Δ_{c}. For instance, at ω/Ω_{k} = 0.589, 2π/T^{av} = 1.0996 ω_{p}, which differs from Δ_{c} by less than 1%.
4.2. Secondary peaks
At every rotation there are other peaks in addition to the main peak caused by the large separation of chaotic modes. At rotations higher than Ω/Ω_{k} = 0.589 some of these peaks raise high above the noise level and we call them secondary peaks. The presence of small amplitude peaks in the autocorrelation can be understood as reflecting the organization of the chaotic spectra on an échelle diagram, illustrated in Fig. 10. Indeed, the spectrum is organized in families of nearly aligned frequencies. Lets us consider first the case of two perfectly aligned families of frequencies. These two families would produce a peak in the autocorrelation at a position given by the distance between consecutive levels in the two series. This is reminiscent of the case of island modes, that show wellaligned families on the échelle diagram and, accordingly, many peaks are seen in the autocorrelation. In this later case, the effect is strong because the separation δ_{ℓ} = ω_{n, ℓ + 1} − ω_{n, ℓ} between consecutive families is fixed. In the case of chaotic modes, the alignments are weaker and consecutive series are not regularly spaced out. Thus, we expect the peaks to be of small amplitude.
The high amplitude peaks observed in the data at Ω/Ω_{k} = 0.589 are caused in part by two series, number 1 and 2 in Fig. 10, that are nearly parallel on the échelle diagram. The spacing between these two series indeed corresponds to the position δ of the leftmost secondary peak in Fig. 9, panel c. If couples of nearly parallel series can occur occasionally, there is however no reason to expect them at every rotation rate. The presence of strong secondary peaks at other rotation rates (and also in the spectra of non axisymmetric modes) is a hint that a more generic mechanism is at play. Moreover, we know that the spectral statistics of the chaotic spectrum at Ω/Ω_{k} = 0.706 indicates that it is divided into independent subspectra, as explained in Sect. 3.2.2. In the following, we will argue that the presence of secondary peaks is a consequence of the presence of partial barriers in the phase space of the ray system which create separate independent subsets of modes.
It is known that transport properties in chaotic phase space can have a significant impact on the wave system spectra (Bohigas et al. 1993). The transport of trajectories from a subregion A to another subregion B can be affected by the presence of partial barriers. These are curves through which classically trajectories can flow, but with a much smaller flux than in other parts of phase space. Thus, contrary to KAM tori that act as complete barriers, ergodic trajectories are able to cross partial barriers after a sufficiently long time. Such partial barriers are typically created when the system is being perturbed, through the destabilization of island orbits. In virtue of the PoincaréBirkhoff theorem (Ott 1993), resonant tori are destroyed by the perturbation and generate new structures in phase space: a new (smaller) island chain around an elliptic central orbit along with an unstable periodic orbit. The unstable orbits created through this process are known to be the source of partial barriers that trap ergodic trajectories (Shim et al. 2011; Bohigas et al. 1993).
Acoustic ray simulations show that two 6period unstable orbit revolves around the main 2period island chain at Ω/Ω_{k} = 0.589. At Ω/Ω_{k} = 0.706, we also find a periodic orbit revolving around the island and strongly suspect the presence of a second one. Moreover, by choosing trajectories with initial conditions around the main island and evolving the system forward in time, we see clearly the contours of the partial barrier reveal themselves (see Fig. 18). For the frequencies considered in our dataset, the partial barriers may act as complete barriers and isolate some modes, as observed in other systems (Shim et al. 2011). Quantitative estimates of the area and outgoing flux have been given in Evano et al. (2019), showing that the barrier grows in size from Ω/Ω_{k} = 0.589 to Ω/Ω_{k} = 0.706 and that, in parallel, it takes longer to go through the partial barrier at Ω/Ω_{k} = 0.706. Thus, the trapping of chaotic trajectories around stable islands seems to be the cause of the additional peaks seen in the autocorrelations. Moreover, as the trapped trajectories revolve around a 6periodic orbit (see Fig. 18), and assuming the modes to quantize like island modes, one can expect that the secondary peaks will be located approximately at Δ_{c}/3. As seen in Fig. 9, this idea is consistent with the data at ω/ω_{k} = 0.706. However it gives only a rough estimate of the peak position at other rotation rates and thus this issue should be analyzed in more depth. The specific rotations where the trapping of trajectories will be efficient can be anticipated only through precise numerical simulations of the ray model.
Fig. 18. Snapshot of the evolution of a bundle of trajectories as they intersect the PSS, represented by black dots, at Ω/Ω_{k} = 0.706. The phase space zones where black dots are dense correspond to regions enclosed by partial barriers. The trajectories are initially in the neighborhood of the two main islands. One of the central periodic orbits of these zones is indicated by colored dots (large gray dots). 
4.3. Amplitude distribution of chaotic modes
Berry proposed that the eigenstates of a classically chaotic system, like quantum billiards, could be modeled by a superposition of random plane waves. Indeed this picture holds locally, in the high frequency regime, where the mode is a superposition of rays of the same magnitude k, but coming from all possibles directions. It is easy to test this idea, as done e.g. in O’Connor et al. (1987), by adding multiple timeharmonic solutions of the Helmholtz equation ∇^{2}ψ + k^{2}ψ = 0. The waves are of the form ψ_{n} = a_{n} cos(k_{n}(α_{n})⋅x + ξ_{n}), where the amplitude a_{n}, the wave vector orientation α_{n} and the phase shift ξ_{n} are random variables but the magnitude k of the wave vector is fixed. The resulting modes indeed exhibit the random ridges characteristic of chaotic modes.
Contrary to quantum billiards, the wavevector magnitude varies strongly within a star, as it is proportional to the inverse of the sound speed. Moreover, each incoming ray approaching the surface is almost aligned in the radial direction (k_{r} ≫ k_{θ}). Thus, the hypothesis that intersecting rays come with a large variety of possible orientations is not valid near the surface. We modified slightly the random wave model to incorporate such a behavior. For simplicity, we treat the x axis as a radial direction. Then we impose a radial increase of the magnitude k and that all waves end up aligned in the radial direction as x → ∞. An example of mode produced in this way is shown in Fig. 19. The mode exhibits random ridges in the center but the outer part is more structured, with nodes regularly spaced radially and irregular in the transverse direction as in the star.
Fig. 19. Superposition of random stationary waves reproducing the kind of amplitude patterns seen in the chaotic pmodes of the polytropic stellar model (see text). 
4.4. Spectrum organization
In the previous sections, we found that the chaotic spectrum can be described as a set of series of modes, where a series corresponds to modes separated by approximately Δ_{c} and having similar amplitude distribution. Sixteen series have been effectively identified at the rotation rate Ω/Ω_{k} = 0.589 and in the frequency range, 25.60 ω_{p} < ω < 33.54 ω_{p}. We also found that whispering gallery modes are present at the lowfrequency end of some chaotic series, and that island mode series appears as highfrequency extensions of some other chaotic series. In this section, we use semiclassical arguments to investigate the origin of the different chaotic series and their link with the whispering gallery and island mode spectra.
We first consider a non rotating star because in this case the modes can be precisely located on the PSS. As the system is integrable, the phase space is foliated with N_{d}dimensional tori, where N_{d} is the number of degrees of freedom. Modes are then constructed on some of the tori, the ones specified by the EinsteinBrillouinKeller quantization rules. The quantization of the norm of the angular momentum L, the invariant associated with the spherical symmetry, reads L = ±(ℓ_{s} + 1/2) (Gough 1993) where ℓ_{s} is the degree of the mode, with the index s denoting the “spherical” case. For the frequencyscaled coordinates used here, we rather use the invariant , where n_{s} and m are the radial order and the azimuthal order of the mode. In the axisymmetric case, L_{z} = 0, the tori imprint the PSS on horizontal lines (Lignières & Georgeot 2009). On Fig. 20, the identification of a few modecarrying tori is displayed. It shows how the position of the tori along the vertical axis depends on ℓ_{s} and n_{s}. For fixed radial order n_{s} and increasing ℓ_{s}, the tori indeed move toward higher or values. Similarly, for fixed ℓ_{s} and increasing n_{s}, the tori approach the axis. In particular, the ℓ_{s} = 0 or ℓ_{s} = 1 modes are already close to the central axis for the smallest radial order.
Fig. 20. Intersection of a few tori with the PSS at Ω/Ω_{k} = 0. Only the part of the PSS where k_{θ} > 0 is shown because of the symmetry with respect to the k_{θ} = 0 axis at zero rotation. The tori correspond to modes of degree ℓ_{s} = 1, ℓ_{s} = 9 or ℓ_{s} = 11. 
When rotation comes into play, the location of the tori at Ω = 0 has strong consequences on their fate. From the evolution of the PSS we know that the phase space becomes rapidly dominated by three main structures: the 2period island chain at low , the large chaotic zone at low and intermediate , and at large , the region of surviving KAM tori corresponding to whispering gallery trajectories. In this context, high tori at Ω = 0 will transform into structures of the surviving KAM tori region, whereas low tori will be destroyed as the 2period island chain and the chaotic zone surrounding it emerge. We thus expect high modes to become whispering gallery modes, while low modes should evolve toward chaotic or island modes.
We can use this phenomenology to predict the fate of a sequence of modes having a fixed degree ℓ_{s} and variable radial orders n_{s}. In such a series, has a maximum value for n_{s} = 1 and it decreases toward zero as n_{s} increases. Generically, we thus expect that, as rotation increases, the lown_{s} (high) modes become whispering gallery modes, the intermediaten_{s} modes become chaotic modes, and the highn_{s} (low) modes become island modes. Combining this picture with our observation that some series of chaotic modes show island modes at their highfrequency end or whispering gallery modes at their lowfrequency end, we are led to interpret the chaotic series together with their island modes and whispering gallery modes extensions as the remnants of series of given ℓ_{s} at zero rotation.
At Ω/Ω_{k} = 0.589, we could indeed relate two chaotic series (the series 1 and 2) with the ℓ = 4 and the ℓ = 5 oddparity island mode series, respectively. Following this interpretation and using the formulas that link ℓ to ℓ_{s} (Reese et al. 2008; Pasek et al. 2012) we can attribute the ℓ_{s} = 9 value to series 1 and the ℓ_{s} = 11 value to series 2. In principle, the ℓ_{s} value of the other chaotic series shown on Fig. 13 could also be determined, either by following them to higher frequency up to the island mode transition or by looking for a whispering gallery mode and its ℓ_{s} value at low frequency.
The generic case of a ℓ_{s} series containing the three type of modes only holds for high enough ℓ_{s}. The computed PSS indeed show that below a specific that depends on the rotation rate, all trajectories are either chaotic or within an island chain. We thus expect that, below some critical ℓ_{s} that depends on the rotation rate, all modes in the series are either of the chaotic or island type. Moreover, from numerical studies (Lignières et al. 2006; Reese et al. 2008; Pasek et al. 2012) where ℓ_{s} = {0, 1, 2, 3} modes have been carefully followed with rotation, we know that, for these lowest degree, modes of all orders behave as island modes. This is coherent with the fact that their initial are all close to the axis.
To summarize, we argued that the chaotic mode spectrum of rapidly rotating stars is organized in series that can be traced back to the fixed ℓ_{s} series of the nonrotating star. While supported by the analysis of our numerical results, this idea needs to be further tested at other rotations and in other frequency domains. It is also important to stress that without the Δ_{c} organization of the chaotic spectrum, the chaotic series would not exist and we could not attribute ℓ_{s} value to them. Indeed, tracing back chaotic modes to their integrable counterpart is not possible for typical mixed systems like quantum billiards. Another property that help recognize the link with a highfrequency island mode series is that Δ_{i} is very close to Δ_{c}.
4.5. Mode identification: chaos versus islands
In the numerically computed spectra, the large separations of island modes and chaotic modes are found to be very close to one another (see Table 2). The asymptotic theory discussed above gives a natural explanation for this apparent coincidence. Indeed we established that the large separation of chaotic modes is defined asymptotically by . On the other hand, the large separation of island modes in the asymptotic regime is related to the acoustic time along the central periodic orbit γ by (Pasek et al. 2012) Δ_{i} = 2δ_{n} with . The closeness of the two peaks is thus due to the closeness of the mean travel time of a chaotic trajectory and the travel time along the central path γ. For the same reason that all chaotic trajectories have almost the same travel time, the acoustic time between two rebounds along the central orbit of the island has to be very close to T_{0}.
Comparison of the large separation of island modes Δ_{i} and chaotic modes Δ_{c}, obtained from the simulated spectra at six rotations.
Despite the proximity of the two peaks, it may be possible to tell them apart by combining odd and even spectra. In this case, the autocorrelation of the island mode spectra shows not only a peak at the large separation but also at half the large separation. This is due to the fact that the island modes are built around a central orbit which is selfretracing, i.e. during a complete period it goes twice through the same points in q with opposite momenta. For such orbit, the semiclassical quantization for an even or odd spectrum uses a twice shorter orbit than for the full spectrum. In contrast, the odd and even spectra of chaotic modes are built on generic orbits with no such property, and are completely independent. Hence, the autocorrelation of the full chaotic spectrum at a given rotation, with both parities, does not show a peak at half separation.
5. Discussion and conclusion
In this paper, we have computed high frequency pmodes in stellar polytropic models for rotation rates between Ω/Ω_{k} = 0.48 and Ω/Ω_{k} = 0.81. Following the methodology of Sect. 2.3, we have then identified chaotic modes and built a dataset of chaotic frequencies. As expected, the nearestneighbor statistics of the chaotic spectra for most rotations follow the WignerDyson distribution, a wellknown generic property of wave chaos systems. The frequency autocorrelations of the chaotic spectra have been computed. All of them exhibit peaks above the noise level. The presence of peaks in the frequency autocorrelation is clearly not generic in wave chaos systems.
Our analysis shows that chaotic modes are organized in series. The frequency difference between consecutive modes being approximately constant and of similar value across all series. By displaying chaotic mode frequencies on échelle diagrams, we showed that this weakly varying frequency interval can be interpreted as a pseudo large separation. We speak about pseudo large separation because contrary to modes in the nonrotating case or to island modes, the interval is slightly irregular and would remain so asymptotically.
The pseudo large separation is responsible for the presence of the socalled main peak in the frequency autocorrelations and we explained it using semiclassical methods. The ray dynamics is indeed peculiar, as the sound speed is strongly inhomogeneous along the radius of the star. We characterized the impact of the sound speed profile on the ray dynamics through two variables, σ_{0} and T_{0} (see Eq. (18)). These two quantities correspond respectively to the mean value and standard deviation of the onechord travel time distribution. Using the formalism of quantum chaos, we then wrote a semiclassical expression of the autocorrelation and showed that knowing σ_{0} and T_{0} is enough to recover the position of the main peak in the numerically computed spectra. This asymptotic analysis also explains the decrease of the peak height as rotation increases, equivalent to a loss of regularity of chaotic modes. There are other peaks in the autocorrelations, which vary in a less predictable way as a function of rotation. We propose that they are created by the presence of phase space structures that develop around the stable island chains, called partial barriers.
The large separation is expressed as Δ_{c} ∼ 2π/T_{0}. Since chaotic trajectories cover the entire meridional plane, at least in the asymmetric case, we found that T_{0} can be estimated approximately, without using ray tracing, by computing the average acoustic time over the meridional plane (see Eq. (19)).
As explained in Sect. 4.5, the small variance of acoustic travel times (σ_{0}/T_{0} ≪ 1) implies the observed quasidegeneracy of Δ_{i} and Δ_{c}. Thus, we expect rapidly rotating stars to be characterized by a unique large separation Δ ∼ Δ_{i} ∼ Δ_{c}. Autocorrelation peaks at the large separation detected in δ scuti stars (García Hernández et al. 2015) could be produced not only by island modes, as previously thought, but also in part by chaotic modes. This would be important for stars that rotate rapidly enough to harbor a significant number of chaotic modes. A way to distinguish the contribution of island modes is to look for a peak at half the large separation as our analysis indicates that it is due to island modes only.
To go further in the comparison with observed spectra, one should construct a database of lowfrequency synthetic spectra as in Reese et al. (2017) but with full mode identification and with a higher sampling in rotation rate. We expect the asymptotic properties described in the present paper to guide the identification of chaotic modes even at low frequency, as it was the case for island modes. Having fully identified synthetic spectra then would help to identify modes in real data.
Calculations of mode visibilities in rapid rotators were performed first in Lignières & Georgeot (2009) and later in Reese et al. (2013) taking more effects into account such as gravity darkening. These calculations showed that the surface structure of chaotic modes should allow them to be visible, especially at high rotation rates. The most direct proof of the occurrence of wave chaos in stars would be to identify a large set of observed chaotic mode frequencies and find that they follow closely the WignerDyson surmise. This is a very difficult task, however a couple of observations might reduce the difficulty by a small amount. First, the spectra of many stars with various rotation rates can be aggregated, as we did in Fig. 7 to construct the nearest neighbor distribution. Also, choosing very fast rotators may help, since chaotic modes are expected to be more present at very high rotation rates. Finally, a poleon configuration can be helpful as avoiding m ≠ 0 subspectra would simplify mode identification.
Looking closely, one can discern a very narrow peak inside each packet. It is created by trajectories trapped near the main island (see Sect. 4.2) and will have no significant impact on the autocorrelation.
Acknowledgments
We thank CALMIP (“CALcul en MIdiPyrénées”) for the use of their supercomputer. We used the code Top developed by D. Reese and made userfriendly by B. Putigny. We thank ISSI (“International Space Science Institute”) through the SoFAR (“Seismology of Fast Rotating Stars”) program for their support.
References
 Aerts, C., ChristensenDalsgaard, J., & Kurtz, D. W. 2010, in Asteroseismology (Netherlands: Springer), Astron. Astrophys. Lib. [Google Scholar]
 Ballot, J., Lignières, F., & Reese, D. R. 2013, in Numerical Exploration of Oscillation Modes in Rapidly Rotating Stars, eds. M. Goupil, K. Belkacem, C. Neiner, F. Lignières, & J. J. Green, 865, 91 [Google Scholar]
 Berry, M. V. 1985, Proc. R. Soc. London A, 400, 229 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Berry, M. V., & Robnik, M. 1984, J. Phys. A: Math. Gen., 17, 2413 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Bogomolny, E., & Hugues, E. 1998, Phys. Rev. E, 57, 5404 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Bogomolny, E. B., & Keating, J. P. 1996, Phys. Rev. Lett., 77, 1472 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Bogomolny, E. B., Georgeot, B., Giannoni, M.J., & Schmit, C. 1992, Phys. Rev. Lett., 69, 1477 [NASA ADS] [CrossRef] [Google Scholar]
 Bogomolny, E. B., Georgeot, B., Giannoni, M.J., & Schmit, C. 1997, Phys. Rep., 291, 219 [NASA ADS] [CrossRef] [Google Scholar]
 Bohigas, O. 1991, in Random Matrix Theories and Chaotic Dynamics, eds. M. J. Giannoni, A. Voros, & J. ZinnJustin (Amsterdam: NorthHolland), Proceedings of the Les Houches Summer School of Theoretical Physics, LII, 87 [Google Scholar]
 Bohigas, O., Tomsovic, S., & Ullmo, D. 1993, Phys. Rep., 223, 43 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Bowman, D. M., & Kurtz, D. W. 2018, MNRAS, 476, 3169 [NASA ADS] [CrossRef] [Google Scholar]
 Chang, S.J., & Shi, K.J. 1986, Phys. Rev. A, 34, 7 [NASA ADS] [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
 CohenTannoudji, C., Diu, B., & Laloe, F. 1977, Mecanique quantique (Paris: Hermann) [Google Scholar]
 Cvitanovic, P., Artuso, R., Mainieri, R., Tanner, G., & Vattay, G. 2017, Chaos: Classical and Quantum, http://chaosbook.org/ [Google Scholar]
 Evano, B., Georgeot, B., & Lignières, F. 2019, EPL, 125, 49002 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 García Hernández, A., Moya, A., Michel, E., et al. 2009, A&A, 506, 79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 García Hernández, A., Moya, A., Michel, E., et al. 2013, A&A, 559, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 García Hernández, A., MartínRuiz, S., Monteiro, M. J. P. F. G., et al. 2015, ApJ, 811, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Gough, D. 1993, in Linear Adiabatic Stellar Pulsation, eds. J. P. Zahn, & J. ZinnJustin (Amsterdam: Elsevier), Proceedings of the Les Houches Summer School of Theoretical Physics, XLVII, 400 [Google Scholar]
 Gutzwiller, M. C. 1990, Chaos in Classical and Quantum Mechanics, Interdisciplinary Applied Mathematics (New York: SpringerVerlag) [CrossRef] [Google Scholar]
 Hannay, J. H., & Ozorio De Almeida, A. M. 1984, J. Phys. A: Math. Gen., 17, 3429 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Hansen, C. J., Kawaler, S. D., & Trimble, V. 2004, in Stellar Interiors: Physical Principles, Structure, and Evolution, 2nd edn. (New York: SpringerVerlag), Astron. Astrophys. Lib. [Google Scholar]
 Kudrolli, A., Sridhar, S., Pandey, A., & Ramaswamy, R. 1994, Phys. Rev. E, 49, R11 [NASA ADS] [CrossRef] [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]
 Lignières, F., Rieutord, M., & Reese, D. 2006, A&A, 455, 607 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mehta, M. 2004, Random Matrices (Amsterdam: Elsevier) [Google Scholar]
 Michel, E., Dupret, M.A., Reese, D., et al. 2017, Eur. Phys. J. Web Conf., 160, 03001 [Google Scholar]
 Mirouh, G. M., Angelou, G. C., Reese, D. R., & Costa, G. 2019, MNRAS, 483, L28 [NASA ADS] [CrossRef] [Google Scholar]
 Ott, E. 1993, Chaos in Dynamical Systems (Cambridge: 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]
 Ouazzani, R.M., Roxburgh, I. W., & Dupret, M.A. 2015, A&A, 579, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 O’Connor, P., Gehlen, J., & Heller, E. J. 1987, Phys. Rev. Lett., 58, 1296 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Paparó, M., Benkö, J. M., Hareter, M., & Guzik, J. A. 2016, ApJS, 224, 41 [NASA ADS] [CrossRef] [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]
 Reese, D., Lignières, F., & Rieutord, M. 2006, A&A, 455, 621 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reese, D., Lignières, F., & Rieutord, M. 2008, A&A, 481, 449 [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]
 Reese, D. R., Prat, V., Barban, C., VeerMenneret, C. V. T., & MacGregor, K. B. 2013, A&A, 550, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reese, D. R., Lignières, F., Ballot, J., et al. 2017, A&A, 601, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schulman, L. S. 1996, Techniques and Applications of Path Integration (New York: Wiley Classics Library) [Google Scholar]
 Shim, J.B., Wiersig, J., & Cao, H. 2011, Phys. Rev. E, 84, 035202 [NASA ADS] [CrossRef] [Google Scholar]
 Sieber, M., & Richter, K. 2001, Phys. Scr., 2001, 128 [CrossRef] [Google Scholar]
 Vidmar, G., Stöckmann, H.J., Robnik, M., et al. 2007, J. Phys. A: Math. Theor., 40, 13883 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Trace formula
In this appendix we will adapt the semiclassical formalism used in quantum mechanics to derive a trace formula for chaotic modes in rotating stars. The general derivation follows the original one due to Gutzwiller, detailed e.g. in (Gutzwiller 1990; Cvitanovic et al. 2017; Ott 1993).
The Hamiltonian system
We denote the canonically conjugate variables (p, q), as is usual in textbooks on Hamiltonian mechanics. In the present subsection, we will consider the system to be onedimensional to avoid carrying indices in the notation. The ray dynamics is governed by the following Hamiltonian:
where the wave vector plays the role of momentum: p ≡ k.
Hamilton’s principal function
We now express p as a function of :
The ray system can thus be seen as analogous to a mechanical system with a varying mass . Then, the Lagrangian L is obtained through the usual Legendre transform:
Hamilton’s principal function, denoted R, is defined as the time integral of the Lagrangian. Its computation involves following a trajectory from (q′,t′) to (q, t). For short times δt, it gives:
using .
Action integral
The action integral S is defined from R as . From Eq. (A.3) we have:
which is the well know expression of the action. Finally, the eikonal equation Eq. (2) gives . Thus the action can be written in terms of the acoustic time:
where s is the curvilinear coordinate along the ray path.
WKB approximation for the semiclassical propagator
We first derive the expression for the semiclassical propagator, adapting the method in (Cvitanovic et al. 2017) for quantum systems to our star model, which has two degrees of freedom and a fourdimensional phase space. Let Λ^{−1} be a small dimensionless parameter. In the wave equation, Eq. (1), we insert the WKB ansatz Ψ(q, t) = A e^{iΛϕ(q, t)} leading to (Gough 1993):
Equation (A.7) can be recognized as the HamiltonJacobi equation ∂R/∂t = ±H. Thus, in the limit of small wavelengths, the phase is simply Hamilton’s principal function: ϕ = R or ϕ = −R (associated to the Hamiltonian −H). The two possible phases lead to two terms in the propagator, with a projector on each subspace P_{1} and P_{2} satisfying P_{1} + P_{2} = I.
Additionally, the substitution yields:
Introducing the density ρ = A^{2} and velocity v = 1/m (∂R/∂q), with (from the mechanical analogy introduced in Sect. A.1.1), Eq. (A.9) appears as a continuity equation. It follows that the ray amplitude A(q, t) can be interpreted as the square root of the density of nearby trajectories. The evolution of this density from q′ to q (variation of volume in coordinate space of a swarm of trajectories) is quantified through the Jacobian determinant .
Having now both the phase and amplitude, we obtain the semiclassical wave function
where A_{1} and A_{2} are the projection of the intial wave function on the two subspaces. This is valid for short time, i.e. when only one classical trajectory connects q′ to q in time t. For longer times, several trajectories which we label by j connect the two points, and the formula becomes:
where the topological index κ is added to account for the phase shift at points where the amplitude becomes singular, such as caustics.
The propagator K is defined by
The propagator is the timedependent Green’s function, it is thus solution of the wave equation with the initial condition lim_{t → 0}K(q, q′,t) = δ(q − q′). Again, we will assume that, for short times δt, the semiclassical propagator is of the form (for the first term, the computation is similar for the second term)
Neglecting the second term ω δt in (A.4) gives
If we impose , then the previous expression is a 2dimensional Gaussian of width . In the limit δt = 0 the condition K(q″, q′,t = 0) = δ(q − q′) is satisfied. Using again the substitution for clarity
using the fact that the expression becomes
This short time expression with the correct limit at t → 0 can be extended to longer times using Eq. (A.11) and the fact that . This gives:
where the sum is over all classical trajectories labeled by j from q′ to q in time t.
Green’s function
To derive the trace formula, the usual procedure necessitates in quantum mechanical systems to go from the propagator to the energydependent Green function G(q, q′,E), which is related to the propagator through the Fourier transform
where ϵ is a small positive number which makes the integral convergent and goes to zero eventually. Then one uses the fact that G(q, q′,E) can be expanded on a basis of eigenvector {ϕ_{j}} of the Hamiltonian with eigenvalues E_{n} as
In the semiclassical approximation this leads to the following equality for E/ℏ → ∞
which is valid at first order in ℏ (acoustic case: first order in Λ^{−1}). One can already forecast the trace formula from this equality, as it involves both the eigenenergies of the quantum system on the lefthand side and classical quantities in the righthand side. We now need to find out an expression similar to Eq. (A.20) for the acoustic waves.
The propagator can be found by taking the matrix element of the evolution operator U(t, 0) between the final and initial states q⟩ and q′⟩, (see e.g. CohenTannoudji et al. 1977):
U(t, 0) satisfies the wave equation: , where is the Hamiltonian operator with . Then:
with P_{1} and P_{2} projectors in two subspaces as above. Let ϕ_{n}⟩ be the eigenfunctions of the Hamiltonian . Then, from the closure relation, we have:
where A_{1} and A_{2} are the projection of the intial wave function on the two subspaces. Using Eq. (A.23) in Eq. (A.21) gives:
The ω dependent Green’s function G(q″, q′,ω) stems from the Fourier transform of K(q, q′,0):
On the other hand, the semiclassical Green’s function is obtained by taking the Fourier transform of the semiclassical propagator (A.16) and evaluating it by stationary phase; one starts from
The phase term of the integrand is the action R(q, q′,t)+ωt = S(q, q′,ω) or −R(q, q′,t)+ωt = S′(q, q′,ω). Stationary points of the first sum are such that which correspond to classical trajectories from q′ to q at frequency ω. As usual one expands the integrand in Eq. (A.26) in powers of t at second order. Then the integral is approximated by the method of stationary phase (Schulman 1996). Again, the computation of the prefactor requires to combine the prefactor of (A.17) with the one coming from the stationary phase. The second sum has stationary points at or . The result is:
where the sum is over all classical trajectories from q′ to q at frequency ω, and are final and initial velocities, and q_{⊥} and are coordinates transverse to the orbit. The index κ′_{j} counts again the singularities along the orbits.
The final formula
To obtain the trace formula, we should compute the trace of the Green’s function from the two formulas we obtained, Eqs. (A.25) and (A.27). We will now keep only the first part in both equations since each term in one equation is equal to its counterpart in the other. Let us first compute the trace of Green’s function from Eq. (A.25). To this aim, let the small imaginary part in the denominator go to zero and get for the imaginary part:
On the other hand, the trace of the semiclassical Green’s function Eq. (A.27) is:
This formula is an integral involving all closed classical paths from q to q at frequency ω. It contains two parts. The first one corresponds to the limit for q → q′ of the short direct trajectories between q and q′, which become of zero length. We call it Tr G_{0} and it should be treated separately. The remaining contains a sum of closed orbits between q and q with nonzero length. We evaluate this sum again by stationary phase; the stationary points in the sum are such that the first derivative of the function in the exponential is cancelled. This implies that for q = q′. This selects closed trajectories with equal initial and final momentum, thus periodic orbits. Again, the prefactors are to be combined correctly, yielding to:
where i labels all the periodic orbits of the system, including repetitions of a primitive orbit, M_{i} is the monodromy matrix describing the linearized motion in the transverse direction to the orbit; the determinant encodes the stability of this orbit i. T_{i} is the geometrical period of the orbit (i.e. without counting the repetitions).
Equation (A.28) connects the trace of the Green’s function to the density of states ∑_{j}δ(ω − ω_{j}). It is known that this sum can be split in two parts . The first term corresponds to the smooth part of the density of states, while the second part contains the fluctuating (oscillatory) part. It turns out that corresponds to Tr G_{0}, while the oscillatory part corresponds to the remaining part of Eq. (A.30). Putting together Eq. (A.28) for positive frequencies and Eq. (A.30) gives the Gutzwiller formula for the acoustic waves:
All Tables
Comparison of the large separation of island modes Δ_{i} and chaotic modes Δ_{c}, obtained from the simulated spectra at six rotations.
All Figures
Fig. 1. PSS at rotation rate Ω/Ω_{k} = 0.589 with , where θ is the colatitude and k_{θ} the projection of the wave vector on the line tangent to the r_{p}(θ) = r_{s}(θ)−z curve. Two stable islands are embedded in a chaotic region, which is itself surrounded by whispering gallery rays. The different types of rays are color coded, blue: whispering gallery, green: 6period island, yellow: 2period island and red: ergodic. 

In the text 
Fig. 2. Chaotic zone of the PSS at increasing values of the rotation rate. The chaotic zone grows monotonically, which is not the case for the main island zone. 

In the text 
Fig. 3. PSS at Ω/Ω_{k} = 0.589 with , where θ is the colatitude and k_{θ} the projection of the wave vector on the line tangent to the r_{p}(θ) = r_{s}(θ)−z curve. The phase space structures are similar to those presented in Fig. 1. 

In the text 
Fig. 4. Four odd axisymmetric modes at rotation Ω/Ω_{k} = 0.589: (a) chaotic mode, (b) whispering gallery mode, (c) 2period island mode (ℓ = 0), with a black line indicating the central periodic orbit, and (d) 6period island mode. The figure shows the scaled pressure amplitude , with d the distance to the rotation axis and ρ_{0} the equilibrium density. 

In the text 
Fig. 5. Histograms showing the distribution of three types of odd modes as a function of their dominant degree l_{max}, at Ω/Ω_{k} = 0.589 in the frequency domain 25.6 ω_{p} to 33.5 ω_{p}. The existence of a gap between the highest value of l_{max} for chaotic modes, marked by a dashed line, and the lowest value of l_{max} for whispering gallery modes is used to automatically separate these two types of modes. 

In the text 
Fig. 6. Three modes at Ω/Ω_{k} = 0.589 belonging to series 1. From left to right: we see a chaotic mode, an interface mode and an island mode with ℓ = 4. The thick white line is the main island central orbit. The modes intensity Ψ^{2} is represented, where Ψ is the scaled pressure amplitude as in Fig. 4. 

In the text 
Fig. 7. Left: nearest neighbors spacing distribution , with 1344 frequency levels obtained from eight independent spectra: Ω/Ω_{k} = 0.481 (206 odd levels), Ω/Ω_{k} = 0.545 (223 odd levels, 105 even levels), Ω/Ω_{k} = 0.589 (217 odd levels, 96 even levels), Ω/Ω_{k} = 0.658 (207 odd levels, 120 even levels) and Ω/Ω_{k} = 0.809 (170 odd levels). Right: integrated distribution N(s) for all eight independent spectra. In both panels the dashed line the Wigner surmise and the dotted line is the prediction for Poissonian spectra. 

In the text 
Fig. 8. Scaled pressure amplitude Ψ shown for twelve odd axisymmetric chaotic modes at Ω/Ω_{k} = 0.589 with quantum number m = 0. 

In the text 
Fig. 9. Autocorrelations R_{2}(ξ), where ξ is a displacement in frequency, of chaotic spectra with odd parity: (a) 206 levels from 28.35 ω_{p} to 46.89 ω_{p}, (b) 223 levels from 28.15 ω_{p} to 44.09 ω_{p}, (c) 217 levels from 26.02 ω_{p} to 40.29 ω_{p}, (d) 283 levels from 23.57 ω_{p} to 36.22 ω_{p} and (e) 170 levels from 24.02 ω_{p} to 30.01 ω_{p}. The autocorrelations have been normalized such that their value at the origin is unity. The main peak’s position is labeled Δ_{c} and marked with a solid line. Secondary peaks are marked in panels c, d and e with a dashed line. 

In the text 
Fig. 10. Échelle diagram of chaotic modes at Ω/Ω_{k} = 0.589 in the range 25.60 ω_{p} to 33.54 ω_{p}, with odd parity. 

In the text 
Fig. 11. Échelle diagram of island modes ℓ = 0, 1, 2, 3, 4 at Ω/Ω_{k} = 0.589, with only odd modes. The lowest point of the ℓ = 4 track corresponds to an islandlike interface mode. 

In the text 
Fig. 12. Mode intensity Ψ^{2} at Ω/Ω_{k} = 0.589, where Ψ is the scaled pressure amplitude, showing the similarity between consecutive modes. Top: two consecutive modes that belong to series 8. Bottom: two consecutive modes that belong to series 3. 

In the text 
Fig. 13. Échelle diagram showing all chaotic frequencies in the range 25.60 ω_{p} to 33.54 ω_{p} modulo the large separation Δ_{c}. Series of modes are labeled by numbers 1 to 16. Black dots correspond to modes that do not fit into a series. 

In the text 
Fig. 14. Top panel: autocorrelation at Ω/ω_{k} = 0.589 for quantum number m = 1, in the frequency domain 30.51 ω_{p} to 38.48 ω_{p}. Bottom panel: autocorrelation at Ω/ω_{k} = 0.589 for quantum number m = 4, in the frequency domain 30.53 ω_{p} to 38.51 ω_{p}. In both cases, the dashed line is the position Δ_{c} of the main peak for axisymmetric modes at the same rotation rate. 

In the text 
Fig. 15. Comparison of the “stellar” autocorrelation R_{2}(ξ) at Ω/Ω_{k} = 0.589, represented in solid line, with the GOE autocorrelation represented by a dashed line. To compare the “stellar” and GOE autocorrelations, two changes have been made. First, the autocorrelation of our model has been rescaled so that the value of the plateau, i.e. the line around which there are fluctuations, is approximately unity. The value of the plateau otherwise depends on the mean density and the Gaussian smoothing of the spectrum. Secondly, the GOE distribution has been rescaled, in the direction of the x axis, by the mean level spacing of the Ω/Ω_{k} = 0.589 spectrum in the frequency domain considered. 

In the text 
Fig. 16. Left panels: number of nchord trajectories, with n = 1, …, 20, vs their travel time T at rotations Ω/Ω_{k} = 0.481 (top) and Ω/Ω_{k} = 0.809 (bottom), with 300 bins for the total distribution. The nchord samples contain ∼4200 chords each. Right panel: Ω/Ω_{k} = 0.481 1chord distribution in more details, with ∼84 000 chords and 100 bins. The mean value T_{0} of the distribution is marked with a dashed line and the standard deviation σ_{0} is shown. 

In the text 
Fig. 17. Upper panel: theoretical large separation (upward triangles) compared to the numerical peak’s position (downward triangles) calculated for axisymmetric modes at six rotation rates. Middle panel: theoretical autocorrelations with quantum number m = 0, from right to left the rotation rate increases. Bottom panel: theoretical autocorrelations at Ω/Ω_{k} = 0.589 with m = 0, 2, 4, 6 and frequency ω = 24.41 ω_{p}. 

In the text 
Fig. 18. Snapshot of the evolution of a bundle of trajectories as they intersect the PSS, represented by black dots, at Ω/Ω_{k} = 0.706. The phase space zones where black dots are dense correspond to regions enclosed by partial barriers. The trajectories are initially in the neighborhood of the two main islands. One of the central periodic orbits of these zones is indicated by colored dots (large gray dots). 

In the text 
Fig. 19. Superposition of random stationary waves reproducing the kind of amplitude patterns seen in the chaotic pmodes of the polytropic stellar model (see text). 

In the text 
Fig. 20. Intersection of a few tori with the PSS at Ω/Ω_{k} = 0. Only the part of the PSS where k_{θ} > 0 is shown because of the symmetry with respect to the k_{θ} = 0 axis at zero rotation. The tori correspond to modes of degree ℓ_{s} = 1, ℓ_{s} = 9 or ℓ_{s} = 11. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.