Regularities in the spectrum of chaotic p-modes in rapidly rotating stars

TL;DR

This paper investigates the spectral properties of chaotic p-modes in rapidly rotating stars, revealing organized structures and correlations in their spectra through wave chaos analysis.

Contribution

It introduces a physical interpretation of spectral peaks in chaotic modes and uncovers a pseudo large separation similar to island modes.

Findings

Chaotic mode spectra exhibit a pseudo large separation.

Chaotic modes are organized in series with regular spacing.

Correlations are explained by acoustic ray trapping near stable islands.

Abstract

Interpreting the oscillations of massive and intermediate mass stars remains a challenging task. In fast rotators, the oscillation spectrum of p-modes is a superposition of sub-spectra which correspond to different types of modes, among which island modes and chaotic modes are expected to be the most visible. In the case of island modes, a semi-analytic formula describing the asymptotic behavior of island modes has been obtained previously. We study the properties of high frequency chaotic p-modes in a polytropic model. Unexpected peaks appear in the frequency autocorrelations of the spectra. Our goal is to find a physical interpretation for these peaks and also to provide an overview of the mode properties. We use the 2D oscillation code TOP to produce the modes and acoustic ray simulations to explore the wave properties in the asymptotic regime. Using the tools developed in the field of quantum chaos (or wave chaos), we derive an expression for the frequency autocorrelation involving the travel time of acoustic rays. Chaotic mode spectra were previously thought to be irregular, i. e. described only through their statistical properties. Our analysis shows the existence, in chaotic mode spectra, of a pseudo large separation. This means that chaotic modes are organized in series, such that the modes in each series follow a nearly regular frequency spacing. The pseudo large separation of chaotic modes is very close to the large separation of island modes. Its value is related to the sound speed averaged over the meridional plane of the star. In addition to the pseudo large separation, other correlations appear in the numerically calculated spectra. We explain their origin by the trapping of acoustic rays near the stable islands.

Tables

Table 1. Table 1: Characteristics of chaotic frequency spectra. We show the symmetry classes, which are denoted by m ± superscript 𝑚 plus-or-minus m^{\pm} (where m 𝑚 m is the azimuthal quantum number, + + is even parity and − - is odd parity), the rotation rate, the number of levels and the frequency domain.

symmetry class	Rotation ($\Omega /{\Omega }_k$)	levels	frequencies ($\omega /{\omega }_p$)
$0^{-}$	$0.481$	206	28.35 - 46.89
$0^{-}$	$0.545$	223	28.15 - 44.09
$0^{-}$	$0.589$	217	26.02 - 40.29
$0^{-}$	$0.658$	207	36.37 - 44.89
$0^{-}$	$0.706$	283	23.57 - 36.22
$0^{-}$	$0.809$	170	24.02 - 30.01
$0^{+}$	$0.545$	105	38.01 - 44.06
$0^{+}$	$0.589$	96	30.52 - 36.60
$0^{+}$	$0.658$	120	36.40 - 41.25
$1^{-}$	$0.589$	125	30.51 - 38.48
$4^{-}$	$0.589$	93	30.53 - 38.51

Table 2. Table 2: Comparison of the large separation of island modes Δ i subscript Δ 𝑖 \Delta_{i} and chaotic modes Δ c subscript Δ 𝑐 \Delta_{c} , obtained from the simulated spectra at six rotations

$\Omega /{\Omega }_k$	${\Delta }_i$	${\Delta }_c$	$\left|{\Delta }_i-{\Delta }_c\right|/{\Delta }_i$
$0.481$	$1.1644{\omega }_p$	$1.1458{\omega }_p$	$1.60\%$
$0.545$	$1.1288{\omega }_p$	$1.0907{\omega }_p$	$3.38\%$
$0.589$	$1.1021{\omega }_p$	$1.1132{\omega }_p$	$1.01\%$
$0.658$	$1.0543{\omega }_p$	$1.0543{\omega }_p$	$0.00\%$
$0.706$	$1.0155{\omega }_p$	$1.0169{\omega }_p$	$0.14\%$
$0.809$	$0.9345{\omega }_p$	$0.9355{\omega }_p$	$0.10\%$

Equations

$$\Delta \Psi +\frac{{\omega }^2-{\omega }_c^2}{c_s^2}\Psi =0,$$

$${\omega }^2=c_s^2{k}^2+{\omega }_c^2.$$

$$H^{\prime }$$

$$W$$

$$\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}={\tilde{\mathbf{k}}}_{\mathbf{p}},\frac{\mathrm{d}{\tilde{\mathbf{k}}}_{\mathbf{p}}}{\mathrm{d}t}=-\nabla W,$$

$$\Psi (r,\theta ,\varphi )=\sum _{l=0}^{N_l}\left(\sum _{j=0}^{N_j}{a}_j^{l,m}{T}_j(2r-1)\right),Y_l^m(\theta ,\varphi ),$$

$$\mathrm{d}(\omega )-\bar{\mathrm{d}}(\omega )=\mathrm{Re}\sum _j{A}_j{e}^{\mathrm{i}{S}_j(\omega )},$$

$$A_j=\frac{1}{\pi }\frac{T_j}{|\det (M_j-I){|}^{1/2}},$$

$$\rho (T)\approx (1/T)e^{\lambda T},$$

$$A(T)\approx \frac{1}{\pi }T{e}^{-(1/2)\lambda T}.$$

$$K(T)=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{\infty }\mathrm{d}\xi \exp (\mathrm{i}\xi T)\mathrm{C}(\xi ),$$

$$C(\xi )$$

$$\displaystyle=\Big{\langle}\Big{(}\mathrm{Re}\sum_{i}A_{i}e^{\mathrm{i}S_{i}(\omega-1/2\xi)}\Big{)}\Big{(}\mathrm{Re}\sum_{j}A_{j}e^{\mathrm{i}S_{j}(\omega+1/2\xi)}\Big{)}\Big{\rangle}.$$

$$C(\xi)=\frac{1}{2}\Big{\langle}\Big{(}\sum_{i}A_{i}e^{\mathrm{i}S_{i}(\omega-1/2\xi)}\Big{)}\Big{(}\sum_{j}A_{j}e^{\mathrm{i}S_{j}(\omega+1/2\xi)}\Big{)}^{*}\Big{\rangle}.$$

$$K(T)$$

$$\displaystyle=\frac{1}{2}\left\langle\sum_{i}\sum_{j}A_{i}A_{j}^{*}\,e^{i(S_{i}-S_{j})}\,\frac{2\pi}{\sqrt{2\pi}}\delta\Big{(}T-\frac{1}{2}(T_{i}+T_{j})\Big{)}\right\rangle.$$

$$\mathrm{K}(T)\approx \sum _j{A}_j^2\delta (T-T_j).$$

$$K(T)\approx A^2(T)\rho (T),$$

$${\mathrm{P}}_{\mathrm{n},\mathrm{\Omega }}(T)=\frac{T_0}{\sqrt{2\pi n}{\sigma }_0}\exp \left(-\frac{(T-n{T}_0{)}^2}{2(\sqrt{n}{\sigma }_0{)}^2}\right),$$

$$F(\xi )=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{\infty }\mathrm{d}T\exp (-i\xi T)T{\mathrm{P}}_{\mathrm{\Omega }}(T),$$

$$T^{\mathrm{av}}=2\left(\frac{2}{\pi }{\int }_0^{\pi /2}\tau (\theta )\mathrm{d}\theta \right),$$

$$H=\sqrt{c_s^2{k}^2+{\omega }_c^2}=\omega ,$$

$$\frac{\partial H}{\partial p}=\dot{q}=2c_s^{2}p/\omega ⟹p=\frac{\dot{q}\omega }{2c_s^2}.$$

$$L=\dot{q}p-H=\frac{{\dot{q}}^2\omega }{2c_s^2}-\omega .$$

$$R(q,q^{\prime },\delta t)=\frac{(q-q^{\prime }{)}^2\omega }{2c_s^2\delta t}-\omega \delta t,$$

$$S(q(t),q^{\prime }(0))$$

$$={\int }_{q^{\prime }(0)}^{q(t)}\frac{{\dot{q}}^2\omega }{2c_s^2}\frac{\mathrm{d}q}{\dot{q}}$$

$$={\int }_{q^{\prime }(0)}^{q(t)}\frac{\dot{q}\omega }{2c_s^2}\mathrm{d}q={\int }_{q^{\prime }(0)}^{q(t)}p\mathrm{d}q,$$

$$S(q,q^{\prime },\omega )=\omega {\int }_{q^{\prime }}^q\frac{\mathrm{d}s}{{\tilde{c}}_s},$$

$$({\partial }_t\varphi {)}^2-\frac{1}{{\Lambda }^2}({\omega }_c^2+c_s^2{k}^2)=0,$$

$$\Lambda {\partial }_t\varphi \pm ({\omega }_c^2+c_s^2{k}^2{)}^{1/2}=0.$$

$$\frac{\partial }{\partial t}(A^2)+\frac{c_s^2}{\omega }\frac{\partial }{\partial q}\left(A^2\frac{\partial R}{\partial q}\right)=0.$$

$${\Psi }_{sc}(q,t)$$

Full text

*Aims. *We study the properties of high frequency chaotic p-modes in a polytropic model. Unexpected peaks appear in the frequency autocorrelations of the spectra. Our goal is to find a physical interpretation for these peaks and also to provide an overview of the mode properties.

*Methods. *We used the 2D oscillation code ”TOP” to produce the modes and acoustic ray simulations to explore the wave properties in the asymptotic regime. Using the tools developed in the field of quantum chaos (or wave chaos), we derived an expression for the frequency autocorrelation involving the travel time of acoustic rays.

*Results. *Chaotic mode spectra were previously thought to be irregular, that is, described only through their statistical properties. Our analysis shows the existence, in chaotic mode spectra, of a pseudo large separation. This means that chaotic modes are organized in series, such that the modes in each series follow a nearly regular frequency spacing. The pseudo large separation of chaotic modes is very close to the large separation of island modes. Its value is related to the sound speed averaged over the meridional plane of the star. In addition to the pseudo large separation, other correlations appear in the numerically calculated spectra. We explain their origin by the trapping of acoustic rays near the stable islands.

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 $\delta$ Scuti stars (García Hernández et al., 2009; García Hernández et al., 2013; García Hernández et al., 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; Ouazzani et al., 2015; Reese et al., 2017). 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 $\delta$ 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 $\delta$ 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 $\sim 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 (Section 3), and develop the semiclassical theory which explains the features observed (Section 4). The Appendix presents detailed derivations needed for Section 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) :

[TABLE]

where $\Psi=\hat{P}/f$, with $\hat{P}$ the complex amplitude associated with the time-harmonic pressure perturbation $P=\mathrm{Re}[\hat{P}\,\exp(-i\omega t)]$, $f$ a function of the background model, $\omega_{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

[TABLE]

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_{s}^{2}k^{2}+\omega_{c}^{2})^{1/2}$. The motion takes place in a meridional plane rotating with the ray at an angular velocity $\mathrm{d}\phi/\mathrm{d}t=\tilde{L}_{z}/\mathrm{d}^{2}$, where $\tilde{L}_{z}=r\sin\theta\,k_{\phi}/\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 :

[TABLE]

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

[TABLE]

where $\mathrm{d}t=c_{s}\mathrm{d}s/\left(1-\frac{\omega_{c}^{2}}{\omega^{2}}\right)^{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 $\boldsymbol{\tilde{k}_{p}}$ as simply $\boldsymbol{\tilde{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}(\theta)$ : $r_{p}(\theta)=r_{s}(\theta)-z$, where $r_{p}(\theta)$ 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).

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 $\Omega$ 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 $\Omega=0$ system. This is represented in Fig. 1. We did not draw the small island structures in the domain where $|k_{\theta}/\omega|$ 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_{k}=(GM/R_{eq}^{3})^{1/2}$, where $G$ is Newton’s constant, $M$ is the mass of the star and $R_{eq}$ the equatorial radius. Near $\Omega/\Omega_{k}\simeq 0.32$ the chaotic zone is comparable in size to the main stable island and continues to grow afterward (see Fig. 2). The ratio $\mathcal{V}_{\mathrm{I}}/\mathcal{V}_{\mathrm{II}}$, where $\mathcal{V}_{\mathrm{I}}$ is the volume occupied by the main stable island and $\mathcal{V}_{\mathrm{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 $\mathcal{A}_{\mathrm{I}}$ and $\mathcal{A}_{\mathrm{II}}$ on the PSS, corresponding respectively to the main island and chaotic zones, we found out that the function $\mathcal{A}_{\mathrm{I}}(\Omega)/\mathcal{A}_{\mathrm{II}}(\Omega)$ is not monotonic. In fact, the ratio of the 2-period island over the chaotic zone shrinks beyond $\Omega/\Omega_{k}=0.658$, reaches a minimum around $\Omega/\Omega_{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 $\Omega/\Omega_{k}=0.682$ were not considered.

The non axisymmetric ray dynamics $\tilde{L}_{z}\neq 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(|\tilde{L}_{z}|/\tilde{L})<\theta<\arcsin(|\tilde{L}_{z}|/\tilde{L})$, $\tilde{L}$ being the frequency-scaled angular momentum norm.

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 $\Gamma=5/3$. We impose that the pressure and density satisfy a polytropic relation $P\propto\rho^{1+1/\mu}$, with $\mu=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}\leq 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 :

[TABLE]

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_{\theta}$ and $k_{r}$. Using the PSS of a given stellar model, we can determine the upper limit of the wave vector components $\tilde{k}_{\theta}^{max}$ or $\tilde{k}_{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 $\Omega_{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 : $\Omega/\Omega_{k}=0.481,0.545,0.589,0.658$, $0.706$ and $0.809$ with frequencies ranging from $\sim 23\,\omega_{p}$ to $\sim 47\,\omega_{p}$, with $\omega_{p}=(GM/R_{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.

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 $\omega_{n,l,m}=n\delta_{n}(m)+\ell\delta_{\ell}(m)+\alpha$, where $n$ (resp. $\ell$) is the number of nodes along (resp. transverse to) the central orbit $\gamma$ of the island, shown in Fig. 4, and $\alpha$ is a constant (Pasek et al., 2012). Thus the frequency spectrum is completely determined by the regular spacings $\delta_{n}$ and $\delta_{\ell}$, with the half large separation $\delta_{n}(m)=2\pi/\oint_{\gamma}(\mathrm{d}s/\tilde{c}_{s})$ 111The large separation is expressed as $\Delta_{i}=2\pi/\int_{\gamma}(\mathrm{d}s/\tilde{c}_{s})$, where the integral is calculated along the path $\gamma$ between two points on the surface, which is half the full period of the orbit.. 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 $\tilde{k}_{\theta}^{max}$ of any whispering gallery ray is higher than the $\tilde{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 $\bar{a}^{l}\equiv\tfrac{1}{N_{j}}\sum_{j=0}^{N_{j}}|a_{j}^{l}|$. Thus the dominant degree $l_{max}$ corresponds to the value of $l$ where $\bar{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}\leq 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.

At finite frequency, here $\omega\leq 47\omega_{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, $\ell=4$ island modes cannot be distinguished from interface modes without ambiguity from the amplitudes.

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(\tilde{s})$, with $\tilde{s}_{n}=(\omega_{n+1}-\omega_{n})/\langle\omega_{n+1}-\omega_{n}\rangle$, where $\langle\omega_{n+1}-\omega_{n}\rangle$ corresponds to the mean frequency spacing, is expected to follow the Wigner surmise of Random Matrix Theory $P(\tilde{s})=\frac{\pi\tilde{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(\tilde{s})\,\mathrm{d}\tilde{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 $\Omega/\Omega_{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 Section 4.2).

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.

3.2 Regularities in the spectra

We computed the autocorrelation $R_{2}(\xi)=\int d(\omega-1/2\xi)d(\omega+1/2\xi)\,\mathrm{d}\omega$ of a variety of chaotic spectra, where the density $d(\omega)$ 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 $\Delta_{c}$, slightly decreases with increasing rotation. Moreover, additional peaks of significant amplitude appear in the autocorrelations, their relative amplitude being large especially at $\Omega/\Omega_{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.

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 $\Delta_{c}$. To verify this idea, we make use of so-called échelle diagrams. On such diagrams, two modes distant in frequency of $\Delta$ 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 $\Omega/\Omega_{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\,\omega_{p}$ and $33.78\,\omega_{p}$ are represented on an échelle diagram, using a folding value equal to the main peak $\Delta_{c}=1.0899\,\omega_{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 $\Delta_{c}$.

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 Section 4.

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 $\Delta_{c}$ 222the systematic search for consecutive modes was performed only for the chaotic spectrum at $\Omega/\Omega_{k}=0.589$ represented in the échelle diagram., 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.

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, $\Delta_{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 $\Delta_{i}=1.0907\omega_{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 $\ell$, 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 $\Delta_{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 $\ell$ 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 $\omega=27\Delta_{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 $\ell$ ; for instance, series 1 is followed by $\ell=4$ island modes and series 2 by $\ell=5$ island modes. This transition is visible by comparing Fig. 13, where series 1 ends just below $27\Delta_{c}$, and Fig. 11 where the $\ell=4$ island modes start just above $27\Delta_{c}$.

3.2.2 Secondary peaks

In addition to the main peak discussed above, we see other peaks that we call secondary peaks. At $\Omega/\Omega_{k}=0.589$, $\Omega/\Omega_{k}=0.706$ and $\Omega/\Omega_{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 $\alpha\delta+\beta\Delta_{c}$, where $\delta$ is the position of the first secondary peak and $\alpha$ and $\beta$ are integers. For instance, the peaks of Fig. 9, panel c, occur at $(\alpha=+1,\,\beta=0)$ and $(\alpha=+1,\,\beta=+1)$. At $\Omega/\Omega_{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 $\Omega/\Omega_{k}=0.706$. We will propose an explanation for these secondary peaks in Section 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 $\tilde{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 $\omega=24.41\omega_{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).

4 Semiclassical interpretation

In section 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.

Using the semiclassical periodic orbit theory, which relates the mode properties to the acoustic ray dynamics, we show in subsection 4.1 that the $\Delta_{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 subsection 4.2. In subsection 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 subsection 4.4 we discuss how the families of chaotic modes evolve when the frequency domain changes. At last in subsection 4.5, we show that the proximity of $\Delta_{c}$ and $\Delta_{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(\omega)$ counts the number of modes below a given frequency $\omega$. From the staircase function one can define the spectral density $\mathrm{d}(\omega)=\mathrm{d}N/\mathrm{d}\omega$. Both the staircase function and the spectral density are often written as the sum of two contributions called the smooth part $\bar{N}(\omega)$ (or $\bar{\mathrm{d}}(\omega)$) and the fluctuating (or oscillating) part $N^{osc}(\omega)$ (or $\mathrm{d}^{osc}(\omega)$), 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 $\gamma_{j}$, the formula reads :

[TABLE]

where $S_{j}=\omega\oint_{\gamma_{j}}\mathrm{d}s/\tilde{c}_{s}$ is the action and $\oint_{\gamma_{j}}\mathrm{d}s/\tilde{c}_{s}=T_{j}$ is the acoustic travel time of $\gamma_{j}$. The amplitude term is expressed as

[TABLE]

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 $\rho(T)$ of orbits with a time period $T$ grows exponentially as

[TABLE]

where $\lambda$ 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

[TABLE]

The exponential growth of the number of long orbits makes it numerically difficult to attain good precision on the spectral density using the trace formula. Moreover in many systems, including the stellar model studied here, there is no automatic procedure to find the periodic orbits. To circumvent these issues, Berry showed (Berry, 1985) that some statistical quantities such as the variance, or more importantly the autocorrelation of the spectrum, can be approximated using the trace formula in systems where individual periodic orbits are not known, provided one knows their distribution.

The autocorrelation $R_{2}(\xi)=\left\langle\mathrm{d}(\omega-1/2\,\xi)\,\mathrm{d}(\omega+1/2\,\xi)\right\rangle$, with the average $\langle f\rangle=\int f\mathrm{d}\omega$, can be re-written as $R_{2}(\xi)=\left<\mathrm{d^{osc}}(\omega-1/2\xi)\,\mathrm{d^{osc}}(\omega+1/2\xi)\right\rangle+\phi(\xi)$, where $\phi(\xi)$ is a smooth function. The theory aims at deriving an approximate expression of the form factor :

[TABLE]

which is the Fourier transform of the autocorrelation

[TABLE]

Getting rid of the real parts through the fact that with $x=\mathrm{Re}(\tilde{x})=A_{i}\cos(\xi T_{i}+S_{i})$ and $y=\mathrm{Re}(\tilde{y})=A_{j}\cos(\xi T_{j}+S_{j})$ then $\langle x\,y\rangle=1/2\,\langle\tilde{x}\,\tilde{y}^{*}\rangle$, gives

[TABLE]

Inserting this expression in Eq. 11 and expanding the action as $S_{j}(\omega\pm\,1/2\xi)\approx S_{j}(\omega)\pm 1/2\,\xi\,(\partial S_{j}/\partial\omega)=S_{j}(\omega)\pm 1/2\,\xi T_{j}(\omega)$, where $T_{j}$ is the acoustic travel time introduced before, leads to a new expression of the form factor :

[TABLE]

Under the frequency average $\langle.\rangle$, the contribution of the off-diagonal terms $i\neq j$ in the double sum can be neglected. This diagonal approximation is valid for ”short times” below the Ehrenfest time $T_{E}\approx(1/\lambda)\ln(\omega)$ (Bogomolny et al., 1997). Beyond $T_{E}$, there are pairs of orbits with very close actions $S_{i}\approx 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

[TABLE]

In this expression $K(T)$, as a Fourier transform of a correlation function in frequency, is a function of time $T$. Eq. 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

[TABLE]

where $\rho(T)\,\mathrm{d}T$ is the the number of periodic orbit with a time period between $T$ and $T+\mathrm{d}T$. A prescription for the amplitude $A(T)$ and the density $\rho(T)$ is given in Eq. 10 and Eq. 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)\propto 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 $\rho(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 $\Omega/\Omega_{k}=0.481$. The distribution is a narrow packet : its standard deviation $\sigma_{0}$ is small compared to its mean value $T_{0}$. We interpret this effect as a consequence of the strong decrease of the sound speed near the surface 333To test this hypothesis, we computed the acoustic time distribution in a domain shaped like a $\Omega/\Omega_{k}=0.706$ rotating star, but with a homogeneous sound speed throughout the interior. It results in a distribution whose dispersion is of the same order as the mean value : $\sigma_{0}\approx 0.32\,T_{0}$.. Indeed, the trajectories travel rapidly through the core of the star and the acoustic time is dominated by the surface behavior, so the actual length of the trajectory has very little impact on its travel time. We expect the characteristic time $T_{0}$ to increase with rotation, because of the increasing volume of the star. This is indeed confirmed, with $T_{0}=5.19/\omega_{p}$ at $\Omega/\Omega_{k}=0.481$ and $T_{0}=6.30/\omega_{p}$ at $\Omega/\Omega_{k}=0.706$.

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 $\leq 20$, with the same number of orbits for all values of n, is represented in Fig. 16, panels (a) and (b) at rotations $\Omega/\Omega_{k}=0.481$ and $\Omega/\Omega_{k}=0.809$. The packets are evenly spaced out ($T_{n}\approx nT_{0}$) 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 $\sigma_{0}/T_{0}$ which is $0.049$ at $\Omega/\Omega_{k}=0.481$ and $0.097$ at $\Omega/\Omega_{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 444Looking closely, one can discern a very narrow peak inside each packet. It is created by trajectories trapped near the main island (see subsection 4.2) and will have no significant impact on the autocorrelation. 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 : $\mathrm{P_{\,\Omega}}(T)=\sum_{n}\mathrm{P_{n,\,\Omega}}(T)$, with $\mathrm{P_{n,\,\Omega}}(T)$ the probability distribution of n-chords travel times :

[TABLE]

where $T_{0}$, in the prefactor, has been added for normalization purposes. It is a crude approximation but it encapsulates in a simple form all the relevant properties of the distributions, namely the presence of regularly spaced packets and their widening with increasing $n$. One has to keep in mind that the dependency on $\Omega$ is not explicit but hidden in the values of $T_{0}$ and $\sigma_{0}$. Finally, we need to take into account the fact that the number of trajectories of acoustic time $T$ grows with $T$. Thus, we write $\rho(T)=(1/T)e^{\lambda T}\times\mathrm{P_{\,\Omega}}(T)$ and the form factor becomes $K(T)\propto A^{2}(T)\;\rho(T)=T\,\mathrm{P_{\,\Omega}}(T)$.

The Fourier transform of this quantity

[TABLE]

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(\xi)$ (see Eq. 12) and is closely related to the autocorrelation $R_{2}(\xi)$ (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 $\Delta_{c}^{\mathrm{th}}$ moves with rotation, as in the numerical mode computations, from $1.1458\,\omega_{p}$ at $\Omega/\Omega_{k}=0.481$ to $0.9254\,\omega_{p}$ at $\Omega/\Omega_{k}=0.809$. In the top panel, $\Delta_{c}^{\mathrm{th}}$ is compared to the peak position of the numerical modes $\Delta_{c}$, with good agreement. The position, height and width of the peaks are completely determined by $T_{0}$ and $\sigma_{0}$. The position is found straightforwardly as $\Delta_{c}^{\mathrm{th}}\approx 2\pi/T_{0}$. The height and width are controlled by the ratio $\sigma_{0}/T_{0}$. Indeed, the condition $\sigma_{0}\ll T_{0}$ is necessary to clearly distinguish the packets in the distribution of travel times, as in Fig. 16. Due to the increase of $\sigma_{0}/T_{0}$ with rotation, the peak gets less visible at high rotation rates.

The theoretical large separation $\Delta_{c}^{\mathrm{th}}$ is not only a function of the rotation rate $\Omega$, but also of the projected angular momentum $\tilde{L}_{z}$. In the bottom panel of Fig. 17 we show the predicted peaks for $\tilde{L}_{z}=0,0.0819/\omega_{p},0.1638/\omega_{p}$ and $0.2458/\omega_{p}$ at $\Omega/\Omega_{k}=0.589$. The shift in position is governed by the value of $T_{0}$ and the change in amplitude by the value of $\sigma_{0}$ as $\tilde{L}_{z}$ increases. The main difference with the axisymmetric case is a reduction of the domain of propagation, resulting in a smaller standard deviation $\sigma_{0}$. As the mean travel time $T_{0}$ between two rebounds depends on the impact of $\tilde{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 $\Delta_{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_{s}(\theta)}\mathrm{d}r/\tilde{c}_{s}$, where $r$ is the radial coordinate, the mean acoustic time between two rebounds can be estimated without using ray dynamics, as

[TABLE]

where the integral is performed over a quarter of the star because of the axial and equatorial symmetries. In the $m\neq 0$ case, the domain of integration must be adapted, since the acoustic ray cavity is reduced in size. We find that $2\pi/T^{\mathrm{av}}$ is a good estimation of $\Delta_{c}$. For instance, at $\omega/\Omega_{k}=0.589$, $2\pi/T^{\mathrm{av}}=1.0996\omega_{p}$, which differs from $\Delta_{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 $\Omega/\Omega_{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 $\delta_{\ell}=\omega_{n,\ell+1}-\omega_{n,\ell}$ 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 $\Omega/\Omega_{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 $\delta$ 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 $\Omega/\Omega_{k}=0.706$ indicates that it is divided into independent subspectra, as explained in section 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 $\Omega/\Omega_{k}=0.589$. At $\Omega/\Omega_{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 $\Omega/\Omega_{k}=0.589$ to $\Omega/\Omega_{k}=0.706$ and that, in parallel, it takes longer to go through the partial barrier at $\Omega/\Omega_{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 $\Delta_{c}/3$. As seen in Fig. 9, this idea is consistent with the data at $\omega/\omega_{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.

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 $\nabla^{2}\psi+k^{2}\psi=0$. The waves are of the form $\psi_{n}=a_{n}\,\cos(\boldsymbol{k_{n}}(\alpha_{n})\cdot\boldsymbol{x}+\xi_{n})$, where the amplitude $a_{n}$, the wave vector orientation $\alpha_{n}$ and the phase shift $\xi_{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}\gg k_{\theta}$). 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\rightarrow\infty$. 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.

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 $\Delta_{c}$ and having similar amplitude distribution. Sixteen series have been effectively identified at the rotation rate $\Omega/\Omega_{k}=0.589$ and in the frequency range, $25.60\,\omega_{p}<\omega<33.54\,\omega_{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=\pm(\ell_{s}+1/2)$ (Gough, 1993) where $\ell_{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 $\tilde{L}=L/\omega=\pm(\ell_{s}+1/2)/\omega_{n_{s},\ell_{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 $\tilde{k}_{\theta}=\pm\tilde{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 $\ell_{s}$ and $n_{s}$. For fixed radial order $n_{s}$ and increasing $\ell_{s}$, the tori indeed move toward higher $\tilde{L}$ or $\tilde{k}_{\theta}$ values. Similarly, for fixed $\ell_{s}$ and increasing $n_{s}$, the tori approach the $\tilde{k}_{\theta}=0$ axis. In particular, the $\ell_{s}=0$ or $\ell_{s}=1$ modes are already close to the central axis for the smallest radial order.

When rotation comes into play, the location of the tori at $\Omega=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 $\tilde{k}_{\theta}$, the large chaotic zone at low and intermediate $\tilde{k}_{\theta}$, and at large $\tilde{k}_{\theta}$, the region of surviving KAM tori corresponding to whispering gallery trajectories. In this context, high-$\tilde{L}$ tori at $\Omega=0$ will transform into structures of the surviving KAM tori region, whereas low-$\tilde{L}$ tori will be destroyed as the 2-period island chain and the chaotic zone surrounding it emerge. We thus expect high-$\tilde{L}$ modes to become whispering gallery modes, while low-$\tilde{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 $\ell_{s}$ and variable radial orders $n_{s}$. In such a series, $\tilde{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-$\tilde{L}$) modes become whispering gallery modes, the intermediate-$n_{s}$ modes become chaotic modes, and the high-$n_{s}$ (low-$\tilde{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 $\ell_{s}$ at zero rotation.

At $\Omega/\Omega_{k}=0.589$, we could indeed relate two chaotic series (the series $1$ and $2$) with the $\ell=4$ and the $\ell=5$ odd-parity island mode series, respectively. Following this interpretation and using the formulas that link $\ell$ to $\ell_{s}$ (Reese et al., 2008; Pasek et al., 2012) we can attribute the $\ell_{s}=9$ value to series $1$ and the $\ell_{s}=11$ value to series $2$. In principle, the $\ell_{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 $\ell_{s}$ value at low frequency.

The generic case of a $\ell_{s}$ series containing the three type of modes only holds for high enough $\ell_{s}$. The computed PSS indeed show that below a specific $\tilde{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 $\ell_{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 $\ell_{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 $\tilde{L}$ are all close to the $\tilde{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 $\ell_{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 $\Delta_{c}$ organization of the chaotic spectrum, the chaotic series would not exist and we could not attribute $\ell_{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 $\Delta_{i}$ is very close to $\Delta_{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 $\Delta_{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 $\gamma$ by (Pasek et al., 2012) $\Delta_{i}=2\delta_{n}$ with $\delta_{n}=2\pi/\oint_{\gamma}(\mathrm{d}s/\tilde{c}_{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 $\gamma$. For the same reason that all chaotic trajectories have almost the same travel time, the acoustic time between two rebounds along the central orbit of the island has to be very close to $T_{0}$.

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 $\Omega/\Omega_{k}=0.48$ and $\Omega/\Omega_{k}=0.81$. Following the methodology of section 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, $\sigma_{0}$ and $T_{0}$ (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 $\sigma_{0}$ and $T_{0}$ is enough to recover the position of the main peak in the numerically computed spectra. This asymptotic analysis also explains the decrease of the peak height as rotation increases, equivalent to a loss of regularity of chaotic modes. There are other peaks in the autocorrelations, which vary in a less predictable way as a function of rotation. We propose that they are created by the presence of phase space structures that develop around the stable island chains, called partial barriers.

The large separation is expressed as $\Delta_{c}\sim 2\pi/T_{0}$. Since chaotic trajectories cover the entire meridional plane, at least in the asymmetric case, we found that $T_{0}$ can be estimated approximately, without using ray tracing, by computing the average acoustic time over the meridional plane (see Eq. 19).

As explained in section 4.5, the small variance of acoustic travel times ($\sigma_{0}/T_{0}\ll 1$) implies the observed quasi-degeneracy of $\Delta_{i}$ and $\Delta_{c}$. Thus, we expect rapidly rotating stars to be characterized by a unique large separation $\Delta\sim\Delta_{i}\sim\Delta_{c}$. Autocorrelation peaks at the large separation detected in $\delta$ 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\neq 0$ sub-spectra would simplify mode identification.

Acknowledgements.

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.

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

A.1 The Hamiltonian system

We denote the canonically conjugate variables $(\boldsymbol{p},\boldsymbol{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 :

[TABLE]

where the wave vector plays the role of momentum : $p\equiv k$.

A.1.1 Hamilton’s principal function

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

[TABLE]

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

[TABLE]

Hamilton’s principal function, denoted $R$, is defined as the time integral of the Lagrangian. Its computation involves following a trajectory from $(q^{\prime},t^{\prime})$ to $(q,t)$. For short times $\delta t$, it gives :

[TABLE]

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

A.1.2 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. 22 we have :

[TABLE]

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

[TABLE]

where $s$ is the curvilinear coordinate along the ray path.

A.2 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 $\Lambda^{-1}$ be a small dimensionless parameter. In the wave equation, Eq. 1, we insert the WKB ansatz $\Psi(q,t)=A\,e^{i\Lambda\phi(q,t)}$ leading to (Gough, 1993) :

[TABLE]

Eq. 26 can be recognized as the Hamilton-Jacobi equation $\partial R/\partial t=\pm H$. Thus, in the limit of small wavelengths, the phase is simply Hamilton’s principal function : $\phi=R$ or $\phi=-R$ (associated to the Hamiltonian $-H$). The two possible phases lead to two terms in the propagator, with a projector on each subspace $P_{1}$ and $P_{2}$ satisfying $P_{1}+P_{2}=I$.

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

[TABLE]

Introducing the density $\rho=A^{2}$ and velocity $v=1/m\,(\partial R/\partial q)$, with $m=\omega/c_{s}^{2}$ (from the mechanical analogy introduced in section A.1.1), Eq. 28 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^{\prime}$ 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

[TABLE]

where $A_{1}$ and $A_{2}$ are the projection of the intial wave function on the two subspaces. This is valid for short time, i.e. when only one classical trajectory connects $q^{\prime}$ to $q$ in time $t$. For longer times, several trajectories which we label by $j$ connect the two points, and the formula becomes:

[TABLE]

where the topological index $\kappa$ is added to account for the phase shift at points where the amplitude becomes singular, such as caustics.

The propagator $K$ is defined by

[TABLE]

The propagator is the time-dependent Green’s function, it is thus solution of the wave equation with the initial condition $\mbox{lim}_{t\rightarrow 0}K(q,q^{\prime},t)=\delta(q-q^{\prime})$. Again, we will assume that, for short times $\delta t$, the semiclassical propagator is of the form (for the first term, the computation is similar for the second term)

[TABLE]

Neglecting the second term $\omega\;\delta t$ in 23 gives

[TABLE]

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

[TABLE]

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

[TABLE]

This short time expression with the correct limit at $t\rightarrow 0$ can be extended to longer times using Eq. 30 and the fact that $\det\frac{\partial p^{\prime}}{\partial q}\det\frac{\partial q}{\partial q^{\prime\prime}}=\det\frac{\partial p^{\prime}}{\partial q^{\prime\prime}}$. This gives:

[TABLE]

where the sum is over all classical trajectories labelled by $j$ from $q^{\prime}$ to $q$ in time $t$.

A.3 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^{\prime},E)$, which is related to the propagator through the Fourier transform

[TABLE]

where $\epsilon$ is a small positive number which makes the integral convergent and goes to zero eventually. Then one uses the fact that $G(q,q^{\prime},E)$ can be expanded on a basis of eigenvector $\{\phi_{j}\}$ of the Hamiltonian with eigenvalues $E_{n}$ as

[TABLE]

In the semiclassical approximation this leads to the following equality for $E/\hbar\rightarrow\infty$

[TABLE]

which is valid at first order in $\hbar$ (acoustic case : first order in $\Lambda^{-1}$). 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. 39 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 $\ket{q}$ and $\ket{q^{\prime}}$, (see e. g. (Cohen-Tannoudji et al., 1977)) :

[TABLE]

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

[TABLE]

with $P_{1}$ and $P_{2}$ projectors in two subspaces as above. Let ${\ket{\phi_{n}}}$ be the eigenfunctions of the Hamiltonian $\hat{H}$. Then, from the closure relation, we have :

[TABLE]

where $A_{1}$ and $A_{2}$ are the projection of the intial wave function on the two subspaces. Using Eq. 42 in Eq. 40 gives :

[TABLE]

The $\omega$ dependent Green’s function $G(q^{\prime\prime},q^{\prime},\omega)$ stems from the Fourier transform of $K(q,q^{\prime},0)$:

[TABLE]

On the other hand, the semiclassical Green’s function is obtained by taking the Fourier transform of the semiclassical propagator 36 and evaluating it by stationary phase; one starts from

[TABLE]

The phase term of the integrand is the action $R(q,q^{\prime},t)+\omega t=S(q,q^{\prime},\omega)$ or $-R(q,q^{\prime},t)+\omega t=S^{\prime}(q,q^{\prime},\omega)$. 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^{\prime}$ to $q$ at frequency $\omega$. As usual one expands the integrand in Eq. 45 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 36 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:

[TABLE]

where the sum is over all classical trajectories from $q^{\prime}$ to $q$ at frequency $\omega$, $\dot{q}$ and $\dot{q}^{\prime}$ are final and initial velocities, and $q_{\perp}$ and $q_{\perp}^{\prime}$ are coordinates transverse to the orbit. The index $\kappa^{\prime}_{j}$ counts again the singularities along the orbits.

A.4 The final formula

To obtain the trace formula, we should compute the trace of the Green’s function from the two formulas we obtained, Eq. Eq. 44 and Eq. 46. 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. 44. To this aim, let the small imaginary part in the denominator go to zero and get for the imaginary part:

[TABLE]

On the other hand, the trace of the semiclassical Green’s function Eq. 46 is:

[TABLE]

This formula is an integral involving all closed classical paths from $q$ to $q$ at frequency $\omega$. It contains two parts. The first one corresponds to the limit for $q\rightarrow q^{\prime}$ of the short direct trajectories between $q$ and $q^{\prime}$, which become of zero length. We call it $\mbox{Tr}\;G_{0}$ and it should be treated separately. The remaining contains a sum of closed orbits between $q$ and $q$ with nonzero length. We evaluate this sum again by stationary phase; the stationary points in the sum are such that the first derivative of the function in the exponential is cancelled. This implies that $\frac{\partial S(q,q^{\prime},\omega)}{\partial q^{\prime}}+\frac{\partial S(q,q^{\prime},\omega)}{\partial q}=0$ for $q=q^{\prime}$. This selects closed trajectories with equal initial and final momentum, thus periodic orbits. Again, the prefactors are to be combined correctly, yielding to:

[TABLE]

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

Eq. 47 connects the trace of the Green’s function to the density of states $\sum_{j}\delta(\omega-\omega_{j})$. It is known that this sum can be split in two parts $d(\omega)=\bar{d}(\omega)+d^{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 $\bar{d}(\omega)$ corresponds to $\mbox{Tr}\;G_{0}$, while the oscillatory part corresponds to the remaining part of Eq. 49. Putting together Eq. 47 for positive frequencies and Eq. 49 gives the Gutzwiller formula for the acoustic waves:

[TABLE]