© The Authors 2025

1. Introduction

The efficiency of angular momentum transport inside stars inferred from standard hydrodynamic evolution modelling is underestimated compared to observations. Cores are rotating slower than they are expected to during a significant portion of the stellar evolution, that is, on the main sequence (MS Aerts et al. 2017; Eggenberger et al. 2019a), as well as in subgiants (SG, Deheuvels et al. 2014; Gehan et al. 2018; Eggenberger et al. 2019b), red giant branch stars (RGB, Deheuvels et al. 2012; Mosser et al. 2012), and red clump stars (RC, Deheuvels et al. 2015; Mosser et al. 2024). As theoretical insights into the dynamical processes at play in stellar interiors are being developed, with candidates such as internal gravity waves (Talon & Charbonnel 2005; Belkacem et al. 2015; Pinçon et al. 2017), meridional circulation (Rieutord 2006; Decressin et al. 2009), shear-driven turbulence (Zahn 1992; Mathis et al. 2018), and internal magnetism (Spruit 1999, 2002; Mathis & Zahn 2005; Fuller et al. 2019; Petitdemange et al. 2023; Barrère et al. 2023), tools needed to accurately probe the rotation rate at different locations inside the star are also required in order to better identify transport mechanisms.

During the MS, gravity oscillation modes (g-modes) are unique probes of the properties of the stellar radiative zone. However, the rotation rate of the convective core remains unconstrained. Nevertheless, measuring convective-core rotation holds valuable insights into stellar dynamics and chemical mixing, as the properties of convection themselves are affected by rotation (Augustson & Mathis 2019), and rotation is a key parameter for core dynamo (Brun et al. 2005; Augustson et al. 2016), which is a progenitor of internal magnetism in evolved phases, and is key to modelling the properties of compact objects (e.g. Bagnulo & Landstreet 2022, for white dwarfs). Measuring the core rotation of intermediate-mass and massive stars would thus allow us to better understand angular-momentum redistribution and magnetic-field generation in convective zones happening over dynamical timescales (Brun et al. 2004, 2022).

Two families of stars on the MS have attracted significant interest, exhibiting clear g-mode pulsations in Kepler (Li et al. 2020; Pápics et al. 2017; Pedersen 2022) and TESS (Antoci et al. 2019; Garcia et al. 2022) data: the massive slowly pulsating B stars (SPB) and the intermediate-mass γ Doradus stars. In the case of these g-mode pulsators, rotation alters the character of the gravity modes propagating in the envelope. The Coriolis force acts as a restoring force alongside the buoyancy and modifies the structure of the modes, becoming gravito-inertial modes. Rotation changes the spatial structure of such modes (Lee & Saio 1997; Dintrans & Rieutord 2000; Mathis et al. 2008), as well as their spacing in period (Ballot et al. 2012; Bouabid et al. 2013), which would be constant in the asymptotic, non-rotating case (Tassoul 1980).

Following this observation, the period spacing between gravity modes of consecutive radial orders and of the same horizontal structure holds valuable information about the radiative interior of the stars. Any deviation from this constant pattern can arise from dynamical processes and mixing. Modulations in the period spacing were proven to arise from sharp discontinuities in the mean molecular weight (Miglio et al. 2008) and can be used as a test for profiles of convective penetration and of the structure of core-envelope boundary (Michielsen et al. 2019). The slope of the period spacings in an inertial frame has been used to measure near-core rotation rates, the location at which g-modes reach the highest sensitivity (Van Reeth et al. 2015; Ouazzani et al. 2017; Christophe et al. 2018). With the surface rotation rate inferred from near-surface acoustic modes (p-modes) in the case of mixed γ Doradus – δ Scuti pulsators (Kurtz et al. 2014; Saio et al. 2015), or inferred from rotational spot modulation (Van Reeth et al. 2018), the differential near-core to surface rotation was measured in γ Doradus stars, and was proven to be limited, with the ratio of surface rotation to core rotation ranging between 0.97 and 1.02 in this latter study. The sample of γ Doradus stars exhibiting gravito-inertial modes is now large, with 611 stars in Kepler data showing prograde Kelvin g-modes, which have the highest visibility, and retrograde r-modes, corresponding to global Rossby waves (see Townsend 2003; Mathis et al. 2008, for gravito-inertial modes typologies). Among those, 58 have measured near-core to surface differential rotation (Li et al. 2020).

Recently, Ouazzani et al. (2020) made a breakthrough in the analysis of the period-spacing pattern, proving that pure inertial modes propagating in the convective core of the star couple through the convective-radiative boundary with gravito-inertial modes. This process results in a characteristic dip in the period-spacing pattern, which was later confirmed by Saio et al. (2021), and observed in Kelvin g-modes series of 16 γ Doradus stars. The location of the dip in period, as well as its width and depth, were proven to depend on stellar parameters and evolution stage, and the prescription of convective penetration in the radiative zone. As a detailed understanding of the formation of the dip was lacking, Tokuno & Takata (2022) described the coupling mechanism theoretically and found a Lorentzian shape for the dip that is specific to this coupling mechanism compared to the dips created by a strong gradient in molecular weight (Kurtz et al. 2014; Saio et al. 2015; Schmid & Aerts 2016; Murphy et al. 2016; Pedersen et al. 2018; Michielsen et al. 2019; Li et al. 2019; Wu et al. 2020). Tokuno & Takata (2022) derives a key coupling parameter controlling the formation of the dip, ϵ, depending on the density stratification as well as the near-core rotation rate of the star. Galoy et al. (2024) further conducted a thorough numerical study of the inertial dip formation using the spectral 2D oscillation code TOP (Reese et al. 2009; Ballot et al. 2012), deriving empirical relations for the width and the location of the dips as a function of the frequencies of the pure inertial modes in an isolated convective core and the near-core stratification gradient. This opens the possibility of future studies allowing structure inversions from the inertial dip structure observed in rapid g-mode pulsators from asteroseismic data. Galoy et al. (2024) also derives an improvement of the model of Tokuno & Takata (2022), taking into account a multi-mode interaction from both sides of the boundary.

The ϵ coupling parameter was estimated by forward modelling methods by Aerts & Mathis (2023) for a sample of 37 γ Doradus and 26 SPB stars, and its dependency on the stellar age and the near-core rotation was established. The SPB sample currently used is not as large as the γ Doradus sample and contains slower rotators, which partly explains the unconfirmed detection of inertial dips in SPBs. However, the γ Doradus stars with confirmed inertial-dip detection do not stand out in terms of their coupling parameter compared to the overall sample. Arguments based on the coupling parameter alone do not appear to be able to fully explain the relatively low proportion of inertial dips findings in observed γ Doradus stars period-spacing patterns (3%). This is mitigated by the fact that no ensemble analysis of the inertial dip has been made so far, with Saio et al. (2021) relying on a visual inspection for their study.

While Ouazzani et al. (2020), Tokuno & Takata (2022), and Galoy et al. (2024) remained in the framework of solid-body rotation, Saio et al. (2021) explored the possibility of allowing a steep differential rotation profile from the convective core to the near-core regions. Their study reveals a minimum ratio of the core rotation rate over the near-core rotation rate of 0.85 for the stars exhibiting inertial dips. Moyano et al. (2024) further used these results to put constraints on the angular momentum redistribution processes at play in the radiative interior of γ Doradus stars. In this respect, it is important to further explore the effect of differential rotation, as we expect it to Doppler shift the location of the dip in the period-spacing pattern, and it might also lead to a modification in its morphology. We thus aim to extend the theoretical explanations of the signature of rotation on the g-modes period spacing given by Tokuno & Takata (2022) to allow for differential rotation, and investigate the effects of differentiality on the coupling mechanism, the subsequent dip formation, and the shape of the dip. We propose an analytical two-zone model, with both the convective core and the radiative envelope rotating as solid bodies with different rotation rates, serving as a first laboratory towards the understanding of the effect of differential rotation on the coupling. We aim to provide predictions for the measurement of such two-zone differential rotation, to investigate the detectability of such differentiality, and to compare its effect with other processes that shift the location of the dip: the density stratification in the core and the gradient of stratification at the interface of the two zones (Ouazzani et al. 2020).

The outline of the present paper is as follows: in Section 2, we describe the theoretical framework used, and summarise the results of Tokuno & Takata (2022). In Section 3, we derive the coupling equation in the presence of differential rotation and solve it both numerically and analytically. We describe the characteristics of the modified Lorentzian profile obtained with differential rotation in the case of Kelvin g-modes. In Section 4, we summarise the effects of differential rotation and discuss our findings in the context of other processes influencing the position of the dip, tackling the detectability of such differential rotation in the analysis of inertial dips. We then conclude in Section 5 on the potentiality of differential rotation measurements and also outline the limitations of our model and the leads that could be pursued by future studies on inertial dips.

2. Framework and derivation of the structure of the modes from both sides of the convective–radiative boundary

2.1. Selection of a γ Doradus star

We are working on a model of g-mode pulsating main sequence star displaying a convective core surrounded by a radiative envelope. Two different stellar pulsators are classically comprised in this analysis: γ Doradus stars, of 1.3 to 2.0 M_⊙, and SPB stars, of 2 to 7 M_⊙. In this work, we will focus on the case of γ Doradus stars, as (1) more oscillation modes are detected in this stellar class, (2) the structural discontinuities exhibited by SPB stars would make the period-spacings modulations due to prominent molecular weight or temperature gradient, making the detection of a dip more difficult, and (3) the sample of γ Doradus stars observed period spacings is more abundant than the one of SPB stars and inertial dips in SPB stars were searched for in Saio et al. (2021) and none were found in the sample studied by Pápics et al. (2017).

The classical instability strip of γ Doradus covers in the low-mass range stars close to their zero-age main sequence (ZAMS), and in the high-mass range stars close to their terminal-age main sequence (TAMS). However, this picture is complicated by the detection of γ Doradus at locations in the colour–magnitude diagram far from the classical instability strip (see e.g. Fig. 12 of Li et al. 2020).

2.2. Physical framework and rotational setup

We now present the relevant assumptions and physical approximations used to study core-to-envelope mode coupling, and discuss them with regard to the state of the art on gravito-inertial modes in the envelope and pure inertial modes in the core.

2.2.1. Gravito-inertial mode typology

Gravito-inertial modes form the predominant features in the spectrum of γ Doradus stars. Their restoring forces are both buoyancy and the Coriolis force and are strongly affected by the latter both in their spatial structure and frequencies, compared to their pure gravity mode counterparts in non-rotating stars. Gravito-inertial modes can be classically separated into four classes based on their geometry and their restoring mechanism, as described for instance in Mathis et al. (2008): Poincaré modes, existing in the non-rotating case, which are internal gravity waves modified by the Coriolis acceleration in rotating stars. Their radial wavenumber rapidly increases with rotation (Townsend 2003), as their damping; r-modes, which are purely retrograde waves existing only for rapid rotators in the sub-inertial regime, driven by the conservation of the vorticity and curvature effect; Yanai modes, mix of the aforementionned modes with smaller radial nodes, hence less damped than the Poincaré modes; and Kelvin g-modes, also driven by the conservation of vorticity, which are purely prograde modes trapped around the equator of the star.

We stay within the framework of Kelvin g-modes, for two reasons: First, Kelvin g-modes have the highest visibility in γ Doradus stars and were by far the most observed in the sample of 611 stars analysed in Li et al. (2020). Second, r-modes that were also observed in this sample would not display dips arising from the interacting mechanism we are describing, as Saio et al. (2021) pointed out. Poincaré modes remain unobserved, while only 7 stars show Yanai modes to this date (Van Reeth et al. 2018; Li et al. 2019).

2.2.2. Frequency regimes and relevant approximations

These gravito-inertial modes have different propagation regions based on the frequency regime at which they appear. As studied by Prat et al. (2016), noting σ_env and σ_core the frequencies of the wave in the frames co-rotating with the envelope and the core, respectively, gravito-inertial waves are evanescent in the core in the super-inertial regime (σ_env > 2Ω_env and σ_core > 2Ω_core), whereas they can penetrate in the core as pure inertial modes in the sub-inertial regime (σ_core < 2Ω_core and σ_env < 2Ω_env). We therefore remain within the sub-inertial regime in both regions in our study, ignoring a trans-inertial regime in which we have (σ_env < 2Ω_env and σ_core > 2Ω_core or σ_env > 2Ω_env and σ_core < 2Ω_core), which results in a chaotic behaviour of the rays and reduced coupling between the core and the envelope (Prat et al. 2018). Gravito-inertial modes are described in solid lines in Fig. 1. They are propagating between two turning points for which their frequency in the co-rotating frame equates the near-core Brunt-Väisälä (N) frequency (r_a) and the minimum between the near-surface Brunt-Väisälä and Lamb (L) frequencies (r_b).

In the sub-inertial regime, waves thus propagate in the core and gain a purely inertial character, because buoyancy can no longer act as a restoring force in an unstratified region. While these waves have been long studied and observed in geophysics, Wu (2005) made a comprehensive study of the properties of pure inertial modes in an astrophysical context.

The equation governing the structure of pure inertial modes is only separable for solid-body rotating spheres with constant density. We thus choose for the density profile a uniform averaged density in the core, for the study to remain analytical. This assumption will be further discussed in Section 4.6. Pure inertial waves are symbolized as dashed oscillations in Fig. 1 and propagate within the convective core r < R_core.

We adopt the traditional approximation of rotation (TAR) within the envelope (e.g. Gerkema et al. 2008). The TAR consists in neglecting the latitudinal component of the rotation vector in the Coriolis acceleration. The TAR is known to hold in highly stratified layers of stellar interiors. For the sub-inertial Kelvin g-modes in which we are interested, we must stay within the regime N > ( ≫ )2Ω, for which the TAR solution approaches the solution obtained with a full treatment of the Coriolis force (Prat et al. 2016). This cannot be satisfied in the convective core, with the Brunt-Väisälä frequency reaching non-positive real values in this zone. Hence the inertial dips could only be seen in numerical studies for which the TAR was lifted (Saio et al. 2018; Ouazzani et al. 2020). As for the envelope, the strong molecular weight gradient near the core ensures high values of the Brunt-Väisälä frequency even in the near-core region. This sharp stratification gradient at the core-envelope boundary ensures that the TAR holds for gravito-inertial modes, and is thus part of our studied cases, as well as in Tokuno & Takata (2022). This allows the region of radial coordinates [R_core; r_a] in which the TAR fails and gravito-inertial modes become evanescent to remain small so that the structure of such modes remain close to the form predicted by the TAR at r_a. Throughout the present study, we neglect as well the wave’s gravific potential fluctuation (Cowling 1941), approximation justified by the strongly radially oscillating character of the modes.

2.2.3. Differential rotation framework

We allow a differential rotation between the convective core and the radiative surrounding envelope. This takes the form of a bi-layer rotation rate, with Ω_env the rotation rate of the radiative envelope and Ω_core the rotation rate of the convective core (see Fig. 1), parametrized as

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\text{env}}={\alpha }_{\text{rot}}{\mathrm{\Omega }}_{\text{core}},\end{array}$$(1)

with α_rot accounting for the degree of differential rotation between the core and the envelope.

In the convective core, the rotation is known to approach cylindricity from the application of the Taylor-Proudmann theorem. The core rotation profile was treated by hydrodynamical (Browning et al. 2004) and magneto-hydrodynamical (Brun et al. 2005; Featherstone et al. 2009; Augustson et al. 2016) simulations, showing an impact of the magnetic field on the limitation of cylindrical differential rotation in the core (5% of radial differential rotation at the equator for model M4 of Augustson et al. 2016). We decide to stay in the framework of solid-body rotation to maintain the analyticity of the inertial mode solutions. Indeed, the introduction of a convective core differential rotation breaks the regularity of the solutions: the wave structure is of an attractor type (Baruteau & Rieutord 2013) and critical layers can form that damp the inertial waves (Astoul et al. 2021). Additionally, seismology allows us to probe the average rotation rate of a spherical zone, thus advocating for a solid-body treatment of the core rotation rate in our study, while a two-zone analytical treatment is useful to get a precise physical understanding of the dip formation in the context of differential rotation.

As for the radiative zone, transport processes can lead to differential rotation (Maeder & Meynet 2000; Rogers 2015). The topic was treated in Van Reeth et al. (2018), where it was shown that differential rotation within the envelope was hard to constrain from period spacings, unless multiple mode series are exhibited, and was less than 5% between the surface and the near-core region in the stars in which a measurement was possible from starspot modulations. This approximation of a solid-body rotating radiative zone is thus done because differential rotation is limited and the key quantity for the coupling is only the near-core rotation. We hereby emphasize that the uniform rotation rate of the envelope Ω_env we take would correspond to the near-core rotation rate in a differentially rotating radiative zone. Differential rotation in our study is then probed between the convective core and the inner radiative envelope near the core and does not correspond to the differential rotation probed by the curvature gravito-inertial modes and rotational spot modulations (see e.g. Van Reeth et al. 2018).

2.3. Formalism and description of the relevant reference frames

We here describe the expansion of the quantities in both the two co-rotating frames relevant to the present study and the inertial one. In the rotational set-up previously described, each scalar and vectorial field will be expanded in Fourier’s series in the azimuthal angle φ and time t:

$$\begin{array}{c}\hfill X(r,\theta ,\phi ,t)={\displaystyle \sum _{{\sigma }_{\text{zone}},m}}X\prime (r,\theta )e^{i(m\phi +{\sigma }_{\text{zone}}t)},\end{array}$$(2)

$$\begin{array}{c}\hfill \mathit{x}(r,\theta ,\phi ,t)={\displaystyle \sum _{{\sigma }_{\text{zone}},m}}\mathit{x}\prime (r,\theta )e^{i(m\phi +{\sigma }_{\text{zone}}t)}.\end{array}$$(3)

We emphasise that with the chosen sign convention, modes with negative m are prograde, and modes with positive m are retrograde.

σ_zone being either σ_env or σ_core, the two local wave angular frequencies related to the frames rotating at respectively Ω_env and Ω_core, linked to the frequency in the inertial frame σ_in as

$$\begin{array}{c}\hfill {\sigma }_{\text{env}}={\sigma }_{\text{in}}+m{\mathrm{\Omega }}_{\text{env}},\end{array}$$(4)

$$\begin{array}{c}\hfill {\sigma }_{\text{core}}={\sigma }_{\text{in}}+m{\mathrm{\Omega }}_{\text{core}}.\end{array}$$(5)

The frequency decomposition in the two co-rotating frames is thus dependent on the zone in which the waves propagate. This is a differing key point from the solid-body rotation case treated by Tokuno & Takata (2022). Since we are considering a differentially rotating case, we derive the expressions for the quantities in the co-rotating frames, and further ensure the matching between the two zones in the inertial frame. Taking into account a Doppler shift related to the two zones on either side of the boundary in the phase term, the relevant quantities to match in the decomposition are therefore

$$\begin{array}{c}\hfill X\prime (r,\theta )e^{i(m\phi +{\sigma }_{\text{zone}}t)}{e}^{-im{\mathrm{\Omega }}_{\text{zone}}t}=X\prime (r,\theta )e^{i(m\phi +{\sigma }_{\text{in}}t)},\end{array}$$(6)

and

$$\begin{array}{c}\hfill \mathit{x}\prime (r,\theta )e^{i(m\phi +{\sigma }_{\text{zone}}t)}{e}^{-im{\mathrm{\Omega }}_{\text{zone}}t}=\mathit{x}\prime (r,\theta )e^{i(m\phi +{\sigma }_{\text{in}}t)}.\end{array}$$(7)

We further define the two spin parameters, which are the relevant parameters in the following analysis, as

$$\begin{array}{c}\hfill s_{\text{env}}=\frac{2{\mathrm{\Omega }}_{\text{env}}}{{\sigma }_{\text{env}}},\end{array}$$(8)

$$\begin{array}{c}\hfill s_{\text{core}}=\frac{2{\mathrm{\Omega }}_{\text{core}}}{{\sigma }_{\text{core}}}.\end{array}$$(9)

We place ourselves in the sub-inertial regime, for which σ_env < 2Ω_env and σ_core < 2Ω_core, thus s_env > 1 and s_env > 1. As argued in Sect. 2.2.2, this is the regime where inertial modes can propagate in the convective core, and we are thus considering low-frequency waves.

2.4. Expression of mode structures from both sides of the boundary

With this bi-layer differentially rotating framework defined, the derivation of the structure of modes from both sides of the boundary is very similar to Sections 3.1.1 and 3.1.2 of Tokuno & Takata (2022), respectively for gravito-inertial modes propagating in the radiative envelope and pure inertial modes propagating in the convective core. The spin parameter the functions depend on in those cases is no longer fixed at a single value s, but rather at s_env and s_core. The full calculations are given respectively in Appendices A and B for readability.

We obtain the expressions for the radial displacement and the Eulerian perturbation of pressure at R_core.

For gravito-inertial modes, corresponding to Hough solutions (see e.g. Lee & Saio 1997), the structure at the boundary depends on the hypothesis made on the continuity or discontinuity of the near-core Brunt-Väisälä profile. In the former continuous case,

$$\begin{array}{c}\hfill \frac{{\xi }_r\prime }{r}{|}_{R_{\text{core}}}={\displaystyle \sum _k}{a}_kϵX_k^m(s_{\text{env}}){\mathrm{\Theta }}_k^m(\mu ;s_{\text{env}}),\end{array}$$(10)

and

$$\begin{array}{c}\hfill \frac{p\prime }{\overline{p}{\mathrm{\Gamma }}_1}{|}_{R_{\text{core}}}=\frac{{\sigma }_{\text{env}}^2{R}_{\text{core}}^2}{c_S^2}{\displaystyle \sum _k}{a}_k{Y}_k^m(s_{\text{env}}){\mathrm{\Theta }}_k^m(\mu ;s_{\text{env}}).\end{array}$$(11)

In the latter discontinuous case, the functions $X_k^m$ and $Y_k^m$ are replaced by respectively ${\stackrel{\sim }{X}}_k^m$ and ${\stackrel{\sim }{Y}}_k^m$ described in Appendix A.

We have introduced μ = cos θ, ${\mathrm{\Theta }}_k^m$, the Hough functions described in the Appendix A and c_S the local sound speed. For pure inertial modes, we obtain Bryan solutions (Bryan 1889):

$$\begin{array}{c}\hfill \frac{{\xi }_r\prime }{r}{|}_{R_{\text{core}}}={\displaystyle \sum _l}{b}_l{C}_l^m(1/s_{\text{core}}){\stackrel{\sim }{P}}_l^m(\mu ),\end{array}$$(12)

and

$$\begin{array}{c}\hfill \frac{p\prime }{\overline{p}{\mathrm{\Gamma }}_1}{|}_{R_{\text{core}}}=\frac{{\sigma }_{\text{core}}^2{R}_{\text{core}}^2}{c_S^2}{\displaystyle \sum _l}{b}_l{P}_l^m(1/s_{\text{core}}){\stackrel{\sim }{P}}_l^m(\mu ),\end{array}$$(13)

where $P_l^m$ are the Legendre polynomials and ${\stackrel{\sim }{P}}_l^m$ their normalized counterparts (see Appendix B). These solutions take into account the Doppler shift induced by the differential rotation.

3. Coupling of the core and envelope modes in a differentially rotating context

With the expressions previously derived at both sides of the convective core boundary, we investigate the coupling of pure inertial modes in the convective core and gravito-inertial modes in the radiative envelope considering differential two-zone rotation. We describe the relevant assumptions allowing us to simplify the coupling equation. We solve the equation numerically and derive an analytical shape for the resulting inertial dip in the period-spacing pattern.

3.1. Matrix formulation of the problem

In this section, we investigate the problem of the continuity of the relevant quantities in the specific context of a continuous Brunt-Väisälä profile, and a two-zone differential rotation described in Fig. 1. The matching must be achieved at fixed azimuthal number m and frequency σ_in in the inertial frame. Thus, compared to a one-zone rotation model, the problem can no longer be equivalent to working at a fixed spin parameter in the two-zone case. The discontinuous case is treated in Section 3.5.

The two quantities that must stay continuous at the two sides of the interface are the radial Lagrangian displacement ξ_r and the Lagrangian perturbation of pressure $\delta p=p\prime +{\xi }_r(\frac{\partial \overline{p}}{\partial r})$. Except for the rotation frequency, background quantities are continuous at the interface in the case of the continuous Brunt-Väisälä profile.

Given the density concentration in the core of γ Doradus stars, and the critical rotation frequency scaling as ${\mathrm{\Omega }}_{\mathrm{crit}}\propto \sqrt{\mathcal{G}\rho }$, even for fast rotators close to Ω_surface = 0.5Ω_{crit, surface}, the rotation frequency at the edge of the convective core will be negligible compared to the critical rotation rate at this location. Thus the contribution of the centrifugal acceleration is considered to be negligible with respect to the local gravity at the core/envelope boundary in the present study (e.g. Ballot et al. 2010).

The continuity of the radial Lagrangian displacement being already ensured, the continuity of the Lagrangian perturbation of pressure is equivalent to the one of the Eulerian pressure p′. Using Eqs. (10) and (12), the matching equations for the Lagrangian displacement writes:

$$\begin{array}{cc}\hfill {\displaystyle \sum _k}& a_kϵX_k^m(s_{\text{env}}){\mathrm{\Theta }}_k^m(\mu ;s_{\text{env}})\hfill \\ \hfill & ={\displaystyle \sum _l}{b}_l{C}_l^m(1/s_{\text{core}}){\stackrel{\sim }{P}}_l^m(\mu ).\hfill \end{array}$$(14)

The matching of the Eulerian perturbation of the pressure writes, using Eqs. (11) and (13):

$$\begin{array}{cc}\hfill {\sigma }_{\text{env}}^2{\displaystyle \sum _k}& a_k{Y}_k^m(s_{\text{env}}){\mathrm{\Theta }}_k^m(\mu ;s_{\text{env}})\hfill \\ \hfill & ={\sigma }_{\text{core}}^2{\displaystyle \sum _l}{b}_l{P}_l^m(1/s_{\text{core}}){\stackrel{\sim }{P}}_l^m(\mu ).\hfill \end{array}$$(15)

Projecting each of the last two equations (⋅) on Hough functions, by taking ${\displaystyle {\int }_{-1}^1}(\mathbf{·}){\mathrm{\Theta }}_k^m(\mu ;s_{\text{env}})\mathrm{d}\mu$, we reduce the dimensionality of the equations, taking profit of the orthogonality of Hough functions.

$$\begin{array}{c}\hfill a_kϵX_k^m(s_{\mathrm{env}})={\displaystyle \sum _l}{b}_l{C}_l^m(1/s_{\text{core}})c_{k,l},\end{array}$$(16)

and

$$\begin{array}{c}\hfill {\sigma }_{\text{env}}^2{a}_k{Y}_k^m(s_{\text{env}})={\sigma }_{\text{core}}^2{\displaystyle \sum _l}{b}_l{P}_l^m(1/s_{\text{core}})c_{k,l},\end{array}$$(17)

with the coefficient c_k, l being defined as

$$\begin{array}{c}\hfill c_{k,l}\equiv {\displaystyle {\int }_{-1}^1}{\mathrm{\Theta }}_k^m(\mu ,s_{\text{env}}){\stackrel{\sim }{P}}_l^m(\mu )\mathrm{d}\mu .\end{array}$$(18)

Injecting the expressions for the a_k given in Eq. (16) in Eq. (17), one can recast the equations in a matrix form:

$$\begin{array}{c}\hfill [\mathcal{M}(s_{\text{core}},s_{\text{env}})-ϵ\mathcal{N}(s_{\text{core}},s_{\text{env}})]\mathit{b}=\mathbf{0},\end{array}$$(19)

b being the column vector of the terms b_l. The matrices are defined as

$$\begin{array}{c}\hfill {[\mathcal{M}]}_{k,l}=c_{k,l}{\sigma }_{\text{env}}^2{Y}_k^m(s_{\text{env}})C_l^m(1/s_{\text{core}}),\end{array}$$(20)

and

$$\begin{array}{c}\hfill {[\mathcal{N}]}_{k,l}=c_{k,l}{\sigma }_{\text{core}}^2{X}_k^m(s_{\text{env}})P_l^m(1/s_{\text{core}}).\end{array}$$(21)

For b to be non-trivial, we have the following condition:

$$\begin{array}{c}\hfill \text{det}[\mathcal{M}(s_{\text{core}},s_{\text{env}})-ϵ\mathcal{N}(s_{\text{core}},s_{\text{env}})]=0.\end{array}$$(22)

To determine the frequencies of the modes, the system has to be closed by an equation linking s_core and s_env. Such an equation arises from the fact that the match must occur at fixed σ_in, the frequency in the inertial frame. We have in each zone

$$\begin{array}{c}\hfill s_{\text{zone}}=\frac{2{\mathrm{\Omega }}_{\text{zone}}}{{\sigma }_{\text{in}}+m{\mathrm{\Omega }}_{\text{zone}}}.\end{array}$$(23)

Thus, with Eq. (1):

$$\begin{array}{c}\hfill s_{\text{env}}=\frac{{\alpha }_{\text{rot}}{s}_{\text{core}}}{1+\frac{m}{2}({\alpha }_{\text{rot}}-1)s_{\text{core}}}.\end{array}$$(24)

For concision, we define G as

$$\begin{array}{c}\hfill G(s)=\frac{{\alpha }_{\text{rot}}s}{1+\frac{m}{2}({\alpha }_{\text{rot}}-1)s},\end{array}$$(25)

so that s_env = G(s_core)⇔s_core = G⁻¹(s_env).

Due to the choice of the bi-layer rotation profile, the function G⁻¹ contains a pole at $1+\frac{m}{2}(1/{\alpha }_{\mathrm{rot}}-1)s_{\mathrm{env}}$, which corresponds to the case of infinitely low frequencies of the core modes: σ_core = 0. The dip model using this function G⁻¹ thus must be used for envelope spin parameters below this limit.

3.2. Synthesized coupling equation

First, putting ϵ = 0, case of $\frac{d{N}^2}{\mathit{dr}}\to \infty$ at R_core, and no avoided crossing, we see that we have either $Y_k^m=0$ (pure inertial mode in the core) or $C_l^m=0$ (pure gravito-inertial mode in the envelope). We then suppose $Y_k^m(s_{\text{env}})=O(ϵ)$ and $C_l^m(s_{\text{core}})=O(ϵ)$.

In a situation of ϵ ≠ 0 ≪ 1, the mode extends now over all radial coordinates of the star and is no longer confined either in the convective core or the radiative envelope. Considering a fixed energy input from the source of excitation (flux-blocking at the bottom of the outermost convective zone in the case of γ Doradus stars), the amplitude of an interacting mode on the two sides of the boundary must be underdominant compared to non-interacting modes. This is why we can consider that both $Y_k^m$ and $C_l^m$ are of the order O(ϵ).

The spectrum of the pure inertial modes in frequency is sparser than the one of the gravito-inertial modes. We then consider the interaction of one mode of the core with several modes of the envelope with a single value of k. The determinant thus reads:

$$\begin{array}{cc}& \text{det}[\mathcal{M}(s_{\text{core}},s_{\text{env}})-ϵ\mathcal{N}(s_{\text{core}},s_{\text{env}})]\hfill \\ \hfill & =\text{det}\left(\begin{array}{ccccccc}& & ⋮& /& ⋮& & \\ & O(1)& ⋮& /& ⋮& O(1)& \\ & & ⋮& /& ⋮& & \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ /& /& /& /& /& /& /\\ \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ & & ⋮& /& ⋮& & \\ & O(1)& ⋮& /& ⋮& O(1)& \\ & & ⋮& /& ⋮& \end{array}\right)\hfill \\ \hfill & =Y_k^m(s_{\text{env}})C_l^m(1/s_{\text{core}})\hfill \\ \hfill & \times \text{det}\left(\begin{array}{ccccccc}& & ⋮& & ⋮& & \\ & O(1)& ⋮& O(1)& ⋮& O(1)& \\ & & & & \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ & O(1)& ⋮& O(1/ϵ)& ⋮& O(1)& \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ & & & & \\ & O(1)& ⋮& O(1)& ⋮& O(1)& \\ & & ⋮& & ⋮& \end{array}\right).\hfill \end{array}$$(26)

The matrix at the right-hand side of this equation has one dominant term, at the crossing of the lth line and the kth column, the coupled modes. Putting this determinant null, each term in its development with ϵ must be null. We have for the leading term of order O(1/ϵ):

$$\begin{array}{c}\hfill c_{k,l}-ϵ\frac{P_l^m(1/s_{\text{core}}){\sigma }_{\text{core}}^2{X}_k^m(s_{\text{env}})}{C_l^m(1/s_{\text{core}}){\sigma }_{\text{env}}^2{Y}_k^m(s_{\text{env}})}{c}_{k,l}=0.\end{array}$$(27)

We thus obtain

$$\begin{array}{c}\hfill \frac{C_l^m(1/s_{\text{core}}){\sigma }_{\text{env}}^2{Y}_k^m(s_{\text{env}})}{P_l^m(1/s_{\text{core}}){\sigma }_{\text{core}}^2{X}_k^m(s_{\text{env}})}\simeq ϵ.\end{array}$$(28)

This reasoning holds only if c_k, l is not negligibly small, of the order of ϵ. This geometric consideration was first used as a criterion for the interaction of modes by Ouazzani et al. (2020) and lies on top of the criterion on the small value of ϵ developed in Tokuno & Takata (2022). Galoy et al. (2024) go further in the development of the determinant by not only taking the most dominant term but by reducing this infinite system of equations to a finite one, truncating it to the mode interactions with appreciable c_k, l values. The complexity is increased by the fact that the number of considered interactions is not known a priori and can vary from one interaction to another. Therefore, using Tokuno & Takata (2022)’s model can be considered as less precise but more general. We let this promising lead of using the Galoy et al. (2024) model in the context of differential rotation for further studies. Defining:

$$\begin{array}{c}\hfill F_l^m(s_{\text{core}})\equiv -\frac{C_l^m(1/s_{\text{core}})}{P_l^m(1/s_{\text{core}})},\end{array}$$(29)

the coupling equation between the Lagrangian displacement and the Lagrangian pressure perturbation in case of differential rotation thus reads, with the expressions of $X_k^m$ and $Y_k^m$ given in Appendix A:

$$\begin{array}{c}\hfill F_l^m(s_{\text{core}})\frac{\sqrt{3}}{2}\frac{\alpha }{{\mathrm{\Lambda }}_k^m(s_{\text{env}})}{s}_{\text{env}}^{2/3}\frac{{\sigma }_{\text{env}}^2}{{\sigma }_{\text{core}}^2}\\ \hfill & \times [cot(\frac{{\pi }^2{s}_{\text{env}}}{{\mathrm{\Omega }}_{\text{env}}{\mathrm{\Pi }}_0}-\frac{\pi }{6})+\frac{1}{\sqrt{3}}]\simeq ϵ.\hfill \end{array}$$(30)

3.3. Solving the coupling system

The solving method differs from the solid-body rotation case. Instead of finding the zeros of a 1D function, we have to find the zeros of a 2D function: g_ϵ(s_core, s_env) = 0, g_ϵ being defined as the LHS of Eq. (30) minus the coupling parameter. We compute numerically such zeros, for various ϵ values, and their location in the (s_core, s_env) plane is given by the blue lines in Fig. 2 for a specific ϵ value. Far from the spin parameter of the pure inertial mode $s_{\text{core}}^{\ast }$ (vertical dashed line in Fig. 2), the gravito-inertial modes in the envelope are not influenced by the interaction and the zeros are regularly spaced in envelope spin parameters (y-axis). On the other hand, the interaction with the pure inertial mode shifts the envelope spin parameters of the zeros close to $s_{\text{core}}^{\ast }$.

In addition, we are evolving on specific coordinates of the surface g_ϵ(s_core, s_env) = 0. Indeed, s_core and s_env are linked via Eq. (24). The solution of the coupling rotation is thus to be found on a characteristic unique of each α_rot. For differential rotations ranging from α_rot = Ω_env/Ω_core = 0.80 to α_rot = 1.05, the characteristics are marked as full lines on Fig. 2. The solid-body rotating case is shown by the red line. Depending on the characteristic, hence the differential rotation value, we can see that the location of the spin parameter of the envelope at which the interaction occurs is shifted. The local slope of the curve at $s_{\text{core}}=s_{\text{core}}^{\ast }$ also varies, which shows that with an increasing differential rotation from the core to the envelope, less gravito-inertial modes will be influenced by the pure inertial mode, and thus the dip will be thinner with increasing differentiality. The curvature of the characteristics introduces a non-symmetry of the dip relative to its minimum, which however stays small for this range of differential rotation. We further interpret this effect in Section 4.2.

3.4. Frequency-dependence of the dips

We recall from Section 3.1 that s_env = G(s_core). The coupling equation is then:

$$\begin{array}{c}{F}_l^m(G^{-1}(s_{\text{env}}))\frac{\sqrt{3}}{2}\frac{\alpha }{{\mathrm{\Lambda }}_k^m(s_{\text{env}})}{s}_{\text{env}}^{2/3}{\alpha }_{\text{rot}}^2\frac{G^{-1}{(s_{\text{env}})}^2}{s_{\text{env}}^2}\hfill \\ \hfill & \times [cot(\frac{{\pi }^2{s}_{\text{env}}}{{\mathrm{\Omega }}_{\text{env}}{\mathrm{\Pi }}_{0,\text{env}}}-\frac{\pi }{6})+\frac{1}{\sqrt{3}}]\simeq ϵ.\hfill \end{array}$$(31)

We now expand the front function near its zero, $s_{\text{core}}^{\ast }$, considering the interaction to happen at close spin parameter values. The values of the zeros are the ones given by Eq. (B.7). Those values are presented in Table 1 of Ouazzani et al. (2020). We obtain:

$$\begin{array}{cc}& \frac{F_l^m(s_{\text{core}})G{(s_{\text{core}})}^{2/3}}{{\mathrm{\Lambda }}_k^m(G(s_{\text{core}}))}{\left(\frac{s_{\text{core}}}{G(s_{\text{core}})}\right)}^2\hfill \\ \hfill & =(s_{\text{core}}-s_{\text{core}}^{\ast })\frac{\mathrm{d}{F}_l^m(s)}{\mathrm{d}s}{|}_{s_{\text{core}}^{\ast }}\frac{G{(s_{\text{core}}^{\ast })}^{2/3}}{{\mathrm{\Lambda }}_k^m(G(s_{\text{core}}^{\ast }))}{\left(\frac{s_{\text{core}}^{\ast }}{G(s_{\text{core}}^{\ast })}\right)}^2.\hfill \end{array}$$(32)

As in Tokuno & Takata (2022), we rewrite the coupling equation using V_diff:

$$\begin{array}{c}\hfill V_{\text{diff}}=-{\left[\frac{\mathrm{d}{F}_l^m(s)}{\mathrm{d}s}\frac{\sqrt{3}}{2}\frac{\alpha G{(s)}^{2/3}}{{\mathrm{\Lambda }}_k^m(G(s))}{\alpha }_{\text{rot}}^2{\left(\frac{s}{G(s)}\right)}^2\right]}_{s_{\text{core}}^{\ast }}.\end{array}$$(33)

The correction of V between the differentially and solid-body rotating cases is then:

$$\begin{array}{c}\hfill \frac{V_{\text{diff}}}{V}={\left[{\alpha }_{\text{rot}}^2(\frac{G(s)}{s}{)}^{-4/3}(\frac{{\mathrm{\Lambda }}_k^m(G(s))}{{\mathrm{\Lambda }}_k^m(s)}{)}^{-1}\right]}_{s_{\text{core}}^{\ast }}.\end{array}$$(34)

With that definition, we retrieve, dropping the env subscript for envelope spin parameters for readability:

$$\begin{array}{c}\hfill cot(\frac{{\pi }^{2}s}{{\mathrm{\Omega }}_{\text{env}}{\mathrm{\Pi }}_0}-\frac{\pi }{6})+\frac{1}{\sqrt{3}}\simeq -\frac{ϵ/V_{\text{diff}}}{G^{-1}(s)-G^{-1}{(s)}^{\ast }}.\end{array}$$(35)

To compute the period spacing expression, we write Eq. (35) for two neighbouring solutions s₁ and s₂ with s₂ > s₁. We assume that the frequency spectrum of the gravito-inertial modes is dense. The details of our computation are given in Appendix C. Using the vicinity in period of the solutions, we show that the Taylor expansion gives:

$$\begin{array}{cc}& [1+(\frac{ϵ/V_{\text{diff}}}{\overline{s}-s^{\ast }}{\left(\frac{\mathrm{d}{G}^{-1}}{\mathrm{d}s}{|}_{\frac{\overline{s}+s^{\ast }}{2}}\right)}^{-1}+\frac{1}{\sqrt{3}}{)}^2]\times [\frac{{\pi }^2(s_1-s_2)}{{\mathrm{\Omega }}_{\text{env}}{\mathrm{\Pi }}_0}-\pi ]\simeq \hfill \\ \hfill & -\frac{ϵ}{V_{\text{diff}}}\frac{s_1-s_2}{{(\overline{s}-s^{\ast })}^2}{\left(\frac{\mathrm{d}{G}^{-1}}{\mathrm{d}s}{|}_{\frac{\overline{s}+s^{\ast }}{2}}\right)}^{-2}\left(\frac{\mathrm{d}{G}^{-1}}{\mathrm{d}s}{|}_{\overline{s}}\right),\hfill \end{array}$$(36)

with $s^{\ast }=G(s_{\mathrm{core}}^{\ast })$ and $\overline{s}=(s_1+s_2)/2$. This equation can thus be easily compared with Equation (63) of Tokuno & Takata (2022). The effect of the differential rotation is comprised in the function G, through the derivative of G⁻¹ and V_diff. To get further insights we derive the approximate expression of ΔP as function of P. We define $\mathrm{P}=\pi \overline{s}/{\mathrm{\Omega }}_{\text{env}}$, ΔP = π(s₁ − s₂)/Ω_env, P^* = πs^*/Ω_env and ${\mathrm{\Gamma }}_{\text{diff}}=\frac{3\pi ϵ}{4{\mathrm{\Omega }}_{\text{env}}{V}_{\text{diff}}}$. This P^* holds the frequency dependence of the location of the dip in the frame co-rotating with the frequency Ω_env. We have:

$$\begin{array}{c}\hfill \frac{1}{\mathrm{\Delta }\mathrm{P}}-\frac{1}{{\mathrm{\Pi }}_0}\simeq \frac{{\mathrm{\Gamma }}_{\text{diff}}/\pi \frac{\mathrm{d}{G}^{-1}}{\mathrm{d}s}{|}_{\overline{s}}}{{((\mathrm{P}-{\mathrm{P}}_{\ast })\frac{\mathrm{d}{G}^{-1}}{\mathrm{d}s}{|}_{\frac{\overline{s}+s_{\ast }}{2}}+\frac{{\mathrm{\Gamma }}_{\text{diff}}}{\sqrt{3}})}^2+{\mathrm{\Gamma }}_{\text{diff}}^2}.\end{array}$$(37)

The dip in the period spacing is thus also approximatively a Lorentzian¹ near the period of the pure inertial mode in the co-rotating frame corresponding to the envelope. Then:

$$\begin{array}{cc}& \mathrm{\Delta }P\simeq {\mathrm{\Pi }}_0\hfill \\ \hfill & (1-\frac{{\mathrm{\Pi }}_0{\mathrm{\Gamma }}_{\text{diff}}/\pi \frac{\mathrm{d}{G}^{-1}}{\mathrm{d}s}{|}_{\overline{s}}}{{((P-P_{\ast })\frac{\mathrm{d}{G}^{-1}}{\mathrm{d}s}{|}_{\frac{\overline{s}+s_{\ast }}{2}}+\frac{{\mathrm{\Gamma }}_{\text{diff}}}{\sqrt{3}})}^2+{\mathrm{\Gamma }}_{\text{diff}}^2+{\mathrm{\Pi }}_0\frac{{\mathrm{\Gamma }}_{\text{diff}}}{\pi }\frac{\mathrm{d}{G}^{-1}}{\mathrm{d}s}{|}_{\overline{s}}}).\hfill \end{array}$$(38)

Naturally, by putting α_rot = 1, we retrieve the canonical inverted Lorentzian profile derived in Tokuno & Takata (2022).

3.5. Case of discontinuous rotation and density profiles

We take the same two-zone rotation profile as previously, and we assume the density at the interface core-radiative envelope to be discontinuous, which produces as well a discontinuous Brunt-Väisälä profile at the core-envelope boundary. This particular case is the one that should be used if the variation of the Brunt-Väisälä frequency near the core happens on a lengthscale that is less extended than the radial wavelength of the mode. First, this is seen in Ouazzani et al. (2020) for models near the Terminal Age Main Sequence at the upper range of mass of the γ Doradus classical instability strip. Second, as the extra-core mixing mechanisms suffer from uncertainties, it is relevant to encompass both continuous and discontinuous cases in our analysis. This problem is thus parametrized as such:

$$\begin{array}{c}\hfill {\displaystyle \underset{r\to R_{\text{core}}^{-}}{lim}}\rho ={\rho }_{\text{core}},\end{array}$$(39)

and

$$\begin{array}{c}\hfill {\displaystyle \underset{r\to R_{\text{core}}^{+}}{lim}}\rho ={\rho }_{\text{core}}+\mathrm{\Delta }\rho ={\rho }_{\text{env}},\end{array}$$(40)

as for the Brunt-Väisälä frequency:

$$\begin{array}{c}\hfill {\displaystyle \underset{r\to R_{\text{core}}^{-}}{lim}}N=0,\end{array}$$(41)

$$\begin{array}{c}\hfill {N|}_{R_{\text{core}}}=+\infty ,\end{array}$$(42)

and

$$\begin{array}{c}\hfill {\displaystyle \underset{r\to R_{\text{core}}^{+}}{lim}}N=N_0.\end{array}$$(43)

Two cases are comprised in this analysis: (1) only the derivative of the density is discontinuous, hence the Brunt-Väisälä frequency. This corresponds to Δρ = 0; (2) both the density and its derivative are discontinuous: Δρ ≠ 0. We emphasize that in this framework, the lower boundary of the gravito-inertial modes cavity and the radial coordinate of the limit of the convective core are assumed to be the same: r_a = R_core.

The derivation of the radial displacement and the Eulerian pressure perturbation from the upper side of the boundary is made in Appendix A.

Core solutions are the same as previously described in the case of constant density in the core. The matching equation for Lagrangian displacement thus becomes:

$$\begin{array}{cc}\hfill {\displaystyle \sum _k}{a}_k& \stackrel{\sim }{ϵ}{\stackrel{\sim }{X}}_k^m(s_{\text{env}}){\mathrm{\Theta }}_k^m(\mu ;s_{\text{env}})\hfill \\ \hfill & ={\displaystyle \sum _l}{b}_l{C}_l^m(1/s_{\text{core}}){\stackrel{\sim }{P}}_l^m(\mu ).\hfill \end{array}$$(44)

The pressure equation has to be precisely taken care of. As Tokuno & Takata (2022) highlighted, the match of the Lagrangian variation of pressure does not equate to the match of the Eulerian pressure anymore, since background quantities are no longer continuous near the interface. We thus derive the matching equation for $\delta p=p\prime +{\xi }_r(\frac{\partial \overline{p}}{\partial r})$. With the same reasoning as the one held in Section 3.1, we neglect the contribution of the centrifugal force with respect to the self-gravity. We rewrite the expression of the Lagrangian pressure perturbation as:

$$\begin{array}{c}\hfill \delta p=\overline{\rho }{c}_S^2\left(\frac{p\prime }{{\mathrm{\Gamma }}_1\overline{p}}\right)-r\overline{\rho }\overline{g}\left(\frac{\xi {\prime }_r}{r}\right),\end{array}$$(45)

using the hydrostatic equilibrium and the definition of the sound speed, $\overline{g}$ being the local background gravity acceleration. The expressions of the Lagrangian pressure perturbation in the two zones are thus, respectively, in the inertial frame

$$\begin{array}{cc}\hfill {\rho }_{\text{env}}& {\displaystyle \sum _k}{a}_k[{\sigma }_{\text{env}}^2{R}_{\text{core}}^2{\stackrel{\sim }{Y}}_k^m(s_{\text{env}})\hfill \\ \hfill & -\left(\frac{\mathcal{G}{M}_{\text{core}}}{R_{\text{core}}}\right)\stackrel{\sim }{ϵ}{\stackrel{\sim }{X}}_k^m(s_{\text{env}})]\times {\mathrm{\Theta }}_k^m(\mu ,s_{\text{env}})\hfill \end{array}$$(46)

in the envelope, and

$$\begin{array}{cc}\hfill {\rho }_{\text{core}}& {\displaystyle \sum _l}{b}_l[{\sigma }_{\text{core}}^2{R}_{\text{core}}^2{P}_l^m(1/s_{\text{core}})\hfill \\ \hfill & -\left(\frac{\mathcal{G}{M}_{\text{core}}}{R_{\text{core}}}\right)C_l^m(1/s_{\text{core}})]\times {\stackrel{\sim }{P}}_l^m(\mu )\hfill \end{array}$$(47)

in the core.

The approximated coupling equation for the pressure reads

$$\begin{array}{cc}& {\rho }_{\text{env}}{\displaystyle \sum _k}{a}_k[{\sigma }_{\text{env}}^2{\stackrel{\sim }{Y}}_k^m(s_{\text{env}})-\frac{\mathcal{G}{M}_{\text{core}}}{R_{\text{core}}^3}\stackrel{\sim }{ϵ}{\stackrel{\sim }{X}}_k^m(s_{\text{env}})]\times {\mathrm{\Theta }}_k^m(\mu ,s_{\text{env}})\hfill \\ \hfill & ={\rho }_{\text{core}}{\displaystyle \sum _l}{b}_l[{\sigma }_{\text{core}}^2{P}_l^m(1/s_{\text{core}})-\frac{\mathcal{G}{M}_{\text{core}}}{R_{\text{core}}^3}{C}_l^m(1/s_{\text{core}})]\times {\stackrel{\sim }{P}}_l^m(\mu ).\hfill \end{array}$$(48)

We follow the same method as previously, projecting both coupling equations on Hough functions. The matrix equation has the same form as Eq. (19):

$$\begin{array}{c}\hfill [\stackrel{\sim }{}\}M(s_{\text{core}},s_{\text{env}})-\stackrel{\sim }{ϵ}\stackrel{\sim }{}\}N(s_{\text{core}},s_{\text{env}})]\mathit{b}=\mathbf{0},\end{array}$$(49)

with

$$\begin{array}{c}\hfill \stackrel{\sim }{}\}M=c_{k,l}{\sigma }_{\text{env}}^2{\rho }_{\text{env}}{\stackrel{\sim }{Y}}_k^m(s_{\text{env}})C_l^m(1/s_{\text{core}})\end{array}$$(50)

and

$$\begin{array}{cc}\hfill \stackrel{\sim }{}\}N=& c_{k,l}({\sigma }_{\text{core}}^2{\rho }_{\text{core}}{\stackrel{\sim }{X}}_k^m(s_{\text{env}})P_l^m(1/s_{\text{core}})\hfill \\ \hfill & +\mathrm{\Delta }\rho \frac{\mathcal{G}{M}_{\text{core}}}{R_{\text{core}}^3}{\stackrel{\sim }{X}}_k^m(s_{\text{env}})C_l^m(1/s_{\text{core}})).\hfill \end{array}$$(51)

The leading term analysis conducted as before leads to the synthetic coupling equation

$$\begin{array}{cc}& \frac{{\sigma }_{\text{env}}^2}{{\sigma }_{\text{core}}^2}(1+\frac{\mathrm{\Delta }\rho }{{\rho }_{\text{core}}})cot(\frac{{\pi }^2{s}_{\text{env}}}{{\mathrm{\Omega }}_{\text{env}}{\mathrm{\Pi }}_0}-\frac{\pi }{4})\times [\frac{s_{\text{env}}}{2{\mathrm{\Lambda }}_k^m{(s_{\text{env}})}^{1/2}}]F_l^m(s_{\text{core}})\hfill \\ \hfill & \simeq \stackrel{\sim }{ϵ}(1-\frac{\mathrm{\Delta }\rho }{{\rho }_{\text{core}}}\frac{\mathcal{G}{M}_{\text{core}}}{{\sigma }_{\text{core}}^2{R}_{\text{core}}^3}{F}_l^m(s_{\text{core}})).\hfill \end{array}$$(52)

This equation, as before, can be computed either numerically or analytically, by an expansion of $F_l^m$ around its zero $s_{\text{core}}^{\ast }$. For the analytical expansion, we rely on the function G previously described. The coupling equation reduces to:

$$\begin{array}{cc}& cot(\frac{{\pi }^{2}s}{{\mathrm{\Omega }}_{\text{env}}{\mathrm{\Pi }}_0}-\frac{\pi }{4})+\stackrel{\sim }{ϵ}{\left(\frac{{\alpha }_{\text{rot}}{G}^{-1}(s)}{2{\mathrm{\Omega }}_{\text{env}}}\right)}^2\frac{\mathcal{G}{M}_{\text{core}}}{R_{\text{core}}^3}\frac{\mathrm{\Delta }\rho }{{\rho }_{\text{env}}}\hfill \\ \hfill & \times [\frac{2{\mathrm{\Lambda }}_k^m{(G(s))}^{1/2}}{{\alpha }_{\text{rot}}^{2}G(s)}(\frac{G(s)}{s}{)}^2{]}_{s_{\text{core}}^{\ast }}\simeq -\frac{\stackrel{\sim }{ϵ}/{\stackrel{\sim }{V}}_{\text{diff}}}{G^{-1}(s)-G^{-1}{(s)}^{\ast }},\hfill \end{array}$$(53)

with

$$\begin{array}{cc}\hfill {\stackrel{\sim }{V}}_{\text{diff}}=& -(1+\frac{\mathrm{\Delta }\rho }{{\rho }_{\text{core}}})\times [\frac{\mathrm{d}{F}_l^m(s)}{\mathrm{d}s}\frac{{\alpha }_{\text{rot}}^{2}G(s)}{2{\mathrm{\Lambda }}_k^m{(G(s))}^{1/2}}(\frac{s}{G(s)}{)}^2{]}_{s_{\text{core}}^{\ast }}.\hfill \end{array}$$(54)

We verify that putting α_rot = 1, and thus G = 𝟙, reduces the coupling equation to equations (88) and (89) of Tokuno & Takata (2022). The correction from the solid-body rotating case is thus

$$\begin{array}{c}\hfill \frac{{\stackrel{\sim }{V}}_{\text{diff}}}{\stackrel{\sim }{V}}=[{\alpha }_{\text{rot}}^2{\left(\frac{G(s)}{s}\right)}^{-1}{\left(\frac{{\mathrm{\Lambda }}_k^m(G(s))}{{\mathrm{\Lambda }}_k^m(s)}\right)}^{-1/2}{]}_{s_{\text{core}}^{\ast }}.\end{array}$$(55)

With the same expansion of the cotangent as before, and ${\stackrel{\sim }{\mathrm{\Gamma }}}_{\text{diff}}=\frac{\pi \stackrel{\sim }{ϵ}}{{\mathrm{\Omega }}_{\text{env}}{\stackrel{\sim }{V}}_{\text{diff}}}$, we get:

$$\begin{array}{c}\hfill \frac{1}{\mathrm{\Delta }\mathrm{P}}-\frac{1}{{\mathrm{\Pi }}_0}\simeq \frac{{\stackrel{\sim }{\mathrm{\Gamma }}}_{\text{diff}}/\pi \frac{\mathrm{d}{G}^{-1}}{\mathrm{d}s}{|}_{\overline{s}}}{{\left((\mathrm{P}-{\mathrm{P}}_{\ast })\frac{\mathrm{d}{G}^{-1}}{\mathrm{d}s}{|}_{(\overline{s}+s_{\ast })/2}\right)}^2+{\stackrel{\sim }{\mathrm{\Gamma }}}_{\text{diff}}^2}.\end{array}$$(56)

4. Results and discussion

We apply in this section our theoretical results to a particular Kelvin g-mode for a typical fast-rotating γ Doradus star and investigate the detectability of the core-near-core differential rotation from asteroseismic data. To do so, we first analyse the effect of differential rotation on the structure and location of the inertial dips in the period-spacing pattern, in the formalism of a continuous near-core Brunt-Väisälä frequency. We then aim to assess the potential uncertainties on the differential rotation rate for different processes: (1) the profile of the near-core Brunt-Väisälä frequency (2) the intrinsic uncertainty coming from the limited time duration of observations (3) the influence of density stratification in the core. We conclude on the detectability of the differential rotation at the core-near-core interface and on the convective core rotation rate from the dip properties. Caveats of the formalism and perspectives are ultimately discussed.

4.1. Case star: Typical values of buoyancy travel time and rotation rate

We select a typical fast-rotating γ Doradus star, with an envelope rotation rate of ${\mathrm{\Omega }}_{\mathrm{env}}/2\pi =2.16{\mathrm{d}}^{-1}$, a buoyancy travel time of Π₀ = 4175 s. Those values are close to the ones observed in Saio et al. (2021) for KIC12066947: (2.159 ± 0.002 d⁻¹, 4175 ± 28 s respectively), which is a good representative of the γ Doradus in their study. These parameters are compatible with stellar modelling prescriptions, such as the ones presented in Ouazzani et al. (2019). As we are taking an assumption of uniform density in the core to maintain the analyticity of the results, the star best fitting our model is a ZAMS star, as the density gradient in the core is increased during evolution (Ouazzani et al. 2020). KIC12066947 was not found to be a ZAMS star by Saio et al. (2021), having a core hydrogen abundance of 0.53 in the model computed without overshooting, yet we keep the values of Π₀ and Ω_env as they remain close to the modelling prescriptions for ZAMS stars. This age dependence will be further discussed in Section 4.6.

To further explore the coupling problem in the presence of differential rotation, we study the interaction between a pure inertial mode (l = 3, m = −1) and a Kelvin g-mode (k = 0, m = −1) in the envelope. In this case,

$$\begin{array}{c}\hfill F_l^m(s)=-\frac{s^2-10s-15}{(s+1)(s^2-5)},\end{array}$$(57)

and the spin parameter of the pure inertial mode is s_core^* = 11.3245. The choice of this particular mode is motivated by the fact that: (1) patterns of Kelvin g-mode (k = 0, m = −1) are the most commonly detected in the 611 γ Doradus sample of Li et al. (2020), forming 62% of the detected period-spacing patterns; (2) as shown by Ouazzani et al. (2020), the pure inertial (l = 3, m = −1) mode is likely to couple with this Kelvin g-mode, having a high geometrical factor of c_0, 3 = 0.5 in the absence of differential rotation, remaining non-negligible with differential rotation, as investigated in Appendix C. Dips arising from this interaction have previously been observed in Saio et al. (2021) in the range s_env ∈ [8; 11]. The high spin parameter of the pure inertial mode ensures our formalism holds in that case, which will be later discussed in Section 4.7.

For this analysis, we rely on Aerts & Mathis (2023) to assess the relevant range of coupling parameters allowed in our study. Among the 37 γ Doradus analysed, the maximum value of coupling parameter is ϵ ≈ 0.25, and $\stackrel{\sim }{ϵ}\approx 0.1$, which corresponds respectively to Γ ≈ 40 h and $\stackrel{\sim }{\mathrm{\Gamma }}\approx 12$ h. Given the reduced amount of analysed stars and the potentiality of further detection of coupling parameters superior to this threshold, we allow Γ and $\stackrel{\sim }{\mathrm{\Gamma }}$ parameters to vary in the range [1; 60] hours.

As for the amount of differential rotation allowed in the present study, we base ourselves on the detection of core-near-core differential rotation made by Saio et al. (2021). In their sample, the highest core-near-core differentiality is reached for KIC 05985441, with α_rot ≈ 0.85. Conversely, in the opposite regime of the core rotating slower than the envelope, the highest value of differentiality is α_rot ≈ 1.01 (KIC 8330056). The regime of α_rot < 1 is favored by realistic 2D models able to reproduce the meridional circulation in the radiative envelope of early-type fast rotating stars and to predict the gravity-darkening observed in such stars (Espinosa Lara & Rieutord 2013). We thus extend the range of analysed differential rotation to α_rot ∈ [0.80; 1.05]. The different ranges of the parameters considered are summarized in Table 1.

4.2. Effects of differential rotation on the morphology of the inertial dip

In Fig. 3, the numerical results described in subsection 3.3, obtained by solving Eq. (28) are overplotted on top of the analytical Lorentzian profile described in Eq. (38). Three different values for the coupling parameter ϵ are used, corresponding to Γ values in the solid-body rotating case of 1, 4, and 12 hours. These values are both taken for the sake of comparison with the results of Tokuno & Takata (2022) in the case of the lower coupling parameter values, as well as to illustrate the case of a moderately high coupling parameter comprised in the present study. To follow the evolution of the morphology of the dip with differential rotation, cases from α_rot = 0.8 to α_rot = 1.05 are considered in the same figure.

From Fig. 3, we first confirm the consistency of the numerical results (dots) with the approximated modified Lorentzian profile (lines), as already predicted by the difference in O(ϵ²) arising from the fact that we take into account the first term in the development of $F_l^m$ near its zero $s_{\mathrm{core}}^{\ast }$ in Section 3.4. The dips, as described in the numerical study held in Section 3.3, arise at an envelope spin parameter varying with the amount of core to near-core differential rotation. From no differential rotation and a spin parameter of $s_{\text{env}}^{\ast }=s_{\text{core}}^{\ast }=11.32$, a differential rotation of α_rot = 0.8 gives $s_{\text{env}}^{\ast }=4.25$. Conversely, a differential rotation of α_rot = 1.05 gives $s_{\text{env}}^{\ast }=16.59$. This corresponds to a shift in the period at which the dip arises in Fig. 3. High core-near-core differentiality corresponds to a low period of the dip, while a faster-rotating envelope increases the period of the dip when compared to the case of uniform rotation. At fixed coupling parameter, with increasing core rotation compared to the envelope, the dip gets also thinner and deeper, being formed of fewer modes. The value of Γ_diff decreases from Γ. On the contrary, in regimes where the core rotates slower than the envelope, the dip gets wider and shallower, with many modes composing the dip. The value of Γ_diff increases from Γ.

The middle and bottom panels of Fig. 3 show the evolution of the dip structure with increasing coupling parameter ϵ and Γ parameter. The results obtained by Tokuno & Takata (2022) in the case of no differential rotation still hold in differentially rotating regimes: at fixed α_rot parameter, with increasing coupling parameter the dip gets wider, and is shifted to lower periods. We see that in the case of Γ = 12 hours, dips in the regime of core rotation inferior to envelope rotation are barely distinguishable from the baseline, which in turn questions the detectability of such feature in data.

The effect of an increasing core-to-envelope differential rotation on the dip shape is thus comparable to a decrease of the coupling parameter. This effect can be physically interpreted using Fig. 2, considering the mode interaction in the reference frame co-rotating with the core. Seen from this frame, the spectrum of gravito-inertial modes gets denser with increasing envelope rotation, sparser with decreasing one. Taking for instance the case α_rot = 1.05, the intersection of the yellow continuous line with blue thin ones forms modes that are separated by a low spacing in terms of core spin parameter compared to the solid-body rotating one (red line). Conversely in the case α_rot = 0.80, the lower inclination of the black line with respect to the blue ones forms modes spaced by an increased spacing in core spin parameter. Thus, in the regime of α_rot > 1, the pure inertial mode is close in period to an increased number of modes, with which it can interact.

We quantify this effect in Table 2. If $s_{\mathrm{core}}^{\min }$ and $s_{\mathrm{core}}^{\max }$ are respectively the lowest and highest value of core spin parameter for which the resulting period spacing is inferior to 98% of the baseline value, Δn is the number of modes in the interval $[s_{\mathrm{core}}^{\min },s_{\mathrm{core}}^{\max }]$. $\mathrm{\Delta }{s}_{\mathrm{core}}=s_{\mathrm{core}}^{\max }-s_{\mathrm{core}}^{\min }$ is the gap in core spin parameter, and Δn/Δs_core a measurement of the density of gravito-inertial modes close to the period of the core pure inertial mode in the frame co-rotating with the core. This density, at fixed ϵ parameter, increases with decreasing core-to-envelope differential rotation, hence increasing α_rot parameter.

This increased coupling efficiency explains the similar effect between a high α_rot and a high coupling parameter ϵ on the dip shape. Conversely, at α_rot < 1 the pure inertial mode is surrounded by fewer modes and the coupling efficiency is less important, resulting in a thin dip.

This effect is captured in Eq. (35) by the parameter ϵ/V_diff, controlling the depth and the width of the inertial dip structure. The coupling parameter ϵ is left unchanged by the existence of differential rotation, but V_diff increases with increasing core-to-envelope differential rotation (or decreasing α_rot). This evolution is due to the different angular structure of the gravito-inertial modes, as at fixed core spin parameter the pure inertial mode keeps the same structure as in the solid-body rotating case while the Hough function ${\mathrm{\Theta }}_k^m(\mu ,s_{\mathrm{env}})$ varies with the considered envelope spin parameter. Overall, we retrieve the influence of increasing core-to-envelope differential rotation, decreasing the coupling efficiency measured by ϵ/V_diff.

4.3. Detectability of the core-near-core differential rotation

In the left panel of Fig. 4, we summarize the results by plotting the period of the inertial dip minimum with regards to its half-width in the co-rotating frame related to the envelope, for the extent of parameter α_rot and Γ described in Table 1. The amount of differential rotation is given by the colour scale. The computations for this left panel are made in the framework of a continuous Brunt-Väisälä profile near-core.

4.3.1. Validity of the theoretical expansion

Two dotted lines are plotted in the left panel of Fig. 4. Below those lines, the ratio ${ϵ}^3/s_{\text{env}}^2$ is respectively less than 0.1 and 0.01. This ratio controls the accuracy of the developments made for the expansion of the gravito-inertial modes to the edge of the convective core, as described by Eq. (A.18), in the framework of a continuous near-core Brunt-Väisälä frequency. We see that this ratio, for the range of analysed parameters, reaches high values in the regime in which the core rotates faster than the envelope, due to the decrease in the envelope spin parameter at which the interaction occurs, $s_{\text{env}}^{\ast }$. We advise against the use of ours and Tokuno & Takata (2022)’s formalisms in the regime in which ${ϵ}^3/s_{\text{env}}^2>0.1$, as the extent of the intermediate region between R_core and r_a would make the structure of the gravito-inertial modes significantly altered from Hough functions. Inertial dips would still be detected but further theoretical developments of the wave propagation in this intermediate-region, or numerical studies (e.g. Dintrans & Rieutord 2000) would be necessary in this case. However, this ratio remains small in the majority of the regions studied, ensuring good validity of our model throughout our analysis.

4.3.2. Effect of noise in the asteroseismic data on the detectability of the dip

To address sources of difficult dip structure detection, and thus uncertainty on the detection of differential rotation from the inertial dip structure study, we aim to trace a zone in which the analysis of an inertial dip in data will lead to an accurate estimate of the core to near-core differential rotation.

For this, we first construct perturbed period-spacing patterns representative of the number of modes usually found in a period-spacing pattern and their periods, taking into account a realistic estimate of the intrinsic observational uncertainty. We artificially introduce a dip in the period-spacing pattern. The generation of realistic period-spacing patterns is thoroughly described in Appendix D. This procedure is reproduced for values in the range 0.9 < α_rot < 1.05 and 0.5h < Γ < 60h. We then infer uncertainties on the parameters (α_rot, Γ) by running a Markov-chain Monte Carlo analysis on these simulated period-spacing patterns.

An open question is then of how many parameters would the fit to data best depend on. In the case of the dips found in Saio et al. (2021), their clear appearance in the period-spacing pattern suggests that a first estimation of (Ω_env, Π₀) can be performed excluding the modes belonging to the inertial dip. This was originally done in Van Reeth et al. (2016) for the star KIC 12066947, where the estimation of Ω_env and Π₀ was done on the low-period part of the spectrum, while the region containing the dip was excluded when fitting the baseline. The global parameters can then be used to move to the co-rotating frame and detect the parameters which the dip shape and location depend on, ϵ (through Γ), and α_rot. This analysis can a priori not be done for all dips, as we saw in Section 4.2 that in the regime of high coupling, or α_rot > 1, the dips can get shallow, hence uneasy to clearly distinguish from non-coupling gravito-inertial modes, and affect the measure of Ω_env and Π₀. We then decide to investigate the detectability of core rotation effects when performing a multidimensional fit, taking into account simultaneously the previously inferred quantities (Ω_env, Π₀) and the ones characterizing the dip (α_rot, Γ) in the vector Θ = (Ω_env, Π₀, α_rot, Γ, P₀), P₀ being the period in the co-rotating frame of the first mode in the period-spacing pattern considered. We consider the period-spacing pattern in the inertial frame, as we aim to fit as well the internal rotation rate.

To probe the parameter space, we build custom period-spacing patterns consisting of a series of periods $[{\mathrm{P}}_{k,\mathrm{in}}^{\mathrm{mod}}]$ in the inertial frame, depending on the five aforementioned parameters, following the differentially rotating model (hereafter “5D model”), in the framework of a continuous near core Brunt-Väisälä profile:

$$\begin{array}{c}\hfill {\mathcal{F}}_{5\mathrm{D},\mathrm{cont}}:({\mathrm{\Omega }}_{\mathrm{env}},{\mathrm{\Pi }}_0,{\mathrm{P}}_0,{\alpha }_{\mathrm{rot}},\mathrm{\Gamma })\to [{\mathrm{P}}_{k,\mathrm{in}}^{\mathrm{mod}}].\end{array}$$(58)

We compare it to the period-spacing pattern $[{\mathrm{P}}_{k,\mathrm{in}}^{\mathrm{pert}}]$ that was originated from a collection of “true” parameters, ${\mathbf{\Theta }}_{\mathrm{true}}=({\mathrm{\Omega }}_{\mathrm{env}}=2\pi \times 2.16{\mathrm{d}}^{-1},{\mathrm{\Pi }}_0=4175\mathrm{s},{\mathrm{P}}_{0,\mathrm{true}},{\alpha }_{\mathrm{rot}},\mathrm{\Gamma }))$ the three former ones being fixed to the values found in KIC12066947, and further perturbed to accurately represent realistic Kepler data following a procedure described in Appendix G. We use the following loglikelihood:

$$\begin{array}{c}\hfill log\mathcal{L}(\mathcal{F})=-\frac{1}{2}{\displaystyle \sum _k}\frac{{({\mathrm{P}}_{k,\mathrm{in}}^{\mathrm{mod}}-{\mathrm{P}}_{k,\mathrm{in}}^{\mathrm{pert}})}^2}{{\sigma }_k^2}.\end{array}$$(59)

The standard deviation ${\sigma }_k^2$ is the observational error we hypothesize, taking a fixed signal-over-noise ratio for all modes contained in the period-spacing pattern. This approach is similar to the one taken in e.g. Moravveji et al. (2015) and adapted to our case in which we model perturbations from the observational noise. We emphasize that even though this error represents modes with a fixed S/N ratio, the values of σ_k can be derived for a realistic case of varying S/N over modes. The form of this likelihood is discussed in Appendix G, as well as an extension to account for baseline variations that are not comprised in this analytical model.

We aim to compare the accuracy of the differentially rotating model with the solid-body rotating one to fit the simulated period-spacing patterns containing an inertial dip. For this, we perform another analysis taking into account the following nested model (hereafter described as “4D model”), fitting it on the same simulated data as the 5D model:

$$\begin{array}{c}\hfill {\mathcal{F}}_{4\mathrm{D},\mathrm{cont}}:({\mathrm{\Omega }}_{\mathrm{env}},{\mathrm{\Pi }}_0,{\mathrm{P}}_0,\mathrm{\Gamma })\to [{\mathrm{P}}_{k,\mathrm{in}}^{\mathrm{mod}}].\end{array}$$(60)

The MCMC is then computed using the emcee package (Foreman-Mackey et al. 2013). Examples of cornerplots resulting from the 5D and 4D models are shown in Appendix E. We first perform a burn-in step of 100 iterations, taking walkers from a uniform prior described in Table 3. The positions of the walkers are then reinitialized, sampling them from a normal distribution around the set of parameters found to maximize the loglikelihood in the burn-in phase. The MCMC is then run until 100 times the auto-correlation time is reached, with a limitation of 20 000 iterations to get the posterior distributions for the parameters. If this limit on the number of iterations has been reached and is inferior to 50 times the auto-correlation time, the detection of the dip is considered “challenging”: the parameters are insufficiently constrained by the fit.

From the posterior distribution, when available, we are interested in the marginal distributions of the parameters α_rot and Γ, which give us uncertainties on the input parameters of our model. Then remains to evaluate for which set of parameters (α_rot, Γ) we can consider the 5D model presented in this work to significantly better fit the data than the 4D solid-body rotating model with only Γ as a parameter. We thus aim to test the following null hypothesis: “The differentially rotating and solid-body rotating model fit the data equally well”. If this hypothesis is true, the detection of differential rotation is not considered significant. To test the hypothesis, we calculate the p-value for the null model using a ${\chi }_{df=1}^2$ cumulative distribution function:

$$\begin{array}{c}\hfill \mathrm{p}=P[{\chi }_{df=1}^2>\lambda ].\end{array}$$(61)

With λ the difference of the loglikelihoods:

$$\begin{array}{c}\hfill \lambda =2[log\mathcal{L}({\mathcal{F}}_{5\mathrm{D},\mathrm{cont}})-log\mathcal{L}({\mathcal{F}}_{4\mathrm{D},\mathrm{cont}})].\end{array}$$(62)

The analysis described here is extended to the case of models derived taking a discontinuous near-core Brunt-Väisälä profile, the 5D and 4D models being, in that case, ℱ_5D, disc and ℱ_4D, disc. We detail hereafter our analysis of the results obtained in the framework of continuous Brunt-Väisälä profile, as the overall conclusions remain qualitatively the same in the discontinuous case, and bridge the two results in Section 4.4.

Results are found in Fig. 5: the zones for which p < 0.05 are marked in blue, indicating that we can reject the null hypothesis with a level of confidence of 95%. This treatment ensures that the detection is not subject to the noise. In other cases, the detection of differential rotation is not significant enough to stand out from the solid-body rotating case, and the region is marked in orange. These zones are then inscribed in Fig. 4. This analysis is done both for the models with continuous or discontinuous Brunt-Väisälä near the boundary.

In the regime close to solid body rotation (0.99 < α_rot < 1.01), we notice the differential rotation to be detectable even for values of α_rot different from 1 by less than 1%, depending on the considered Γ parameter. The dip being more extended and shallower for high Γ values, this threshold of detection is reduced with increasing values of Γ. For values inferior to 15 hours, this translates to a zone of 0.1 days around the minimum period of the dip in the corotating frame for which differential rotation would be impossible to disentangle from the solid-body rotating case. We however point out that a precision of less than a percent in this regime shows the great potential of taking the dip structure as a probe for core-to-envelope differential rotation.

When the core rotates slower than the envelope (α_rot  >  1), the detection is limited by the shallowness of the dip, as well as by the important observational noise at high periods. For values of Γ > 12 hours (in the case of n = 33), the null hypothesis cannot be rejected at a 95% level of confidence, indicating an impossibility of reliable detection of differential rotation from the analysis of a single inertial dip in data. In our particular case, above α_rot = 1.02 and Γ > 4 hours, the fitting shows no satisfactory results with dip parameters unconstrained in a reasonable time due to the high observational noise in the high-period part of the pattern in which the dip is situated. We thus consider the fit impossible in realistic data in the regime of a core rotating 2% slower than the envelope, all the more if we consider slower rotators. This analysis provides the upper limit of the detectability region at large periods reported on the left panel of Fig. 4, by choosing the case where n = 33.

In the regime in which the core rotates faster than the envelope (α_rot < 1), the detection is eased by the depth of the dip and the low impact of the noise in the region of the pattern in which the dip appears. We find the dip to be detectable for high values of the coupling parameter in this regime.

To give further insights, we consider the only study of the distribution of coupling parameter on a sample of γ Doradus: 37 stars analysed in Aerts & Mathis (2023). In their sample of Kelvin g-mode series (k = 0, m = −1), the mean value of radial order separation is 33, and the mean number of modes 29, not far from the number of modes found in the sample of Saio et al. (2021) (33). The maximum value of the coupling parameter found is 0.25. The majority of the stars in the sample do not exhibit clear inertial dips, with only two common to the analysis of Saio et al. (2021). The authors state that those stars do not exhibit a value of coupling parameter different from the rest of the sample. We thus take this distribution of coupling parameters among “regular” γ Doradus stars to be representative of a distribution of coupling parameters among γ Doradus stars on which inertial dips can be fitted and a measurement of differential rotation made. According to Fig. 5, if all of the stars had a detected core to near-core differential rotation of more than 6% (or α_rot ≤ 0.94) and n = 33 modes in their period-spacing pattern, the differentially rotating model would significantly better fit their period-spacing patterns than the solid-body rotating model. The detection of differential rotation would be considered firm, with a confidence level of 95%. If the stars had a detected differential rotation of 1% (α_rot = 0.99), still with n = 33 modes in their period-spacing patterns, only stars with ϵ < 0.16 would fall in the region. In the sample of Aerts & Mathis (2023), this would translate to 25 stars having firm detection of differential rotation.

Interestingly, the limit of detectability of the dip in terms of Γ parameter increases with increasing core-to-envelope rotation. Indeed, at high coupling parameter, the dip is shallow and an increase of the rotation of the core with fixed envelope rotation deepens the inertial dip, as seen in the bottom panel of Fig. 3: the dip in the regime of core rotating faster than the envelope (in black for α_rot = 0.8) is deepened compared to the solid-body rotating case (in red). This is an interesting property, as Aerts & Mathis (2023) stated that among the fast rotators analysed sample, γ Doradus stars exhibiting dips do not stand out in term of their coupling parameter. We see here that in the regime of moderate to high coupling parameters, the presence of core-near-core differential rotation is one of the potential mechanisms rendering the dip clearer in the period-spacing pattern.

Detectability limits taking hypotheses of lower (n = 21) or higher (n = 50) number of modes in the period-spacing pattern are also shown in Fig. 5, as respectively solid and dashed-dotted lines. As expected, a higher number of modes in the period-spacing pattern allows for a better constraint on the relevant parameters of the fit, related to the properties of the radiative envelope (Ω_env, Π₀) or the inertial dip (α_rot, Γ). This pushes the detectability limit towards high coupling values, and conversely low coupling values for fewer modes.

We emphasize that even if the dip is detectable, one should be careful when studying dips in the regime ${ϵ}^3/s_{\mathrm{env}}^2>0.1$, as explained in Section 4.3.1. Lines under which ϵ³/s_env² < 0.1 and ϵ³/s_env² < 0.01 are shown on the left panel of Fig. 6. This should encourage to be cautious using the analytical model described in the present study for dips situated at low period values in the co-rotating frame related to the envelope, in the framework of a continuous Brunt-Väisälä frequency near the boundary.

We emphasize that this analysis assumes a constant value for the period spacing outside of the inertial dip. However, chemical and thermal structure features leave modulations in the pattern (Miglio et al. 2008; Cunha et al. 2015; Degroote et al. 2010). The detectable zones described hereby thus correspond to the ideal case in which those modulations are accurately described from stellar modelling. We discuss the impact of chemical stratification on the detectability of the dip in Appendix G.

Additionally, due to the extended propagation region in the convective core, modes appearing in the dip are likely to have a low signal-over-noise ratio and thus can be unidentified in the period-spacing pattern. Such a phenomenon, with the extreme case of having dips appearing as gaps in the period-spacing pattern, will reduce the detectability region of differential rotation from the case studied here. This effect will be most prominent in the regions in which the dip is composed of a few modes, i.e. at low α_rot and coupling parameters. We expect on the contrary dips in the α_rot > 1 detectability region to be still detectable with a few modes missing, as the dip extends over tens of modes in this regime even for low coupling parameter values.

For hypothetically strong core-to-envelope differential rotation below α_rot = 0.9, we chose not to pursue the MCMC analysis otherwise performed. In our routine of building perturbed period-spacing patterns, the detectability in this regime is not limited by the shallowness of the dip or the intrinsic observational noise hypothesized but by the presence or absence of the dip in the period-spacing pattern. Indeed the spin parameter value at which the dip appears, s_env^* is shifted towards values uncomprised in the period-spacing pattern analysed. This is readily found in the case of α_rot = 0.90 for which the MCMC fitting does not give constraints on α_rot and Γ, as the period of the first mode in the period-spacing pattern is too high compared to periods at which the effect of the interaction is sizeable. Moreover, in this regime, the dip is formed of a few modes and the results would be vastly different if the signal-to-noise ratio of the few modes detected in the dip would significantly depart from the hypothesis we made for this analysis. Thus as both the extent of the period-spacing pattern in terms of spin parameter (see Appendix D) and the Signal-to-Noise ratio of the modes in the dip are varying in each particular case, we find our analysis not to be applicable for this regime. Conversely, we point out that high core-to-envelope differential rotation is only to be found in stars displaying a low spin parameter for their mode of lowest radial order displayed in the period-spacing pattern.

Our analysis tackled the impact of the uncertainty intrinsic to the limited time duration of the observations of the Kepler mission. In the following sections, we investigate the influences of other physical processes or uncertainties that could hinder the detection of differential rotation, which are more difficult to estimate because of their variability across γ Doradus stars.

4.4. Degeneracy with the treatment of the stratification gradient near the core: Evolution constraints

We saw in Section 3.4 that by taking a continuous Brunt-Väisälä profile between the core and the envelope, the Lorentzian is shifted from the location of the pure inertial mode in the co-rotating frame related to the envelope and gets its minimum at ${\mathrm{P}}_{\min }={\mathrm{P}}_{\ast }-{\mathrm{\Gamma }}_{\mathrm{diff}}/\sqrt{3}{(\frac{\mathrm{d}{G}^{-1}}{\mathrm{d}s}{|}_{\frac{s^{\ast }+\overline{s}}{2}})}^{-1}$, whereas the minimum is reached at P_min = P_* in the discontinuous case. The two continuous and discontinuous cases, represented in Figure 2 of Tokuno & Takata (2022) can be considered as two extreme cases and thus represent a source of uncertainty on differential rotation, which also shifts the period at which the pure inertial mode forms a dip in the period-spacing pattern. Such a dependency is represented in Fig. 6, in the case of the model star exhibiting a dip with parameters $\mathrm{\Gamma }=\stackrel{\sim }{\mathrm{\Gamma }}=1.0$h, and a parameter α_rot = 0.95. We see that compared to the continuous treatment (in blue), the dip arising from the discontinuous treatment (in red) is shifted to high periods, as hinted from the theoretical model. This effect is due to the presence of an intermediate region between R_core and r_a in the case of a continuous Brunt-Väisälä frequency, and is thus the same as a decrease in core rotation compared to the envelope rotation. The shift in terms of envelope spin parameter, from the expression of P_min in the continuous Brunt-Väisälä frequency framework, can be estimated as:

$$\begin{array}{cc}\hfill \delta s_{\mathrm{cont}\to \mathrm{disc}}=& 6.7\times {10}^{-4}\left(\frac{{\mathrm{\Omega }}_{\text{env}}}{2\pi \times 2.16{\mathrm{d}}^{-1}}\right)\left(\frac{\mathrm{\Gamma }}{1\mathrm{h}}\right)\hfill \\ \hfill & \times {\alpha }_{\mathrm{rot}}{(1+m(1/{\alpha }_{\mathrm{rot}}-1)s_{\mathrm{core}}^{\ast })}^2.\hfill \end{array}$$(63)

In the regime of typical fast rotators, with deep dips and core rotating faster than the envelope, the uncertainty coming from the discontinuity of the Brunt-Väisälä profile is thus under-dominant compared to the shift of the spin parameter due to differential rotation. In the precise case of Fig. 6, the shift in the spin parameter of the location of the dip that is due to the discontinuous treatment of the Brunt-Väisälä is equivalent to a decreased value of α_rot by 3 × 10⁻⁶ from α_rot = 0.95. However, the influence of this source of uncertainty becomes appreciable for high values of Γ (hence high coupling). The right panel of Fig. 4 illustrates the different values of α_rot forming an equivalent dip (same location of the minimum in period and half-width in the co-rotating frame related to the envelope): the colour map shows the value α_rot in the discontinuous case in these coordinates, whereas the contour overplotted shows the value of α_rot in the continuous case as fully described in the left panel. For instance, for the model star exhibiting a dip at P_min = 1.93 d, with a half-width of 1.38 d, the continuous formalism gives a differential rotation of α_rot = 0.94, whereas the discontinuous formalism gives α_rot = 0.98, hence a relative error of 4%.

In the right panel of Fig. 4, we superimpose detectability zones inferred from the same treatment as the one described in Appendix D, in the framework of a discontinuous near-core Brunt-Väisälä frequency, and n = 33 modes visible in the period-spacing pattern. Results are consistent with the continuous case, in the α_rot < 1 regime, with differential rotation being detectable for higher values of Γ as the dip increases in amplitude and is shifted to low periods. The situation is slightly different in the regime of α_rot > 1, with a differential rotation only found to be detectable for α_rot ≤ 1.02 and a low value of Γ. We interpret it as an incidence of the previously highlighted shift of the period of the minimum of the dip with increasing Γ in the framework of the continuous near-core Brunt-Väisälä frequency, resulting in a detection of dip less attained at fixed α_rot between the two regimes.

The detectability zone, considering the degeneracy of the treatment of the near-core Brunt-Väisälä profile, can thus be considered as the intersection between the detectability zones inferred in the two cases of continuous and discontinuous Brunt-Väisälä profile. This effectively strengthens the limits in terms of accessible coupling parameter values for which differential rotation is considered detectable.

It is worth noting that the building of chemical gradient near-core during the evolution would hint towards a better accuracy of the discontinuous case for evolving main-sequence stars after the ZAMS. However, given the complex dependence of the Brunt-Väisälä profile on age or mixing parameters, one cannot state the most realistic case and reduce the uncertainty related to the Brunt-Väisälä profile near-core, unless a thorough forward-modelling study is led, taking into account the most up-to-date stellar modelling prescriptions. As it is not the aim of this present paper, we do not consider any bias towards one case or the other in the present study.

4.5. From differentiality to the measure of convective core rotation

We tackled two distinct internal structures: a continuous treatment of the near-core Brunt-Väisälä profile and a discontinuous one. In the discontinuous case, the period of the inertial dip is directly related to the eigenmode condition for the inertial dip ($s_{\mathrm{core}}^{\ast }=11.32$ in the case we used). The sole measure of the period of the dip would then directly give an estimate of the convective core rotation rate, independently of a measure of the near-core rotation rate. However, when the Brunt-Väisälä profile is continuous at the core-near-core interface, the period of the inertial dip not only depends on $s_{\mathrm{core}}^{\ast }$ but also on the Γ parameter controlled by the sharpness of the Brunt-Väisälä gradient, and on level of the core to near-core differentiality through the function G. Therefore, an accurate measure of the convective core rotation depends on a reliable estimate of the stellar stratification and the near-core rotation rate. Those two treatments can be considered as two extreme cases, and stars to be studied will be attributed one or the other formalism depending on the considered mode: the JWKB treatment made for the continuous case in Section 3.4 is only valid if the radial wavelength of the gravito-inertial modes is much shorter than the scale height related to the Brunt-Väisälä frequency, and the discontinuous Brunt-Väisälä profile can be attributed to modes having wavelength much longer than this scale-height.

For the shape of the Brunt-Väisälä profile to be constrained simultaneously as the differential rotation between the core and the envelope, we saw that one has to characterize the shape of the inertial dip in the co-rotating frame. The near-core rotation has first to be constrained to perform such a measurement. The method currently used to constrain near-core rotation is based on the study of the slope of period-spacing patterns in the inertial frame. When a clear dip is seen, it can be artificially removed from the period-spacing pattern to get a slope measurement. However, the impact of the dip spans several radial orders outside of what observers consider a dip (see for instance Fig. D.2). When no clear dip is detected, for instance in the regime of high coupling parameters, or core rotation inferior to the one of the envelope, the fitting of the slope in the period-spacing pattern does not include any information on the dip structure, even though it can be present but faint. Both removing coupled modes and ignoring the presence of a dip can lead to an error in estimating the near-core rotation rate, and therefore to errors in the convective core rotation rate estimate. In any case, especially when mechanisms of potential dip structure suppression will be known, it is of interest to integrate a dip model into the period-spacing pattern fitting routine to increase the precision of near-core rotation measurement.

4.6. Effect of core density stratification along the evolution and need for forward modeling

Along the evolution, the density profile gets farther from the approximation of uniform density in the core (see Fig. 5 of Ouazzani et al. 2020). As highlighted by Wu (2005), this shifts the eigenfrequency of the pure inertial mode. Their work showed that with a density profile modeled by a power-law:

$$\begin{array}{c}\hfill \rho \propto {(1-{\left(\frac{r}{R_{\text{core}}}\right)}^2)}^{\beta },\end{array}$$(64)

with r the radial coordinate within the core, an increasing parameter β would lower the spin parameter $s_{\text{core}}^{\ast }$, shifting the dip towards low periods. This effect is backed numerically by both Ouazzani et al. (2020) and Saio et al. (2021): the mode studied in this work (l = 3; m = −1), for which the spin parameter in a core of uniform density is at $s_{\text{core}}^{\ast }=11.32$ will exhibit a lower spin parameter with density stratification. As Ouazzani et al. (2020) previously pointed out, this effect is strongest for models close to the TAMS, reaching a value of 8.83 for their most evolved model (Table 6), explaining the range of [8, 11] for dips found in Saio et al. (2021) with limited differential rotation.

To approach this problem semi-analytically, one has to solve the eigenfrequency problem of pure inertial modes in the presence of density stratification with the power-law profile or approximate the profile as described in Section 2.3 of Wu (2005). An updated function $F_l^m$ can be derived and the full impact of density stratification estimated. In the case of more complex profiles, or different core-to-envelope mixing prescriptions such as convective penetration or overshooting, a numerical study is necessary (Galoy et al. 2024). It is worth noting that as highlighted by Tokuno & Takata (2022) and seen in Saio et al. (2021), an overshooting zone will result in a multiplicity of dips, which confirms the high potential of inertial dips to probe core to envelope processes happening in deep stellar interiors.

Nevertheless, especially in the regime of high coupling parameters, and therefore shallow dips, forward modelling is necessary to be able to disentangle the effect of differential rotation from the effect of density stratification. But the present study, as well as a semi-analytical one on density stratification, holds promising information on the prediction of the location of the dips and thus an easier detection of this feature from gravity-mode period spacings.

4.7. Limitations of the scope

In this section, we discuss the impact of choices and assumptions we made on the results and evaluate the accuracy of the assumptions taken.

4.7.1. Choice of the modes

We took into account only one pure inertial mode in the analysis, mainly analysing the interaction of the (l = 3, m = −1) pure inertial mode with the (k = 0, m = −1) Kelvin g-modes. Another significant interaction exists with a non-negligible geometric factor: (l = 5, m = −1). Yet, this interaction occurs at the spin parameter s_core^* = 29.33. Even with an effect of density stratification, and differential rotation, it is very unlikely that this interaction results in a dip measurable in the range of observed period-spacing patterns. This interaction at a high spin parameter in turn justifies the assumption of the sparsity of the pure inertial modes spectra in period, which was used to derive the coupling equation.

4.7.2. Impact of the LTE eigenvalue

We have neglected the variation of the LTE horizontal eigenvalue ${\mathrm{\Lambda }}_k^m$ with the spin parameter s_env in the analysis in Section 4.2, which would lead to a non-constant buoyancy travel time Π₀ and an additional curvature in the period-spacing pattern otherwise. This is motivated by the fact that asymptotically, the eigenvalue of the LTE for Kelvin modes is constant. In the range that we analysed, from 4.25 (α_rot = 0.8) to 16.59, (α_rot = 1.05), the eigenvalue ${\mathrm{\Lambda }}_0^{-1}$ of the LTE varies for the mode (k = 0, m = −1) from 1.12 to 1.03, leading to a variation of the baseline of 4%. If the variation of the eigenvalue has to be taken into account, for instance, if the differential rotation considered is more important, or modes of different azimuthal orders have to be compared, the formalism here developed holds if we take a stretched period $\sqrt{{\mathrm{\Lambda }}_k^m}\mathrm{P}$ in the co-rotating frame. Indeed, modes in the asymptotic theory are equally spaced in this stretched period: $\sqrt{{\mathrm{\Lambda }}_k^m}\mathrm{P}(n)={\mathrm{\Pi }}_0(n+{\alpha }_g)$, with n the radial order (e.g. Christophe et al. 2018.) On the other hand, more dependencies have to be taken into account in the derivation of the period spacing, such as the decreased radial coordinate of the upper turning point r_b with increasing frequency, resulting in a variation of Π₀, or modulations due to glitches in the Brunt-Väisälä frequency. These sources of potential deviation from constant baselines in period-spacing patterns are listed in Appendix A of Tokuno & Takata (2022) and would require a systematic study on age, and mixing prescriptions on top of the sole influence of differential rotation described in this work. We nevertheless highlight that due to the very low frequency of the analysed sub-inertial Kelvin g-mode series, the latter effect of variation of r_b with frequency will likely be unimportant as for this regime N < L and the decrease of N close to the upper boundary is very sharp.

4.7.3. Underlying assumptions in the coupling mechanism

The reasoning held in Section 3.2 allowing us to simplify an infinite matrix problem into a simple coupling equation involving only one interacting mode relies on both the low value of ϵ and the non-negligible value of the geometric parameter c_k, l. Regarding the coupling parameter ϵ, it remains unclear as to why stars with dips do not show particularly low values of coupling parameter compared to stars with no detected inertial dips, as highlighted by Aerts & Mathis (2023). As already argued in Section 4.2, differential rotation leads to the formation of a clearer dip with high values of ϵ, compared to the solid-body rotating case. However, a clear and definite picture will only be obtained with more realistic situations on the rotation profile or the density stratification, accessible with numerical calculations. The geometric factor, regarded as the main parameter controlling the interaction in Ouazzani et al. (2020), was not found to vary appreciably for the mode interaction considered (see Appendix F). Still, its dependence on differential rotation holds valuable information for newer theoretical developments, such as the one developed in Galoy et al. (2024). This work derived an improvement of the Tokuno & Takata (2022) model, taking into account more modes from both sides of the boundary, rather than only considering the dominant ones, as done in this present work and Tokuno & Takata (2022)’s one. The improved model shows that both the Γ parameter and the dip’s central spin parameter depend on the geometrical factors related to the dominant interactions from both sides of the boundary. As the projection of Hough functions on Legendre polynomials also varies with the differential rotation rate considered, it is also yet unclear if the number of interactions considered for a solid-body rotating star in Galoy et al. (2024) will stay identical when considering differential rotation. A combination of the model derived in this work, accounting for differential rotation, and the model used in Galoy et al. (2024), taking into account the variation of the geometrical factors with differential rotation holds promising potential for an even more accurate modelling of the dip, and thus a more accurate detection of core and near-core properties by the dip study.

4.7.4. Validity of the traditional approximation of rotation

The traditional approximation is used in this work in the radiative zone, allowing us to stay analytical throughout the present study. This solution of the gravito-inertial modes under this approximation, as described in Prat et al. (2016), can be seen as the limit of the full system of calculations in the regime N/σ_env ≫ 1, or N/(2Ω_env)≫1 in the sub-inertial regime we are interested in. As highlighted by Gerkema et al. (2008), even for values such as N = 2Ω, the solutions under the TAR match reasonably with the full calculations. In the case of the radiative zone of the γ Doradus stars, the high stratification of molecular weight ensures a high value of the Brunt-Väisälä frequency near the core. We thus expect calculations outside of the TAR not to diverge significantly from the results of the present analysis. However, the picture is blurrier with different core-envelope mixing prescriptions. If one is to consider a core-to-envelope mixing mechanism taking into account a non-radiative near-core gradient, the step overshoot with an adiabatic gradient being the most extreme case, the traditional approximation in the envelope will no longer hold in this case, and the matching of the solutions will have to be revisited with (1) an extended propagating zone of pure inertial modes and (2) an intermediate zone in which the gravito-inertial modes structure will significantly diverge from the one given by the TAR. We hereby point out that the study of the morphology of the inertial dip can be considered also as a probe of the core-to-envelope mixing mechanisms, processes that are otherwise difficult to constrain accurately. Analyses of the potentiality of gravity modes to probe the convective core boundary and different mixing prescriptions were made in Pedersen et al. (2018) and Michielsen et al. (2019). In this regard, complete numerical calculations of gravito-inertial modes (such as the ones pursued in Dintrans & Rieutord 2000; Ballot et al. 2012; Galoy et al. 2024) are of a prime importance.

5. Conclusion and perspectives

We studied the coupling between gravito-inertial modes in the radiative envelope and pure inertial modes in the convective core of a γ Doradus star in the sub-inertial regime, allowing for a two-zones differential rotation from the two sides of the core–envelope boundary. We provide an analytical model for the interaction, and discuss the possibility of constraining core to near-core differential rotation from the study of inertial dips in the period-spacing pattern of fast-rotating γ Doradus stars.

We demonstrate that core to near-core differential rotation has an effect on the location in period of the inertial dip, and also influences the depth and the width of this characteristic feature. With envelope rotation taken as a constant, increasing core rotation will shift the inertial dip to low periods, deepening and tightening it. On the contrary, decreasing core rotation will shift the inertial dip to high periods, widening it and making it more shallow. The wider the dip, the higher the number of modes affected by the coupling. We also retrieve the influence of the coupling parameter described in Tokuno & Takata (2022), making the dip shallower for a decreasing Brunt-Väisälä near-core gradient. We also highlight the influence of considering two different models, differing on their hypothesis on the Brunt-Väisälä frequency, being either continuous or discontinuous near the boundary.

Taking realistic values of the coupling parameter from the recent study of Aerts & Mathis (2023), and differential rotation rates from Saio et al. (2021), we then discuss the detectability of core to near-core differential rotation from the study of the inertial dip. We isolate two dip parameter regions in which we find differential rotation to be detectable and demonstrate that both core-to-envelope positive and negative differential rotation can be detectable in asteroseismic data from the Kepler mission, even for core rotation differing from that of the envelope at the 1% level.

As a result of our study, the rotation rate of the convective core alone can be inferred from the dip location, with an increasing uncertainty for decreasing near-core stratification gradient. In turn, we foresee that taking into account a dip model in the period-spacing pattern will also lower the uncertainty on the measurement of the near-core rotation rate. Both measurements can be enriched by an accurate modelling of the inertial dip.

We isolate two main caveats preventing detections from this model in real stars: the unconstrained shape of the Brunt-Väisälä near-core profile and the existence of a density gradient in the convective core. For the former, we find that the impact of the choice of one of the two formalisms for the Brunt-Väisälä profile has an appreciable impact on relatively high coupling parameters, introducing an additional uncertainty on the core-to-near-core differential rotation. For the latter, the existence of a density gradient is known to lower the period of the inertial dip from the constant density picture, as shown in Ouazzani et al. (2020). The effects of a core rotating faster than the envelope and a gradient of the density of the core are therefore degenerate. Further efforts must be made to lift these degeneracies, for example with an investigation into the stratification gradient in the core and near it for the most up-to-date stellar models, taking into account the impact of age and core-to-envelope mixing prescriptions. Such a study will lower the uncertainties related to the measure of the differential rotation throughout evolution.

We decided to remain within analytical framework for the present study, by taking a bi-layer profile for the stellar rotation. This allowed us to isolate the effect of differential rotation and demonstrate its detectability. Going further, it would be of great interest to investigate a smoother rotation gradient in both the convective core and the radiative zone. Properties of inertial modes propagating in a differentially rotating shell, either with a shellular or cylindrical profile, have been proven to differ from the solid-body rotating case: the spectrum, ray paths, and damping are altered and waves tend to form attractors, both in the case of pure inertial (Baruteau & Rieutord 2013; Guenel et al. 2016) and gravito-inertial modes (Mirouh et al. 2016). Given the increased technical complexity of the mode computation, the use of a spectral code to investigate the coupling problem would be required for such a refinement (Ouazzani et al. 2020; Galoy et al. 2024).

The model hereby presented remains in the purely hydrodynamic framework. Nevertheless, internal stellar magnetism is ubiquitous in stars and can potentially be intense in convective zones and near-core regions in MS stars due to the dynamo processes taking place in the core, which could be strengthened by a potential fossil field in the radiative zone arising from the past evolution of the stars (Brun et al. 2005; Featherstone et al. 2009; Augustson et al. 2016). Moreover, stars displaying no evidence of differential rotation as derived from the formalism developed in this work could be potential candidates for strong internal magnetism, as the dynamo field can contribute to a locked rotation between the convective core and the radiative zone. It would therefore be of prime interest to revisit the coupling problem in a magneto-hydrodynamical framework, derive the potential effects on the inertial dip formation picture and the degeneracies between differential rotation, and to explore the effect of density stratification in the core. This will be the scope of a follow-up paper.

---

Acknowledgments

We thank the referee for very constructive and detailed comments that led to an improvement of the quality of our study. L.B. and L.B. gratefully acknowledge support from the European Research Council (ERC) under the Horizon Europe programme (Calcifer; Starting Grant agreement N°101165631). S. Mathis acknowledges support from the PLATO CNES grant at CEA/DAp and from the European Research Council through HORIZON ERC SyG Grant 4D-STAR 101071505. While partially funded by the European Union, views and opinions expressed are however those of the author only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. L. Barrault thanks the members of the asteroseismology group of the Institute of Astronomy (IvS) of KU Leuven, in particular T. Van Reeth, M. Vanrespaille, Z. Guo and C. Aerts, for their warm welcome during a work visit in Spring 2024, and very insightful input on the present study. The authors thank also the members of the Asteroseismology and Stellar Dynamics group of the Institute of Science and Technology Austria (ISTA) for very useful discussion: K. M. Smith, L. Einramhof, S. Torres and A. Cristea.

References

Appendix A: Derivation of gravito-inertial modes structure in the radiative envelope and expansion near the core

This Appendix describes the structure of the gravito-inertial modes propagating in the radiative envelope of the star, along with an expansion of the Lagrangian displacement and pressure perturbation at the inner boundary of the radiative zone.

Under the Traditional Approximation of Rotation (TAR) and the Cowling approximation, the oscillation variables are separable in the spherical coordinates (r,θ,φ), and can be expended as:

$$\begin{array}{c}\hfill {\xi }_r^{\prime }(r,\theta )={\displaystyle \sum _k}{\xi }_{r;k,m}^{\prime }(r){\mathrm{\Theta }}_k^m(\mu ;s_{\text{env}})\end{array}$$(A.1)

and

$$\begin{array}{c}\hfill p^{\prime }(r,\theta )={\displaystyle \sum _k}{p}_{k,m}^{\prime }(r){\mathrm{\Theta }}_k^m(\mu ;s_{\text{env}}),\end{array}$$(A.2)

where ${\mathrm{\Theta }}_k^m(\mu ;s_{\text{env}})$ are the Hough functions, defined as the eigenfunctions of the Laplace tidal equation (LTE). Hough functions depend both on θ through μ = cos θ and σ_env through s_env = 2Ω_env/σ_env the spin parameter in the envelope. The LTE describes the horizontal variations of physical quantities within the TAR (Lee & Saio 1997), and reads:

$$\begin{array}{c}\hfill {\mathcal{L}}_{s_{\text{env}}}[{\mathrm{\Theta }}_k^m(\mu ;s_{\text{env}})]=-{\mathrm{\Lambda }}_k^m(s_{\text{env}}){\mathrm{\Theta }}_k^m(\mu ;s_{\text{env}}),\end{array}$$(A.3)

with the Laplace tidal operator:

$$\begin{array}{cc}\hfill {\mathcal{L}}_{s_{\text{env}}}=& \frac{\mathrm{d}}{\mathrm{d}\mu }\left(\frac{1-{\mu }^2}{1-s_{\text{env}}^2{\mu }^2}\frac{\mathrm{d}}{\mathrm{d}\mu }\right)\hfill \\ \hfill & -\frac{1}{1-s_{\text{env}}^2{\mu }^2}(\frac{m^2}{1-{\mu }^2}+m{s}_{\text{env}}\frac{1+s_{\text{env}}^2{\mu }^2}{1-s_{\text{env}}^2{\mu }^2}).\hfill \end{array}$$(A.4)

The radial functions of ξ_r′ and p′ satisfy the following system:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}}{\mathrm{d}r}({\xi }_{r;k,m}^{\prime })& =-(\frac{2}{r}-\frac{1}{{\mathrm{\Gamma }}_1{H}_p}){\xi }_{r;k,m}^{\prime }\hfill \\ \hfill & +\frac{1}{\overline{\rho }{c}_S^2}(\frac{{\mathrm{\Lambda }}_k^m(s_{\text{env}})c_S^2}{r^2{\sigma }_{S,\text{env}}^2}-1)p_{k,m}^{\prime }(r),\hfill \end{array}$$(A.5)

$$\begin{array}{c}\hfill \frac{\mathrm{d}{p}_{k,m}^{\prime }(r)}{\mathrm{d}r}=\overline{\rho }({\sigma }_{\text{env}}^2-N^2){\xi }_{r;k,m}^{\prime }-\frac{p_{k,m}^{\prime }(r)}{{\mathrm{\Gamma }}_1{H}_p},\end{array}$$(A.6)

with H_p, N, Γ₁, $\overline{\rho }$, and c_S being the pressure scale height, the Brunt-Väisälä frequency, the first adiabatic index, the mean density, and the sound speed respectively. ${\mathrm{\Lambda }}_k^m(s_{\text{env}})$ is the eigenvalue of the LTE, quantized by k (introduced by Lee & Saio 1997) and m, depending on the spin parameter of the envelope s_env. For Kelvin modes, the eigenvalue reaches a constant value in the regime of high spin parameters. The asymptotic development has been derived by Townsend (2003), later extended in Townsend (2020). Tokuno & Takata (2022) derive, using a Jeffreys-Wentzel-Kramers-Brillouin (JWKB) analysis, the radial wavenumber (see also e.g. Unno et al. 1989):

$$\begin{array}{c}\hfill k_r^2=\frac{{\sigma }_{\text{env}}^2-N^2}{c_S^2}(1-\frac{{\mathrm{\Lambda }}_k^m(s_{\text{env}})c_S^2}{r^2{\sigma }_{\text{env}}^2}).\end{array}$$(A.7)

The integral over the whole radiative zone reads, in the low-frequency regime:

$$\begin{array}{c}\hfill {\int }_{r_a}^{r_b}{k}_r\mathrm{d}r\simeq {\displaystyle {\int }_{r_a}^{r_b}}\frac{N}{{\sigma }_{\mathrm{env}}}\frac{\sqrt{{\mathrm{\Lambda }}_k^m(s_{\text{env}})}}{r}\mathrm{d}r=\frac{{\pi }^2{s}_{\text{env}}}{{\mathrm{\Omega }}_{\text{env}}{\mathrm{\Pi }}_0},\end{array}$$(A.8)

r_a and r_b being respectively the lower and the upper turning points of the cavity (see Fig. 1), and Π₀ the asymptotic period spacing:

$$\begin{array}{c}\hfill {\mathrm{\Pi }}_0\simeq \frac{2{\pi }^2}{\sqrt{{\mathrm{\Lambda }}_k^m(s_{\text{env}})}}{\left({\displaystyle {\int }_{r_a}^{r_b}}\frac{N}{r}\mathrm{d}r\right)}^{-1}.\end{array}$$(A.9)

This period spacing is approximately constant for Kelvin modes, neglecting the variation of the eigenvalue of the LTE. A discussion of this approximation can be found in Appendix A of Tokuno & Takata (2022), complemented by a discussion taking into account the differentially rotating situation in our subsection 4.7.2.

A.1. The continuous case

The expressions derived using this JWKB analysis for the Eulerian pressure and the Lagrangian displacement remain valid in the range [r_a, r_b]. Yet, the matching of the inertial waves in the core and the gravito-inertial waves in the envelope occurs at the edge of the convective core, at the radial coordinate R_core. In the case of a continuous Brunt-Väisälä frequency near-core, the two radial coordinates r_a and R_core differ. The next key step of the analysis is thus to derive an approximate expression for the radial displacement and the pressure at the edge of the convective core, from the expressions of the Eulerian pressure perturbation and the Lagrangian displacement at r_a. We perform a Taylor expansion near r = r_a, given that the distance between R_core and r_a is small.

The vertical wave number can be written as

$$\begin{array}{c}\hfill k_r^2\simeq \frac{{\mathrm{\Lambda }}_k^m(s_{\text{env}})s_{\text{env}}^2}{{ϵ}^3{r}_a^3}(r-r_a),\end{array}$$(A.10)

where we have introduced the small parameter

$$\begin{array}{c}\hfill ϵ={\left(\frac{r_a}{4{\mathrm{\Omega }}_{\text{env}}^2}\frac{\text{d}{N}^2}{\text{d}r}{|}_{r_a}\right)}^{-1/3}\end{array}$$(A.11)

assuming ϵ ≪ 1 because of the steepness of the gradient of the Brunt-Väisälä frequency at the lower turning point. Expanding all quantities, Tokuno & Takata (2022) show that the radial Lagrangian displacement and the Eulerian pressure perturbation can be rewritten synthetically as such:

$$\begin{array}{c}\hfill \frac{{\xi }_{r;k,m}^{\prime }}{r}{|}_{r_a}\simeq QϵX_k^m(s_{\text{env}})\end{array}$$(A.12)

and

$$\begin{array}{c}\hfill \frac{p_{k,m}^{\prime }}{\overline{p}{\mathrm{\Gamma }}_1}{|}_{r_a}\simeq Q\frac{r_a^2{\sigma }_{\text{env}}^2}{c_S^2}{Y}_k^m(s_{\text{env}}),\end{array}$$(A.13)

where we have introduced, very similarly to Tokuno & Takata (2022), the following functions:

$$\begin{array}{c}\hfill X_k^m(s_{\text{env}})={\mathrm{\Lambda }}_k^m{(s_{\text{env}})}^{1/6}{s}_{\text{env}}^{2/3}sin(\frac{{\pi }^2{s}_{\text{env}}}{{\mathrm{\Omega }}_{\mathrm{env}}{\mathrm{\Pi }}_0}-\frac{\pi }{6})\end{array}$$(A.14)

and

$$\begin{array}{c}\hfill Y_k^m(s_{\text{env}})=\alpha {\mathrm{\Lambda }}_k^m{(s_{\text{env}})}^{-1/2}{s}_{\text{env}}^{4/3}sin(\frac{{\pi }^2{s}_{\text{env}}}{{\mathrm{\Omega }}_{\text{env}}{\mathrm{\Pi }}_0}-\frac{5\pi }{6}).\end{array}$$(A.15)

In those equations, the dimensionless parameters α and Q are defined as:

$$\begin{array}{c}\hfill \alpha \equiv \frac{3^{1/3}\mathrm{\Gamma }(2/3)}{\mathrm{\Gamma }(1/3)}\simeq 0.73,\end{array}$$(A.16)

$$\begin{array}{c}\hfill Q\equiv -\frac{c}{3^{2/3}\mathrm{\Gamma }(2/3)}\frac{{\mathrm{\Lambda }}_k^m{(s_{\text{env}})}^{1/6}}{{ϵ}^{1/2}{\mathrm{\Omega }}_{\text{env}}{r}_a^{5/2}}\frac{1}{{\overline{\rho }}^{1/2}}{|}_{r_a},\end{array}$$(A.17)

c being an arbitrary constant and Γ(x) the Gamma function.

For this expansion to be accurate, and to match the approximate expressions at the first order with the respective quantities at the edge of the convective core, we need to assert that the values of ξ_r and p′ at R_core are close to the ones at r_a. From a first order expansion of N² from r_a to R_core, we obtain taking Eq. (A.11):

$$\begin{array}{c}\hfill \frac{R_{\text{core}}}{r_{\text{a}}}\sim 1-\frac{{ϵ}^3}{s_{\text{env}}^2}.\end{array}$$(A.18)

We thus highlight here that the formalism of Tokuno & Takata (2022) only holds for small values of ϵ and high values of the spin parameter. This limitation will be further discussed in subsection 4.7. The expressions of the pressure perturbation and the radial Lagrangian displacement are, in the inertial frame:

$$\begin{array}{c}\hfill \frac{{\xi }_r^{\prime }}{r}{|}_{R_{\text{core}}}={\displaystyle \sum _k}{a}_kϵX_k^m(s_{\text{env}}){\mathrm{\Theta }}_k^m(\mu ;s_{\text{env}})\end{array}$$(A.19)

and

$$\begin{array}{c}\hfill \frac{p^{\prime }}{\overline{p}{\mathrm{\Gamma }}_1}{|}_{R_{\text{core}}}=\frac{{\sigma }_{\text{env}}^2{R}_{\text{core}}^2}{c_S^2}{\displaystyle \sum _k}{a}_k{Y}_k^m(s_{\text{env}}){\mathrm{\Theta }}_k^m(\mu ;s_{\text{env}}).\end{array}$$(A.20)

A.2. The discontinous case

The vertical wave number can be written at the upper edge of the boundary as

$$\begin{array}{c}\hfill k_r^2=\frac{{\mathrm{\Lambda }}_k^m(s_{\mathrm{env}})s_{\mathrm{env}}^2}{4{\stackrel{\sim }{ϵ}}^2{R}_{\mathrm{core}}^2},\end{array}$$(A.21)

with $\stackrel{\sim }{ϵ}=\frac{{\mathrm{\Omega }}_{\mathrm{env}}}{N_0}$. Very similarly to the continuous case, the terms appearing in the decomposition of the Lagrangian displacement and Eulerian pressure perturbation are:

$$\begin{array}{c}\hfill \frac{{\xi }_{r;k,m}^{\prime }}{r}{|}_{R_{\text{core}}^{+}}\simeq Qϵ{\stackrel{\sim }{X}}_k^m(s_{\text{env}})\end{array}$$(A.22)

and

$$\begin{array}{c}\hfill \frac{p_{k,m}^{\prime }}{\overline{p}{\mathrm{\Gamma }}_1}{|}_{R_{\text{core}}^{+}}\simeq Q\frac{R_{\mathrm{core}}^2{\sigma }_{\text{env}}^2}{c_S^2}{Y}_k^m(s_{\text{env}}).\end{array}$$(A.23)

The JWKB analysis described in Tokuno & Takata (2022) is providing:

$$\begin{array}{c}\hfill {\stackrel{\sim }{X}}_k^m(s_{\text{env}})=2{\mathrm{\Lambda }}_k^m{(s_{\text{env}})}^{1/4}{s}_{\text{env}}^{1/2}sin(\frac{{\pi }^2{s}_{\text{env}}}{{\mathrm{\Omega }}_{\text{env}}{\mathrm{\Pi }}_0}-\frac{\pi }{4})\end{array}$$(A.24)

and:

$$\begin{array}{c}\hfill {\stackrel{\sim }{Y}}_k^m(s_{\text{env}})=-{\mathrm{\Lambda }}_k^m{(s_{\text{env}})}^{-1/4}{s}_{\text{env}}^{3/2}cos(\frac{{\pi }^2{s}_{\text{env}}}{{\mathrm{\Omega }}_{\text{env}}{\mathrm{\Pi }}_0}-\frac{\pi }{4}).\end{array}$$(A.25)

We end up with the following expression for the Lagrangian displacement and the Eulerian pressure perturbation:

$$\begin{array}{c}\hfill \frac{{\xi }_r^{\prime }}{r}{|}_{R_{\text{core}}^{+}}={\displaystyle \sum _k}{a}_k\stackrel{\sim }{ϵ}{\stackrel{\sim }{X}}_k^m(s_{\text{env}}){\mathrm{\Theta }}_k^m(\mu ;s_{\text{env}})\end{array}$$(A.26)

and:

$$\begin{array}{cc}\hfill \frac{p^{\prime }}{{\mathrm{\Gamma }}_1\overline{p}}{|}_{R_{\text{core}}^{+}}& =\frac{{\sigma }_{\text{env}}^2{R}_{\text{core}}^2}{c_S^2}\hfill \\ \hfill & \times {\displaystyle \sum _k}{a}_k{\stackrel{\sim }{Y}}_k^m(s_{\text{env}}){\mathrm{\Theta }}_k^m(\mu ;s_{\text{env}})\hfill \end{array}$$(A.27)

Appendix B: Deriving Bryan solutions for pure inertial modes in the convective core

We hereby derive the structure of pure inertial modes propagating in the convective core of the star and give an expression for the Lagrangian displacement and pressure perturbation at the outer boundary of the convective core.

In the framework of solid-body rotation and uniform averaged density in the core, under the Cowling approximation, using the scalar variable Ψ defined as:

$$\begin{array}{c}\hfill \mathrm{\Psi }=\frac{1}{{\sigma }_{\text{core}}^2}\frac{p^{\prime }}{\overline{\rho }},\end{array}$$(B.1)

the momentum and mass conservation equations reduce to the Poincaré equation:

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2\mathrm{\Psi }-s_{\text{core}}^2\frac{{\partial }^2\mathrm{\Psi }}{\partial z^2}=0,\end{array}$$(B.2)

z being the axial coordinate in the cylindrical coordinate system, the cylindrical axis being along the rotation vector.

The radial Lagrangian displacement vector projected on ${\widehat{\mathit{e}}}_r$ in the formalism of prograde modes having m< 0 is:

$$\begin{array}{c}\hfill {\xi }_r^{\prime }=\frac{1}{1-s_{\text{core}}^2}[\frac{\partial }{\partial r}+\frac{m{s}_{\text{core}}}{r}-\mu s_{\text{core}}^2\frac{\partial }{\partial z}]\mathrm{\Psi }\end{array}$$(B.3)

Eq. (B.1) can also be written as:

$$\begin{array}{c}\hfill p^{\prime }=\overline{p}{\mathrm{\Gamma }}_1\frac{{\sigma }_{\text{core}}^2}{c_S^2}\mathrm{\Psi },\end{array}$$(B.4)

where we have used the sound speed velocity definition: $c_S^2={\mathrm{\Gamma }}_1\overline{p}/\overline{\rho }$. Using the ellipsoidal coordinates, Wu (2005) shows that the expressions of the relevant quantities near the core are:

$$\begin{array}{c}\hfill \frac{{\xi }_{r;l,m}^{\prime }(r)}{r}{|}_{R_{\text{core}}^{-}}\propto C_l^m(1/s_{\text{core}})P_l^m(\mu ),\end{array}$$(B.5)

$$\begin{array}{c}\hfill \frac{p_{l,m}^{\prime }(r)}{\overline{p}{\mathrm{\Gamma }}_1}{|}_{R_{\text{core}}^{-}}\propto \frac{R_{\text{core}}^2{\sigma }_{\text{core}}^2}{c_S^2}{P}_l^m(1/s_{\text{core}})P_l^m(\mu )\end{array}$$(B.6)

with $C_l^m$ being as defined in Ouazzani et al. (2020):

$$\begin{array}{c}\hfill C_l^m(x)=x(\frac{\mathrm{d}{P}_l^m(x)}{\mathrm{d}x}-\frac{m}{1-x^2}{P}_l^m(x)).\end{array}$$(B.7)

These analytical solutions are well-known as Bryan solutions, having been first described by Bryan (1889).

The fixed boundary condition ξ′(r) = 0 at the edge of the core translates into $C_l^m(1/s_{\text{core}})=0$, which is an eigenvalue problem setting the spin parameters of the pure inertial modes in the core. Table 1 of Ouazzani et al. (2020) reports the spin parameters of the prograde m = −1 modes in the core. Finally, the solution of the variables of interest near the core is the linear combination of such solutions, in the inertial frame taking into account the Doppler shift:

$$\begin{array}{c}\hfill \frac{{\xi }_r^{\prime }}{r}{|}_{R_{\text{core}}^{-}}={\displaystyle \sum _l}{b}_l{C}_l^m(1/s_{\text{core}}){\stackrel{\sim }{P}}_l^m(\mu ),\end{array}$$(B.8)

and

$$\begin{array}{c}\hfill \frac{p^{\prime }}{\overline{p}{\mathrm{\Gamma }}_1}{|}_{R_{\text{core}}^{-}}=\frac{{\sigma }_{\text{core}}^2{R}_{\text{core}}^2}{c_S^2}{\displaystyle \sum _l}{b}_l{P}_l^m(1/s_{\text{core}}){\stackrel{\sim }{P}}_l^m(\mu ).\end{array}$$(B.9)

where b_l are constant factors, and ${\stackrel{\sim }{P}}_l^m$ are the normalized Legendre polynomials:

$$\begin{array}{c}\hfill {\stackrel{\sim }{P}}_l^m\equiv \sqrt{\frac{(2l+1)(l-m)!}{2(l+m)!}}{P}_l^m(x).\end{array}$$(B.10)

Appendix C: Derivation of the modified Lorentzian profile

To compute the period spacing expression, we write Eq. (35) for two neighbouring solutions s₁ and s₂ with G⁻¹(s₂) > G⁻¹(s₁). We assume that the frequency spectrum of the gravito-inertial modes is dense. The two solutions being on two different branches of the cotangent, and the function 1/(G⁻¹(s)−G⁻¹(s)^*) not varying appreciably due to the density of the spectrum, we have $\frac{{\pi }^2(s_2-s_1)}{{\mathrm{\Omega }}_{\text{env}}{\mathrm{\Pi }}_0}-\pi =2x\ll 1$. In the case of $G^{-1}(s_1)<s_{\text{core}}^{\ast }<G^{-1}(s_2)$ the solutions belong to the same branch, as pointed out by Tokuno & Takata (2022). We thus rewrite Eq. (35) as:

$$\begin{array}{c}\hfill tan\left(\frac{{\pi }^{2}s}{{\mathrm{\Omega }}_{\text{env}}{\mathrm{\Pi }}_0}\right)\simeq -\frac{1}{\frac{1}{\sqrt{3}}+\frac{ϵ/V_{\text{diff}}}{G^{-1}(s)-G^{-1}(s^{\ast })}}\end{array}$$(C.1)

and the solutions are on two branches of the tangent. Thus, in both cases:

$$\begin{array}{c}\hfill \frac{{\pi }^2{s}_{1,2}}{{\mathrm{\Omega }}_{\text{env}}{\mathrm{\Pi }}_0}-\frac{\pi }{6}=\frac{{\pi }^2\overline{s}}{{\mathrm{\Omega }}_{\text{env}}{\mathrm{\Pi }}_0}-\frac{\pi }{6}\pm \frac{\pi }{2}\mp x\end{array}$$(C.2)

with $\overline{s}=\frac{s_1+s_2}{2}$. The resulting system is thus, noting $\overline{S}=\frac{{\pi }^2\overline{s}}{{\mathrm{\Omega }}_{\text{env}}{\mathrm{\Pi }}_0}-\frac{\pi }{6}$ and using a first order Taylor expansion:

$$\begin{array}{c}\hfill \{\begin{array}{cc}cot(\overline{S})-x(1+{cot}^2(\overline{S}))+1/\sqrt{3}\hfill & \simeq -\frac{ϵ/V_{\text{diff}}}{G^{-1}(s_1)-G^{-1}{(s)}^{\ast }}\hfill \\ cot(\overline{S})+x(1+{cot}^2(\overline{S}))+1/\sqrt{3}\hfill & \simeq -\frac{ϵ/V_{\text{diff}}}{G^{-1}(s_2)-G^{-1}{(s)}^{\ast }}\hfill \end{array}.\end{array}$$(C.3)

We then take the sum and the difference. Considering the right-hand side terms, with $\overline{G^{-1}(s)}=(G^{-1}(s_1)+G^{-1}(s_2))/2$:

$$\begin{array}{cc}& \frac{1}{G^{-1}(s_2)-G^{-1}{(s)}^{\ast }}+\frac{1}{G^{-1}(s_1)-G^{-1}{(s)}^{\ast }}\hfill \\ \hfill & =\frac{2(\overline{G^{-1}(s)}-G^{-1}{(s)}^{\ast })}{{(\overline{G^{-1}(s)}-G^{-1}{(s)}^{\ast })}^2-{(\frac{G^{-1}(s_2)-G^{-1}(s_1)}{2})}^2}\hfill \end{array}$$(C.4)

and

$$\begin{array}{cc}& \frac{1}{G^{-1}(s_2)-G^{-1}{(s)}^{\ast }}-\frac{1}{G^{-1}(s_1)-G^{-1}{(s)}^{\ast }}\hfill \\ \hfill & =\frac{G^{-1}(s_2)-G^{-1}(s_1)}{{(\overline{G^{-1}(s)}-G^{-1}{(s)}^{\ast })}^2-{(\frac{G^{-1}(s_2)-G^{-1}(s_1)}{2})}^2}.\hfill \end{array}$$(C.5)

We approximate the denominator considering: $\overline{G^{-1}(s)}-G^{-1}{(s)}^{\ast }\gg \frac{G^{-1}(s_2)-G^{-1}(s_1)}{2}$. This approximation is accurate moving away from the location of the pure inertial mode. With such approximation, we retrieve:

$$\begin{array}{cc}& [1+(\frac{ϵ/V_{\text{diff}}}{\overline{G^{-1}(s)}-G^{-1}{(s)}^{\ast }}+\frac{1}{\sqrt{3}}{)}^2]\hfill \\ \hfill & \times [\frac{{\pi }^2(s_1-s_2)}{{\mathrm{\Omega }}_{\text{env}}{\mathrm{\Pi }}_0}-\pi ]\simeq -\frac{ϵ}{V_{\text{diff}}}\frac{G^{-1}(s_1)-G^{-1}(s_2)}{(\overline{G^{-1}(s)})-G^{-1}{(s)}^{\ast }{)}^2}.\hfill \end{array}$$(C.6)

We further manipulate this equation to give an analytical insight. Provided that $s_1-s_2\ll \overline{s}$, and defining $s_{\text{env}}^{\ast }=G(s_{\text{core}}^{\ast })=G(G^{-1}{(s)}^{\ast })$, we rewrite the equation as such:

$$\begin{array}{cc}& [1+(\frac{ϵ/V_{\text{diff}}}{G^{-1}(\overline{s})-G^{-1}(s^{\ast })}+\frac{1}{\sqrt{3}}{)}^2]\hfill \\ \hfill & \times [\frac{{\pi }^2(s_1-s_2)}{{\mathrm{\Omega }}_{\text{env}}{\mathrm{\Pi }}_0}-\pi ]\simeq -\frac{ϵ}{V_{\text{diff}}}\frac{G^{-1}(s_1)-G^{-1}(s_2)}{{(G^{-1}(\overline{s})-G^{-1}(s^{\ast }))}^2}.\hfill \end{array}$$(C.7)

We then perform a first order expansion of G⁻¹, both around $\overline{s}$ and $\frac{\overline{s}+s^{\ast }}{2}$. The smoothness of G⁻¹ ensures that the latter expansion remains accurate near the dip. This gives Eq. (36). Reasoning on Eq. (C.1), the expansion holds the same result.

Appendix D: Building perturbed period-spacing patterns

We investigate the detectability of the differential rotation from the inertial dip analysis, in real data, by considering a typical noise from the observations.

The question of the uncertainty to attach to the period-spacing patterns has been extensively analysed by Bowman & Michielsen (2021) in the case of a well-studied star, KIC 7760680, a slowly pulsating B star. It is highly non-trivial and has potential sources of uncertainty that can come from (1) the light-curve used, (2) the iterative pre-whitening process, (3) the building of the period-spacing pattern, which has still a part of subjectivity from the scientist. The excitation and damping mechanisms of g-modes are still unsufficiently known and the lowest the signal of a peak, and the highest the uncertainty on the period obtained. Thus, no accurate model of noise can be derived from first principles, such as the finite observation times of space missions. We take our noise model from Kjeldsen (2003) and Montgomery & O’Donoghue (1999), which inferred the following structure for the uncertainty on the frequency:

$$\begin{array}{c}\hfill \delta f=0.44\frac{⟨a⟩}{a}\frac{1}{\mathrm{\Delta }T},\end{array}$$(D.1)

where ⟨a⟩ is the noise level in the amplitude spectrum, a is the amplitude of peak, so that a/⟨a⟩ is the signal over noise ratio. ΔT is the observing time of the considered space mission. In our case, we take ΔT to be of 4 years. Without prior information on the signal-over-noise ratio for each peak, the ratio ⟨a⟩/a is not possible to compute. We then take a hypothesis of S/N = 4.6 for every signal composing the period-spacing pattern. This value of signal over noise ratio was used in the study of Bowman & Michielsen (2021) as a stopping criterion for their analysis. It is worth noting that as mentioned by Van Reeth et al. (2015), the value of S/N to be taken is controversial and a priori distinct for every star. Li et al. (2019) for instance used a S/N stopping criterion of 3. Converting this relation in the period domain, with the adopted S/N:

$$\begin{array}{c}\hfill \delta {\mathrm{P}}_{\mathrm{in}}=\frac{{\mathrm{P}}_{\mathrm{in}}^2}{10.45\mathrm{\Delta }T}.\end{array}$$(D.2)

Then remains the issue of the extent of the generated period-spacing pattern in terms of accessible spin parameters. For this, we analyse the period-spacing patterns with a confirmed analysis by Saio et al. (2021). The results are shown in Fig. D.1. First, the minimum, mean and maximum amount of modes in the analysis of Saio et al. (2021) are respectively 21, 33, and 50. We thus take these values as three distinct cases for our detectability analysis.

Second, we aim to build a representative period-spacing pattern in term of the range of periods of the gravito-inertial modes. Sticking to the case of KIC12066947 for which the extremal values of periods in the inertial frame are P_in, min = 0.33 d and P_in, max = 0.405 d, corresponding to values P_co, min and P_co, max in the co-rotating frame, we take n/2 modes around the mean value (P_co, min + P_co, max)/2 to build the period-spacing pattern, then convert it in the inertial frame for the MCMC analysis described in the main text to fit simultaneously the internal rotation rate Ω_env and the buoyancy travel time Π₀ as well as the parameters (α_rot and Γ specific to the inertial dip. Two examples of mock period-spacing patterns are shown in Fig. D.2 in the co-rotating frame for a better vizualisation of the dip structure.

Appendix E: Bayesian fitting

This appendix presents the cornerplots issued following the procedure described in subsection 4.3.2, using the emcee package (Foreman-Mackey et al. 2013), for the particular situation of input parameters (Ω_env/2π, Π₀, P₀, α_rot,Γ) = (2.16 c.d⁻¹, 4175 s, 1.419 d, 0.98, 12 hours).

Appendix F: Variation of the geometric factor with the amount of differential rotation

The geometric factor c_k, l measures the geometrical overlap between Hough functions ${\mathrm{\Theta }}_k^m(\mu ,s_{\text{env}})$ and Bryan solutions, represented by a Legendre polynomial $P_l^m(\mu )$, as defined in Eq. (18). The two definitions given both in Tokuno & Takata (2022) and in this present work on the one hand, and in Ouazzani et al. (2020) on the other hand, are consistent if both the Hough functions and the Legendre polynomials are duly normalized. For the reasoning held in subsection 3.2 to hold, this geometric factor needs not to be negligibly small. We thus compute the geometric factors corresponding to the different s_env^* at which the interaction occurs with varying differential rotation rate α_rot. Hough functions are computed using a code described in Appendix A of Prat et al. (2019), using a Chebyshev collocation scheme described in Wang et al. (2016).

We find that at the differential rotation rates taken into account in this work, the geometric factor does not reach negligibly small values. With increasing differential rotation rates, the geometric factor is gradually increasing. This comes from the fact that Hough functions, which reduce to a simple Legendre polynomial in the non-rotating case, progressively diverge from their non-rotating counterparts, being composed of further harmonics in their Legendre polynomial decomposition. In the case of the considered interaction, the spin parameter at which the interaction occurs is progressively shifted towards low values with increasing core rotation over constant envelope rotation. The geometric factor thus decreases in this regime.

Appendix G: Impact of modulations due to chemical or thermal stratification on the inertial dip detectability

The impact of buoyancy glitches is enscribed in period-spacing patterns through modulations, whose shape depend on the particular hypothesized glitch profile (Miglio et al. 2008; Cunha et al. 2015, 2019, 2024). The location along the stellar radius, number and shape of those glitches depends on the age and particular choices of stellar modelling parameters, which renders the sole analysis of chemical or thermal stratification challenging in itself, and often relying on forward modelling (see e.g. Mombarg et al. 2022). We thus chose not to consider these modulations in our fitting analysis, and to assume a knowledge of the baseline only limited by the intrinsic observational errors. We aim to discuss this approach in this Appendix. Rather than taking an arbitrary glitch location and shape in this analysis, we aim to base ourselves on observed modulations in real data. We take the case of KIC12066947, for which a periodic modulation is observed at low periods in the period-spacing pattern (Fig. G.1), as well as a clear inertial dip (highlighted with a red dashed inverse Lorentzian profile).

We introduce this periodic modulation in custom period-spacing patterns reproducing the behaviour of the baseline shown in Fig. G.1, with a same number of modes, while introducing an inertial dip whose amplitude is made varying, allowing us to study the impact of oscillatory features on the dip fitting. We apply our MCMC analysis on two different period-spacing patterns exhibiting (1) an inertial dip of a significant amplitude compared to the chemical modulation (Bottom panel of Fig. G.2, green points) (2) a shallow inertial dip of a similar amplitude compared to the modulation (Top panel of Fig. G.2, green points).

By its treatment of the sole errors related to the intrinsic finite observing time, the likelihood described in Eq. (59) would not be fitted for this analysis, as it gives increasing weight to low-period modes, which are primarily perturbed in this present analysis. We thus investigate the results given by three different implementations of the likelihood. Alongside with the likelihood used in the main text (case 1), we adapt our likelihood to take into account a free parameter in our error estimations (case 2). This implementation is related to the one of Aerts et al. (2018), without taking the Mahalanobis distance in our simplified case. The log-likelihood takes the following form:

$$\begin{array}{c}\hfill log\mathcal{L}(\mathcal{F})=-\frac{1}{2}{\displaystyle \sum _k}[\frac{{({\mathrm{P}}_{k,\mathrm{in}}^{\mathrm{mod}}-{\mathrm{P}}_{k,\mathrm{in}}^{\mathrm{pert}})}^2}{s_k^2}+ln(2\pi s_k^2)],\end{array}$$(G.1)

with $s_k^2={\sigma }_k^2+f^2{({\mathrm{P}}_{k,\mathrm{in}}^{\mathrm{mod}})}^2$, ${\sigma }_k^2$ being the error previously taken for the likelihood related to intrinsic observationnal uncertainty and f an additional free parameter. This choice makes the model flexible to modulations not comprised in the observationnal error. We also consider fitting period spacings instead of period with a similar treatment of the error (case 3) for completeness, as results can differ taking one or the other formulation (see e.g. Cunha et al. 2024).

Results are shown in Fig. G.2, for each of the considered cases. The precise estimated distribution of the parameters is given in Table G.1 for the shallowest dip. We find that in the first configuration, the parameters related to the dip (α_rot, Γ) are retrieved accurately irrespective of the likelihood used. This in turn justifies our use of the first likelihood for the dips studied in the main text, without modulation. The rotation rate is constrained as well, but Π₀, the baseline is overestimated. In the second configuration, a difference emerges between the values retrieved using the different likelihoods. We notice our method used so far (case 1, blue points) to be biased towards high core-to-envelope differential rotation in this configuration, as we were only considering the intrinsic observationnal error in the calculations of the likelihood. Modes of low periods being less affected by this noise, they take a significant weight in the obtained result. We find that fitting period spacings (case 3) provides the closest estimation of dip parameters (α_rot, Γ), all cases introducing a bias on the fitted baseline as well in this case. We nevertheless highlight first that this case shows a relative robustness of the differential rotation (α_rot) and near-core sratification (Γ) detection even in the hypothesis of a shallow dip relative to the chemical modulations. Second, this fitting routine does not comprise any information on the chemical modulations and a more educated fitting method would show better constraints on both the inertial dips and the chemical modulations. A dedicated data-oriented study would be necessary in this case but is outside of the scope of our present proof-of-concept study.

All Tables

All Figures