Issue 
A&A
Volume 667, November 2022



Article Number  A43  
Number of page(s)  21  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202244247  
Published online  04 November 2022 
Stochastic excitation of internal gravity waves in rotating late Ftype stars: A 3D simulation approach
^{1}
Université ParisCité, Université ParisSaclay, CEA, CNRS, AIM, 91191 GifsurYvette, France
email: sylvain.breton@cea.fr
^{2}
Université ParisSaclay, Université ParisCité, CEA, CNRS, AIM, 91191 GifsurYvette, France
Received:
11
June
2022
Accepted:
22
August
2022
Context. There are no strong constraints placed thus far on the amplitude of internal gravity waves (IGWs) that are stochastically excited in the radiative interiors of solartype stars. Late Ftype stars have relatively thin convective envelopes with fast convective flows and tend to be fast rotators compared to solartype stars of later spectral types. These two elements are expected to directly impact the IGW excitation rates and properties.
Aims. We want to estimate the amplitude of stochastically excited gravity modes (gmodes) in Ftype stars for different rotational regimes.
Methods. We used the ASH code to perform 3D simulations of deepshell models of 1.3 M_{⊙} Ftype solartype stars, including the radiative interior and the shallow convective envelope.
Results. We found different differential rotation regimes in the convective zone, depending on the rotation rate we imposed on the stellar models. We find that the convective structures and the overshoot properties are affected by rotation. The IGWs are excited by interface interactions between convective plumes and the top of the radiative interior. We were able to characterise the IGWs and gmode properties in the radiative interior, and we compared these properties using the computation from the 1D oscillation code GYRE. The amplitude of lowfrequency modes is significantly higher in fastrotating models and the evolution of the period spacing of consecutive modes exhibits evidence of a behaviour that is modified by the influence of the Coriolis force. For our fastest rotating model, we were able to detect the intermediate degree gmode signature near the top of the simulation domain. Nevertheless, the predicted luminosity perturbations from individual modes still remain at small amplitudes.
Conclusions. We obtained mode amplitudes that are several orders of magnitude higher than those of prior 3D simulations of solar models. Our simulations suggest that gmode signatures could be detectable in late Ftype stars, which are the hottest mainsequence solartype pulsating stars. We therefore emphasise that they constitute object of primary importance for improving our understanding of internal stellar dynamics.
Key words: asteroseismology / stars: rotation / stars: solartype / waves / methods: numerical / hydrodynamics
© S. N. Breton et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1. Introduction
The propagation of internal gravity waves (IGW) in stellar interiors is an expected and wellknown phenomenon (see e.g. ChristensenDalsgaard 2008; Maeder 2009; Kippenhahn et al. 2012, and references therein). These waves are propagative in stably stratified resonant cavities located in radiative regions and evanescent in unstable convective regions. Standing IGWs are usually referred as gravity modes (gmodes). In the stellar interior, they may be excited by different mechanisms. In massive stars, small perturbations of the nuclear reaction rate can lead to growing instability in temperature and may therefore excite oscillation modes, an effect known as the ϵmechanism (e.g. Moravveji et al. 2012). In intermediatemass and massive stars, the κmechanism that is related to opacitybump in ionization regions is also able to drive mode excitation (e.g. Unno et al. 1989; Bowman 2020); whereas for early Ftype stars pulsating as γ Doradus pulsators, a flux blocking mechanism at the bottom of the shallow convective zone is also invoked (e.g. Guzik et al. 2000; Dupret et al. 2005). In lowmass stars, IGWs are stochastically excited by turbulent convective motions at the interface between the radiative zone and the convective zone. When considering these stochastic mechanisms, the relative importance between the large scale convective plumes penetrating the overshoot regions (e.g. Press 1981; Hurlburt et al. 1986; Garcia Lopez & Spruit 1991; Zahn 1991; Zahn et al. 1997; Browning et al. 2004; Rogers & Glatzmaier 2005; Brun et al. 2011; Alvan et al. 2014, 2015; Pinçon et al. 2016) and the Reynold stresses in the bulk region above the radiative zone (e.g. Goldreich & Kumar 1990; Belkacem et al. 2009; Samadi et al. 2010; Lecoanet & Quataert 2013) has been discussed over the years (Talon & Charbonnel 2003; Rogers et al. 2013; Perrard et al. 2013; Lecoanet et al. 2015; Augustson et al. 2020; Le Saux et al. 2022).
Because their properties are intrinsically related to the structure and dynamics of the innermost stellar regions, the question of the gmode surface amplitude in the Sun and the possibility to observe such modes has therefore been a heavily debated topic since the advent of helioseismology (see e.g. Appourchaux et al. 2013). With the introduction of spaceborne asteroseismology and the observation of hundreds of mainsequence solartype pulsating stars (e.g. Chaplin et al. 2011; Mathur et al. 2022), it is interesting to consider the case of solartype pulsators other than the Sun.
Among mainsequence solartype pulsators, Ftype stars are probably the most promising for gmode detections. In such stars, fast convective flows (e.g. Brun et al. 2017) develop inside a thin surface convective zone, which tend to favour the tunneling of gmodes towards the surface. Most of observed Ftype stars also tend to be fast rotators: in the Kepler (Borucki et al. 2010) catalog provided by Santos et al. (2021), the median surface rotation period for Ftype stars with measured rotation rate is close to 6 days. Through the action of the Coriolis force, rotation is expected to have an effect on IGW behaviour (Lee & Saio 1997; Townsend 2003b) by turning them into gravitoinertial waves, which remain propagative in an unstable layer if their frequency, ω, is below twice the local rotation frequency (Mathis et al. 2014). In recent years, the asteroseismology of fast rotators has presented spectacular developments both with regard to theoretical (e.g. Prat et al. 2019, 2020; Dhouib et al. 2021a,b, 2022; Ouazzani et al. 2017, 2019, 2020) and observational (e.g. Van Reeth et al. 2015, 2016; Pápics et al. 2017; Szewczuk et al. 2021; Li et al. 2020; Aerts et al. 2021; Pedersen et al. 2021) aspects, which has enabled probes of the internal dynamics of intermediatemass and massive stars at the interface between the convective core and the radiative interior. Stochastically excited gravitoinertial modes have also been detected in rapidly rotating massive stars (Neiner et al. 2012, 2020). Due to the difficulties in observing gmodes in solartype pulsators discussed above, it appears that performing 3D deepshell simulations of late Ftype stars constitutes a first approach to understanding the interplay between IGWs and rotation in their radiative interiors. While 3D simulations of Ftype stars have already been performed, they have only been aimed at studying the convectiveenvelope dynamics (Augustson et al. 2012) or dynamo properties (Augustson et al. 2013).
To this day, the range of stellar masses explored by IGW 3D hydrodynamic simulations remains limited. Brun et al. (2011) and Alvan et al. (2014; 2015, referred as A14 and A15, respectively), performed solarmodel simulations to evaluate the possibility to detect individual gmodes with helioseismic instruments. Browning et al. (2004) presented a simulation of the central regions of a 2 M_{⊙} A star, where IGWs are excited by convective motions in the core. Edelmann et al. (2019) extended the 2D simulations wave excitation analysis of 3 M_{⊙} stars from Rogers et al. (2013) into the 3D space. Augustson et al. (2016) studied the coreradiative interior interplay in 10 M_{⊙} O stars, while André (2019) examined the breaking of waves in 15 M_{⊙} O stars. In this work, we intend to extend the study of IGWs in lowmass stars to nonsolar models with varying rotation rates. We present the first deepshell hydrodynamical simulations of rotating Ftype stars. The layout of the paper is as follows. In Sect. 2, we present the model and cases for which we solve the hydrodynamical equations. In Sect. 3, we describe the global dynamics of our Ftype star model for the different rotation regimes we considered. The IGWs and gmode properties are extensively described and analysed in Sect. 4. In particular, we estimate the evolution with frequency of rotational splittings and period spacing. We then estimate the surface mode visibility of g modes in Sect. 5. We discuss the theoretical and observational perspectives opened up by this work in Sect. 6.
2. Numerical setup
2.1. Model equations
We used the ASH code (Clune et al. 1999; Brun et al. 2004) to solve the 3D hydrodynamic equations in the anelastic approximation. We consider a system of spherical coordinates (r, θ, ϕ) in a frame rotating at constant angular velocity Ω_{0} = Ω_{0}e_{z}. The reference density, pressure, temperature, and specific entropy are denoted as , , , and , while the fluctuations about this reference state are ρ, P, T, S. Following the prescription from Brown et al. (2012) and using the implementation presented in A14, we used the LantzBraginskyRoberts (Lantz 1992; Braginsky & Roberts 1995, LBR) equations to define the momentum equation. Indeed, the wave energy may be overestimated when using the traditional anelastic approximation. In the LBR formulation, the reduced pressure is considered instead of the fluctuating pressure P. The nonlinear momentum equation is therefore:
where v = (v_{r}, v_{θ}, v_{ϕ}) is the local velocity, g is the gravitational acceleration, and 𝒟 is the viscous stress tensor:
with e_{ij} = 1/2(∂_{j}v_{i}+∂_{i}v_{j}) the strain rate tensor and δ_{ij} the Kronecker symbol. In the anelastic approximation, the continuity equation is expressed as:
We assume a linearised equation of state and the zerothorder ideal gas law:
where γ is the adiabatic exponent, c_{p} is the specific heat per unit mass at constant pressure, and ℛ is the gas constant. Finally, the equation of conservation of internal energy is:
where κ_{r} is the radiative diffusivity and the volumeheating term related to the energy generation by nuclear burning. The ϵ profile is parametrised as a power law of , . The parameters ϵ_{0} and k are computed in order to have the integrated heating equal to the stellar luminosity at the top of the radiative zone. The parameters ν and κ are the effective diffusivities representing the unresolved momentum and heattransport by subgridscale (SGS) motions, while the diffusivity, κ_{0}, is set to carry the unresolved entropy eddy flux in the convective zone near the surface. We ensure that this flux does not play any role in the radiative zone by chosing a κ_{0} profile that decreases exponentially with depth (Miesch et al. 2000).
2.2. Models
We used the Modules for Experiments in Stellar Astrophysics (MESA, Paxton et al. 2011, 2013, 2015, 2018, 2019) to generate a 1D model of a 1.3 M_{⊙} star. We considered the model at the evolutionary stage when its hydrogen center mass fraction is 0.35. At this stage, the model luminosity, L_{⋆}, is 3.31 L_{⊙}, the effective temperature, T_{eff}, is 6353 K and the logarithm of the surface gravity, log g, is 4.2. The stellar radius R_{⋆} is 1.0465 × 10^{11} cm (1.5 R_{⊙}) and for the simulation, we considered the region spanning from 0.07 to 0.98 R_{⋆} (purposely omitting the relatively small convective core). We denote r_{bottom} and r_{top} the innermost and outermost radius of our simulation region, respectively. The upper edge of the convective core is treated as an impenetrable boundary and we consider the following boundary conditions at the top and bottom of the domain:
The interface between the radiative zone and the upper convective zone is located at r_{CZ} = 0.88 R_{⋆} (9.18 × 10^{10} cm). We note that the convective zone absolute thickness of our model is approximately 60% the thickness of the solar convective zone, while the total volume of the convective shell is approximately 68% larger than for the solar case. Concerning the energy generation parameter, we performed a fit on the MESA luminosity profile to obtain ϵ_{0} = 3.23 × 10^{−7} and k = 7.414. The dS/dr and g profiles provided to ASH are obtained from a polynomial fit of the MESA model. These two profiles yield the BruntVäisäla frequency profile N (given in Hz) according to the relation (e.g. Maeder 2009):
where we take a uniform c_{p} = 3.42 × 10^{8} erg g^{−1} K^{−1} for our input ASH profile. We know that IGWs are propagative in regions where N^{2} > 0 and evanescent in regions where N^{2} < 0. Close to the convective core, the dS/dr and N profiles are affected by the chemical gradient as it can be seen in Fig. 1, where we show the structural contribution, N_{t, MESA}, and the chemical contribution, N_{μ, MESA}, to the MESA BruntVäisälä profile, N_{MESA}. The extent of the convective envelope and the areas from the MESA models excluded from the 3D simulation domain are represented in the figure as well. We also compared N_{MESA} to the N_{ASH} profile obtained from the dS/dr and g fits. Although we are not able to reproduce the stiff frequency bump at the bottom of the radiative zone, we find only a 3.7% discrepancy when we integrate N_{ASH} and N_{MESA} along r. When we compute the , which is directly related to the gmode period spacing, we find a 11.6% discrepancy, which means that the gmode properties obtained from the 1D MESA model cannot be compared to the 3D simulations in a straightforward way. The trends and global properties could certainly be compared between the two set of profiles but in order to be as consistent as possible, we will compare the outcome of the nonlinear 3D simulations with the gmode frequencies predicted by a 1D oscillation code using the 1D ASH reference profile and Sect. 4.3. The and are computed using a NewtonRaphson algorithm in order to set the reference state at the hydrostatic equilibrium. In Fig. 2, we compare the profiles obtained with the NewtonRaphson algorithm and the reference profiles from MESA. Despite small deviations in the convective zone, we note that there is a good overall agreement.
Fig. 1. Radial profiles of the BruntVäisälä frequency N in the ASH model (black) and the MESA model (dashed grey). The structural and chemical contributions to N_{MESA}, N_{T, MESA}, and N_{μ, MESA}, are represented in dottedorange and dotteddarkblue lines, respectively. The dashed green vertical line shows the boundary between the radiative zone and the convective envelope, while the grey hatched regions are excluded from the simulation domain. The frequency corresponding to the ray paths of Fig. 15 are represented by the solid horizontal red, blue, and yellow lines. 
Fig. 2. Density (blue) and temperature (orange) profiles used as reference profiles for the simulation are represented with straight lines. MESA model profiles are shown in dashed lines for comparison. 
To study the impact of rotation on the model differential rotation in the convective zone and on gravity waves dynamics in the radiative interior, we ran cases with Ω_{0} = 1 (F1a and F1b cases), 3 (F3 case), and 5 Ω_{⊙} (F5 case), where we take Ω_{⊙} = 2.6 × 10^{−6} rad s^{−1}. All the considered cases verify that rad s^{−1}, with Ω_{K} the Keplerian critical breakup angular velocity and G the universal gravitational constant. Indeed, 5 Ω_{⊙} corresponds to ∼3.4% of Ω_{K} for the considered stellar model. The centrifugal deformation of the star can therefore be neglected (see Zahn 1992). We also emphasise the fact that, for this range of rotation rate, the relative difference in radius is below 0.1%, and that the position of the convective core differs by no more than 1% of the relative radius (Amard et al. 2019, and Amard, private communication), justifying our choice to use the same reference structure for all cases. As shown in Fig. 4, where we represent the distribution of photometric surface rotation periods from the Santos et al. (2021)Kepler catalog for stars in the range 6000 < T_{eff} < 6600 K and log g > 4 dex (using the values from Berger et al. 2020), the F1a and F1b cases are close to the slowrotator tail of the distribution for this population of stars. We emphasis that the F3 case is close to the median of the distribution (9.4 days) and the F5 case to the distribution maximum. Our choice of rotation rates is therefore a good representation of what is observed in the Kepler sample.
Fig. 3. Radial profile of the diffusivities, ν (straight lines) and κ (dashed lines), for the F1a (orange), F3 (blue), F1b, and F5 (red) cases. 
Fig. 4. Distribution of surface rotation periods P_{rot} measured by Santos et al. (2021) in the Kepler data, considering stars with 6000 < T_{eff} < 6600 K and log g > 4 dex (using values from Berger et al. 2020). Vertical colour lines indicate the periods corresponding to the rotation rates selected for our simulations. 
In order to simulate both the convective and the radiative zone, we take the diffusivities, ν and κ, as function of the radius and we use the following profiles:
with , , , σ_{t} = 0.08 the profile stiffness parameter, β_{ν} = 10^{−4}, and β_{κ} = 2 × 10^{−5}. The Prandtl number, Pr = ν/κ, is therefore not uniform along the stellar radius, from 0.25 at the top of the domain to 1.25 at the bottom. We show in Fig. 3, the diffusivity profiles adopted for the F1a, F1b, F3, and F5 cases. We emphasise that, with this choice of profile, we are able to resolve the motions in the convective envelope, while the abrupt drop of ν and κ in the tachocline limits the viscous and radiative damping of the IGWs in the radiative interior. We remind here that 3D stellar models must be compared with each other with the required level of caution. Indeed, as numerical limitations of simulations prevent us from reaching actual stellar regimes, it is important to take into account and discuss the role of the different characteristic fluid numbers exhibited by the different cases. Supercriticality considerations related to the choice of ν_{top} and κ_{top} will be discussed in Sect. 3.1 while the influence of the Rossby number Ro on the dynamic of the convective envelope will be exposed in Sect. 3.2. The properties of the four different simulations are summarised in Table 1.
Global properties for the considered cases.
2.3. Numerical resolution
In order to solve the hydrodynamic equations, following A14, the horizontal structure of the velocity and thermodynamic variables are expanded in spherical harmonics Y_{ℓ, m}(θ, ϕ), with ℓ the spherical degree and m the azimuthal number (see Table 1 for N_{θ} × N_{ϕ} resolution), while for the radial structure we use a nonuniformgrid with a finite difference approach. The grid we use has N_{r} = 1205 radial points. The radiative zone is resolved with 764 points while the convective zone is resolved with 441 points. As shown in Fig. 5, the radial resolution is significantly finer in the convective zone and in the overshoot region, where the diffusivity values, ν and κ, drop abruptly. The equations are solved using an explicit AdamsBashforth time integration scheme for the advection and Coriolis terms, while the diffusive and buoyancy terms are treated through a semiimplicit CrankNicolson method (Glatzmaier 1984; Clune et al. 1999).
Fig. 5. Spacing in the nonuniform radial grid used in the 3D ASH simulations. 
Once the simulation has been evolved over several tens of convective overturning times (4 < τ_{conv} < 5 days in our simulations), we obtain the flux balance between the different energy transport processes represented in Fig. 6. In our hydrodynamic setup, the total luminosity, L_{tot}, can be decomposed as
Fig. 6. Flux balance between radiative luminosity L_{rad} (dashed orange), kinetic energy luminosity L_{ke} (dashed red), enthalpy luminosity L_{en} (dashed grey), diffusivity luminosity L_{ν} (dashed green), and subgridscale eddy luminosity L_{ed} (dashed blue) for the F5 case. The total luminosity, L_{tot}, is represented in black. The top panel shows the flux balance for the whole simulation domain, while the bottom panel shows an enlargement in the convective zone and the overshoot region. 
where L_{rad} is the radiative flux, L_{ke} the kinetic energy flux, L_{ν} the diffusive processes energy flux, L_{en} the enthalpy flux, and L_{ed} the unresolved eddy flux (as described in Sect. 2.1). As expected, the energy transport is purely radiative below the tachocline. The unresolved energy flux becomes dominant near the top of the simulation domain (as explained in Sect. 2.1). In the middle of the convective zone, the enthalpy flux excess compensates for the inwards kinetic energy flux. At the interface between the radiative and the convective zone, we note that a significant amount of enthalpy is transported towards the interior due to overshoot mechanisms in the tachocline. As the timescales we consider in the simulation are much smaller than the thermal relaxation timescale (about 1 × 10^{5} yr) required for the system to reach a new equilibrium (Zahn 1991), we modified the κ_{r} value at the interface in order to increase the radiative flux and balance the enthalpy flux excess in this region (Miesch et al. 2000; Brun et al. 2011), thus easing the relaxation time.
3. Ftype star dynamics in 3D
To illustrate the different behaviours of the radiative zone and the convective zone, we represent a 3D volume of the F5 case in Fig. 7. We show the radial velocity, v_{r}, normalised by its rootmean square (rms) value at each given radius, . In the convective zone, large convective structures are shaped by the stratification of the reference state, with wide upward flows and thinner downwards flows. These downflow plumes act as an excitation piston when they interact with the stratified radiative zone. Below the interface, excited IGWs propagate towards the convective core and are reflected at the bottom of the radiative zone. We note the change of convection scales as we reach the base of the convective zone and the presence of concentric rings typical of lowfrequency IGWs (see Sect. 4). We must draw attention to the fact that, even if the hydrodynamical equations we solve are deterministic, everything occurs as if the IGWs are stochastically excited, since the properties of the plumes interacting at any given moment with the top of the radiative zone are not known a priori.
Fig. 7. 3D volume of the radial velocity, v_{r}, normalised by the rms velocity, (denoted v_rms in the colorbar) for the F5 case, from r = 0.07 to 0.94 R_{⋆}. Upward flows are in represented in yellow and downward flows in blue. 
3.1. Rotation and convection
In order to assess the degree of turbulence of the flows in the convective zone, we estimate the Reynolds number, Re, as:
where is the rootmeansquare (rms) radial velocity and L a characteristic length. We consider L = d_{CZ}, where d_{CZ} is the thickness of the convective zone. As expected from the chosen ν profiles, the most turbulent case is the F1b case, with a rms Re of 73 in the middle of the convective zone, followed by the F5 case with Re = 64. The F1a and F3 cases have Re = 20 and Re = 39, respectively. We underline that the action of the turbulent Reynolds stress is likely to be enhanced with increasing Re.
Convection is only possible above a critical value of the Rayleigh number Ra (e.g. Chandrasekhar 1961; Jones et al. 2009) and, therefore, a sufficient level of supercriticality needs to be sustained in the simulations for convection to dominate diffusive phenomena. Takehiro et al. (2020) confirmed that, in the anelastic approximation the critical value of Ra scaled as Ta^{2/3}, similarly to the Boussinesq case, where is the Taylor number. The diffusivity profiles, ν and κ, should scale as in order to keep the same level of supercriticality in each simulation but the limitation of available computing resources prevents from actually adopting this scaling (e.g. Augustson et al. 2012). Hence, we adopted, for the diffusivities, a scaling of between the F1a and F3 cases and between the F1a and the F5 cases, as a best compromise between numerical costs and physical constraints. In order to estimate the supercriticality level of each of our models, we computed the modified Rayleigh number defined by Takehiro et al. (2020), Ra^{*}, and we compared it with the critical value Ra_{c} we obtain when taking their M11R5 model as reference for the scaling (for their model Ta = 5.9 × 10^{6} and Ra_{c} = 4 × 10^{5}). We obtained Ra^{*}/Ra_{c} = 14.3, 23.9, 4.5, and 2.8 for the F1a, F1b, F3, and F5 cases, respectively, confirming that all of them are in a supercritical state. The supercriticality level achieved in the different cases is however significantly different, mainly because of the different rotation rates we impose. In order to assess how this influences the convection power spectrum, we therefore represent in Fig. 8 the spherical harmonic decomposition of the time average of the rms radial velocity, , for each case. The decomposition is summed over m in the top panel and over ℓ in the bottom panel. At low ℓ, the F3 case exhibit the largest values, followed by the F5 and F1b cases. It appears that the transition between the inertial range and the viscousdominated domain happens in the ℓ = 20 − 50 range. Beyond ℓ = 100, the spectrum is completely dominated by viscous diffusion. The inhibiting effect of rotation can be distinguished around the peak at ℓ = 30, where we note that the F1b decomposition peaks significantly higher than for the F5 case, although the two cases have identical κ and ν profiles. As expected due to its fastest rotation rate, the F5 case peaks at higher ℓ (Takehiro et al. 2020). Concerning the m decomposition, the velocity spectrum is flat for for the F1a and F1b cases at low and intermediate m. It is flat only at low m for the F3 and F5 cases, then increases and peaks between m = 10 and m = 30. Around the m ≈ 30 threshold, the velocity drastically decreases for all cases as m increases. The mspectrum reaches its maximal value for m = 17, 8, 20, 25 for the F1a, F1b, F3, and F5 cases, respectively. We note that these maximal values increase with rotation, following the trend identified by Takehiro et al. (2020) for critical azimutal numbers. We recall that due to the apparent mismatch with helioseismic solar convective velocity (Hanasoge et al. 2012), the socalled convective conundrum, absolute values for convective velocities obtained from simulations must be considered with care although the general trend they follow is consistent (Hanasoge et al. 2016; Hotta & Kusano 2021; Brun et al. 2022). The values of Re, Ta, Ra^{*}/Ra_{c}, and τ_{conv} are summarised in Table 2.
Fig. 8. Spherical harmonic decomposition of the time average of the rms radial velocity in the middle of the convective zone (r = 0.93 R_{⋆}) for the F1a (orange), F1b (cyan), F3 (blue), and F5 (red) cases. In the top panel, the rms components are summed over m and shown as a function of ℓ. In the bottom panel, the rms components are summed over ℓ and shown as a function of m. 
Dimensionless numbers in the middle of the convective zone, τ_{conv}, and r_{c} − r_{0} penetration depth.
Figure 9 shows the impact of rotation on the convective structures, represented through radial velocity, v_{r}, and temperature perturbation maps at r = 0.95 R_{⋆}. At this depth, v_{r}, and are clearly correlated: downward flows carry negative temperature perturbations, while upward flows carry positive temperature perturbations. At 1 Ω_{⊙}, the structures are only marginally affected by rotation, their shape does not depend much of the latitudinal position. Smaller scales are visible for the F1b case as it is more turbulent than F1a. On the contrary, when considering the F3 and F5 cases, the increase of Ω_{0} yields sharper patterns for the velocity field. Bananashape cells akin to Busse’s columns appear at this rotation rate (Busse 1970), showing evidence that the convection dynamics is strongly affected by the rotation rate of the model and the strength of the Coriolis acceleration. We notice that the structures tend to align perpendiculary to the equator. We note that in the F3 and F5 cases, largescale structures with large temperature perturbations (positive or negative) appear at high latitude. As it can be seen in Fig. 10, the radial profile is similar in the convective zone for the four cases. In the middle of the convective zone, we have , 3.7 × 10^{4}, 3.3 × 10^{4}, and 3.2 × 10^{4} cm s^{−1} in the F1a, F1b, F3, and F5 cases, respectively. We therefore confirm that the most turbulent model is the one with the highest convective velocities. The value in the radiative zone is associated with the amplitude of IGWs propagating in the stellar interior. Larger are reached in the radiative zone for the F3 and F5 cases. F3 and F5 have lower κ value in the radiative zone. Therefore, waves are less damped by thermal dissipation. We also notice that, due to the deeper position of the transition radius for the diffusivity drop, the mean stay at a level equivalent to the F3 case until r/R_{⋆} ≈ 0.6. At the bottom of the radiative zone, however, values from the F3 case are comparable to the F1b case.
Fig. 9. Radial velocity v_{r} (left) and temperature perturbation (right) at depth r = 0.95 R_{⋆} for the F1a, F1b, F3, and F5 cases (top to bottom). The latitudinal mean temperature perturbation, m = 0, is subtracted from T. 
Fig. 10. profiles for the F1a (orange), F1b (cyan), F3 (blue), and F5 (red) cases. The vertical dashedgrey line correspond to the depth where we examine the mode signatures for the F5 case in Sect. 5. 
Finally, it should be reminded that convection is intriscally a dynamical process, with cells evolving over time. The interaction and combination between turbulent flows shape upwards and downwards travelling convective structure. Downwards flows with the largest amount of power give rise to the convective plumes interacting with the radiative zone. The plumes act similarly to a piston as they overshoot in the radiative and inject power in the inner regions in the form of IGWs.
3.2. Differential rotation in the convective zone and the sharp tachocline
We now turn to considering the differential rotation regimes achieved in the simulations. Following Brun et al. (2017), we define the fluid Rossby number as:
where ∇ × v is the vorticity of the flow. Brun et al. (2017) expect a transition between the solar (fast equator, slow poles) and antisolar rotation (slow equator, fast poles) regime at Ro_{f} ≈ 1 (see also Gastine et al. 2014; Guerrero et al. 2019; Warnecke & Käpylä 2020). At very low Ro_{f}, as a consequence of the TaylorProudman theorem, the rotational regime becomes cylindrical.
We also compute the stellar Rossby number Ro_{s} (e.g. Noyes et al. 1984; Corsaro et al. 2021) and the convective Rossby number Ro_{c} (Gilman & Glatzmaier 1981). Here, Ro_{s} is given by
where the rotation period is P_{rot} = 2π/Ω_{0}, and τ_{conv} is the convective turnover time that we compute as the ratio of the thickness of the convective zone in the simulation domain, d_{CZ}, and the mean in the convective zone, . For Ro_{c}, we use the following definition
where Ra = ( − ∂ρ/∂S)(∂S_{tot}/∂r)gL^{4}/ρκν, with . Here, again, for the characteristic length, L, we consider d_{CZ}. As expected by Brun et al. (2017), Ro_{f}, Ro_{s}, and Ro_{c} scale in a similar way, except for the Ro_{c} value of the F1b case, which is smaller than for the F1a case. These values are between 3 and 6 for the 1 Ω_{⊙} cases and approach unity for the fastest rotating cases.
The values of Ro_{f}, Ro_{s}, Ro_{c}, and Ra in the middle of the convective zone are summarised for each case in Table 2.
We present the Ro_{f} profile in the convective zone in Fig. 11. Considering the v_{ϕ} component of the velocity field averaged over longitude and time, we also computed the differential rotation state for the different cases. The (r, θ) differential rotation profiles are represented in Fig. 12. For the F1a and F1b cases, the simulation exhibits a shellular rotation profile which is consistent with the fact that Ro_{f} in the middle of the convective zone is significantly larger than 1, as seen in Fig. 11. The F3 case behaviour is still clearly antisolar, with a Ro_{f} of 1.9 in the middle of the convective zone. Interestingly, the transition from prograde to retrograde flows (in the corotational frame), at every latitude for F1a and F1b and close to the equator for F3, intervene close to r_{CZ}. The F5 case has R_{of} ≈ 1.4 in the middle of the convective zone and the structure of the flows observed in Fig. 12 suggest that the model is in a transitional regime from an antisolar to a solar differential rotation regime. In our F5 case, we notice an asymmetry between behaviours at high latitudes, with slow flows close to the north pole and fast flows close to the south pole. Faster rotating Fstar models published by Augustson et al. (2012) confirm the change of rotation regime towards TaylorProudman constrained states for low Ro. We recall that their 1.3 M_{⊙} model rotating at 10 Ω_{⊙} exhibits a Ro_{f} value of 0.84 and a solar differential rotation regime, which is consistent with the predictions from Brun et al. (2017).
Fig. 11. Fluid Rossby number Ro_{f} profiles in the convective zone for the F1a (orange), F1b (cyan), F3 (blue), and F5 (red) cases. 
Fig. 12. Differential rotation for F1a (top left), F1b (top left), F3 (bottom left), and F5 (bottom right) cases. Retrograde and prograde (relatively to the corotational frame) flows are shown in blue and red, respectively. The dashed black line corresponds to r_{CZ} and the dashed grey lines are isocontours. 
3.3. Overshoot
The convective motions overshoot into the radiative interior, leading to mixing and IGWs generation. In Fig. 13, we illustrate how the correlation between v_{r} and change in the overshoot region relatively to the convective zone. We show the v_{r} and maps at depth r = 0.86 R_{⋆} for the F1a case. In the convective zone, v_{r} and were correlated, that is, the upwards flows were related to positive temperature perturbations and reciprocally. On the contrary, in the overshoot region, the two quantities are now anticorrelated: downwards travelling plumes penetrating the top of the stablystratified radiative zone are associated with positive temperature perturbations. This is expected from our understanding of penetrative convection in stellar interiors (Zahn 1991; Brummell et al. 2002).
Fig. 13. Radial velocity v_{r} (left) and temperature perturbation (right) at depth r = 0.86 R_{⋆} for the F1a case. The latitudinal mean temperature perturbation was substracted from T. 
Figure 14 shows the enthalpy flux at the bottom of the convective zone and in the overshoot region for the F1a, F1b, F3, and F5 cases. Following Brun et al. (2011, 2017), we computed, for every latitude, the boundaries r_{c} and r_{0} where the enthalpy flux becomes negative and where the enthalpy flux is only 10 % of the maximal absolute value in the overshoot region at this latitude, respectively. We note that the shape of the r_{0} profile does not change significantly with the rotation, however, the r_{c} profile differs significantly in the F3 and F5 cases compared to the F1b case. We also note for fast rotating models, the enthalpy flux intensity in the overshoot region is concentrated towards the equator likely due to the stronger alignment of convective rolls for these stars. By computing the difference of the mean value of r_{c} and r_{0} at latitudes below 55^{o} for each case, we obtain the mean penetration depth of the overshooting flows, , in the same way as Brun et al. (2017). We obtain , 0.23, 0.19 and 0.15 H_{p} where H_{p} ≈ 8.8 × 10^{9} cm is the pressure scale height at the base of the convective zone. The comparison between F1a and F1b shows that, at constant Ω_{0}, the overshoot depth increases with Re but, as expected, the penetration depth of the overshooting flow is significantly reduced as rotation increases. The values are summarised together with the dimensionless numbers in Table 2. These overshooting motions contribute to the excitation of the IGWs. In particular, (Pratt et al. 2017) showed that two characteristic layers of penetration could be distinguished for overshooting plumes, the deepest one corresponding to the excitation region of the IGWs.
Fig. 14. Overshoot region (blue) and bottom of the convective zone (red) in the F1a (top left), F1b (top right), F3 (bottom left), and F5 (bottom right) models, for θ spanning from 30 to 150 ^{o}. The white line show the radius r_{c} at which the enthalpy flux crosses zero at a given latitude while the orange line signals the radius r_{0} where the local enthalpy flux is equal to a tenth of the maximal enthalpy flux at the corresponding latitude. 
4. IGW properties
In this section, we study in detail the properties of the IGWs that are generated by the interaction of the convective motions with the top of the radiative zone. The dispersion relation for IGWs is (Press 1981):
where k is the norm of the wave vector k = k_{r} + k_{h}. The norm of the horizontal wave vector k_{h} = k_{θ} + k_{ϕ} is:
Concerning gravitoinertial waves, when the Coriolis parameter f = 2Ω_{0} is such that f ≪ N, the effect of the Coriolis acceleration on the wave behaviour is negligible and we retrieve the puregravity case. Hence, we expect highfrequency waves propagation not to be affected by the rotation regime of the different cases. However, as we increase Ω_{0}, the range of waves significantly affected by inertial effects expands.
4.1. Raytracing
In order to better understand the expected IGWs behaviour at different frequencies for our Ftype model, we start by computing the propagation path of IGWs at different frequencies following a raytracing Hamiltonian method (Gough 1993). As k_{h} depends only on ℓ and r, we set θ = 0 and k_{θ} = 0 and we consider only a 2D problem in the equatorial plane. However, as θ and ϕ can be freely interchanged in this specific case, we could have chosen to place ourselves in any meridian plane without any modification to our subsequent analysis. To perform the raytracing, we must numerically solve the following set of equations:
where W(r, θ, k_{r}, k_{ϕ}) = ω is the Hamiltonian of the considered system. Using the dispersion relation of Eq. (15), we compute the propagation path for three IGWs with different properties, ℓ = 1 at ω = 10 μHz, ℓ = 1 at ω = 140 μHz, ℓ = 2 at ω = 319 μHz. The integration times, Δt, are (respectively) 15, 5, and 1 days. The radial resonant cavities corresponding to the three chosen frequencies are represented along with the N profile in Fig. 1. The obtained propagation paths are shown in Fig. 15. As expected from the N profile shown in Fig. 1, highfrequency IGWs are trapped in the deepest regions of the radiative zone. Highfrequency IGWs propagate much faster in the radiative interior. Lowfrequency IGWs have a characteristic spiralling trajectory and take more time to reach the edge of the convective core.
Fig. 15. Propagation path computed with raytracing for three IGWs: ℓ = 1, ω = 10 μHz (top), ℓ = 1, ω = 140 μHz (middle), ℓ = 2, ω = 319 μHz (bottom). The integration time, Δt, considered for each ray is specified. From the surface to the center, the black circles represent the position of the stellar surface, the bottom of the convective envelope and the edge of the convective core. The grey area corresponds to the part of the radiative zone where the wave is evanescent. 
In order to present more evidence of IGWs propagating in the radiative zone of our simulations, we considered the v_{r} temporal evolution in an equatorial cut of the F1a case. We applied a finite impulse response (FIR) filter with a passband centered at 140 μHz in order to isolate the n = 2, ℓ = 1 mode. We compare in Fig. 16 the result of the filtering with the propagation path obtained by solving the Hamiltonian system. Considering the sign of , the dipolar structure of the mode oscillation is clearly visible after the filtering, as well as the position of the two radial nodes of the mode. What we see in this filtering can be interpreted as the pattern of a ray interfering with itself (Gough 1993). Due to the diffusive effects related to the κ and ν profiles, the characteristic shape of the mode spreads when compared to the path obtained with the raytracing. However, having such a good qualitative agreement between our filtered model and the shape of the waves obtained with the raypath theory gives us confidence in the realism of internal waves excited in our simulations.
Fig. 16. Equatorial cut for for the radiative zone in the F1a case, with application of a temporal filtering centered at 140 μHz. The 2D propagation ray obtained by numerically solving the system of Eq. (17) and shown in the middle panel of Fig. 15 is overplotted in white. We note that there is a good qualitative agreement. 
4.2. IGWs power spectrum
To obtain the IGW power spectrum, we produced outputs of v_{r}(r, θ, ϕ, t) at a mean sampling dt low enough to have N_{max} < ω_{N}, with ω_{N} the Nyquist frequency of the output signal and N_{max} the maximal value of the N profile. In our case, this means that we are careful to have dt ≤ 1250 s. The length of the time series we use for each case are given in Table 1. We expand v_{r}(r, θ, ϕ, t) into a spherical harmonics representation to obtain v_{r}(r, ℓ, m, t). This representation is more suitable to study the modes as their properties directly depend of ℓ and m (as for example mode period spacings or rotation splittings). From the Fourier transform of v_{r}, we then obtain the power spectrum E_{ℓ, m} as
and, by summing over the m for each ℓ, we have the power spectrum, E_{ℓ},
Some additional practical details concerning the way we obtain E_{ℓ, m} and E_{ℓ} are provided in Appendix A.
The E_{ℓ} power spectra at r = 0.18 R_{⋆} obtained for the F1a, F1b, F3, and F5 cases are represented in Fig. 17. The IGWs spatial damping rate in a nonadiabatic medium increases with ℓ and decreases with ω (e.g. Zahn et al. 1997). As pointed out by Alvan et al. (2015), two regions can therefore be distinguished in each spectrum: the bottom right part of the spectrum (high degree and/or low frequency) shows progressive waves which are significantly damped before reaching the edge of the convective core while in the upper left part of the spectrum (low degree or high frequency, or both), IGWs reflect on the core edge with enough amplitude to form standing waves, or g modes. In this region, gmodes are distributed along ridges in relation to their radial order, n. Following Ahuir et al. (2021), we characterise the frequency cutoff between the two regions, ω_{c}, through the following relation
Fig. 17. E_{ℓ} power spectrum at r = 0.18 R_{⋆} for the F1a (top left), F1b (top right), F3 (bottom left), and F5 (bottom) case. The mode frequencies computed with GYRE are overplotted (black dots) for ℓ = 1 to ℓ = 40. In each panel, the white line highlights the separation between the gmode region (top left) and the progressive waves region (bottom right). 
where we find that adopting a critical damping parameter τ_{c} = 0.02 correctly describe the cutoff profile observed in Fig. 17.
As expected from Fig. 10 and the various level of in the radiative interior, the IGWs mean amplitude is significantly larger in F5 and F3 than F1a. In particular, lowfrequencymode amplitudes are significantly larger in the F3 and F5 cases. Being a more turbulent version of F1a, F1b exhibits highfrequency modes of larger amplitude than in any other case, but the excitation of lowfrequency modes is similar to what is observed in F1a.
In order to compare the mode excitation rate degree by degree, we compute the power index E_{ℓ, tot}, taken as the summation of the E_{ℓ} component over the frequency bins, ω_{i},
For 1 ≤ ℓ ≤ 10, we represent in Fig. 18 the values obtained for E_{ℓ, tot} at depth r = 0.8 R_{⋆}, in the top of the radiative zone, as well as the ratio α = 100 × E_{ℓ, tot}(r = 0.8 R_{⋆})/E_{ℓ, tot}(r = 0.9 R_{⋆}). This allows us to compare on one hand the relative excitation level of each degree ℓ, and to evaluate on the other hand the efficiency of power injection from the convective motions to the IGWs. As expected from it being more turbulent than F1a, the excitation level of the waves is more important in F1b than in F1a. The excitation level decreases with ℓ in F1a while it is relatively flat in ℓ for F1b, suggesting that intermediate degree ℓ have been more efficiently excited by the more turbulent flows of F1b (see Fig. 8). As already shown by Fig. 17, the excitation level of the IGWs is the highest in F5, with a peak at ℓ = 4. It is remarkable to note that despite the F3 case being more dissipative than the F1b case, the transmission of power and the wave excitation is more important. It is interesting to note that, for the F5 case, the α ratio is the highest for ℓ = 1 (with α = 0.12 %), suggesting an efficient transfer from convection, but remains the degree with the lowest excitation level in this case. For each case, α tends to decrease as ℓ increases. In the F5 case, the power injection peaks at ℓ = 4. The enhanced power injection in the IGWs as Ro_{c} (see Table 2) decreases has been predicted by Augustson et al. (2020) and the result we present in Figs. 17 and 18 are in agreement with their theoretical considerations. However, the fact that intermediate ℓ are more excited than low ℓ might be considered surprising as modes of higher ℓ have increased inertia and require more energy from convection to be excited (Provost et al. 2000). As Fig. 8 shows, the intermediate ℓ have increased convective velocities compared to low ℓ. The absolute power in convection is nevertheless not a sufficient explanation alone as the slowly rotating cases have a similar convective spectrum in the low to intermediate ℓ. The mechanism enhancing power injection from rotation into the modes seems particularly efficient for intermediate ℓ.
Fig. 18. Amount of power injected in IGW. Top: total power E_{ℓ, tot} for 1 ≤ ℓ ≤ 10 at r = 0.8 R_{⋆} for the F1a (orange diamonds), F1b (cyan circles) F3 (blue crosses), and F5 (red pentagons) cases. Bottom: ratio between E_{ℓ, tot}(r = 0.8 R_{⋆}), the radiative zone, and E_{ℓ, tot}(r = 0.9 R_{⋆}), convective zone. The same symbols and colour coding are used as in previous figures. 
We will discuss in more detail the surface amplitudes and signatures of the modes in Sect. 5.
4.3. gmodes frequencies and eigenfunctions comparison with outputs from a 1D oscillation code
We use the GYRE code (Townsend & Teitler 2013; Townsend et al. 2018; Goldstein & Townsend 2020) to compute the expected oscillation frequencies from the 1D input ASH profiles, for ℓ = 1 to ℓ = 40, in a case without rotation (referred later as the Ω_{⋆} = 0 GYRE run). As shown by Fig. 17, the GYRE computed frequencies are globally in good agreement with what we observe in the 3D simulations. In particular, we can clearly see that the ridge structure for modes of different degrees ℓ but same orders n coincide in the 3D simulation and in the GYRE computation.
We also compare for some modes the ξ_{r} displacement eigenfunctions with the outcomes of the 3D simulations. The eigenfunctions for ℓ = 5, n = 7 and n = 12 are shown in Fig. 19. The node position for these modes is in good agreement between GYRE and the 3D simulation, confirming the type of the waves excited in the 3D simulations as being gravity waves. This is very satisfactory and allows us to advance further our analysis of their properties.
Fig. 19. GYRE ξ_{r} eigenfunctions (blue) compared with ASH outputs for ℓ = 5, n = 7 and ℓ = 5, n = 12. The depths where ASH outputs are sampled are signaled by orange squares. 
4.4. Period spacing
For n ≫ 1, asymptotic m = 0 g modes are evenly spaced in period. The asymptotic period spacing for consecutive highorder gmodes is (Tassoul 1980; Provost & Berthomieu 1986):
where, for a mode of frequency ω, r_{1} and r_{2} are the inner and outer boundary of the corresponding resonant cavity, respectively, defined by N(r_{1}) = N(r_{2}) = ω.
In Fig. 20, for each case, we represent the E_{1} (0.38 R_{⋆}, ω) power spectra and the corresponding E_{1} (0.38 R_{⋆}, P = 1/ω) power spectra scaled in period. The middle panels of Fig. 20 show the modes are almost equidistant in period. Therefore, we consider 50 < P < 520 min for ℓ = 1 and 50 < P < 300 min for ℓ = 2, and, in order to highlight this ΔP_{ℓ}periodicity in E_{ℓ} (0.38 R_{⋆}, P), we use an approach similar to the one presented by García et al. (2007), by computing the LombScargle periodogram (Lomb 1976; Scargle 1982) of E_{ℓ} (0.38 R_{⋆}, P) for periods in the LombScargle periodogram spanning from 5 to 60 min. Each LombScargle periodograms is normalised by its standard deviation σ. We perform the same operation for the E_{2} power spectra. The obtained LombScargle periodograms are shown in the bottom panels of Fig. 20 for E_{1} and in Fig. 21 for E_{2}. We use the periodograms to identify the period spacings, ΔP_{1} and ΔP_{2}, which we compare to the asymptotic value, min and min. We find ΔP_{1} = 45.4, 42.3, 45.4, and 45.8 min and ΔP_{1} = 24.3, 25.8, 25.5, and 26 min for the F1a, F1b, F3, and F5 cases, respectively. These values are all below the asymptotic references, which is expected as we considered loworder modes for which period spacings are usually smaller than for higher order modes. It is nevertheless satisfactory that we are able to identify without ambiguity the periodicity related to the mode pattern in the LombScargle periodograms. We also note that the shape of the peaks yielded by the periodogram significantly changes with rotation, as mode splittings is considerably more visible on the frequency range we consider for fast rotating models and the split components are not equidistant in period.
Fig. 20. Top panels: power spectra E_{1} (0.38 R_{⋆}, ω). Middle panels: power spectra E_{1} (0.38 R_{⋆}, P). Bottom panels: LombScargle periodograms computed from E_{1} (0.38 R_{⋆}, P). Each periodogram is normalised with its standard deviation, σ and the obtained ΔP_{l} value is shown by a vertical dashed blue line. From left to right, the represented cases are F1a, F1b, F3, and F5. 
The difference in amplitudes for lowfrequency modes between slowrotating cases F1a and F1b, and fast rotating cases F3 and F5 clearly appears in Figs. 20 and 21. On these two figures, we also note that highfrequency modes are less excited in the F3 and F5 cases than in the F1a and F1b, suggesting that the mode excitation efficiency is shifted towards low frequencies.
Fig. 21. LombScargle periodograms computed from E_{2} (0.38 R_{⋆}, P). Each periodogram is normalised with its standard deviation, σ and the obtained ΔP_{2} value is shown by a vertical dashed red line. From left to right, the represented cases are F1a, F1b, F3, and F5. 
We also compare the evolution of ΔP_{ℓ} along the period with the outputs of GYRE computations. In order to take into account the Coriolis effect in our comparison, we use GYRE to compute the predicted mode frequencies for ℓ = 1; m = { − 1, 0, 1} with the traditional approximation of rotation (TAR), assuming a solid body rotation at 5 Ω_{⊙}. We compare the outputs with the F5 case. The TAR treatment is useful when the ratio ω/2Ω becomes to small to take rotation into account using only a perturbative approach (Lee & Saio 1997). In the TAR, the Coriolis force is included in the equations to solve but its radial component is neglected. The TAR is therefore wellsuited to study wave behaviour in stably stratified radiative interiors but breaks in convective zones. The GYRE run using the TAR will be referred in what follows as the 5 Ω_{⊙} GYRE TAR run.
In order to measure the period spacings in the simulation, we considered a 280day long time series of the F5 case (see Table 1) and we compute the corresponding E_{ℓ, m}(r = 0.13 R_{⋆}, ω) power spectrum (see Eq. (18)) to fit a Lorentzian profile for each m component of the ℓ = 1 modes visible at this depth. The Lorentzian profile is fitted with a Markov chain Monte Carlo (MCMC) approach, assuming a χ^{2} statistics with two degrees of freedom, implemented with the emcee module (ForemanMackey et al. 2013). We compute the corresponding periods, P_{n, 1, m}, and period spacings, ΔP_{n, 1, m}, relative to the azimutal number, m, taken as:
The results of this analysis are shown in Fig. 22. We compare them with the ΔP_{n, 1, m} obtained with the 5 Ω_{⊙} GYRE TAR run and the perturbative approach applied to the Ω_{⋆} = 0 GYRE runs. At long periods, the ΔP_{n, 1, m}, obtained considering the TAR significantly differs from what the perturbative approach predicts. Since the Coriolis force is fully taken into account in the equations solved by the ASH code, the ΔP_{n, 1, m} evolution observed in the 3D simulation follows much better the 5 Ω_{⊙} GYRE TAR run. To assess this more quantitatively, we define ℳ, a leastsquare distance metric:
Fig. 22. ΔP_{1} as a function of modes period in the F5 case, for azimuthal numbers m = −1 (yellow), m = 0 (blue), and m = 1 (red). The ΔP_{1} obtained from the Ω_{⋆} = 0 GYRE run and the perturbative approach are represented as grey hexagons while the white crosses shows the ΔP_{1} obtained for m = { − 1, 0, 1} with the 5 Ω_{⊙} GYRE TAR run. The black crosses signal the ΔP_{1} position for m = { − 1, 1} in the asymptotic approximation of the perturbative method. The measured uncertainties on period and ΔP_{1} are represented. We show only fitted modes for which we are able to measure ΔP_{1} with an uncertainty below 1.75 min. 
where the are the uncertainties measured for each value of ΔP_{n, 1, m (ASH)}. As a leastsquare distance metric, ℳ is equivalent to a loglikelihood with variables following a Gaussian distribution. Considering both approaches, we can therefore determine the most likely by computing δ = ℳ_{p.a} − ℳ_{TAR}, where ℳ_{TAR} is the value obtained with the TAR and ℳ_{p.a} is the value obtained with the perturbative approach. A positive value of δ would favour the TAR hypothesis while a negative value would favour the perturbative approach hypothesis. We obtain δ = 416, which strongly favours the TAR hypothesis. This provides another evidence that the period spacings that we measure in the simulation are more compatible with the values obtained using the TAR.
4.5. Rotational splittings
In the case of uniform rotation and in the slow rotation limit, ω ≫ 2Ω_{0}, the rotational splitting for a mode component with frequency ω_{nℓm} is given by (e.g. ChristensenDalsgaard 2008)
with, in the inertial frame,
and in the corotating frame,
The parameter β_{nℓ} is:
where ξ_{r} and ξ_{h} are the radial and horizontal displacements, respectively, and . In the case of highorder gmodes, we can consider
Considering a ℓ = 2 mode, we provide an example of how rotational splittings appear in the power spectrum of the F5 case in Fig. 23. In the corotating frame, as expected from Eq. (27), the mcomponent frequency decreases as m increase. In order to measure the rotational splittings in the simulation, we use the mode frequencies measured in Sect. 4.4 for ℓ = 1 modes. For each fitted component m, we take the Lorentzian central frequency as the m component frequency and we use Eqs. (27) and (29) to compute the corresponding Ω value. Figure 24 shows the result of this analysis. We show the value inferred by using the m = −1 or the m = 1 alone, along with the value inferred by considering the averaged mean of the two measurements. We are able to fit almost every order with period P below 1800 min. The one for which we do not provide a measurement are those who are close to a node ridge at this depth of the simulation and, therefore, we have therefore very low amplitude in the considered power spectrum. In most of the cases, the uncertainties over the frequency we measure are significantly below the 10% interval around Ω_{0}. As we are constrained by the resolution of the time series we use, we note that the uncertainties on the inferred Ω increase as we go towards long periods (low frequencies). Obviously, the asymptotic approximation does not hold for modes of lowest orders and it is not possible to retrieve the correct rotation rate just by using the β_{nℓ} approximation yielded by Eq. (29) with such lown modes.
Fig. 23. Example of rotational splitting observed in the F5 case for a ℓ = 2 mode, with the m = −2 shown in brown, the m = −1 in yellow, the m = 0 in blue, the m = 1 in red, and the m = 2 in orange. 
Fig. 24. Rotation rate Ω inferred from the asymptotic approximation of Eqs. (27) and (29) with δ_{n, 1, m} measurements, in the F5 case. The position of the true rotation value of the model is highlighted by the dashed grey line, with the 10% error interval in grey. Rotation rate measured only with δ_{n, 1, −1} and δ_{n, 1, 1} are represented in yellow and red, respectively, while the averaged mean of the two values is shown in blue. The error bar for each Ω measurement is also represented. 
5. Mode tunneling and surface visibility
We go on to study the radial dependency of the mode amplitudes in the power spectrum. In Figs. 25 and 26, we represent for F1b and F5 the radial evolution of E_{ℓ} for ℓ = 1 and ℓ = 4, respectively, and for P < 1440 min. The boundary between the resonant cavity and the evanescent region is represented by the white line. As expected, the power level is significantly higher in the convective envelope because of the contribution of convective motions. The location of the mode nodes draw dark ridges of low amplitudes in the spectrum. As already underlined in Sect. 4.2, the F3 and F5 cases exhibit more power in highperiod (lowfrequency) modes. The splitted character of the ℓ = 1 modes at high period (low frequency) appears for the F3 and F5 cases. As expected from Eq. (27), in the corotational frame, mode splittings appear less clearly for ℓ = 4 modes. For the E_{4} power spectrum, the separation between gmodes and progressive waves is visible in Fig. 26 and located close to P = 800 min. Therefore, modes with a radial order below n = 50 are visible in the figure, while the range presented for the ℓ = 1 modes in Fig. 25 covers 1 ≤ n ≤ 30. As underlined in Sect. 4.2 (see Fig. 18), modes with ℓ = 4 are particularly excited in the F5 case. The bottom right panel of Fig. 26 shows that, for 200 < P < 800 min and above r = 0.88 R_{⋆}, the mode signature is clearly visible above the white line and among the convective flows contribution.
Fig. 25. Power spectrum for ℓ = 1, E_{1}, radial evolution in the F1b (left) and F5 (right) cases. In each panel, the white line represents the boundary of the mode resonant cavity. The range of radial orders visible in the figure is 1 ≤ n ≤ 30. 
Fig. 26. Power spectrum for ℓ = 4, E_{4}, radial evolution in the F1b (left) and F5 (right) cases. In each panel, the white line represents the boundary of the mode resonant cavity. The cutoff between gmodes and progressive waves is visible around 800 min. Below this cutoff, the range of visible radial orders is 1 ≤ n ≤ 50. 
As illustrated with ℓ = 4 modes, the F5 case presents evidence of low and middegree modes tunneling through the convective zone and therefore we look for gmode signatures in the convective velocity signal near the top of the domain, where we still have of the order of 1.5 × 10^{4} cm s^{−1} (see Fig. 10). The F5 case is the only model where clear gmode patterns are observable in the convective envelope. Using the same periodogram method presented in Sect. 4.4, we look for periodicity in the v_{r} signal at r = 0.97 R_{⋆} for 1 ≤ ℓ ≤ 10, considering the period range 250 to 450 min. We are able to detect the gmode signature of ℓ = 3, 4, 5, 6, and 7 modes. This is particularly interesting for ℓ = 3 modes which are still observable in diskintegrated observations. As an illustration, we show in Fig. 27 the LombScargle periodograms computed for ℓ = 3 and ℓ = 5 modes, for the F5 case. The clearest detection is the ΔP_{5}, with a peak height at 7.1 σ against 4.4 σ for ΔP_{3}, and 5.8 for ΔP_{4}. We report in Table 3 the measured ΔP_{ℓ} and corresponding peak heights. They are compared with the asymptotic computed with Eq. (22).
Fig. 27. LombScargle periodograms computed from E_{3} (0.97 R_{⋆}, P) (top) and E_{5} (0.97 R_{⋆}, P) (bottom). The periodogram is normalised with its standard deviation σ and the obtained ΔP_{3} value is shown by a vertical dashed red line. 
Asymptotic , ΔP_{ℓ} measured at r = 0.97 R_{⋆} and corresponding peak heights for the F5 case.
Finally, we estimate the bolometric luminosity perturbation for individual g modes. In what follows, all the considered perturbations, δL_{⋆}, δR_{⋆}, and δT_{eff} are related to the action of an individual mode. From the StefanBoltzmann law we have (e.g. Samadi et al. 2010)
and, following Townsend (2003a)
where ∇_{ad} = 0.4 is the adiabatic temperature gradient, and is a dimensionless frequency defined as
To compute the displacement δR_{⋆} related to a single mode, we simply consider:
where is the radial rms velocity (measured in the E_{ℓ} power spectrum) of the considered mode.
For ℓ = 3, 4, and 5, that is the degrees with the most excited individual modes, we compute the luminosity perturbation δL_{⋆}/L_{⋆} as a function of for the modes with the largest amplitude at the top of the radiative zone. These three modes happen to have close eigenfrequencies, 41.5, 41.5, and 39.4 μHz, respectively, and correspond to n = 19, 23, and 28. We consider both the mode radial velocity at the top of the radiative zone, r = 0.87 R_{⋆}, and near the top of the domain, r = 0.97 R_{⋆}. The results of this analysis are represented in Fig. 28. When using the estimate at the top of the radiative zone, we find corresponding δL_{⋆}/L_{⋆} of a few partpermillions (ppm) for the ℓ = 4 and ℓ = 5 mode. The value for the ℓ = 3 mode is significantly smaller, close to 0.3 ppm. To maintain a propagative behaviour above the upper limit of the BruntVäisala resonant cavity, gravitoinertial waves related to the modes would have to be subinertial (Mathis et al. 2014), which is not the case here, as for the F5 case, we have 2Ω_{⋆} = 4.138 μHz. Nevertheless, as their eigenfrequencies are close to 10Ω_{⋆} where the inertial contribution still has a significant influence on the wave behaviour (Fig. 22 and Mathis, private communication), it is likely that the superinertial character of these modes increases the characteristic length of their evanescent tail in the convective envelope, compared to slower rotating case. The luminosity perturbation for mode velocities near the top of the domain are all below 1 ppm. The δL_{⋆}/L_{⋆} estimated in our analysis are several orders of magnitudes below what is typically observed for unstable modes in γ Doradus stars. This underlines that the stochastic excitation of stable modes remains much less efficient than unstable mechanisms. We remind that our model has been chosen to be compatible with the existence of solartype pmode oscillations, but has T_{eff} and log g values that locate it close to the γ Doradus instability strip lower boundary. However, evidence were presented above that gmode may efficiently tunnel through the convective zone and may be detected through their ΔP_{ℓ} signature. This leaves open the perspective to detect stochastically excited g modes and gravitoinertial modes in late Ftype stars.
Fig. 28. Luminosity perturbation estimation as a function of v_{r} for a ℓ = 3 mode at 41.5 μHz (red), a ℓ = 4 mode at 41.5 μHz (blue), and a ℓ = 5 mode at 39.4 μHz (yellow). The diamonds correspond to the radial velocity we measure in the F5 case at r = 0.87 R_{⋆} and the dots to the same quantity at r = 0.97 R_{⋆}. 
6. Discussion and conclusion
In this work, we presented and studied the first deepshell 3D simulation of IGWs excited in the radiative interior of a Ftype solartype star by convective motions. We considered a 1.3 M_{⊙} model for which we ran simulations at 1, 3, and 5 Ω_{⊙}, which are rotation rates representative of Ftype stars surface rotation rates observed by Kepler (see Fig. 4). We described the effect of rotation on the convective structure and on the differential rotation regime in the convective zone. The properties of the stochastically excited IGWs were extensively studied. In particular, we showed that the excitation rate of lowfrequency waves was significantly higher for fast rotating cases. We compared the eigenfrequencies of the modes in the simulations with the result of computations performed with the 1D oscillation code GYRE and we found a good agreement between the two approaches. We verified that the internal rotation rate could be inferred from rotational splittings in the asymptotic regime. Taking advantage of the regular period spacing of g modes, we computed the LombScargle periodograms of the power spectra to detect this periodicity and we compared the measured values with asymptotic predictions. Studying the periodspacing evolution at long periods in the F5 case, we were also able to observe the effect of the Coriolis force on the behaviour of highorder modes. Finally, we tried to quantify the possibility of observing gmode signatures in actual Ftype solartype pulsators. We showed (see Fig. 27 and Table 3) that for intermediate degrees ℓ = 3, 4, 5, 6, and 7, the mode signature could still be detected near the top of the simulation domain. As highlighted by Fig. 28, the most efficiently excited modes have a surface velocity of ∼10^{−1} cm s^{−1}, corresponding to a bolometric luminosity perturbation of ∼10^{−1} ppm, for the frequency range of these modes (∼40 μHz). It should be underlined that these values are much larger than what was observed in the solar models from A14, where the maximal surface velocity (obtained with the most turbulent model) was ∼1 × 10^{−6} cm s^{−1}. Near the top of the radiative zone (before tunneling through the evanescent region), we measure some mode velocities above 1 cm s^{−1} in the F5 case, corresponding to a luminosity perturbation of ∼3 ppm. In the F1a and F1b cases, considering the most excited modes, we find near the top of the radiative zone typical velocities of ∼1 × 10^{−3} cm s^{−1} which is more than one order of magnitude higher than the maximal value observed by A14 near the tachocline (∼5 × 10^{−5} cm s^{−1}). These differences between the solar model and our Ftype star models at 1 Ω_{⊙} (F1a and F1b) can be explained by considering that the power injection from convection to the waves find a kinetic energy flux proportional to (see e.g. Press 1981; Pinçon et al. 2016) where ρ_{b} is the density at the base of the convective zone and v_{b} is a characteristic velocity at the base of the convective zone. For a given rotation rate, we therefore expect the rms velocity of the modes, to be subject to the following scaling formula:
where the ⊙ index denotes the solar values. By comparing F1a and F1b with the solar model from A14, we find on one hand that convective motions at the bottom of the convective zone are significantly faster for our Ftype model, by a factor of approximately 20 for the rms velocity, and 10 for the maximal velocity of the downward flows. On the other hand, the density at the interface is ∼100 times lower in the Ftype case, thus meaning that the rms velocity of the modes should be ∼10 to 40 times larger, which is roughly consistent with what we see in the simulations. This considered, we therefore strongly advocate for the realisation of dedicated parametric studies of plume properties at the interface between the convective envelope and the radiative zone, including the effect of rotation. Such simulations would help to constrain the analytic models of IGWs excitation by convection, and would allow predicting the mode surface amplitude for a large range of models. Indeed, the inclusion of rotation in our work shows that a parameterisation from 3D simulations of the form:
with a an unknown scaling parameter, could be compared with the theoretical predictions from Augustson et al. (2020) and would prove extremely useful for future studies of IGWs in rotating stars.
We underline that some caveat must be kept in mind when considering these simulations. Due to the numerical constraints, we remind that the fluid regimes that are considered here are far from the actual regime in stellar interiors. The convective envelopes are in reality significantly more turbulent than what is currently achieved in any 2D or 3D simulation. In order to model a radiative zone as realistic as possible, we took κ and ν values that are five and four orders of magnitudes lower to the values used in the convective zone. In actual stellar interiors with Pr thought to be comprised between 1 × 10^{−5} and 1 × 10^{−9} (e.g. Garaud 2021), radiative diffusion dominates the damping of the waves as they travel through the radiative regions, in a quasiadiabatic regime (Zahn et al. 1997), which should allow modes to form below a much lower frequency cutoff than observed in our simulations. The values in the convective zone are constrained by L_{⋆}. Therefore, it is mainly the profile of the power transfer function from convection to the wave and the ability of the waves to propagate through the radiative interior to form standing modes that will determine the shape of the gmode power spectrum.
This work opens some observational perspectives for Ftype solartypepulsating stars. For our fastest rotating case, F5, we were able to detect the signature of ℓ = 3 modes at the top of the domain, a spherical harmonic degree which is still accessible to observations from diskintegrated instruments like the Kepler satellite, the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015), or the PLAnetary Transits and Oscillations of stars (PLATO, Rauer et al. 2014) satellite. For targets bright enough to consider this type of followup, ground observations with échelle spectrograph as the ones from the Stellar Observations Network Group (SONG, Grundahl et al. 2007) could also be of greatest interest in the perspective of improving the characterisation of lowfrequency regions of the brightest mainsequence Ftype stars with solartype oscillations observed by Kepler. As our simulations show that ℓ = 4, 5, 6, and 7 modes are also particularly excited, longterm observations of some Ftype stars from the Kepler LEGACY sample (Lund et al. 2017; Silva Aguirre et al. 2017), with a stellar imager instrument dedicated to asteroseismology represent an interesting perspective (e.g. ChristensenDalsgaard et al. 2011). The estimated luminosity perturbations induced by an individual mode remains small (below 1 ppm) from what we compute with the mode velocity near the top of the domain, especially with the fact that we deal with modes laying in frequency regions where the power contribution from the convective signal is important. However, this suggests that modes excited by a more efficient mechanism like tidal forcing (e.g. Fuller 2017) could reach the surface with an amplitude large enough to be detectable in these stars. In the future, we plan to include the convective core in the simulated domain in order to study the behaviour of IGWs excited simultaneously at the internal and the external interfaces of the radiative zone. The convective motions inside the core are also susceptible to play a significant role from a magnetohydrodynamic point of view. Indeed, in this work, we did not take the effect of the magnetic field into account. Lecoanet et al. (2022) recently showed that highorder gmodes can be suppressed through the action of a strong internal magnetic field generated by convective motions in the core.
Acknowledgments
The authors thank the referee for useful comments that helped improve the paper. S.N.B., A.S.B and R.A.G acknowledge the support from PLATO CNES grant. S.N.B. and A.S.B acknowledge the support from Solar Orbiter CNES grant and financial support by ERC Whole Sun Synergy grant #810218. S.N.B and R.A.G acknowledge the support from GOLF CNES grant. The authors acknowledge support from the GENCI project 1623. S.N.B reserves a special acknowledgment to L. Amard for his precious knowledge in stellar structure and evolution and is also grateful to Q. Noraz and A. Strugarek for advice relative to the ASH code and its outputs. He also thanks F. Grundahl for insights concerning the capabilities of the SONG telescopes. Finally, the authors thank K. Augustson, M. Delorme, A. Le Saux, S. Mathis, and C. Pinçon for fruitful discussions. Software: ASH (Clune et al. 1999; Brun et al. 2004), Python (Van Rossum & Drake 2009), numpy (Oliphant 2006; Harris et al. 2020), matplotlib (Hunter 2007), scipy (Virtanen et al. 2020), emcee (ForemanMackey et al. 2013).
References
 Aerts, C., Augustson, K., Mathis, S., et al. 2021, A&A, 656, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ahuir, J., Mathis, S., & Amard, L. 2021, A&A, 651, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Alvan, L., Brun, A. S., & Mathis, S. 2014, A&A, 565, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Alvan, L., Strugarek, A., Brun, A. S., Mathis, S., & Garcia, R. A. 2015, A&A, 581, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Amard, L., Palacios, A., Charbonnel, C., et al. 2019, A&A, 631, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 André, Q. 2019, Ph.D. Thesis, Université ParisCité, France [Google Scholar]
 Appourchaux, T., & Pallé, P. L. 2013, in Fifty Years of Seismology of the Sun and Stars, eds. K. Jain, S. C. Tripathy, F. Hill, J. W. Leibacher, & A. A. Pevtsov, ASP Conf. Ser., 478, 125 [NASA ADS] [Google Scholar]
 Augustson, K. C., Brown, B. P., Brun, A. S., Miesch, M. S., & Toomre, J. 2012, ApJ, 756, 169 [Google Scholar]
 Augustson, K. C., Brun, A. S., & Toomre, J. 2013, ApJ, 777, 153 [NASA ADS] [CrossRef] [Google Scholar]
 Augustson, K. C., Brun, A. S., & Toomre, J. 2016, ApJ, 829, 92 [Google Scholar]
 Augustson, K. C., Mathis, S., & Astoul, A. 2020, ApJ, 903, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Belkacem, K., Samadi, R., Goupil, M. J., et al. 2009, A&A, 494, 191 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Berger, T. A., Huber, D., van Saders, J. L., et al. 2020, AJ, 159, 280 [Google Scholar]
 Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977 [Google Scholar]
 Bowman, D. M. 2020, Front. Astron. Space Sci., 7, 70 [Google Scholar]
 Braginsky, S. I., & Roberts, P. H. 1995, Geophys. Astrophys. Fluid Dyn., 79, 1 [Google Scholar]
 Brown, B. P., Vasil, G. M., & Zweibel, E. G. 2012, ApJ, 756, 109 [Google Scholar]
 Browning, M. K., Brun, A. S., & Toomre, J. 2004, ApJ, 601, 512 [Google Scholar]
 Brummell, N. H., Clune, T. L., & Toomre, J. 2002, ApJ, 570, 825 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., Miesch, M. S., & Toomre, J. 2004, ApJ, 614, 1073 [Google Scholar]
 Brun, A. S., Miesch, M. S., & Toomre, J. 2011, ApJ, 742, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., Strugarek, A., Varela, J., et al. 2017, ApJ, 836, 192 [Google Scholar]
 Brun, A. S., Strugarek, A., Noraz, Q., et al. 2022, ApJ, 926, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Busse, F. H. 1970, J. Fluid Mech., 44, 441 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1961, Hydrodynamic and Hydromagnetic Stability (Oxford University Press) [Google Scholar]
 Chaplin, W. J., Kjeldsen, H., ChristensenDalsgaard, J., et al. 2011, Science, 332, 213 [Google Scholar]
 ChristensenDalsgaard, J. 2008, Lecture Notes on Stellar Oscillations, 5th edn. (Institut for Fysik og Astronomi, Aarhus Universitet) [Google Scholar]
 ChristensenDalsgaard, J., Carpenter, K. G., Schrijver, C. J., Karovska, M., & Si Team. 2011, J. Phys. Conf. Ser., 271, 012085 [NASA ADS] [CrossRef] [Google Scholar]
 Clune, T., Elliott, J., Miesch, M., Toomre, J., & Glatzmaier, G. 1999, Parallel Comput., 25, 361 [CrossRef] [Google Scholar]
 Corsaro, E., Bonanno, A., Mathur, S., et al. 2021, A&A, 652, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dhouib, H., Prat, V., Van Reeth, T., & Mathis, S. 2021a, A&A, 652, A154 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dhouib, H., Prat, V., Van Reeth, T., & Mathis, S. 2021b, A&A, 656, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dhouib, H., Mathis, S., Bugnet, L., Van Reeth, T., & Aerts, C. 2022, A&A, 661, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dupret, M. A., Grigahcène, A., Garrido, R., Gabriel, M., & Scuflaire, R. 2005, A&A, 435, 927 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Edelmann, P. V. F., Ratnasingam, R. P., Pedersen, M. G., et al. 2019, ApJ, 876, 4 [NASA ADS] [CrossRef] [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
 Fuller, J. 2017, MNRAS, 472, 1538 [Google Scholar]
 Garaud, P. 2021, Phys. Rev. Fluids, 6, 030501 [NASA ADS] [CrossRef] [Google Scholar]
 García, R. A., TurckChièze, S., JiménezReyes, S. J., et al. 2007, Science, 316, 1591 [Google Scholar]
 Garcia Lopez, R. J., & Spruit, H. C. 1991, ApJ, 377, 268 [NASA ADS] [CrossRef] [Google Scholar]
 Gastine, T., Yadav, R. K., Morin, J., Reiners, A., & Wicht, J. 2014, MNRAS, 438, L76 [CrossRef] [Google Scholar]
 Gilman, P. A., & Glatzmaier, G. A. 1981, ApJS, 45, 335 [CrossRef] [Google Scholar]
 Glatzmaier, G. A. 1984, J. Comput. Phys., 55, 461 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Kumar, P. 1990, ApJ, 363, 694 [NASA ADS] [CrossRef] [Google Scholar]
 Goldstein, J., & Townsend, R. H. D. 2020, ApJ, 899, 116 [NASA ADS] [CrossRef] [Google Scholar]
 Gough, D. O. 1993, in Astrophysical Fluid Dynamics  Les Houches 1987, 399 [Google Scholar]
 Grundahl, F., Kjeldsen, H., ChristensenDalsgaard, J., Arentoft, T., & Frandsen, S. 2007, Commun. Asteroseismol., 150, 300 [NASA ADS] [CrossRef] [Google Scholar]
 Guerrero, G., Zaire, B., Smolarkiewicz, P. K., et al. 2019, ApJ, 880, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Guzik, J. A., Kaye, A. B., Bradley, P. A., Cox, A. N., & Neuforge, C. 2000, ApJ, 542, L57 [Google Scholar]
 Hanasoge, S. M., Duvall, T. L., & Sreenivasan, K. R. 2012, Proc. Natl. Acad. Sci., 109, 11928 [Google Scholar]
 Hanasoge, S., Gizon, L., & Sreenivasan, K. R. 2016, Annu. Rev. Fluid Mech., 48, 191 [Google Scholar]
 Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
 Hotta, H., & Kusano, K. 2021, Nat. Astron., 5, 1100 [NASA ADS] [CrossRef] [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Hurlburt, N. E., Toomre, J., & Massaguer, J. M. 1986, ApJ, 311, 563 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, C. A., Kuzanyan, K. M., & Mitchell, R. H. 2009, J. Fluid Mech., 634, 291 [NASA ADS] [CrossRef] [Google Scholar]
 Kippenhahn, R., Weigert, A., & Weiss, A. 2012, Stellar Structure and Evolution (Springer) [Google Scholar]
 Lantz, S. R. 1992, Ph.D. Thesis, Cornell University, USA [Google Scholar]
 Le Saux, A., Guillet, T., Baraffe, I., et al. 2022, A&A, 660, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lecoanet, D., & Quataert, E. 2013, MNRAS, 430, 2363 [Google Scholar]
 Lecoanet, D., Le Bars, M., Burns, K. J., et al. 2015, Phys. Rev. E, 91 [CrossRef] [Google Scholar]
 Lecoanet, D., Bowman, D. M., & Van Reeth, T. 2022, MNRAS, 512, L16 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, U., & Saio, H. 1997, ApJ, 491, 839 [Google Scholar]
 Li, G., Van Reeth, T., Bedding, T. R., et al. 2020, MNRAS, 491, 3586 [Google Scholar]
 Lomb, N. R. 1976, Ap&SS, 39, 447 [Google Scholar]
 Lund, M. N., Silva Aguirre, V., Davies, G. R., et al. 2017, ApJ, 835, 172 [Google Scholar]
 Maeder, A. 2009, Physics, Formation and Evolution of Rotating Stars (Springer) [Google Scholar]
 Mathis, S., Neiner, C., & Tran Minh, N. 2014, A&A, 565, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathur, S., García, R. A., Breton, S., et al. 2022, A&A, 657, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Miesch, M. S., Elliott, J. R., Toomre, J., et al. 2000, ApJ, 532, 593 [CrossRef] [Google Scholar]
 Moravveji, E., Moya, A., & Guinan, E. F. 2012, ApJ, 749, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Neiner, C., Floquet, M., Samadi, R., et al. 2012, A&A, 546, A47 [CrossRef] [EDP Sciences] [Google Scholar]
 Neiner, C., Lee, U., Mathis, S., et al. 2020, A&A, 644, A9 [EDP Sciences] [Google Scholar]
 Noyes, R. W., Hartmann, L. W., Baliunas, S. L., Duncan, D. K., & Vaughan, A. H. 1984, ApJ, 279, 763 [Google Scholar]
 Oliphant, T. 2006, NumPy: A guide to NumPy (Trelgol Publishing) [Google Scholar]
 Ouazzani, R.M., Salmon, S. J. A. J., Antoci, V., et al. 2017, MNRAS, 465, 2294 [Google Scholar]
 Ouazzani, R. M., Marques, J. P., Goupil, M. J., et al. 2019, A&A, 626, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ouazzani, R. M., Lignières, F., Dupret, M. A., et al. 2020, A&A, 640, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pápics, P. I., Tkachenko, A., Van Reeth, T., et al. 2017, A&A, 598, A74 [Google Scholar]
 Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
 Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [Google Scholar]
 Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
 Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [Google Scholar]
 Pedersen, M. G., Aerts, C., Pápics, P. I., et al. 2021, Nat. Astron., 5, 715 [NASA ADS] [CrossRef] [Google Scholar]
 Perrard, S., Le Bars, M., & Le Gal, P. 2013, in Experimental and Numerical Investigation of Internal Gravity Waves Excited by Turbulent Penetrative Convection in Water Around Its Density Maximum, eds. M. Goupil, K. Belkacem, C. Neiner, F. Lignières, & J. J. Green, 865, 239 [NASA ADS] [Google Scholar]
 Pinçon, C., Belkacem, K., & Goupil, M. J. 2016, A&A, 588, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pratt, J., Baraffe, I., Goffrey, T., et al. 2017, A&A, 604, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Prat, V., Mathis, S., Buysschaert, B., et al. 2019, A&A, 627, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Prat, V., Mathis, S., Neiner, C., et al. 2020, A&A, 636, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Press, W. H. 1981, ApJ, 245, 286 [Google Scholar]
 Provost, J., & Berthomieu, G. 1986, A&A, 165, 218 [NASA ADS] [Google Scholar]
 Provost, J., Berthomieu, G., & Morel, P. 2000, A&A, 353, 775 [NASA ADS] [Google Scholar]
 Rauer, H., Catala, C., Aerts, C., et al. 2014, Exp. Astron., 38, 249 [Google Scholar]
 Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, J. Astron. Telesc. Instrum. Syst., 1, 014003 [Google Scholar]
 Rogers, T. M., & Glatzmaier, G. A. 2005, ApJ, 620, 432 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, T. M., Lin, D. N. C., McElwaine, J. N., & Lau, H. H. B. 2013, ApJ, 772, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Samadi, R., Belkacem, K., Goupil, M. J., et al. 2010, Ap&SS, 328, 253 [NASA ADS] [CrossRef] [Google Scholar]
 Santos, A. R. G., Breton, S. N., Mathur, S., & García, R. A. 2021, ApJS, 255, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Scargle, J. D. 1982, ApJ, 263, 835 [Google Scholar]
 Silva Aguirre, V., Lund, M. N., Antia, H. M., et al. 2017, ApJ, 835, 173 [Google Scholar]
 Szewczuk, W., Walczak, P., & DaszyńskaDaszkiewicz, J. 2021, MNRAS, 503, 5894 [NASA ADS] [CrossRef] [Google Scholar]
 Takehiro, S.I., Brun, A. S., & Yamada, M. 2020, ApJ, 893, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Talon, S., & Charbonnel, C. 2003, A&A, 405, 1025 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tassoul, M. 1980, ApJS, 43, 469 [Google Scholar]
 Townsend, R. H. D. 2003a, MNRAS, 343, 125 [Google Scholar]
 Townsend, R. H. D. 2003b, MNRAS, 340, 1020 [Google Scholar]
 Townsend, R. H. D., & Teitler, S. A. 2013, MNRAS, 435, 3406 [Google Scholar]
 Townsend, R. H. D., Goldstein, J., & Zweibel, E. G. 2018, MNRAS, 475, 879 [Google Scholar]
 Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial Oscillations of Stars (University of Tokyo Press) [Google Scholar]
 Van Reeth, T., Tkachenko, A., Aerts, C., et al. 2015, ApJS, 218, 27 [Google Scholar]
 Van Reeth, T., Tkachenko, A., & Aerts, C. 2016, A&A, 593, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Van Rossum, G., & Drake, F. L. 2009, Python 3 Reference Manual (Scotts Valley, CA: Create Space) [Google Scholar]
 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
 Warnecke, J., & Käpylä, M. J. 2020, A&A, 642, A66 [EDP Sciences] [Google Scholar]
 Zahn, J. P. 1991, A&A, 252, 179 [Google Scholar]
 Zahn, J. P. 1992, A&A, 265, 115 [NASA ADS] [Google Scholar]
 Zahn, J. P., Talon, S., & Matias, J. 1997, A&A, 322, 320 [NASA ADS] [Google Scholar]
Appendix A: Spherical harmonic expansion and power spectrum computation
This section provides some practical details on the techniques we used to obtain the power spectra E_{ℓ}. In practice, we only compute the development for m ≥ 0. We then consider the Fourier transform of the complex quantity v_{r}(r, ℓ, m, t), with −ω_{N} < ω < ω_{N}. From the spherical harmonic properties, we have . The power contribution of the retrogade modes, m < 0, therefore lies in −ω_{N} < ω < 0. The power contribution of the prograde modes, m > 0 is in 0 < ω < ω_{N}. The power contribution of a zonal mode, m = 0, at frequency ω, is split between ω and −ω. For ω ≥ 0, the power in each ℓ, m component is:
To compute the power spectrum E_{ℓ}(r, ℓ, ω) and restrict the domain to the positive frequencies ω ≥ 0, we consider for each degree ℓ
All Tables
Dimensionless numbers in the middle of the convective zone, τ_{conv}, and r_{c} − r_{0} penetration depth.
Asymptotic , ΔP_{ℓ} measured at r = 0.97 R_{⋆} and corresponding peak heights for the F5 case.
All Figures
Fig. 1. Radial profiles of the BruntVäisälä frequency N in the ASH model (black) and the MESA model (dashed grey). The structural and chemical contributions to N_{MESA}, N_{T, MESA}, and N_{μ, MESA}, are represented in dottedorange and dotteddarkblue lines, respectively. The dashed green vertical line shows the boundary between the radiative zone and the convective envelope, while the grey hatched regions are excluded from the simulation domain. The frequency corresponding to the ray paths of Fig. 15 are represented by the solid horizontal red, blue, and yellow lines. 

In the text 
Fig. 2. Density (blue) and temperature (orange) profiles used as reference profiles for the simulation are represented with straight lines. MESA model profiles are shown in dashed lines for comparison. 

In the text 
Fig. 3. Radial profile of the diffusivities, ν (straight lines) and κ (dashed lines), for the F1a (orange), F3 (blue), F1b, and F5 (red) cases. 

In the text 
Fig. 4. Distribution of surface rotation periods P_{rot} measured by Santos et al. (2021) in the Kepler data, considering stars with 6000 < T_{eff} < 6600 K and log g > 4 dex (using values from Berger et al. 2020). Vertical colour lines indicate the periods corresponding to the rotation rates selected for our simulations. 

In the text 
Fig. 5. Spacing in the nonuniform radial grid used in the 3D ASH simulations. 

In the text 
Fig. 6. Flux balance between radiative luminosity L_{rad} (dashed orange), kinetic energy luminosity L_{ke} (dashed red), enthalpy luminosity L_{en} (dashed grey), diffusivity luminosity L_{ν} (dashed green), and subgridscale eddy luminosity L_{ed} (dashed blue) for the F5 case. The total luminosity, L_{tot}, is represented in black. The top panel shows the flux balance for the whole simulation domain, while the bottom panel shows an enlargement in the convective zone and the overshoot region. 

In the text 
Fig. 7. 3D volume of the radial velocity, v_{r}, normalised by the rms velocity, (denoted v_rms in the colorbar) for the F5 case, from r = 0.07 to 0.94 R_{⋆}. Upward flows are in represented in yellow and downward flows in blue. 

In the text 
Fig. 8. Spherical harmonic decomposition of the time average of the rms radial velocity in the middle of the convective zone (r = 0.93 R_{⋆}) for the F1a (orange), F1b (cyan), F3 (blue), and F5 (red) cases. In the top panel, the rms components are summed over m and shown as a function of ℓ. In the bottom panel, the rms components are summed over ℓ and shown as a function of m. 

In the text 
Fig. 9. Radial velocity v_{r} (left) and temperature perturbation (right) at depth r = 0.95 R_{⋆} for the F1a, F1b, F3, and F5 cases (top to bottom). The latitudinal mean temperature perturbation, m = 0, is subtracted from T. 

In the text 
Fig. 10. profiles for the F1a (orange), F1b (cyan), F3 (blue), and F5 (red) cases. The vertical dashedgrey line correspond to the depth where we examine the mode signatures for the F5 case in Sect. 5. 

In the text 
Fig. 11. Fluid Rossby number Ro_{f} profiles in the convective zone for the F1a (orange), F1b (cyan), F3 (blue), and F5 (red) cases. 

In the text 
Fig. 12. Differential rotation for F1a (top left), F1b (top left), F3 (bottom left), and F5 (bottom right) cases. Retrograde and prograde (relatively to the corotational frame) flows are shown in blue and red, respectively. The dashed black line corresponds to r_{CZ} and the dashed grey lines are isocontours. 

In the text 
Fig. 13. Radial velocity v_{r} (left) and temperature perturbation (right) at depth r = 0.86 R_{⋆} for the F1a case. The latitudinal mean temperature perturbation was substracted from T. 

In the text 
Fig. 14. Overshoot region (blue) and bottom of the convective zone (red) in the F1a (top left), F1b (top right), F3 (bottom left), and F5 (bottom right) models, for θ spanning from 30 to 150 ^{o}. The white line show the radius r_{c} at which the enthalpy flux crosses zero at a given latitude while the orange line signals the radius r_{0} where the local enthalpy flux is equal to a tenth of the maximal enthalpy flux at the corresponding latitude. 

In the text 
Fig. 15. Propagation path computed with raytracing for three IGWs: ℓ = 1, ω = 10 μHz (top), ℓ = 1, ω = 140 μHz (middle), ℓ = 2, ω = 319 μHz (bottom). The integration time, Δt, considered for each ray is specified. From the surface to the center, the black circles represent the position of the stellar surface, the bottom of the convective envelope and the edge of the convective core. The grey area corresponds to the part of the radiative zone where the wave is evanescent. 

In the text 
Fig. 16. Equatorial cut for for the radiative zone in the F1a case, with application of a temporal filtering centered at 140 μHz. The 2D propagation ray obtained by numerically solving the system of Eq. (17) and shown in the middle panel of Fig. 15 is overplotted in white. We note that there is a good qualitative agreement. 

In the text 
Fig. 17. E_{ℓ} power spectrum at r = 0.18 R_{⋆} for the F1a (top left), F1b (top right), F3 (bottom left), and F5 (bottom) case. The mode frequencies computed with GYRE are overplotted (black dots) for ℓ = 1 to ℓ = 40. In each panel, the white line highlights the separation between the gmode region (top left) and the progressive waves region (bottom right). 

In the text 
Fig. 18. Amount of power injected in IGW. Top: total power E_{ℓ, tot} for 1 ≤ ℓ ≤ 10 at r = 0.8 R_{⋆} for the F1a (orange diamonds), F1b (cyan circles) F3 (blue crosses), and F5 (red pentagons) cases. Bottom: ratio between E_{ℓ, tot}(r = 0.8 R_{⋆}), the radiative zone, and E_{ℓ, tot}(r = 0.9 R_{⋆}), convective zone. The same symbols and colour coding are used as in previous figures. 

In the text 
Fig. 19. GYRE ξ_{r} eigenfunctions (blue) compared with ASH outputs for ℓ = 5, n = 7 and ℓ = 5, n = 12. The depths where ASH outputs are sampled are signaled by orange squares. 

In the text 
Fig. 20. Top panels: power spectra E_{1} (0.38 R_{⋆}, ω). Middle panels: power spectra E_{1} (0.38 R_{⋆}, P). Bottom panels: LombScargle periodograms computed from E_{1} (0.38 R_{⋆}, P). Each periodogram is normalised with its standard deviation, σ and the obtained ΔP_{l} value is shown by a vertical dashed blue line. From left to right, the represented cases are F1a, F1b, F3, and F5. 

In the text 
Fig. 21. LombScargle periodograms computed from E_{2} (0.38 R_{⋆}, P). Each periodogram is normalised with its standard deviation, σ and the obtained ΔP_{2} value is shown by a vertical dashed red line. From left to right, the represented cases are F1a, F1b, F3, and F5. 

In the text 
Fig. 22. ΔP_{1} as a function of modes period in the F5 case, for azimuthal numbers m = −1 (yellow), m = 0 (blue), and m = 1 (red). The ΔP_{1} obtained from the Ω_{⋆} = 0 GYRE run and the perturbative approach are represented as grey hexagons while the white crosses shows the ΔP_{1} obtained for m = { − 1, 0, 1} with the 5 Ω_{⊙} GYRE TAR run. The black crosses signal the ΔP_{1} position for m = { − 1, 1} in the asymptotic approximation of the perturbative method. The measured uncertainties on period and ΔP_{1} are represented. We show only fitted modes for which we are able to measure ΔP_{1} with an uncertainty below 1.75 min. 

In the text 
Fig. 23. Example of rotational splitting observed in the F5 case for a ℓ = 2 mode, with the m = −2 shown in brown, the m = −1 in yellow, the m = 0 in blue, the m = 1 in red, and the m = 2 in orange. 

In the text 
Fig. 24. Rotation rate Ω inferred from the asymptotic approximation of Eqs. (27) and (29) with δ_{n, 1, m} measurements, in the F5 case. The position of the true rotation value of the model is highlighted by the dashed grey line, with the 10% error interval in grey. Rotation rate measured only with δ_{n, 1, −1} and δ_{n, 1, 1} are represented in yellow and red, respectively, while the averaged mean of the two values is shown in blue. The error bar for each Ω measurement is also represented. 

In the text 
Fig. 25. Power spectrum for ℓ = 1, E_{1}, radial evolution in the F1b (left) and F5 (right) cases. In each panel, the white line represents the boundary of the mode resonant cavity. The range of radial orders visible in the figure is 1 ≤ n ≤ 30. 

In the text 
Fig. 26. Power spectrum for ℓ = 4, E_{4}, radial evolution in the F1b (left) and F5 (right) cases. In each panel, the white line represents the boundary of the mode resonant cavity. The cutoff between gmodes and progressive waves is visible around 800 min. Below this cutoff, the range of visible radial orders is 1 ≤ n ≤ 50. 

In the text 
Fig. 27. LombScargle periodograms computed from E_{3} (0.97 R_{⋆}, P) (top) and E_{5} (0.97 R_{⋆}, P) (bottom). The periodogram is normalised with its standard deviation σ and the obtained ΔP_{3} value is shown by a vertical dashed red line. 

In the text 
Fig. 28. Luminosity perturbation estimation as a function of v_{r} for a ℓ = 3 mode at 41.5 μHz (red), a ℓ = 4 mode at 41.5 μHz (blue), and a ℓ = 5 mode at 39.4 μHz (yellow). The diamonds correspond to the radial velocity we measure in the F5 case at r = 0.87 R_{⋆} and the dots to the same quantity at r = 0.97 R_{⋆}. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.