© The Authors 2026

1. Introduction

With their peculiar structural properties, F-type stars populate the region between the Sun and intermediate-mass stars (see e.g. Fig. 1 from Breton et al. 2023). Their radiative interior is surrounded both by a shallow convective envelope and a small convective core, which makes them fascinating test-beds to explore the non-linear excitation of internal waves by decoupled convective sources. When they are able to form standing modes, these internal waves give rise to stellar oscillations (e.g. Unno et al. 1989; Aerts et al. 2010). In particular, direct constraints on the core properties and on the stellar deep dynamics can only be obtained through the observation of gravity modes (g modes, e.g. Christensen-Dalsgaard et al. 1985; García et al. 2007), which are the standing modes of oscillation related to the propagation of internal gravity waves (IGWs) in the stellar interior. However, due to their low surface amplitude when a convective envelope is present, the detectability of such g modes is unclear (Belkacem et al. 2009, 2022; Breton et al. 2023) below the γ Dor instability strip (e.g. Kaye et al. 1999; Warner et al. 2003; Dupret et al. 2004).

F-type main-sequence stars with an effective temperature below 6500 K exhibit stochastically excited acoustic modes of oscillation (p modes, e.g. Belkacem et al. 2008; Chaplin et al. 2011; Mathur et al. 2014; Lund et al. 2017) and are therefore still considered solar-like from this point of view. They are of paramount importance for our understanding of stellar objects below 1.5 M_⊙, since a significant part of the upcoming Planetary Transit and Oscillation of Stars mission (PLATO, Rauer et al. 2025) core sample will be made up of such main-sequence F-type solar-like pulsators (Montalto et al. 2021; Goupil et al. 2024). Although p-mode asteroseismology enables us to put some constraints on the extent of the convective core of F-type stars, this requires the implementation of model-dependent fitting of observational frequencies on a stellar evolutionary track (e.g. Deheuvels et al. 2016), or the application of complex inversion techniques (e.g. Bétrisey & Buldgen 2022). Moreover, p modes alone are unable to provide constraints on the rotational dynamics of the deep radiative layers (e.g. Thompson et al. 1996, 2003; García et al. 2007; Benomar et al. 2015). On the contrary, the core extent is directly encoded in the g-mode asymptotic period spacing (Tassoul 1980), and g modes are extremely sensitive to the rotation of such regions (Deheuvels et al. 2012; Li et al. 2020). The main issue for the application of g-mode asteroseismology on main-sequence solar-type stars stems from the fact that, unlike γ Dor stars that exhibit g modes with large amplitudes due to the existence of non-adiabatic destabilising mechanisms (e.g. Dupret et al. 2005), g modes in solar-type stars are stable and can only be excited by the stochastic energy input from turbulent convective motions (e.g. Kumar et al. 1996; Lecoanet & Quataert 2013; Pinçon et al. 2016; Augustson et al. 2020). Previous 2D and 3D hydrodynamical simulations in polar and spherical configurations demonstrated that IGWs are more efficiently excited in F-type stars (Breton et al. 2022a) than in the Sun (e.g. Rogers & Glatzmaier 2005; Brun et al. 2011; Alvan et al. 2014, 2015).

In main-sequence solar-type stars, g modes are low-frequency (ν ≲ 500 μHz) modes with an absolute radial order, n, that increases with decreasing frequency. In F-type stars, low-order high-frequency modes are in principle able to mix with low-order p modes (Breton et al. 2023) to form mixed modes, while high-order low-frequency g modes will become gravito-inertial modes under the effect of the Coriolis force (Lee & Saio 1997; Townsend 2003; Van Reeth et al. 2015, 2018), with possible additional couplings between the eigenmodes of the core, the radiative interior, and the envelope (Ouazzani et al. 2020; Tokuno & Takata 2022; Barrault et al. 2025; Breton et al. 2026). Unambiguously observing and characterising g modes in late F-type stars, under any of their manifestations, would bring new, strong structural constraints on stellar evolution processes in solar-type stars. Moreover, measuring their rotational splittings would open a direct window onto the rotation of the deep radiative layers. Both of these aspects would represent an invaluable asset for our understanding of stellar structure and evolution.

In this work, we present the first 3D simulations of a stellar interior including both a convective core and a convective envelope. We show that the g modes encode the structural extent of the core in their properties. The modes maintain their coherence up to the highest regions of the convective envelope. With respect to models that only account for envelope excitation, we find that the core plays a crucial role by filtering out high-frequency, high-degree modes. By constructing observational time series from our simulations, we are able to unambiguously identify the signature of the modes at the top of the convective envelope. This should justify and encourage future efforts for their identification and characterisation in observational time series. The structure of the article is as follows. In Sect. 2, we present the anelastic numerical setup of our simulations. In Sect. 3, we compute the power spectrum of the g modes excited in the simulation in order to analyse the different behaviour exhibited by the setups depending on the mode excitation source. Section 4 is dedicated to the computation of pseudo-observational time series and the discussion of the mode detectability. Finally, we draw the conclusions of our work in Sect. 5.

2. Simulation setups

In order to assess the importance of the convective core in the internal wave dynamics of an F-type stars, we developed and we analysed three 1.3 M_⊙ spherical models with the Anelastic Spherical Harmonics code (ASH, Clune et al. 1999; Brun et al. 2004). The first one includes the full extent of the stellar structure (convective envelope surrounding the radiative interior and the convective core; the model is named ‘Full’ hereafter). The second one only includes the convective envelope and the radiative interior (model named ‘EnvCZ’ hereafter), and the third one only includes the radiative interior and the convective core (model named ‘CoreCZ’ hereafter). We note that, despite their similar configuration, the background structure of the EnvCZ model slightly differs from the model presented in Breton et al. (2022a). Each simulation is initialised with small perturbations in the unstable convective layer to allow for the growth of convective instability, and evolved until it statistically reaches a relaxed steady state. More details on the numerical setup used for the simulations are provided in Appendix A.

A volume representation of the radial velocity field, v_r, for the Full model is shown in Fig. 1. To enable a visualisation of both the convective structure and the IGW pattern, v_r is normalised by the root-mean-square value, v_rms, computed on the radial shell to which each cell belongs. In the convective envelope, large-scale upwards and downwards flows are visible. In particular, downwards motions can be seen interacting with the top of the radiative interior through plumes penetrating the stably stratified layers and exciting the IGWs. Given the more limited extent of the core, the typical scale of the convective motions close to the centre of the star is smaller than for the envelope. As shown in Fig. A.3 and expected from a mixing-length argument (Böhm-Vitense 1958), the convective motions are also slower in the core than in the envelope. In the upper layers of the radiative interior, the signal is constituted by the superposition of progressive IGWs and standing g modes of all degrees. Close to the convective core, the wave signal is dominated by high-amplitude, high-frequency g modes. For comparison with the Full model, the volume display of the EnvCZ and CoreCZ models are shown in Fig. 2. We note that, due to the absence of the convective core, the superposition of mode patterns in the lower half of the radiative interior differs between the Full and the EnvCZ cases. Finally, the wave structure we distinguish in the CoreCZ case has a large scale even close to the top of the radiative interior, which suggests that the signal is dominated by low-order modes with few radial nodes.

Fig. 1. Left: Volume display of the radial velocity field in the Full model, normalised by the shell root-mean square value, vr/vrms. The inset on the left displays a zoom on the convective core. Top right: Volume display of vr/vrms after filtering the signal to isolate the signal at 121 μHz ≤ ν ≤ 123 μHz. Bottom right: Same as above, but isolating the signal at ν ≤ 5 μHz.

Fig. 2. Volume display of the radial velocity field in the EnvCZ model (left) and the CoreCZ model (right), normalised by the shell root-mean square value, vr/vrms.

In Fig. 3 we show the profile of the Brunt-Väisäla frequency, N, of the reference stellar structure. We have N = 0 in the convective layer and N > 0 in the radiative layers. Given that N imposes the maximal frequency of IGWs propagating in the local medium, waves in the radiative interior are expected to have frequencies up to 400 μHz. As the ray trajectory of an IGW is frequency-dependent (Gough 1993), it is possible to filter the temporal signal to disentangle the signature of waves at distinct frequencies (Alvan et al. 2015). We illustrate this in the right panel of Fig. 1 where we show the result of this filtering in the radiative interior for a frequency of 121 μHz ≤ ν ≤ 123 μHz and ν ≤ 5 μHz. In the former case, the ray travels from inner caustic to outer caustic, with a reflection geometry that creates the typical rosetta pattern of IGWs in this range of frequency. In the latter case, at very low frequency, the rays converge towards the centre in a spiral pattern, with winding angles so low that they appear as concentric circles.

Fig. 3. Profile of the Brunt-Väisäla frequency of the reference stellar structure, N. The extent of the convective core, the radiative interior, and the convective envelope are indicated in colour. The value of N was set to zero in the convective regions where N2 < 0. The frequency limit where intermediate- and high-degree modes are unable to form in the Full model while visible in the EnvCZ model is indicated with the dotted grey line.

3. Mode power spectrum

Once the simulations are relaxed, we computed the power spectrum of the radial velocity at given depths for our different setups. We considered time series of lengths 104 days.

3.1. Power spectrum computation

For each 3D model, the power spectrum was computed by constructing time series with the simulation outputs (Breton et al. 2022a). Considering the radial grid of the simulations, a sub-sample of spherical shells with radial coordinates, r, were selected. The corresponding (θ, ϕ) velocity maps, with θ the latitude and ϕ the longitude, were output with a regular sampling dt = 600 s, which corresponds to a Nyquist frequency of approximately 833 μHz, which is larger than the maximal value of our structural N profile. We computed the spherical harmonic transform of each spherical map in order to obtain the ℓ, m coefficients of the velocity field, with ℓ the angular degree and m the azimuthal number. Given that these high-cadence time series represent a very large amount of data, we computed this transform up to ℓ = 100, as a trade-off between preservation of the signal at a low-angular scale and data storage availability. We then applied a discrete Fourier transform (DFT). The transformation process can then be summarised as

$$\begin{array}{c}\hfill v_r(r,\theta ,\varphi ,t)\to {\stackrel{\sim }{v}}_r(r,\ell ,m,t)\to {\widehat{v}}_r(r,\ell ,m,\nu ),\end{array}$$(1)

where the complex quantities ${\stackrel{\sim }{v}}_r$ and ${\widehat{v}}_r$ are computed from the radial velocity v_r. Given the 104-day length of the time series, the frequency bin width is ∼0.11 μHz. To obtain the energy coefficients with respect of each angular degree, E_ℓ, that we represent in Fig. 4, we computed the quadratic sum of the modulus of the m-components for each ℓ:

$$\begin{array}{c}\hfill E_{\ell }(r)={\displaystyle \sum _{m=-\ell }^{m=\ell }}{|{\widehat{v}}_r(r,\ell ,m,\nu )|}^2.\end{array}$$(2)

Fig. 4. Power spectra Eℓ at depths r = [0.243, 0.485, 0.728] R★ for the Full model (first row), the EnvCZ model (second row), and the CoreCZ model (third row). The degree-dependent cut-off frequency between the standing modes and the progressive waves is shown on every panel with a dashed black line. In the first two rows, the limit frequency where intermediate- and high-degree modes are unable to form in the Full model while visible in the EnvCZ model is indicated with the dotted black line. The border colour of the panels is different for each row in order to guide the eye (each colour corresponds to a model: orange, blue, and yellow for the Full, EnvCZ, and CoreCZ models, respectively).

3.2. Impact of the excitation regions on the power spectrum

The power spectrum of the signal at three depths, r = [0.243, 0.485, 0.728] R_★, expanded on the spherical harmonics with degrees up to ℓ = 100, is shown in Fig. 4. For each ℓ, the power contribution of each azimuthal number, m, has been summed quadratically. The spectra from the three models are extremely different, which emphasises that neither the presence of the envelope nor the presence of the core can be omitted in order to model the behaviour of IGWs in the radiative interior.

Internal gravity waves propagate inside the resonant cavity located inside the radiative interior. The location of the top and bottom boundaries of these cavities only depends on the wave frequency and corresponds to the local value of the Brunt-Väisälä frequency, N (see Fig. 3). The IGWs in each spectrum are separated into two distinct regions (Alvan et al. 2014; Ahuir et al. 2021; Breton et al. 2022a), with a boundary materialised by the degree-dependent cut-off frequency represented as the dashed grey lines in Fig. 4. Below this cut-off the progressive IGWs are unable to form standing modes because, excited at one end of the resonant cavity, they are damped by radiative diffusion before reaching the opposite boundary (Zahn et al. 1997). The spectrum of progressive waves is similar between the Full model and the EnvCZ model, while no progressive IGW is visible in the CoreCZ model. This underlines the fact that the convective envelope alone is able to excite these progressives waves for a large range of frequencies and degrees.

Above the cut-off frequency, standing modes are visible; launched from one boundary, the IGWs travel through the entirety of the resonant cavity and reach the opposite boundary with sufficient energy. There, a fraction of the wave energy is reflected towards the cavity, allowing the standing mode to form, while another fraction is transmitted outside the cavity as an evanescent tail. This reasoning corresponds to an asymptotic representation of the wave behaviour (n ≫ 1), where the boundary of the resonant cavity can be validly identified as the radius where ν_nℓ = N, with ν_nℓ the mode frequency. In the following, we see that non-asymptotic modes have a more complex behaviour. Indeed, it is very striking that the rich pattern of high-degree, high-frequency modes exhibited by the EnvCZ model is not present for the Full model. We also note that, at a high frequency and intermediate degrees (5 ≲ ℓ ≲ 20), the CoreCZ model also exhibits a slightly richer mode spectrum. We further illustrate these findings with Figs. 5 and 6, where we show the evolution of E_ℓ as a function of depth for ℓ = 5 and ℓ = 15, respectively. For ℓ = 5 modes, low-order modes are visible in the three models. We note that the signature of the mode extends beyond the N = ν boundary. On the contrary, for ℓ = 15 modes, the mode behaviour above 150 μHz is very different from one model to another. As expected from Fig. 4, the EnvCZ model has a rich spectrum of excited eigenmodes above 150 μHz, while only a few are visible in the CoreCZ case. Eigenmodes above this frequency are barely visible in the Full case, and they are orders of magnitude below the EnvCZ model in terms of power. Before discussing the condition that allows the modes to form, we emphasise that, despite this depletion in high-degree modes, the Full model exhibits large-amplitude low-degree modes, and these are the ones most susceptible to be detected in disc-resolved observations (García & Ballot 2019).

Fig. 5. Power spectra Eℓ as a function of depth for ℓ = 5. The thick grey lines correspond to the Brunt-Väisälä frequency, while the dotted vertical black line indicates the 150 μHz frequency limit discussed in the text. The hashed areas in grey correspond to the stellar regions that are not included in the EnvCZ and CoreCZ models.

Fig. 6. Same as Fig. 5 for ℓ = 15.

3.3. Interpreting the convective core’s role in mode formation

As highlighted by the dotted line in Fig. 4, the modes that are unable to form in the Full case, while being excited to large amplitudes in the EnvCZ case, have a frequency ν ≳ 150 μHz. Below this frequency, the lower boundary of the mode resonant cavity approximately coincides with the interface between the core and the radiative interior (see Fig. 3). We suggest that the distinct behaviour observed for these modes in the Full and EnvCZ cases can be interpreted in terms of mode asymptoticity. We argue that the waves that correspond to an asymptotic mode will perceive the ν_nℓ = N boundary (in this case, the top of the convective core) as a transmission-reflection interface. It can be demonstrated through the use of the Werner, Jeffreys, Brillouin, and Kramers (WKBJ, e.g. Unno et al. 1989) method that their radial eigenfunction is oscillatory within the bounds of the cavity and is thus defined and evanescent outside. Under this asymptotic configuration, in the EnvCZ case, the inner boundary condition imposes complete reflection at r = 0.07 R_★ (it is the bottom of the simulation domain), while in the Full case a sufficient fraction of the wave energy is reflected and a standing mode is still able to form.

By contrast, the WKBJ approximation is no longer valid for non-asymptotic modes. A travelling wave will not perceive the ν_nℓ = N boundary as a transmission-reflection interface and it will propagate beyond this point. In the EnvCZ case, the wave will therefore reach the bottom of the simulation domain and will be reflected there. The situation is different when the dynamical picture includes the convective core. The wave would need to be refracted close to the centre of the star to form a coherent standing mode. However, as it enters the overshoot region, the wave starts experiencing viscous damping from the convective eddies, which adds to the radiative damping. In this picture, the wave is damped while travelling in the convective core, and is therefore unable to form a mode. In the asymptotic picture, the damping efficiency increases with the degree, ℓ, and the distribution we obtain in the power spectrum for the Full case supports that this should still be valid as only low-degree non-asymptotic modes are able to form. We illustrate our argument with Fig. 7, where we compare the eigenfunctions computed with the GYRE code (Townsend & Teitler 2013, see Appendix B) for the simulation reference states, for two ℓ = 15 modes. The first mode, with n = 5 and ν = 290.8 μHz, is excited in the EnvCZ case but not in the Full case. The second mode, with n = 29 and ν = 100.1 μHz, is excited in both cases. In particular, as visible from the location of the vertical dashed lines, the oscillatory behaviour of the eigenfunction in the asymptotic case is confined inside the whole extent of the radiative interior, while there is no clear distinction between an oscillatory component and an exponential one in the non-asymptotic case.

Fig. 7. Radial displacement (ξr) computed with the Full reference state (orange) and the EnvCZ reference state (dotted blue). The top panel shows a non-asymptotic ℓ = 15, n = 5 mode with ν = 290.8 μHz, while the left panel shows a ℓ = 15, n = 5 mode with ν = 100.1 μHz, which can be interpreted as asymptotic. The vertical dashed grey lines correspond to the bounds of the radiative interior while the vertical dashed yellow lines highlight the radial location of the νnℓ = N boundary. On both panels, ξr was normalised by its maximal value.

4. Mode detectability

We now turn to the issue of mode detectability. Observationally, the main issue stems from the fact that the low-surface amplitude of g modes makes it difficult to disentangle their contribution from the incoherent turbulent noise that dominates the Fourier spectrum in this range of frequency (e.g. García et al. 2007; Breton et al. 2023). In order to examine this issue with our simulation, we developed a process to compute pseudo-observational time series.

4.1. Observational time series reconstruction

In order to produce the pseudo-observational power spectral densities (PSD), we applied the following process to our data. We denote θ₀ and ϕ₀ as the co-latitude and the longitude corresponding to the observer line of sight. The resultant projected velocity, v_proj, in this direction is therefore the sum of the three components v_r, proj, v_θ, proj, and v_ϕ, proj

$$\begin{array}{c}\hfill v_{\mathrm{proj}}(\theta ,\varphi )=v_{r,\mathrm{proj}}(\theta ,\varphi )+v_{\theta ,\mathrm{proj}}(\theta ,\varphi )+v_{\varphi ,\mathrm{proj}}(\theta ,\varphi ),\end{array}$$(3)

where

$$\begin{array}{cc}\hfill v_{r,\mathrm{proj}}(\theta ,\varphi )& =v_{r}cos(\theta -{\theta }_0)cos(\varphi -{\varphi }_0),\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill v_{\theta ,\mathrm{proj}}(\theta ,\varphi )& =v_{\theta }sin({\theta }_0-\theta )cos(\varphi -{\varphi }_0),\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill v_{\varphi ,\mathrm{proj}}(\theta ,\varphi )& =v_{\varphi }sin({\varphi }_0-\varphi )sin{\theta }_0,\hfill \end{array}$$(6)

where v_θ and v_ϕ are the latitudinal and longitudinal component of the fluid velocity, respectively. The angle ψ between the normal vector of a given surface element of the sphere and the observer line of sight follows the relation, denoting μ = cos ψ (Lanza 2016; Breton et al. 2024a)

$$\begin{array}{c}\hfill \mu (\theta ,\varphi )=sin{\theta }_{0}sin\theta cos(\varphi -{\varphi }_0)+cos{\theta }_{0}cos\theta ,\end{array}$$(7)

and only the elements with μ > 0 are visible to the observer. In what follows, we consider for simplicity ϕ₀ = 0 and θ₀ = π/2 in the co-rotating frame. This means that our fictitious observer is rotating around the star with the same angular velocity. Nevertheless, the advantage of this approach is that the same surface elements continually remain in the observer’s field of view, and we do not introduce a discrete spurious modulation related to the appearance and disappearance of surface elements when they cross the μ = 0 threshold. We model the limb-darkening contribution, ℒ_d, with the Eddington law (García & Ballot 2019) so that

$$\begin{array}{c}\hfill {\mathcal{L}}_d(\mu )=(0.4+0.6\mu )\mu .\end{array}$$(8)

Finally, the observed velocity, v_obs, is written

$$\begin{array}{c}\hfill v_{\mathrm{obs}}=\frac{{\displaystyle {\displaystyle \sum _{i,j}}{a}_i{v}_{\mathrm{proj}}({\theta }_i,{\varphi }_j){\mathcal{L}}_d({\mu }_{i,j}){\mu }_{i,j}}}{{\displaystyle {\displaystyle \sum _{i,j}}{a}_i{\mu }_{i,j}}},\end{array}$$(9)

where the indices (i, j) identify each surface element and a_i = sin θ_i is the corresponding relative surface. In what follows, we reconstruct v_obs considering angular degrees ℓ from 0 to 100, which is up to the largest angular degree we stored for the time series. Nevertheless, the additional tests we performed showed that the reconstruction did not significantly differ considering a truncation at a significantly lower ℓ, such as ℓ = 10.

In our simulations, the convective noise contribution at a high frequency mostly results from a power leakage of the low-frequency signal associated with convection. This means that, by applying a high-pass filter on our reconstructed time series, we are able to increase the signal-to-noise ratio of the periodic signatures at a higher frequency, which corresponds, in our case, to the g-mode signal.

In order to remove the contamination from the convective noise in the high-frequency range of the PSD, a high-pass finite-impulse-response (FIR) filter with a 12-hour (∼23 μHz) cut-off is therefore applied on the time series (Breton et al. 2022b, 2024b) before computing the Fast Fourier Transform (FFT). In Fig. 8, we show the time series we obtain before applying the FIR filter. As shown, due to the contribution of the convective motions, the amplitude of the v_obs signal significantly varies as a function of the depth, with a maximum amplitude of about 1 × 10³ cm s⁻¹ at r/R_★ = 0.97 and only 10 cm s⁻¹ at r/R_★ = 0.80, in the radiative interior. It also clearly appears that the long-period variations of the time series at r/R_★ = 0.97 are strongly correlated with the average variation in the convective envelope, which reflects the impact of the large-scale structures that are visible in Fig. 1. In Fig. 9, we show a six-day chunk of the same time series after applying the 12-hour cut-off FIR filter. With the long-period trend removed, the typical amplitude of the signal is now much smaller. As can be seen by comparing the three rows of the figure, some contribution from the convection in the period range shorter that 12 hours remains, and enhances the total amplitude of the signal observed in the convective envelope with respect to the radiative interior.

Fig. 8. Time series obtained for vobs before filtering, at r/R★ = 0.97 (top row), averaged in the convective envelope (0.88 ≤ r/R★ ≤ 0.97, middle row), and r/R★ = 0.80 (bottom row), for the Full (left column) and the EnvCZ (right column) models, respectively.

Fig. 9. Six-day chunk of the time series obtained for vobs after filtering, at r/R★ = 0.97 (top row), averaged in the convective envelope (0.88 ≤ r/R★ ≤ 0.97, middle row), and r/R★ = 0.80 (bottom row), for the Full (left column) and the EnvCZ (right column) models, respectively.

4.2. Observational PSD

By integrating the projected velocity signal on the stellar disc, as would be seen by an observer, we were therefore able to construct a set of time series that allows us to model how the disc-integrated observable g-mode signature evolves with the stellar depth. In Fig. 10, we show the corresponding PSD obtained after applying a DFT on these time series. It is expected that, in this observational PSD, the contribution from increasing degree ℓ will be gradually filtered out by the disc-integrated spatial integration (Dziembowski 1977). Moreover, by temporally filtering the low-frequency contribution to the time series, we are able to unambiguously isolate the g-mode signal in the PSD, even in the upper layers of the convective envelope. Here again, the signature differs importantly between the Full case and the EnvCZ model. In both cases, the main g-mode signature contribution comes from the ℓ = 1 and ℓ = 2 modes. As visible, the EnvCZ case exhibits a few peaks related to the contribution from modes at higher degrees (this is also the case for the mode with the lowest frequency in the Full spectrum). The existence of this signature demonstrates without ambiguity that the evanescent tail of the g modes maintains its coherence up to the top of the convective envelope and should be able to reach the stellar photosphere. Removing the incoherent contribution from convective noise allows us to reveal the presence of the g-mode pattern, regularly spaced in periods for each degree ℓ (Tassoul 1980). The convective core has an imprint not only on this period spacing, but also by limiting the number of g modes. In particular, our results highlight the need to look for the combined signature of ℓ = 1 and ℓ = 2 modes in observational data, as they should reach the surface of the star with comparable amplitude.

Fig. 10. PSD of the disc-integrated projected velocity signal, vobs, obtained from the simulation in the upper radiative interior (r/R★ = 0.80, left column), averaged in the convective envelope (0.88 ≤ r/R★ ≤ 0.97, middle column), and close to the top of the simulation domain (r/R★ = 0.97, right column). The PSD obtained for the Full model is shown at the top and the EnvCZ model PSD is at the bottom. The regular pattern of the ℓ = 1 modes is highlighted with vertical dashed grey lines and that of the ℓ = 2 modes with the vertical dotted blue lines.

5. Conclusion

In this work, we presented the first 3D simulations of the dynamics of a star with two convective zones. F-type stars possess both an external convective envelope and a small convective core that surrounds a stably stratified radiative interior. In order to investigate the respective contribution from the core and the envelope in the excitation of g modes, we compared a full-extent model (Full) with a coreless model (EnvCZ) and a model with no convective envelope model (CoreCZ). We also highlighted the crucial role that the convective core plays regarding the excitation of the modes, as it prevents high-degree modes from forming in the non-asymptotic high-frequency range. By integrating the signal on the stellar disc in order to compute pseudo-observational time series, we demonstrated that, despite their evanescence, g modes maintain their coherent integrity up to the top of the convective turbulent layer, which supports the claim that they have a detectable observable signature at the stellar photosphere.

This work paves the way for a systematic search for g-modes signatures in the F-type solar-like pulsators that will be observed by the PLATO mission. A first attempt at such an analysis was performed by Breton et al. (2023) considering the 34 F-type solar pulsators with the best characterisation in the Kepler mission (Borucki et al. 2010; Mathur et al. 2014; Davies et al. 2016; Lund et al. 2017; Silva Aguirre et al. 2017), but definitive conclusions could not be drawn given the limited size of the sample. Given the low amplitude of the modes, they will have to be looked for in a large sample of stars and using analysis techniques built on solid statistical ground.

Acknowledgments

The authors thank the anonymous referee for comments and suggestions that helped improving the manuscript. S.N.B. acknowledges support from PLATO ASI-INAF agreement no. 2022-28-HH.0 “PLATO Fase D”. SNB acknowledges support from the INAF grant MASTODINT. A.S.B. acknowledges support from Solar Orbiter CNES grant and financial support by ERC Whole Sun Synergy grant #810218. While partially funded by the European Union, views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. A.S.B. acknowledges CNES support for PLATO WP123400. R.A.G. acknowledges support from PLATO and GOLF CNES grants.

References

Appendix A: Hydrodynamical setup

The ASH numerical setups presented in Sect. 2 are defined in the spherical co-rotating frame (r, θ, ϕ), with unit vectors (e_r, e_θ, e_ϕ). The reference density, pressure, temperature, and specific entropy are denoted as $\overline{\rho }$, $\overline{P}$, $\overline{T}$, and $\overline{S}$, with corresponding fluctuations ρ, P, T, S. They are connected through the equation of state and the zeroth-order ideal gas law

$$\begin{array}{cc}\hfill \frac{\rho }{\overline{\rho }}& =\frac{P}{\overline{P}}-\frac{T}{\overline{T}}=\frac{P}{\gamma \overline{P}}-\frac{S}{c_p},\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill \overline{P}& =\mathcal{R}\overline{\rho }\overline{T},\hfill \end{array}$$(A.2)

where c_p is the specific heat per unit mass at constant pressure, γ the adiabatic exponent, and R is the gas constant. The quantities S and c_p can be combined with the gravitational acceleration g = −ge_r to define the Brunt-Väisälä frequency, N, through the relation

$$\begin{array}{c}\hfill N^2=\frac{g}{c_p}\frac{\mathrm{d}S}{\mathrm{d}r}·\end{array}$$(A.3)

In the WKBJ approximation, N > 0 defines the resonant cavities of IGWs. It is also directly related to the asymptotic period spacing of consecutive g modes of same degree, ΔP_ℓ, through (Tassoul 1980)

$$\begin{array}{c}\hfill {(\mathrm{\Delta }{P}_{\ell })}^{-1}=\frac{\sqrt{\ell (\ell +1)}}{2{\pi }^2}{\displaystyle {\int }_{r_1}^{r_2}}\frac{N}{r}\mathrm{d}r,\end{array}$$(A.4)

where r₁ and r₂ correspond to the boundaries N² = 0. N is represented in Fig. 3. We note that the maximal value of N is reached in the deep radiative interior, close to r/R_★ = 0.25.

The hydrodynamic equations are solved in the Lantz-Braginsky-Roberts (LBR, Lantz 1992; Braginsky & Roberts 1995) formulation of the anelastic approximation (Gough 1969). By filtering out acoustic waves, the anelastic approximation allows for larger integration time steps with respect to fully compressible setups. The LBR formulation is specifically chosen because it was shown that it preserves better the IGWs energy (Brown et al. 2012). The momentum equation is

$$\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}$$(A.5)

where $\stackrel{\sim }{\omega }=P/\overline{\rho }$ is the reduced pressure fluctuation, v = (v_r, v_θ, v_ϕ) is the local velocity, Ω₀ = Ω₀e_z is the rotation frequency of the reference frame. The viscous stress tensor, [BOLD]𝒟, is

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

with ν_diff the effective viscous diffusivity, 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}$$(A.7)

Our equation of conservation of internal energy writes

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

where κ_r is the radiative diffusivity (Kippenhahn et al. 2012). The ϵ term accounts for the energy generation by nuclear burning, modelled as a temperature power-law, $ϵ={ϵ}_0{\overline{T}}^k$. Similarly to ν_diff, κ_diff is an effective diffusivity while the diffusivity κ₀ is set to carry the unresolved entropy eddy flux in the convective zone near the surface (see Brun et al. 2004). We set κ_diff and ν_diff as a trade-off between the desired turbulent regime and numerical tractability. In particular, we decrease both κ_diff and ν_diff in the radiative interior in order to minimise the IGWs damping. The ν_diff and κ_diff are described by an equation of the form

$$\begin{array}{cc}& {\nu }_{\mathrm{diff}}={\nu }_1[\beta +\frac{1-\beta }{2}{\left(\frac{{\overline{\rho }}_{\mathrm{top}}}{\overline{\rho }}\right)}^{\frac{1}{2}}(tanh\frac{r-r_{t1}}{{\sigma }_t}+1)]\hfill \\ \hfill & +\frac{{\nu }_2}{2}(tanh\frac{r-r_{t2}}{{\sigma }_t}+1).\hfill \end{array}$$(A.9)

with, in both case, the transition radii r_t1 = 8.9 × 10¹⁰ cm and r_t2 = 7.3 × 10⁹ cm, the stiffness σ_t = 8 × 10⁸ cm. This way, at the top of the simulation domain ν_diff, top ≈ ν₁, in the radiative interior, ν_{diff, radiative} ≈ βν₁, and in the convective core, ν_diff, core ≈ ν₂ + βν₁. We use ν₁ = 1.05 × 10¹³ cm² s⁻¹, ν₂ = 3.395 × 10¹⁰ cm² s⁻¹, κ₁ = 4.2 × 10¹³ cm² s⁻¹, and κ₂ = 2.716 × 10¹⁰ cm² s⁻¹. For ν_diff, β = 1 × 10⁻⁴, and for κ_diff, β = 2 × 10⁻⁵. Both profiles are represented in Fig. A.1.

Fig. A.1. νdiff (solid line) and κdiff (dashed line) diffusivity profiles used in the simulations. The extent of the convective core, the radiative interior, and the convective envelope are indicated in colour as in Fig. 3.

We finally have to define the boundary conditions. For the Full and the CoreCZ models, the singularity at r = 0 is treated with a regularisation procedure that was previously validated for the study of IGWs (Alvan et al. 2014). For every other boundary, we adopt the following conditions: rigid, stress-free impenetrable boundaries with constant mean-entropy gradient (Breton et al. 2022a).

From the problem formulation, the total energy balance can be expressed as the sum of the flux related to the different transport processes (Brun et al. 2004)

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

where L_rad is the radiative luminosity, L_ke the kinetic energy luminosity, L_ν the diffusive processes energy luminosity, L_en the enthalpy luminosity, and L_ed the unresolved eddy luminosity. As explained in the main body of the text, small perturbations introduced in the convective unstable layer at the start of the simulation allow for the convective instability to grow and build the L_ke luminosity. The simulation is evolved for several times the convective turnover time, τ_CZ, in order to allow it to reach a dynamical steady state. In the convective envelope, τ_CZ ≈ 4 days while in the core, τ_CZ ≈ 160 days. We therefore evolve the Full and CoreCZ simulations for about 1250 days in order to reach the dynamical steady state. In Fig. A.2, we show the flux balance achieved by the Full model once it has reached this steady state and is relaxed. As the thermal relaxation timescale (Zahn 1991) of the background state (about 1 × 10⁵ yr) is beyond reach, we adjust L_rad in the overshoot region to maintain the flux balance (Miesch et al. 2000; Breton et al. 2022a). We also show in Fig. A.3, the rms velocity regime that develop in the relaxed simulations for the three cases. Finally, the grid resolution of each model is indicated in Table A.1, with the number of radial (N_r), latitudinal (N_θ), and azimuthal (N_ϕ) mesh points.

Fig. A.2. Flux balance obtained for the Full setup. Lrad (dashed orange), Lke (dashed red), Led (dashed blue), Lν (dashed green), Len (dashed grey), and Ltot are shown. See Eq. (A.10) and subsequent explanations.

Fig. A.3. Root-mean-square (rms) radial velocity profile, vrms, when the steady state is reached for each of the three models. The Full profile is shown in orange, the EnvCZ one in dotted blue, and the CoreCZ one in dashed yellow. The dashed vertical grey lines indicate the boundaries of the convective core and of the convective envelope.

Appendix B: Power spectrum comparison with linear eigenfunctions

The linear eigenfunction shown in Fig. 7 were computed with the GYRE code (Townsend & Teitler 2013), considering the 1D reference profile of our 3D setup. The implemented GYRE formalism closest to the anelastic approximation used in ASH is the γ formalism (Ong & Basu 2020), as it isolates the g-mode buoyancy cavity from the acoustic cavity, similarly to the anelastic approximation where acoustic waves are filtered out. The oscillation equations are solved both with and without the Cowling approximation, which neglects the perturbation of the gravitational potential (Cowling 1941).

In Fig. B.1, we compare the GYRE frequencies for modes with ν_nℓ > 100 μHz and ℓ ≤ 20 to the E_ℓ distribution computed with the ASH simulations. As it is visible in Fig. B.1, the agreement between the GYRE eigenfunctions and the 3D-simulation power spectrum is very good. We note that the first ℓ = 1 mode of the ASH spectra, visible at ν_nℓ ≈ 190 μHz, is present in the GYRE predictions including the Cowling approximation, but absent otherwise. It corresponds to a n = 0 dipolar mode that is unphysical for a self-gravitating spherical object (that would mean the object is oscillating around its centre of mass) but permitted in the context of the Cowling approximation and in ASH setups, where the gravitational potential is defined by a fixed background term with neglected fluctuations.

Fig. B.1. Power spectrum for the Full (top), the EnvCZ (middle), and the CoreCZ (bottom) cases at r = 0.117 R★. The eigenfrequencies computed with GYRE are shown as orange dots for the computations without the Cowling approximation, and as grey dots for the computation with the Cowling approximation. At high degree and high frequency, grey and orange dots are superposed.

All Tables

All Figures

	Fig. 1. Left: Volume display of the radial velocity field in the Full model, normalised by the shell root-mean square value, vr/vrms. The inset on the left displays a zoom on the convective core. Top right: Volume display of vr/vrms after filtering the signal to isolate the signal at 121 μHz ≤ ν ≤ 123 μHz. Bottom right: Same as above, but isolating the signal at ν ≤ 5 μHz.
In the text

	Fig. 2. Volume display of the radial velocity field in the EnvCZ model (left) and the CoreCZ model (right), normalised by the shell root-mean square value, vr/vrms.
In the text

Fig. 3. Profile of the Brunt-Väisäla frequency of the reference stellar structure, N. The extent of the convective core, the radiative interior, and the convective envelope are indicated in colour. The value of N was set to zero in the convective regions where N2 < 0. The frequency limit where intermediate- and high-degree modes are unable to form in the Full model while visible in the EnvCZ model is indicated with the dotted grey line.

In the text

Fig. 4. Power spectra Eℓ at depths r = [0.243, 0.485, 0.728] R★ for the Full model (first row), the EnvCZ model (second row), and the CoreCZ model (third row). The degree-dependent cut-off frequency between the standing modes and the progressive waves is shown on every panel with a dashed black line. In the first two rows, the limit frequency where intermediate- and high-degree modes are unable to form in the Full model while visible in the EnvCZ model is indicated with the dotted black line. The border colour of the panels is different for each row in order to guide the eye (each colour corresponds to a model: orange, blue, and yellow for the Full, EnvCZ, and CoreCZ models, respectively).

In the text

	Fig. 5. Power spectra Eℓ as a function of depth for ℓ = 5. The thick grey lines correspond to the Brunt-Väisälä frequency, while the dotted vertical black line indicates the 150 μHz frequency limit discussed in the text. The hashed areas in grey correspond to the stellar regions that are not included in the EnvCZ and CoreCZ models.
In the text

	Fig. 6. Same as Fig. 5 for ℓ = 15.
In the text

Fig. 7. Radial displacement (ξr) computed with the Full reference state (orange) and the EnvCZ reference state (dotted blue). The top panel shows a non-asymptotic ℓ = 15, n = 5 mode with ν = 290.8 μHz, while the left panel shows a ℓ = 15, n = 5 mode with ν = 100.1 μHz, which can be interpreted as asymptotic. The vertical dashed grey lines correspond to the bounds of the radiative interior while the vertical dashed yellow lines highlight the radial location of the νnℓ = N boundary. On both panels, ξr was normalised by its maximal value.

In the text

	Fig. 8. Time series obtained for vobs before filtering, at r/R★ = 0.97 (top row), averaged in the convective envelope (0.88 ≤ r/R★ ≤ 0.97, middle row), and r/R★ = 0.80 (bottom row), for the Full (left column) and the EnvCZ (right column) models, respectively.
In the text

	Fig. 9. Six-day chunk of the time series obtained for vobs after filtering, at r/R★ = 0.97 (top row), averaged in the convective envelope (0.88 ≤ r/R★ ≤ 0.97, middle row), and r/R★ = 0.80 (bottom row), for the Full (left column) and the EnvCZ (right column) models, respectively.
In the text

Fig. 10. PSD of the disc-integrated projected velocity signal, vobs, obtained from the simulation in the upper radiative interior (r/R★ = 0.80, left column), averaged in the convective envelope (0.88 ≤ r/R★ ≤ 0.97, middle column), and close to the top of the simulation domain (r/R★ = 0.97, right column). The PSD obtained for the Full model is shown at the top and the EnvCZ model PSD is at the bottom. The regular pattern of the ℓ = 1 modes is highlighted with vertical dashed grey lines and that of the ℓ = 2 modes with the vertical dotted blue lines.

In the text

	Fig. A.1. νdiff (solid line) and κdiff (dashed line) diffusivity profiles used in the simulations. The extent of the convective core, the radiative interior, and the convective envelope are indicated in colour as in Fig. 3.
In the text

	Fig. A.2. Flux balance obtained for the Full setup. Lrad (dashed orange), Lke (dashed red), Led (dashed blue), Lν (dashed green), Len (dashed grey), and Ltot are shown. See Eq. (A.10) and subsequent explanations.
In the text

	Fig. A.3. Root-mean-square (rms) radial velocity profile, vrms, when the steady state is reached for each of the three models. The Full profile is shown in orange, the EnvCZ one in dotted blue, and the CoreCZ one in dashed yellow. The dashed vertical grey lines indicate the boundaries of the convective core and of the convective envelope.
In the text

	Fig. B.1. Power spectrum for the Full (top), the EnvCZ (middle), and the CoreCZ (bottom) cases at r = 0.117 R★. The eigenfrequencies computed with GYRE are shown as orange dots for the computations without the Cowling approximation, and as grey dots for the computation with the Cowling approximation. At high degree and high frequency, grey and orange dots are superposed.
In the text