# Effects of stratification on overshooting and waves atop the convective core of M main-sequence stars

A. Morison, A. Le Saux, I. Baraffe, J. Morton, T. Guillet, D. G. Vlaykov, T. Goffrey, J. Pratt

► To cite this version:

A. Morison, A. Le Saux, I. Baraffe, J. Morton, T. Guillet, et al. Effects of stratification on overshooting and waves atop the convective core of M main-sequence stars. Monthly Notices of the Royal Astronomical Society, 2024, 532, pp.2778-2787. 10.1093/mnras/stae1678. insu-04729841

HAL Id: insu-04729841

https://insu.hal.science/insu-04729841v1

Submitted on 10 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Effects of stratification on overshooting and waves atop the convective core of 5 M$_{\odot}$ main-sequence stars

A. Morison<sup>1,★</sup> A. Le Saux<sup>1,†</sup> I. Baraffe,<sup>1,3</sup> J. Morton<sup>1,‡</sup> T. Guillet<sup>1,§</sup> D. G. Vlaykov<sup>1,¶</sup> T. Goffrey<sup>5,6</sup>
and J. Pratt<sup>7</sup>

Accepted 2024 July 6. Received 2024 July 5; in original form 2024 March 21

## ABSTRACT

As a massive star evolves along the main sequence, its core contracts, leaving behind a stable stratification in helium. We simulate two-dimensional convection in the core at three different stages of evolution of a 5 M$_{\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.

**Key words:** asteroseismology – hydrodynamics – waves – software: simulations – stars: interiors.

## 1 INTRODUCTION

One-dimensional (1D) formalisms for mixing due to convective overshooting and internal gravity waves (IGWs) 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 colour–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). This

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) 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 IGWs generated at the convective core boundary and whether they can be observed (Bowman et al. 2019, 2020; Lecoanet et al. 2019, 2021). Recently, Vanon et al. (2023) showed simulations of evolved 7 M$_{\odot}$ stars on the main sequence where the propagation of IGWs is greatly affected by the Brunt–Väisälä peak above the convective core. Ratnasingam et al. (2023) showed simulations of 3–13 M$_{\odot}$ 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 M$_{\odot}$ stars along the main sequence. Two key differences of our simulations compared to previous ones are that ours are fully compressible instead of anelastic, and use temperature and helium mass fraction profiles from stellar evolution models instead of prescribing the N$^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 set-up for the 2D numerical simulations with the fully compressible time-

<footnote n="1"><sup>1</sup>Physics and Astronomy, University of Exeter, Exeter EX4 4QL, UK</footnote>

<footnote n="2"><sup>2</sup>Laboratoire de Météorologie Dynamique (IPSL), Sorbonne Université, CNRS, Ecole Polytechnique, Ecole Normale Supérieure, 4 place Jussieu, 75252 Paris cedex 05, France</footnote>

<footnote n="3"><sup>3</sup>Ecole Normale Supérieure de Lyon, CRAL (UMR CNRS 5574), Université de Lyon, F-69007 Lyon, France</footnote>

<footnote n="4"><sup>4</sup>Mathematics and Statistics, University of Exeter, Exeter EX4 4QF, UK</footnote>

<footnote n="5"><sup>5</sup>Centre for Fusion, Space and Astrophysics, Department of Physics, University of Warwick, Coventry CV4 7AL, UK</footnote>

<footnote n="6"><sup>6</sup>Healthcare Technology Institute, School of Chemical Engineering, The University of Birmingham, Birmingham B15 2TT, UK</footnote>

<footnote n="7"><sup>7</sup>Lawrence Livermore National Laboratory, 7000 East Avenue, Livermore, CA 94550, USA</footnote>

<footnote>* E-mail: a.morison@exeter.ac.uk</footnote>
<table><caption>**Table 1.** Physical parameters of the three evolutionary stages. R\* is the total radius of the star, R\*conv is the position of the Schwarzschild boundary of the convective core, H\*p,conv is the pressure scale height at that boundary, L is the total luminosity of the star, and Y\*core is the helium mass fraction in the core.</caption><thead><tr><th>Model</th><th>R<sub>*</sub> (cm)</th><th>R<sub>conv</sub>/R<sub>*</sub></th><th>H<sub>p,conv</sub>/R<sub>*</sub></th><th>L (L<sub>☉</sub>)</th><th>Y<sub>core</sub></th></tr></thead><tbody><tr><td>ZAMS</td><td>1.84 × 10<sup>11</sup></td><td>0.18</td><td>0.099</td><td>522</td><td>0.28</td></tr><tr><td>MidMS</td><td>2.51 × 10<sup>11</sup></td><td>0.12</td><td>0.070</td><td>680</td><td>0.56</td></tr><tr><td>TAMS</td><td>2.92 × 10<sup>11</sup></td><td>0.095</td><td>0.056</td><td>758</td><td>0.70</td></tr></tbody></table>

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 IGWs that propagate in the radiative region. Section 3 groups our results, showing the consequences of the helium stratification on the dynamics, penetration length, and IGW propagation. We end this paper with discussions and conclusive remarks in Section 4.

## 2 METHODS

We focus on three evolutionary stages of a 5 M☉ star along the main sequence: ZAMS, MidMS, and TAMS. All three models are part of the same evolutionary path of a 5 M☉ star, with an initial helium mass fraction $Y = 0.28$ and metallicity $Z = 0.02$. Table 1 summarizes the three cases.

### 2.1 1D stellar models

We use the 1D stellar evolution code presented in Baraffe & El Eid (1991) and Baraffe et al. (1998) to compute 1D models for each of those three stages.

Fig. 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 $Y = 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 $N^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 $N^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 M☉ model only marginally increases the radial extent of the $N^2$-peak layer but not its amplitude; they report however more drastic differences for 1.6 M☉ stars. Therefore, we expect that for 5 M☉ stars, the impacts of the $N^2$-peak layer on convective penetration and waves that we report

<figure><img src="image_2.png" /><figcaption>**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 equation (18) for the definition of N and equation (25) for the definition of $N_t$.</figcaption></figure>

in this work will remain qualitatively the same independently of the assumptions used to construct the 1D models.

### 2.2 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,
<table><caption>**Table 2.** Set-up 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 IGWs in the other models. $\tau_{\text{conv}}$ is the convective timescale equation (9), $f_{\text{conv}}$ the convective frequency equation (11), $R_{\text{in}}$ the inner radius of the simulated domain, $R_{\text{out}}$ the outer radius of the simulated domain, $R_{\text{conv}}$ the Schwarzschild radius of the convective core, $\tau_{\text{steady}}$ the duration of the steady-state of the simulation, $n_r$, $n_\theta$ the number of grid cells along the radial and latitudinal directions, $H_{p,\text{conv}}$ the pressure scale height at $R_{\text{conv}}$, $R_*$ the radius of the star, and $\Delta r$ the grid cell size along the radial direction.</caption><thead><tr><th>Model</th><th>τ<sub>conv</sub> (s)</th><th>f<sub>conv</sub> (μHz)</th><th>R<sub>in</sub>/R<sub>*</sub></th><th>R<sub>conv</sub>/R<sub>*</sub></th><th>R<sub>out</sub>/R<sub>*</sub></th><th>τ<sub>steady</sub>/τ<sub>conv</sub></th><th>n<sub>r</sub> × n<sub>θ</sub></th><th>H<sub>p,conv</sub>/Δr</th></tr></thead><tbody><tr><td>ZAMS</td><td>1.6 × 10<sup>6</sup></td><td>0.625</td><td>0.02</td><td>0.18</td><td>0.3</td><td>1188</td><td>400 × 200</td><td>141</td></tr><tr><td>MidMS</td><td>8.7 × 10<sup>5</sup></td><td>1.15</td><td>0.02</td><td>0.12</td><td>0.21</td><td>1368</td><td>384 × 192</td><td>142</td></tr><tr><td>TAMS</td><td>6.7 × 10<sup>5</sup></td><td>1.49</td><td>0.02</td><td>0.095</td><td>0.21</td><td>910</td><td>472 × 236</td><td>140</td></tr><tr><td>TAMS ext</td><td>6.4 × 10<sup>5</sup></td><td>1.56</td><td>0.02</td><td>0.095</td><td>0.59</td><td>31</td><td>1416 × 236</td><td>140</td></tr></tbody></table>

internal energy, and helium conservation equations for a fully compressible inviscid fluid:

$$ \frac{\partial \rho}{\partial t} = -\nabla \cdot (\rho \mathbf{v}), \qquad (1) $$

$$ \frac{\partial \rho \mathbf{v}}{\partial t} = -\nabla \cdot (\rho \mathbf{v} \otimes \mathbf{v}) - \nabla p + \rho g, \qquad (2) $$

$$ \frac{\partial \rho e}{\partial t} = -\nabla \cdot (\rho e \mathbf{v}) - p \nabla \cdot \mathbf{v} + \nabla \cdot (\chi \nabla T) + \rho \varepsilon_{\text{nuc}}, \qquad (3) $$

$$ \frac{\partial \rho Y}{\partial t} = -\nabla \cdot (\rho Y \mathbf{v}), \qquad (4) $$

where $\rho$ is the fluid density, $e$ the specific internal energy, $\mathbf{v}$ the fluid velocity, and $Y$ the helium mass fraction. The gravitational acceleration $g = -g(r)\hat{r}$ is purely radial and updated after each time-step from the instantaneous density profile $\bar{\rho}(r)$:

$$ g(r) = 4\pi \frac{G}{r^2} \int_0^r \bar{\rho}(u)u^2 du. \qquad (5) $$

The density profile $\bar{\rho}$ is calculated as

$$ \bar{\rho} = \langle \rho \rangle_{\theta}, \qquad (6) $$

where $\langle \bullet \rangle_\theta$ denotes an angular average across the domain accounting for sphericity

$$ \langle f \rangle_\theta = \frac{1}{2} \int_0^\pi f \sin \theta \, d\theta. \qquad (7) $$

We recover the pressure $p$ and temperature $T$ from $\rho$ and $e$ using tabulated equation of state, which for the simulated regime of 5 M$_\odot$ main-sequence stars are the OPAL tables from Rogers & Nayfonov (2002). Radiative transfer is accounted for via a diffusive formulation involving the heat conductivity $\chi$. The latter is derived from the Rosseland mean opacity $\kappa$ given by the OPAL opacity tables (Iglesias & Rogers 1996):

$$ \chi = \frac{16\sigma T^3}{3\kappa\rho} \qquad (8) $$

(with $\sigma$ 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, $\varepsilon_{\text{nuc}}$ is the specific energy released by nuclear burning. Its radial profile is taken from the 1D model (see Fig. 1), and considered constant in time throughout the 2D run since nuclear burning evolves on much longer time-scales than the duration of the simulations. This is also why no variation of $Y$ due to nuclear reactions appears in equation (4).

All the simulations performed in this study are 2D in $(r, \theta)$ space, assuming $v_\phi = 0$ and invariance along $\phi$ (see Table 2 and hereafter for details). For all three cases, we solve in a domain that spans the entire latitudinal range $(\theta \in [0, \pi])$, and truncate along the radial direction $(r \in [R_{\text{in}}, R_{\text{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 $R_{\text{conv}}$): $H_{p,\text{conv}} \sim 140\Delta r$, with $\Delta r$ the grid size in radial direction. With this resolution, the jump in helium mass fraction above the core spans across $\sim 40\Delta r$ in the MidMS model and $\sim 80\Delta r$ in the TAMS model. We use reflective boundary conditions at $\theta = 0$ and $\theta = \pi$. At $R_{\text{in}}$ and $R_{\text{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); time-steps are therefore not limited by stability considerations. The time-step is set based on limits on both the acoustic and advective CFL. We pick the maximum time-step $\Delta t$ that satisfies both $\Delta t \frac{c_s+|v|}{\Delta r} \le 100$ and $\Delta t \frac{|v|}{\Delta r} \le 0.1$ for every velocity component $v_i$ everywhere in the simulation domain, where $c_s$ 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 time-steps 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 time-scale $\tau_{\text{conv}}$ defined as

$$ \tau_{\text{conv}} = \left\langle \int_{R_{\text{in}}}^{R_{\text{conv}}} \frac{dr}{v_{\text{rms}}(r, t)} \right\rangle_t \qquad (9) $$

where $\langle \bullet \rangle_t$ denotes a time average across the steady state. $v_{\text{rms}}$ is the instantaneous root-mean-square velocity profile:

$$ v_{\text{rms}} = \sqrt{\langle \mathbf{v} \cdot \mathbf{v} \rangle_{\theta}}. \qquad (10) $$

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 time-scale is the convective frequency

$$ f_{\text{conv}} = \frac{1}{\tau_{\text{conv}}}. \qquad (11) $$

Notice in 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 IGWs 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.
## 2.3 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. (2021, 2023), of which an overview is done here. This method is built around two different fluxes: the radial kinetic energy flux $f_k$

$$f_k(r, \theta, t) = \frac{1}{2} \rho v^2 v_r \quad (12)$$

and the radial heat flux $f_{\delta T}$

$$f_{\delta T}(r, \theta, t) = \rho c_p \delta T v_r. \quad (13)$$

The temperature fluctuations $\delta T$ are the departure from the mean temperature profile of the steady state:

$$\delta T = T - \langle\langle T \rangle\rangle_t \quad (14)$$

The overshooting distance $l_k$ with respect to $f_k$ is then defined as

$$l_k(\theta, t) = r_k - R_{\text{conv}}, \quad (15)$$

where $r_k$ is the first radius above $R_{\text{conv}}$ at which $f_k = 0$. Similarly, $l_{\delta T}(\theta, t)$ is defined from the zero-crossing of $f_{\delta T}$. For such low-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 $R_{\text{conv}}$ in equation (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 $l_k(\theta, t)$:

$$l_k^{\max}(t) = \max_{\theta} l_k(\theta, t) \quad (16)$$

and, following Baraffe et al. (2023), its time average defines the overshooting length

$$l_k^{\text{ov}} = \langle l_k^{\max} \rangle_t \quad (17)$$

and similarly for $l_{\delta T}^{\max}$ and $l_{\delta T}^{\text{ov}}$. Several hundreds of overturn times (see Table 2) are needed for $l_k^{\text{ov}}$ and $l_{\delta T}^{\text{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) and Baraffe et al. (2023).

## 2.4 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 IGWs in the envelope. This follows the work presented in Le Saux et al. (2022, 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-of-equilibrium (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 $N$ (in rad s⁻¹) is the maximum frequency of such waves. It is directly related to the density stratification of the medium:

$$N^2 = g \left( \frac{1}{\Gamma_1} \frac{\mathrm{d}\ln p}{\mathrm{d}r} - \frac{\mathrm{d}\ln \rho}{\mathrm{d}r} \right). \quad (18)$$

$\Gamma_1$ is the first adiabatic exponent:

$$\Gamma_1 = \left. \frac{\partial \ln p}{\partial \ln \rho} \right|_{\text{ad}}. \quad (19)$$

A fluid is unstable against convection if $N^2 \le 0$, in which case IGW are evanescent. Conversely, a fluid is stably stratified if $N^2 > 0$, and IGW of angular frequency $\omega \le N$ can propagate in the medium. Note that in this paper, $\omega$ denotes the angular frequency in rad s⁻¹, and $f = \omega/(2\pi)$ is the frequency in hertz. In the absence of rotation, the dispersion relation of IGW verifies (Press 1981):

$$\frac{\omega^2}{N^2} = \frac{k_h^2}{k_r^2 + k_h^2} \quad (20)$$

where $k_r$ is the radial wavenumber and $k_h$ the horizontal wavenumber. The latter is related to the spherical harmonic degree $l$ of the considered wave:

$$k_h^2 = \frac{l(l+1)}{r^2}. \quad (21)$$

The linear theory developed in Press (1981) and extended to variable molecular weight $\mu$ in Zahn, Talon & Matias (1997) uses the Wentzel–Kramers–Brillouin (WKB) approximation to predict the amplitude of travelling IGW in the radiative region. This yields the following expression for the radial velocity:

$$\hat{v}_r(r, l, \omega) = \hat{v}_e(l, \omega) \left(\frac{\rho_e}{\rho}\right)^{1/2} \left(\frac{r_e}{r}\right)^{3/2} \left(\frac{N_e^2 - \omega^2}{N^2 - \omega^2}\right)^{1/4} \times \exp(-\tau/2), \quad (22)$$

where $e$-indices denote quantities evaluated at the excitation radius. Note that $\hat{v}_e(l, \omega)$ is an arbitrary amplitude for the linear theory. In equation (22), $\tau$ 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. $\tau$ verifies (Zahn et al. 1997):

$$\tau(r, l, \omega) = \int_{r_e}^{r} \kappa_{\text{rad}} k_h^3 \frac{N N_t^2}{\omega^4} \left( \frac{N^2}{N^2 - \omega^2} \right)^{1/2} dr. \quad (23)$$

$\kappa_{\text{rad}}$ is the radiative diffusivity

$$\kappa_{\text{rad}} = \frac{\chi}{\rho c_p}, \quad (24)$$

where $c_p$ is the specific heat capacity as constant pressure. $N_t$ is the thermal contribution to the angular Brunt–Väisälä frequency

$$N_t^2 = - \left. \frac{\partial \ln \rho}{\partial \ln T} \right|_{p,\mu} \cdot \frac{g}{H_p} \cdot \left( \left. \frac{\partial \ln T}{\partial \ln p} \right|_{\text{ad}} - \frac{\mathrm{d}\ln T}{\mathrm{d}\ln p} \right). \quad (25)$$

A related quantity is the wave luminosity, defined as (Press 1981; Zahn et al. 1997):

$$L_{\text{wave}} = 4\pi r^2 \rho v_{\text{rms}}^2 V_{g,r}, \quad (26)$$

where $V_{g,r}$ is the radial group velocity

$$V_{g,r} = \frac{\omega^2 \sqrt{N^2 - \omega^2}}{N^2 k_h}. \quad (27)$$

This luminosity is only affected by damping as waves travel in the radiative region:

$$L_{\text{wave}}(r, l, \omega) = L_{\text{wave}}(r_e, l, \omega) \exp(-\tau). \quad (28)$$

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. $P[\hat{v}_r](r, l, f)$ 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 $v_h/v_r \simeq k_r/k_h$ and the dispersion relation equation (20), the
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Temperature perturbation maps for the ZAMS (left) and MidMS (right) cases. The solid black line marks the position of the convective boundary.</figcaption></figure>

wave luminosity equation (26) can then be expressed from this power spectrum:

$$L_{\text{wave}} \approx 4\pi r^2 \rho \frac{\sqrt{N^2 - (2\pi f)^2}}{k_{\text{h}}} P[\hat{v}_r]. \quad (29)$$

The linear expression on $\hat{v}_r$ equation (22) also applies to $\sqrt{P[\hat{v}_r]}$ since $\hat{v}_r = \sqrt{P[\hat{v}_r]}/2$. For simplicity, we test the linear prediction on profiles of $\sqrt{P[\hat{v}_r]}$ directly. The excitation radius $r_e$ in equation (22) is taken slightly above the convective core boundary, at twice the overshooting length equation (17):

$$r_e = R_{\text{conv}} + 2l_k^{\text{ov}}. \quad (30)$$

Note that taking $r_e$ at $l_k^{\text{ov}}$ or $2l_k^{\text{ov}}$ above the convective boundary does not change our results significantly. To improve the readability of figures, we choose $\hat{v}_e(\ell, \omega)$ in equation (22) as the local maximum of $\sqrt{P[\hat{v}_r]}$ evaluated around $r_e$. 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.

# 3 RESULTS

## 3.1 General observations on the dynamics

Fig. 2 shows the instantaneous temperature perturbation

$$\epsilon_T = \frac{T - \langle T \rangle_{\theta}}{\langle T \rangle_{\theta}} \quad (31)$$

for the ZAMS and MidMS cases. The TAMS case exhibits similar behaviour to the MidMS case. The convective part is extremely well mixed with $|\epsilon_T| \lesssim 10^{-6}$ while the stably stratified region shows stronger deviations from the mean profile ($|\epsilon_T| \sim 10^{-4}$) that organize in well-defined patterns. Those patterns are the signature of IGWs. 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 $\epsilon_T$ are clearly visible across the stable region of the MidMS model.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Kinetic energy density profiles for the three cases, total $\langle\langle\frac{1}{2}\rho v^2\rangle\rangle_{\theta}$ and radial $\langle\langle\frac{1}{2}\rho v_r^2\rangle\rangle_{\theta}$. The solid vertical lines mark the position of the convective core boundary according to the Schwarzschild criterion.</figcaption></figure>

Sitting atop the core is the $N^2$-peak region, showing patterns of $\epsilon_T$ with larger angular and radial extent than in the rest of the stable region at lower $N^2$. Section 3.3 presents an in-depth analysis of the wave content in those simulations.

As can be seen in Fig. 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 $N^2$-peak layer of evolved models reduces the transfer of kinetic energy from the convective core to the radiative envelope. Fig. 3 provides an additional insight on the effects of the $N^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 $N^2$-peak layer of evolved models compared to the rest of the stable region. This is readily explained by the dispersion relation equation (20) which implies that $v^2/v_r^2 \simeq N^2/\omega^2$; as $N^2$ increases, the horizontal component of the kinetic energy density will therefore carry more of the energy.

## 3.2 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 $l_k^{\text{ov}}$ and $l_{\delta T}^{\text{ov}}$ are shown in Table 3 for the three models. For the ZAMS case, one recovers the value found in Baraffe et al. (2023) for the 5L0 case ($l_k \sim 0.09 H_p$, $l_{\delta T} \sim 0.06 H_p$). For both evolved cases MidMS and TAMS, the penetration length associated with either flux is much lower $l \sim 0.025 H_p$. It should be noted that $l(\theta, t)$ is only one grid cell $\Delta r \sim 0.007 H_p$ in many areas, meaning the penetration length is likely underresolved and possibly overestimated for those evolved
<table><caption>**Table 3.** Overshooting length measured for the three evolution stages.</caption><thead><tr><th>Model</th><th><i>l</i><sub>k</sub><sup>ov</sup> (cm)</th><th><i>l</i><sub>δT</sub><sup>ov</sup> (cm)</th><th><i>l</i><sub>k</sub><sup>ov</sup>/<i>H</i><sub>p,conv</sub></th><th><i>l</i><sub>δT</sub><sup>ov</sup>/<i>H</i><sub>p,conv</sub></th></tr></thead><tbody><tr><td>ZAMS</td><td>1.53 × 10<sup>9</sup></td><td>1.51 × 10<sup>9</sup></td><td>0.084</td><td>0.083</td></tr><tr><td>MidMS</td><td>3.91 × 10<sup>8</sup></td><td>4.86 × 10<sup>8</sup></td><td>0.022</td><td>0.028</td></tr><tr><td>TAMS</td><td>3.65 × 10<sup>8</sup></td><td>4.64 × 10<sup>8</sup></td><td>0.022</td><td>0.028</td></tr></tbody></table>

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.

### 3.3 Wave study

The previous subsection shows that the $N^2$-peak layer hinders radial motions in the evolved models. This subsection focuses on the consequences on IGW in the radiative region. Fig. 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, Christensen-Dalsgaard & Kurtz 2010). A notable feature of both evolved models compared to the ZAMS model is that nodes are narrowed in the $N^2$-peak layer in evolved models. This is, again, explained by the dispersion relation equation (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. Fig. 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 Fig. 6).

Fig. 5 compares power spectra at ($ℓ = 5, f = 7.6$ µHz) and ($ℓ = 10, f = 1.2$ µHz) with the linear theory equation (22) for all three models, with damping following equation (23) and without damping ($τ = 0$). Note that since we use realistic values for the heat conductivity $χ$ (see equation (8)), we expect our 2D models to exhibit realistic damping profiles. The ($ℓ = 5, f = 7.6$ µHz) mode undergoes very limited damping in the simulated region, and as can be seen in Fig. 5, the simulated spectrum matches well the linear theory for all three models. On the other hand, the ($ℓ = 10, f = 1.2$ µHz) mode is predicted to be severely damped (which is a consequence of the higher $ℓ$ and lower $f$ 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 $N^2$-peak. The TAMS model is the one that matches the most the linear prediction for this mode, with a good match in the $N^2$-peak region and roughly following the predicted slope in the rest of the stratified region. The increase in amplitude when exiting the $N^2$-peak layer is systematically higher in the simulations than what is predicted by the linear theory. This linear theory is developed using the WKB approximation, and therefore assumes that spatial derivatives of various parameters (such as $N^2$) are small;

<figure><img src="image_5.png" /><figcaption>**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.</figcaption></figure>

this assumption is potentially not valid at the top of the $N^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 Fig. 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 $λ_r$ from the dispersion relation equation (20) is well above the cell width $Δr$ for the mode ($ℓ = 5, f = 7.6$ µHz) across all models ($6Δr < λ_r < 60Δr$), but the radial wavelength for the mode ($ℓ = 10, f = 1.2$ µHz) is comparable to $2Δr$ for all three models and even smaller in the $N^2$-peak region: the latter mode is therefore underresolved. It should finally be noted that the linear theory equation (22) concerns travelling waves, while
<figure><img src="image_6.png" /><figcaption>**Figure 5.** 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 equation (22), and the dotted line is the linear prediction without radiative damping. These two lines are superimposed for the ($l = 5, f = 7.6$ µHz) mode as damping is small. The black dashed line marks the convective boundary $R_{conv}$.</figcaption></figure>

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 $N^2$ at the top of the $N^2$-peak region. Reflection of waves in the overall domain and/or within the $N^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.

Fig. 6 compares the wave power spectra at $l = 3$ for the three stellar models at two radii: at the excitation radius equation (30), and $0.07R_*$ above it (outside of the $N^2$-peak region in evolved models). The damping at low $f \lesssim 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 around 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.

Fig. 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 (Fig. 3).

As a check, we performed an extended simulation of the TAMS model (with $R_{out} = 0.59R_*$); Fig. 7 shows the obtained wave luminosity spectrum and comparison with the linear theory for the ($l = 15, f = 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 $r \lesssim 0.2R_*$ (with, again, a stronger increase in amplitude than predicted when exiting the $N^2$-peak layer at $r \simeq 0.13R_*$), but strong oscillations due to g-modes renders the comparison to the theory difficult at higher radius.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Wave luminosity spectrum profiles at $ℓ = 3$ for the three cases, at two different radii above the convective boundary. $r_e$ is slightly above the convective boundary (at twice the penetration length $l_k^{\text{conv}}$ found in Section 3.2), and $r_e + 0.07R_*$ is outside of the $N^2$-peak region of the MidMS and TAMS models. The thin light-coloured solid lines are the spectra obtained from the simulations, and the thick dark-coloured solid lines are the same profiles smoothed for readability purposes. The vertical grey dash-dotted line is the convective frequency. Dotted lines are gaussian models of excitation by plume penetration (Pinçon, Belkacem & Goupil 2016), and dashed lines are scalings considering Reynolds stress (Lecoanet & Quataert 2013).</figcaption></figure>

Fig. 8 shows the theoretical wave power spectrum for the three evolution stages, using the 1D initial profiles which extend up to $r = R_*$. 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 Anders et al. (2023) and Le Saux 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.

## 4 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; IGWs 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 $N^2$-peak region, reflections at the top of the $N^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 strengthens the conclusions from Anders et al. (2023) and Le Saux et al. (2023). Our results clearly highlight that the amplitude of IGWs propagating in the radiative envelope of intermediate mass stars is expected to be small. It is thus even less likely to detect IGWs excited by core convection in these stars, casting doubts on the interpretation of Bowman et al. (2019, 2020).

Finally, Baraffe et al. (2023) show that 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 $M \gtrsim 10 M_\odot$. 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 Ratnasingam et al. (2023), Vanon et al. (2023), and Varghese et al. (2023). They present simulations in similar set-ups, 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 a helium mass fraction field, and with modified heat fluxes, diffusivities, and temperature profiles compared to stellar evolution models to mimic the $N^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 $N^2$-peak via a modified temperature profile leads to a more realistic damping profile in the $N^2$-peak region as we do not conflate $N_t$ and $N$ in equation (23).

There are some differences between our findings and those of previous studies. Ratnasingam et al. (2023) and Vanon et al. (2023) find that their power spectra have a better agreement with plume excitation, while we find a better agreement with Reynolds stress excitation (Fig. 6). Vanon et al. (2023) report discrepancies between their TAMS simulations and theoretical predictions of g modes with GYRE. They suggest that 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, Vanon et al. (2023) and Varghese et al. (2023) find that IGWs 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 studies. A comparison with the
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Wave luminosity spectrum for the extended TAMS model (top), and comparison with the linear theory for the ($\ell = 15, f = 2.5$ µHz) mode (bottom). Top: as in Fig. 4, the dashed lines denote the g modes predicted by GYRE. Labels for nodes −8 to −15 are omitted for readability. Bottom: as in Fig. 5, the light solid line is the spectrum obtained from the simulation, the dark solid line is the linear prediction with radiative damping equation (22), and the dotted line is the linear prediction without radiative damping.</figcaption></figure>

recent work by Thompson et al. (2024) based on 3D simulations of a 25 M⊙ is also difficult since it does not include radiative diffusion. As we show in this 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 to the 2D ones from Rogers et al. (2013), suggesting that the power spectra we show in this 2D study should be largely unaffected by dimensionality. Convective velocities in 2D simulations are usually larger than those in 3D simulations (Meakin & Arnett 2007; Pratt et al. 2020), suggesting larger overshooting lengths predicted by 2D simulations compared to 3D models. However, the plume filling

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Wave luminosity spectrum calculated from the linear theory equation (28) for all three models at $\ell = 3$. The dashed line marks the convective frequency $f_{conv}$.</figcaption></figure>

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, Clune & Toomre 2002; Rogers, Glatzmaier & Jones 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 those of 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., in preparation). 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 M⊙ 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.

## ACKNOWLEDGEMENTS

The authors thank an anonymous reviewer for valuable comments that helped clarify this paper. This work was supported by the ERC grant no. 787361-COBOM and the STFC consolidated grant ST/V000721/1. The authors would like to acknowledge the use of the University of Exeter High-Performance Computing (HPC) facility ISCA and of the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility. The equipment was funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1, and STFC DiRAC operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure. Part of this work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under contract no. DE-AC52-07NA27344 (LLNL-JRNL-859808).

## DATA AVAILABILITY

Data backing this study will be shared on reasonable request to the corresponding author.

## REFERENCES

* Aerts C., Christensen-Dalsgaard J., Kurtz D. W., 2010, Asteroseismology. Springer, Dordrecht, p. 866

* Anders E. H. et al., 2023, Nat. Astron., 7, 1228

* Baraffe I., El Eid M. F., 1991, A&A, 245, 548

* Baraffe I., Chabrier G., Allard F., Hauschildt P. H., 1998, A&A, 337, 403

* Baraffe I., Pratt J., Vlaykov D. G., Guillet T., Goffrey T., Le Saux A., Constantino T., 2021, A&A, 654, A126

* Baraffe I. et al., 2023, MNRAS, 519, 5333

* Bossini D. et al., 2015, MNRAS, 453, 2290

* Bowman D. M. et al., 2019, A&A, 621, A135

* Bowman D. M., Burssens S., Simón-Díaz S., Edelmann P. V. F., Rogers T. M., Horst L., Röpke F. K., Aerts C., 2020, A&A, 640, A36

* Brummell N. H., Clune T. L., Toomre J., 2002, ApJ, 570, 825

* Claret A., Torres G., 2016, A&A, 592, A15

* Edelmann P. V. F., Ratnasingam R. P., Pedersen M. G., Bowman D. M., Prat V., Rogers T. M., 2019, ApJ, 876, 4

* Goffrey T. et al., 2017, A&A, 600, A7

* Iglesias C. A., Rogers F. J., 1996, ApJ, 464, 943

* Le Saux A. et al., 2022, A&A, 660, A51

* Le Saux A., Baraffe I., Guillet T., Vlaykov D. G., Morison A., Pratt J., Constantino T., Goffrey T., 2023, MNRAS, 522, 2835

* Lecoanet D., Quataert E., 2013, MNRAS, 430, 2363

* Lecoanet D. et al., 2019, ApJ, 886, L15

* Lecoanet D., Cantiello M., Anders E. H., Quataert E., Couston L.-A., Bouffard M., Favier B., Le Bars M., 2021, MNRAS, 508, 132

* Meakin C. A., Arnett D., 2007, ApJ, 667, 448

* Miglio A., Montalbán J., Noels A., Eggenberger P., 2008, MNRAS, 386, 1487

* Pinçon C., Belkacem K., Goupil M. J., 2016, A&A, 588, A122

* Pratt J., Baraffe I., Goffrey T., Constantino T., Viallet M., Popov M. V., Walder R., Folini D., 2017, A&A, 604, A125

* Pratt J., Baraffe I., Goffrey T., Geroux C., Constantino T., Folini D., Walder R., 2020, A&A, 638, A15

* Press W. H., 1981, ApJ, 245, 286

* Ratnasingam R. P., Rogers T. M., Chowdhury S., Handler G., Vanon R., Varghese A., Edelmann P. V. F., 2023, A&A, 674, A134

* Rogers F. J., Nayfonov A., 2002, ApJ, 576, 1064

* Rogers T. M., Glatzmaier G. A., Jones C. A., 2006, ApJ, 653, 765

* Rogers T. M., Lin D. N. C., McElwaine J. N., Lau H. H. B., 2013, ApJ, 772, 21

* Rosenfield P. et al., 2017, ApJ, 841, 69

* Schatzman E., 1993, A&A, 279, 431

* Shaviv G., Salpeter E. E., 1973, ApJ, 184, 191

* Stancliffe R. J., Fossati L., Passy J.-C., Schneider F. R. N., 2016, A&A, 586, A119

* Takata M., 2006, Publ. Astron. Soc. Jpn., 58, 893

* Thompson W., Herwig F., Woodward P. R., Mao H., Denissenkov P., Bowman D. M., Blouin S., 2024, MNRAS, 531, 1316

* Townsend R. H. D., Teitler S. A., 2013, MNRAS, 435, 3406

* Vanon R., Edelmann P. V. F., Ratnasingam R. P., Varghese A., Rogers T. M., 2023, ApJ, 954, 171

* Varghese A., Ratnasingam R. P., Vanon R., Edelmann P. V. F., Rogers T. M., 2023, ApJ, 942, 53

* Viallet M., Goffrey T., Baraffe I., Folini D., Geroux C., Popov M. V., Pratt J., Walder R., 2016, A&A, 586, A153

* Zahn J. P., 1991, A&A, 252, 179

* Zahn J. P., 1992, A&A, 265, 115

* Zahn J.-P., Talon S., Matias J., 1997, A&A, 322, 320

This paper has been typeset from a TeX/LATEX file prepared by the author.