11institutetext: IRAP, Université de Toulouse, CNRS, CNES, UPS, 14 avenue
Edouard Belin, 31400 Toulouse, France
11email<EMAIL_ADDRESS>
# Properties of observable mixed inertial and gravito-inertial modes in
$\gamma$ Doradus stars

Marion Galoy  François Lignières  Jérôme Ballot

(Received 13 Mars 2024 / Accepted 3 July 2024)

###### Abstract

Context. The space missions Kepler and TESS provided a large number of highly
detailed time series for main-sequence stars, including $\gamma$ Doradus
stars. Additionally, numerous $\gamma$ Doradus stars are to be observed in the
near future thanks to the upcoming PLATO mission. In $\gamma$ Doradus stars,
gravito-inertial modes in the radiative zone and inertial modes in the
convective core can interact resonantly, which translates into the appearance
of dip structures in the period spacing of modes. Those dips are information-
rich, as they are related to the star core characteristics.

Aims. Our aim is to characterise these dips according to stellar properties
and thus to develop new seismic diagnostic tools to constrain the internal
structure of $\gamma$ Doradus stars, especially their cores.

Methods. We used the two-dimensional oscillation code TOP to compute sectoral
prograde and axisymmetric dipolar modes in $\gamma$ Doradus stars at different
rotation rates and evolutionary stages. We then characterised the dips we
obtained by their width and location on the period spacing diagram.

Results. We found that the width and the location of the dips depend quasi-
linearly on the ratio of the rotation rate and the Brunt-Väisälä frequency at
the core interface. This allowed us to determine empirical relations between
the width and location of dips as well as the resonant inertial mode frequency
in the core and the Brunt-Väisälä frequency at the interface between the
convective core and the radiative zone. We propose an approximate theoretical
model to support and discuss these empirical relations.

Conclusions. The empirical relations we established could be applied to dips
observed in data, which would allow for the estimation of frequencies of
resonant inertial modes in the core and of the Brunt-Väisälä jump at the
interface between the core and the radiative zone. As those two parameters are
both related to the evolutionary stage of the star, their determination could
lead to more accurate estimations of stellar ages.

###### Key Words.:

Asteroseismology – Stars: oscillations – Stars: rotation

##  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, $\gamma$ Doradus (hereafter $\gamma$ Dor) stars are
main-sequence stars that typically have a mass of 1.4–2.0 M⊙ and a radius of
1.4–2.3 R⊙ (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 3 d (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 $\gamma$ Dor stars.
Gravito-inertial modes have been identified in 611 $\gamma$ Dor stars observed
by _Kepler_ (Li et al., 2020) and in 106 $\gamma$ 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 $\gamma$ Dor stars to be determined (see Van Reeth et al., 2015,
2016; Ouazzani et al., 2017; Van Reeth et al., 2018; Saio et al., 2018;
Christophe et al., 2018; Mombarg et al., 2019; Ouazzani 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 $\gamma$ 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 $\Delta 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 $\gamma$ Dor stars. Saio et al. (2021)
analysed dips in 16 $\gamma$ 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 $\Delta 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 $\rm{}^{14}N(p,\gamma){}^{15}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
$\alpha_{\mathrm{MLT}}=1.70$, close to a solar calibration. We included
turbulent diffusion with a constant diffusion coefficient of
$D_{\mathrm{t}}=700\>\mathrm{cm^{2}\,s^{-1}}$. 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 $\gamma$
Dor stars (Bouabid et al., 2013). The model labelled ‘1z’ describes a young
star of $1.40\>\mathrm{M_{\odot}}$ close to the zero age main sequence (ZAMS;
age of $180$ My from the ZAMS, $X_{\rm c}=0.68$, with $X_{\rm c}$ being the
mass fraction of hydrogen in the core). The model labelled ‘2m’ describes an
evolved mid-main-sequence star of $1.60\>\mathrm{M_{\odot}}$ (age of $1800$ My
from the ZAMS, $X_{\rm c}=0.38$). The model labelled ‘3t’ describes a star of
$1.86\>\mathrm{M_{\odot}}$ at the end of the main sequence (age of $1480$ My
from the ZAMS, $X_{\rm c}=0.05$).

Table 1: Characteristics of the 1z, 2m, and 3t star models computed by CESAM.

Model name | 1z | 2m | 3t
---|---|---|---
$M/\rm M_{\odot}$ | 1.40 | 1.60 | 1.86
$T_{\rm eff}$ (K) | 6935 | 7289 | 6694
$\log L/\rm L_{\odot}$ | 0.692 | 1.052 | 1.379
$\log g$ | 4.31 | 3.99 | 3.59
$R/\rm R_{\odot}$ | 1.38 | 2.10 | 3.64
Age (Myr) | 180 | 1800 | 1480
$X_{\rm c}$ | 0.68 | 0.38 | 0.05
111Mass, effective temperature, luminosity, surface gravity, radius, age and
central hydrogen mass fraction _(from top to bottom)_

Figure 1: Squared Brunt-Väisälä frequency $N^{2}$ 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.

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_{\rm t}$ mentioned above.

###  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\phi+\omega t)}$, where $m$ is the azimuthal number and
$\omega$ is the mode frequency. The governing equations are

$\displaystyle\rho+\vec{\xi}\cdot\vec{\nabla}\rho_{0}+\rho_{0}\vec{\nabla}\cdot\vec{\xi}=0,$
(1) $\displaystyle
2i\rho_{0}\vec{\Omega}\times[\omega+m\Omega]\vec{\xi}-[\omega+m\Omega]^{2}\rho_{0}\vec{\xi}=-\vec{\nabla}p+\frac{\vec{\nabla}P_{0}}{\rho_{0}}\rho-\rho_{0}\vec{\nabla}\psi,$
(2) $\displaystyle
p-c_{0}^{2}\rho+\vec{\xi}\cdot\left(\vec{\nabla}P_{0}-c_{0}^{2}\vec{\nabla}\rho_{0}\right)=0,$
(3) $\displaystyle\vec{\Delta}\psi-4\pi G\rho=0,$ (4)

where $\vec{\xi}$, $p$, $\rho$, and $\psi$ are the amplitude of the Lagrangian
displacement, the Eulerian perturbations of pressure, density, and
gravitational potential, respectively, and where $\rho_{0}$, $P_{0}$, and
$c_{0}$ are the density, the pressure, and the sound velocity of the star
model, respectively, with $c_{0}=\sqrt{\Gamma_{1}p_{0}/\rho_{0}}$,
$\Gamma_{1}$ being the first adiabatic index. The star rotation rate is
$\Omega$, and $G$ is the gravitational constant. The mode frequency in the
inertial frame $\omega$ is related to the mode frequency in the co-rotating
frame $\omega_{\mathrm{co}}$ through

$\omega=\omega_{\mathrm{co}}-m\Omega.$ (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 $\ell$ and azimuthal number $m$ but different radial
orders $n$. The degree $\ell$ corresponds to the associated Legendre
polynomial $P_{\ell}^{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 ($\ell,m$) gravito-inertial mode can be
estimated by

$n+\epsilon_{g}\approx\frac{P_{\mathrm{co}}\sqrt{\Lambda_{\ell}^{m}}}{\Pi_{0}},$
(6)

with $P_{\mathrm{co}}$ as the period of the mode in the co-rotating frame,
$\Pi_{0}$ as the buoyancy radius, $\Lambda_{\ell}^{m}$ as the eigenvalue
associated with the Hough function $H_{\ell}^{m}$ describing the latitudinal
part of the mode in the TAR (see Appendix D.2), and $\epsilon_{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_{\theta}$ corresponds to the number of
associated Legendre polynomials $P_{\ell}^{m}$ on which the modes are
decomposed. The degree $\ell$ spans from $\ell_{\mathrm{min}}=|m|+i_{p}$ to
$\ell_{\mathrm{max}}=\ell_{\rm min}+2(n_{\theta}-1)$, where $i_{p}=0$ or 1
when $\ell+m$ is even or odd, respectively. For example, a mode of the series
($\ell=1$, $m=-1$) calculated with $n_{\theta}=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\Omega/\omega_{\mathrm{co}}\sim 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_{\theta}=6$ is needed for this mode.

Figure 2: Projection of the Hough function describing the latitudinal profile
of the mode $n=66$ ($\ell=1$, $m=-1$) in the radiative zone (TAR) on the
associated Legendre polynomials $P_{\ell_{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
($\ell=1,m=-1$) onto the associated Legendre polynomials for different
latitudinal resolutions. We first noticed that the amplitude decreases rapidly
with the degree $\ell_{i}$. Nevertheless, when we used too few polynomials, we
missed non-negligible contributions of polynomials of degree $\ell_{i}$ larger
than $\ell_{\mathrm{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 $\ell_{\mathrm{max}}$. For
this mode, the frequency varies by less than 0.01% for an
$\ell_{\mathrm{max}}$ higher than nine ($n_{\theta}=5$).

Figure 3: Effect of the latitudinal resolution on the mode $n=66$ ($s\sim
16.4$) of the series ($\ell=1$, $m=-1$). Upper panel: Projections of the mean
latitudinal profile in the radiative zone of this mode at different
latitudinal resolution $n_{\theta}$ on the associated Legendre polynomials for
a rotation rate $\Omega=0.4\Omega_{\mathrm{K}}$. Here,
$\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
$\ell_{\mathrm{max}}$. The reference frequency $\omega_{\mathrm{ref}}$
corresponds to the frequency obtained using the highest latitudinal resolution
($n_{\theta}=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 $\Delta P_{\mathrm{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 $\Omega=0.4\Omega_{\mathrm{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 $\Delta P_{\mathrm{co}}$ increases for the higher $s$ as the radial
resolution decreases, which increases the general slope of the curve. When the
latitudinal resolution decreases, $\Delta P_{\mathrm{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_{\theta}>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\times n$ points is enough to reach convergence.

Figure 4: Period spacing $\Delta P_{\mathrm{co}}$, normalised by $\Pi_{0}$
(buoyancy radius), plotted as a function of the spin parameter $s$ with
different radial resolution $n_{r}$ for the series of modes ($\ell=1$, $m=-1$)
and for a rotation rate $\Omega=0.4\Omega_{\mathrm{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. Figure 5: Period spacing $\Delta P_{\mathrm{co}}/\Pi_{0}$
plotted as a function of the spin parameter $s$ with different latitudinal
resolutions $n_{\theta}$ for the series of modes ($\ell=1$, $m=-1$) and for a
rotation rate $\Omega=0.4\Omega_{\mathrm{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
$\delta\Delta P_{\mathrm{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
$\delta\omega$ using

$\delta\omega=\omega_{\mathrm{co}}^{2}\frac{\delta\Delta
P_{\mathrm{co}}}{2\pi\sqrt{2}}.$ (7)

To ensure that $\delta\omega<1\%$, we showed that it is sufficient to choose a
radial resolution $n_{r}$ ranging from 5000 to 10000 for the ($\ell=1$,
$m=-1$) series, from 3000 to 10000 for the ($\ell=2$, $m=-2$) series, and from
3000 to 10000 for the ($\ell=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_{\theta}=6$ to $8$ for the series ($\ell=1$, $m=-1$) and
($\ell=2$, $m=-2$), and $n_{\theta}=4$ to $6$ for the series ($\ell=1$,
$m=0$).

Figure 6: Example of a resonant mode. Left panel: Period spacing $\Delta
P_{\mathrm{co}}$ as a function of the spin parameter for the series ($\ell=1$,
$m=-1$) and for the 1z model at a rotation rate of
$\Omega=0.1\Omega_{\mathrm{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 $\pi$/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.

##  3 Results

In this section, we present the results of our numerical calculations. We
studied the series of modes ($\ell=1$, $m=-1$), ($\ell=2$, $m=-2$), and
($\ell=1$, $m=0$) at rotation rates ranging from $0.1$ to $0.5$
$\Omega/\Omega_{\mathrm{K}}$, with $\Omega_{\mathrm{K}}=\sqrt{GM/R^{3}}$, for
the three stellar models described in Table 1. Among these three series, the
($\ell=1$, $m=-1$) is by far the most frequently observed (Li et al., 2020).
For each series, we computed the period spacing $\Delta P_{\mathrm{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 $\Delta P_{\mathrm{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
($\ell=1$, $m=-1$), ($\ell=2$, $m=-2$), and ($\ell=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 ($\ell=1$,
$m=-1$). The dip is located around $s=10.7$ and was calculated with the 1z
model at rotation rate $\Omega=0.1\Omega_{\mathrm{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 $\Delta
P_{\mathrm{co}}-s$ curve around $s=10.7$ in our calculations, the TAR is
nearly constant. For prograde sectoral modes ($\ell=-m$), the TAR indeed leads
to a constant $\Delta P_{\mathrm{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

$\frac{E_{\mathrm{kin}}(\mathrm{core})}{E_{\mathrm{kin}}(\mathrm{core+RZ})}=\frac{\int_{0}^{r_{c}}\int_{0}^{2\pi}\int_{0}^{\pi}\rho_{0}\left|\vec{v}\right|^{2}r^{2}\sin{\theta}\,\mathrm{d}\theta\,\mathrm{d}\phi\,\mathrm{d}r}{\int_{0}^{R}\int_{0}^{2\pi}\int_{0}^{\pi}\rho_{0}\left|\vec{v}\right|^{2}r^{2}\sin{\theta}\,\mathrm{d}\theta\,\mathrm{d}\phi\,\mathrm{d}r},$
(8)

with $E_{\mathrm{kin}}(\mathrm{core})$ being the mode kinetic energy in the
core and $E_{\mathrm{kin}}(\mathrm{core+RZ})$ as the mode kinetic energy in
the core and the radiative zone. The mode perturbation velocity is denoted as
$\vec{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 $\Delta
P_{\mathrm{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 $(\ell,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):

$\frac{1}{\Delta
P_{\mathrm{co}}}=x_{1}+x_{2}s+\frac{\Omega\sigma/\pi^{2}}{(s-s_{c})^{2}+\sigma^{2}},$
(9)

where $s$ is the spin parameter of the gravito-inertial modes in the co-
rotating frame and $x_{1}$, $x_{2}$, $s_{c}$, and $\sigma$ are free fitting
parameters. With this parametrisation, $x_{1}$ and $x_{2}$ account for the
general slope of the curve, $s_{c}$ is the spin parameter on which the dip is
centred, and $\sigma$ 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 $\sigma$ is always less
than 1% for the series ($\ell=1$, $m=-1$) and less than 2% for the series
($\ell=2$, $m=-2$), but it goes up to 10% with the rotation rate for the
series ($\ell=1$, $m=0$). The latter case could be explained by the fact that
$x_{1}$ and $x_{2}$ account for a linear slope of $1/\Delta P_{\mathrm{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 $\gamma$ 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

Figure 7: Period spacing in the co-rotating frame as a function of the spin
parameter for the modes of the series ($\ell=1$, $m=-1$) at different rotation
rates using the 1z model. The period spacing is normalised by the buoyancy
radius $\Pi_{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 $\ell_{i}=3$, $m=-1$ calculated
analytically with the uniform density model. Figure 8: Same as Fig. 7 but for
the model 2m. The low-amplitude oscillations in $\Delta 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).  Figure 9: Same as Fig. 7 but for the model
3t.  Figure 10: Evolution of the centre of the dip $s_{\mathrm{c}}$ 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.
Figure 11: Width of the dips $\sigma$ as a function of the rotation rate for
the series ($\ell=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 show in Figs. 7, 8, and 9 the evolution of the first dip of the series
($\ell=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 ($\ell=1$, $m=-1$) series are typically of a radial order $30<n<70$.
Analogous figures for the other two series, ($\ell=2$ $m=-2$) and ($\ell=1$,
$m=0$), are shown in Appendix C.

We fit the dips using Eq. (9) in order to determine their width $\sigma$ and
central spin parameter $s_{c}$. Figures 10 and 11 show the evolution of
$s_{c}$ and $\sigma$ 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
$\Omega=0.1\Omega_{\mathrm{K}}$ and $\Omega=0.15\Omega_{\mathrm{K}}$. We also
observed that $s_{c}$ decreases as the star evolves. The width $\sigma$
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.
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 $\sigma$ and $s_{c}$ and derive
approximate empirical relations usable on observational data.

###  4.1 Empirical expression of $\sigma$

According to the model of Tokuno & Takata (2022), $\sigma$ evolves linearly
with the parameter $\epsilon=\Omega/N_{0}$, where $N_{0}$ is the jump of the
Brunt-Väisälä frequency at the interface between the convective core and the
radiative region. We tested whether $\sigma$ is proportional to $\Omega$ by
plotting its evolution with $\Omega$ alongside the theoretical linear
evolution predicted in Tokuno & Takata (2022). Figure 12 shows the evolution
of $\sigma/\sigma_{\mathrm{ref}}$ with the rotation rate, where
$\sigma_{\mathrm{ref}}$ is the width computed at
$\Omega=\Omega_{\mathrm{ref}}=0.2\Omega_{\mathrm{K}}$. This figure shows that
$\sigma$ 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 $\sigma$. The origin of this particular
behaviour remains unclear at the moment, and numerical effects cannot be
excluded. Analogous figures for the series ($\ell=2$, $m=-2$) and ($\ell=1$,
$m=0$) are shown in Appendix C.

Figure 12: Evolution of $\sigma/\sigma_{\mathrm{ref}}$ with the rotation rate
for the dip of the series ($\ell=1$, $m=-1$) for three different evolutionary
stages (1z, 2m, 3t). The associated error bars are too small to be visible.
The reference value $\sigma_{\mathrm{ref}}$ is the width for
$\Omega=\Omega_{\mathrm{ref}}=0.2\Omega_{\mathrm{K}}$. If, as predicted by
Tokuno & Takata (2022), $\sigma\propto\Omega$, then
$\sigma/\sigma_{\mathrm{ref}}=\Omega/\Omega_{\mathrm{ref}}=5\Omega$ (the
dashed line).

Then, assuming a linear relation $\sigma=\beta\epsilon$, we performed a linear
regression to determine the slope $\beta$. We used the least-squares method.
Table 2 gathers the values of $\beta$ obtained for each dip and stellar model.
While $\beta$ 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 ($\ell=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 $\sigma=\beta\epsilon$, the
1z-model $\sigma$ values are weighted more in the determination of the free
parameter $\beta$ because, due to the smaller $N_{0}$ of the 1z model, they
extend over a larger $\epsilon$ range. To counter that, we actually fitted the
relation $\sigma N_{0}=\beta\Omega$. Analogous figures for the dips of the
($\ell=2$, $m=-2$) and ($\ell=1$, $m=0$) series are shown in Appendix C. Using
our empirical relation, the mean value of $\beta$ could thus be used to
estimate $N_{0}$ from $\sigma$ via

$N_{0}^{\mathrm{calc}}=\frac{\Omega\beta}{\sigma}.$ (10)

Table 2: Values of $\beta$.

$\ell$ | $m$ | 1z | 2m | 3t | All
---|---|---|---|---|---
1 | -1 | $5.03\pm 0.06$ | $4.63\pm 0.05$ | $4.33\pm 0.37$ | $4.66\pm 0.19$
1 | 0 | $0.59\pm 0.01$ | $0.64\pm 0.03$ | $0.60\pm 0.08$ | $0.61\pm 0.03$
2 | -2 | $2.89\pm 0.04$ | $2.83\pm 0.02$ | $2.55\pm 0.30$ | $2.76\pm 0.11$
222Values of $\beta$ are obtained by fitting the relation $\sigma
N_{0}=\beta\Omega$ for the series of modes ($\ell=1$, $m=-1$), ($\ell=1$,
$m=0$), and ($\ell=2$, $m=-2$) for the models 1z, 2m, and 3t separately and
together.

Figure 13: Width $\sigma$ as a function of $\epsilon$ for the series
($\ell=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 $\xi_{r}=0$, with $\xi_{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 $\Omega$ as long as an acoustic term can be neglected
in the governing equations (see Eq. (23)).

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
($\ell=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 ($\ell=2$, $m=-2$) and
($\ell=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 $\Omega=0.1\Omega_{\mathrm{K}}$ and
$\Omega=0.15\Omega_{\mathrm{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_{0}$ grows as the star evolves on the main sequence, which
makes the radial displacement $\xi_{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 $\epsilon$
usable on observation data, we noticed that though $s_{\rm c}$ does not vary
linearly with the rotation for all models (especially 3t), the difference
$s_{c}-s_{*}$ appears to be proportional to $\Omega$ (figure not shown). We
thus fit $s_{c}-s_{*}$ as a linear function of $\epsilon$:
$s_{c}-s_{*}=\alpha\epsilon$. We estimated $\alpha$ for the series of modes
($\ell=1$, $m=-1$), ($\ell=1$, $m=0$), and ($\ell=2$, $m=-2$) and for the
models 1z, 2m, and 3t separately. Table 3 provides the estimate of $\alpha$
for each case.

Table 3: Values of $\alpha$.

$\ell$ | $m$ | 1z | 2m | 3t
---|---|---|---|---
1 | -1 | $8.36\pm 0.24$ | $6.29\pm 0.58$ | $6.75\pm 0.49$
1 | 0 | $0.90\pm 0.01$ | $1.04\pm 0.03$ | $1.16\pm 0.05$
2 | -2 | $5.16\pm 0.06$ | $3.97\pm 0.40$ | $4.29\pm 0.49$
333Values of $\alpha$ are obtained by fitting the relation
$s_{c}(\epsilon)-s_{*}(\epsilon)=\alpha\epsilon$ for the series of modes
($\ell=1$, $m=-1$), ($\ell=1$, $m=0$), and ($\ell=2$, $m=-2$) for the models
1z, 2m, and 3t.

One can see that $\alpha$ varies significantly (up to $\sim$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 $\epsilon$ 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\epsilon+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
$\epsilon$ 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 $\beta$, we fitted
$(s_{c}/s_{*}-1)N_{0}=A\Omega$ so that all the models were weighted equally.
The result of the fit for $A$ is shown in Table 4.

Table 4: Values of $A$.

$\ell$ | $m$ | 1z | 2m | 3t | All
---|---|---|---|---|---
1 | -1 | $0.81\pm 0.02$ | $0.70\pm 0.06$ | $0.78\pm 0.06$ | $0.76\pm 0.04$
1 | 0 | $0.41\pm 0.01$ | $0.48\pm 0.01$ | $0.54\pm 0.02$ | $0.48\pm 0.03$
2 | -2 | $0.65\pm 0.01$ | $0.57\pm 0.05$ | $0.63\pm 0.07$ | $0.62\pm 0.03$
444Values of $A$ are obtained by fitting the relation
$(s_{c}/s_{*}-1)N_{0}=A\Omega$ for the series of modes ($\ell=1$, $m=-1$),
($\ell=1$, $m=0$), and ($\ell=2$, $m=-2$) for the models 1z, 2m, and 3t,
separately and together.

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

The global fit of all models together is plotted in Fig. 14 for the series
($\ell=1$, $m=-1$). Analogous figures for the ($\ell=2$, $m=-2$) and
($\ell=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:

$s_{*}^{\mathrm{calc}}=\frac{N_{0}^{\mathrm{calc}}s_{c}}{A\Omega+N_{0}^{\mathrm{calc}}}.$
(11)

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

$s_{*}^{\mathrm{calc}}=\frac{\beta s_{c}}{A\sigma+\beta}.$ (12)

Since $s_{c}$ and $\sigma$ are obtained from the $\Delta P_{\mathrm{co}}$-$s$
relation in the co-rotating frame, the rotation rate $\Omega$ is needed. It
can be determined by fitting the TAR outside the dip (e.g. Li et al., 2020).

Figure 14: Ratio $s_{c}/s_{*}$ as a function of $\Omega/N_{0}$ ($=\epsilon$)
for the series ($\ell=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.

###  4.3 Relative errors made on the estimations of $N_{0}$ 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_{0}$ and
$s_{*}$ using Eqs. (10) and (12), denoted as $N_{0}^{\mathrm{calc}}$ and
$s_{*}^{\mathrm{calc}}$. We then calculated the errors relative to the real
known values $N_{0}^{\mathrm{true}}$ and $s_{*}^{\mathrm{true}}$. We estimated
$s_{c}$ and $\sigma$ 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_{0}$ and $s_{*}$ as a function of
$\epsilon^{\mathrm{calc}}=\Omega/N_{0}^{\mathrm{calc}}$ for the series
($\ell=1$, $m=-1$). Analogous figures for the series ($\ell=2$, $m=-2$) and
($\ell=1$, $m=0$) are shown in Appendix C. One can see that for the ($\ell=1$,
$m=-1$) series, the relative error on $N_{0}$ 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 $\Omega=0.4,0.5\Omega_{\mathrm{K}}$, between 10 and 20% for the rotation
rates $\Omega=0.2,0.3\Omega_{\mathrm{K}}$, and between 30 and 50% for the
rotation rates $\Omega=0.1,0.15\Omega_{\mathrm{K}}$. Such behaviour is due to
the direct dependency between $N_{0}^{\mathrm{calc}}$ and $\sigma$. Indeed, as
$\epsilon$ gets smaller, $\sigma$ goes to zero, making relative errors explode
when we evaluate the difference between $\sigma$ and $\beta\epsilon$, despite
decreasing absolute errors. Nevertheless, it has to be noted that as typical
radial orders for observable modes of the ($\ell=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\approx 90$ at $\Omega=0.5\Omega_{\mathrm{K}}$ to
$n\approx 650$ at $\Omega=0.1\Omega_{\mathrm{K}}$. The same applies to the
($\ell=2$, $m=-2$) and ($\ell=1$, $m=0$) series, as the observed modes are
typically of about the same radial order as the ($\ell=1$, $m=-1$) series, but
the radial order of the modes calculated with the 3t model never goes lower
than $n\approx 170$ for the ($\ell=2$, $m=-2$) series and
$n\approx\leavevmode\nobreak\ 90$ for the ($\ell=1$, $m=0$) series.
The relative error of $s_{*}$ is far smaller and does not become greater than
$\sim 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 $\beta$ using only the 1z and 2m models (see the results in
Table 5). The relative error made on $N_{0}$ and $s_{*}$ for the series
($\ell=1$, $m=-1$) using these new estimates of $A$ and $\beta$ are shown in
Fig. 16. The relative error of $N_{0}$ is always inferior to 6%, and the
relative error on $s_{*}$ is inferior to 0.2%, except for the case of 1z and
$\Omega/\Omega_{\mathrm{K}}=0.5$, where it reaches $\sim$ 1%. The same
approach was used for the series ($\ell=1$, $m=0$) and ($\ell=2$, $m=-2$), and
the results are reported in Appendix C.

Figure 15: Relative errors for the series ($\ell=1$, $m=-1$) on the estimated
$N_{0}^{\mathrm{calc}}$ (upper panel) and $s_{*}^{\mathrm{calc}}$ (lower
panel) using the three models 1z, 2m, and 3t. Table 5: Values of $A$ and
$\beta$ when only taking into account the models 1z and 2m together.

$\ell$ | $m$ | $A$ | $\beta$
---|---|---|---
1 | -1 | $0.76\pm 0.05$ | $4.83\pm 0.14$
1 | 0 | $0.44\pm 0.03$ | $0.61\pm 0.02$
2 | -2 | $0.61\pm 0.04$ | $2.86\pm 0.03$
555Values of $A$ and $\beta$ are obtained by fitting the relations
$(s_{c}/s_{*}-1)N_{0}=A\Omega$ and $\sigma N_{0}=\beta\Omega$ for the series
of modes ($\ell=1$, $m=-1$), ($\ell=1$, $m=0$), and ($\ell=2$, $m=-2$).

Figure 16: Relative errors for the series ($\ell=1$, $m=-1$) on the estimated
$N_{0}^{\mathrm{calc}}$ (upper panel) and $s_{*}^{\mathrm{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_{0}$ 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 $\Delta
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, $(\ell_{i},m)$, of the Legendre polynomial involved in the solution
(see Appendix A). Under the TAR, the latitudinal part of gravito-inertial
$(\ell,m)$ modes is characterised by the Hough function $H_{k}^{m}(\theta,s)$
associated with the eigenvalue $\Lambda_{k}^{m}(s)$, where $k=\ell-|m|$ (see
Eq. 33). 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 $\gamma$ 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_{0}$ on the radiative side increases with age. According
to Tokuno & Takata (2022), the coupling between the $(\ell,m)$ series of the
gravito-inertial modes and a $(\ell_{i},m)$ inertial mode produces a
Lorentzian-shaped dip in the $\Delta P=f(s)$ that depends on two parameters,
namely its width,666The expressions of $\sigma$ and $s_{c}$ in the case
considered here are not explicitly written in Tokuno & Takata (2022) but can
easily be derived from their equations (89), (90), and (91) with
$\Delta\rho=0$ and replacing $\tilde{F}$ by $F$ given by their equation (58).

$\sigma=\frac{2P_{\ell_{i}}^{m}(1/s_{*})}{{\gamma^{\prime}}_{\ell_{i}}^{m}(s_{*})}\sqrt{\Lambda_{k}^{m}(s_{*})}\frac{\Omega}{N_{0}},$
(13)

and its central spin parameter,

$s_{c}=s_{*},$ (14)

where ${\gamma^{\prime}}_{\ell_{i}}^{m}$ is the derivative with respect to $s$
of the function

$\gamma_{\ell_{i}}^{m}(s)={P^{\prime}}_{\ell_{i}}^{m}(1/s)-\frac{m}{1-1/s^{2}}P_{\ell_{i}}^{m}(1/s),$
(15)

where ${P^{\prime}}_{\ell_{i}}^{m}$ is the derivative with respect to $x$ of
the Legendre polynomial ${P}_{\ell_{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
$\sigma$ and $s_{c}$ indicate that the dip parameters depend on the ratio
$\Omega/N_{0}$, on the gravito-inertial oscillations through
$\Lambda_{k}^{m}$, and on the inertial oscillation through $s_{*}$ and
$R(s_{*})=2P_{\ell_{i}}^{m}(1/s_{*})/{\gamma^{\prime}}_{\ell_{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 $\sigma$ 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 $\sigma\propto\Omega$ can be tested by considering the
evolution of dip widths with rotation for a given star model and $(\ell,m)$
series. Figure 12 shows that $\sigma\propto\Omega$ 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
$\sigma$ corresponding to the discontinuous case. A questionable property of
the theoretical $\sigma$ is that each inertial mode with the same $m$ and
equatorial parity as the gravito-inertial modes of a $(\ell,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 $\Delta 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 $\Omega/N_{0}\ll 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 $(\ell=1,m=-1)$ versus $(\ell_{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.

Figure 17: Latitudinal profiles of the $\ell_{i}=3,m=-1,s_{*}=11.32$ inertial
mode ($P_{3}^{-1}(\cos\theta)$, blue) and $\ell=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
$(\ell=1,m=-1),(\ell=2,m=-2),(\ell=1,m=0)$ gravito-inertial modes and
respectively the $(\ell_{i}=3,m=-1),(\ell_{i}=4,m=-2),(\ell_{i}=3,m=0)$
inertial modes produces Lorentzian-shaped dips in the $\Delta P=f(s)$ relation
whose widths and central spin parameters are described by

$\displaystyle\sigma=$
$\displaystyle\;\left(\frac{f_{k,\ell_{i}}^{m}(s_{*})}{1+f_{k,\ell_{i}}^{m}(s_{*})}\frac{2P_{\ell_{i}}^{m}(1/s_{*})}{{\gamma^{\prime}}_{\ell_{i}}^{m}(s_{*})}\sqrt{\Lambda_{k}^{m}(s_{*})}\right)\frac{\Omega}{N_{0}},$
(16) $\displaystyle s_{c}=$ $\displaystyle
s_{*}+\;\left(\frac{1}{1+f_{k,\ell_{i}}^{m}(s_{*})}\frac{2P_{\ell_{i}}^{m}(1/s_{*})}{{\gamma^{\prime}}_{\ell_{i}}^{m}(s_{*})}\sqrt{\Lambda_{k-2}^{m}(s_{*})}\right)\frac{\Omega}{N_{0}},$
(17)

where $f_{k,\ell_{i}}^{m}$ is defined by

$f_{k,\ell_{i}}^{m}(s_{*})=-\frac{c_{k,\ell_{i}}^{m}(s_{*})c_{k-2,\ell_{i}-2}^{m}(s_{*})}{c_{k,\ell_{i}-2}^{m}(s_{*})c_{k-2,\ell_{i}}^{m}(s_{*})}.$
(18)

Here, $c^{m}_{k,\ell_{i}}$ characterises the geometrical coupling between
inertial and gravito-inertial oscillations along the convective-radiative
interface:

$c^{m}_{k,\ell_{i}}(s)=\int_{0}^{\pi}\tilde{P}_{\ell_{i}}^{m}(\cos\theta)H_{k}^{m}(\theta,s)\sin\theta\;\mathrm{d}\theta,$
(19)

with $\tilde{P}_{\ell_{i}}^{m}(\cos\theta)$ and $H_{k}^{m}(\theta,s)$ being
the normalised Legendre polynomial and Hough functions. In addition to
$P_{\ell_{i}}^{m}$ and $H_{k}^{m}$, these mixed modes involve
$P_{\ell_{i}-2}^{m}$ and $H_{k-2}^{m}$. In Appendix D, we argue this is a good
approximation for the first dip of the $(\ell=1,m=-1)$, $(\ell=2,m=-2)$, and
$(\ell=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 $\sigma$ is multiplied by the factor
$f_{k,\ell_{i}}^{m}/(1+f_{k,\ell_{i}}^{m})$. As $f_{k,\ell_{i}}^{m}$ is
proportional to the coupling coefficient $c_{k,\ell_{i}}(s_{*})$ introduced in
Ouazzani et al. (2020), this new expression of $\sigma$ 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_{*}\propto\Omega$.

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

$\displaystyle\sigma=$ $\displaystyle\beta\frac{\Omega}{N_{0}},$ (20)
$\displaystyle s_{c}=$ $\displaystyle s_{*}+\alpha\frac{\Omega}{N_{0}},$ (21)

where $\beta$, $s_{*}$, and $\alpha$ depend on both the gravito-inertial modes
and on the resonant inertial mode. For a fixed $(\ell,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 $\beta$ remains roughly constant,
while $\alpha$ is approximately proportional to $1/s_{*}$. An open question
for future and more detailed empirical relations is whether $\beta$ and
$\alpha$ 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 $(\ell,m)$ mode
series.

Because of the uniform density assumption, a quantitative agreement between
the theoretical and numerical values of $\sigma$ and $s_{c}$ is not expected.
We nevertheless compared them by calculating the theoretical values of
$\beta_{\rm th}(s_{*})$, $\alpha_{\rm th}(s_{*})$, and $A_{\rm th}=\alpha_{\rm
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 $\alpha$, the deviations from the uniform density values
increase for more evolved star models with higher density contrasts.

Table 6: Parameters of the empirical relations for $\sigma$ and $s_{c}$ proposed in Sect. 4 using a uniform density theoretical model. $\ell$ | $m$ | $\alpha_{\rm th}$ | $A_{\rm th}$ | $\beta_{\rm th}$
---|---|---|---|---
1 | -1 | $8.99$ | $0.794$ | $4.78$
1 | 0 | $0.769$ | $0.344$ | $0.599$
2 | -2 | $5.77$ | $0.67$ | $2.74$

##  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
$\gamma$ Dor stars. As in Ouazzani et al. (2020), we obtained dips in the
$\Delta P_{\mathrm{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 $\gamma$ 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
$\gamma$ Dor stars, and similar dip structures could arise from the $\Delta 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
$\gamma$ 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 $\gamma$ 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 $\beta$ 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 $\gamma$ 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 $\Delta 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. (23), 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 $\gamma$ Dor stars as scientific calibrators.

###### Acknowledgements.

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

  * Aerts (2021) Aerts, C. 2021, Reviews of Modern Physics, 93, 015001
  * Aerts & Mathis (2023) Aerts, C. & Mathis, S. 2023, A&A, 677, A68
  * Aerts et al. (2019) Aerts, C., Mathis, S., & Rogers, T. M. 2019, ARA&A, 57, 35
  * Angulo et al. (1999) Angulo, C., Arnould, M., Rayet, M., et al. 1999, Nucl. Phys. A, 656, 3
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481
  * Ballot et al. (2012) Ballot, J., Lignières, F., Prat, V., Reese, D. R., & Rieutord, M. 2012, in Astronomical Society of the Pacific Conference Series, Vol. 462, Progress in Solar/Stellar Physics with Helio- and Asteroseismology, ed. H. Shibahashi, M. Takata, & A. E. Lynas-Gray, 389
  * Ballot et al. (2010) Ballot, J., Lignières, F., Reese, D. R., & Rieutord, M. 2010, A&A, 518, A30
  * Böhm-Vitense (1958) Böhm-Vitense, E. 1958, ZAp, 46, 108
  * Borucki et al. (2010) Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977
  * Bouabid et al. (2013) Bouabid, M.-P., Dupret, M.-A., Salmon, S., et al. 2013, Monthly Notices of the Royal Astronomical Society, 429, 2500
  * Chatelin (1988) Chatelin, F. 1988, Valeurs propres de matrices, Collection Mathématiques appliquées pour la maîtrise (Masson)
  * Christophe et al. (2018) Christophe, S., Ballot, J., Ouazzani, R. M., Antoci, V., & Salmon, S. J. A. J. 2018, A&A, 618, A47
  * Dupret et al. (2002) Dupret, M. A., De Ridder, J., Neuforge, C., Aerts, C., & Scuflaire, R. 2002, A&A, 385, 563
  * Dupret et al. (2005) Dupret, M. A., Grigahcène, A., Garrido, R., Gabriel, M., & Scuflaire, R. 2005, A&A, 435, 927
  * Ferguson et al. (2005) Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585
  * Garcia et al. (2022) Garcia, S., Van Reeth, T., De Ridder, J., & Aerts, C. 2022, A&A, 668, A137
  * Guzik et al. (2000) Guzik, J. A., Kaye, A. B., Bradley, P. A., Cox, A. N., & Neuforge, C. 2000, ApJ, 542, L57
  * Iglesias & Rogers (1996) Iglesias, C. A. & Rogers, F. J. 1996, ApJ, 464, 943
  * Imbriani et al. (2004) Imbriani, G., Costantini, H., Formicola, A., et al. 2004, A&A, 420, 625
  * Kaye et al. (1999) Kaye, A. B., Handler, G., Krisciunas, K., Poretti, E., & Zerbi, F. M. 1999, PASP, 111, 840
  * Lebreton et al. (2014) Lebreton, Y., Goupil, M. J., & Montalbán, J. 2014, in EAS Publications Series, Vol. 65, EAS Publications Series, ed. Y. Lebreton, D. Valls-Gabaud, & C. Charbonnel, 177–223
  * Lee & Saio (1989) Lee, U. & Saio, H. 1989, MNRAS, 237, 875
  * Lee & Saio (1997) Lee, U. & Saio, H. 1997, ApJ, 491, 839
  * Li et al. (2019) Li, G., Van Reeth, T., Bedding, T. R., Murphy, S. J., & Antoci, V. 2019, MNRAS, 487, 782
  * Li et al. (2020) Li, G., Van Reeth, T., Bedding, T. R., et al. 2020, Monthly Notices of the Royal Astronomical Society, 491, 3586
  * Miglio et al. (2008) Miglio, A., Montalbán, J., Noels, A., & Eggenberger, P. 2008, MNRAS, 386, 1487
  * Mombarg et al. (2021) Mombarg, J. S. G., Van Reeth, T., & Aerts, C. 2021, A&A, 650, A58
  * Mombarg et al. (2019) Mombarg, J. S. G., Van Reeth, T., Pedersen, M. G., et al. 2019, MNRAS, 485, 3248
  * Morel (1997) Morel, P. 1997, A&AS, 124, 597
  * Morel & Lebreton (2008) Morel, P. & Lebreton, Y. 2008, Ap&SS, 316, 61
  * Ouazzani et al. (2020) Ouazzani, R.-M., Lignières, F., Dupret, M.-A., et al. 2020, A&A, 640, A49
  * Ouazzani et al. (2019) Ouazzani, R. M., Marques, J. P., Goupil, M. J., et al. 2019, A&A, 626, A121
  * Ouazzani et al. (2017) Ouazzani, R.-M., Salmon, S. J. A. J., Antoci, V., et al. 2017, MNRAS, 465, 2294
  * Rauer et al. (2014) Rauer, H., Catala, C., Aerts, C., et al. 2014, Experimental Astronomy, 38, 249
  * Reese (2006) Reese, D. 2006, PhD thesis, Universite de Toulouse Paul Sabatier, France
  * Reese (2013) Reese, D. R. 2013, A&A, 555, A148
  * Reese et al. (2009) Reese, D. R., MacGregor, K. B., Jackson, S., Skumanich, A., & Metcalfe, T. S. 2009, A&A, 506, 189
  * Ricker et al. (2014) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9143, Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave, ed. J. Oschmann, Jacobus M., M. Clampin, G. G. Fazio, & H. A. MacEwen, 914320
  * Rogers & Nayfonov (2002) Rogers, F. J. & Nayfonov, A. 2002, ApJ, 576, 1064
  * Saio et al. (2018) Saio, H., Kurtz, D. W., Murphy, S. J., Antoci, V. L., & Lee, U. 2018, MNRAS, 474, 2774
  * Saio et al. (2021) Saio, H., Takata, M., Lee, U., Li, G., & Van Reeth, T. 2021, Monthly Notices of the Royal Astronomical Society, 502, 5856
  * Takata et al. (2020) Takata, M., Ouazzani, R. M., Saio, H., et al. 2020, A&A, 635, A106
  * Tokuno & Takata (2022) Tokuno, T. & Takata, M. 2022, MNRAS, 514, 4140
  * Townsend (2003) Townsend, R. 2003, Monthly Notices of the Royal Astronomical Society, 340, 1020
  * Unno et al. (1989) Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial oscillations of stars
  * Uytterhoeven et al. (2011) Uytterhoeven, K., Moya, A., Grigahcène, A., et al. 2011, A&A, 534, A125
  * Van Reeth et al. (2018) Van Reeth, T., Mombarg, J. S. G., Mathis, S., et al. 2018, A&A, 618, A24
  * Van Reeth et al. (2016) Van Reeth, T., Tkachenko, A., & Aerts, C. 2016, A&A, 593, A120
  * Van Reeth et al. (2015) Van Reeth, T., Tkachenko, A., Aerts, C., et al. 2015, ApJS, 218, 27
  * Wu (2005) Wu, Y. 2005, ApJ, 635, 674

##  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 $\nabla^{2}\Psi-s^{2}\frac{\partial^{2}\Psi}{\partial
s^{2}}$=0, the Poincaré equation. It is separable using the ellipsoidal
coordinates $(x_{1},x_{2},\phi)$ (see Eqs. 24, 25, 26).

Using regularity conditions and the boundary condition $\xi_{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

$\frac{dP_{\ell_{i}}^{m}(\mu)}{d\mu}=\frac{m}{1-\mu^{2}}P_{\ell_{i}}^{m}(\mu),$
(22)

where $\mu=1/s=\omega_{\mathrm{co}}/2\Omega$ and $\ell_{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 7. 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 7.
The modes $(\ell_{i}=3$ $m=-1$), $(\ell_{i}=3$ $m=0$), and $(\ell_{i}=4$
$m=-2$) with the highest correlation coefficients are shown in Table 7.

Table 7: Spin parameters of inertial modes obtained using an analytical uniform density core model. The correlation coefficient between the inertial mode and the Hough function of the same symmetry class and spin parameter is shown in parentheses. $\ell_{i}$ | $m=0$ |  |  |
---|---|---|---|---
1 | $\emptyset$ |  |  |
3 |  $2.236$ $(0.62)$ |  |  |
5 |  $3.506$ $(0.51)$ |  $1.307$ $(0.014)$ |  |
7 |  $4.778$ $(0.44)$ |  $1.690$ $(3.3\hskip 3.0pt10^{-3})$ |  $1.147$ $(1.6\hskip 3.0pt10^{-4})$ |
$\ell_{i}$ | $m=-1$ |  |  |
1 | $\emptyset$ |  |  |
3 |  $11.325$ $(0.50)$ |  |  |
5 |  $29.330$ $(0.39)$ |  $1.690$ $(4.3\hskip 3.0pt10^{-3})$ |  |
7 |  $55.332$ $(0.34)$ |  $2.318$ $(3.8\hskip 3.0pt10^{-5})$ |  $1.293$ $(1.89\hskip 3.0pt10^{-5})$ |
$\ell_{i}$ | $m=-2$ |  |  |
2 | $\emptyset$ |  |  |
4 |  $8.624$ $(0.48)$ |  |  |
6 |  $19.648$ $(0.38)$ |  $1.831$ $(1.6\hskip 3.0pt10^{-2})$ |  |
8 |  $34.662$ $(0.32)$ |  $2.482$ $(3.3\hskip 3.0pt10^{-3})$ |  $1.374$ $(2.4\hskip 3.0pt10^{-4})$ |

##  Appendix B Core-only model

We run TOP with a truncated version of our models to obtain the resonant
frequencies. Table 8 shows the spin parameters of the inertial modes coupling
with the series of gravito-inertial modes ($\ell=1$ $m=-1$), ($\ell=1$ $m=0$),
and ($\ell=2$ $m=-2$) for the 3 models. They are labelled as ($\ell_{i}=3$
$m=-1$), ($\ell_{i}=3$ $m=0$), and ($\ell_{i}=4$ $m=-2$) respectively as they
remain similar to the modes calculated with the uniform density model
presented in Appendix A (see Table 7). 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. 18.

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

$\nabla^{2}\Psi-s^{2}\frac{\partial^{2}\Psi}{\partial
z^{2}}=-\frac{1}{\rho_{0}}\frac{\mathrm{d}\rho_{0}}{\mathrm{d}r}\left(\frac{\partial\Psi}{\partial
r}-s^{2}\cos{\theta}\frac{\partial\Psi}{\partial
z}+\frac{ms}{r}\Psi\right)-(1-s^{2})\frac{\omega_{\mathrm{co}}^{2}}{c_{0}^{2}}\Psi$
(23)

where $\Psi=p^{\prime}/(\rho_{0}\omega_{\mathrm{co}}^{2})$ and
$z=r\cos{\theta}$ 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 $\gamma$ Dor stars.

Table 8: Spin parameters $s_{*}$ of the resonant inertial modes coupling with the series of modes ($\ell=1$ $m=-1$), ($\ell=1$ $m=0$), and ($\ell=2$ $m=-2$), for the models 1z, 2m, and 3t using a core-only model at $\Omega=0.1\Omega_{K}$ and $\Omega=0.5\Omega_{K}$.  $\ell_{i}$ | $m$ | 1z | 2m | 3t
---|---|---|---|---
3 | -1 | 10.289 - 10.290 | 8.966 - 8.964 | 8.712 - 8.660
3 | 0 | 2.207 - 2.207 | 2.163 - 2.163 | 2.153 - 2.152
4 | -2 | 7.920 - 7.921 | 6.980 - 6.978 | 6.805 - 6.755
Table 9: Inertial modes in the core coupling with the ($\ell=1$ $m=-1$),
($\ell=1$ $m=0$), and ($\ell=2$ $m=-2$) series of gravito-inertial modes in
the envelope for the uniform density and 3t models.

$\ell_{i}$ | $m$ |  Uniform density | 3t
---|---|---|---
3 | -1 |  |
3 | 0 |  |
4 | -2 |  |
777The white line shows the critical latitude $\theta_{c}=\arccos(1/s)$ in the
case of the uniform density model. $\Omega=0.3\Omega_{K}$.

Figure 18: Eulerian pressure perturbations normalised by their value at the
core radius $R_{\mathrm{c}}$ of the ($\ell_{i}=3$, $m=-1$) mode, as a function
of the radius at the colatitude $\theta=\pi/2$ (equator), for the uniform
density model (blue) and for the 3t model (red).

##  Appendix C Results and analysis for the ($\ell=1$, $m=0$) and ($\ell=2$,
$m=-2$) series

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

###  C.1 ($\ell=1$, $m=0$) series

The dips in the series of modes ($\ell=1$, $m=0$) computed for the models 1z,
2m, and 3t at different rotations are presented in Figs. 19, 20, and 21.

Figure 19: Period spacing in the co-rotating frame as a function of the spin
parameter for the modes of the series ($\ell=1$, $m=0$) at different rotation
rates using the 1z model. The period spacing is normalised by the buoyancy
radius $\Pi_{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). Figure 20: Same as Fig. 19 for the model 2m. Figure 21: Same as
Fig. 19 for the model 3t.

The evolution of the width $\sigma$ and location $s_{c}$ of dips as a function
of the rotation are shown in Figs. 22 and 23. 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.

Figure 22: Width of the dips $\sigma$ as a function of the rotation rate for
the series ($\ell=1$, $m=0$), for the three stellar models described in Table
1. Figure 23: Evolution of the centre of the dip $s_{\mathrm{c}}$ 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 $\sigma$ and $s_{c}/s_{*}$ as a function of $\epsilon=\Omega/N_{0}$
for the three models in Figs. 24 and 25, along with linear fits of the form
$\sigma=\beta\epsilon$ and $s_{c}/s_{*}=A\epsilon+1$. We can see that $\sigma$
is a lot smaller for the ($\ell=1$, $m=0$) series than it is for the
($\ell=1$, $m=-1$), which is consistent with our analytical model (see the
$\beta_{\mathrm{th}}$ column of Table 6). The ratio $s_{c}/s_{*}$ is also
smaller, which is in agreement with our analytical model ($A_{\mathrm{th}}$
column of Table 6).

Figure 24: Width $\sigma$ as a function of $\epsilon$ for the series ($\ell=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.

.

Figure 25: Ratio $s_{c}/s_{*}$ as a function of $\Omega/N_{0}$ ($\epsilon$)
for the series ($\ell=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. 26, we compare the relative evolution of $\sigma$ with rotation in our
numerical computations and in the theoretical model of Tokuno & Takata (2022).
Compared to the ($\ell=1$, $m=-1$) series, the non-linearity of $\sigma$ 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.

Figure 26: Evolution of $\sigma/\sigma_{\mathrm{ref}}$ with the rotation rate,
for the dip of the series ($\ell=1$, $m=0$) for three different evolutionary
stages (1z, 2m, 3t), with the associated error bars. The reference value
$\sigma_{\mathrm{ref}}$ is the width for
$\Omega=\Omega_{\mathrm{ref}}=0.2\Omega_{\mathrm{K}}$. If, as predicted by
Tokuno & Takata (2022), $\sigma\propto\Omega$,
$\sigma/\sigma_{\mathrm{ref}}=\Omega/\Omega_{\mathrm{ref}}=5\Omega$ (the
dashed line).

We used the empirical relations (10) and (12) on our numerical calculations to
estimate $s_{*}$ and $N_{0}$ and compared them to their true values. The
errors made are shown in Figs. 27 and 28. 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 $\ell=|m|$ (Kelvin) modes. The estimate
for $N_{0}$ is bad for the 3t model regardless of the rotation rate, making
the estimate of $s_{*}$ also worse for the 3t model.

Figure 27: Relative errors for the series ($\ell=1$ $m=0$) on the estimated
$N_{0}^{\mathrm{calc}}$ (upper panel) and $s_{*}^{\mathrm{calc}}$ (lower
panel) using the three models 1z, 2m, and 3t. Figure 28: Relative errors for
the series ($\ell=1$ $m=0$) on the estimated $N_{0}^{\mathrm{calc}}$ (upper
panel) and $s_{*}^{\mathrm{calc}}$ (lower panel) using the two models 1z and
2m.

###  C.2 ($\ell=2$, $m=-2$) series

The dips in the series of modes ($\ell=2$, $m=-2$) computed for the models 1z,
2m, and 3t at different rotations are presented in Figs. 29, 30 and 31.

Figure 29: Same as Fig. 19 for the series ($\ell=2$, $m=-2$). Figure 30: Same
as Fig. 20 for the series ($\ell=2$, $m=-2$). Figure 31: Same as Fig. 21 for
the series ($\ell=2$, $m=-2$).

The evolution of the width $\sigma$ and location $s_{c}$ of dips as a function
of the rotation are shown in Figs. 32 and 33. 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.

Figure 32: Same as Fig. 22 for the series ($\ell=2$, $m=-2$). Figure 33: Same
as Fig. 23 for the series ($\ell=2$, $m=-2$).

We plot $\sigma$ and $s_{c}/s_{*}$ as a function of $\epsilon=\Omega/N_{0}$
for the three models in Figs. 34 and 35, along with linear fits of the form
$\sigma=\beta\epsilon$ and $s_{c}/s_{*}=A\epsilon+1$. At a given $\epsilon$,
both $\sigma$ and $s_{c}/s_{*}$ are slightly smaller for the ($\ell=2$,
$m=-2$) series compared to the ($\ell=1$, $m=-1$) series. This behaviour is
consistent with our analytical model (see the $\beta_{\mathrm{th}}$ and
$A_{\mathrm{th}}$ column of Table 6).

Figure 34: Same as Fig. 24 for the series ($\ell=2$, $m=-2$). Figure 35: Same
as Fig. 25 for the series ($\ell=2$, $m=-2$).

In Fig. 36, we compare the relative evolution of $\sigma$ with rotation in our
numerical computations and in the theoretical model of Tokuno & Takata (2022).
As for the series ($\ell=1$, $m=-1$, while $\sigma$ 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 ($\ell=2$, $m=-2$) series than it is for
the ($\ell=1$, $m=-1$) series.

Figure 36: Same as Fig. 26 for the series ($\ell=2$, $m=-2$).

We used the empirical relations (10) and (12) on our numerical calculations to
estimate $s_{*}$ and $N_{0}$ and compared them to their true values. The
errors made are shown in Figs. 37 and 38. 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_{0}$ and $s_{*}$ follow the same kind of trend as for
the ($\ell=1$, $m=-1$) series.

Figure 37: Same as Fig. 27 for the series ($\ell=2$, $m=-2$). Figure 38: Same
as Fig. 28 for the series ($\ell=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 $(\ell,m)=(1,-1)$ gravito-inertial modes
interacting with the $(\ell_{i},m)=(3,-1),s_{*}=11.32$ inertial mode. The
effect of this mode coupling on the $\Delta 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\phi+\omega t)}$ where $m$ is the azimuthal number. The
radial displacement and the Eulerian pressure perturbations are thus expressed
as $\xi^{\rm conv}_{r}(r,\theta)e^{i(m\phi+\omega t)}$ and ${p^{\prime}}^{\rm
conv}(r,\theta)e^{i(m\phi+\omega t)}$ respectively. In an isentropic
convective core of uniform density, perturbations of low frequency are
governed by the Poincaré equation
$\Delta\Psi=s^{2}\frac{\partial\Psi}{\partial z}$ where $\Psi$ is proportional
to ${p^{\prime}}^{\rm conv}$ (Wu 2005). In a full sphere, separable solutions
are obtained using the ellipsoidal coordinates ($x_{1},x_{2},\phi$) defined in
terms of the Cartesian coordinates as

$\displaystyle x/r_{c}$ $\displaystyle=$
$\displaystyle\left(\frac{(1-x_{1}^{2})(1-x_{2}^{2})}{1-\mu_{s}^{2}}\right)^{1/2}\cos(\phi)$
(24) $\displaystyle y/r_{c}$ $\displaystyle=$
$\displaystyle\left(\frac{(1-x_{1}^{2})(1-x_{2}^{2})}{1-\mu_{s}^{2}}\right)^{1/2}\sin(\phi)$
(25) $\displaystyle z/r_{c}$ $\displaystyle=$
$\displaystyle\frac{x_{1}x_{2}}{\mu_{s}}$ (26)

where $\mu_{s}=1/s$, $x_{1}\in[\mu_{s},1]$, and $x_{2}\in[-\mu_{s},\mu_{s}]$.
With these coordinates, the surface of the sphere is given by
$\\{x_{1}\in[\mu_{s},1],x_{2}=\mu_{s}\\}\cup\\{x_{1}=\mu_{s},x_{2}\in[-\mu_{s},\mu_{s}]\\}$.

The separable solutions read:

$\displaystyle\Psi=r_{c}^{2}P_{\ell_{i}}^{m}(x_{1})P_{\ell_{i}}^{m}(x_{2})$
(27)

where $P_{\ell_{i}}^{m}$ is the associated Legendre polynomial of degree
$\ell_{i}$. The radial displacement $\xi^{\rm conv}_{r}$ and the Eulerian
pressure perturbations $p^{\prime\rm conv}$ are then obtained from the
relations (Wu 2005):

$\displaystyle\xi^{\rm conv}_{r}=$
$\displaystyle\frac{1}{1-s^{2}}\;\left(\frac{\partial\Psi}{\partial
r}+\frac{ms}{r}\Psi-s^{2}\cos\theta\frac{\partial\Psi}{\partial z}\right)$
(28) $\displaystyle p^{\prime\rm conv}=$
$\displaystyle\rho_{0}\omega^{2}\Psi.$ (29)

At the outer radius of the convective core $r=r_{c}$, it reads (Tokuno &
Takata 2022):

$\displaystyle\xi^{\rm conv}_{r}(r=r_{c})=$ $\displaystyle
r_{c}\frac{1}{s}\gamma_{\ell_{i}}^{m}(s)P_{\ell_{i}}^{m}(\cos\theta)\;\;\;{\rm
with}$ (30) $\displaystyle\gamma_{\ell_{i}}^{m}(s)=$
$\displaystyle{P^{\prime}}_{\ell_{i}}^{m}(1/s)-\frac{m}{1-1/s^{2}}P_{\ell_{i}}^{m}(1/s)$
(31)

where ${P^{\prime}}_{\ell_{i}}^{m}$ is the derivative of the polynomial
${P}_{\ell_{i}}^{m}$, and

$\displaystyle p^{\prime\rm conv}(r=r_{c})=$
$\displaystyle\rho_{0}r_{c}^{2}\omega^{2}P_{\ell_{i}}^{m}(1/s)P_{\ell_{i}}^{m}(\cos\theta).$
(32)

These solutions oscillate freely at the convective-radiative interface because
no boundary condition is imposed there. If the $\xi^{\rm conv}_{r}(r=r_{c})=0$
boundary condition is enforced, one recovers the eigenfrequencies of pure
inertial modes presented in Appendix A which satisfy Eq. 22 or equivalently
$\gamma_{\ell_{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^{\rm rad}_{r}(r,\theta)=\zeta(r)H(\theta)$ and ${p^{\prime}}^{\rm
rad}(r,\theta)=p^{\prime}(r)H(\theta)$, where $\zeta$ and $p^{\prime}$ denote
the radial part of the radial displacement and Eulerian pressure perturbation
respectively, and $H(\theta)$ their latitudinal part.

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

${\cal L_{\rm T}}(H)+\Lambda H=0$ (33)

with

$\displaystyle{\cal L_{\rm T}}(H)=$
$\displaystyle\frac{d}{d\mu}\left[\frac{1-\mu^{2}}{(1-s^{2}\mu^{2})}\frac{d}{d\mu}H\right]$
$\displaystyle-\frac{1}{1-s^{2}\mu^{2}}\left(\frac{m^{2}}{1-\mu^{2}}+ms\frac{1+s^{2}\mu^{2}}{1-s^{2}\mu^{2}}\right)H,$

where $\mu=\cos\theta$ and $\Lambda$ 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
$\Lambda^{m}_{k}(s)$ and $H^{m}_{k}(\theta,s)$ respectively. They form an
orthogonal set of functions as
$\int_{0}^{\pi}H^{m}_{k}(\theta,s)H^{m}_{k^{\prime}}(\theta,s)\sin\theta\;d\theta$=0
if $k\neq k^{\prime}$. In the following, we shall use normalised Hough
functions such that
$\int_{0}^{\pi}\;\left(H^{m}_{k}(\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 $\ell$
of the spherical harmonics of the gravity mode in a non-rotating star through
$\ell=|m|+k$. Thus, for example, the prograde sectoral modes, also called
Kelvin modes, can be equivalently labelled $(\ell=|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 $\Lambda<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 $\Lambda<0$ non
propagative modes.

Within the Cowling approximation, the equations governing the radial part of
the oscillation, $\zeta$ and $p^{\prime}$, are:

$\displaystyle\frac{d\zeta}{dr}=$
$\displaystyle-\;\left(\frac{2}{r}-\frac{1}{\Gamma_{1}H_{P}}\right)\zeta+\frac{1}{\rho_{0}c_{0}^{2}}\;\left(\frac{\Lambda
c_{0}^{2}}{r^{2}\omega^{2}}-1\right)p^{\prime}$ (34)
$\displaystyle\frac{dp^{\prime}}{dr}=$
$\displaystyle\rho_{0}(\omega^{2}-N^{2})\xi_{r}-\frac{1}{\Gamma_{1}H_{P}}p^{\prime}$
(35)

where $\Lambda^{m}_{k}(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}}\;\left(1-\frac{\Lambda_{k}^{m}(s)c_{0}^{2}}{r^{2}\omega^{2}}\right)$
:

$\displaystyle{\rm if}\;k_{r}^{2}>0$
$\displaystyle\;\;\zeta_{k}^{m}(r,s)=-c\frac{|1-\frac{\Lambda
c_{0}^{2}}{r^{2}\omega^{2}}|^{1/2}}{rc_{0}(\pi\rho_{0}k_{r})^{1/2}}\sin\\!\;\left(\int_{r}^{r_{b}}k_{r}\;dr-\frac{\pi}{4}\right)$
(36)
$\displaystyle\;\;\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\\!\;\left(\int_{r}^{r_{b}}k_{r}\;dr-\frac{\pi}{4}\right)$
(37) $\displaystyle{\rm if}\;k_{r}^{2}<0$
$\displaystyle\;\;\zeta_{k}^{m}(r,s)=c^{\prime}\frac{\;\left(1-\frac{\Lambda
c_{0}^{2}}{r^{2}\omega^{2}}\right)^{1/2}}{rc_{0}(\pi\rho_{0}|k_{r}|)^{1/2}}\exp\\!\;\left(-\int_{r_{c}}^{r}|k_{r}|\;dr\right)$
(38)
$\displaystyle\;\;\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\\!\;\left(-\int_{r_{c}}^{r}|k_{r}|\;dr\right)$
(39)

where $c$ and $c^{\prime}$ are arbitrary constants and the labels ${k,m,s}$ of
$\Lambda$ 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{\Lambda_{k}^{m}(s)c_{0}(r_{b})^{2}}{r_{b}^{2}\omega^{2}}\approx 1$ as
$\omega\ll N$. The solution exponentially decreasing from $r_{c}$ corresponds
to Hough modes with negative $\Lambda$.

These oscillations are free at the convective/radiative interface. If a
boundary condition $\xi_{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. 36, this leads to the quantisation
condition $\int_{r_{c}}^{r_{b}}k_{r}\;dr-\pi/4=n\pi$. As $\Lambda
c_{0}^{2}/(r^{2}\omega^{2})\gg 1$ away from $r_{b}$ and $\omega\ll N$, the
condition can be simplified into
$\frac{\pi^{2}s\sqrt{\Lambda}}{\Omega\Pi_{0}}=\pi(n+1/4)$ with
$\Pi_{0}=\frac{2\pi^{2}}{\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 $\xi_{r}$ and
$\delta p$ have to be continuous there. However, as illustrated in Fig. 17 for
the $(\ell=1,m=-1)$ vs $(\ell_{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 :

$\displaystyle{\xi}^{\rm rad}_{r}=$
$\displaystyle\sum_{k=-\infty}^{k=+\infty}a_{k}\zeta_{k}^{m}(r,s)H_{k}^{m}(\theta,s)$
(40) $\displaystyle p^{\prime\rm rad}=$
$\displaystyle\sum_{k=-\infty}^{k=+\infty}a_{k}{p^{\prime}}_{k}^{m}(r,s)H_{k}^{m}(\theta,s)$
(41)

and in the convective zone :

$\displaystyle\Psi\propto$ $\displaystyle
r_{c}^{2}\sum_{\ell_{i}=1}^{+\infty}b_{\ell_{i}}P_{\ell_{i}}^{m}(x_{1})P_{\ell_{i}}^{m}(x_{1}).$
(42)

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 $\xi_{r}$ and $\delta p$ is equivalent to
the continuity of $\xi_{r}$ and $p^{\prime}$, thus to ${\xi}^{\rm
conv}_{r}(r=r_{c}^{-})={\xi}^{\rm rad}_{r}(r=r_{c}^{+})$ and $p^{\prime\rm
conv}(r=r_{c}^{-})=p^{\prime\rm rad}(r=r_{c}^{+})$.

The expressions of ${\xi}^{\rm rad}_{r}(r_{c}^{+})$ and $p^{\prime\rm
rad}(r_{c}^{+})$ are

$\displaystyle{\xi}^{\rm rad}_{r}(r_{c}^{+})=$
$\displaystyle\sum_{k=-\infty}^{k=+\infty}a_{k}\zeta_{k}^{m}(r_{c}^{+},s)H_{k}^{m}(\theta,s)\;\;\;{\rm
with}$ (43) $\displaystyle\zeta_{k}^{m}(r_{c}^{+},s)=$
$\displaystyle-\frac{1}{\rho_{0}^{1/2}\omega
r_{c}^{3}}\;\left(\frac{\Lambda}{\pi k_{r}}\right)^{1/2}\sin(\psi)\;\;\;{\rm
if}\;\;\;k_{r}^{2}>0$ (44) $\displaystyle\zeta_{k}^{m}(r_{c}^{+},s)=$
$\displaystyle\frac{1}{\rho_{0}^{1/2}\omega
r_{c}^{3}}\;\left(\frac{|\Lambda|}{\pi|k_{r}|}\right)^{1/2}\;\;\;{\rm
if}\;\;\;k_{r}^{2}<0$ (45)

and

$\displaystyle p^{\prime\rm rad}_{r}(r_{c}^{+})=$
$\displaystyle\sum_{k=-\infty}^{k=+\infty}a_{k}{p^{\prime}}_{k}^{m}(r_{c}^{+},s)H_{k}^{m}(\theta,s)\;\;\;{\rm
with}$ (46) $\displaystyle{p^{\prime}}_{k}^{m}(r_{c}^{+},s)=$
$\displaystyle\frac{\rho_{0}^{1/2}}{r_{c}}\;\left(\frac{N_{0}^{2}-\omega^{2}}{\pi
k_{r}}\right)^{1/2}\cos(\psi)$ $\displaystyle\;\;\;{\rm if}\;\;\;k_{r}^{2}>0$
(47) $\displaystyle{p^{\prime}}_{k}^{m}(r_{c}^{+},s)=$
$\displaystyle\frac{\rho_{0}^{1/2}}{r_{c}}\;\left(\frac{N_{0}^{2}-\omega^{2}}{\pi|k_{r}|}\right)^{1/2}$
$\displaystyle\;\;\;{\rm if}\;\;\;k_{r}^{2}<0$ (48)

where $N_{0}=N(r=r_{c}^{+})$ and
$\psi_{k}^{m}(s)=\int_{r_{c}}^{r_{b}}k_{r}\;dr-\pi/4\approx\frac{\pi^{2}s\sqrt{\Lambda}}{\Omega\Pi_{0}}-\pi/4$
because $\Lambda c_{0}^{2}/(r^{2}\omega^{2})\gg 1$ away from $r_{b}$ and
$\omega\ll N$.

The expressions of ${\xi}^{\rm conv}_{r}(r_{c}^{-})$ and $p^{\prime\rm
conv}(r_{c}^{-})$ are

$\displaystyle{\xi}^{\rm conv}_{r}(r_{c}^{-})=$
$\displaystyle\frac{r_{c}}{s}\sum_{\ell_{i}=1}^{+\infty}b_{\ell_{i}}\gamma_{\ell_{i}}^{m}(s)\tilde{P}_{\ell_{i}}^{m}(\cos\theta)$
(49) $\displaystyle p^{\prime\rm conv}(r_{c}^{-})=$
$\displaystyle\rho_{0}r_{c}^{2}\omega^{2}\sum_{\ell_{i}=1}^{+\infty}b_{\ell_{i}}P_{\ell_{i}}^{m}(1/s)\tilde{P}_{\ell_{i}}^{m}(\cos\theta)$
(50)

where, for convenience, we used $\tilde{P}_{\ell_{i}}^{m}$, the normalised
form of the associated Legendre polynomial, i.e.
$\int_{0}^{\pi}\;\left(\tilde{P}_{\ell_{i}}^{m}(\cos\theta)\right)^{2}\sin\theta\;d\theta=1$.

The continuity conditions then read:

$\displaystyle\frac{r_{c}}{s}\sum_{\ell_{i}=1}^{+\infty}b_{\ell_{i}}\gamma_{\ell_{i}}^{m}(s)\tilde{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)$
(51)
$\displaystyle\rho_{0}r_{c}^{2}\omega^{2}\sum_{\ell_{i}=1}^{+\infty}b_{\ell_{i}}P_{\ell_{i}}^{m}(1/s)\tilde{P}_{\ell_{i}}^{m}(\cos\theta)=$
$\displaystyle\sum_{k=-\infty}^{k=+\infty}a_{k}{p^{\prime}}_{k}^{m}(r_{c}^{+},s)H_{k}^{m}(\theta,s)$
(52)

Multiplying both equations by $H_{k}^{m}(\theta,s)$ and taking the integral
$\int_{0}^{\pi}...\sin\theta\;d\theta$ leads to :

$\displaystyle\frac{r_{c}}{s}\sum_{\ell_{i}=1}^{+\infty}b_{\ell_{i}}\gamma_{\ell_{i}}^{m}(s)c^{m}_{k,\ell_{i}}=$
$\displaystyle a_{k}\zeta_{k}^{m}(r_{c}^{+},s)$ (53)
$\displaystyle\rho_{0}r_{c}^{2}\omega^{2}\sum_{\ell_{i}=1}^{+\infty}b_{\ell_{i}}P_{\ell_{i}}^{m}(1/s)c^{m}_{k,\ell_{i}}=$
$\displaystyle a_{k}{p^{\prime}}_{k}^{m}(r_{c}^{+},s)$ (54)

for each integer $k$, where
$c^{m}_{k,\ell_{i}}(s)=\int_{0}^{\pi}\tilde{P}_{\ell_{i}}^{m}(\cos\theta)H_{k}^{m}(\theta,s)\sin\theta\;d\theta$
is a coupling coefficient between $\tilde{P}_{\ell_{i}}^{m}$ and
$H_{k}^{m}(\theta,s)$. Eliminating $a_{k}$ we obtain:

$\sum_{\ell_{i}=1}^{+\infty}\;\left(\frac{r_{c}}{s}\gamma_{\ell_{i}}^{m}(s){p^{\prime}}_{k}^{m}(r_{c}^{+},s)-\rho_{0}r_{c}^{2}\omega^{2}P_{\ell_{i}}^{m}(1/s)\zeta_{k}^{m}(r_{c}^{+},s)\right)c^{m}_{k,\ell_{i}}(s)b_{\ell_{i}}=0$
(55)

for each integer $k$. Introducing the expressions of
${p^{\prime}}_{k}^{m}(r_{c}^{+},s)$ and $\zeta_{k}^{m}(r_{c}^{+},s)$, we
obtain for $k\geq 0$:

$\displaystyle\sum_{\ell_{i}=1}^{+\infty}\;\left(\gamma_{\ell_{i}}^{m}(s)\frac{N_{0}}{2\Omega}\cos(\psi_{k}^{m}(s))+P_{\ell_{i}}^{m}(1/s)\Lambda_{k}^{m}(s)^{1/2}\sin(\psi_{k}^{m}(s))\right)c^{m}_{k,\ell_{i}}(s)b_{\ell_{i}}$
$\displaystyle=0$ (56)

and for $k<0$:

$\displaystyle\sum_{\ell_{i}=1}^{+\infty}\;\left(\gamma_{\ell_{i}}^{m}(s)\frac{N_{0}}{2\Omega}-P_{\ell_{i}}^{m}(1/s)|\Lambda_{k}^{m}(s)|^{1/2}\right)c^{m}_{k,\ell_{i}}(s)b_{\ell_{i}}=0$
(57)

where we used the approximation $(1-\omega^{2}/N_{0}^{2})^{1/2}\approx 1$.

To get approximate analytical solutions, we may try to truncate this infinite
system of equations to a finite system. We first consider the case of the
resonance between the Kelvin modes $(\ell=1,m=-1)$ and the $(\ell_{i}=3,m=-1)$
inertial mode near $s_{*}=11.32$. To analyse the geometrical matching between
these modes at $r=r_{c}$, the upper panel of Fig. 39 shows how
$H_{0}^{-1}(\theta,s_{*})$, projects onto Legendre polynomials of same $m$ and
equatorial parity. We observe that the two main contributions come from
$P_{1}^{-1}$ and $P_{3}^{-1}$. Conversely, the lower panel of Fig. 39 shows
that the inertial mode $P_{3}^{-1}$ mostly projects onto $H_{0}^{-1}(s_{*})$
and $H_{-2}^{-1}(s_{*})$. This suggests to truncate the infinite system to
linear combinations of $P_{3}^{-1}$ and $P_{1}^{-1}$ on the convective side
and of $H_{0}^{-1}(s)$ and $H_{-2}^{-1}(s)$ on the radiative side. For
consistency one should also consider the modes that are excited by
$H_{-2}^{-1}(s_{*})$ and $P_{1}^{-1}$ at the interface. The projections of
$H_{-2}^{-1}(s)$ and $P_{1}^{-1}$ onto Legendre polynomials and Hough
functions respectively are displayed in Fig. 40. It shows that, except for the
contribution of $P_{5}^{-1}$ to $H_{-2}^{-1}(s)$, the two functions
$H_{-2}^{-1}(s)$ and $P_{1}^{-1}$ are dominantly projected onto respectively
$P_{3}^{-1}$ and $P_{1}^{-1}$, and $H_{0}^{-1}(s_{*})$ and $H_{-2}^{-1}(s)$.
We thus expect that the truncation of the infinite system to $P_{1}^{-1}$,
$P_{3}^{-1}$, $H_{0}^{-1}(s)$, and $H_{-2}^{-1}(s)$ provides an approximate
representation of the mixed modes near $s_{*}=11.32$.

Figure 39: Upper panel: Projection of the Hough function
$H_{k=0}^{m=-1}(s_{*}=11.32)$ on the associated Legendre polynomial of the
same $m$ and parity. Lower panel: Projection of the associated Legendre
polynomial $P_{\ell_{i}=3}^{m=-1}$ on the Hough functions of the same $m$ and
parity determined at $s_{*}=11.32$.

Figure 40: Upper panel: Projection of the Hough function
$H_{k=-2}^{m=-1}(s_{*}=11.32)$ on the associated Legendre polynomial of the
same $m$ and parity. Lower panel: Projection of the associated Legendre
polynomial $P_{\ell_{i}=1}^{m=-1}$ on the Hough functions of the same $m$ and
parity determined at $s_{*}=11.32$.

With this assumption, the system of Eqs. D.3 and 57 reduces to the two
equations corresponding to $k=0$ and $k=-2$:

$\displaystyle\;\left(\gamma_{1}\frac{N_{0}}{2\Omega}\cos(\psi_{0})+P_{1}\Lambda_{0}^{1/2}\sin(\psi_{0})\right)c_{0,1}b_{1}+$
$\displaystyle\;\left(\gamma_{3}\frac{N_{0}}{2\Omega}\cos(\psi_{0})+P_{3}\Lambda_{0}^{1/2}\sin(\psi_{0})\right)c_{0,3}b_{3}$
$\displaystyle=0$ (58)
$\displaystyle\;\left(\gamma_{1}\frac{N_{0}}{2\Omega}-P_{1}|\Lambda_{-2}^{1/2}|\right)c_{-2,1}b_{1}+$
$\displaystyle\;\left(\gamma_{3}\frac{N_{0}}{2\Omega}-P_{3}|\Lambda_{-2}^{1/2}|\right)c_{-2,3}b_{3}$
$\displaystyle=0$ (59)

where the $m=-1$ indices and the $s$ dependence have been omitted for clarity.

Having a non-trivial solution for this system requires that its determinant
vanishes. This can be written as:

${\cal A}\cos(\psi_{0})+{\cal B}\sin(\psi_{0})=0$ (60)

with

$\displaystyle{\cal A}=$ $\displaystyle
c_{0,1}c_{-2,3}\gamma_{1}P_{3}\Lambda_{-2}^{1/2}\frac{2\Omega}{N_{0}}$
$\displaystyle+\left[(c_{0,3}c_{-2,1}-c_{0,1}c_{-2,3})\gamma_{1}-c_{0,3}c_{-2,1}P_{1}\Lambda_{-2}^{1/2}\frac{2\Omega}{N_{0}}\right]\gamma_{3}$
(61) $\displaystyle{\cal B}=$
$\displaystyle\Biggl{[}c_{0,3}c_{-2,1}\gamma_{1}P_{3}-(c_{0,3}c_{-2,1}-c_{0,1}c_{-2,3})\Lambda_{-2}^{1/2}P_{1}P_{3}\frac{2\Omega}{N_{0}}\Biggr{.}$
$\displaystyle-\Biggl{.}c_{0,1}c_{-2,3}\gamma_{3}P_{1}\Biggr{]}\Lambda_{0}^{1/2}\frac{2\Omega}{N_{0}}$
(62)

We can simplify this equation using the fact that $\Omega/N_{0}\ll 1$ in
$\gamma$ Dor stars. However, one should pay attention that $\gamma_{3}$ is
also small as we consider spin parameters close to $s_{*}$ and
$\gamma_{3}(s_{*})=0$. We obtain:

$\displaystyle{\cal A}\approx$ $\displaystyle
c_{0,1}c_{-2,3}\gamma_{1}P_{3}\Lambda_{-2}^{1/2}\frac{2\Omega}{N_{0}}+(c_{0,3}c_{-2,1}-c_{0,1}c_{-2,3})\gamma_{1}\gamma_{3}$
(63) $\displaystyle{\cal B}\approx$
$\displaystyle\;\left(c_{0,3}c_{-2,1}\gamma_{1}P_{3}-c_{0,1}c_{-2,3}\gamma_{3}P_{1}\right)\Lambda_{0}^{1/2}\frac{2\Omega}{N_{0}}$
(64)

Close to $s_{*}$ we can write $\gamma_{3}\approx\gamma^{\prime}_{3}(s-s_{*})$,
while the other terms in ${\cal A}$ and ${\cal B}$ remain approximately equal
to their value at $s_{*}$. Keeping the dominant terms, we get:

$\displaystyle{\cal A}\approx$ $\displaystyle
c_{0,1}c_{-2,3}\gamma_{1}P_{3}\Lambda_{-2}^{1/2}\frac{2\Omega}{N_{0}}+(c_{0,3}c_{-2,1}-c_{0,1}c_{-2,3})\gamma_{1}\gamma^{\prime}_{3}(s-s_{*})$
(65) $\displaystyle{\cal B}\approx$ $\displaystyle
c_{0,3}c_{-2,1}\gamma_{1}P_{3}\Lambda_{0}^{1/2}\frac{2\Omega}{N_{0}}$ (66)

so that Eq. 60 becomes:

$\cot(\psi_{0}(s))=-\frac{\epsilon/V}{s-s_{c}}$ (67)

where

$\displaystyle\psi_{0}=$
$\displaystyle\frac{\pi^{2}s\sqrt{\Lambda_{0}}}{\Omega\Pi_{0}}-\frac{\pi}{4}$
(68) $\displaystyle\epsilon=$ $\displaystyle\Omega/N_{0}$ (69) $\displaystyle
s_{c}=$ $\displaystyle
s_{*}+\;\left(\frac{1}{1+f(s_{*})}\frac{2P_{3}(1/s_{*})}{\gamma^{\prime}_{3}(s_{*})}\sqrt{\Lambda_{-2}(s_{*})}\right)\epsilon$
(70) $\displaystyle\frac{1}{V}=$
$\displaystyle\frac{f(s_{*})}{1+f(s_{*})}\frac{2P_{3}(1/s_{*})}{\gamma^{\prime}_{3}(s_{*})}\sqrt{\Lambda_{0}(s_{*})}$
(71) $\displaystyle f(s_{*})=$
$\displaystyle-\frac{c_{0,3}(s_{*})c_{-2,1}(s_{*})}{c_{0,1}(s_{*})c_{-2,3}(s_{*})}$
(72)

Following Tokuno & Takata (2022), the relation Eq. 67 translates into a
Lorentzian-shaped dip in the $\Delta P=f(s)$, the dip being centred on $s_{c}$
and of width $\sigma=\epsilon/V$.

We obtain similar dip models for the $(\ell=1,m=0)$ vs $(\ell_{i}=3,m=0)$
resonance near $s_{*}=\sqrt{5}$ and the $(\ell=2,m=-2)$ vs $(\ell_{i}=4,m=-2)$
resonance near $s_{*}=8.624$. Indeed, as for the $(\ell=1,m=-1)$ vs
$(\ell_{i}=3,m=-1)$ resonance, the projection of $H_{k}^{m}$ on the Legendre
polynomials and the projection of $P_{\ell_{i}}^{m}$ on the Hough functions
suggest it is sufficient to add two modes to get approximate continuity
conditions. These are $H_{k-2}^{m}$, the Hough function associated with the
smallest exponential decrease, and $P_{\ell_{i}-2}^{m}$. When expressed for
these three dips, the width and location of the dips are given by Eqs. (16)
and (17).

To conclude, we note that up to the (infinite) system of equations Eq. 55
describing the continuity conditions, the present model is identical to that
of Tokuno & Takata (2022). From there, we improved the continuity conditions
by adding the contribution of the Hough function $H_{k-2}^{m}$ on the
radiative side and of the Legendre polynomial $P_{\ell_{i}-2}^{m}$ mode on the
convective side. While we argue this is enough to approximate the resonances
considered in this paper, we stress that this model cannot be generalised to
all dips. Indeed, for other resonances (e.g. the $(\ell=2,m=0)$ vs
$(\ell_{i}=4,m=0)$ resonance), the projections of the $H_{k}^{m}$ Hough
function on the Legendre polynomials and of the $P_{\ell_{i}}^{m}$ Legendre
polynomial on the Hough functions strongly suggest that many more free
oscillation modes will be involved leading to much more complex mixed modes
than the ones considered here.