© The Authors 2024

1. Introduction

Asteroseismology has proven to be a crucial tool in unravelling the mysteries of stellar evolution and determining the ages and evolutionary stages of stars (e.g. Lebreton et al. 2014; Aerts 2021, and references therein). Estimating the evolutionary stages of stars is particularly significant to addressing the complex question of angular momentum transport within stars (Aerts et al. 2019).

Among pulsating stars, γ Doradus (hereafter γ Dor) stars are main-sequence stars that typically have a mass of 1.4–2.0M_⊙ and a radius of 1.4–2.3R_⊙ (Kaye et al. 1999; Uytterhoeven et al. 2011). Their structure is characterised by a convective core, a large radiative zone, and a small sub-surface convective envelope. They exhibit non-radial gravity modes oscillating with a period of 0.3 to 3d (Kaye et al. 1999). Thanks to high-precision photometry provided by space-borne instruments, especially by the Kepler mission (Borucki et al. 2010) and the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014), a large quantity of good quality data is now available, and it has enabled the study of numerous γ Dor stars.

Gravito-inertial modes have been identified in 611 γ Dor stars observed by Kepler (Li et al. 2020) and in 106 γ Dor stars observed by TESS (Garcia et al. 2022). These observed modes have enabled both the rotation rate near the bottom of the radiative zone and the buoyancy radius for a large number of γ Dor stars to be determined (see Van Reeth et al. 2015, 2016, 2018; Ouazzani et al. 2017, 2019; Saio et al. 2018; Christophe et al. 2018; Mombarg et al. 2019; Li et al. 2019, 2020; Takata et al. 2020), which has led to unprecedented constraints on the models of angular momentum transport (Ouazzani et al. 2019; Li et al. 2020). These seismic analyses used an approximation of the rotation effects on gravity modes called the traditional approximation of rotation (TAR; see Unno et al. 1989; Lee & Saio 1989, 1997; Townsend 2003; Bouabid et al. 2013). While generally convenient for modes confined in radiative zones (Ballot et al. 2012; Ouazzani et al. 2017, 2020), the TAR cannot describe modes oscillating in the convective cores of γ Dor stars.

Using an oscillation code that includes a full treatment of the Coriolis force, Ouazzani et al. (2020) showed the existence of gravito-inertial modes resonantly interacting with pure inertial modes confined in the convective core. In the otherwise smooth curve representing the period spacing between modes of consecutive order ΔP as a function of the period P, this interaction produces observable dips localised near the period of the inertial mode. This discovery constitutes a major breakthrough since such a phenomenon can help probe the convective core of γ Dor stars. Saio et al. (2021) analysed dips in 16 γ Dor stars observed by Kepler and estimated their core rotation rate. Moreover, Tokuno & Takata (2022) proposed a theoretical model of the shape of resonant dips in the ΔP–P curve.

In this article, we numerically study the resonance between gravito-inertial and pure inertial modes of the convective core with the aim of extracting as much information as possible on the convective core from seismic data. For the different series of observed gravito-inertial modes (e.g. Li et al. 2020), we studied the evolution of the dips with the rotation rate and the evolutionary stage of the star. This led us to develop seismic diagnostic tools that are usable on observational data. We also tested the model of Tokuno & Takata (2022) against our numerical results, which led us to construct a new model that better compares with these numerical results.

In the following, we first present the stellar models and the oscillation code we used and discuss the mode identification in Sect. 2. We then present the results of our numerical computations for three series of modes at different rotation rates and evolutionary stages in Sect. 3. In Sect. 4, we provide an empirical model that relates the shape and position of dips to the inner structure of the stars. We discuss the model of Tokuno & Takata (2022) and present our new model in Sect. 5. We then move to a discussion and the conclusion of paper in Sect. 6.

2. Method

In this section we describe the oscillation code and the stellar models we used to study properties of the dip. We then explain how oscillation modes are identified and selected, and finally we treat the question of numerical resolution.

2.1. Stellar models

In this subsection, we present the three main-sequence star models we used for our study. They were computed using the CESAM code (Morel 1997; Morel & Lebreton 2008). Their main characteristics are shown in Table 1. They are named after the CLES models used in Ouazzani et al. (2020), as they share similar characteristics. For these models, we adopted a solar metal mixture (Asplund et al. 2009), an initial helium fraction Y = 0.27, and a initial metallicity Z = 0.02. We used opacity tables from OPAL (Iglesias & Rogers 1996) completed at low temperature with tables from Ferguson et al. (2005). We used the OPAL2005 equation of state (Rogers & Nayfonov 2002) and the nuclear reaction rates from the NACRE collaboration (Angulo et al. 1999) except for the ¹⁴N(p,γ)¹⁵ O reaction, for which we used the LUNA reaction rate given in Imbriani et al. (2004). Convection was treated using the mixing-length theory (Böhm-Vitense 1958) with a mixing-length parameter α_MLT = 1.70, close to a solar calibration. We included turbulent diffusion with a constant diffusion coefficient of D_t = 700 cm² s⁻¹. The physical prescriptions adopted for these models are thus very similar to the ones used in Ouazzani et al. (2020).

The three main-sequence models seek to cover the instability strip of γ Dor stars (Bouabid et al. 2013). The model labelled ‘1z’ describes a young star of 1.40 M_⊙ close to the zero age main sequence (ZAMS; age of 180 My from the ZAMS, X_c = 0.68, with X_c being the mass fraction of hydrogen in the core). The model labelled ‘2m’ describes an evolved mid-main-sequence star of 1.60 M_⊙ (age of 1800 My from the ZAMS, X_c = 0.38). The model labelled ‘3t’ describes a star of 1.86 M_⊙ at the end of the main sequence (age of 1480 My from the ZAMS, X_c = 0.05).

The Brunt-Väisälä frequency N is an important parameter in our study since gravity modes depend directly on it. Moreover, the resonance with inertial modes in the convective core is greatly affected by the behaviour of N at the bottom of the radiative zone, as we demonstrate later. As shown in Fig. 1, the three models have different radial profiles of the Brunt-Väisälä frequency. They all present a discontinuity at the convective core-radiative zone interface, even the model that has just left the ZAMS (model 1z). This is due to the chemical composition gradient that develops on the outer edge of the convective core as the star evolves on the main sequence. This gradient, and thus the jump of the Brunt-Väisälä frequency, is sensitive to mixing processes such as rotational mixing, overshooting and microscopic diffusion. These mixing processes are modelled here through the turbulent diffusivity D_t mentioned above.

Fig. 1. Squared Brunt-Väisälä frequency N2 of the models presented in Table 1 as a function of the relative radius. Model 1z is drawn in blue, 2m in orange, and 3t in red.

2.2. Oscillation code

To calculate the oscillation modes and frequencies, we used the code TOP (Two-dimensional Oscillation Program, Reese 2006; Reese et al. 2009), which we describe in this subsection.

We considered time-harmonic adiabatic small perturbations of a uniformly rotating spherical star. Thanks to the axial symmetry, the solutions are proportional to e^{i(mϕ + ωt)}, where m is the azimuthal number and ω is the mode frequency. The governing equations are

$$\begin{array}{cc}& \rho +\mathit{\xi }·\mathbf{\nabla }{\rho }_0+{\rho }_0\mathbf{\nabla }·\mathit{\xi }=0,\hfill \end{array}$$(1)

$$\begin{array}{cc}& 2i{\rho }_0\mathbf{\Omega }\times [\omega +m\mathrm{\Omega }]\mathit{\xi }-{[\omega +m\mathrm{\Omega }]}^2{\rho }_0\mathit{\xi }=-\mathbf{\nabla }p+\frac{\mathbf{\nabla }{P}_0}{{\rho }_0}\rho -{\rho }_0\mathbf{\nabla }\psi ,\hfill \end{array}$$(2)

$$\begin{array}{cc}& p-c_0^2\rho +\mathit{\xi }·(\mathbf{\nabla }{P}_0-c_0^2\mathbf{\nabla }{\rho }_0)=0,\hfill \end{array}$$(3)

$$\begin{array}{cc}& \mathbf{\Delta }\psi -4\pi G\rho =0,\hfill \end{array}$$(4)

where ξ, p, ρ, and ψ are the amplitude of the Lagrangian displacement, the Eulerian perturbations of pressure, density, and gravitational potential, respectively, and where ρ₀, P₀, and c₀ are the density, the pressure, and the sound velocity of the star model, respectively, with $c_0=\sqrt{{\mathrm{\Gamma }}_1{p}_0/{\rho }_0}$, Γ₁ being the first adiabatic index. The star rotation rate is Ω, and G is the gravitational constant. The mode frequency in the inertial frame ω is related to the mode frequency in the co-rotating frame ω_co through

$$\begin{array}{c}\hfill \omega ={\omega }_{\mathrm{co}}-m\mathrm{\Omega }.\end{array}$$(5)

The governing equations were discretised using a fourth-order finite difference method in the radial direction (for details on the scheme, see Reese 2013) and spectral decomposition on spherical harmonics in the horizontal directions. After discretisation, an algebraic eigenvalue problem was obtained and solved using the Arnoldi-Chebyshev algorithm (Chatelin 1988; Reese et al. 2009). Due the axial and equatorial symmetries of the problem, independent eigenvalue problems were solved for a given azimuthal order m and a given parity with respect to the equator.

2.3. Mode identification

In practice, TOP finds a certain number of modes (typically two to eight) around a given frequency for a prescribed resolution. These modes have the same equatorial parity and azimuthal order, but modes associated with different degrees and radial orders are computed simultaneously. Unlike with 1D calculations, identification is thus non-trivial. Additionally, modes unresolved at the chosen resolution can be among the computed ones.

The series of gravito-inertial modes we computed for this work are the modes with the same degree ℓ and azimuthal number m but different radial orders n. The degree ℓ corresponds to the associated Legendre polynomial $P_{\mathcal{l}}^m$ that characterises the mode at zero rotation but remains relevant to label these same modes at higher rotation rates. For identification, we used the fact that the modes of the series are recognisable by their radial and latitudinal structure. We applied three criteria on the perturbed pressure p to select the desired modes among the computed ones; two are related to the radial and latitudinal profiles, while the third criterion takes care of miscalculated modes. Using the asymptotic Wentzel-Kramers-Brillouin (WKB) formulation of the TAR (see Unno et al. 1989), the number of radial nodes n of a (ℓ, m) gravito-inertial mode can be estimated by

$$\begin{array}{c}\hfill n+{ϵ}_g\approx \frac{P_{\mathrm{co}}\sqrt{{\mathrm{\Lambda }}_{\ell }^m}}{{\mathrm{\Pi }}_0},\end{array}$$(6)

with P_co as the period of the mode in the co-rotating frame, Π₀ as the buoyancy radius, ${\mathrm{\Lambda }}_{\mathcal{l}}^m$ as the eigenvalue associated with the Hough function $H_{\mathcal{l}}^m$ describing the latitudinal part of the mode in the TAR (see Appendix D.2), and ϵ_g as a small offset. We used this estimate of n to remove modes with radial orders in the radiative zone noticeably higher than expected. Additionally, we computed the correlation of the mode latitudinal profiles in the radiative zone with the corresponding Hough function and rejected modes with weak correlation coefficients. Finally, we excluded badly resolved and spurious modes. This included modes with an extremely high amplitude in the core or at the surface as well as modes showing strong variations between consecutive radial grid points.

2.4. Numerical resolution

In this subsection, we estimate which radial and latitudinal resolutions are needed for our analysis. The radial resolution was determined by the number of radial points n_r and needed to be adjusted depending on the number of radial nodes of the considered mode.

The latitudinal resolution n_θ corresponds to the number of associated Legendre polynomials $P_{\mathcal{l}}^m$ on which the modes are decomposed. The degree ℓ spans from ℓ_min = |m|+i_p to ℓ_max = ℓ_min + 2(n_θ − 1), where i_p = 0 or 1 when ℓ + m is even or odd, respectively. For example, a mode of the series (ℓ = 1, m = −1) calculated with n_θ = 3 is decomposed latitudinally on the three first associated Legendre polynomials with the same azimuthal number and equatorial parity $P_1^{-1}$, $P_3^{-1}$ and $P_5^{-1}$. The minimal latitudinal resolution needed to compute a given mode can be estimated by looking for the number of associated Legendre polynomials needed to represent the associated Hough function. Figure 2 shows that for the Hough function $H_{\ell =1}^{m=-1}$ calculated at the spin parameter s = 2Ω/ω_co ∼ 16.4, the projection on Legendre polynomials peaks at $P_1^{-1}$, while the contributions of polynomials above $P_{11}^{-1}$ are very small. This suggests that a latitudinal resolution of at least n_θ = 6 is needed for this mode.

Fig. 2. Projection of the Hough function describing the latitudinal profile of the mode n = 66 (ℓ = 1, m = −1) in the radiative zone (TAR) on the associated Legendre polynomials$P_{{\mathcal{l}}_i}^{m=-1}$.

For a more accurate estimation, we applied the same process with modes calculated by TOP. For example, in the upper panel of Fig. 3, one can see the decomposition of the mean latitudinal profile of a mode of the series (ℓ = 1, m = −1) onto the associated Legendre polynomials for different latitudinal resolutions. We first noticed that the amplitude decreases rapidly with the degree ℓ_i. Nevertheless, when we used too few polynomials, we missed non-negligible contributions of polynomials of degree ℓ_i larger than ℓ_max, leading to inaccuracies. This translated into errors on the frequency determination. In the lower panel of the same figure, one can see how the computed frequency evolves with the latitudinal resolution and converges towards a stable value at a large ℓ_max. For this mode, the frequency varies by less than 0.01% for an ℓ_max higher than nine (n_θ = 5).

Fig. 3. Effect of the latitudinal resolution on the mode n = 66 (s ∼ 16.4) of the series (ℓ = 1, m = −1). Upper panel: Projections of the mean latitudinal profile in the radiative zone of this mode at different latitudinal resolution nθ on the associated Legendre polynomials for a rotation rate Ω = 0.4ΩK. Here, ${\mathrm{\Omega }}_{\mathrm{K}}=\sqrt{GM/R^3}$ is the critical rotation rate of the star and M and R are its mass and radius, respectively. Lower panel: Evolution of the relative frequency of the same mode as a function of ℓmax. The reference frequency ωref corresponds to the frequency obtained using the highest latitudinal resolution (nθ = 12). The model used is 1z.

Because we studied the dips, the resolution was chosen such that the dips are computed with enough accuracy. Figures 4 and 5 show the impact of the radial and latitudinal resolutions on a series of modes in diagrams depicting period spacing ΔP_co (period differences of two modes with consecutive n in the co-rotating frame) as a function of the spin parameter s for the rotation rate Ω = 0.4Ω_K. The vertical line indicates the spin parameter of the expected resonant mode in the case of a core of uniform density (see Sect. 3.2). We observed that the dip deepens and that ΔP_co increases for the higher s as the radial resolution decreases, which increases the general slope of the curve. When the latitudinal resolution decreases, ΔP_co decreases for the higher s, which leads to the opposite change in the general slope of the curve. The curve barely evolves for n_r > 1000 and n_θ > 6, which confirms the primary guess we made using the projection of the Hough function on the associated Legendre polynomials. Based on Fig. 4 and knowing that the mode with the highest number of radial nodes is n = 72, we deduced that about 15 × n points is enough to reach convergence.

Fig. 4. Period spacing ΔPco, normalised by Π0 (buoyancy radius), plotted as a function of the spin parameter s with different radial resolution nr for the series of modes (ℓ = 1, m = −1) and for a rotation rate Ω = 0.4ΩK. Their radial order n goes from 32 to 72. The model used is 1z. The pink line shows the traditional approximation of rotation. The vertical line shows the spin parameter of the resonant inertial mode calculated with an analytical uniform density model.

Fig. 5. Period spacing ΔPco/Π0 plotted as a function of the spin parameter s with different latitudinal resolutions nθ for the series of modes (ℓ = 1, m = −1) and for a rotation rate Ω = 0.4ΩK. Their radial order n goes from 32 to 72. The model used is 1z. The pink line shows the traditional approximation of rotation. The vertical line shows the spin parameter of the resonant inertial mode calculated with an analytical uniform density model.

To describe the dips accurately, we wanted the error on the period spacing δΔP_co not to exceed 1% of the total depth of the dip. We determined the radial and latitudinal resolutions for our numerical calculations in order to meet this criterion on the frequency error δω using

$$\begin{array}{c}\hfill \delta \omega ={\omega }_{\mathrm{co}}^2\frac{\delta \mathrm{\Delta }{P}_{\mathrm{co}}}{2\pi \sqrt{2}}.\end{array}$$(7)

To ensure that δω < 1%, we showed that it is sufficient to choose a radial resolution n_r ranging from 5000 to 10 000 for the (ℓ = 1, m = −1) series, from 3000 to 10 000 for the (ℓ = 2, m = −2) series, and from 3000 to 10 000 for the (ℓ = 1, m = 0) series. These ranges cover the different resolutions that are needed since the radial order increases with the spin parameter and evolutionary status. For the latitudinal resolution, we typically used n_θ = 6 to 8 for the series (ℓ = 1, m = −1) and (ℓ = 2, m = −2), and n_θ = 4 to 6 for the series (ℓ = 1, m = 0).

3. Results

In this section, we present the results of our numerical calculations. We studied the series of modes (ℓ = 1, m = −1), (ℓ = 2, m = −2), and (ℓ = 1, m = 0) at rotation rates ranging from 0.1 to 0.5 Ω/Ω_K, with ${\mathrm{\Omega }}_{\mathrm{K}}=\sqrt{GM/R^3}$, for the three stellar models described in Table 1. Among these three series, the (ℓ = 1, m = −1) is by far the most frequently observed (Li et al. 2020). For each series, we computed the period spacing ΔP_co in a frequency interval where the gravito-inertial modes strongly interact with an inertial mode of the convective core. We then described the properties of the dips obtained in the ΔP_co − s diagram. In the following, we first present the basic characteristics of a typical dip in Sect. 3.1. The simple model we used to determine the approximate spin parameter of the dips is described in Sect. 3.2. The results of our investigation of dips of the (ℓ = 1, m = −1), (ℓ = 2, m = −2), and (ℓ = 1, m = 0) mode series are then reported in Sect. 3.3.

3.1. Basic dip properties

We first describe the properties of a typical dip for the series (ℓ = 1, m = −1). The dip is located around s = 10.7 and was calculated with the 1z model at rotation rate Ω = 0.1Ω_K (shown in the left panel of Fig. 6). One can see the difference between our calculations and the TAR shown in pink. While there is a clear dip in the ΔP_co − s curve around s = 10.7 in our calculations, the TAR is nearly constant. For prograde sectoral modes (ℓ = −m), the TAR indeed leads to a constant ΔP_co in the limit of high spin parameters. The middle and right panels of Fig. 6 show the Eulerian pressure perturbation p for the mode located at the centre of the dip. It exhibits a very high amplitude in the convective core compared to the radiative zone. The presence of significant oscillations in both the convective core and the radiative zone indicates the mixed nature of the mode. To confirm the mixed nature of the modes in the dip, we calculated the portion of kinetic energy in the core for each mode. It reads

$$\begin{array}{c}\hfill \frac{E_{\mathrm{kin}}(\mathrm{core})}{E_{\mathrm{kin}}(\mathrm{core}+\mathrm{RZ})}=\frac{{\displaystyle {\int }_0^{r_c}}{\displaystyle {\int }_0^{2\pi }}{\displaystyle {\int }_0^{\pi }}{\rho }_0{\left|\mathit{v}\right|}^2{r}^{2}sin\theta \mathrm{d}\theta \mathrm{d}\varphi \mathrm{d}r}{{\displaystyle {\int }_0^R}{\displaystyle {\int }_0^{2\pi }}{\displaystyle {\int }_0^{\pi }}{\rho }_0{\left|\mathit{v}\right|}^2{r}^{2}sin\theta \mathrm{d}\theta \mathrm{d}\varphi \mathrm{d}r},\end{array}$$(8)

Fig. 6. Example of a resonant mode. Left panel: Period spacing ΔPco as a function of the spin parameter for the series (ℓ = 1, m = −1) and for the 1z model at a rotation rate of Ω = 0.1ΩK. Red dots represent computed modes. The pink line shows the TAR. The blue line shows the proportion of kinetic energy in the core for each mode. Middle panel: Quantity $p/{\rho }_0^{0.8}$ plotted in a meridional plane for the mode marked by a black star in the left panel. Right panel: Radial profile of $p/{\rho }_0^{0.8}$ at the latitude π/2 (equator). The factor ${\rho }_0^{0.8}$ is a scaling factor used to make the visualisation easier. The red vertical lines show the radiative zone limits.

with E_kin(core) being the mode kinetic energy in the core and E_kin(core + RZ) as the mode kinetic energy in the core and the radiative zone. The mode perturbation velocity is denoted as v, the radius of the core as r_c, and the radius of the star as R. As observed in the left panel of Fig. 6, the portion of kinetic energy in the core increases for the modes in the dip, reaching its maximum when ΔP_co is minimum. Dips are formed by mixed (inertial and gravito-inertial) modes around the spin parameter of an inertial mode of the core that couples with the gravito-inertial modes of the (ℓ,m) series. In the dip, a frequency is added to the series.

To study the dips and quantify their evolution, we modelled them as inverse Lorentzian profiles by following Tokuno & Takata (2022):

$$\begin{array}{c}\hfill \frac{1}{\mathrm{\Delta }{P}_{\mathrm{co}}}=x_1+x_{2}s+\frac{\mathrm{\Omega }\sigma /{\pi }^2}{{(s-s_c)}^2+{\sigma }^2},\end{array}$$(9)

where s is the spin parameter of the gravito-inertial modes in the co-rotating frame and x₁, x₂, s_c, and σ are free fitting parameters. With this parametrisation, x₁ and x₂ account for the general slope of the curve, s_c is the spin parameter on which the dip is centred, and σ is the width of the dip.

Adding one frequency locally modifies the period spacing but not the whole period gap between unaffected modes. This constraint imposes that the integral of the Lorenztian term over s is constant and thus that the width and the depth are anti-correlated (Tokuno & Takata 2022). We fit this model to our numerically computed dip using a least-squares method. The errors associated with the fit tend to increase with the rotation rate, as the dip implies less modes when the rotation increases, decreasing the quality of the fit. The relative error of s_c is negligible; the one of σ is always less than 1% for the series (ℓ = 1, m = −1) and less than 2% for the series (ℓ = 2, m = −2), but it goes up to 10% with the rotation rate for the series (ℓ = 1, m = 0). The latter case could be explained by the fact that x₁ and x₂ account for a linear slope of 1/ΔP_co, which is only a good local approximation of the TAR. As the rotation rate increases, the dip widens, and this approximation becomes less valid, explaining the growing relative error.

3.2. Model of Ouazzani et al. (2020) as an initial guess

Assuming a convective core of uniform density, Ouazzani et al. (2020) proposed a simple analytical model to identify the inertial modes that produce significant dips in the period spacings of a gravito-inertial mode series. Ouazzani et al. (2020) showed that this model provides a useful approximation of the dip spin parameters observed in full numerical computations.

We thus used this model to choose the dips we wanted to study and to get an initial guess of their spin parameter. The model details are given in Appendix A. Once a dip was found numerically, we could easily check that it is due to the expected inertial mode because the spatial structure of the inertial mode in our complete calculations (see for example the core region in the middle panel of Fig. 6) is very similar to that of the inertial mode in the sphere of uniform density (see comparisons in Appendix A).

For each series of gravito-inertial modes, we limited our study to resonances that could occur in spin-parameter ranges actually observed in γ Dor stars. This condition was verified for one resonance in each series, which happens to be the one occurring at the smallest spin parameter.

3.3. Evolution of the width and central spin parameter of the dips

We show in Figs. 7, 8, and 9 the evolution of the first dip of the series (ℓ = 1, m = −1) with the star rotation rate for the 1z, 2m, and 3t model, respectively. The square dots indicate modes within an observable range of radial orders. Indeed, according to Li et al. (2020), the observed modes of the (ℓ = 1, m = −1) series are typically of a radial order 30 < n < 70. Analogous figures for the other two series, (ℓ = 2 m = −2) and (ℓ = 1, m = 0), are shown in Appendix C.

Fig. 7. Period spacing in the co-rotating frame as a function of the spin parameter for the modes of the series (ℓ = 1, m = −1) at different rotation rates using the 1z model. The period spacing is normalised by the buoyancy radius Π0. The pink line shows the traditional approximation of rotation. The square dots indicate the modes that are within the range of typically observed radial orders 30 < n < 70. The dashed line corresponds to the frequency of the resonant inertial mode ℓi = 3, m = −1 calculated analytically with the uniform density model.

Fig. 8. Same as Fig. 7 but for the model 2m. The low-amplitude oscillations in ΔP are due to a glitch generated by a small numerical discontinuity in its Brunt-Väisälä frequency profile (see Miglio et al. 2008).

Fig. 9. Same as Fig. 7 but for the model 3t.

We fit the dips using Eq. (9) in order to determine their width σ and central spin parameter s_c. Figures 10 and 11 show the evolution of s_c and σ with the star rotation for the three stellar models. The label ‘core-only model’ refers to a model that we discuss later in Sect. 4.2 (some oscillation modes calculated with this model are shown in Appendix B). One can see in Fig. 10 that s_c varies with the rotation rate but becomes less sensitive to it as the star evolves. It generally increases with the rotation rate, with the exception of the 3t model between Ω = 0.1Ω_K and Ω = 0.15Ω_K. We also observed that s_c decreases as the star evolves. The width σ follows the same trend: It increases with the rotation rate but becomes smaller and less sensitive to the rotation rate as the star evolves (see Fig. 11). In summary, the rotation rate of the star tends to make the dips occur at a higher s and widens them, whereas the dips become narrower and occur at a lower s as the star ages.

Fig. 10. Evolution of the centre of the dip sc obtained using a full model (full line) and evolution of the spin parameter of the inertial mode obtained using a core-only model (dotted line). The dash-and-dot blue line shows the spin parameter of the inertial mode in the case of a uniform density core (analytical model). The light blue curves refer to the model 1z, orange curves to the model 2m, and red curves to the model 3t.

Fig. 11. Width of the dips σ as a function of the rotation rate for the series (ℓ = 1, m = −1) for the three stellar models described in Table 1. The error bars are too small to be visible with the scale used.

We note that outside of the dip and for the three models we considered, the curve remains further and further from the TAR as the rotation rate grows. This behaviour has already been observed in previous studies (e.g. Ouazzani et al. 2017).

4. Analysis

In this section we analyse the evolution on σ and s_c and derive approximate empirical relations usable on observational data.

4.1. Empirical expression of σ

According to the model of Tokuno & Takata (2022), σ evolves linearly with the parameter ϵ = Ω/N₀, where N₀ is the jump of the Brunt-Väisälä frequency at the interface between the convective core and the radiative region. We tested whether σ is proportional to Ω by plotting its evolution with Ω alongside the theoretical linear evolution predicted in Tokuno & Takata (2022). Figure 12 shows the evolution of σ/σ_ref with the rotation rate, where σ_ref is the width computed at Ω = Ω_ref = 0.2Ω_K. This figure shows that σ is indeed proportional to the rotation rate of the star, with the exception of the 3t model for which one can see a deviation that is significant compared to the errors of σ. The origin of this particular behaviour remains unclear at the moment, and numerical effects cannot be excluded. Analogous figures for the series (ℓ = 2, m = −2) and (ℓ = 1, m = 0) are shown in Appendix C.

Fig. 12. Evolution of σ/σref with the rotation rate for the dip of the series (ℓ = 1, m = −1) for three different evolutionary stages (1z, 2m, 3t). The associated error bars are too small to be visible. The reference value σref is the width for Ω = Ωref = 0.2ΩK. If, as predicted by Tokuno & Takata (2022), σ ∝ Ω, then σ/σref = Ω/Ωref = 5Ω (the dashed line).

Then, assuming a linear relation σ = βϵ, we performed a linear regression to determine the slope β. We used the least-squares method. Table 2 gathers the values of β obtained for each dip and stellar model. While β is expected to depend on the core structure (Tokuno & Takata 2022), we find that it varies by less than 15% for the three stellar models considered. For the dip of the (ℓ = 1, m = −1) series, we show in Fig. 13 the result of the fit obtained by considering the three stellar models together. When directly fitting with relation σ = βϵ, the 1z-model σ values are weighted more in the determination of the free parameter β because, due to the smaller N₀ of the 1z model, they extend over a larger ϵ range. To counter that, we actually fitted the relation σN₀ = βΩ. Analogous figures for the dips of the (ℓ = 2, m = −2) and (ℓ = 1, m = 0) series are shown in Appendix C. Using our empirical relation, the mean value of β could thus be used to estimate N₀ from σ via

$$\begin{array}{c}\hfill N_0^{\mathrm{calc}}=\frac{\mathrm{\Omega }\beta }{\sigma }.\end{array}$$(10)

Fig. 13. Width σ as a function of ϵ for the series (ℓ = 1, m = −1) for three different evolutionary stages (1z, 2m, 3t). The purple dashed line shows the fit presented in Table 2 for the three models together, and the purple area shows the related error. The x-axis is logarithmic for visualisation purposes.

4.2. Empirical expression of s_c

In order to understand the variations of s_c, we first computed the spin parameters s_* of inertial modes of truncated versions of our models limited to their convective cores, as done in Ouazzani et al. (2020). It allowed us to see the impact of varying the density as well as the effect of the rotation rate on the frequency of the resonant inertial modes. The modes were calculated with the boundary condition ξ_r = 0, with ξ_r being the radial displacement, and identified as the resonant modes using their number of latitudinal and radial nodes. The inertial modes and their spin parameter obtained for this core-only model are reported in Appendix B. As already observed in Ouazzani et al. (2020), the spin parameter of these modes remains independent of Ω as long as an acoustic term can be neglected in the governing equations (see Eq. (B.1)).

We could then compare the central spin parameter of the dip s_c obtained in the previous section with the spin parameter s_* of the resonant inertial mode of the core-only model. We show the comparison for the series (ℓ = 1, m = −1) for the three stellar models in Fig. 10. The spin parameter of the inertial modes obtained using the uniform density model presented in Appendix A is also indicated. Analogous figures for the (ℓ = 2, m = −2) and (ℓ = 1, m = 0) series are shown in Appendix C. We note that for the same inertial mode, the spin parameter decreases with the star evolution. This behaviour depends on the density stratification of the core, which is related to the stellar mass and the stellar evolution. We observed that while the spin parameters of the eigenmodes of the uniform-density model and the core-only model remain constant with the rotation (with the exception of s_* for the 3t model between Ω = 0.1Ω_K and Ω = 0.15Ω_K), the central spin parameter of the dip shows a significant evolution with the rotation rate. For the three models, s_c tends towards s_* as the rotation rate decreases. Additionally, the differences between the full model and the core-only model significantly decrease as the star ages. This is consistent with the fact that the Brunt-Väisälä jump N₀ grows as the star evolves on the main sequence, which makes the radial displacement ξ_r at the core interface increase and thus makes the core more isolated and closer to the core-only model. This also explains the decrease of the dip width with evolutionary stage, as a less deformable interface weakens the coupling, which thus impacts fewer and fewer modes and leads to narrower dips.

While looking for an empirical relationship between s_c and ϵ usable on observation data, we noticed that though s_c does not vary linearly with the rotation for all models (especially 3t), the difference s_c − s_* appears to be proportional to Ω (figure not shown). We thus fit s_c − s_* as a linear function of ϵ: s_c − s_* = αϵ. We estimated α for the series of modes (ℓ = 1, m = −1), (ℓ = 1, m = 0), and (ℓ = 2, m = −2) and for the models 1z, 2m, and 3t separately. Table 3 provides the estimate of α for each case.

One can see that α varies significantly (up to ∼30%) with the star model for a same series of modes. However, we observed that the ratio s_c/s_* plotted as a function of ϵ appears to almost gather on a same line for the three models (see Fig. 14). Since this quantity depends less on the model, we fitted it as a linear function s_c/s_* = Aϵ + B for the three models, separately and together. We obtained values of B compatible with 1, and this confirms that s_c tends to s_* when ϵ tends to 0. Therefore, we redid fits by fixing B = 1, letting only A be a free parameter. As we did for the estimation of β, we fitted (s_c/s_* − 1)N₀ = AΩ so that all the models were weighted equally. The result of the fit for A is shown in Table 4.

Fig. 14. Ratio sc/s* as a function of Ω/N0 (=ϵ) for the series (ℓ = 1, m = −1) for the 1z, 2m, and 3t models. The purple dashed line shows the fit presented in Table 4 for the three models together, and the purple area shows the related error. The x-axis is logarithmic for visualisation purposes.

The fitting parameter A varies by less than 15% for the (ℓ = 1, m = −1) and the (ℓ = 2, m = −2) series between the three models, while it goes up to about 30% for the series (ℓ = 1, m = 0).

The global fit of all models together is plotted in Fig. 14 for the series (ℓ = 1, m = −1). Analogous figures for the (ℓ = 2, m = −2) and (ℓ = 1, m = 0) series are shown in Appendix C. As A does not vary much with the stellar structure, we assumed it to be constant in order to retrieve the frequency of the resonant inertial mode in the core, s_*, from the location of an observed dip through the following relation:

$$\begin{array}{c}\hfill s_{\ast }^{\mathrm{calc}}=\frac{N_0^{\mathrm{calc}}{s}_c}{A\mathrm{\Omega }+N_0^{\mathrm{calc}}}.\end{array}$$(11)

Using Eq. (10), we rewrote this relation as

$$\begin{array}{c}\hfill s_{\ast }^{\mathrm{calc}}=\frac{\beta s_c}{A\sigma +\beta }.\end{array}$$(12)

Since s_c and σ are obtained from the ΔP_co–s relation in the co-rotating frame, the rotation rate Ω is needed. It can be determined by fitting the TAR outside the dip (e.g. Li et al. 2020).

4.3. Relative errors made on the estimations of N₀ and s_*

To estimate the errors made by our empirical method, we applied it to dips calculated numerically, and we obtained empirical measurements of N₀ and s_* using Eqs. (10) and (12), denoted as $N_0^{\text{calc}}$ and $s_{*}^{\text{calc}}$. We then calculated the errors relative to the real known values $N_0^{\text{true}}$ and $s_{*}^{\text{true}}$. We estimated s_c and σ by fitting the dips with the Eq. (9), and the rotation rates were assumed to be perfectly known. Figure 15 shows the relative error made on N₀ and s_* as a function of ${ϵ}^{\text{calc}}=\mathrm{\Omega }/N_0^{\text{calc}}$ for the series (ℓ = 1, m = −1). Analogous figures for the series (ℓ = 2, m = −2) and (ℓ = 1, m = 0) are shown in Appendix C. One can see that for the (ℓ = 1, m = −1) series, the relative error on N₀ is always inferior or equal to 10% for the model 1z and 2m, while it increases up to 50% for slowly rotating 3t models. The relative error for the 3t model is below 10% for the rotation rates Ω = 0.4, 0.5Ω_K, between 10 and 20% for the rotation rates Ω = 0.2, 0.3Ω_K, and between 30 and 50% for the rotation rates Ω = 0.1, 0.15Ω_K. Such behaviour is due to the direct dependency between $N_0^{\text{calc}}$ and σ. Indeed, as ϵ gets smaller, σ goes to zero, making relative errors explode when we evaluate the difference between σ and βϵ, despite decreasing absolute errors. Nevertheless, it has to be noted that as typical radial orders for observable modes of the (ℓ = 1, m = −1) series range between 30 and 70 (see Li et al. 2020), none of the dips obtained with the 3t model are observable. Indeed, they exhibit radial orders lying out of the [30,70] range, going from n ≈ 90 at Ω = 0.5Ω_K to n ≈ 650 at Ω = 0.1Ω_K. The same applies to the (ℓ = 2, m = −2) and (ℓ = 1, m = 0) series, as the observed modes are typically of about the same radial order as the (ℓ = 1, m = −1) series, but the radial order of the modes calculated with the 3t model never goes lower than n ≈ 170 for the (ℓ = 2, m = −2) series and n ≈  90 for the (ℓ = 1, m = 0) series.

Fig. 15. Relative errors for the series (ℓ = 1, m = −1) on the estimated$N_0^{\text{calc}}$ (upper panel) and $s_{*}^{\text{calc}}$ (lower panel) using the three models 1z, 2m, and 3t.

The relative error of s_* is far smaller and does not become greater than ∼0.7%. It tends to decrease as the stars evolves. As dips are not currently observed in stars as evolved as the 3t model, we also estimated the parameters A and β using only the 1z and 2m models (see the results in Table 5). The relative error made on N₀ and s_* for the series (ℓ = 1, m = −1) using these new estimates of A and β are shown in Fig. 16. The relative error of N₀ is always inferior to 6%, and the relative error on s_* is inferior to 0.2%, except for the case of 1z and Ω/Ω_K = 0.5, where it reaches ∼1%. The same approach was used for the series (ℓ = 1, m = 0) and (ℓ = 2, m = −2), and the results are reported in Appendix C.

Fig. 16. Relative errors for the series (ℓ = 1, m = −1) on the estimated $N_0^{\text{calc}}$ (upper panel) and $s_{*}^{\text{calc}}$ (lower panel) using the two models 1z and 2m.

We showed that for observable dips, our empirical method can accurately predict the frequency of the resonant inertial mode in the core, and it provides a reasonable estimate of the jump in the Brunt-Väisälä frequency at the interface between the core and the envelope. Both N₀ and s_* are correlated to the evolutionary stage of the star. In particular, the frequency of the resonant inertial mode is only sensitive to the structure of the core. It is thus independent of the mixing and diffusion processes occurring in the overlying radiative layers, which, in contrast, may deeply affect the shape of the Brunt-Väisälä frequency profile.

Further investigations are needed in order to characterise how the frequencies of inertial modes depend on the characteristics of the convective core. Inertial modes could be a powerful tool in obtaining constraints on convective cores and thus precise determinations of stellar ages.

5. Approximate analytical models

We present the dip model of Tokuno & Takata (2022) as well as a newly improved version of the model and compare them with the properties of the dip obtained by numerical calculations in the previous sections.

5.1. Model of Tokuno & Takata (2022)

Tokuno & Takata (2022) proposed an analytical model of the dips in ΔP = f(s) provoked by the coupling between gravito-inertial and inertial oscillations. It is based on three main assumptions: the density is uniform in the convective core, the gravito-inertial oscillations are described by the TAR, and the ratio between the (uniform) rotation rate and the Brunt-Väisälä frequency at the bottom of the radiative zone is small. The assumption of uniform density is obviously unrealistic, but it enables one to get analytical solutions of the inertial oscillations in the core. In the absence of coupling with the radiative envelope, inertial modes in the uniform density core are identified by their spin parameter, s_*, and the degree and azimuthal order, (ℓ_i, m), of the Legendre polynomial involved in the solution (see Appendix A). Under the TAR, the latitudinal part of gravito-inertial (ℓ,m) modes is characterised by the Hough function $H_k^m(\theta ,s)$ associated with the eigenvalue ${\mathrm{\Lambda }}_k^m(s)$, where k = ℓ − |m| (see Eq. D.10). The model results depend on the continuity of the density and the Brunt-Väisälä frequency at the interface between the convective core and the radiative zone. For all stellar models of γ Dor stars presented in Sect. 2.1, the density is continuous at the convective-radiative interface, while the Brunt-Väisälä is discontinuous. Indeed, the jump from zero on the convective side to N₀ on the radiative side increases with age. According to Tokuno & Takata (2022), the coupling between the (ℓ,m) series of the gravito-inertial modes and a (ℓ_i, m) inertial mode produces a Lorentzian-shaped dip in the ΔP = f(s) that depends on two parameters, namely its width¹

$$\begin{array}{c}\hfill \sigma =\frac{2P_{{\ell }_i}^m(1/s_{\ast })}{{\gamma \prime }_{{\ell }_i}^m(s_{\ast })}\sqrt{{\mathrm{\Lambda }}_k^m(s_{\ast })}\frac{\mathrm{\Omega }}{N_0},\end{array}$$(13)

and its central spin parameter,

$$\begin{array}{c}\hfill s_c=s_{\ast },\end{array}$$(14)

where ${\gamma }^{\prime }{}_{{\mathcal{l}}_i}^m$ is the derivative with respect to s of the function

$$\begin{array}{c}\hfill {\gamma }_{{\ell }_i}^m(s)={P\prime }_{{\ell }_i}^m(1/s)-\frac{m}{1-1/s^2}{P}_{{\ell }_i}^m(1/s),\end{array}$$(15)

where $P^{\prime }{}_{{\mathcal{l}}_i}^m$ is the derivative with respect to x of the Legendre polynomial $P_{{\mathcal{l}}_i}^m(x)$.

A first success of this model is that the Lorentzian form of the dips is fully consistent with our numerical results (see Sect. 3.1). The expressions of σ and s_c indicate that the dip parameters depend on the ratio Ω/N₀, on the gravito-inertial oscillations through ${\mathrm{\Lambda }}_k^m$, and on the inertial oscillation through s_* and $R(s_{*})=2P_{{\mathcal{l}}_i}^m(1/s_{*})/{\gamma }^{\prime }{}_{{\mathcal{l}}_i}^m(s_{*})$. Tokuno & Takata (2022) argued that these expressions remain relevant in the case of a realistic core with a non-uniform density, although the quantities related to the inertial oscillation, that is s_* and the function R(s), are no longer known analytically. We can nevertheless determine s_* numerically by computing the inertial modes of stellar models truncated at the boundary of the convective core where the radial displacements vanish (see Appendix B). This allowed us to test the s_c = s_* prediction against the numerical determination of s_c for the dips studied in this paper. As shown in Fig. 10, our numerical results are not consistent with s_c = s_*.

The theoretical expression of σ cannot be fully tested with our numerical results because the function R(s), which is derived from the analytical form of free inertial oscillations at the convective-radiative interface, is unknown for cores of variable density. Nevertheless, the proportionality σ ∝ Ω can be tested by considering the evolution of dip widths with rotation for a given star model and (ℓ,m) series. Figure 12 shows that σ ∝ Ω is indeed in agreement with our numerical results despite the small deviations observed in the case of the most evolved star model. Tokuno & Takata (2022) also treated the case of a continuous Brunt-Väisälä frequency profile at the core interface and found a different scaling of the dip width with the rotation. We note that even the small jump of the Brunt-Väisälä in the 1z model produces an evolution of σ corresponding to the discontinuous case. A questionable property of the theoretical σ is that each inertial mode with the same m and equatorial parity as the gravito-inertial modes of a (ℓ,m) series produces a dip of non-negligible width. As already discussed in Tokuno & Takata (2022), this property is at odds with the Ouazzani et al. (2020) phenomenological model, where an efficient resonance requires a non-negligible geometrical matching between the gravito-inertial and inertial modes at the convective-radiative interface. This matching selects the significantly resonant inertial modes and thus the dips in the ΔP = f(s), a prediction that was successfully tested with full numerical calculations by Ouazzani et al. (2020).

5.2. The improved model

In this subsection, we present a newly improved analytical model of the dips studied in this paper. This model is similar to the Tokuno & Takata (2022) model in that it uses the same three basic assumptions: uniform density core, TAR gravito-inertial oscillations, and Ω/N₀ ≪ 1. However, it goes a step further by improving on the geometrical matching of the inertial and gravito-inertial oscillations at the convective-radiative interface. This matching is indeed required by the continuity condition on the radial displacement and the pressure perturbations at this interface. Figure 17 shows the latitudinal profiles of the gravito-inertial mode and the inertial mode involved in the (ℓ = 1, m = −1) versus (ℓ_i = 3, m = −1) resonance. It is clear that the continuity condition cannot be realised by only considering these two modes, as is done in Tokuno & Takata (2022). We improved the matching by allowing more complex eigenfunctions on both sides of the interface. The detail of the derivation is given in Appendix D.

Fig. 17. Latitudinal profiles of the ℓi = 3, m = −1, s* = 11.32 inertial mode ($P_3^{-1}$(cos θ), blue) and ℓ = 1, m = −1 gravito-inertial modes near s* ($H_0^{-1}(\theta ,s_{*})$, orange) at the convective-radiative interface. Mixed modes that are continuous at the interface can hardly be constructed with these two modes only.

We find that the coupling between the (ℓ = 1, m = −1),(ℓ = 2, m = −2),(ℓ = 1, m = 0) gravito-inertial modes and respectively the (ℓ_i = 3, m = −1),(ℓ_i = 4, m = −2),(ℓ_i = 3, m = 0) inertial modes produces Lorentzian-shaped dips in the ΔP = f(s) relation whose widths and central spin parameters are described by

$$\begin{array}{cc}\hfill \sigma =& \left(\frac{f_{k,{\ell }_i}^m(s_{\ast })}{1+f_{k,{\ell }_i}^m(s_{\ast })}\frac{2P_{{\ell }_i}^m(1/s_{\ast })}{{{\gamma }^{\prime }}_{{\ell }_i}^m(s_{\ast })}\sqrt{{\mathrm{\Lambda }}_k^m(s_{\ast })}\right)\frac{\mathrm{\Omega }}{N_0},\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill s_c=& s_{\ast }+\left(\frac{1}{1+f_{k,{\ell }_i}^m(s_{\ast })}\frac{2P_{{\ell }_i}^m(1/s_{\ast })}{{\gamma \prime }_{{\ell }_i}^m(s_{\ast })}\sqrt{{\mathrm{\Lambda }}_{k-2}^m(s_{\ast })}\right)\frac{\mathrm{\Omega }}{N_0},\hfill \end{array}$$(17)

where $f_{k,{\mathcal{l}}_i}^m$ is defined by

$$\begin{array}{c}\hfill f_{k,{\ell }_i}^m(s_{\ast })=-\frac{c_{k,{\ell }_i}^m(s_{\ast })c_{k-2,{\ell }_i-2}^m(s_{\ast })}{c_{k,{\ell }_i-2}^m(s_{\ast })c_{k-2,{\ell }_i}^m(s_{\ast })}.\end{array}$$(18)

Here, $c_{k,{\mathcal{l}}_i}^m$ characterises the geometrical coupling between inertial and gravito-inertial oscillations along the convective-radiative interface:

$$\begin{array}{c}\hfill c_{k,{\ell }_i}^m(s)={\displaystyle {\int }_0^{\pi }}{\stackrel{\sim }{P}}_{{\ell }_i}^m(cos\theta )H_k^m(\theta ,s)sin\theta \mathrm{d}\theta ,\end{array}$$(19)

with ${\stackrel{\sim }{P}}_{{\ell }_i}^m(cos\theta )$ and $H_k^m(\theta ,s)$ being the normalised Legendre polynomial and Hough functions. In addition to $P_{{\mathcal{l}}_i}^m$ and $H_k^m$, these mixed modes involve $P_{{\mathcal{l}}_i-2}^m$ and $H_{k-2}^m$. In Appendix D, we argue this is a good approximation for the first dip of the (ℓ = 1, m = −1), (ℓ = 2, m = −2), and (ℓ = 1, m = 0) series studied in this paper. This model cannot be generalised to all dips, as we expect some resonances to produce more complex mixed modes.

The dip properties predicted by this model have important differences with Tokuno & Takata (2022), and these differences appear to improve the model. First, the width σ is multiplied by the factor $f_{k,{\mathcal{l}}_i}^m/(1+f_{k,{\mathcal{l}}_i}^m)$. As $f_{k,{\mathcal{l}}_i}^m$ is proportional to the coupling coefficient $c_{k,{\mathcal{l}}_i}(s_{*})$ introduced in Ouazzani et al. (2020), this new expression of σ can account for the fact that significant dips only occur when the resonance involves modes with a non-negligible geometrical matching at the convective-radiative interface. Second, the new form of the central spin parameter is now consistent with the numerical results as s_c − s_* ∝ Ω.

When extrapolated to a realistic core with a non-uniform density, our model suggests that the first dip of the (ℓ = 1, m = −1), (ℓ = 2, m = −2), and (ℓ = 1, m = 0) series can be expressed as

$$\begin{array}{cc}\hfill \sigma =& \beta \frac{\mathrm{\Omega }}{N_0},\hfill \end{array}$$(20)

$$\begin{array}{cc}\hfill s_c=& s_{\ast }+\alpha \frac{\mathrm{\Omega }}{N_0},\hfill \end{array}$$(21)

where β, s_*, and α depend on both the gravito-inertial modes and on the resonant inertial mode. For a fixed (ℓ,m), these parameters only depend on the properties of the inertial mode, that is, on the radial gradient of the density profile in the convective core (see Eq. (1) in Ouazzani et al. 2020). In the absence of analytical expressions, they can be derived by fitting the dips computed for different star models. In Sect. 4, using three star models, it was found that β remains roughly constant, while α is approximately proportional to 1/s_*. An open question for future and more detailed empirical relations is whether β and α can indeed be described by functions of s_* only. In that is the case, investigating the link between the seismic observable s_* and the star properties (age, mass, overshooting, etc.) would be much facilitated. Computing s_* for a given stellar model is indeed much easier than computing the frequency pattern and the dip parameters of a (ℓ,m) mode series.

Because of the uniform density assumption, a quantitative agreement between the theoretical and numerical values of σ and s_c is not expected. We nevertheless compared them by calculating the theoretical values of β_th(s_*), α_th(s_*), and A_th = α_th(s_*)s_*. Table 6 presents these parameters for the three dips we studied. It shows that despite the uniform density assumption, the theoretical values are in the range of those obtained by numerical calculations (see Tables 2, 3, 4, and 5). This is consistent with the fact that the spin parameter and the surface distribution of the inertial modes are not strongly modified by the core density stratification (see Appendix B). We also note that as expected for α, the deviations from the uniform density values increase for more evolved star models with higher density contrasts.

6. Discussion and conclusion

We have used the code TOP to study the resonances between inertial modes in the convective core and gravito-inertial modes in the radiative zone in γ Dor stars. As in Ouazzani et al. (2020), we obtained dips in the ΔP_co − s relation around spin parameters of gravito-inertial modes approximately predicted by a simple analytical model. We computed those dips for different series of modes and different rotation rates and evolutionary stages across the main sequence. We quantified the dips using a Lorentzian function as suggested in Tokuno & Takata (2022), which allowed us to determine the evolution of their basic properties with the rotation and the evolutionary stage of the star. We then proposed an empirical model for the evolution of the dips.

We also compared our numerical results with the model of Tokuno & Takata (2022), showing that it correctly predicts that the dip width is proportional to the rotation rate, but it does not account for the variation of the dip central spin parameter with rotation. We have presented a new model that accounts for this variation and selects the significant dips by including the geometrical matching between the inertial mode and the gravito-inertial modes. The model is not fully predictive, but it is useful to support and discuss the proposed empirical relations.

Our empirical model gives access to the Brunt-Väisälä frequency jump at the core interface as well as the resonant inertial mode frequency for stars presenting a clear dip in their period spacing patterns. Such dips have been observed, as shown in Saio et al. (2021), and further research could be led in order to find dips in a larger number of γ Dor stars in the TESS and Kepler data. This could be challenging, as modes in the dips may be harder to observe due to their higher inertia, as their energy is more concentrated in their core. As our empirical model relies only on two different stellar models, it needs to be tested on more stellar models in order to estimate its accuracy. Such a study is not trivial, as computing and identifying series of gravito-inertial modes with 2D oscillation codes can be a tedious process. The resonance phenomenon between the core and the envelope may not be exclusive to γ Dor stars, and similar dip structures could arise from the ΔP – P diagram of other types of stars. Slow pulsating B (SPB) stars are good candidates even though no dips have been observed yet (Aerts & Mathis 2023). A thorough investigation taking into account the expected spin parameters of the resonant inertial modes would be needed to establish whether resonances occur in the frequency ranges observed in SPB stars.

The main limitation of our present numerical study is related to the modelling of the rotation: We neglected the distortion induced by the centrifugal force and assumed a uniform rotation profile. For rotation rates observed in γ Dor stars, we expect a weak impact of the centrifugal distortion on our results since it mainly affects the external layers of the stars. Indeed, the structure of convective cores and the surrounding regions where g modes have large amplitudes remains quasi-spherical (see, e.g., Ballot et al. 2010, 2012). However, the influence of the differential rotation should be carefully investigated in a further study, especially when the core spins at a rotation that is different from the radiative zone. That kind of differential rotation has already been investigated and identified through the study of dips observed in 16 γ Dor stars in Saio et al. (2021). Although their differential rotation is shown to be small to non-significant in most cases, the core of one of the 16 stars was found to rotate 20% faster than the surrounding radiative zone. Such differential rotation would cause the frequency of the resonant inertial mode to be misestimated, as its estimate relies on the use of the inner radiative zone rotation rate obtained through the TAR. Nevertheless, how such differential rotation would affect the width of the dips and the estimate of the Brunt-Väisälä jump at the core interface remains unclear. If β remains constant, the width is only related to the radiative zone rotation rate rather than the core rotation rate. It implies that the estimate of the Brunt-Väisälä frequency jump through the dip width and the rotation rate of the inner radiative zone should remain correct despite the presence of such differential rotation.

We neglected non-adiabatic effects in our analysis. This approximation is supported by non-adiabatic calculations of high-order g modes, which have been performed to investigate the excitation mechanism of γ Dor pulsations (Guzik et al. 2000; Dupret et al. 2005; Bouabid et al. 2013). These calculations showed that the adiabatic analysis provides precise enough oscillation frequencies (Dupret et al. 2005). The spatial distribution of the g modes is affected by non-adiabatic effects, but those are only significant in the upper layers of the star, where the thermal relaxation time is smaller or of the same order as the oscillation period (Dupret et al. 2002, 2005). We thus expect that deep in the star, at the interface of the convective core and the radiative zone, the mode geometry is not affected by non-adiabatic effects. Thus, the adiabatic approximation we used to model the coupling between modes in the convective core and in the radiative envelope should be precise enough.

Our stellar structure models include turbulent diffusion to reproduce the different mixing processes that occur in stars. The value of the diffusion coefficient we used is the one proposed by Ouazzani et al. (2020), which is based on a previous calibration to Genova models (Miglio et al. 2008). The mixing is thus very efficient and prevents the formation of sharp features in the Brunt-Väisälä frequency profile in the g-mode cavity. Such features are known to generate glitches, visible in ΔP − s diagrams as oscillations or periodic dips (see for example Miglio et al. 2008). While absent in our computations, glitches may be present in real data and perturb the characterisation or even the detection of dips. Conversely, dips caused by the resonance phenomenon can potentially be mistaken for a glitch signature and could lead to wrong estimations of stellar parameters (see Mombarg et al. 2021).

The frequency and the width of the dips depends on the structure of the convective core through the spin parameter and the spatial distribution of the resonant core inertial modes. As shown by Eq. (B.1), it is the density stratification of the convective core that controls the inertial mode characteristics. Thus, the evolutionary stage or the stellar mass affects the dips as far as they affect the convective core density stratification. It is worth noting that the core size itself has no direct impact on the inertial modes. Further investigations are needed to understand the relation between the frequencies of the inertial modes in the core and the stellar properties. Such studies would thus allow one to use the said frequencies extracted by our empirical model to constrain stellar parameters. In particular, coupled with the jump in the Brunt-Väisälä frequency at the core interface we extracted from our model and the buoyancy radius deduced from the TAR, the inertial mode frequencies could be a powerful tool to determine stellar ages. Determining precise and accurate ages will be crucial for the PLATO mission (Rauer et al. 2014), which will use γ Dor stars as scientific calibrators.

Acknowledgments

We thank the anonymous referee for carefully reading our manuscript and providing helpful and constructive comments. We thank R.M. Ouazzani, M. Takata and D.R. Reese for very useful discussions and for performing some tests and comparisons, which were crucial to validate our computations with TOP. We acknowledge support from the Centre National d’Etudes Spatiales (CNES). This work was supported by the “Programme National de Physique Stellaire” (PNPS) of CNRS/INSU co-funded by CEA and CNES.

References

Appendix A: Uniform density model

Here we briefly present the uniform-density resonance model proposed in Ouazzani et al. (2020) and apply it to the gravito-inertial mode series studied in this paper. When there is no density stratification, the wave equation reduces ${\mathrm{\nabla }}^2\mathrm{\Psi }-s^2\frac{{\partial }^2\mathrm{\Psi }}{\partial s^2}$ = 0, the Poincaré equation. It is separable using the ellipsoidal coordinates (x₁, x₂, ϕ) (see Eqs. D.1, D.2, D.3).

Using regularity conditions and the boundary condition ξ_r = 0 applied at the sphere surface, it becomes possible to calculate the eigenfrequencies of the inertial modes as they are the non-trivial positive roots of

$$\begin{array}{c}\hfill \frac{d{P}_{{\ell }_i}^m(\mu )}{d\mu }=\frac{m}{1-{\mu }^2}{P}_{{\ell }_i}^m(\mu ),\end{array}$$(A.1)

where μ = 1/s = ω_co/2Ω and ℓ_i is the degree of the inertial mode. The spin parameters of inertial modes that could couple with gravito-inertial series studied here are presented in Table A.1. For an inertial mode in the core to be significantly coupled with a gravito-inertial mode in the radiative zone, there needs to be spatial and temporal correspondences. The two coupling modes must share a similar latitudinal profile at the interface between the core and the radiative zone and must oscillate with a similar frequency. The condition on the frequency is easily fulfilled as the frequency spectrum of gravito-inertial modes is dense. To quantify the spatial correspondence between an inertial and a gravito-inertial mode, Ouazzani et al. (2020) uses a correlation coefficient between the associated Legendre polynomial (which describes the inertial mode in the core in the analytical model) and the Hough function (which describes the gravito-inertial mode in the envelope in the TAR approximation). Those coefficient are presented in parentheses next to the spin parameter of the modes in Table A.1. The modes (ℓ_i = 3 m = −1), (ℓ_i = 3 m = 0), and (ℓ_i = 4 m = −2) with the highest correlation coefficients are shown in Table B.2.

Appendix B: Core-only model

We run TOP with a truncated version of our models to obtain the resonant frequencies. Table B.1 shows the spin parameters of the inertial modes coupling with the series of gravito-inertial modes (ℓ = 1 m = −1), (ℓ = 1 m = 0), and (ℓ = 2 m = −2) for the 3 models. They are labelled as (ℓ_i = 3 m = −1), (ℓ_i = 3 m = 0), and (ℓ_i = 4 m = −2) respectively as they remain similar to the modes calculated with the uniform density model presented in Appendix A (see Table B.2). Despite their apparent resemblance, it has to be noted that their spin parameters vary significantly with the density stratification as the relative differences between the 2m and 3t models’ spin parameters and the uniform-density model exceed 20%. To stress their differences, we show the mode profiles along the equator for both the uniform-density and 3t models in Fig. B.1.

Fig. B.1. Eulerian pressure perturbations normalised by their value at the core radius Rc of the (ℓi = 3, m = −1) mode, as a function of the radius at the colatitude θ = π/2 (equator), for the uniform density model (blue) and for the 3t model (red).

The rotation rate has little to no impact on the spin parameter and spatial distribution of inertial modes, as expected from the equation governing adiabatic waves in an isentropic convective core under the Cowling approximation (Wu 2005):

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2\mathrm{\Psi }-s^2\frac{{\partial }^2\mathrm{\Psi }}{\partial z^2}=-\frac{1}{{\rho }_0}\frac{\mathrm{d}{\rho }_0}{\mathrm{d}r}(\frac{\partial \mathrm{\Psi }}{\partial r}-s^{2}cos\theta \frac{\partial \mathrm{\Psi }}{\partial z}+\frac{\mathit{ms}}{r}\mathrm{\Psi })-(1-s^2)\frac{{\omega }_{\mathrm{co}}^2}{c_0^2}\mathrm{\Psi }\end{array}$$(B.1)

where $\mathrm{\Psi }=p^{\prime }/({\rho }_0{\omega }_{\text{co}}^2)$) and z = r cos θ is the coordinate alongside the rotation axis. Indeed, the acoustic term, the last on the right-hand side, which could introduce an explicit dependence on the rotation rate, is negligible with respect to the first term on the left-hand side in the convective core of γ Dor stars.

Appendix C: Results and analysis for the (ℓ = 1, m = 0) and (ℓ = 2, m = −2) series

In this appendix, we present our results for (ℓ = 1, m = 0) and (ℓ = 2, m = −2) series. We computed and analysed these series of modes in the same way as the (ℓ = 1, m = −1) series, which is extensively described in the main text.

C.1. (ℓ = 1, m = 0) series

The dips in the series of modes (ℓ = 1, m = 0) computed for the models 1z, 2m, and 3t at different rotations are presented in Figs. C.1, C.2, and C.3.

Fig. C.1. Period spacing in the co-rotating frame as a function of the spin parameter for the modes of the series (ℓ = 1, m = 0) at different rotation rates using the 1z model. The period spacing is normalised by the buoyancy radius Π0. The pink line shows the traditional approximation of rotation. Modes that are observable are marked with a square-shaped dot (Li et al. 2020).

Fig. C.2. Same as Fig. C.1 for the model 2m.

Fig. C.3. Same as Fig. C.1 for the model 3t.

The evolution of the width σ and location s_c of dips as a function of the rotation are shown in Figs. C.4 and C.5. In the latter, we also show the evolution of the spin parameter s_* of the resonant inertial mode of core-only models and the spin parameter of the same mode in a uniform density sphere.

Fig. C.4. Width of the dips σ as a function of the rotation rate for the series (ℓ = 1, m = 0), for the three stellar models described in Table 1.

Fig. C.5. Evolution of the centre of the dip sc obtained using a full model (full line) and evolution of the spin parameter of the inertial mode obtained using a core-only model (dotted line). The point-and-dot blue line shows the spin parameter of the inertial mode in the case of a uniform-density core (analytical model). Light blue curves refer to the model 1z, orange curves to the model 2m, and red curves to the model 3t.

We plot σ and s_c/s_* as a function of ϵ = Ω/N₀ for the three models in Figs. C.6 and C.7, along with linear fits of the form σ = βϵ and s_c/s_* = Aϵ + 1. We can see that σ is a lot smaller for the (ℓ = 1, m = 0) series than it is for the (ℓ = 1, m = −1), which is consistent with our analytical model (see the β_th column of Table 6). The ratio s_c/s_* is also smaller, which is in agreement with our analytical model (A_th column of Table 6).

Fig. C.6. Width σ as a function of ϵ for the series (ℓ = 1 m = 0) for 3 different evolutionary stages (1z, 2m, 3t). The purple dashed line shows the fit presented in Table 2 for the three models together and the purple area shows the related error. The x-axis is logarithmic for visualisation purposes.

Fig. C.7. Ratio sc/s* as a function of Ω/N0 (ϵ) for the series (ℓ = 1 m = 0), for the 1z, 2m, and 3t models. The purple dashed line shows the fit presented in Table 4 for the three models together and the purple area shows the related error. The x-axis is logarithmic for visualisation purposes.

In Fig. C.8, we compare the relative evolution of σ with rotation in our numerical computations and in the theoretical model of Tokuno & Takata (2022). Compared to the (ℓ = 1, m = −1) series, the non-linearity of σ for the 3t model is less clear. As the dip were less well described by a Lorentzian function in this case, the fit errors are significant.

Fig. C.8. Evolution of σ/σref with the rotation rate, for the dip of the series (ℓ = 1, m = 0) for three different evolutionary stages (1z, 2m, 3t), with the associated error bars. The reference value σref is the width for Ω = Ωref = 0.2ΩK. If, as predicted by Tokuno & Takata (2022), σ ∝ Ω, σ/σref = Ω/Ωref = 5Ω (the dashed line).

We used the empirical relations (10) and (12) on our numerical calculations to estimate s_* and N₀ and compared them to their true values. The errors made are shown in Figs. C.9 and C.10. In the first one, we used the three models and in the latter only the models 1z and 2m. Relative errors behave quite differently when compared to the ℓ = |m| (Kelvin) modes. The estimate for N₀ is bad for the 3t model regardless of the rotation rate, making the estimate of s_* also worse for the 3t model.

Fig. C.9. Relative errors for the series (ℓ = 1 m = 0) on the estimated $N_0^{\text{calc}}$ (upper panel) and $s_{*}^{\text{calc}}$ (lower panel) using the three models 1z, 2m, and 3t.

Fig. C.10. Relative errors for the series (ℓ = 1 m = 0) on the estimated $N_0^{\text{calc}}$ (upper panel) and $s_{*}^{\text{calc}}$ (lower panel) using the two models 1z and 2m.

C.2. (ℓ = 2, m = −2) series

The dips in the series of modes (ℓ = 2, m = −2) computed for the models 1z, 2m, and 3t at different rotations are presented in Figs. C.11, C.12 and C.13.

Fig. C.11. Same as Fig. C.1 for the series (ℓ = 2, m = −2).

Fig. C.12. Same as Fig. C.2 for the series (ℓ = 2, m = −2).

Fig. C.13. Same as Fig. C.3 for the series (ℓ = 2, m = −2).

The evolution of the width σ and location s_c of dips as a function of the rotation are shown in Figs. C.14 and C.15. In the latter, we also show the evolution of the spin parameter s_* of the resonant inertial mode of core-only models and the spin parameter of the same mode in a uniform density sphere.

Fig. C.14. Same as Fig. C.4 for the series (ℓ = 2, m = −2).

Fig. C.15. Same as Fig. C.5 for the series (ℓ = 2, m = −2).

We plot σ and s_c/s_* as a function of ϵ = Ω/N₀ for the three models in Figs. C.16 and C.17, along with linear fits of the form σ = βϵ and s_c/s_* = Aϵ + 1. At a given ϵ, both σ and s_c/s_* are slightly smaller for the (ℓ = 2, m = −2) series compared to the (ℓ = 1, m = −1) series. This behaviour is consistent with our analytical model (see the β_th and A_th column of Table 6).

Fig. C.16. Same as Fig. C.6 for the series (ℓ = 2, m = −2).

Fig. C.17. Same as Fig. C.7 for the series (ℓ = 2, m = −2).

In Fig. C.18, we compare the relative evolution of σ with rotation in our numerical computations and in the theoretical model of Tokuno & Takata (2022). As for the series (ℓ = 1, m = −1, while σ is clearly linear for the 1z and 2m models, it is not the case for the 3t model. One can see that this tendency is more pronounced for the (ℓ = 2, m = −2) series than it is for the (ℓ = 1, m = −1) series.

Fig. C.18. Same as Fig. C.8 for the series (ℓ = 2, m = −2).

We used the empirical relations (10) and (12) on our numerical calculations to estimate s_* and N₀ and compared them to their true values. The errors made are shown in Figs. C.19 and C.20. In the first one, we used the three models and in the latter only the models 1z and 2m. One can see that the relative errors on N₀ and s_* follow the same kind of trend as for the (ℓ = 1, m = −1) series.

Fig. C.19. Same as Fig. C.9 for the series (ℓ = 2, m = −2).

Fig. C.20. Same as Fig. C.10 for the series (ℓ = 2, m = −2).

Appendix D: Analytical model of the resonant coupling between inertial and gravito-inertial oscillations.

We present an analytical model of mixed (inertial and gravito-inertial modes) when the density is assumed to be uniform in the convective core. This model is similar to the one proposed by Tokuno & Takata (2022) in that its uses the same hypothesis to derive analytical solutions in the convective core and in the radiative zone respectively. However, it goes a step further by improving the continuity of the mixed modes at the convective-radiative interface.

In the following, we determine analytical solutions of free oscillations in the convective core (Sect. D.1) and then in the radiative envelope (Sect. D.2). In Sect. D.3, the continuity conditions at the interface are written and approximated in the case of (ℓ,m) = (1, −1) gravito-inertial modes interacting with the (ℓ_i, m) = (3, −1),s_* = 11.32 inertial mode. The effect of this mode coupling on the ΔP = f(s) relation is then derived.

D.1. Analytical solutions in the convective core

In this subsection, we describe small amplitude inertial oscillations in the convective core with no specified boundary conditions at the interface with the radiative zone. In an axially symmetric star, we search for solutions proportional to e^{i(mϕ + ωt)} where m is the azimuthal number. The radial displacement and the Eulerian pressure perturbations are thus expressed as ${\xi }_r^{\text{conv}}(r,\theta )e^{i(m\varphi +\omega t)}$ and p′^conv(r, θ)e^{i(mϕ + ωt)} respectively. In an isentropic convective core of uniform density, perturbations of low frequency are governed by the Poincaré equation $\mathrm{\Delta }\mathrm{\Psi }=s^2\frac{\partial \mathrm{\Psi }}{\partial z}$ where Ψ is proportional to p′^conv (Wu 2005). In a full sphere, separable solutions are obtained using the ellipsoidal coordinates (x₁, x₂, ϕ) defined in terms of the Cartesian coordinates as

$$\begin{array}{cc}\hfill x/r_c& ={\left(\frac{(1-x_1^2)(1-x_2^2)}{1-{\mu }_s^2}\right)}^{1/2}cos(\varphi )\hfill \end{array}$$(D.1)

$$\begin{array}{cc}\hfill y/r_c& ={\left(\frac{(1-x_1^2)(1-x_2^2)}{1-{\mu }_s^2}\right)}^{1/2}sin(\varphi )\hfill \end{array}$$(D.2)

$$\begin{array}{cc}\hfill z/r_c& =\frac{x_1{x}_2}{{\mu }_s}\hfill \end{array}$$(D.3)

where μ_s = 1/s, x₁ ∈ [μ_s, 1], and x₂ ∈ [ − μ_s, μ_s]. With these coordinates, the surface of the sphere is given by {x₁ ∈ [μ_s, 1],x₂ = μ_s}∪{x₁ = μ_s, x₂ ∈ [ − μ_s, μ_s]}.

The separable solutions read:

$$\begin{array}{c}\hfill \mathrm{\Psi }=r_c^2{P}_{{\ell }_i}^m(x_1)P_{{\ell }_i}^m(x_2)\end{array}$$(D.4)

where $P_{{\mathcal{l}}_i}^m$ is the associated Legendre polynomial of degree ℓ_i. The radial displacement ${\xi }_r^{\text{conv}}$ and the Eulerian pressure perturbations p′^conv are then obtained from the relations (Wu 2005):

$$\begin{array}{cc}\hfill {\xi }_r^{\mathrm{conv}}=& \frac{1}{1-s^2}(\frac{\partial \mathrm{\Psi }}{\partial r}+\frac{\mathit{ms}}{r}\mathrm{\Psi }-s^{2}cos\theta \frac{\partial \mathrm{\Psi }}{\partial z})\hfill \end{array}$$(D.5)

$$\begin{array}{cc}\hfill p^{\prime \mathrm{conv}}=& {\rho }_0{\omega }^2\mathrm{\Psi }.\hfill \end{array}$$(D.6)

At the outer radius of the convective core r = r_c, it reads (Tokuno & Takata 2022):

$$\begin{array}{cc}\hfill {\xi }_r^{\mathrm{conv}}(r=r_c)=& r_c\frac{1}{s}{\gamma }_{{\ell }_i}^m(s)P_{{\ell }_i}^m(cos\theta )\mathrm{with}\hfill \end{array}$$(D.7)

$$\begin{array}{cc}\hfill {\gamma }_{{\ell }_i}^m(s)=& {P^{\prime }}_{{\ell }_i}^m(1/s)-\frac{m}{1-1/s^2}{P}_{{\ell }_i}^m(1/s)\hfill \end{array}$$(D.8)

where $P^{\prime }{}_{{\mathcal{l}}_i}^m$ is the derivative of the polynomial $P_{{\mathcal{l}}_i}^m$, and

$$\begin{array}{cc}\hfill p^{\prime \mathrm{conv}}(r=r_c)=& {\rho }_0{r}_c^2{\omega }^2{P}_{{\ell }_i}^m(1/s)P_{{\ell }_i}^m(cos\theta ).\hfill \end{array}$$(D.9)

These solutions oscillate freely at the convective-radiative interface because no boundary condition is imposed there. If the ${\xi }_r^{\text{conv}}(r=r_c)=0$ boundary condition is enforced, one recovers the eigenfrequencies of pure inertial modes presented in Appendix A which satisfy Eq. A.1 or equivalently ${\gamma }_{{\mathcal{l}}_i}^m(s)=0$.

D.2. Analytical solutions in the radiative zone

In the radiative zone, analytical gravito-inertial oscillations with no prescribed boundary conditions at the convective interface can also be found in the framework of the TAR and in the limit of short-radial wavelengths. These solutions are separable in the spherical coordinates and can thus be written ${\xi }_r^{\text{rad}}(r,\theta )=\zeta (r)H(\theta )$ and p′^rad(r, θ) = p′(r)H(θ), where ζ and p′ denote the radial part of the radial displacement and Eulerian pressure perturbation respectively, and H(θ) their latitudinal part.

The equation governing the latitudinal part of the solution is the Laplace’s tidal equation:

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{T}}(H)+\mathrm{\Lambda }H=0\end{array}$$(D.10)

with

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\mathrm{T}}(H)=& \frac{d}{d\mu }\left[\frac{1-{\mu }^2}{(1-s^2{\mu }^2)}\frac{d}{d\mu }H\right]\hfill \\ \hfill & -\frac{1}{1-s^2{\mu }^2}(\frac{m^2}{1-{\mu }^2}+ms\frac{1+s^2{\mu }^2}{1-s^2{\mu }^2})H,\hfill \end{array}$$

where μ = cos θ and Λ is a constant. Applying regularity conditions at the polar axis amounts to solve an eigenvalue problem. For each doublet (m, s), the eigenfunctions and their associated eigenvalues can be labelled with an integer k (Lee & Saio 1997). The eigenvalues and the eigenfunctions, known as the Hough functions, are thus denoted ${\mathrm{\Lambda }}_k^m(s)$ and $H_k^m(\theta ,s)$ respectively. They form an orthogonal set of functions as ${\displaystyle {\int }_0^{\pi }}{H}_k^m(\theta ,s)H_{k\prime }^m(\theta ,s)sin\theta d\theta$ = 0 if k ≠ k′. In the following, we shall use normalised Hough functions such that ${\displaystyle {\int }_0^{\pi }}{\left(H_k^m(\theta ,s)\right)}^{2}sin\theta d\theta =1$. Physically, a positive or null k corresponds to a gravity mode modified by the Coriolis force. In this case, the index k relates to the degree ℓ of the spherical harmonics of the gravity mode in a non-rotating star through ℓ = |m|+k. Thus, for example, the prograde sectoral modes, also called Kelvin modes, can be equivalently labelled (ℓ=|m|,m < 0) or (k = 0, m < 0). Negative k in turn are associated to solutions that only exist in rotating stars either the (retrograde) Rossby modes or the Λ < 0 modes that do not propagate in the radial direction. For the prograde or axisymmetric modes that we shall consider below, the k < 0 modes correspond to Λ < 0 non propagative modes.

Within the Cowling approximation, the equations governing the radial part of the oscillation, ζ and p′, are:

$$\begin{array}{cc}\hfill \frac{d\zeta }{\mathit{dr}}=& -(\frac{2}{r}-\frac{1}{{\mathrm{\Gamma }}_1{H}_P})\zeta +\frac{1}{{\rho }_0{c}_0^2}(\frac{\mathrm{\Lambda }{c}_0^2}{r^2{\omega }^2}-1)p^{\prime }\hfill \end{array}$$(D.11)

$$\begin{array}{cc}\hfill \frac{d{p}^{\prime }}{\mathit{dr}}=& {\rho }_0({\omega }^2-N^2){\xi }_r-\frac{1}{{\mathrm{\Gamma }}_1{H}_P}{p}^{\prime }\hfill \end{array}$$(D.12)

where ${\mathrm{\Lambda }}_k^m(s)$ is the eigenvalue of the Laplace’s tidal equation.

In the limit of short radial wavelengths, the WKB solutions of these equations depends on the sign of $k_r^2=\frac{{\omega }^2-N^2}{c_0^2}(1-\frac{{\mathrm{\Lambda }}_k^m(s)c_0^2}{r^2{\omega }^2})$ :

$$\begin{array}{cc}\hfill \mathrm{if}{k}_r^2>0& {\zeta }_k^m(r,s)=-c\frac{|1-\frac{\mathrm{\Lambda }{c}_0^2}{r^2{\omega }^2}{|}^{1/2}}{r{c}_0{(\pi {\rho }_0{k}_r)}^{1/2}}sin({\displaystyle {\int }_r^{r_b}}{k}_{r}dr-\frac{\pi }{4})\hfill \end{array}$$(D.13)

$$\begin{array}{cc}\hfill & \frac{{p^{\prime }}_k^m(r,s)}{{\rho }_0}=c\frac{{(N^2-{\omega }^2)}^{1/2}}{r{(\pi {\rho }_0{k}_r)}^{1/2}}cos({\displaystyle {\int }_r^{r_b}}{k}_{r}dr-\frac{\pi }{4})\hfill \end{array}$$(D.14)

$$\begin{array}{cc}\hfill \mathrm{if}{k}_r^2<0& {\zeta }_k^m(r,s)=c^{\prime }\frac{{(1-\frac{\mathrm{\Lambda }{c}_0^2}{r^2{\omega }^2})}^{1/2}}{r{c}_0(\pi {\rho }_0|k_r{|)}^{1/2}}exp(-{\displaystyle {\int }_{r_c}^r}|k_r|dr)\hfill \end{array}$$(D.15)

$$\begin{array}{cc}\hfill & \frac{{p^{\prime }}_k^m(r,s)}{{\rho }_0}=c^{\prime }\frac{{(N^2-{\omega }^2)}^{1/2}}{r(\pi {\rho }_0|k_r{|)}^{1/2}}exp(-{\displaystyle {\int }_{r_c}^r}|k_r|dr)\hfill \end{array}$$(D.16)

where c and c′ are arbitrary constants and the labels k, m, s of Λ and k_r have been omitted for clarity. In these expressions, r_b is the outer radius of the oscillation cavity which satisfies k_r(r_b) = 0 that is $\frac{{\mathrm{\Lambda }}_k^m(s)c_0{(r_b)}^2}{r_b^2{\omega }^2}\approx 1$ as ω ≪ N. The solution exponentially decreasing from r_c corresponds to Hough modes with negative Λ.

These oscillations are free at the convective/radiative interface. If a boundary condition ξ_r(r = r_c) = 0 were to be enforced, solutions ${\zeta }_k^m(r,s)H_k^m(\theta ,s)$ with ${\zeta }_k^m(r=r_c,s)=0$ would be the eigenmodes. According to Eq. D.13, this leads to the quantisation condition ${\displaystyle {\int }_{r_c}^{r_b}}{k}_r\hspace{0.17em}dr-\pi /4=n\pi$. As $\mathrm{\Lambda }{c}_0^2/(r^2{\omega }^2)\gg 1$ away from r_b and ω ≪ N, the condition can be simplified into $\frac{{\pi }^{2}s\sqrt{\mathrm{\Lambda }}}{\mathrm{\Omega }{\mathrm{\Pi }}_0}=\pi (n+1/4)$ with ${\mathrm{\Pi }}_0=\frac{2{\pi }^2}{{\displaystyle {\int }_{r_i}^{r_o}}\frac{N}{r}\mathrm{d}r}$.

D.3. Modes of mixed inertial and gravito-inertial nature

To produce mixed modes, the inertial and gravito-inertial oscillations have to match at the convective-radiative interface in the sense that ξ_r and δp have to be continuous there. However, as illustrated in Fig. 17 for the (ℓ = 1, m = −1) vs (ℓ_i = 3, m = −1) resonance, these matching conditions can not be fulfilled with only one separable solutions in each domain. We thus need more general solutions on both sides of the interface. This can be done by considering linear combinations of the separable solutions derived in the two previous sections. Namely, in the radiative zone :

$$\begin{array}{cc}\hfill {\xi }_r^{\mathrm{rad}}=& {\displaystyle \sum _{k=-\infty }^{k=+\infty }}{a}_k{\zeta }_k^m(r,s)H_k^m(\theta ,s)\hfill \end{array}$$(D.17)

$$\begin{array}{cc}\hfill p^{\prime \mathrm{rad}}=& {\displaystyle \sum _{k=-\infty }^{k=+\infty }}{a}_k{p^{\prime }}_k^m(r,s)H_k^m(\theta ,s)\hfill \end{array}$$(D.18)

and in the convective zone :

$$\begin{array}{cc}\hfill \mathrm{\Psi }\propto & r_c^2{\displaystyle \sum _{{\ell }_i=1}^{+\infty }}{b}_{{\ell }_i}{P}_{{\ell }_i}^m(x_1)P_{{\ell }_i}^m(x_1).\hfill \end{array}$$(D.19)

At the interface between the convective core and the radiative zone, the density is continuous while its radial derivative is discontinuous. In such a case, the continuity condition on ξ_r and δp is equivalent to the continuity of ξ_r and p′, thus to ${\xi }_r^{\text{conv}}(r=r_c^{-})={\xi }_r^{\text{rad}}(r=r_c^{+})$ and $p^{\prime }{}^{\text{conv}}(r=r_c^{-})=p^{\prime }{}^{\text{rad}}(r=r_c^{+})$.

The expressions of ${\xi }_r^{\text{rad}}(r_c^{+})$ and $p^{\prime }{}^{\text{rad}}(r_c^{+})$ are

$$\begin{array}{cc}\hfill {\xi }_r^{\mathrm{rad}}(r_c^{+})=& {\displaystyle \sum _{k=-\infty }^{k=+\infty }}{a}_k{\zeta }_k^m(r_c^{+},s)H_k^m(\theta ,s)\mathrm{with}\hfill \end{array}$$(D.20)

$$\begin{array}{cc}\hfill {\zeta }_k^m(r_c^{+},s)=& -\frac{1}{{\rho }_0^{1/2}\omega r_c^3}{\left(\frac{\mathrm{\Lambda }}{\pi k_r}\right)}^{1/2}sin(\psi )\mathrm{if}{k}_r^2>0\hfill \end{array}$$(D.21)

$$\begin{array}{cc}\hfill {\zeta }_k^m(r_c^{+},s)=& \frac{1}{{\rho }_0^{1/2}\omega r_c^3}{\left(\frac{|\mathrm{\Lambda }|}{\pi |k_r|}\right)}^{1/2}\mathrm{if}{k}_r^2<0\hfill \end{array}$$(D.22)

and

$$\begin{array}{cc}\hfill p_r^{\prime \mathrm{rad}}(r_c^{+})=& {\displaystyle \sum _{k=-\infty }^{k=+\infty }}{a}_k{p^{\prime }}_k^m(r_c^{+},s)H_k^m(\theta ,s)\mathrm{with}\hfill \end{array}$$(D.23)

$$\begin{array}{ccc}\hfill {p^{\prime }}_k^m(r_c^{+},s)=& \frac{{\rho }_0^{1/2}}{r_c}{\left(\frac{N_0^2-{\omega }^2}{\pi k_r}\right)}^{1/2}cos(\psi )\hfill & \hfill \mathrm{if}{k}_r^2>0\end{array}$$(D.24)

$$\begin{array}{ccc}\hfill {p^{\prime }}_k^m(r_c^{+},s)=& \frac{{\rho }_0^{1/2}}{r_c}{\left(\frac{N_0^2-{\omega }^2}{\pi |k_r|}\right)}^{1/2}\hfill & \hfill \mathrm{if}{k}_r^2<0\end{array}$$(D.25)

where $N_0=N(r=r_c^{+})$ and ${\psi }_k^m(s)={\displaystyle {\int }_{r_c}^{r_b}}{k}_{r}dr-\pi /4\approx \frac{{\pi }^{2}s\sqrt{\mathrm{\Lambda }}}{\mathrm{\Omega }{\mathrm{\Pi }}_0}-\pi /4$ because $\Lambda c_0^2/(r^2{\omega }^2)\gg 1$ away from r_b and ω ≪ N.

The expressions of ${\xi }_r^{\text{conv}}(r_c^{-})$ and $p^{\prime }{}^{\text{conv}}(r_c^{-})$ are

$$\begin{array}{cc}\hfill {\xi }_r^{\mathrm{conv}}(r_c^{-})=& \frac{r_c}{s}{\displaystyle \sum _{{\ell }_i=1}^{+\infty }}{b}_{{\ell }_i}{\gamma }_{{\ell }_i}^m(s){\stackrel{\sim }{P}}_{{\ell }_i}^m(cos\theta )\hfill \end{array}$$(D.26)

$$\begin{array}{cc}\hfill p^{\prime \mathrm{conv}}(r_c^{-})=& {\rho }_0{r}_c^2{\omega }^2{\displaystyle \sum _{{\ell }_i=1}^{+\infty }}{b}_{{\ell }_i}{P}_{{\ell }_i}^m(1/s){\stackrel{\sim }{P}}_{{\ell }_i}^m(cos\theta )\hfill \end{array}$$(D.27)

where, for convenience, we used ${\stackrel{\sim }{P}}_{{\ell }_i}^m$, the normalised form of the associated Legendre polynomial, i.e. ${\displaystyle {\int }_0^{\pi }}{\left({\stackrel{\sim }{P}}_{{\ell }_i}^m(cos\theta )\right)}^{2}sin\theta d\theta =1$.

The continuity conditions then read:

$$\begin{array}{cc}\hfill \frac{r_c}{s}{\displaystyle \sum _{{\ell }_i=1}^{+\infty }}{b}_{{\ell }_i}{\gamma }_{{\ell }_i}^m(s){\stackrel{\sim }{P}}_{{\ell }_i}^m(cos\theta )=& {\displaystyle \sum _{k=-\infty }^{k=+\infty }}{a}_k{\zeta }_k^m(r_c^{+},s)H_k^m(\theta ,s)\hfill \end{array}$$(D.28)