© B. Evano et al. 2019

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 solar-like 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 spherical-harmonic 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 high-frequencies for acoustic modes. This motivates high-frequency 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 2-period 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 semi-analytical 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 high-frequency 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 auto-correlation 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 high-frequency pressure modes and finally the method used to classify the modes and obtain a set of high-frequency 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 self-gravitating gas. We are focusing on p-modes in the high-frequency 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 high-frequency 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 high-frequency regime, see Reese et al. (2006). Finally, we can discard the gravity waves since the Brunt-Väisälä frequency is very small compared to the high p-mode frequencies. With all these approximations taken into account, the adiabatic pressure perturbations obey a Helmholtz-like equation derived in Lignières & Georgeot (2009):

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{\Psi }+\frac{{\omega }^2-{\omega }_{\mathrm{c}}^2}{c_{\mathrm{s}}^2}\mathrm{\Psi }=0,\end{array}$$(1)

where $\mathrm{\Psi }=\widehat{P}/f$, with $\widehat{P}$ the complex amplitude associated with the time-harmonic pressure perturbation $P=\mathit{Re}[\widehat{P}exp(-i\omega t)]$, f a function of the background model, ω_c is a cut-off frequency of the model and c_s is the inhomogeneous sound velocity.

The short-wavelength approximation of Eq. (1) leads to the eikonal equation

$$\begin{array}{c}\hfill {\omega }^2=c_{\mathrm{s}}^2{k}^2+{\omega }_{\mathrm{c}}^2.\end{array}$$(2)

This equation can be put in the form of a Hamiltonian system describing the propagation of acoustic rays (Ott 1993), by setting $H=\omega ={(c_{\text{s}}^2{k}^2+{\omega }_{\text{c}}^2)}^{1/2}$. The motion takes place in a meridional plane rotating with the ray at an angular velocity $\mathrm{d}\varphi /\mathrm{d}t={\stackrel{\sim }{L}}_z/d^2$, where ${\stackrel{\sim }{L}}_z=rsin\theta k_{\varphi }/\omega$ 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:

$$\begin{array}{cc}\hfill H\prime & =\frac{{\stackrel{\mathbf{\sim }}{\mathit{k}}}_{\mathit{p}}^2}{2}+W,\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill W& =-\frac{1}{2c_{\mathrm{s}}^2}(1-\frac{{\omega }_{\mathrm{c}}^2}{{\omega }^2})+\frac{{\stackrel{\sim }{L}}_z^2}{2d^2},\hfill \end{array}$$(4)

where ${\stackrel{\mathbf{\sim }}{\mathit{k}}}_{\mathit{p}}$ is the frequency-scaled wave vector $\stackrel{\mathbf{\sim }}{\mathit{k}}=\mathit{k}/\omega$ projected onto the meridional plane and W is the potential. The corresponding dynamical system is

$$\begin{array}{c}\hfill \frac{\mathrm{d}\mathit{x}}{\mathrm{d}t}={\stackrel{\mathbf{\sim }}{\mathit{k}}}_{\mathit{p}},\frac{\mathrm{d}{\stackrel{\mathbf{\sim }}{\mathit{k}}}_{\mathit{p}}}{\mathrm{d}t}=-\mathbf{\nabla }W,\end{array}$$(5)

where $\mathrm{d}t=c_{\mathrm{s}}\mathrm{d}s/{(1-\frac{{\omega }_{\mathrm{c}}^2}{{\omega }^2})}^{1/2}$, is a time-like variable, s being the curvilinear coordinate along the ray. This system is then integrated using a fifth order Runge-Kutta method. To simplify the notation, we now refer to ${\stackrel{\mathbf{\sim }}{\mathit{k}}}_{\mathit{p}}$ as simply $\stackrel{\mathbf{\sim }}{\mathit{k}}$.

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 ${\stackrel{\sim }{L}}_z=0$, where θ is the colatitude and kθ the projection of the wave vector on the line tangent to the rp(θ) = rs(θ)−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: 6-period island, yellow: 2-period island and red: ergodic.

The system undergoes a KAM-like 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 break-up rotation rate ${\Omega }_{\text{k}}={(GM/R_{\text{eq}}^3)}^{1/2}$, 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 2-period 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 non-monotonic 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 ${\stackrel{\sim }{L}}_z\ne 0$, 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) $-arcsin(|{\stackrel{\sim }{L}}_z|/\stackrel{\sim }{L})<\theta <arcsin$ $(|{\stackrel{\sim }{L}}_z|/\stackrel{\sim }{L})$, $\stackrel{\sim }{L}$ being the frequency-scaled angular momentum norm.

Fig. 3. PSS at Ω/Ωk = 0.589 with ${\stackrel{\sim }{L}}_z=0.16/{\omega }_{\mathrm{p}}$, where θ is the colatitude and kθ the projection of the wave vector on the line tangent to the rp(θ) = rs(θ)−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 high-frequency modes and classifying them according to the phase space structure.

2.2. The set of numerically computed high-frequency p-modes

As for the ray dynamics, the star is modeled by a uniformly rotating, self-gravitating 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 high-frequency, and thus small-wavelength, p-modes, 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 $Y_l^m$, through the decomposition:

$$\begin{array}{c}\hfill \mathrm{\Psi }(r,\theta ,\varphi )={\displaystyle \sum _{l=0}^{N_l}}\left({\displaystyle \sum _{j=0}^{N_j}}{a}_j^{l,m}{T}_j(2r-1)\right)Y_l^m(\theta ,\varphi ),\end{array}$$(6)

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 ${\stackrel{\sim }{k}}_{\theta }^{\max }$ or ${\stackrel{\sim }{k}}_{\mathrm{r}}^{\max }$ 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 high-frequency 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 ${\omega }_{\text{p}}={(GM/R_{\text{p}}^3)}^{1/2}$, 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 6-period island mode corresponding to the 6-period island chain, a 2-period island mode corresponding to the 2-period 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) 2-period island mode (ℓ = 0), with a black line indicating the central periodic orbit, and (d) 6-period island mode. The figure shows the scaled pressure amplitude $\mathrm{\Psi }=P\sqrt{d/{\rho }_0}$, 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 2-period 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 ${\delta }_n(m)=2\pi /{\oint }_{\gamma }(\mathrm{d}s/{\stackrel{\sim }{c}}_{\mathrm{s}})$ ¹. However, the regularity of the 2-period 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, 6-period 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 ${\stackrel{\sim }{k}}_{\theta }^{\max }$ of any whispering gallery ray is higher than the ${\stackrel{\sim }{k}}_{\theta }^{\max }$ 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 ${\overline{a}}^l\equiv \frac{1}{N_j}\sum _{j=0}^{N_j}|a_j^l|$. Thus the dominant degree l_max corresponds to the value of l where ${\overline{a}}^l$ 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 lmax, 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 lmax for chaotic modes, marked by a dashed line, and the lowest value of lmax for whispering gallery modes is used to automatically separate these two types of modes.

At finite frequency, here ω ≤ 47ω_p, the Berry-Robnik 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 gallery-like and island-like. An example of island-like 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 clear-cut, 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 gallery-like interface mode. The difference between chaotic, island-like 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 $P(\stackrel{\sim }{s})$, with ${\stackrel{\sim }{s}}_n=({\omega }_{n+1}-{\omega }_n)/⟨{\omega }_{n+1}-{\omega }_n⟩$, where ⟨ω_n + 1 − ω_n⟩ corresponds to the mean frequency spacing, is expected to follow the Wigner surmise of Random Matrix Theory $P(\stackrel{\sim }{s})=\frac{\pi \stackrel{\sim }{s}}{2}{e}^{-\pi s^2/4}$ (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 $N(s)={\int }_0^{s}P(\stackrel{\sim }{s})\mathrm{d}\stackrel{\sim }{s}$, 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).

Fig. 7. Left: nearest neighbors spacing distribution $P(\stackrel{\sim }{s})$, 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 high-frequency 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 p-modes 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₂(ξ) = ∫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 sub-sections) and then extend our analysis to the non-axisymmetric modes.

Fig. 9. Autocorrelations R2(ξ), 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 so-called é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 well-known for low-degree modes in non-rotating stars as well as island modes in rotating ones (see Fig. 11). The frequency spacing between two consecutive modes is the so-called 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 so-called Einstein-Brillouin-Keller 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 island-like 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 ², we were able to label the vertical series of chaotic mode échelle diagram, ending-up 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 non-rotating and island modes. First the frequency spacing is much more regular for non-rotating 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 2-period 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 gallery-like modes or at the end of a series as island-like 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. Non-axisymmetric 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 ${\stackrel{\sim }{L}}_z=L_z/\omega$ 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 time-reversal 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 R2(ξ) 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 re-scaled 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 $\overline{d}(\omega )=\mathrm{d}\overline{N}/\mathrm{d}\omega$ 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 $\overline{N}(\omega )$ (or $\overline{d}(\omega )$) 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 re-derive the trace formula for the system considered here in Appendix A. Using j as a label for the periodic orbits γ_j, the formula reads:

$$\begin{array}{c}\hfill d(\omega )-\overline{d}(\omega )=\mathit{Re}{\displaystyle \sum _j}{A}_j{e}^{i{S}_j(\omega )},\end{array}$$(7)

where $S_j=\omega {\oint }_{{\gamma }_j}\mathrm{d}s/{\stackrel{\sim }{c}}_{\mathrm{s}}$ is the action and ${\oint }_{{\gamma }_j}\mathrm{d}s/{\stackrel{\sim }{c}}_{\mathrm{s}}=T_j$ is the acoustic travel time of γ_j. The amplitude term is expressed as

$$\begin{array}{c}\hfill A_j=\frac{1}{\pi }\frac{T_j}{|det(M_j-I){|}^{1/2}},\end{array}$$(8)

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 so-called primitive orbits or repetitions of them. In many systems, the density ρ(T) of orbits with a time period T grows exponentially as

$$\begin{array}{c}\hfill \rho (T)\approx (1/T)e^{\lambda T},\end{array}$$(9)

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

$$\begin{array}{c}\hfill A(T)\approx \frac{1}{\pi }T{e}^{-(1/2)\lambda T}.\end{array}$$(10)

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₂(ξ) = ⟨d(ω−1/2 ξ) d(ω+1/2 ξ)⟩, with the average ⟨f⟩ = ∫fdω, can be re-written as $R_2(\xi )=\langle d^{\mathrm{osc}}(\omega -1/2\xi )d^{\mathrm{osc}}(\omega +1/2\xi )\rangle +\varphi (\xi )$, where ϕ(ξ) is a smooth function. The theory aims at deriving an approximate expression of the form factor:

$$\begin{array}{c}\hfill K(T)=\frac{1}{\sqrt{2\pi }}{{\displaystyle \int }}_{-\infty }^{\infty }d\xi \exp (i\xi T)\mathrm{C}(\xi ),\end{array}$$(11)

which is the Fourier transform of the autocorrelation

$$\begin{array}{cc}\hfill C(\xi )& =\langle d^{\mathrm{osc}}(\omega -1/2\xi )d^{\mathrm{osc}}(\omega +1/2\xi )\rangle ,\hfill \end{array}$$(12)

$$\begin{array}{cc}& =\langle (\mathit{Re}{\displaystyle \sum _i}{A}_i{e}^{i{S}_i(\omega -1/2\xi )})(\mathit{Re}{\displaystyle \sum _j}{A}_j{e}^{i{S}_j(\omega +1/2\xi )})\rangle .\hfill \end{array}$$(13)

Getting rid of the real parts through the fact that with $x=\mathit{Re}(\stackrel{\sim }{x})=A_{i}cos(\xi T_i+S_i)$ and $y=\mathit{Re}(\stackrel{\sim }{y})=A_{j}cos(\xi T_j+S_j)$ then $⟨xy⟩=1/2⟨\stackrel{\sim }{x}{\stackrel{\sim }{y}}^{\ast }⟩$, gives

$$\begin{array}{c}\hfill C(\xi )=\frac{1}{2}\langle ({\displaystyle \sum _i}{A}_i{e}^{i{S}_i(\omega -1/2\xi )})({\displaystyle \sum _j}{A}_j{e}^{i{S}_j(\omega +1/2\xi )}{)}^{\ast }\rangle ·\end{array}$$(14)

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:

$$\begin{array}{cc}\hfill K(T)& =\frac{1}{2}\langle {\displaystyle \sum _i\sum _j}{A}_i{A}_j^{\ast }{e}^{i(S_i-S_j)}\frac{1}{\sqrt{2\pi }}{{\displaystyle \int }}_{-\infty }^{\infty }d\xi e^{i\xi (T-\frac{1}{2}(T_i+T_j))}\rangle ,\hfill \\ \hfill & =\frac{1}{2}\langle {\displaystyle \sum _i\sum _j}{A}_i{A}_j^{\ast }{e}^{i(S_i-S_j)}\frac{2\pi }{\sqrt{2\pi }}\delta (T-\frac{1}{2}(T_i+T_j))\rangle ·\hfill \end{array}$$(15)

Under the frequency average ⟨.⟩, the contribution of the off-diagonal 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

$$\begin{array}{c}\hfill \mathrm{K}(T)\approx {\displaystyle \sum _{\mathrm{j}}}{A}_{\mathrm{j}}^2\delta (T-T_j).\end{array}$$(16)

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

$$\begin{array}{c}\hfill K(T)\approx A^2(T)\rho (T),\end{array}$$(17)

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 n-chord. For any integer n, periodic trajectories with n rebounds are a subset of n-chord 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 n-chord trajectories.

At a given rotation rate, we compute the acoustic travel times of thousands of n-chord trajectories. The 1-chord distribution is shown in the right panel of Fig. 16, at Ω/Ω_k = 0.481. The distribution is a narrow packet: its standard deviation σ₀ is small compared to its mean value T₀. We interpret this effect as a consequence of the strong decrease of the sound speed near the surface³. 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₀ to increase with rotation, because of the increasing volume of the star. This is indeed confirmed, with T₀ = 5.19/ω_p at Ω/Ω_k = 0.481 and T₀ = 6.30/ω_p at Ω/Ω_k = 0.706.

Fig. 16. Left panels: number of n-chord 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 n-chord samples contain ∼4200 chords each. Right panel: Ω/Ωk = 0.481 1-chord distribution in more details, with ∼84 000 chords and 100 bins. The mean value T0 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 n-chords 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₀) but get thicker as n increases (${\sigma }_n\approx \sqrt{n}{\sigma }_0$), 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 σ₀/T₀ which is 0.049 at Ω/Ω_k = 0.481 and 0.097 at Ω/Ω_k = 0.809.

The constraints imposed by the n-chord 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⁴. 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) = ∑_nP_n,  Ω(T), with P_n,  Ω(T) the probability distribution of n-chords travel times:

$$\begin{array}{c}\hfill P_{n,\mathrm{\Omega }}(T)=\frac{T_0}{\sqrt{2\pi n}{\sigma }_0}\exp (-\frac{{(T-n{T}_0)}^2}{2{(\sqrt{n}{\sigma }_0)}^2}),\end{array}$$(18)

where T₀, 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₀ and σ₀. 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²(T) ρ(T) = T P_Ω(T).

The Fourier transform of this quantity

$$\begin{array}{c}\hfill F(\xi )=\frac{1}{\sqrt{2\pi }}{{\displaystyle \int }}_{-\infty }^{\infty }\mathrm{d}Texp(-i\xi T)T{P}_{\mathrm{\Omega }}(T),\end{array}$$

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₂(ξ) (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 ${\mathrm{\Delta }}_{\mathrm{c}}^{\mathrm{th}}$ 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, ${\mathrm{\Delta }}_{\mathrm{c}}^{\mathrm{th}}$ 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₀ and σ₀. The position is found straightforwardly as ${\mathrm{\Delta }}_{\mathrm{c}}^{\mathrm{th}}\approx 2\pi /T_0$. The height and width are controlled by the ratio σ₀/T₀. Indeed, the condition σ₀ ≪ T₀ is necessary to clearly distinguish the packets in the distribution of travel times, as in Fig. 16. Due to the increase of σ₀/T₀ with rotation, the peak gets less visible at high rotation rates.

Fig. 17. Upper panel: theoretical large separation ${\mathrm{\Delta }}_{\mathrm{c}}^{\mathrm{th}}/{\omega }_{\mathrm{p}}$ (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 ${\mathrm{\Delta }}_{\mathrm{c}}^{\mathrm{th}}$ is not only a function of the rotation rate Ω, but also of the projected angular momentum ${\stackrel{\sim }{L}}_z$. In the bottom panel of Fig. 17 we show the predicted peaks for ${\stackrel{\sim }{L}}_z=0,0.0819/{\omega }_{\mathrm{p}},0.1638/{\omega }_{\mathrm{p}}$ and 0.2458/ω_p at Ω/Ω_k = 0.589. The shift in position is governed by the value of T₀ and the change in amplitude by the value of σ₀ as ${\stackrel{\sim }{L}}_z$ increases. The main difference with the axisymmetric case is a reduction of the domain of propagation, resulting in a smaller standard deviation σ₀. As the mean travel time T₀ between two rebounds depends on the impact of ${\stackrel{\sim }{L}}_z$ 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 ${\mathrm{\Delta }}_{\mathrm{c}}^{\mathrm{th}}$ 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 $\tau (\theta )={\int }_0^{r_{\mathrm{s}}(\theta )}\mathrm{d}r/{\stackrel{\sim }{c}}_{\mathrm{s}}$, where r is the radial coordinate, the mean acoustic time between two rebounds can be estimated without using ray dynamics, as

$$\begin{array}{c}\hfill T^{\mathrm{av}}=2{\displaystyle \left(\frac{2}{\pi }{\int }_0^{\pi /2}\tau (\theta )\mathrm{d}\theta \right)},\end{array}$$(19)

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 well-aligned 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 6-period unstable orbit revolves around the main 2-period 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 6-periodic 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 time-harmonic solutions of the Helmholtz equation ∇²ψ + k²ψ = 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 p-modes 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 low-frequency end of some chaotic series, and that island mode series appears as high-frequency 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 Einstein-Brillouin-Keller 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 frequency-scaled coordinates used here, we rather use the invariant $\stackrel{\sim }{L}=L/\omega =\pm ({\ell }_{\mathrm{s}}+1/2)/{\omega }_{n_{\mathrm{s}},{\ell }_{\mathrm{s}},m}$, 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 ${\stackrel{\sim }{k}}_{\theta }=\pm \stackrel{\sim }{L}$ (Lignières & Georgeot 2009). On Fig. 20, the identification of a few mode-carrying 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 $\stackrel{\sim }{L}$ or ${\stackrel{\sim }{k}}_{\theta }$ values. Similarly, for fixed ℓ_s and increasing n_s, the tori approach the ${\stackrel{\sim }{k}}_{\theta }=0$ 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 2-period island chain at low ${\stackrel{\sim }{k}}_{\theta }$, the large chaotic zone at low and intermediate ${\stackrel{\sim }{k}}_{\theta }$, and at large ${\stackrel{\sim }{k}}_{\theta }$, the region of surviving KAM tori corresponding to whispering gallery trajectories. In this context, high-$\stackrel{\sim }{L}$ tori at Ω = 0 will transform into structures of the surviving KAM tori region, whereas low-$\stackrel{\sim }{L}$ tori will be destroyed as the 2-period island chain and the chaotic zone surrounding it emerge. We thus expect high-$\stackrel{\sim }{L}$ modes to become whispering gallery modes, while low-$\stackrel{\sim }{L}$ 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, $\stackrel{\sim }{L}$ 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 low-n_s (high-$\stackrel{\sim }{L}$) modes become whispering gallery modes, the intermediate-n_s modes become chaotic modes, and the high-n_s (low-$\stackrel{\sim }{L}$) modes become island modes. Combining this picture with our observation that some series of chaotic modes show island modes at their high-frequency end or whispering gallery modes at their low-frequency 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 odd-parity 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 ${\stackrel{\sim }{k}}_{\theta }$ 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 $\stackrel{\sim }{L}$ are all close to the ${\stackrel{\sim }{k}}_{\theta }=0$ 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 non-rotating 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 high-frequency 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 ${\mathrm{\Delta }}_{\mathrm{c}}^{\mathrm{th}}=2\pi /T_0$. 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 ${\delta }_n=2\pi /{\oint }_{\gamma }(\mathrm{d}s/{\stackrel{\sim }{c}}_{\mathrm{s}})$. 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₀.

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 self-retracing, 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 p-modes 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 nearest-neighbor statistics of the chaotic spectra for most rotations follow the Wigner-Dyson distribution, a well-known 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 non-rotating 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 so-called 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, σ₀ and T₀ (see Eq. (18)). These two quantities correspond respectively to the mean value and standard deviation of the one-chord travel time distribution. Using the formalism of quantum chaos, we then wrote a semiclassical expression of the autocorrelation and showed that knowing σ₀ and T₀ 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₀. Since chaotic trajectories cover the entire meridional plane, at least in the asymmetric case, we found that T₀ 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 (σ₀/T₀ ≪ 1) implies the observed quasi-degeneracy 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 low-frequency 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 Wigner-Dyson 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 pole-on configuration can be helpful as avoiding m ≠ 0 sub-spectra would simplify mode identification.

Acknowledgments

We thank CALMIP (“CALcul en MIdi-Pyrénées”) for the use of their supercomputer. We used the code Top developed by D. Reese and made user-friendly by B. Putigny. We thank ISSI (“International Space Science Institute”) through the SoFAR (“Seismology of Fast Rotating Stars”) program for their support.

References

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 one-dimensional to avoid carrying indices in the notation. The ray dynamics is governed by the following Hamiltonian:

$$\begin{array}{c}\hfill H=\sqrt{c_{\mathrm{s}}^2{k}^2+{\omega }_c^2}=\omega ,\end{array}$$(A.1)

where the wave vector plays the role of momentum: p ≡ k.

Hamilton’s principal function

We now express p as a function of $\dot{q}$:

$$\begin{array}{c}\hfill \frac{\partial H}{\partial p}=\dot{q}=2c_{\mathrm{s}}^{2}p/\omega \Rightarrow p=\frac{\dot{q}\omega }{2c_{\mathrm{s}}^2}·\end{array}$$(A.2)

The ray system can thus be seen as analogous to a mechanical system with a varying mass $\omega /c_{\mathrm{s}}^2$. Then, the Lagrangian L is obtained through the usual Legendre transform:

$$\begin{array}{c}\hfill L=\dot{q}p-H=\frac{{\dot{q}}^2\omega }{2c_{\mathrm{s}}^2}-\omega .\end{array}$$(A.3)

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:

$$\begin{array}{c}\hfill R(q,q\prime ,\delta t)=\frac{{(q-q\prime )}^2\omega }{2c_{\mathrm{s}}^2\delta t}-\omega \delta t,\end{array}$$(A.4)

using $\dot{q}=(q-q\prime )/\delta t$.

Action integral

The action integral S is defined from R as $S(q,q\prime ,\omega )=R(q,q\prime ,t)+\omega t={\int }_0^t(L+H)\mathrm{d}t$. From Eq. (A.3) we have:

$$\begin{array}{cc}\hfill S(q(t),q\prime (0))& ={{\displaystyle \int }}_0^t\frac{{\dot{q}}^2\omega }{2c_{\mathrm{s}}^2}\mathrm{d}t\hfill \\ \hfill & ={{\displaystyle \int }}_{q\prime (0)}^{q(t)}\frac{{\dot{q}}^2\omega }{2c_{\mathrm{s}}^2}\frac{\mathrm{d}q}{\dot{q}}\hfill \\ \hfill & ={{\displaystyle \int }}_{q\prime (0)}^{q(t)}\frac{\dot{q}\omega }{2c_{\mathrm{s}}^2}\mathrm{d}q={{\displaystyle \int }}_{q\prime (0)}^{q(t)}p\mathrm{d}q,\hfill \end{array}$$(A.5)

which is the well know expression of the action. Finally, the eikonal equation Eq. (2) gives $\stackrel{\sim }{k}=p/\omega =1/{\stackrel{\sim }{c}}_{\mathrm{s}}$. Thus the action can be written in terms of the acoustic time:

$$\begin{array}{c}\hfill S(q,q\prime ,\omega )=\omega {{\displaystyle \int }}_{q\prime }^q\frac{\mathrm{d}s}{{\stackrel{\sim }{c}}_{\mathrm{s}}},\end{array}$$(A.6)

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 four-dimensional phase space. Let Λ⁻¹ 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):

$$\begin{array}{cc}& {({\partial }_t\varphi )}^2-\frac{1}{{\mathrm{\Lambda }}^2}({\omega }_{\mathrm{c}}^2+c_{\mathrm{s}}^2{k}^2)=0,\hfill \end{array}$$(A.7)

$$\begin{array}{cc}& \mathrm{\Lambda }{\partial }_t\varphi \pm {({\omega }_c^2+c_{\mathrm{s}}^2{k}^2)}^{1/2}=0.\hfill \end{array}$$(A.8)

Equation (A.7) can be recognized as the Hamilton-Jacobi 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₁ and P₂ satisfying P₁ + P₂ = I.

Additionally, the substitution $\dot{\varphi }=-H=-\omega$ yields:

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}(A^2)+\frac{c_{\mathrm{s}}^2}{\omega }\frac{\partial }{\partial q}\left(A^2\frac{\partial R}{\partial q}\right)=0.\end{array}$$(A.9)

Introducing the density ρ = A² and velocity v = 1/m (∂R/∂q), with $m=\omega /c_{\mathrm{s}}^2$ (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 $det\left(\frac{\partial q\prime }{\partial q}\right)$.

Having now both the phase and amplitude, we obtain the semiclassical wave function

$$\begin{array}{cc}\hfill {\mathrm{\Psi }}_{\mathrm{sc}}(q,t)=& A_1{\displaystyle \int }\mathrm{d}q\prime {|det\frac{\partial q\prime }{\partial q}|}^{1/2}{e}^{iR(q,q\prime ,t)}{\mathrm{\Psi }}_{\mathrm{sc}}(q\prime ,0)\hfill \\ \hfill & +A_2{\displaystyle \int }\mathrm{d}q\prime {|det\frac{\partial q\prime }{\partial q}|}^{1/2}{e}^{-iR(q,q\prime ,t)}{\mathrm{\Psi }}_{\mathrm{sc}}(q\prime ,0),\hfill \end{array}$$(A.10)

where A₁ and A₂ 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:

$$\begin{array}{cc}\hfill {\mathrm{\Psi }}_{\mathrm{sc}}(q,t)=& A_1{\displaystyle \int }\mathrm{d}q\prime {\displaystyle \sum _j}{|det\frac{\partial q\prime }{\partial q}|}^{1/2}{e}^{i{R}_j(q,q\prime ,t)-i{\kappa }_j\pi /2}{\mathrm{\Psi }}_{\mathrm{sc}}(q\prime ,0)\hfill \\ \hfill & +A_2{\displaystyle \int }\mathrm{d}q\prime {\displaystyle \sum _j}{|det\frac{\partial q\prime }{\partial q}|}^{1/2}{e}^{-i{R}_j(q,q\prime ,t)-i{\kappa }_j\pi /2}{\mathrm{\Psi }}_{\mathrm{sc}}(q\prime ,0),\hfill \end{array}$$(A.11)

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

$$\begin{array}{c}\hfill \psi (q,t)={\displaystyle \int }\mathrm{d}q\prime K(q,q\prime ,t)\psi (q\prime ,0).\end{array}$$(A.12)

The propagator is the time-dependent Green’s function, it is thus solution of the wave equation with the initial condition lim_t → 0K(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)

$$\begin{array}{c}\hfill K_{\mathrm{sc}}(q,q\prime ,\delta t)=A(q,q\prime ,\delta t)e^{iR(q,q\prime ,\delta t)}.\end{array}$$(A.13)

Neglecting the second term ω δt in (A.4) gives

$$\begin{array}{c}\hfill K_{\mathrm{sc}}(q,q\prime ,\delta t)=A(q,q\prime ,\delta t)e^{i\frac{\omega }{2c_{\mathrm{s}}^2{\delta }_t}{(q-q\prime )}^2}.\end{array}$$(A.14)

If we impose $A(q,q\prime ,\delta t)=(\frac{\omega }{2\pi i{c}_{\mathrm{s}}^2\delta t})$, then the previous expression is a 2-dimensional Gaussian of width $\sigma ={(\delta t{c}_{\mathrm{s}}^2/\omega )}^{1/2}$. In the limit δt = 0 the condition K(q″, q′,t = 0) = δ(q − q′) is satisfied. Using again the substitution $m=\omega /c_{\mathrm{s}}^2$ for clarity

$$\begin{array}{c}\hfill K_{\mathrm{sc}}(q,q\prime ,\delta t)=\left(\frac{m}{2\pi i\delta t}\right)e^{iR(q,q\prime ,\delta t)},\end{array}$$(A.15)

using the fact that $m/\delta t=det(-\frac{\partial p\prime }{\partial q})=det(-\frac{{\partial }^{2}R}{\partial q\partial q\prime })$ the expression becomes

$$\begin{array}{c}\hfill K_{\mathrm{sc}}(q,q\prime ,t)={(2\pi i)}^{-1}{|det\frac{{\partial }^{2}R}{\partial q\partial q\prime }|}^{1/2}{e}^{i{R}_j(q,q\prime ,t)}.\end{array}$$(A.16)

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 $det\frac{\partial p\prime }{\partial q}det\frac{\partial q}{\partial q″}=det\frac{\partial p\prime }{\partial q″}$. This gives:

$$\begin{array}{cc}\hfill K_{\mathrm{sc}}(q,q\prime ,t)=& A_1{\displaystyle \sum _j}{(2\pi i)}^{-1}{|det\frac{{\partial }^{2}R}{\partial q\partial q\prime }|}^{1/2}{e}^{i{R}_j(q,q\prime ,t)-i{\kappa }_j\pi /2}\hfill \\ \hfill & +A_2{\displaystyle \sum _j}{(2\pi i)}^{-1}{|det\frac{{\partial }^{2}R}{\partial q\partial q\prime }|}^{1/2}{e}^{-i{R}_j(q,q\prime ,t)-i{\kappa }_j\pi /2},\hfill \end{array}$$(A.17)

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 energy-dependent Green function G(q, q′,E), which is related to the propagator through the Fourier transform

$$\begin{array}{c}\hfill G(q,q\prime ,E)=\frac{1}{iħ}{{\displaystyle \int }}_0^{\infty }K(q,q\prime ,t)e^{\frac{i}{ħ}Et+\frac{ϵt}{ħ}}\mathrm{d}t,\end{array}$$(A.18)

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

$$\begin{array}{c}\hfill G(q,q\prime ,E)={\displaystyle \sum _n}\frac{{\varphi }_n^{\ast }(q){\varphi }_n(q\prime )}{E-E_n+iϵ}·\end{array}$$(A.19)

In the semiclassical approximation this leads to the following equality for E/ℏ → ∞

$$\begin{array}{c}\hfill {\displaystyle \sum _n}\frac{{\varphi }_n^{\ast }(q){\varphi }_n(q\prime )}{E-E_n+iϵ}=\frac{1}{iħ}{{\displaystyle \int }}_0^{\infty }{K}_{\mathrm{sc}}(q,q\prime ,t\prime )e^{\frac{i}{ħ}Et}\mathrm{d}t,\end{array}$$(A.20)

which is valid at first order in ℏ (acoustic case: first order in Λ⁻¹). One can already forecast the trace formula from this equality, as it involves both the eigenenergies of the quantum system on the left-hand side and classical quantities in the right-hand 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. Cohen-Tannoudji et al. 1977):

$$\begin{array}{c}\hfill K(q,q\prime ,t)=⟨q|U(t,0)|q\prime ⟩.\end{array}$$(A.21)

U(t, 0) satisfies the wave equation: ${\partial }^2/\partial t^{2}U(t)=-{\widehat{H}}^{2}U(t,0)$, where $\widehat{H}$ is the Hamiltonian operator with ${\widehat{H}}^2=c_{\mathrm{s}}^2{\mathrm{\nabla }}^2+{\omega }_c^2$. Then:

$$\begin{array}{c}\hfill U(t,0)=P_1{e}^{-i\widehat{H}t}+P_2{e}^{i\widehat{H}t},\end{array}$$(A.22)

with P₁ and P₂ projectors in two subspaces as above. Let |ϕ_n⟩ be the eigenfunctions of the Hamiltonian $\widehat{H}$. Then, from the closure relation, we have:

$$\begin{array}{c}\hfill U(t,0)={\displaystyle \sum _n}|{\varphi }_n⟩⟨{\varphi }_n|(A_1{e}^{-i{\omega }_{n}t}+A_2{e}^{i{\omega }_{n}t}),\end{array}$$(A.23)

where A₁ and A₂ are the projection of the intial wave function on the two subspaces. Using Eq. (A.23) in Eq. (A.21) gives:

$$\begin{array}{c}\hfill K(q,q\prime ,t)={\displaystyle \sum _n}{\varphi }_n^{\ast }(q){\varphi }_n(q\prime )(A_1{e}^{-i{\omega }_{n}t}+A_2{e}^{i{\omega }_{n}t}),\end{array}$$(A.24)

The ω dependent Green’s function G(q″, q′,ω) stems from the Fourier transform of K(q, q′,0):

$$\begin{array}{cc}\hfill G_{\pm }(q,q\prime ,\omega )& ={\displaystyle \sum _n}{\varphi }_n^{\ast }(q){\varphi }_n(q\prime )\frac{1}{i}{{\displaystyle \int }}_0^{\infty }{e}^{-i(\omega \pm {\omega }_n)t+ϵt}\mathrm{d}t\hfill \\ \hfill & =A_1{\displaystyle \sum _n}\frac{{\varphi }_n^{\ast }(q){\varphi }_n(q\prime )}{\omega -{\omega }_n+iϵ}+A_2{\displaystyle \sum _n}\frac{{\varphi }_n^{\ast }(q){\varphi }_n(q\prime )}{\omega +{\omega }_n+iϵ}·\hfill \end{array}$$(A.25)

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

$$\begin{array}{c}\hfill G_{\mathrm{sc}}(q,q\prime ,\omega )=\frac{1}{i}{{\displaystyle \int }}_0^{\infty }{K}_{\mathrm{sc}}(q,q\prime ,t)e^{i\omega t}\mathrm{d}t.\end{array}$$(A.26)

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 $\frac{\partial R(q,q\prime ,t)}{\partial t}+\omega =0$ 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 $\frac{\partial R(q,q\prime ,t)}{\partial t}-\omega =0$. The result is:

$$\begin{array}{cc}\hfill G_{\mathrm{sc}}(q,q\prime ,\omega )=& A_1\frac{1}{i\sqrt{2\pi }}{\displaystyle \sum _j}{{\displaystyle \left|\frac{1}{\dot{q}\dot{q}\prime }\frac{{\partial }^{2}S}{\partial q_{\perp }\partial q_{\perp }\prime }\right|}}^{1/2}{e}^{i{S}_j(q,q\prime ,\omega )-i\kappa {\prime }_j\pi /2}\hfill \\ \hfill & +A_2\frac{1}{i\sqrt{2\pi }}{\displaystyle \sum _j}{{\displaystyle \left|\frac{1}{\dot{q}\dot{q}\prime }\frac{{\partial }^{2}S}{\partial q_{\perp }\partial q_{\perp }\prime }\right|}}^{1/2}{e}^{iS{\prime }_j(q,q\prime ,\omega )-i\kappa {\prime }_j\pi /2},\hfill \end{array}$$(A.27)

where the sum is over all classical trajectories from q′ to q at frequency ω, $\dot{q}$ and $\dot{q}\prime$ are final and initial velocities, and q_⊥ and $q_{\perp }^{\text{'}}$ 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:

$$\begin{array}{c}\hfill \mathrm{Im}\mathrm{Tr}G(q,q\prime ,\omega )=-\frac{1}{\pi }{\displaystyle \sum _j}\delta (\omega -{\omega }_j).\end{array}$$(A.28)

On the other hand, the trace of the semiclassical Green’s function Eq. (A.27) is:

$$\begin{array}{cc}\hfill \mathrm{Tr}{G}_{\mathrm{sc}}& ={\displaystyle \int }\mathrm{d}q{G}_{\mathrm{sc}}(q,q,\omega )\hfill \\ \hfill & ={\displaystyle \int }\mathrm{d}q\frac{1}{i\sqrt{2\pi }}{\displaystyle \sum _j}{{\displaystyle \left|\frac{1}{\dot{q}\dot{q}\prime }\frac{{\partial }^{2}S}{\partial q_{\perp }\partial q_{\perp }\prime }\right|}}^{1/2}{e}^{i{S}_j(q,q\prime ,\omega )-i\kappa {\prime }_j\pi /2},\hfill \end{array}$$(A.29)

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₀ 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 $\frac{\partial S(q,q\prime ,\omega )}{\partial q\prime }+\frac{\partial S(q,q\prime ,\omega )}{\partial q}=0$ 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:

$$\begin{array}{c}\hfill \mathrm{Tr}{G}_{\mathrm{sc}}(q,q\prime ,\omega )=\mathrm{Tr}{G}_0+{\displaystyle \sum _j}\frac{1}{i}{T}_j{|det(M_i-I)|}^{-1/2}{e}^{i{S}_j(\omega )-i\kappa {\prime }_j\pi /2},\end{array}$$(A.30)

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 $d(\omega )=\overline{d}(\omega )+d^{\mathrm{osc}}(\omega )$. 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 $\overline{d}(\omega )$ corresponds to Tr G₀, 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:

$$\begin{array}{c}\hfill d^{\mathrm{osc}}(\omega )=\frac{1}{\pi }\mathit{Re}{\displaystyle \sum _j}{T}_j{|det(M_i-I)|}^{-1/2}{e}^{i{S}_j(\omega )-i\kappa {\prime }_j\pi /2}.\end{array}$$(A.31)

All Tables

All Figures

Fig. 1. PSS at rotation rate Ω/Ωk = 0.589 with ${\stackrel{\sim }{L}}_z=0$, where θ is the colatitude and kθ the projection of the wave vector on the line tangent to the rp(θ) = rs(θ)−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: 6-period island, yellow: 2-period 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 ${\stackrel{\sim }{L}}_z=0.16/{\omega }_{\mathrm{p}}$, where θ is the colatitude and kθ the projection of the wave vector on the line tangent to the rp(θ) = rs(θ)−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) 2-period island mode (ℓ = 0), with a black line indicating the central periodic orbit, and (d) 6-period island mode. The figure shows the scaled pressure amplitude $\mathrm{\Psi }=P\sqrt{d/{\rho }_0}$, 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 lmax, 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 lmax for chaotic modes, marked by a dashed line, and the lowest value of lmax 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 $P(\stackrel{\sim }{s})$, 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