Correlations in the chaotic spectrum of pressure modes in rapidly rotating stars

Most of the information that we can obtain from stars stems from the light that they emit. Variations of this light can be monitored for many stars, enabling to detect periodic patterns produced by oscillation modes of stars. In the case of the Sun, it has been possible to theoretically construct these modes and compare them with observations, giving crucial pieces of information on the internal structure [1]. Ultra-precise photometric data from the recent space missions COROT [2] and Kepler [3] include many rapidly rotating stars, for which a theory of oscillation modes is needed to infer physical properties of their internal structure [4].

Stellar oscillations are divided into two categories that are separated in frequencies, with gravity modes below the Brunt-Väisälä frequency and pressure (acoustic) modes in the high part of the frequency spectrum. Recently, ray dynamics and semiclassical techniques have been used to describe oscillation modes in rapidly rotating stars [5,6]. Indeed, these oscillations have a short-wavelength limit in the same way as quantum or electromagnetic waves do, and the ray dynamics is also governed by Hamiltonian equations of motion [7]. In the case of acoustic waves, which are the focus of this letter, numerical simulations of this dynamics for a polytropic stellar model showed that when the rotation rate increases, stable and chaotic regions coexist in the phase space. As a consequence, the stationary acoustic modes can be divided into modes localized in the regular zones or in the chaotic zones, with markedly different properties [5]. For regular island modes an asymptotic formula was built in [8], showing that they are characterized by regular spacings. Chaotic modes, which have been studied in the context of quantum mechanics by the field of quantum chaos [9], are expected to be distributed in accordance with the predictions of random matrix theory (the relevant matrix ensemble for our system is the Gaussian Orthogonal Ensemble (GOE)) [10,11], and should not display regularities.

In this letter we show that chaotic modes in models of rapidly rotating stars display pseudo-regularities which can be clearly seen from peaks in the autocorrelation of the spectra. They could correspond to peaks extracted from the observed frequency spectra of δ Scuti stars [12], a class of rapidly rotating pulsating stars. Understanding their origin would then be key to derive physical constraints on star interiors. These peaks are also of interest from a theoretical standpoint, since they are unseen in the autocorrelation of GOE spectra [11]. Here, we explain the pseudo-regularities in the chaotic mode spectrum from a theory based on the general properties of acoustic rays in stars, coupled with the semiclassical theory of correlations in chaotic spectra [13].

To model the star, we study a self-gravitating monoatomic perfect gas of adiabatic exponent $\Gamma = 5/3$ in solid body rotation. The equilibrium model is a simplified model, where the pressure and density satisfy the relation $P \propto \rho^{1+1/N}$ with an effective polytropic index N = 3. The acoustic ray dynamics is derived making standard approximations valid for high frequencies: we neglect the Coriolis force, the perturbations of the gravitational potential and the effect of viscosity and thermal diffusivity. In the linear limit, one obtains a Helmholtz-type equation. The pressure perturbations $\Psi(\bm{x},t)$ satisfy

$$\begin{equation} \left(\frac{\partial^2}{\partial t^2} + \omega_c^2 \right) \Psi(\bm{x},t) - c_s^2 \, \nabla^2 \Psi(\bm{x},t) = 0, \end{equation} \tag{ 1 }$$

where c_s is the sound speed, that decreases from the center to the surface as the square root of the temperature and $\omega_c$ is a cutoff frequency which increases sharply close to the boundary, confining the wave inside the star [14].

We seek solutions of the form $\Psi = e^{i \Lambda \phi(\bm{x},t)}$ where $\Lambda^{-1}$ is a small dimensionless parameter and get an equation for the phase $\phi(\bm{x},t)$ [7]:

$$\begin{equation} -\Lambda \dot{\phi}(\bm{x},t) = (\omega_c^2 + c_s^2 k^2)^{1/2} = H, \end{equation} \tag{ 2 }$$

where $\bm{k} = \bm{\nabla}{\phi}$ is the wave vector and H is the ray dynamics hamiltonian. As is usual in stellar physics [4], rotation frequencies will be given in units of the Keplerian rotation rate $\Omega_k = (GM/R_{eq}^3)^{1/2}$ where the centrifugal acceleration equates gravity at the equator, G being Newton's constant, M the mass and R_eq the equatorial radius of the star. We will express acoustic frequencies in terms of $\omega_p = (GM/R_p^3)^{1/2}$, where R_p is the polar radius. Because R_p varies slowly with rotation, the choice of $\omega_p$ is well suited to compare acoustic frequencies at different rotation rates.

The equations of motion are solved using a 5th-order Runge-Kutta method. As the system is cylindrically symmetric, the ray dynamics can be reduced to a two-dimensional problem [14]. In the particular case of axisymmetric modes, acoustic rays propagate in a meridional plane of the star. The study of the dynamics through the Poincaré Surface of Section (PSS) reveals three main types of structures for a wide range of rotation rates (see fig. 1): stable islands are built around periodic orbits, whispering gallery rays remain close to the surface and chaotic ergodic trajectories fill all available space. Stationary modes are associated to these different phase space regions. As presented in [14], modes associated with the stable islands have simple spectra of the form $\omega_{n \ell} = n \delta n + \ell \delta \ell + \alpha$ and in the case of the most important 2-period island modes this formula was explicited by an asymptotic theory [8].

Zoom In Zoom Out Reset image size

By contrast, chaotic modes associated with the ergodic phase space region are not predicted in general to follow any simple asymptotic formula. According to the Bohigas-Giannoni-Schmit conjecture [10] chaotic mode spectra should have correlations given by GOE, that have the property of repelling each other at short distance. Conversely, generic regular modes should follow a Poisson distribution with no level repulsion [15]. Using a two-dimensional code that computes the stationary oscillations of a rotating star [16] we produced the frequency spectra of the polytropic model at five different rotations (namely the four rotations shown in fig. 1 and additionally $\Omega/\Omega_k$ = 0.658, whose PSS is very similar to panel c but with a disappearing 6-periodic island chain). All modes are axisymmetric and either symmetric (even) or antisymmetric (odd) with respect to the equator. We then selected the chaotic modes by removing from the full spectra the island modes and whispering gallery modes¹ . In fig. 2 we display the distribution of the ratio of consecutive level spacings $r_n = (\omega_{n+1}-\omega_n)/(\omega_n-\omega_{n-1})$ of chaotic modes and in table 1 we indicate the average value of $\tilde{r}_n = \mathrm{min}(r_n, 1/r_n)$ [17,18]. The main panel of fig. 2 shows the aggregated ratio distributions of four rotations. This procedure is standard to smooth out the statistical fluctuations due to small datasets [19]. The resulting distribution is in good agreement with the GOE prediction for chaotic systems. This is not the case for the $\Omega / \Omega_k = 0.706$ spectrum, shown in the inset, as the ratio distribution does not show level repulsion at small r values and the averaged $\tilde{r}$ falls between the GOE and Poisson values. Such anomalous statistics is often associated to the fact that the spectrum contains independent subspectra. We will explain this peculiarity below by the presence of partial barriers in the chaotic zone.

Table 1:.  Average value of $\tilde{r}$.

Poisson	$\approx 0.38629$
GOE	$\approx 0.53590$
$\Omega / \Omega_k = 0.706$	$\approx 0.4743$
Other rotations	$\approx 0.5618$

Zoom In Zoom Out Reset image size

The ratio of consecutive level spacings is a short-range quantity; to investigate the correlations at longer range we computed the two-point autocorrelation function $\text{R_2}(\xi) = \left<\text{d}(\omega-1/2 \xi) \, \text{d}(\omega+1/2 \xi) \right>$, where $\text{d}(\omega)$ is the spectral density and $\left<.\right>$ is a frequency average defined by $\left< f \right> = \int f \text{d}\omega$. As shown in fig. 3, a deviation from GOE statistics [11] appears in the form of a peak centered at value $\Delta \simeq 1.15 \, \omega_p$ (panel (a)). This peak is robust in the sense that it appears at every rotation rate, though its position shifts slightly (see fig. 3, panels (b), (c), (d)). Other peaks emerge from the noise, especially as rotation increases (fig. 3, panels (c) and (d)) around $1/3 \Delta$ and $4/3 \Delta$. In the following, we will explain the origin of the most robust peaks and their evolution with respect to rotation, using semiclassical methods. We will also link the presence of the additional peaks to the existence of partial barriers.

Zoom In Zoom Out Reset image size

Semiclassics (here large-ω limit) is built on the Gutzwiller trace formula, which relates the spectral density to a sum over periodic orbits of the system [9], and was introduced in the framework of quantum mechanics. We have adapted the derivation to the case of acoustic modes; it leads to the formula $\text{d}(\omega)-\text{\bar{d}}(\omega) = \mathrm{Re}\sum_i A_i(\omega) \, e^{\text{i} S_i (\omega)}$ where i labels periodic orbits, S is the action, A is an amplitude, supposed to be slowly varying with ω, which depends on the length and stability of the orbit and $\text{\bar{d}}(\omega)$ is the smoothed density. In our system the action can be expressed as $S = \omega \int \text{d}s \, \sqrt{1-\omega_c^2 / \omega^2} / c_s$, where s is the curvilinear coordinate along a given ray path. Therefore $T = \partial S / \partial \omega$ corresponds to the travel time (or acoustic time) of a ray. In order to relate the autocorrelation of the chaotic spectra to the acoustic ray dynamics, we follow the method of Berry [13], which uses the so-called diagonal approximation. The idea is to consider the form factor:

$$\begin{equation} \text{K}(T) = \frac{1}{\sqrt{2 \pi}}\int^{\infty}_{-\infty} \text{d} \xi \, \text{exp}(\text{i} \xi T) \, \text{C}(\xi) \end{equation} \tag{ 3 }$$

which is the Fourier transform of the autocorrelation:

$$\begin{equation} \begin{array}{rcl} \text{C}(\xi)&=& \Big\langle[\text{d}(\omega-1/2 \xi)-\text{\bar{d}}(\omega-1/2 \xi)] \\ && \times [\text{d}(\omega+1/2 \xi)-\text{\bar{d}}(\omega+1/2 \xi)] \Big\rangle \\ &=& \Big \langle \Big(\text{Re} \sum\limits_i A_i e^{\text{i} S_i(\omega - 1/2 \xi)}\Big)\\ && \times \Big(\text{Re} \sum\limits_j A_j e^{\text{i} S_j(\omega + 1/2 \xi)}\Big) \Big \rangle. \end{array} \end{equation} \tag{ 4 }$$

For short times, the frequency average eliminates the off-diagonal terms $i \neq j$ by phase incoherence, allowing to approximate the form factor as $\text{K}(T) \approx \sum_j A_j^2 \delta(T-T_j)$, i.e., the density of periodic orbits weighted with intensities $A_j^2$ [13,20]. This approximation fails for long times, typically higher than the Ehrenfest time [21], because there are pairs of orbits with very close actions $S_i \approx S_j$ and higher-order terms need to be computed [22–25].

The next step of Berry's method is to use the general result of Hannay and Ozorio de Almeida [26] which implies that $\sum_j A_j^2 \delta(T-T_j) \propto T$ for sufficiently large T. To make this result explicit we express the amplitude and the density as functions of the metric and topological entropies h_M, h_T: $A (T) \propto T e^{-(1/2) h_M T}$ and $\sum_j \delta(T-T_j) \rightarrow \rho(T) \approx (1/T) \, e^{h_T T}$ and suppose these two entropies to be equal [20,27]. These formulas give $\text{K}(T) \rightarrow A^2(T) \; \rho(T) \propto 2 T$, which is the GOE prediction at leading order (the factor of two comes from the fact that orbits are twice degenerate due to time-reversal symmetry [13]). In our system, however, the peculiar distribution of periodic orbits $\rho(T)$ modifies the form factor behavior.

Like classical trajectories in billiards, acoustic rays bounce on the reflective caustic close to the surface. Thus any periodic orbit can be divided into chords, each chord connecting two successive surface bouncing points. Periodic orbits thus belong to a class of trajectories that we call n-chords, whose endpoints have to lie on the surface. While it is extremely difficult to find all periodic orbits of our system, we can nevertheless infer some of their properties from the study of large samples of n-chord trajectories.

Figure 4 shows the distribution of travel times for n-chords in the chaotic region of phase space at rotation $\Omega / \Omega_k = 0.589$. In stars, the travel time is not proportional to the geometric length, as, e.g., in billiards. Indeed, since the sound velocity is much smaller near the surface of the star, all trajectories spend much more time near the surface than in the core. This results in a travel time distribution that concentrates around specific values evenly spaced out (see fig. 4). As the number of chords increases the individual packets grow wider (similarly to the law of large numbers), until the width becomes much larger than the interpeak distance, and the distribution turns into a smooth curve. This is quantified by the ratio $\sigma_0/T_0$, $\sigma_0$ being the standard deviation of the first packet and T₀ the average distance between the mean values of consecutive packets. In the case of a billiard whose boundary is shaped like our star surface, $\sigma_0/T_0 \approx 0.32$ meaning that the packet structure disappears after a few rebounds. In the stellar case $\sigma_0/T_0$ grows with rotation but remains below 0.08. The packets are thus discernable for much longer times.

Zoom In Zoom Out Reset image size

We model the distribution of fig. 4 as $\text{P}(T) = \sum_n \text{P}_n(T)$, with $\text{P}_n(T)$ the probability distribution of n-chords travel times:

$$\begin{equation} \text{P}_n(T) = \frac{1}{\sqrt{2 \pi n} \, \sigma_0} \, \exp \left(- \frac{(T-n \, T_0)^2}{2 (\sqrt{n} \, \sigma_0)^2} \right), \end{equation} \tag{ 5 }$$

where the values of T₀ and $\sigma_0$ are taken from the numerical simulations of acoustic rays. We checked that the packets have indeed a mean value close to nT₀ and that their standard deviation grows as $\sqrt{n}\sigma_0$. Physically, T₀ is the characteristic travel time of a single chord. Thus, T₀ and $\sigma_0$ are directly related to the star volume and the sound speed profile in the star interior. In particular, the small variance of the travel times is a consequence of the strong decrease of the sound speed near the surface. As mentioned before we cannot find directly the periodic orbits, but we claim that the constraints imposed by $\text{P}(T)$ are strong enough to explain the correlations seen in the spectra.

To take into account the exponential growth rate of the number of periodic orbits with T, the density of periodic orbits in our system is thus modeled as $\rho(T) \propto (1/T) \, e^{h_T T} \times \text{P}(T)$. Assuming $A(T) \propto T e^{-(1/2)h_M T}$ gives the result

$$\begin{equation} \text{K}(T) \rightarrow A^2(T) \; \rho(T) \propto T \; \text{P}(T). \end{equation} \tag{ 6 }$$

As the variance of the $\text{P}_n(T)$ grows with n, $\text{P}(T)$ becomes flat for sufficiently long times. For such long times we recover the sum rule [26] $\sum_j A_j^2 \delta(T-T_j) \propto T$, which is consistent with the existence of GOE statistics at short spectral distance (fig. 2). For shorter times $\text{P}(T)$ creates a specific regime departing from GOE.

The Fourier transform

$$\begin{equation*} F(\xi) = 1/(\sqrt{2 \pi}) \int_{-\infty}^\infty \text{d}T \, \exp(-i \xi T) \, T\text{P}(T) \end{equation*}$$

is shown in the left inset of fig. 5 for different rotations. It shows a peak, as in the spectral data (fig. 3), that moves towards low frequencies for increasing rotation. The peaks extracted from the numerical spectral data and obtained from semiclassical theory are in good agreement (fig. 5, main panel). The discrepancy of about $5 \%$ is similar to what was obtained for island modes in [8], and can be attributed to the relatively low values of numerically computed frequencies, as the theoretical values are derived under the semiclassical (high-frequency) approximation. The model also predicts that the peak width should show a slow increase with rotation, as can be seen in the left panel of fig. 5. In the right panel, we show that the peak height decreases rapidly with $\sigma_0/T_0$; this explains why the peak is not usally visible in chaotic systems with large $\sigma_0 / T_0$.

Zoom In Zoom Out Reset image size

The semiclassical theory outlined above explains the origin of the main peak in the autocorrelation. Other peaks are sometimes visible in the autocorrelations of fig. 3, especially for the rotation $\Omega / \Omega_k = 0.706$ where P(r) deviates strongly from GOE (see fig. 2). We attribute these peaks to the presence of partial barriers in the chaotic zone around the main stable island. Partial barriers isolate some zones of phase space, from which trajectories escape more slowly than in the rest of the chaotic zone [28,29]. To find them we compute trajectories with initial condition near the main island and monitor their evolution on the PSS. Without barrier all points would spread ergodically, however fig. 6 shows that a subset of points remains near the island for a long time. Close to $\Omega / \Omega_k = 0.706$, the region enclosed by partial barriers grows in size and traps trajectories more efficiently. For the relatively low-frequency waves here considered, they will act as barriers and quantize independently a subset of modes. This will weaken the level repulsion at short spectral distance as seen in fig. 2. Besides, as trapped trajectories revolve around a 6-periodic orbit, we expect the modes to quantize like island modes [8], leading to a $\Delta/3$ regularity.

Zoom In Zoom Out Reset image size

In this letter, we have shown that specificities in the distribution of the actions of periodic orbits can create a peak in the autocorrelation of chaotic spectra. Such a peak can potentially be detected in the frequency spectrum of chaotic pressure modes in rapidly rotating stars. Recent data have confirmed the existence of peaks in the autocorrelation spectrum of the rapidly rotating δ Scuti stars [12]. They have been attributed so far to regular island modes but our results indicate that chaotic modes should also produce such peaks. These two kinds of peaks should be close² but still distinct at most rotation rates. Whether these two peaks could be discernable in observed data remains to be investigated. In addition to astrophysical observations, the phenomenon described here could be tested with an experimental setup such as, e.g., electromagnetic waves in a cavity with a strong gradient of the refractive index along the radial direction.

We thank CALMIP for providing computational resources and the Université Paul Sabatier for funding. This work has made use of the TOP code developed by D. Reese and of a user-friendly version of this code developed by B. Putigny.