Effects of stratification on overshooting and waves atop the convective core of $5M_{\odot}$ main-sequence stars

As a massive star evolves along the main sequence, its core contracts, leaving behind a stable stratification in helium. We simulate 2D convection in the core at three different stages of evolution of a $5M_{\odot}$ star, with three different stratifications in helium atop the core. We study the propagation of internal gravity waves in the stably-stratified envelope, along with the overshooting length of convective plumes above the convective boundary. We find that the stratification in helium in evolved stars hinders radial motions and effectively shields the radiative envelope against plume penetration. This prevents convective overshooting from being an efficient mixing process in the radiative envelope. In addition, internal gravity waves are less excited in evolved models compared to the zero-age-main-sequence model, and are also more damped in the stratified region above the core. As a result, the wave power is several orders of magnitude lower in mid- and terminal-main-sequence models compared to zero-age-main-sequence stars.


INTRODUCTION
One-dimensional formalisms for mixing due to convective overshooting and internal gravity waves have the potential to improve the comparison between stellar evolution models and observations.Several mechanisms have been proposed to contribute to the radial transport of chemical elements: mainly overshooting of plumes from a convective region in a surrounding radiative region (Shaviv & Salpeter 1973;Stancliffe et al. 2016), transport via internal waves (Schatzman 1993), or erosion of compositional gradient via shearing (possibly enhanced by rotation, Zahn 1992).The mixing induced by these mechanisms is however typically too weak to convincingly explain several observations.Such observations include evolutionary tracks of eclipsing binaries (Claret & Torres 2016), constraints from color-magnitude diagrams of stellar clusters in the Large Magellanic Cloud (Rosenfield et al. 2017), and constraints from asteroseismology (Bossini et al. 2015).Recently, Baraffe et al. (2023) showed that convective overshooting above the cores of massive stars determined from 2D hydrodynamical simulations is too weak to explain the observed width of the main sequence in the Hertzsprung-Russell diagram.Convective overshooting was measured in simulations at the Zero-Age Main Sequence (ZAMS).The present study is a follow-up of the latter, where we also measure the extent of convective overshooting at later times along the main sequence: mid main sequence (MidMS), ★ E-mail: a.morison@exeter.ac.uk and terminal main-sequence (TAMS).This study provides hints at how convective overshooting changes as stars evolve along the main sequence and whether this mechanism is a viable candidate for radial mixing.
Another open question in the field of stellar dynamics is that of the propagation of internal gravity waves (IGW) generated at the convective core boundary and whether they can be observed (Bowman et al. 2019(Bowman et al. , 2020;;Lecoanet et al. 2019Lecoanet et al. , 2021)).Recently Vanon et al. (2023) showed simulations of evolved 7 ⊙ stars on the main sequence where the propagation of IGW is greatly affected by the Brunt-Väisälä peak above the convective core.Ratnasingam et al. (2023) showed simulations of 3 − 13 ⊙ at a mid-main-sequence age, also pointing out the effects of the Brunt-Väisälä peak on the propagation of waves.We perform a similar study for our three models of 5 ⊙ stars along the main sequence.Two key differences of our simulations compared to previous ones is that ours are fully compressible instead of anelastic, and use temperature and helium mass fractions profiles from stellar evolution models instead of prescribing the  2 -peak region via a modified temperature profile.
In Section 2, we detail the models corresponding to the three evolutionary stages along the main sequence.We present our setup for the 2D numerical simulations with the fully compressible time-implicit code MUSIC, and how we assess the penetration length above the convective core in those simulations.We then present the theoretical framework we use to study internal gravity waves that propagate in the radiative region.Section 3 groups our results, showing the con-  1. Physical parameters of the three evolutionary stages. ★ is the total radius of the star,  conv is the position of the Schwarzschild boundary of the convective core,  ,conv is the pressure scale height at that boundary,  is the total luminosity of the star, and  core is the helium mass fraction in the core.
sequences of the helium stratification on the dynamics, penetration length, and IGW propagation.We end this paper with discussions and conclusive remarks in Section 4.

METHODS
We focus on three evolutionary stages of a 5 ⊙ star along the main sequence: zero age (ZAMS), mid main sequence (MidMS) and terminal main sequence (TAMS).All three models are part of the same evolutionary path of a 5 ⊙ -star, with an initial helium mass fraction  = 0.28 and metallicity  = 0.02.Table 1 summarizes the three cases.

1D stellar models
We use the one-dimensional (1D) stellar evolution code presented in Baraffe & El Eid (1991); Baraffe et al. (1998) to compute 1D models for each of those three stages.
Figure 1 shows profiles from these three 1D models.One interesting feature is the helium mass fraction profiles that vary widely above the convective core across the three evolution stages considered here.The core retracts and its helium content increases as a result of hydrogen burning as the star evolves along the main sequence.The convective core is assumed to be well-mixed and therefore to have a radially constant helium fraction, and most of the stable envelope stays at the initial helium mass fraction  = 0.28.Due to the core retracting, a layer therefore develops atop of the convective core with a strong stratification in helium.This gradient of concentration in helium directly translates into a strong density gradient above the convective core, and therefore a region where the Brunt-Väisälä frequency is much higher than in the rest of the radiative envelope (around 400 µHz versus 100 µHz).This region is hence referred to as the  2 -peak region in this paper.
The 1D stellar models used in this study rely on the Schwarzschild criterion for the onset of convective instability and do not account for overshooting at the convective core boundary.Because the convective core retracts during core H-burning, the presence of this  2 -peak layer and of a molecular weight gradient above the convective core is robust.Assumptions used in the 1D models can certainly affect the properties of this layer (width, radial profile of the molecular weight gradient), but not its existence.Miglio et al. (2008) find for instance that accounting for convective overshooting in a 6 ⊙ model only marginally increases the radial extent of the  2 -peak layer but not its amplitude; they report however more drastic differences for 1.6 ⊙ stars.Therefore we expect that for 5 ⊙ stars, the impacts of the  2 -peak layer on convective penetration and waves that we report in this work will remain qualitatively the same independently of the assumptions used to construct the 1D models.The dotted box shows the additional radial domain in the "TAMS ext" model (identical to the TAMS model but with a higher external radius).See eq. ( 18) for the definition of  and eq.( 25) for the definition of   .

2D numerical simulations
The 1D profiles described previously are fed as initial condition for the fully compressible stellar hydrodynamics code MUSIC (Viallet et al. 2016;Goffrey et al. 2017).We solve the mass, momentum, internal energy, and helium conservation equations for a fully com-pressible inviscid fluid: where  is the fluid density,  the specific internal energy,  the fluid velocity, and  the helium mass fraction.The gravitational acceleration  = −() r is purely radial and updated after each time step from the instantaneous density profile ρ(): The density profile ρ is calculated as where ⟨•⟩  denotes an angular average across the domain accounting for sphericity We recover the pressure  and temperature  from  and  using tabulated equation of state, which for the simulated regime of 5 ⊙ main-sequence stars are the OPAL tables from Rogers & Nayfonov (2002).Radiative transfer is accounted for via a diffusive formulation involving the heat conductivity .The latter is derived from the Rosseland mean opacity  given by the OPAL opacity tables (Iglesias & Rogers 1996): (with  the Stefan-Boltzmann constant).Note that the 1D stellar evolution code and MUSIC use the same equation of state tables and opacities for consistency.Finally,  nuc is the specific energy released by nuclear burning.Its radial profile is taken from the 1D model (see Figure 1), and considered constant in time throughout the 2D run since nuclear burning evolves on much longer timescales than the duration of the simulations.This is also why no variation of  due to nuclear reactions appears in eq. ( 4).
All the simulations performed in this study are 2D in (, ) space, assuming   = 0 and invariance along  (see Table 2 and hereafter for details).For all three cases, we solve in a domain that spans the entire latitudinal range ( ∈ [0, ]), and truncate along the radial direction ( ∈ [ in ,  out ], see Table 2).All simulations are performed on a uniform mesh at a similar resolution per pressure scale height at the Schwarzschild boundary (denoted  conv ):  ,conv ∼ 140Δ, with Δ the grid size in radial direction.With this resolution, the jump in helium mass fraction above the core spans across ∼ 40Δ in the MidMS model and ∼ 80Δ in the TAMS model.We use reflective boundary conditions at  = 0 and  = .At  in and  out , we impose reflective boundary conditions on velocity, a constant radial derivative in density (following e.g.Baraffe et al. 2023), and the energy fluxes taken from the 1D model at those radii.
MUSIC uses a fully implicit time scheme (Viallet et al. 2016); timesteps are therefore not limited by stability considerations.The timestep is set based on limits on both the acoustic and advective CFL.We pick the maximum timestep Δ that satisfies both Δ   +|  | Δ ⩽ 100 and Δ |  | Δ ⩽ 0.1 for every velocity component   everywhere in the simulation domain, where   is the local sound speed.This corresponds to an acoustic and advective Courant number of 100 and 0.1, respectively.In practice, this leads to timesteps of the order of a minute.
As the initial conditions for the 2D solver are built from 1D models, the 2D models go through an initial transient as the 2D convective motions and waves develop.In this study, we focus on the statistical steady state reached by the models after this transient stage.
To make sure our simulation time in steady state is long enough to capture meaningful statistics, we use the convective timescale  conv defined as where ⟨•⟩  denotes a time average across the steady state. rms is the instantaneous root-mean-square velocity profile: The models being 2D and of relatively modest resolution allows us to run the simulations for a long time, achieving about a thousand of convective times for all three models (see Table 2).Such simulation times are necessary for the diagnostics related to overshooting length based on extreme events to converge (Pratt et al. 2017;Baraffe et al. 2023).A quantity associated with the convective timescale is the convective frequency Notice on Table 2 that the convective frequency is higher in models that are further along on the main sequence (due to higher convective velocities), leading to higher frequency waves being excited."TAMS ext" in Table 2 is a check for the impact of the low cut-off radius on internal gravity waves in evolved models: it is identical to the "TAMS" model but with a larger external radius, to verify that the features observed with the small cut-off radius are also observed in the larger domain.

Penetration length
A first aspect that this study focuses on is the impact of the stratification in helium atop the core on overshooting length.This follows the method developed in Pratt et al. (2017) and used in Baraffe et al. (2021Baraffe et al. ( , 2023)), of which an overview is done here.This method is built around two different fluxes: the radial kinetic energy flux and the radial heat flux The temperature fluctuations  are the departure from the mean temperature profile of the steady state: The overshooting distance   with respect to   is then defined as where   is the first radius above  conv at which   = 0. Similarly,   (, ) is defined from the zero-crossing of   .For such low Table 2. Setup for 2D simulations of the three evolution stages."TAMS ext" is an extended simulation that uses the same initial condition as the TAMS model, but with a higher external radius.This is to check the impact of the fairly low cut-off radius on internal gravity waves in the other models. conv is the convective timescale eq. ( 9),  conv the convective frequency eq. ( 11),  in the inner radius of the simulated domain,  out the outer radius of the simulated domain,  conv the Schwarzschild radius of the convective core,  steady the duration of the steady-state of the simulation,   ,   the number of grid cells along the radial and latitudinal directions,  ,conv the pressure scale height at  conv ,  ★ the radius of the star, and Δ the grid cell size along the radial direction.
mass models, we do not observe any significant penetration (significant modification of the temperature gradient above the convective core was found in Baraffe et al. (2023) for more luminous models); the position of the convective boundary  conv in eq. ( 15) is therefore taken from the 1D model (see Table 1).Following Pratt et al. (2017), we assess the overshooting length via the angular maximum of   (, ): and, following Baraffe et al. (2023), its time average defines the overshooting length and similarly for  max  and  ov  .Several hundreds of overturn times (see Table 2) are needed for  ov  and  ov  to converge to a steady value.Once that steady state reached, both diagnostics give a similar value for the overshooting length, as discussed in Pratt et al. (2017); Baraffe et al. (2023).

Waves study
The second aspect on which this study focuses is the impact of the stably stratified layer atop of the core on the excitation and propagation of internal gravity waves in the envelope.This follows the work presented in Le Saux et al. (2022Saux et al. ( , 2023)).The reader is referred to those studies for more details, an overview of the method is presented here for convenience.
In a stably-stratified fluid, a particle of fluid that is out-ofequilibrium (i.e.denser or lighter than its surroundings) is brought back to its equilibrium position by the buoyancy force.This mechanism allows IGW to propagate in a stably-stratified fluid.The angular Brunt-Väisälä frequency  (in rad/s) is the maximum frequency of such waves.It is directly related to the density stratification of the medium: Γ 1 is the first adiabatic exponent: A fluid is unstable against convection if  2 ⩽ 0, in which case IGW are evanescent.Conversely, a fluid is stably-stratified if  2 > 0, and IGW of angular frequency  ⩽  can propagate in the medium.Note that in this paper,  denotes the angular frequency in rad/s, and  = /(2) is the frequency in hertz.In the absence of rotation, the dispersion relation of IGW verifies (Press 1981): where   is the radial wavenumber and  ℎ the horizontal wavenumber.The latter is related to the spherical harmonic degree ℓ of the considered wave: The linear theory developed in Press (1981) and extended to variable molecular weight  in Zahn et al. (1997) uses the WKB approximation to predict the amplitude of travelling IGW in the radiative region.This yields the following expression for the radial velocity: where -indices denote quantities evaluated at the excitation radius.Note that v (ℓ, ) is an arbitrary amplitude for the linear theory.In eq. ( 22),  represents the damping of waves due to radiative effects.Radiative diffusion is the main damping mechanism of IGW in stellar interiors, and the only one considered in this linear model. verifies (Zahn et al. 1997): rad is the radiative diffusivity where   is the specific heat capacity as constant pressure.  is the thermal contribution to the angular Brunt-Väisälä frequency A related quantity is the wave luminosity, defined as (Press 1981;Zahn et al. 1997): where  , is the radial group velocity This luminosity is only affected by damping as waves travel in the radiative region: To identify and study IGW in the radiative envelope of our simulations, we use a temporal Fourier transform and decomposition in spherical harmonics of the radial velocity field of the 2D models.[ v ] (, ℓ,  ) denotes the power spectrum obtained from such a spectral analysis (see Appendix A of Le Saux et al. 2022, for more details).This power spectrum is directly the rms of the radial velocity.Using  ℎ /  ≃   / ℎ and the dispersion relation eq. ( 20), the wave luminosity eq. ( 26) can then be expressed from this power spectrum: The linear expression on v eq. ( 22) also applies to For simplicity, we test the linear prediction on profiles of √︁ [ v ] directly.The excitation radius   in eq. ( 22) is taken slightly above the convective core boundary, at twice the overshooting length eq. ( 17): Note that taking   at  ov  or 2 ov  above the convective boundary does not change our results significantly.To improve the readability of figures, we choose v (ℓ, ) in eq. ( 22) as the local maximum of √︁ [ v ] evaluated around   .Comparing the wave amplitudes in our 2D models with this linear theory constitutes a diagnostic of whether our 2D solver is capable to properly solve for IGW propagation, at least for waves that are indeed in the linear regime considered here.

General observations on the dynamics
Figure 2 shows the instantaneous temperature perturbation for the ZAMS and MidMS cases.The TAMS case exhibit similar behavior to the MidMS case.The convective part is extremely well mixed with |  | ≲ 10 −6 while the stably stratified region shows stronger deviations from the mean profile (|  | ∼ 10 −4 ) that organize in well-defined patterns.Those patterns are the signature of internal gravity waves.These plots already exhibit a striking difference between the ZAMS model and evolved models: patterns of temperature perturbation have a similar angular extent and amplitude throughout the stable layer in the ZAMS case, while two shells with a different organization of   are clearly visible across the stable region of the MidMS model.Sitting atop the core is the  2 -peak region, showing patterns of   with larger angular and radial extent than in the rest of the stable region at lower  2 .Section 3.3 presents an in-depth analysis of the wave content in those simulations.
As can be seen on Figure 3, the kinetic energy density is constant in the convective core, decreases sharply at the convective-core boundary, and slightly decreases with radius in the stably stratified region for all three models.The kinetic energy density in the convective region is higher in the evolved models compared to the ZAMS model by about a factor of three.This is expected as the luminosity of the star (and therefore the energy feeding convection) increases with age along the main sequence.However, the drop of kinetic energy at the convective boundary is much deeper for evolved models (4 orders of magnitude) than for the ZAMS model (about 3 orders of magnitude), resulting in an about 10-fold lower kinetic energy density in the stable region of evolved models.The  2 -peak layer of evolved models reduces the transfer of kinetic energy from the convective core to the radiative envelope.Figure 3 provides an additional insight on the effects of the  2 -peak layer.While the radial component of the kinetic energy density is roughly monotonic in the ZAMS model and most of the MidMS and TAMS radiative envelopes, it dips by more than an order of magnitude in the  2 -peak layer of evolved models compared to the rest of the stable region.This is readily explained by the dispersion relation eq. ( 20) which implies that  2 / 2  ≃  2 / 2 ; as  2 increases, the horizontal component of the kinetic energy density will therefore carry more of the energy.

Penetration length
This subsection focuses specifically on the transfer of matter across the convective boundary via the estimation of an overshooting length as described in section 2.3.
Average values of the overshooting lengths  ov  and  ov  are shown in Table 3 for the three models.For the ZAMS case, one recovers the value found in Baraffe et al. (2023)    ∼ 0.06  ).For both evolved cases MidMS and TAMS, the penetration length associated with either flux is much lower  ∼ 0.025  .It should be noted that  (, ) is only one grid cell Δ ∼ 0.007  in many areas, meaning the penetration length is likely under-resolved and possibly over-estimated for those evolved cases.
Overshooting is therefore severely limited in the MidMS and TAMS models (likely due to the stabilizing stratification in helium that builds up), and is unlikely to be an efficient mixing and wave excitation mechanism above the convective cores of massive stars which have evolved along the main sequence.

Wave study
The previous subsection shows that the  2 -peak layer hinders radial motions in the evolved models.This subsection focuses on the consequences on IGW in the radiative region.Figure 4 shows the wave luminosity spectrum of all three models for ℓ = 3. Vertical bright lines are the signature of g-modes, i.e. standing IGW that form in the radiative zone.Dark lines follow the radial nodes of the standing modes.As expected for IGW, the number of nodes increases with decreasing  (Aerts et al. 2010).A notable feature of both evolved models compared to the ZAMS model is that nodes are narrowed in the  2 -peak layer in evolved models.This is, again, explained by the dispersion relation eq. ( 20): for a given frequency , increasing the Brunt-Väisälä frequency leads to a more horizontal wavefront.We used the stellar oscillation code GYRE (Townsend & Teitler 2013) to compute the frequency of g-modes in the simulated cavity.Figure 4 shows a good match between the predicted frequencies (white dashed lines) and the position of the peaks in the power spectra, especially for the MidMS and TAMS model.The match is less convincing for the ZAMS model at higher frequency, in particular for the −2 and −3 g-modes.This could be due to the stronger reflection that results in noisier signal at high frequencies for this model (see Figure 6).Figure 5 compares power spectra at (ℓ = 5,  = 7.6 µHz) and (ℓ = 10,  = 1.2 µHz) with the linear theory eq. ( 22) for all three models, with damping following eq.( 23) and without damping ( = 0).Note that since we use realistic values for the heat conductivity  (see eq. ( 8)), we expect our 2D models to exhibit realistic damping profiles.The (ℓ = 5,  = 7.6 µHz) mode undergoes very limited damping in the simulated region, and as can be seen on Figure 5, the simulated spectrum matches well the linear theory for all three models.On the other hand, the (ℓ = 10,  = 1.2 µHz) mode is predicted to be severely damped (which is a consequence of the higher ℓ and lower  compared to the previous mode).The match between the linear prediction and the simulated spectrum is much less convincing for this mode.The ZAMS model does not match the prediction.The MidMS model departs from the linear prediction in many places, with slopes that are different from the linear prediction everywhere except in a small region above the  2 -peak.The TAMS model is the one that matches the most the linear prediction for this mode, with a good match in the  2 -peak region and roughly following the predicted slope in the rest of the stratified region.The increase in amplitude when exiting the  2 -peak layer is systematically higher in the simulations than what is predicted by the linear theory.This From top to bottom: ZAMS, MidMS, and TAMS.The dashed lines locate the g-modes predicted by GYRE (above 10 µHz), with their radial order within the Eckart-Scuflaire-Osaki-Takata scheme (Takata 2006).Note that the labels for modes -8 to -10 (MidMS) and modes -8 to -16 (TAMS) are omitted for readability purposes.The dotted line on the TAMS spectrum marks the position of an alias of the -2 mode (predicted at 56.7 µHz) due to spectral folding around the Nyquist frequency 50 µHz.The plot on the right is the Brunt-Väisälä frequency profile of each model, note that the frequency scale is different from the spectra.linear theory is developed using the WKB approximation, and therefore assumes that spatial derivatives of various parameters (such as  2 ) are small; this assumption is potentially not valid at the top of the  2 -peak region, and the linear theory may thus incorrectly predict the change of amplitude as waves exit that region.
The two modes discussed above were chosen as they are in different regimes: one is only very weakly damped while the other undergoes strong damping.Similar observations can be made on other modes: high-frequency modes are subject to low damping and match well the theoretical prediction, while modes at low frequency are more strongly damped and the simulated spectrum does not always match the linear prediction.The poorer match of low-frequency modes could be due to their higher amplitude than high-frequency modes (see Figure 4): they are further away from the linear regime in which the linear theory is valid.Moreover, low-frequency IGW have a shorter radial wavelength and therefore are not as well resolved radially, which could partly explain the poorer match for those modes.Indeed, the radial wavelength   from the dispersion relation eq. ( 20) is well above the cell width Δ for the mode (ℓ = 5,  = 7.6 µHz) across all models (6Δ <   < 60Δ), but the radial wavelength for the mode (ℓ = 10,  = 1.2 µHz) is comparable to 2Δ for all three models and even smaller in the  2 -peak region: the latter mode is therefore under-resolved.It should finally be noted that the linear theory eq. ( 22) concerns travelling waves, while waves in the simulated domain are subject to reflection due to the reflective boundary condition at the outer boundary, and to the strong gradient of  2 at the top of the  2 -peak region.Reflection of waves in the overall domain and/or within the  2 -peak region could explain the higher amplitude of some modes in the 2D models compared to the linear prediction.Moreover, large oscillations characteristic of g-modes make it difficult to compare to the linear theory.
Figure 6 compares the wave power spectra at ℓ = 3 for the three stellar models at two radii: at the excitation radius eq. ( 30), and 0.07 ★ above it (outside of the  2 -peak region in evolved models).The damping at low  ≲ 2 µHz is clearly visible when comparing both radial positions.The damping is stronger in the evolved models (three orders of magnitude) compared to the ZAMS model (two orders of magnitude).The gaussian shape aroud the peak of power spectra (around 2 µHz to 5 µHz), usually suggested to be related to the excitation of IGW by convective penetration (e.g.Pinçon et al. 2016;Edelmann et al. 2019;Le Saux et al. 2022), is less clearly marked in the evolved models than in the ZAMS model.This could be a signature of the much lower penetration length in evolved models (see section 3.2).The power spectra align along the slope predicted for IGW excited by Reynolds stress at higher frequencies, although the spread in amplitude due to high-amplitude g-modes in the simulated spectra makes it impossible to discriminate which model is the most relevant in the present case.
Figure 6 exhibits extremely marked peaks in wave luminosity above 10 µHz for the ZAMS model, that have much lower amplitude in the more extended but otherwise identical simulation shown in Le Saux et al. (2023).Those strong peaks are probably the result of wave reflections in the truncated domain.The spectra for the evolved case do not seem to be affected as much by reflection, which could be a result of the overall lower kinetic energy in evolved models compared to ZAMS (Figure 3).
As a check, we performed an extended simulation of the TAMS model (with  out = 0.59 ★ ); Figure 7 shows the obtained wave luminosity spectrum and comparison with the linear theory for the (ℓ = 15,  = 2.5 µHz) mode.The wave luminosity spectrum shows similar features to the ones previously discussed for the truncated model.Note that the g-modes are at different frequencies owing to the different aspect ratio of the cavity in the extended case.The simulated spectrum follows roughly the slopes predicted by the linear theory for  ≲ 0.2 ★ (with, again, a stronger increase in amplitude than predicted when exiting the  2 -peak layer at  ≃ 0.13 ★ ), but strong oscillations due to g-modes renders the comparison to the theory difficult at higher radius.
Figure 8 shows the theoretical wave power spectrum for the three evolution stages, using the 1D initial profiles which extend up to  =  ★ .This shows damping is very efficient at frequencies lower than about 20 µHz, likely preventing most IGW from reaching the surface.IGW that are less damped and susceptible of reaching the surface are at higher frequencies, but those frequencies are weakly excited.Their predicted amplitude at the surface is then extremely small.This confirms the results of Le Saux et al. (2023); Anders et al. (2023) for ZAMS models and extends their conclusion to models evolved on the main sequence: IGW excited by core convection likely have a very low amplitude at the surface.

CONCLUSIONS
The three simulations presented in this paper show several consequences of the stratification in helium that builds up on top of the convective core of massive stars during the evolution on the main sequence.The stratification in helium is mechanically stable and strongly limits the penetration of plumes in the radiative envelope.This hinders mixing by overshooting in evolved stars compared to ZAMS stars.Moreover, the convective frequency increases as stars evolve along the main sequence, IGW excited by convection are then at higher frequency in older stars.IGW propagation is limited in evolved stars due to several effects: damping in the  2 -peak region, reflections at the top of the  2 -peak region due to the steep change in stratification, and/or weaker excitation by plume penetration.This reduces the overall power of waves susceptible to reach the surface of an evolved star.This strengthen the conclusions from Le Saux et al. (2023); Anders et al. (2023).Our results clearly highlight that the amplitude of IGW propagating in the radiative envelope of intermediate mass stars is expected to be small.It is thus even less likely to detect IGW excited by core convection in these stars, casting doubts on the interpretation of Bowman et al. (2019Bowman et al. ( , 2020)).
Finally, Baraffe et al. (2023) shows mixing by overshooting as calculated for ZAMS models is already too weak to explain the observed width of the main sequence in the Hertzsprung-Russell diagram for stellar masses  ≳ 10 ⊙ .Our study further strengthens this conclusion by showing that overshooting only gets weaker compared to ZAMS as stars evolve along the main sequence and a stable stratification in helium builds up above the convective core.
Our results complement recent studies by Vanon et al. (2023); Varghese et al. (2023); Ratnasingam et al. (2023).They present simulations in similar setups, studying the interaction between the convective core and radiative zone of a massive star along the main sequence.It should be underlined that those studies present 2D/3D simulations under the anelastic approximation, without an helium mass fraction field, and with modified heat fluxes, diffusivities and temperature profiles compared to stellar evolution models to mimic the  2 -peak region above the core of evolved models.Our simulations on the other hand are 2D, fully compressible, solving for helium mass fraction conservation, and with diffusivities and temperature profiles consistent with 1D stellar evolution models.In particular, including a realistic helium mass fraction profile rather than mimicking the  2 -peak via a modified temperature profile leads to a more realistic damping profile in the  2 -peak region as we do not conflate   and  in eq. ( 23).
There are some differences between our findings and those previous studies.Ratnasingam et al. (2023); Vanon et al. (2023) find their power spectra have a better agreement with plume excitation, while we find a better agreement with Reynolds stress excitation (Figure 6).Vanon et al. (2023) report discrepancies between their TAMS simulations and theoretical predictions of g-modes with GYRE.They suggest some compressible effects that are not resolved by their anelastic solver could be at play here, for example mixed modes.We however find a good agreement between GYRE predictions and our TAMS models.Finally, Varghese et al. (2023);Vanon et al. (2023) find that IGW become evanescent in radial domain for TAMS models, and suggest that the WKB approximation does not apply in those models.This is in contrast with our results, where we find a good agreement between the WKB linear theory and our simulations (particularly for weakly-damped modes where this approximation should be relevant), and do not observe the waves to become evanescent.The reasons for all those discrepancies are unclear, and might stem from the unrealistic temperature and diffusivity profiles used in these stud-  Wave power spectrum and linear prediction for the three models, for two different modes.On each plot, the light solid line is the spectrum obtained from the simulation, the dark solid line is the linear prediction with radiative damping eq. ( 22), and the dotted line is the linear prediction without radiative damping.These two lines are superimposed for the ( = 5,  = 7.6 µHz) mode as damping is small.The black dashed line marks the convective boundary  conv .Wave luminosity spectrum profiles at ℓ = 3 for the three cases, at two different radii above the convective boundary.  is slightly above the convective boundary (at twice the penetration length  ov  found in section 3.2), and   + 0.07 ★ is outside of the  2 -peak region of the MidMS and TAMS models.The thin light-colored solid lines are the spectra obtained from the simulations, and the thick dark-colored solid lines are the same profiles smoothed for readability purposes.The vertical gray dash-dotted line is the convective frequency.Dotted lines are gaussian models of excitation by plume penetration (Pinçon et al. 2016), and dashed lines are scalings considering Reynolds stress (Lecoanet & Quataert 2013).4, the dashed lines denote the g-modes predicted by GYRE.Labels for nodes -8 to -15 are omitted for readability.Bottom: as in Figure 5, the light solid line is the spectrum obtained from the simulation, the dark solid line is the linear prediction with radiative damping eq. ( 22), and the dotted line is the linear prediction without radiative damping.
ies.A comparison with the recent work by Thompson et al. (2024) based on 3D simulations of a 25 ⊙ is also difficult since it does not include radiative diffusion.As we show in present work and in Le Saux et al. (2023), the propagation of low frequency waves in the radiative zone is strongly impacted by radiative effects.
Our main findings are based on 2D numerical simulations and do not consider the effect of rotation on convection and wave properties.Vanon et al. (2023) find similar power spectra in their 3D simulations as the 2D ones from Rogers et al. (2013), suggesting the power spectra we show in this 2D study should be largely unaffected by dimensionality.Convective velocities in 2D simulations are usually larger than in 3D simulations (Meakin & Arnett 2007;Pratt et al. 2020), suggesting larger overshooting lengths predicted by 2D sim- ulations compared to 3D models.However, the plume filling factor, which is linked to the convective flow structure, also plays a role in the determination of an effective overshooting length and is expected to be different between 2D and 3D simulations, with smaller filling factors producing smaller overshooting lengths (Zahn 1991;Brummell et al. 2002;Rogers et al. 2006;Pratt et al. 2020).How the plume structure, which affects the determination of an overshooting length, varies between 2D and 3D descriptions is still an open question.One can however note that our results are qualitatively consistent with Vanon et al. (2023) who also report a decrease of overshooting as stars evolve along the main sequence.We are currently performing a detailed study of the filling factor and plume shape at the convective boundary of stellar cores and envelopes using the same simulation framework in 2D and 3D to address this question (Dethero et al., 2024, in prep.).Regarding rotation, it can impact convective structures as they align with the rotation axis, affecting both convective velocities and wave excitation.As the Coriolis force deflects radial convective plumes at a convective-radiative interface, one expects decreasing overshooting lengths as the rotation rate increases (Brummell et al. 2002;Rogers et al. 2013).The effects of rotation and of the presence of a helium stratification above the core for stars on the Main Sequence would thus combine to limit convective penetration, making the problem found in ZAMS stellar stars to reproduce observations (Baraffe et al. 2023) even worse for rotating MS stars.
Regarding the impact of rotation on waves, recent 3D simulations by Anders et al. (2023) of a 15 ⊙ ZAMS model suggest that moderate rotation rate can slightly boost the luminosity of waves excited at the convective core boundary.This boosting effect on the wave luminosity may be significant for fast rotators, but it will still be limited by the presence of a strongly stratified layer above the convective core as the star evolves on the main sequence.While the power excess observed at the surface of intermediate mass stars and interpreted as due to internal waves excited by the convective core (Bowman et al. 2019) is ubiquitous, all these stars are not fast rotators and not on the ZAMS.
To conclude, our work suggests that overshooting becomes weaker along the main sequence as also found by Vanon et al. (2023).This aggravates the problem pointed out in Baraffe et al. (2023) to reproduce observations in the Hertzsprung-Russell diagram.In addition, the stabilizing helium stratification that builds up above the core of massive stars on the main sequence strongly affects the excitation and propagation of IGW.This limits the amplitude of IGW reaching the surface in evolved models compared to ZAMS models.Understanding IGW and overshooting signatures in the context of long-term stellar evolution therefore requires further efforts on the modelling and observational fronts.

Figure 1 .
Figure 1.Radial profiles for the three stages of evolution from the 1D evolution model.Those profiles are then interpolated to use them as initial condition for the 2D simulations.Dashed lines mark the position of the convective core boundary for each model according to the Schwarzschild criterion.Shaded areas mark the radial extent of the simulated domain for each evolution stage.The dotted box shows the additional radial domain in the "TAMS ext" model (identical to the TAMS model but with a higher external radius).See eq.(18) for the definition of  and eq.(25) for the definition of   .

Figure 2 .Figure 3 .
Figure 2. Temperature perturbation maps for the ZAMS (left) and MidMS (right) cases.The solid black line marks the position of the convective boundary.

Figure 4 .
Figure 4. Wave luminosity spectrum for all three models, along with the Brunt-Väisälä frequency profiles and g-modes.From top to bottom: ZAMS, MidMS, and TAMS.The dashed lines locate the g-modes predicted by GYRE (above 10 µHz), with their radial order within the Eckart-Scuflaire-Osaki-Takata scheme(Takata 2006).Note that the labels for modes -8 to -10 (MidMS) and modes -8 to -16 (TAMS) are omitted for readability purposes.The dotted line on the TAMS spectrum marks the position of an alias of the -2 mode (predicted at 56.7 µHz) due to spectral folding around the Nyquist frequency 50 µHz.The plot on the right is the Brunt-Väisälä frequency profile of each model, note that the frequency scale is different from the spectra.

Figure 5 .
Figure5.Wave power spectrum and linear prediction for the three models, for two different modes.On each plot, the light solid line is the spectrum obtained from the simulation, the dark solid line is the linear prediction with radiative damping eq.(22), and the dotted line is the linear prediction without radiative damping.These two lines are superimposed for the ( = 5,  = 7.6 µHz) mode as damping is small.The black dashed line marks the convective boundary  conv .

Figure 7 .
Figure7.Wave luminosity spectrum for the extended TAMS model (top), and comparison with the linear theory for the (ℓ = 15,  = 2.5 µHz) mode (bottom).Top: as in Figure4, the dashed lines denote the g-modes predicted by GYRE.Labels for nodes -8 to -15 are omitted for readability.Bottom: as in Figure5, the light solid line is the spectrum obtained from the simulation, the dark solid line is the linear prediction with radiative damping eq.(22), and the dotted line is the linear prediction without radiative damping.

Figure 8 .
Figure 8. Wave luminosity spectrum calculated from the linear theory eq.(28) for all three models at ℓ = 3.The dashed line marks the convective frequency  conv .

Table 3 .
for the 5L0 case (  ∼ 0.09  , Overshooting length measured for the three evolution stages.