Deciphering the oscillation spectrum of $\gamma$ Doradus and SPB stars

The space-based Kepler mission provided four years of highly precise and almost uninterrupted photometry for hundreds of $\gamma$ Doradus stars and tens of SPB stars, finally allowing us to apply asteroseismology to these gravity mode pulsators. Without rotation, gravity modes are equally spaced in period. This simple structure does not hold in rotating stars for which rotation needs to be taken into account to accurately interpret the oscillation spectrum. We aim to develop a stellar-model-independent method to analyse and interpret the oscillation spectrum of $\gamma$ Dor and SPB stars. Within the traditional approximation of rotation, we highlight the possibility of recovering the equidistance of period spacings by stretching the pulsation periods. The stretching function depends on the degree and azimuthal order of gravity modes and the rotation rate of the star. In this new stretched space, the pulsation modes are regularly spaced by the stellar buoyancy radius. On the basis of this property, we implemented a method to search for these new regularities and simultaneously infer the rotation frequency and buoyancy radius. Tests on synthetic spectra computed with a non-perturbative approach show that we can retrieve these two parameters with reasonable accuracy along with the mode identification. In uniformly rotating models of a typical $\gamma$ Dor star, and for the most observed prograde dipole modes, we show that the accuracy on the derived parameters is better than 5% on both the internal rotation rate and the buoyancy radius. Finally, we apply the method to two stars of the Kepler field, a $\gamma$ Dor and an SPB, and compare our results with those of other existing methods. We provide a stellar-model-independent method to obtain the near-core rotation rate, the buoyancy radius and mode identification from g-mode spectra of $\gamma$ Dor and SPB stars.


Introduction
In stellar physics, the treatment of interfaces between convective and radiative regions remains rather simplistic and subject to large uncertainties despite the critical effects on the stellar structure and evolution. This is especially true for main-sequence stars that possess a convective core. Indeed, mixing processes at the interface can extend the central mixed zone and increase the amount of hydrogen fuel available for nuclear fusion reactions, that is, they can increase the time spent on the main sequence. A better treatment of these interfaces will thereby result in more precise determinations of stellar ages. In standard models, the convective core boundary is set by the Schwarzschild criterion that is equivalent to the location where the acceleration of convective motions cancels out. Comparison of standard stellar models with observations of eclipsing binaries (Andersen et al. 1990) or open clusters (Maeder & Mermilliod 1981) showed that this criterion actually underestimates the size of the convective core and proved the need for including extra-mixing processes in stellar models. If this fact is well-established, the nature of these processes (overshooting, rotationally induced mixing, internal gravity waves), the extent of the extra-mixed region, and their dependence on stellar parameters are still not fully understood.
Another major uncertainty of stellar physics is the distribution and evolution of angular momentum (e.g. Meynet et al. 2013). Rotation distorts stars and triggers supplementary hydrodynamical instabilities, which mixes stellar interiors and feeds back into rotation by transporting angular momentum. In particular, meridional circulation and shear instabilities are commonly accepted as operating mechanisms in stellar interiors. However, several studies (Chaboyer et al. 1995;Maeder & Zahn 1998;Mathis & Zahn 2004) pointed out that these mechanisms are insufficient to reproduce the solar rotation profile as measured by helioseismology (e.g. Schou et al. 1998;García et al. 2007;Fossat et al. 2017). To explain the discrepancy between observations and rotating models, several additional processes of angular momentum transport have been proposed that either involve internal gravity waves (Charbonnel & Talon 2005) or magnetic fields (Eggenberger et al. 2005), but whether one is dominating over the others remains somewhat unclear. Similar uncertainties hold in evolved low-mass stars. The internal rotation rates of red giants, as measured by asteroseismology, are a few orders of magnitude smaller than what is predicted by rotating stellar models (Mosser et al. 2012), indicating that one or several missing angular momentum transport processes are at work in these evolved stars (e.g. Fuller et al. 2014;Rüdiger et al. 2015;Belkacem et al. 2015;Eggenberger et al. 2017).
Among the variety of pulsating stars, γ Doradus (γ Dor) and slowly pulsating B-type (SPB) stars are promising targets for obtaining constraints on both convective boundary mixing and angular momentum transport processes. These stars pulsate in high radial order gravity modes that probe the innermost radiative layers close to their convective core, where they have larger amplitude. Notably, Miglio et al. (2008) and Bouabid et al. (2013) demonstrated that the properties of gravity modes (g modes) in these stars is particularly sensitive to the shape of the chemical gradient at the edge of the convective core. Stars of γ Dor type are late A-to early F-type main-sequence stars with masses between roughly 1.3 and 2.0 M . Their g-mode pulsations are thought to be excited by the modulation of the radiative flux at the base of their thin convective envelope, or so-called convective blocking mechanism (Guzik et al. 2000;Dupret et al. 2005). It is worth noting that a more specific interest of γ Dors lies in the fact that they are progenitors of red giants. SPB stars are main-sequence stars that have spectral types B3-B9 and masses between 2.5 and 8 M . In these stars, the pulsations are driven by the κ-mechanism due to the metal opacity bump at T ∼ 2 × 10 5 K (Dziembowski et al. 1993).
From the point of view of seismology, ground-based observations of these pulsators are especially impractical as their g-mode pulsations typically have periods of around one day. This seriously hindered the application of asteroseismology to γ Dor and SPB stars until the advent of space missions dedicated to high precision photometry. In particular, Kepler provided four years of highly precise and nearly uninterrupted photometry for many of these stars allowing the detection of g-mode series nearly equally spaced in periods as predicted by theory (Tassoul 1980;Ledoux 1951). Thus, in recent years, detailed seismic studies of slowly rotating γ Dor (Kurtz et al. 2014;Saio et al. 2015;Keen et al. 2015;Bedding et al. 2015;Murphy et al. 2016) and SPB stars (Pápics et al. 2014;Moravveji et al. 2015) targeted by Kepler could be undertaken. However, this regular structure does not hold for moderately and rapidly rotating stars which, nevertheless, constitute the majority of these pulsating variables. Indeed, the projected equatorial velocities of B-and A-type stars is commonly around 100 km s −1 but this can reach 250 km s −1 in some stars (Abt et al. 2002;Royer et al. 2007).
Moderate or rapid rotation significantly affects the g-mode oscillation spectrum because of the Coriolis force. Observationally, the Fourier spectrum of the photometric time series often contains hundreds of peaks with no apparent structure, at least not as predicted by classical perturbative treatments of rotation. This leads Van Reeth et al. (2015a) to develop a new algorithm to search for non-equally spaced g-mode series during the pre-whitening process. Successfully, they reported the detection of period spacing patterns in 50 γ Dor stars (Van Reeth et al. 2015b). To identify the g modes and estimate their interior rotation rates, Van Reeth et al. (2016) subsequently modelled the patterns following a direct approach. Using the asymptotic formulation of the traditional approximation that accounts for the main effects of the Coriolis force, the authors computed the gravity mode periods in a grid of stellar models representative of the γ Dor instability strip and then fitted the observed patterns with a least-square minimisation. Following a similar approach, Pápics et al. (2017) conducted a detailed seismic study of five rotating SPB stars in the Kepler field.
Given the uncertainties contained in stellar models, we propose a different method that is model-independent and solely relies on the asymptotic traditional approximation. Section 2 revisits the theoretical background of the traditional approximation and introduces the details of this method. In Sect. 3, we validate it on synthetic spectra of representative γ Dor models in which oscillations are computed with a two-dimensional (2D) non-perturbative approach. We also assess the biases on our estimates of the buoyancy radius and near-core rotation frequency. In Sect. 4, we confirm the potential of such a method on two already studied Kepler targets, the γ Dor star KIC12066947 and the SPB star KIC3459297. Finally, we discuss the limitations of the method in its current implementation and conclude in the last section.

Theoretical background
Formally, the pulsation periods are obtained by solving the equations of hydrodynamics perturbed by small fluctuations around an equilibrium structure. Two approaches may be adopted: solving the whole system of equations numerically or conducting an asymptotic analysis assuming reasonable approximations. With this latter approach, in a non-rotating spherical star, Tassoul (1980) demonstrated that the pulsation periods of high-order gravity modes (|n| ) can be well-approximated to first order by the expression where n is the radial order, the angular degree, and m the azimuthal order of the pulsation mode. We note that, by convention, radial orders of g modes are negative. We adopt the convention that m > 0 are prograde modes and m < 0 retrograde modes. The phase term is near-constant and linked to the star's structure. The buoyancy radius P 0 can be expressed as, where N BV is the Brunt-Vaisala frequency and R represents the resonant cavity of g modes. Hence, modes of same degree and consecutive radial orders n are equally spaced in period, which is characterised by constant period spacings defined as, This is a particularly interesting property of high-order g modes when searching for mode series in an observed spectrum. The pulsation periods do not depend on m.
In a rotating star, pulsations are described by a coupled set of equations that is both computationally expensive and numerically complex to solve. First introduced in the context of geophysics (Eckart 1960) and later applied to stellar pulsations (e.g. Lee & Saio 1987, 1997Townsend 2003), the traditional approximation of rotation (TAR) is a non-perturbative, but simplified, treatment of rotation that allows the separability of these equations while still acknowledging the main effects of the Coriolis force. Assuming the star is in solid-body rotation, the TAR neglects the horizontal component of the rotation vector. In other words, the radial motions due to the Coriolis force and the radial component of the Coriolis force associated with horizontal motions are discarded. This approximation is well-justified for low-frequency pulsations as it is discussed here (Berthomieu et al. 1978), except near the stellar centre where radial and horizontal motions become comparable. To obtain the separability of the system, the Cowling (1941) approximation, which neglects the perturbation of the gravitational potential, is further made. Finally, the centrifugal distortion is neglected ensuring the spherical symmetry of the star's structure. This latter hypothesis can be justified by the fact that the g modes are mostly sensitive to the innermost radiative layers, where centrifugal distortion remains minimal.
Thus, in the reference frame in co-rotation with the star, the equations of pulsations become separable in the three spherical coordinates (r, θ, ϕ). The dependence in latitudinal coordinate is governed by the Hough functions, θ ,m , which are the eigenfunctions of the well-known Laplace's tidal eigenvalue problem, where L s is a linear operator, λ ,m are the eigenvalues of the problem, and µ = cos θ. A complete derivation of the problem is derived for example in Lee & Saio (1997). Both the Hough functions and associated eigenvalues depend on , m, and the spin parameter s = 2P co n, ,m /P rot , where P co n, ,m is the mode period in the co-rotating frame and P rot is the rotation period of the star. This latter reduced parameter appears naturally in the frame of the TAR and is indicative of the influence of the Coriolis force on pulsation modes. The equation for the radial component is similar to that of the non-rotating case, except that ( + 1) are replaced by the new eigenvalues λ ,m . The azimuthal dependence stays identical.
As their form is akin to the non-rotating case, the asymptotic analysis of Tassoul (1980) may also be applied to the pulsation equations obtained within the TAR. It follows a more general asymptotic formula for periods of high-order g modes in the corotating frame, The mode periods in the inertial frame are subsequently deduced from P in n, ,m = 1 ν co n, ,m + mν rot = P co n, ,m 1 + mν rot P co n, ,m , where ν co n, ,m is the mode frequency in the co-rotating frame and ν rot is the cyclic rotation frequency. Therefore, in rotating stars, the degeneracy in azimuthal order is lifted and the regular structure of the spectrum no longer holds as period spacings now also depend on the spin parameter.

Behaviour of gravity modes in rotating stars
We recall here the behaviour of high-order g modes in rotating stars. We refer the reader to the works of Miglio et al. (2008), Bouabid et al. (2013), and Ouazzani et al. (2017) for a more complete and detailed picture. In slowly-rotating stars (s 1), the spectrum is organised in well-separated multiplets (n, ). Information on the rotation of the star is carried by the rotational splittings, which are expressed by δν n, = ν n, ,m − ν n, ,0 in terms of the observed mode frequencies. At moderate rotation, when the rotation period approaches that of the modes (s ∼ 0.1), these multiplets start to overlap and it becomes difficult to disentangle them visually. Nonetheless, period spacings remain almost constant, which makes the detection of a regular structure still possible (see e.g. Bedding et al. 2015). For  Fig. 1. Illustration of the rotational shift of observed pulsation periods, in the inertial frame, for four rotation rates. The oscillation spectra were computed within the asymptotic TAR (see Eq. (5)) using P 0 = 4320 s, a value typical of γ Dors. Black, orange, and blue bars are dipole prograde (m = 1), zonal (m = 0), and retrograde (m = −1) modes, respectively. They are shifted vertically for clarity. rapid rotation, when the rotation period is of the same order (or greater) than the pulsation periods, the structure of the spectrum is fundamentally different. Indeed, rotation tends to separate modes according to their geometry, that is, according to their values of ( , m). In the inertial frame, prograde and zonal modes are shifted towards shorter periods whereas retrograde modes tend to drift towards longer periods. If the rotation frequency is sufficiently high, ( , m) modes form clusters of modes distinct from each other. These clusters partially overlay otherwise. This evolution of the g-mode spectrum with rotation is illustrated on Fig. 1.
Rotation also leaves its imprint on the period spacings. Indeed, in the inertial frame, the period spacings of prograde and zonal modes decrease almost linearly with increasing mode periods. At fixed rotation frequency, the slope of this linear trend is smaller for prograde than for zonal modes. Also, for a given star, this slope get steeper with increasing rotation rate. As for retrograde modes, the change of reference frame and rotation have competitive effects on the mode periods, and so, on the period spacings, resulting in more complicated behaviour. However, the spacings tend to have an upper trend in the ∆P in − P in plan for the most part.
As γ Dor and SPB stars evolve on the main sequence, a chemical composition gradient develops near the convective core boundary. If no mixing processes smooth this gradient, this leads to a local sharp variation of the Brunt-Vaisala frequency. Such a buoyancy glitch periodically traps the g modes in the gradient region. As a result, the period spacings slightly oscillate around the linear trend. The period of this oscillatory component is related to the location of the glitch while the amplitude gives information about its abruptness (Miglio et al. 2008;Bouabid et al. 2013).

The method
In the frame of the asymptotic TAR, we highlight the possibility of recovering the equidistance of the period spacings by stretching the pulsation periods. Indeed, rearranging the terms of Eq. (5),  Fig. 2. Stretching of the pulsation periods in the example of a series of dipole prograde modes. Mode periods were generated using the asymptotic TAR with ν rot = 20 µHz and P 0 = 4320 s (or 0.05 d). Top panel: before the stretching as they would be seen by an observer. Bottom panel: after the change of reference frame and the stretching, assuming the correct mode identity and rotation rate.
shows that, by multiplying the period scale by the square root of the Laplace's eigenvalue λ ,m (s), where s matches the star's rotation frequency, the associated ( , m) modes become regularly spaced of P 0 . This is illustrated in Fig. 2 where we applied the stretching to a synthetic series of prograde dipole modes ( = 1, m = 1). On the basis of this property, we implemented a method to search for these new regularities in the g-mode spectrum of γ Dor and SPB stars, and simultaneously infer ν rot and P 0 . Here, we describe the methodology used in detail. Similar stretching techniques have been developed to interpret the mixed mode oscillation patterns of red giant stars (Mosser et al. 2015).
First, while not absolutely necessary, we carry out the frequency analysis of the oscillation spectrum (e.g. by prewhitening the periodogram). This simplifies greatly the subsequent steps and allows us to identify possible combination frequencies that arise from non-linear processes occurring in the star. The periodogram of γ Dor and SPB stars is typically dense and often contains a large amount of frequency peaks, where several series of ( , m) modes may coexist. The equidistance of the period spacings can only be recovered for a given series of ( , m) modes at once. In other words, for the method to work correctly, the ( , m) modes searched for need to be dominant in the list of frequencies used for the analysis. The strategy adopted here consists in selecting a part of the spectrum where we expect modes of given ( , m) to be prevailing, relying on the behaviour of g modes in rotating stars (Sect. 2.2). In the current implementation, this step is performed by visual inspection of the periodogram or of the period distribution. Typically, we look for frequency groupings or mode density variations. Such filtering is obviously subjective and may be a matter of trial and error when the g-mode spectrum appears especially disordered. Automation of this step is possible and is currently under study. At this point, we consider that a given mode series is dominant over other series possibly present in the spectrum.
Once the list of mode periods to analyse is established, we pick a guess for and m and choose a range of rotation frequency to test. In practice, surface cancellation effects limit greatly the visibilities of high degree modes ( 4). As for the rotation rate, we restrain the interval to 0-35 µHz, which largely contains all the rotation frequencies measured in γ Dors and SPBs until now. Then, for each rotation rate, we compute the pulsation periods in the co-rotating frame using Eq. (6) and stretch the spectrum according to Eq. (7).
To detect a regularity, we subsequently compute the power spectral density (PSD) of each stretched spectrum from the Discrete Fourier Transform (DFT). For a frequency f , this is given by where P co i are the mode periods in the co-rotating frame and N the total number of modes in the list. The DFT spectra obtained are then stacked vertically by increasing rotation rate to build a DFT map representative of the space of parameters explored. Figure 3 shows an example of such a map. The detection of a regularity materialises into a characteristic ridge of high PSD. Such a ridge is indicative of the correlation between P 0 and ν rot . Indeed, because ν rot and P 0 are related through Eq. (5), the effect of a small change in ν rot can be nearly, but not exactly, counterbalanced by a small change in P 0 .
Finally, to ensure the detection is not coincidental, we compare the maximum of PSD to a threshold value computed from a false-alarm probability p = 0.01 of having a peak generated by random noise. Further details are provided in Sect. 2.4. If greater than this threshold, the maximum of PSD is used as an estimator of ν rot and P 0 (or more exactly 1/P 0 ). The period échelle diagram of the stretched spectrum may then be plotted to have an overview of the modes actually participating in the regularity detected. Otherwise, the trial and error process may be continued by changing either the filtering of mode periods, the guess for and m, or the interval of rotation rate tested. Briefly, the algorithm proceeds as follows: 1. Establishing a list of frequencies (e.g. from Fourier analyses). 2. (Optional) Filtering of pulsation modes according to their period values.
A47, page 4 of 14 3. Pick a guess for ( , m) and choose the interval of rotation frequencies to test. 4. For each rotational frequency: (a) Switch from the inertial to the co-rotating frame.
(b) Stretch the spectrum. (c) Compute the DFT. 5. Stack the DFT spectra obtained on top of one another by increasing rotation rate (DFT map). 6. Check if the maximum of PSD is significant.
(a) If significant: formally identify modes and estimate ν rot and P 0 from the maximum of PSD. (b) If not: continue the trial and error process by returning to step 2 or 3.

Detection threshold
We derive here the detection threshold needed to achieve a specified probability of false alarm p for a peak found in a DFT map. This p-value corresponds to the probability for that peak to be generated by pure noise. Each line of the map is computed as the norm square of the DFT of a stretched spectrum. From Eq. (8), we can show that the noise in the power spectrum asymptotically (i.e. for large N) follows a two-degrees-of-freedom χ 2 statistics (e.g. Abramowitz & Stegun 1964) when the period set P co i is purely random 1 . These statistics are correct as long as N, the number of periods in the set, is large enough. We have numerically verified that a dozen is sufficient in practice, which is often the case for real observations. The normalisation used for the DFT in Eq. (8) ensures that the mean and variance of the statistics is 1. Thus, the probability for pure noise to generate, in a given DFT bin, a peak larger than a threshold T is p = exp −T . As a consequence, given M independent bins in the DFT, the detection threshold is The number of independent bins is then M = ∆ f /δ f , where δ f is the spectral resolution and ∆ f is the frequency window of interest, that is, the range in which we look for 1/P 0 . The spectral resolution expresses δ f = λ ,m P co max − λ ,m P co min −1 where P co max and P co min are the largest and smallest periods in the set. For p = 0.01 and a typical value of M = 50, the detection threshold is then T = 8.5.

Estimating the uncertainties
Due to the rotation-pulsation coupling, the mode identity, the rotation frequency, and the buoyancy radius are intrinsically related, which can give rise to degeneracies in the process of mode identification. Indeed, several mode identities can be found for a given ( 1 , m 1 ) mode series. This is due to the fact that, assuming another geometry ( 2 , m 2 ), it is sometimes possible to mimic the ( 1 , m 1 ) pattern by adjusting the values of ν rot and P 0 accordingly. This problem has already been noticed by Van Reeth et al. (2016). To break the degeneracy, the authors compare the found value of the asymptotic period spacing ∆P (or equivalently P 0 ) to expected values in models, as a consistent check. We used the same strategy in this work. The buoy-1 Since P co i are random set, both real and imaginary parts of the DFT follow normal distributions according to the Central Limit theorem when N is large enough. Hence, by definition, the norm square follows a two-degrees-of-freedom χ 2 distribution. ancy radius adopts a sufficiently narrow range of values in γ Dor (roughly between 3500 and 5000 s) and SPB stars (approximatively from 5000 to 11 500 s) to infer the correct mode identity in most cases (Miglio et al. 2008). If that is insufficient, other expected properties of g modes can be considered such as the pulsation periods of excited modes or the slope of the pattern in the ∆P − P plan.
The rotation frequency and buoyancy radius as determined by our method is affected by two main sources of uncertainties: errors on oscillation mode periods and the adequacy of the asymptotic TAR to assess the properties of gravity modes. The former uncertainties can be taken into account by simply propagating the errors throughout the analysis. On the other hand, the evaluation of errors caused by the use of the asymptotic TAR is more complicated as it requires knowledge of the "true" parameter values. It is nonetheless possible to assess the effect of a sharp feature of the stellar structure, such as a buoyancy glitch, or any other effect that does not affect the global trend of the period spacing patterns. This can be treated in the same way as a statistical error on the mode periods would be. Here, we describe the procedure used to evaluate the impact of these two types of uncertainties on the estimate of ν rot and P 0 . Other potential error sources or biases will be addressed in Sect. 3.
As a first step, it is necessary to quantify the relative contributions from each source. This is achieved by performing a first analysis of the oscillation spectrum. From the values of ν rot and P 0 obtained, we compute the mode periods in the asymptotic TAR using Eq. (5) and identify each observed mode by comparing the two period values. If the standard deviation of the residuals P n, ,m − P TAR n, ,m is significantly larger than the mean error on the observed mode periods, we consider it to be the dominant source of errors. Otherwise, we keep the individual uncertainty on each observed mode period that was determined during the mode extraction process.
Secondly, we propagate these errors by means of a Monte-Carlo simulation. To this end, we draw 500 random samples of the observed pulsation periods assuming the errors on them follow a normal distribution and are not correlated. In case the uncertainty on pulsation periods is not the dominant source of error, we adopt the standard deviation of the residual P n, ,m − P TAR n, ,m as being the error bar on all mode periods. Then we apply the method to each sample to get the distributions of ν rot and P 0 . The final values of the rotation rate and buoyancy radius are determined from the maximum of PSD during the first analysis and their error bars are evaluated as the 1-σ deviation of the distributions.

Tests on synthetic spectra
As validation, we tested the method on a set of synthetic oscillation spectra computed for representative models of γ Doradus stars, which are presented in Sect. 3.1. In particular, we investigate the ability of the method to detect regularities in situations that are likely to be encountered in observations. This also lets us gauge potential biases introduced by the use of the asymptotic TAR as a prescription of the rotation-pulsation coupling. To this end, we selected three test cases. Firstly, we examine the simple case of a uniformly rotating star with no structural glitch (Sect. 3.2). With the second model, we assess the method's robustness against the effect of a buoyancy glitch (Sect. 3.3). Finally, the effect of differential rotation is addressed with the last model (Sect. 3.4).
We restrained our study to dipole modes only as it is mostly those which are observed in high-order g-mode pulsators. For each synthetic spectrum, we analysed the three mode series independently ensuring that the maximum of PSD found is above the detection threshold defined in Sect. 2.4. The DFT maps obtained are compiled under the form of a contour map, where, for each of three DFT map, we depict the contours at 50% and 95% of the maximum of PSD.

Stellar structure and oscillation models
Here we explain the choices made for the calculation of the stellar structure and oscillation used to generate the synthetic oscillation spectra. Let us first discuss the hypotheses made for the stellar structure calculation. Throughout this study, we opted for one-dimensional (1D) spherical stellar models. Indeed, spherically symmetric 1D calculations have been tested against complete 2D calculations for g modes in polytropic models (Ballot, et al. 2012), as well as for a model of a γ Dor star (Ouazzani et al. 2017). According to them, the 1D nonperturbative approach, which presents the advantage of requiring less numerical resources, gives satisfactory results compared to the full 2D approach.
Under this assumption, stellar models were computed with the stellar evolution code CLES (Scuflaire et al. 2008) for a mass of 1.86 M , and with initial helium mass fraction Y = 0.27 and metallicity Z = 0.014. We adopted the AGSS09 metal mixture (Asplund et al. 2009) and corresponding opacity tables obtained with OPAL opacities (Iglesias & Rogers 1996), completed at low temperature (log T < 4.1) with Ferguson et al. (2005) opacity tables. We used the OPAL2001 equation of state (Rogers & Nayfonov 2002) and the nuclear reaction rates from NACRE compilation (Angulo et al. 1999), except for the 14 N(p, γ) 15 O nuclear reaction, for which we adopted the crosssection from Formicola et al. (2004). Surface boundary conditions at T = T eff were provided by ATLAS model atmospheres (Kurucz 1998). Convection was treated using the mixing-length theory (MLT) formalism (Böhm-Vitense 1958) with a parameter α MLT = 1.70.
We considered models with and without turbulent diffusion. Since the CLES code does not include effects of rotation on transport of angular momentum or chemical species, we instead introduced mixing by turbulent diffusion, following the approach of Miglio et al. (2008). This reproduces an effect of rotationallyinduced mixing that is quite similar to overshooting, but in addition tends to smooth chemical composition gradients inside the star. In models including this type of mixing, the coefficient of turbulent diffusion was set to D t = 700 cm −2 s −1 and kept constant to this value during evolution and in every layer of the models. This value was selected from a previous calibration to Geneva models with similar masses. Since the evolution code does not generate rotation profiles, those are then added ad hoc after evolution calculations. When a differential rotation profile is considered, it is adapted to the structure by the mean of an error function, centred right above the convective core, with a width depending on the profile of chemical gradient.
The oscillation modes of such 1D models were computed with a non-perturbative method using the ACOR oscillation code (Ouazzani et al. 2012(Ouazzani et al. , 2015, which accounts for both the Coriolis and the centrifugal force. The ACOR code solves the hydrodynamics equations perturbed by Eulerian fluctuations, performing direct integration of the problem. The numerical method is based on a spectral multi-domain method that expands the angular dependence of eigenfunctions into spherical harmonics series, and whose radial treatment is particularly well adapted to the behaviour of equilibrium quantities in evolved models (at the interface of convective and radiative regions, and at the stellar surface). In order to determine the range of radial orders to investigate, we have relied on the non-adiabatic stability calculations provided in Bouabid et al. (2013). All the pulsations spectra studied here are calculated for g-mode radial orders between n = −50 and −20.

Simple case: Solid-body rotation and smooth period spacing patterns
We started with the simplest case and computed a model of a γ Dor star in solid-body rotation with no glitches in the structure (Model A in Table 1). Initially, we set the rotation frequency of the model at 7 µHz, which is included in the range observed in γ Dor stars (Van Reeth et al. 2016).
The top left pannel of Fig. 5 shows the contours of the three DFT maps obtained. While the three contours seem to agree on the values of ν rot and P 0 , there are still small departures from the true parameter values as measured directly from the model. Table 2 compiles the results of the three analyses. We found that relative systematic errors on the parameters do not exceed a few percent but vary according to the mode geometry considered. The stretched periodéchelle diagram for prograde dipole modes ( = 1, m = 1) is plotted on the top right panel of Fig. 5 as an example. The ridge slightly weaves along a vertical line that would represent the ridge if the asymptotic TAR was perfectly suited. Periodéchelle diagrams for other series present identical features. Although Model A has turbulent diffusion mixing, the chemical gradient at the convective core boundary is still quite substantial and gives rise to a light undulation of the period spacings.
In order to compare the systematic errors against other sources, we estimated the uncertainties on ν rot and P 0 by proceeding as described in Sect. 2.5 where we assumed a typical uncertainty on the mode periods of 1×10 −4 d. The standard deviations on the residuals (P n, ,m − P TAR n, ,m ) are 6 × 10 −5 , 1 × 10 −4 and 5 × 10 −4 d for the prograde, zonal, and retrograde modes, respectively. The uncertainties derived in this manner stay one order of magnitude smaller than the biases highlighted in Table 2, suggesting the latter is actually the critical source of errors on the parameter values.
To investigate the provenance of these systematic errors, we compared the period spacings that can be obtained with the TAR using either the parameter values determined from the method analysis (Table 2) or from the model, to those of the complete calculations. From Eq. (5), the asymptotic period spacings within the TAR can be expressed as in the co-rotating frame (see appendix of Bouabid et al. 2013, for a derivation). The spacings in the inertial frame are then derived from the expression These functions were interpolated on the period values of the complete calculations to work out the differences δ (∆P in ) = ∆P TAR in − ∆P acor in . Figure 4 displays the results of each comparison. The parameter values output by the method are those that minimise the differences δ (∆P in ), which turn out to be different from the true parameter values. In other words, the stretched spectrum using the true values is less regular than if we used the output values of the method. This means that the asymptotic TAR does not perfectly model the spacings of complete calculations.
The approximations made in the derivation of the asymptotic TAR, that is either the asymptotic approximation or the TAR (or both), then cause the disparity between the true parameter values and those recovered by the method. To investigate this further, we computed additional synthetic spectra. We set three rotation rates (7, 15, and 23 µHz) in Model A and computed the oscillation frequencies either within the non-asymptotic TAR or with ACOR. In this way, we are able to compare the relative contribution of each approximation to the bias and how it may vary with the rotation frequency. We limited this study to prograde modes ( = 1, m = 1) since other dipole modes are affected by an avoided crossing for model at 15 and 23 µHz in the range of radial orders calculated with ACOR. We note also that prograde modes are mostly those detected in γ Dor and SPB stars. The results of the analyses are listed in Table 3. The buoyancy radius is even more mis-estimated when the model is rotating quickly, while no clear trend is visible for the recovered rotation frequencies. Moreover, the retrieved values are biased whether we take the non-asymptotic TAR or complete calculations. As could be expected, the errors obtained for the ACOR calculations are more important than for those in the non-asymptotic TAR. This indicates that both the asymptotic treatment and the TAR contribute to the bias.

Effect of a buoyancy glitch
The effect of a buoyancy glitch on the pulsation periods and therefore on the period spacings is not taken into account within the asymptotic TAR. Although the global trend of the period  Fig. 4. Differences between the period spacings computed within the asymptotic TAR and those derived from complete calculations with ACOR. Orange squares: period spacings within the asymptotic TAR were computed using the parameter values of the model (ν rot = 7 µHz and P 0 = 4579 s). Blue squares: using the parameter values found by the analysis of the spectrum (see Table 2).
spacing pattern is not affected, it should be made clear if this impacts upon our estimates of ν rot and P 0 . Model B was computed without turbulent diffusion mixing so that a sharp chemical composition gradient grows at the convective core boundary forming a buoyancy glitch. We set the rotation frequency of this model to 7 µHz to compute the mode periods. Our results (Table 4 and middle left panel of Fig. 5) show significant disagreements between the three dipole mode series. In addition, substantial deviations from the true parameter values are found. These are one order greater than for the smooth model used in Sect. 3.2, reaching up to a maximum relative difference of ∼17% for ν rot and ∼6% for P 0 . The middle right panel of Fig. 5 displays the stretched periodéchelle diagram for prograde modes where the oscillatory component due to the buoyancy glitch is clearly visible. Otheréchelle diagrams show identical behaviour.
In the general case, the oscillatory component, which is superposed to the smooth period spacing pattern in the presence of a buoyancy glitch, is not symmetrical. By using the location of the maximum of PSD as an estimator of P 0 , we actually find the mean spacing of the stretched pattern. The systematic errors Notes. The non-asymptotic TAR and the complete treatment of rotation (ACOR) are compared for three different rotation rates. (a) For this specific model, an accumulation of higher modes occurs around the pulsation period of the radial order n = −45, which causes an avoided crossing. underlined here can be, for the most part, attributed to the difference between this mean spacing and the true value of P 0 . We note, however, that the biases found in Sect. 3.2 also contribute here.

Differential rotation
The TAR is a simplified treatment of rotation. One of its major and limiting hypotheses is the assumption of solid-body rotation, which is more invoked for the sake of mathematical simplification rather than from physical considerations. As differentially rotating stars will almost surely be encountered in observational data, it is interesting to test if we are able to find regularities in the stretched spectrum for such stars and, if so, investigate the impact of differential rotation on the estimate of the near-core rotation rate and buoyancy radius as obtained with our method. With this aim, we modelled a synthetic spectrum from Model A (Table 1) on which we applied a two-zone rotation profile. The convective core rotates at ν rot,core = 15 µHz while the rotation rate is ν rot,env = 7 µHz in the envelope. In Model A, the gradient of chemical composition at the convective core boundary is smoothed, which allows us to dismiss the effect of a buoyancy glitch hereafter. The bottom left panel of Fig. 5 shows the contours for each DFT map. Table 5 lists the results of the analysis. Our results suggest that we should be in a position to find regularities in the stretched spectrum of a differentially rotating star. Moreover, the values of ν rot and P 0 obtained from each mode series differ from each other at a substantial level. These differences are superior to the expected biases and errors in an equivalent model in solid-body rotation (Model A and B, see Sects. 3.2 and 3.3). The bottom right panel of Fig. 5 represents the stretched period echelle diagram in the case of the prograde dipole modes. The ridge is almost vertical but also slightly curved. Moreover, two avoided crossings occur at ∼15.5 and ∼17.5 µHz due to an accumulation of higher degree modes ( = 5 and = 7, respectively) at these frequencies. This illustrates the extent of the TAR limitations.
In a rotating star, the resonant cavity of gravity modes varies from one pulsation mode to another, which in turn impacts on their properties like the pulsation period. While this variation is small (but still noticeable) within modes of the same ( , m), it can be quite large between modes of distinct ( , m) couples. Then, we have a signature of differential rotation. In particular, prograde modes probe deeper into the star than zonal and retrograde modes, which is consistent with much faster rotation. This is less clear for zonal and retrograde modes as the two rotation frequencies obtained are quite similar and this slight difference could indeed be explained by other factors. In particular, the corresponding ridges are not vertical in the stretched periodéchelle diagram. For this same reason, it is also difficult to reliably interpret the estimate of the buoyancy radii.

Applications on Kepler targets
As a proof of concept, we applied our method to two stars targeted by Kepler, the γ Dor star KIC12066947 and the SPB star KIC3459297. Both were previously studied in the literature hence giving us a point of comparison with other existing methods.

Frequency analysis
For the two stars, the frequency analysis was performed using the classical iterative pre-whitening process, where, at each step, the peak with the highest amplitude in the periodogram was subtracted from the light curve. The statistical significance of each peak was derived on the basis of the false alarm probability (Scargle 1982), which gives robust results for Kepler data. To make sure that the peak extracted was not introduced during the pre-whitening process, we compare the amplitude of the extracted peak with the value in the original data. If these deviate more than 25%, we disregard this peak. A similar procedure was also used by Van Reeth et al. (2015a). A detailed description of the code used here will be given in Antoci et al. (in prep.).
Despite these precautions, noise peaks or spurious frequencies created by the pre-whitening process may still be extracted. Robustness of the stretching method is ensured as long as these frequencies remain in the minority compared to the modes that are looked for, which justifies being conservative in the frequency analysis. On the contrary, actual pulsation mode frequencies of a series may not be extracted because of overly cautious criteria. This would introduce windowing in the DFT map or, in the worst case, prevent us from identifying the mode series. In this case, accepting frequencies with a worst agreement between extracted and original amplitude may improve the detection of the regularity. We tested the effect of missing orders in a modes series with simulated spectra. We found that a regularity could be detected as long as more than approximately half of the modes were present in the frequency list.
Line-of-sight motions of stars introduce Doppler shifts of the pulsation frequencies, which can be larger than their associated uncertainties in Kepler data (Davies et al. 2014). Here, we neglected these Doppler frequency shifts as we could not find Left panels: contour maps at 95% (solid) and 50% (dashed) of the maximum of power spectral density for Model A in solid-body rotation (top), Model B in solid-body rotation (middle), and Model A in differential rotation (bottom). Colours are indicative of the type of modes ( , m) on which was applied the method. Red dots indicates the input value of ν rot and 1/P 0 as measured from the models, when relevant. Right panels: example of stretched periodéchelle diagrams for prograde modes, which are plotted twice for clarity. line-of-sight velocity measurements in the literature for these two stars. The impact on the measured ν rot and P 0 is expected to be contained within our internal uncertainties, even assuming large radial velocities (±150 km s −1 ).

The γ Dor star KIC12066947
The γ Dor star KIC12066947 was observed by the Kepler spacecraft in long-cadence mode during quarters 0-1 and 10-17 thus gathering almost 670 days of high-quality photometric data. We could extract 22 peak frequencies from the Lomb-Scargle periodogram (Table A.1). As shown in Fig. 6, most extracted peaks cluster in two groups in frequency roughly between 20 and 23 µHz (hereafter referred as Cluster A), and between 27 and 35 µHz (Cluster B). We analysed each cluster of peaks independently ensuring that we explored a sufficiently broad parameter space. The analysis of Cluster B frequencies leads to a detection for several trial values of ( , m). However, unless we suppose they are (1,1) modes, the buoyancy radius given by the maximum of PSD is much larger than generally expected values for γ Dor stars. We therefore identify these peaks as prograde dipole modes. The left panel of Fig. 7 shows the DFT map obtained in this case. We found ν rot = 24.95 ± 0.05 µHz (2.156 ± 0.004 c/d) and P 0 = 4181 ± 63 s (0.0484 ± 0.0008 d). The quoted uncertainties were computed from a Monte-Carlo simulation as explained in Sect. 2.5 and do not account for systematic errors. The stretchedéchelle diagram built using these values is shown on the right panel of Fig. 7. The ridge shows important curvatures indicating a substantial deviation from the asymptotic TAR. Such a feature is, however, difficult to explain without a proper detailed study of the star. Furthermore, there is a great gap between the frequency at ∼27.8 µHz and the rest of the ridge, causing a severe window effect in the DFT map. We checked if this has any influence on the derived values of ν rot and P 0 by discarding the frequency and re-analysing this part of the spectrum. In this manner, we derived ν rot = 24.94 ± 0.04 µHz and P 0 = 4157 ± 51 s in good agreement with the previously determined values.
Despite searching for modes up to = 4 and trying rotation frequencies from 2 to 35 µHz, we could not identify Cluster A peaks. Given their high amplitudes, these peaks are more likely to be dipole modes. Using Eq. (5) and the parameter values derived from the analysis of Cluster B peaks, we computed the pulsation periods of dipole retrograde and zonal modes within the asymptotic TAR and then checked if they were consistent with the range of periods observed. Unfortunately, they do not agree well suggesting that these peaks do not originate from gravito-inertial mode pulsations.
The oscillation spectrum shows an excess of power around 55-65 µHz with peaks of very low amplitude that are not extracted. Although they are located in the expected frequency range for (2,2) modes, the majority of these peak frequencies are simple combinations (ν i + ν j ) or first harmonics (2ν i ) of the highest amplitude peaks.
KIC12066947 was studied by Van Reeth et al. (2015b as part of their sample. Using the method of Van Reeth et al. (2015a) to analyse the oscillation spectrum, they detected two non-equidistant period spacing patterns. The first mode series consists of pulsation modes with frequencies in the range 27-37 µHz and a downwards slope in the ∆P − P plane. Modelling this pattern, they identified them as (1, 1) pulsation modes and found ν rot = 25.00 ± 0.10 µHz (2.160 ± 0.008 c/d). From their best-fit model, they could also estimate the asymptotic period spacing ∆P = 2950 ± 70 s, which translates into P 0 = 4171 ± 99 s (0.0484 ± 0.0012 d) using Eq. (1). These results agree very well with our determination from Cluster B peaks. The second pattern has an upward slope consisting of pulsation frequencies in the range 19-22 µHz, corresponding to our Cluster A. The authors could not identify them as g modes. The upward slope in the ∆P − P diagram suggest their retrograde character, however, assuming the found ν rot is accurate, their pulsation periods are much shorter than what would be expected from retrograde g modes. On the other hand, the authors were able to show that the period spacing pattern is consistent with retrograde Rossby modes (or r modes) using the asymptotic expression of Townsend (2003). This seems also supported by theory and numerical computations that predict the excitation of r modes in fast-rotating stars due to an interaction between toroidal motion and rotation (Townsend 2003;Salmon et al. 2014;Saio et al. 2018). Interestingly, this would also explain why we could not identify these modes with our method. The TAR also predicts the r-mode periods with the same expression as Eq. (5). Our stretching method can then also be applied to these modes by providing the λ functions. First attempts are encouraging.

The SPB star KIC3459297
The SPB star KIC3459297 was observed by Kepler in long cadence mode during quarters Q0-5, 7-9, 11-13, and 15-17, which represents 1069 days of high-quality photometric data. The Lomb-Scargle periodogram computed from the light curve is plotted in Fig. 8. We extracted 27 peak frequencies, of which six are found to be combinations and are discarded for the rest of the analysis (Table A.2). Two groups of peaks clearly stand out. The cluster at very low frequency is most likely due to instrumental effects or may have been introduced by the detrending of the light curve. We analysed the other group in the 9-14 µHz interval.
From a first analysis, we could detect a series of modes, but due to numerous missing orders and consequent windowing, we could not reliably obtain ν rot and P 0 . In the following, we therefore included peaks with a lower agreement between the amplitude extracted and that of the original data (50%). The 14 independent frequencies added are indicated by yellow bars on Fig. 8. The detection threshold was reached for several trial ( , m) values. However, similarly to the case of KIC12066947, only (1, 1) gives a buoyancy radius consistent with expected values for SPB stars. The peaks are thereby identified as prograde dipole modes. We estimated ν rot = 6.85 ± 0.07 µHz and P 0 = 7018 ± 190 s (see left panel of Fig. 9). To obtain this estimate, we discarded the two isolated peaks at 9.25 and 9.48 µHz to avoid the effect of windowing. Furthermore, as can be seen on the stretched periodéchelle diagram (right panel of Fig. 9), the belonging of these two frequencies to the prograde series is debatable. Remarkably, the period spacing pattern shows an oscillatory behaviour, which is a typical signature of a buoyancy glitch (Miglio et al. 2008). The power excess between 23 and 28 µHz is dominated by simple combination frequencies of high amplitude peaks (ν i + ν j ). The couple of remaining frequencies may be actual (2, 2) modes but cannot be formally identified with our method due to their restricted number.
This SPB star was studied in Pápics et al. (2017). The authors found a long period spacing pattern consisting of 43 modes in the 9-14 µHz domain with a downward slope in the ∆P − P plane. Following the method of Van Reeth et al. (2016), they identified the mode series as dipole prograde modes. Moreover, they estimated ν rot = 7.3 ± 0.5 µHz (0.63 ± 0.04 c/d) and ∆P = 5840 +950 −860 s, which translates into a buoyancy radius of P 0 = 8260 +1400 −1300 s (0.096 +0.016 −0.014 d). These findings are compatible at a 1-σ level with ours although are somewhat greater. We made the choice to stay more conservative in the frequency analysis of the light curve. This results in a lesser number of extracted frequencies than Pápics et al. (2017). We analysed the list of prograde mode periods provided in their Table 6 using our method. The results (ν rot = 6.81 ± 0.02 µHz, P 0 = 6857 ± 66 s) are not significantly different from our first estimate showing that the difference does not come from the period list. Pápics et al. (2017) found a second series of eight modes in the 25-32 µHz range that they do not identify as combination frequencies. They suggested these are (2,2) modes.

Discussion and conclusions
We have developed and tested a stellar-model-independent method to disentangle the g-mode spectrum of moderately to rapidly rotating γ Dor and SPB stars. Using this method, we are able to simultaneously obtain the mode identity and estimate the near-core rotation frequency and buoyancy radius of these stars. Moreover, we successfully apply this method to two Kepler targets. Our results are compatible with preceding determinations that used a different method (Van Reeth et al. 2016;Pápics et al. 2017), but the uncertainties on ν rot and P 0 intrinsic to our method are substantially lower. It should be noted, however, that the dominant source of errors comes from the inadequacies of the asymptotic TAR on which Van Reeth et al. (2016) also relied. It is also expected that the non-asymptotic formulation that was used by Pápics et al. (2017)  biases. While the method works well in most cases, precautions need to be taken during the analysis or in the interpretation of results in some specific cases. Because we use the DFT to detect regularities in the stretched spectrum, modes in a given series must be in sufficient number so that a peak in the DFT map can be unambiguously attributed to a regularity, that is, with an amplitude that is above the detection threshold (see Sect. 2.4). We are also subject to the effects of windowing especially in the case of missing orders in a ridge. Therefore, particular attention should be paid to the results when signatures of these effects are visible in the DFT map.
As a result of the non-linear coupling of pulsation modes, combination frequencies are common in the oscillation spectrum of γ Dor and SPB stars. These combinations are generally considered to be harmonics or simple combinations of the frequencies of highest amplitude but this has been subject to debate as the underlying physical mechanism is not understood (see Kurtz et al. 2015, and references therein). Possible combination frequencies must be identified. The analysis of combination frequencies, as if there were independent frequencies, can lead to wrong interpretations. For instance, the first harmonics of (1, 1) may be misidentified as (2, 2) modes with the same rotation rate and buoyancy radius that can be derived from the fundamental frequencies. In this work, to avoid these misinterpretations, we searched for frequencies that could be explained as linear combinations ( c i ν i where c i ∈ Z) and we discarded those of order 3 or less ( |c i | ≤ 3) from the analyses.
Numerous γ Dor and SPB stars are in the Kepler field of view. Few of them are slow rotators with rotational splittings that can be accurately interpreted with classical perturbative treatments of rotation. In this regard, we emphasise that the method presented in this paper complements the perturbative approach by allowing the interpretation of the spectrum of moderate to fast rotators. Thus, for rapid rotators, the analysis is straightforward as each mode series ( , m) is well separated from the others. In moderate rotators, several g-mode series may overlap, which slightly complicates the identification process as the stretching function is relevant only for a given ( , m) couple. This can be easily overcome by using the expected distribution of g modes in rotating stars as an indicator to guide the analysis (see Sect. 2.2). As long as some ( , m) modes are dominant in number in the portion of the spectrum studied, the associated regularity should be detectable. Then, proceeding iteratively (identifying the dominant mode series, removing these modes from the list of frequencies, and re-analysing the remaining frequencies) allows the identification of overlapping series.
In the current implementation of the method, we assess the rotation-pulsation coupling in the framework of the asymptotic TAR. While the hypotheses behind this approximation are rather crude, it seems to offer reasonably accurate results. For a typical model of a γ Dor star in uniform rotation and no structural glitch, we evaluated that the derived rotation rate and buoyancy radius are biased by only a few per cent. The TAR and the asymptotic approximations contribute to this bias in approximately equal proportions. In an equivalent model with a buoyancy glitch, the systematic errors can reach up to ∼17% for ν rot and 6% for P 0 because the effects of the buoyancy glitch on the pulsation periods are not accounted for in the asymptotic approximation. For most observed dipole prograde modes, these biases remain limited to 5% for both ν rot and P 0 . In a differentially rotating model, we show a clear signature of differential rotation although the full interpretation of the g-mode spectrum is not straightforward. The TAR assumes solid-body rotation and describes the pulsations in the co-rotating frame, which cannot be properly defined when there is differential rotation. Secondly, as mentioned before, the g-mode cavity changes from mode to mode. This questions the definition of a unique rotation frequency for a given mode series.
We stress that the limitations and potential biases shown here are not specifically linked to our stretching method but are inherited from the asymptotic TAR, thus any method based on it would suffer the same issues. Furthermore, such a method can be adapted to better prescriptions for pulsations in rotating stars. As long as the expression of pulsation periods adopts a suitable form, it is possible to define a stretching function such that the pulsation periods are regularly spaced once stretched. In which case, the determination of ν rot and P 0 can be refined without modifying the core principle of the method. Recent theoretical studies (Prat et al. 2016(Prat et al. , 2017 have taken an interest in developing new and more accurate asymptotic theories on the basis of ray theory, which provides an interesting lead for future improvements of this method.
Kepler provided high-quality photometric data for large numbers of γ Dor and SPB stars. Thanks to recent progress in the theory of pulsations in rotating stars as well as in the interpretation of their oscillation spectra, we are now in a good position to exploit the Kepler data to their full potential. That promises to put tight constraints on convective-core boundary mixing and angular momentum transport.