© S. N. Breton et al. 2022

1. Introduction

The propagation of internal gravity waves (IGW) in stellar interiors is an expected and well-known phenomenon (see e.g. Christensen-Dalsgaard 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 (g-modes). 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 intermediate-mass and massive stars, the κ-mechanism that is related to opacity-bump in ionization regions is also able to drive mode excitation (e.g. Unno et al. 1989; Bowman 2020); whereas for early F-type 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 low-mass 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 g-mode 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 main-sequence solar-type pulsating stars (e.g. Chaplin et al. 2011; Mathur et al. 2022), it is interesting to consider the case of solar-type pulsators other than the Sun.

Among main-sequence solar-type pulsators, F-type stars are probably the most promising for g-mode 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 g-modes towards the surface. Most of observed F-type 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 F-type 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 gravito-inertial 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 intermediate-mass and massive stars at the interface between the convective core and the radiative interior. Stochastically excited gravito-inertial modes have also been detected in rapidly rotating massive stars (Neiner et al. 2012, 2020). Due to the difficulties in observing g-modes in solar-type pulsators discussed above, it appears that performing 3D deep-shell simulations of late F-type stars constitutes a first approach to understanding the interplay between IGWs and rotation in their radiative interiors. While 3D simulations of F-type stars have already been performed, they have only been aimed at studying the convective-envelope 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 solar-model simulations to evaluate the possibility to detect individual g-modes 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 core-radiative 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 low-mass stars to non-solar models with varying rotation rates. We present the first deep-shell hydrodynamical simulations of rotating F-type 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 F-type star model for the different rotation regimes we considered. The IGWs and g-mode 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 Ω₀ = Ω₀e_z. The reference density, pressure, temperature, and specific entropy are denoted as $\overline{\rho }$, $\overline{P}$, $\overline{T}$, and $\overline{S}$, 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 Lantz-Braginsky-Roberts (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 $\stackrel{\sim }{\omega }=P/\overline{\rho }$ is considered instead of the fluctuating pressure P. The non-linear momentum equation is therefore:

$$\begin{array}{c}\hfill \overline{\rho }(\frac{\partial \mathit{v}}{\partial t}+(\mathit{v}.\mathrm{\nabla })\mathit{v})=-\overline{\rho }\mathrm{\nabla }\stackrel{\sim }{\omega }-\overline{\rho }\frac{S}{c_p}\mathit{g}-2\overline{\rho }{\mathbf{\Omega }}_0\times \mathit{v}-\mathrm{\nabla }.\mathcal{D},\end{array}$$(1)

where v = (v_r, v_θ, v_ϕ) is the local velocity, g is the gravitational acceleration, and 𝒟 is the viscous stress tensor:

$$\begin{array}{c}\hfill {\mathcal{D}}_{\mathit{ij}}=-2\overline{\rho }v(e_{\mathit{ij}}-\frac{1}{3}(\mathrm{\nabla }.\mathit{v}){\delta }_{\mathit{ij}}),\end{array}$$(2)

with e_ij = 1/2(∂_jv_i+∂_iv_j) the strain rate tensor and δ_ij the Kronecker symbol. In the anelastic approximation, the continuity equation is expressed as:

$$\begin{array}{c}\hfill \mathrm{\nabla }.(\overline{\rho }\mathit{v})=0.\end{array}$$(3)

We assume a linearised equation of state and the zeroth-order ideal gas law:

$$\begin{array}{c}\hfill \frac{\rho }{\overline{\rho }}=\frac{P}{\overline{P}}-\frac{T}{\overline{T}}=\frac{P}{\gamma \overline{P}}-\frac{S}{c_p},\end{array}$$(4)

$$\begin{array}{c}\hfill \overline{P}=\mathcal{R}\overline{\rho }\overline{T},\end{array}$$(5)

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:

$$\begin{array}{cc}\hfill \overline{\rho }\overline{T}\frac{\partial S}{\partial t}+\overline{\rho }\overline{T}\mathit{v}.\mathrm{\nabla }(S+\overline{S})& =\overline{\rho }ϵ+\mathrm{\nabla }.[{\kappa }_r\overline{\rho }{c}_p\mathrm{\nabla }(T+\overline{T})+\kappa \overline{\rho }\overline{T}\mathrm{\nabla }S\hfill \\ \hfill & +{\kappa }_0\overline{\rho }\overline{T}\mathrm{\nabla }\overline{S}]+2\overline{\rho }\nu [e_{\mathit{ij}}{e}_{\mathit{ij}}-\frac{1}{3}{(\mathrm{\nabla }.\mathit{v})}^2],\hfill \end{array}$$(6)

where κ_r is the radiative diffusivity and $\overline{\rho }ϵ$ the volume-heating term related to the energy generation by nuclear burning. The ϵ profile is parametrised as a power law of $\overline{T}$, $ϵ={ϵ}_0{\overline{T}}^k$. The parameters ϵ₀ 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 heat-transport by subgrid-scale (SGS) motions, while the diffusivity, κ₀, 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 κ₀ 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¹¹ 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:

1. rigid: v_r|_rtop = v_r|_rbottom = 0 ;

2. stress-free: $\frac{\partial }{\partial r}\left(\frac{v_{\theta }}{r}\right){|}_{r_{\mathrm{top}},r_{\mathrm{bottom}}}=\frac{\partial }{\partial r}\left(\frac{v_{\varphi }}{r}\right){|}_{r_{\mathrm{top}},r_{\mathrm{bottom}}}=0$ ;

3. constant mean entropy gradient:

   $$\begin{array}{cc}\hfill \frac{\partial \overline{S}}{\partial r}{|}_{r_{\mathrm{bottom}}}& =0,\hfill \\ \hfill \frac{\partial \overline{S}}{\partial r}{|}_{r_{\mathrm{top}}}& =-1.5\times {10}^{-6}\mathrm{cm}{\mathrm{K}}^{-1}{\mathrm{s}}^{-2}.\hfill \end{array}$$

The interface between the radiative zone and the upper convective zone is located at r_CZ = 0.88 R_⋆ (9.18 × 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 ϵ₀ = 3.23 × 10⁻⁷ 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 Brunt-Väisäla frequency profile N (given in Hz) according to the relation (e.g. Maeder 2009):

$$\begin{array}{c}\hfill N^2(r)=\frac{1}{2\pi }\frac{g}{c_p}\frac{\mathrm{d}S}{\mathrm{d}r},\end{array}$$(7)

where we take a uniform c_p = 3.42 × 10⁸ erg g⁻¹ K⁻¹ for our input ASH profile. We know that IGWs are propagative in regions where N² > 0 and evanescent in regions where N² < 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 Brunt-Vä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 ${\int }_{r_{\text{bottom}}}^{r_{\text{top}}}N/r\text{d}r$, which is directly related to the g-mode period spacing, we find a 11.6% discrepancy, which means that the g-mode 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 g-mode frequencies predicted by a 1D oscillation code using the 1D ASH reference profile and Sect. 4.3. The $\overline{\rho }$ and $\overline{T}$ are computed using a Newton-Raphson algorithm in order to set the reference state at the hydrostatic equilibrium. In Fig. 2, we compare the profiles obtained with the Newton-Raphson 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 Brunt-Väisälä frequency N in the ASH model (black) and the MESA model (dashed grey). The structural and chemical contributions to NMESA, NT, MESA, and Nμ, MESA, are represented in dotted-orange and dotted-dark-blue 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 $\overline{\rho }$ (blue) and temperature $\overline{T}$ (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 Ω₀ = 1 (F1a and F1b cases), 3 (F3 case), and 5 Ω_⊙ (F5 case), where we take Ω_⊙ = 2.6 × 10⁻⁶ rad s⁻¹. All the considered cases verify that ${\mathrm{\Omega }}_0\ll {\mathrm{\Omega }}_K=\sqrt{\mathrm{G}{M}_{\star }/R_{\star }^3}\sim 3.9\times {10}^{-4}$ rad s⁻¹, 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 slow-rotator 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 Prot measured by Santos et al. (2021) in the Kepler data, considering stars with 6000 < Teff < 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:

$$\begin{array}{cc}\hfill \nu (r)& ={\nu }_{\mathrm{top}}[{\beta }_{\nu }+\frac{1-{\beta }_{\nu }}{2}{\left(\frac{{\overline{\rho }}_{\mathrm{top}}}{\overline{\rho }}\right)}^{\frac{1}{2}}(tanh\frac{r-r_t}{{\sigma }_t}+1)],\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill \kappa (r)& ={\kappa }_{\mathrm{top}}[{\beta }_{\kappa }+\frac{1-{\beta }_{\kappa }}{2}{\left(\frac{{\overline{\rho }}_{\mathrm{top}}}{\overline{\rho }}\right)}^{\frac{1}{2}}(tanh\frac{r-r_t}{{\sigma }_t}+1)],\hfill \end{array}$$(9)

with ${\overline{\rho }}_{\mathrm{top}}=\overline{\rho }(r_{\mathrm{top}})$, ${\overline{\nu }}_{\mathrm{top}}=\overline{\nu }(r_{\mathrm{top}})$, ${\overline{\kappa }}_{\mathrm{top}}=\overline{\kappa }(r_{\mathrm{top}})$, σ_t = 0.08 the profile stiffness parameter, β_ν = 10⁻⁴, and β_κ = 2 × 10⁻⁵. 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.

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 non-uniform-grid 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 Adams-Bashforth time integration scheme for the advection and Coriolis terms, while the diffusive and buoyancy terms are treated through a semi-implicit Crank-Nicolson method (Glatzmaier 1984; Clune et al. 1999).

Fig. 5. Spacing in the non-uniform 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

$$\begin{array}{c}\hfill L_{\mathrm{tot}}=L_{\mathrm{rad}}+L_{\mathrm{ke}}+L_{\nu }+L_{\mathrm{en}}+L_{\mathrm{ed}},\end{array}$$(10)

Fig. 6. Flux balance between radiative luminosity Lrad (dashed orange), kinetic energy luminosity Lke (dashed red), enthalpy luminosity Len (dashed grey), diffusivity luminosity Lν (dashed green), and subgrid-scale eddy luminosity Led (dashed blue) for the F5 case. The total luminosity, Ltot, 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⁵ 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. F-type 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 root-mean square (rms) value at each given radius, ${\stackrel{\sim }{v}}_r$. 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 low-frequency 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, vr, normalised by the rms velocity, ${\stackrel{\sim }{v}}_r$ (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:

$$\begin{array}{c}\hfill Re=\frac{{\stackrel{\sim }{v}}_{r}L}{\nu },\end{array}$$(11)

where ${\stackrel{\sim }{v}}_r$ is the root-mean-square (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 $Ta=4{\mathrm{\Omega }}_0^2{L}^4/{\nu }^2$ is the Taylor number. The diffusivity profiles, ν and κ, should scale as ${\mathrm{\Omega }}_0^{-2}$ 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 ${\mathrm{\Omega }}_0^{-0.63}$ between the F1a and F3 cases and ${\mathrm{\Omega }}_0^{-0.75}$ 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⁶ and Ra_c = 4 × 10⁵). 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, ${\stackrel{\sim }{v}}_r$, 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 viscous-dominated 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 m-spectrum 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 so-called 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.

Figure 9 shows the impact of rotation on the convective structures, represented through radial velocity, v_r, and temperature perturbation $T-\overline{T}$ maps at r = 0.95 R_⋆. At this depth, v_r, and $T-\overline{T}$ 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 Ω₀ yields sharper patterns for the velocity field. Banana-shape 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, large-scale structures with large temperature perturbations (positive or negative) appear at high latitude. As it can be seen in Fig. 10, the ${\stackrel{\sim }{v}}_r$ radial profile is similar in the convective zone for the four cases. In the middle of the convective zone, we have ${\stackrel{\sim }{v}}_r\approx 3.4\times {10}^4$, 3.7 × 10⁴, 3.3 × 10⁴, and 3.2 × 10⁴ cm s⁻¹ 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 ${\stackrel{\sim }{v}}_r$ value in the radiative zone is associated with the amplitude of IGWs propagating in the stellar interior. Larger ${\stackrel{\sim }{v}}_r$ 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 ${\stackrel{\sim }{v}}_r$ stay at a level equivalent to the F3 case until r/R_⋆ ≈ 0.6. At the bottom of the radiative zone, however, ${\stackrel{\sim }{v}}_r$ values from the F3 case are comparable to the F1b case.

Fig. 9. Radial velocity vr (left) and temperature perturbation $T-\overline{T}$ (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. ${\stackrel{\sim }{v}}_r$ profiles for the F1a (orange), F1b (cyan), F3 (blue), and F5 (red) cases. The vertical dashed-grey 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:

$$\begin{array}{c}\hfill R{o}_f=\frac{|\mathrm{\nabla }\times \mathit{v}|}{2{\mathrm{\Omega }}_0},\end{array}$$(12)

where ∇ × v is the vorticity of the flow. Brun et al. (2017) expect a transition between the solar (fast equator, slow poles) and anti-solar 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 Taylor-Proudman 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

$$\begin{array}{c}\hfill R{o}_s=\frac{P_{\mathrm{rot}}}{{\tau }_{\mathrm{conv}}},\end{array}$$(13)

where the rotation period is P_rot = 2π/Ω₀, 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 ${\stackrel{\sim }{v}}_r$ in the convective zone, ${⟨{\stackrel{\sim }{v}}_r⟩}_{\mathrm{CZ}}$. For Ro_c, we use the following definition

$$\begin{array}{c}\hfill R{o}_c=\sqrt{\frac{\mathit{Ra}}{\mathit{TaPr}}},\end{array}$$(14)

where Ra = ( − ∂ρ/∂S)(∂S_tot/∂r)gL⁴/ρκν, with $S_{\mathrm{tot}}=\overline{S}+S$. 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 anti-solar, with a Ro_f of 1.9 in the middle of the convective zone. Interestingly, the transition from prograde to retrograde flows (in the co-rotational 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 anti-solar 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 F-star models published by Augustson et al. (2012) confirm the change of rotation regime towards Taylor-Proudman 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 Rof 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 co-rotational frame) flows are shown in blue and red, respectively. The dashed black line corresponds to rCZ 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 $T-\overline{T}$ change in the overshoot region relatively to the convective zone. We show the v_r and $T-\overline{T}$ maps at depth r = 0.86 R_⋆ for the F1a case. In the convective zone, v_r and $T-\overline{T}$ 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 anti-correlated: downwards travelling plumes penetrating the top of the stably-stratified 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 vr (left) and temperature perturbation $T-\overline{T}$ (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₀ 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₀ 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₀ at latitudes below 55^o for each case, we obtain the mean penetration depth of the overshooting flows, ${\overline{d}}_{\mathrm{ov}}$, in the same way as Brun et al. (2017). We obtain ${\overline{d}}_{\mathrm{ov}}=0.22$, 0.23, 0.19 and 0.15 H_p where H_p ≈ 8.8 × 10⁹ cm is the pressure scale height at the base of the convective zone. The comparison between F1a and F1b shows that, at constant Ω₀, the overshoot depth increases with Re but, as expected, the penetration depth of the overshooting flow is significantly reduced as rotation increases. The ${\overline{d}}_{\mathrm{ov}}$ 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 rc at which the enthalpy flux crosses zero at a given latitude while the orange line signals the radius r0 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):

$$\begin{array}{c}\hfill {\omega }^2=N^2\frac{k_h^2}{k^2},\end{array}$$(15)

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:

$$\begin{array}{c}\hfill k_h=\frac{\sqrt{\ell (\ell +1)}}{r}.\end{array}$$(16)

Concerning gravito-inertial waves, when the Coriolis parameter f = 2Ω₀ is such that f ≪ N, the effect of the Coriolis acceleration on the wave behaviour is negligible and we retrieve the pure-gravity case. Hence, we expect high-frequency waves propagation not to be affected by the rotation regime of the different cases. However, as we increase Ω₀, 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 F-type model, we start by computing the propagation path of IGWs at different frequencies following a ray-tracing 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 ray-tracing, we must numerically solve the following set of equations:

$$\begin{array}{c}\hfill \{\begin{array}{cc}\hfill \frac{\mathrm{d}r}{\mathrm{d}t}& =\frac{\partial W}{\partial r},\hfill \\ \hfill \frac{\mathrm{d}\varphi }{\mathit{vt}}& =\frac{1}{r}\frac{\partial W}{\partial \varphi },\hfill \\ \hfill \frac{\mathrm{d}{k}_r}{\mathrm{d}t}& =-\frac{\partial W}{\partial k_r}+\frac{1}{r}\frac{\partial W}{\partial k_{\varphi }}{k}_{\varphi },\hfill \\ \hfill \frac{\mathrm{d}{k}_{\varphi }}{\mathrm{d}t}& =-\frac{1}{r}\frac{\partial W}{\partial k_{\varphi }}-\frac{1}{r}\frac{\partial W}{\partial k_r}{k}_{\varphi },\hfill \\ \hfill \end{array}\end{array}$$(17)

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, high-frequency IGWs are trapped in the deepest regions of the radiative zone. High-frequency IGWs propagate much faster in the radiative interior. Low-frequency IGWs have a characteristic spiralling trajectory and take more time to reach the edge of the convective core.

Fig. 15. Propagation path computed with ray-tracing 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 $v_r/{\stackrel{\sim }{v}}_r$, 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 ray-tracing. 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 $v_r/{\stackrel{\sim }{v}}_r$ 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 ${\stackrel{\sim }{v}}_r(r,\ell ,m,\omega )$ of v_r, we then obtain the power spectrum E_ℓ, m as

$$\begin{array}{c}\hfill E_{\ell ,m}(r,\omega )={|\stackrel{\sim }{v_r}(r,\ell ,m,\omega )|}^2,\end{array}$$(18)

and, by summing over the m for each ℓ, we have the power spectrum, E_ℓ,

$$\begin{array}{c}\hfill E_{\ell }(r,\omega )={\displaystyle \sum _{m=-\ell }^{\ell }}{E}_{\ell ,m}(r,\omega ).\end{array}$$(19)

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 non-adiabatic 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, g-modes 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

$$\begin{array}{c}\hfill {\omega }_c(\ell )={[\ell (\ell +1)]}^{3/8}{{\displaystyle \left[{\int }_{r_{\mathrm{bottom}}}^{r_{\mathrm{CZ}}}\kappa \frac{N^3}{{\tau }_c{r}^3}dr\right]}}^{1/4},\end{array}$$(20)

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 g-mode 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 ${\stackrel{\sim }{v}}_r$ in the radiative interior, the IGWs mean amplitude is significantly larger in F5 and F3 than F1a. In particular, low-frequency-mode amplitudes are significantly larger in the F3 and F5 cases. Being a more turbulent version of F1a, F1b exhibits high-frequency modes of larger amplitude than in any other case, but the excitation of low-frequency 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,

$$\begin{array}{c}\hfill E_{\ell ,\mathrm{tot}}(r)={\displaystyle \sum _{{\omega }_i}}{E}_{\ell }(r,{\omega }_i).\end{array}$$(21)

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. g-modes 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 high-order g-modes is (Tassoul 1980; Provost & Berthomieu 1986):

$$\begin{array}{c}\hfill {\overline{\mathrm{\Delta }P}}_{\ell }=\frac{\pi }{\sqrt{\ell (\ell +1)}{\int }_{r_1}^{r_2}\frac{N}{r}\mathrm{d}r},\end{array}$$(22)

where, for a mode of frequency ω, r₁ and r₂ are the inner and outer boundary of the corresponding resonant cavity, respectively, defined by N(r₁) = N(r₂) = ω.

In Fig. 20, for each case, we represent the E₁ (0.38 R_⋆, ω) power spectra and the corresponding E₁ (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 Lomb-Scargle periodogram (Lomb 1976; Scargle 1982) of E_ℓ (0.38 R_⋆, P) for periods in the Lomb-Scargle periodogram spanning from 5 to 60 min. Each Lomb-Scargle periodograms is normalised by its standard deviation σ. We perform the same operation for the E₂ power spectra. The obtained Lomb-Scargle periodograms are shown in the bottom panels of Fig. 20 for E₁ and in Fig. 21 for E₂. We use the periodograms to identify the period spacings, ΔP₁ and ΔP₂, which we compare to the asymptotic value, ${\overline{\mathrm{\Delta }P}}_1=47$ min and ${\overline{\mathrm{\Delta }P}}_2=27.1$ min. We find ΔP₁ = 45.4, 42.3, 45.4, and 45.8 min and ΔP₁ = 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 low-order 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 Lomb-Scargle 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 E1 (0.38 R⋆, ω). Middle panels: power spectra E1 (0.38 R⋆, P). Bottom panels: Lomb-Scargle periodograms computed from E1 (0.38 R⋆, P). Each periodogram is normalised with its standard deviation, σ and the obtained ΔPl 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 low-frequency modes between slow-rotating 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 high-frequency 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. Lomb-Scargle periodograms computed from E2 (0.38 R⋆, P). Each periodogram is normalised with its standard deviation, σ and the obtained ΔP2 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 well-suited 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 280-day 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 χ² statistics with two degrees of freedom, implemented with the emcee module (Foreman-Mackey 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:

$$\begin{array}{c}\hfill \mathrm{\Delta }{P}_{n,1,m}=P_{n+1,1,m}-P_{n,1,m}.\end{array}$$(23)

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 least-square distance metric:

$$\begin{array}{c}\hfill \mathcal{M}={\displaystyle \sum _{n,m}}{|\mathrm{\Delta }{P}_{n,1,m(\mathrm{ASH})}-\mathrm{\Delta }{P}_{n,1,m(\mathrm{GYRE})}|}^2/{\sigma }_{n,1,m}^2,\end{array}$$(24)

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

where the ${\sigma }_{n,1,m}^2$ are the uncertainties measured for each value of ΔP_{n, 1, m (ASH)}. As a least-square distance metric, ℳ is equivalent to a log-likelihood 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Ω₀, the rotational splitting for a mode component with frequency ω_nℓm is given by (e.g. Christensen-Dalsgaard 2008)

$$\begin{array}{c}\hfill {\omega }_{n\ell m}={\omega }_{n\ell }+{\delta }_{n\ell m}\end{array}$$(25)

with, in the inertial frame,

$$\begin{array}{c}\hfill {\delta }_{n\ell m}=m{\beta }_{n\ell }{\mathrm{\Omega }}_0,\end{array}$$(26)

and in the co-rotating frame,

$$\begin{array}{c}\hfill {\delta }_{n\ell m}=-m(1-{\beta }_{n\ell })\mathrm{\Omega }.\end{array}$$(27)

The parameter β_nℓ is:

$$\begin{array}{c}\hfill {\beta }_{n\ell }=\frac{{\int }_0^{R_{\star }}({\xi }_r^2+{\mathcal{L}}^2{\xi }_h^2-2{\xi }_r{\xi }_h-{\xi }_h^2)r^2\overline{\rho }\mathrm{d}r}{{\int }_0^{R_{\star }}({\xi }_r^2+{\mathcal{L}}^2{\xi }_h^2-2{\xi }_r{\xi }_h)r^2\overline{\rho }\mathrm{d}r},\end{array}$$(28)

where ξ_r and ξ_h are the radial and horizontal displacements, respectively, and $\mathcal{L}=\sqrt{\ell (\ell +1)}$. In the case of high-order g-modes, we can consider

$$\begin{array}{c}\hfill {\beta }_{n\ell }\approx 1-\frac{1}{\ell (\ell +1)}.\end{array}$$(29)

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 co-rotating frame, as expected from Eq. (27), the m-component 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 Ω₀. 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 low-n 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 high-period (low-frequency) 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 co-rotational frame, mode splittings appear less clearly for ℓ = 4 modes. For the E₄ power spectrum, the separation between g-modes 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, E1, 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, E4, 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 g-modes 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 mid-degree modes tunneling through the convective zone and therefore we look for g-mode signatures in the convective velocity signal near the top of the domain, where we still have ${\stackrel{\sim }{v}}_r$ of the order of 1.5 × 10⁴ cm s⁻¹ (see Fig. 10). The F5 case is the only model where clear g-mode 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 g-mode signature of ℓ = 3, 4, 5, 6, and 7 modes. This is particularly interesting for ℓ = 3 modes which are still observable in disk-integrated observations. As an illustration, we show in Fig. 27 the Lomb-Scargle periodograms computed for ℓ = 3 and ℓ = 5 modes, for the F5 case. The clearest detection is the ΔP₅, with a peak height at 7.1 σ against 4.4 σ for ΔP₃, and 5.8 for ΔP₄. We report in Table 3 the measured ΔP_ℓ and corresponding peak heights. They are compared with the asymptotic ${\overline{\mathrm{\Delta }P}}_{\ell }$ computed with Eq. (22).

Fig. 27. Lomb-Scargle periodograms computed from E3 (0.97 R⋆, P) (top) and E5 (0.97 R⋆, P) (bottom). The periodogram is normalised with its standard deviation σ and the obtained ΔP3 value is shown by a vertical dashed red line.

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 Stefan-Boltzmann law we have (e.g. Samadi et al. 2010)

$$\begin{array}{c}\hfill \frac{\delta L_{\star }}{L_{\star }}=4\frac{\delta T_{\mathrm{eff}}}{T_{\mathrm{eff}}}+2\frac{\delta R_{\star }}{R_{\star }},\end{array}$$(30)

and, following Townsend (2003a)

$$\begin{array}{c}\hfill \frac{\delta T_{\mathrm{eff}}}{T_{\mathrm{eff}}}={\mathrm{\nabla }}_{\mathrm{ad}}|\frac{\ell (\ell +1)}{{\overline{\omega }}^2}-4-{\overline{\omega }}^2|\frac{\delta R_{\star }}{R_{\star }},\end{array}$$(31)

where ∇_ad = 0.4 is the adiabatic temperature gradient, and ${\overline{\omega }}^2$ is a dimensionless frequency defined as

$$\begin{array}{c}\hfill {\overline{\omega }}^2=\frac{4{\pi }^2{\omega }^2{R}_{\star }^3}{G{M}_{\star }}.\end{array}$$(32)

To compute the displacement δR_⋆ related to a single mode, we simply consider:

$$\begin{array}{c}\hfill \delta R_{\star }=\frac{{\stackrel{\sim }{v}}_{r,\mathrm{mode}}}{2\pi \omega },\end{array}$$(33)

where ${\stackrel{\sim }{v}}_{r,\mathrm{mode}}$ 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 ${\stackrel{\sim }{v}}_{r,\mathrm{mode}}$ 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 part-per-millions (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 Brunt-Väisala resonant cavity, gravito-inertial waves related to the modes would have to be sub-inertial (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 super-inertial 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 solar-type p-mode 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 g-mode 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 gravito-inertial modes in late F-type stars.

Fig. 28. Luminosity perturbation estimation as a function of vr 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 deep-shell 3D simulation of IGWs excited in the radiative interior of a F-type solar-type 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 F-type 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 low-frequency 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 Lomb-Scargle periodograms of the power spectra to detect this periodicity and we compared the measured values with asymptotic predictions. Studying the period-spacing evolution at long periods in the F5 case, we were also able to observe the effect of the Coriolis force on the behaviour of high-order modes. Finally, we tried to quantify the possibility of observing g-mode signatures in actual F-type solar-type 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⁻¹ cm s⁻¹, corresponding to a bolometric luminosity perturbation of ∼10⁻¹ 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⁻⁶ cm s⁻¹. Near the top of the radiative zone (before tunneling through the evanescent region), we measure some mode velocities above 1 cm s⁻¹ 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⁻³ cm s⁻¹ which is more than one order of magnitude higher than the maximal value observed by A14 near the tachocline (∼5 × 10⁻⁵ cm s⁻¹). These differences between the solar model and our F-type 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 ${\rho }_b{v}_b^4$ (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, ${\stackrel{\sim }{v}}_{\mathrm{mode}}$ to be subject to the following scaling formula:

$$\begin{array}{c}\hfill {\stackrel{\sim }{v}}_{\mathrm{mode}}\propto {\left(\frac{{\rho }_b{v}_b^4}{{\rho }_{b,\odot }{v}_{b,\odot }^4}\right)}^{1/2},\end{array}$$(34)

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 F-type 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 F-type 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:

$$\begin{array}{c}\hfill {\stackrel{\sim }{v}}_{\mathrm{mode}}\propto {\left(\frac{{\rho }_b{v}^4}{{\rho }_{b,\odot }{v}_{b,\odot }^4}\right)}^{1/2}{\left(\frac{\mathit{Ro}}{{\mathrm{Ro}}_{\odot }}\right)}^a,\end{array}$$(35)

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⁻⁵ and 1 × 10⁻⁹ (e.g. Garaud 2021), radiative diffusion dominates the damping of the waves as they travel through the radiative regions, in a quasi-adiabatic regime (Zahn et al. 1997), which should allow modes to form below a much lower frequency cutoff than observed in our simulations. The ${\stackrel{\sim }{v}}_r$ 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 g-mode power spectrum.

This work opens some observational perspectives for F-type solar-type-pulsating 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 disk-integrated 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 follow-up, 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 low-frequency regions of the brightest main-sequence F-type stars with solar-type oscillations observed by Kepler. As our simulations show that ℓ = 4, 5, 6, and 7 modes are also particularly excited, long-term observations of some F-type 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. Christensen-Dalsgaard 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 magneto-hydrodynamic 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 high-order g-modes 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 (Foreman-Mackey et al. 2013).

References

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 $\widehat{v_r}(r,\ell ,m,\omega )$ of the complex quantity v_r(r, ℓ, m, t), with −ω_N < ω < ω_N. From the spherical harmonic properties, we have $\widehat{v_r}(r,\ell ,-m,\omega )=\widehat{v_r}(r,\ell ,m,-\omega )$. 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:

$$\begin{array}{c}\hfill \{\begin{array}{ccc}\hfill E_{\ell ,m}(r,\omega )& =|\widehat{v_r}{(r,\ell ,m,\omega )|}^2;\hfill & \hfill m>0,\\ \hfill E_{\ell ,m}(r,\omega )& =|\widehat{v_r}{(r,\ell ,-m,-\omega )|}^2;\hfill & \hfill m<0,\\ \hfill E_{\ell ,0}(r,\omega )& =|\widehat{v_r}{(r,\ell ,0,\omega )|}^2+{|\widehat{v_r}(r,\ell ,0,-\omega )|}^2.\hfill \end{array}\end{array}$$(A.1)

To compute the power spectrum E_ℓ(r, ℓ, ω) and restrict the domain to the positive frequencies ω ≥ 0, we consider for each degree ℓ

$$\begin{array}{c}\hfill E_{\ell }(r,\omega )={\displaystyle \sum _{m=-\ell }^{\ell }}{E}_{\ell ,m}(r,\omega );\omega \ge 0.\end{array}$$(A.2)

All Tables

All Figures

Fig. 1. Radial profiles of the Brunt-Väisälä frequency N in the ASH model (black) and the MESA model (dashed grey). The structural and chemical contributions to NMESA, NT, MESA, and Nμ, MESA, are represented in dotted-orange and dotted-dark-blue 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 $\overline{\rho }$ (blue) and temperature $\overline{T}$ (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 Prot measured by Santos et al. (2021) in the Kepler data, considering stars with 6000 < Teff < 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 non-uniform radial grid used in the 3D ASH simulations.
In the text

Fig. 6. Flux balance between radiative luminosity Lrad (dashed orange), kinetic energy luminosity Lke (dashed red), enthalpy luminosity Len (dashed grey), diffusivity luminosity Lν (dashed green), and subgrid-scale eddy luminosity Led (dashed blue) for the F5 case. The total luminosity, Ltot, 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, vr, normalised by the rms velocity, ${\stackrel{\sim }{v}}_r$ (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 vr (left) and temperature perturbation $T-\overline{T}$ (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. ${\stackrel{\sim }{v}}_r$ profiles for the F1a (orange), F1b (cyan), F3 (blue), and F5 (red) cases. The vertical dashed-grey 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 Rof 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 co-rotational frame) flows are shown in blue and red, respectively. The dashed black line corresponds to rCZ and the dashed grey lines are isocontours.
In the text

	Fig. 13. Radial velocity vr (left) and temperature perturbation $T-\overline{T}$ (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 rc at which the enthalpy flux crosses zero at a given latitude while the orange line signals the radius r0 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 ray-tracing 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 $v_r/{\stackrel{\sim }{v}}_r$ 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 g-mode 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 E1 (0.38 R⋆, ω). Middle panels: power spectra E1 (0.38 R⋆, P). Bottom panels: Lomb-Scargle periodograms computed from E1 (0.38 R⋆, P). Each periodogram is normalised with its standard deviation, σ and the obtained ΔPl 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. Lomb-Scargle periodograms computed from E2 (0.38 R⋆, P). Each periodogram is normalised with its standard deviation, σ and the obtained ΔP2 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. ΔP1 as a function of modes period in the F5 case, for azimuthal numbers m = −1 (yellow), m = 0 (blue), and m = 1 (red). The ΔP1 obtained from the Ω⋆ = 0 GYRE run and the perturbative approach are represented as grey hexagons while the white crosses shows the ΔP1 obtained for m = { − 1, 0, 1} with the 5 Ω⊙ GYRE TAR run. The black crosses signal the ΔP1 position for m = { − 1, 1} in the asymptotic approximation of the perturbative method. The measured uncertainties on period and ΔP1 are represented. We show only fitted modes for which we are able to measure ΔP1 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, E1, 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, E4, 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 g-modes 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. Lomb-Scargle periodograms computed from E3 (0.97 R⋆, P) (top) and E5 (0.97 R⋆, P) (bottom). The periodogram is normalised with its standard deviation σ and the obtained ΔP3 value is shown by a vertical dashed red line.
In the text

	Fig. 28. Luminosity perturbation estimation as a function of vr 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