# Two-dimensional simulations of internal gravity waves in a 5 $M_\odot$ Zero-Age-Main-Sequence model

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

Constantino, T Goffrey

► To cite this version:

A. Le Saux, I Baraffe, T Guillet, D.G Vlaykov, A Morison, et al. Two-dimensional simulations of internal gravity waves in a 5 $M_\odot$ Zero-Age-Main-Sequence model. Monthly Notices of the Royal Astronomical Society, 2023, 522 (2), pp.2835-2849. 10.1093/mnras/stad1067. hal-04084071

HAL Id: hal-04084071

https://hal.science/hal-04084071

Submitted on 24 May 2024

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
<figure><img src="image_1.png" /></figure>

# Two-dimensional simulations of internal gravity waves in a 5 M⊙ zero-age-main-sequence model

A. Le Saux<sup>1,2*</sup> I. Baraffe,<sup>1,2</sup> T. Guillet<sup>1</sup>,<sup>1</sup> D. G. Vlaykov<sup>1</sup>,<sup>1</sup> A. Morison,<sup>1</sup> J. Pratt,<sup>3</sup> T. Constantino<sup>1</sup>
and T. Goffrey<sup>4</sup>

Accepted 2023 April 4. Received 2023 April 4; in original form 2023 January 31

## ABSTRACT

Main-sequence intermediate-mass stars present a radiative envelope that supports internal gravity waves (IGWs). Excited at the boundary with the convective core, IGWs propagate towards the stellar surface and are suspected to impact physical processes such as rotation and chemical mixing. Using the fully compressible time-implicit code MUSIC, we study IGWs in two-dimensional simulations of a zero-age-main-sequence 5 solar mass star model up to 91 per cent of the stellar radius with different luminosity and radiative diffusivity enhancements. Our results show that low-frequency waves excited by core convection are strongly impacted by radiative effects as they propagate. This impact depends on the radial profile of radiative diffusivity which increases by almost 5 orders of magnitude between the centre of the star and the top of the simulation domain. In the upper layers of the simulation domain, we observe an increase of the temperature. Our study suggests that this is due to heat added in these layers by IGWs damped by radiative diffusion. We show that non-linear effects linked to large amplitude IGWs may be relevant just above the convective core. Both these effects are intensified by the artificial enhancement of the luminosity and radiative diffusivity, with enhancement factors up to $10^4$ times the realistic values. Our results also highlight that direct comparison between numerical simulations with enhanced luminosity and observations must be made with caution. Finally, our work suggests that thermal effects linked to the damping of IGWs could have a non-negligible impact on stellar structure.

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

## 1 INTRODUCTION

Studying the properties of internal gravity waves (IGWs) in stars and planets is of crucial importance as they are supposed to be involved in many physical processes linked to rotation, mixing and magnetism. In the Earth atmosphere, it is generally accepted that they are particularly efficient at transporting angular momentum (AM), energy and chemical elements (see Sutherland 2010, for a review). Therefore, stellar physicists have considered their existence in stars since the 1960s (Stein 1967) and studied their properties and how they could be excited in this context. Stein (1967) extended the pioneering work on acoustic waves by Lighthill (1952) to a stratified medium, and concluded that turbulent convection can generate IGWs that will propagate in an adjacent stably stratified region. In this work, we consider IGWs excited by this mechanism.

Since then, IGWs have been considered to be involved in various phenomena in stellar interiors. For instance, they have been suggested as a possible explanation for the solid rotation of the solar core inferred from helioseismology (see e.g. Charbonnel & Talon 2005),

the mixing profile in the radiative zones of intermediate-mass stars induced from asteroseismology (Pedersen et al. 2021) or the slowing down of the core rotation of subgiant stars (Pinçon et al. 2017) as observed by the *Kepler* mission (Borucki et al. 2010).

However, the properties of IGWs excited by turbulent convection in stellar interiors remain poorly constrained due to the challenge of their observations. There have been claims of detection in the Sun (see e.g. Appourchaux et al. 2010) and in more massive stars (see e.g. Aerts & Rogers 2015), but these are still a matter of debate. Observations of these waves are essential to probe the deeper layers of stars. In the Sun, they will allow the measurement of the rotation profile of the inner 20 per cent of the solar radius that remain inaccessible to acoustic modes (García et al. 2007). In more massive stars, they would bring new constraints on internal rotation, mixing, and convective boundary mixing, which has been a long-standing puzzle of stellar structure (see e.g. Renzini 1987). These more massive stars with a convective core are probably better targets for observations of IGWs.

Unlike main-sequence low-mass stars where the radiative zone is located in the inner part of the star, more massive main-sequence stars with $M \gtrsim 2M_\odot$ present a radiative envelope. Internal gravity

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

<footnote n="2"><sup>2</sup>École Normale Supérieure, Lyon, CRAL (UMR CNRS 5574), Université de Lyon, 69007 Lyon, France</footnote>

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

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

<footnote>*E-mail: al598@exeter.ac.uk</footnote>
<figure><img src="image_2.png" /></figure>

# Two-dimensional simulations of internal gravity waves in a 5 M⊙ zero-age-main-sequence model

A. Le Saux<sup>1,2*</sup> I. Baraffe,<sup>1,2</sup> T. Guillet<sup>1</sup>,<sup>1</sup> D. G. Vlaykov<sup>1</sup>,<sup>1</sup> A. Morison,<sup>1</sup> J. Pratt,<sup>3</sup> T. Constantino<sup>1</sup> and T. Goffrey<sup>4</sup>

Accepted 2023 April 4. Received 2023 April 4; in original form 2023 January 31

## ABSTRACT

Main-sequence intermediate-mass stars present a radiative envelope that supports internal gravity waves (IGWs). Excited at the boundary with the convective core, IGWs propagate towards the stellar surface and are suspected to impact physical processes such as rotation and chemical mixing. Using the fully compressible time-implicit code MUSIC, we study IGWs in two-dimensional simulations of a zero-age-main-sequence 5 solar mass star model up to 91 per cent of the stellar radius with different luminosity and radiative diffusivity enhancements. Our results show that low-frequency waves excited by core convection are strongly impacted by radiative effects as they propagate. This impact depends on the radial profile of radiative diffusivity which increases by almost 5 orders of magnitude between the centre of the star and the top of the simulation domain. In the upper layers of the simulation domain, we observe an increase of the temperature. Our study suggests that this is due to heat added in these layers by IGWs damped by radiative diffusion. We show that non-linear effects linked to large amplitude IGWs may be relevant just above the convective core. Both these effects are intensified by the artificial enhancement of the luminosity and radiative diffusivity, with enhancement factors up to $10^4$ times the realistic values. Our results also highlight that direct comparison between numerical simulations with enhanced luminosity and observations must be made with caution. Finally, our work suggests that thermal effects linked to the damping of IGWs could have a non-negligible impact on stellar structure.

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

## 1 INTRODUCTION

Studying the properties of internal gravity waves (IGWs) in stars and planets is of crucial importance as they are supposed to be involved in many physical processes linked to rotation, mixing and magnetism. In the Earth atmosphere, it is generally accepted that they are particularly efficient at transporting angular momentum (AM), energy and chemical elements (see Sutherland 2010, for a review). Therefore, stellar physicists have considered their existence in stars since the 1960s (Stein 1967) and studied their properties and how they could be excited in this context. Stein (1967) extended the pioneering work on acoustic waves by Lighthill (1952) to a stratified medium, and concluded that turbulent convection can generate IGWs that will propagate in an adjacent stably stratified region. In this work, we consider IGWs excited by this mechanism.

Since then, IGWs have been considered to be involved in various phenomena in stellar interiors. For instance, they have been suggested as a possible explanation for the solid rotation of the solar core inferred from helioseismology (see e.g. Charbonnel & Talon 2005),

the mixing profile in the radiative zones of intermediate-mass stars induced from asteroseismology (Pedersen et al. 2021) or the slowing down of the core rotation of subgiant stars (Pinçon et al. 2017) as observed by the *Kepler* mission (Borucki et al. 2010).

However, the properties of IGWs excited by turbulent convection in stellar interiors remain poorly constrained due to the challenge of their observations. There have been claims of detection in the Sun (see e.g. Appourchaux et al. 2010) and in more massive stars (see e.g. Aerts & Rogers 2015), but these are still a matter of debate. Observations of these waves are essential to probe the deeper layers of stars. In the Sun, they will allow the measurement of the rotation profile of the inner 20 per cent of the solar radius that remain inaccessible to acoustic modes (García et al. 2007). In more massive stars, they would bring new constraints on internal rotation, mixing, and convective boundary mixing, which has been a long-standing puzzle of stellar structure (see e.g. Renzini 1987). These more massive stars with a convective core are probably better targets for observations of IGWs.

Unlike main-sequence low-mass stars where the radiative zone is located in the inner part of the star, more massive main-sequence stars with $M \gtrsim 2M_\odot$ present a radiative envelope. Internal gravity

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

<footnote n="2"><sup>2</sup>École Normale Supérieure, Lyon, CRAL (UMR CNRS 5574), Université de Lyon, 69007 Lyon, France</footnote>

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

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

<footnote>*E-mail: al598@exeter.ac.uk</footnote>
waves generated at the edge of the convective core propagate towards the surface through a medium of decreasing density, which tend to increase their amplitude. These waves are also damped by radiative diffusion as they travel, and are suggested to transport energy and angular momentum and mix chemical elements through this mechanism (Schatzman 1993). Thus, the evolution of IGWs amplitude will depend on the interplay between growth due to decreasing density and decay due to radiative damping. It remains unclear if these waves should be able to propagate up to the surface. More than a decade ago, Blomme et al. (2011) observed a low frequency power excess in the spectra of O type stars observed by the CoRoT (Auvergne et al. 2009) mission. They concluded that the physical origin of this power excess was unclear. In recent studies, (Bowman et al. 2019, 2020) claim that this low frequency power excess in the spectrum of O and B type stars observed by the CoRoT and TESS (Ricker et al. 2015) satellites is due to IGWs excited by turbulent core convection. This hypothesis is supported by hydrodynamical simulations in 3D (Edelmann et al. 2019) and 2D (Ratnasingam, Edelmann & Rogers 2020). However, theoretical work by Lecoanet et al. (2019) and numerical simulations by Lecoanet et al. (2021) do not agree with this conclusion and state that the origin of this power excess is more probably due to a near surface convection zone (Cantiello et al. 2021). This question will remain difficult to answer while the properties of IGWs in such stars remain poorly known.

Hydrodynamical simulations offer a great opportunity to test theoretical models and guide observations. Particularly numerical modelling of internal waves, is a good way to get constraints on the spectrum generated by convection, their amplitude and damping rate within the star. This has already proven to be efficient in simulations of solar-like stars (see e.g. Rogers & Glatzmaier 2005; Alvan, Brun & Mathis 2014) as well as stars with convective cores (see e.g. Rogers et al. 2013; Horst et al. 2020). In this work, we present two-dimensional (2D) simulations of a 5 solar mass star with different enhancement factors for the luminosity and the radiative diffusivity, as it is widely used for hydrodynamical simulations. This is usually referred to as *boosting* a numerical model. In this study, we focus mainly on the damping of IGWs and the impact of boosting on these waves, and we neglect the effect of rotation, even though most OB stars may be moderate to fast rotators. In Section 2, we present the numerical set-up of our simulations. Then Section 3 and 4 focus on amplitude related properties of IGWs, whereas Section 5 studies the non-linearity of waves and Section 6 studies an increase in the temperature observed in the upper layers. Finally, in Section 7 we summarize and discuss our results.

## 2 NUMERICAL SIMULATIONS

We performed 2D simulations of the interior of a 5 M$_{\odot}$ star model at zero-age-main-sequence (ZAMS) with the fully compressible time-implicit code MUSIC (Viallet, Baraffe & Walder 2011, 2013; Geroux et al. 2016; Pratt et al. 2016; Viallet et al. 2016; Goffrey et al. 2017; Pratt et al. 2017, 2020). The code solves the inviscid hydrodynamical equations in a fully compressible medium

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

$$ \frac{\partial \rho e}{\partial t} = -\nabla \cdot (\rho e \mathbf{v}) - p \nabla \cdot \mathbf{v} - \nabla \cdot \mathbf{F}_r + Q_{\text{nuc}}, \qquad (2) $$

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

where $\rho$ is the density, $e$ is the specific internal energy, $\mathbf{v}$ is the velocity field, $p$ the gas pressure, $Q_{\text{nuc}}$ is the nuclear energy rate,

<table><caption>**Table 1.** Properties of the initial stellar model used for the 2D hydrodynamical simulations.</caption><thead><tr><th>M/M<sub>☉</sub></th><th>L<sub>star</sub>/L<sub>☉</sub><sup>a</sup></th><th>R<sub>star</sub> (cm)</th><th>r<sub>conv</sub>/R<sub>tot</sub></th><th>H<sub>P, conv</sub> (cm)</th></tr></thead><tbody><tr><td>5</td><td>523</td><td>1.8424 × 10<sup>11</sup></td><td>0.1814</td><td>1.828 × 10<sup>10</sup></td></tr></tbody></table>

Notes. With mass $M$, luminosity $L_{\text{star}}$, radius $R_{\text{star}}$, size of the convective core $r_{\text{conv}}$, and pressure scale height at the convective boundary $H_{\text{P, conv}}$.

$^{a}$We use $L_{\odot} = 3.839 \times 10^{33}$ erg s$^{-1}$.

and $\mathbf{g}$ is the gravitational acceleration, which we assume is radial. For the stellar simulations considered in this work, the major heat transport that contributes to thermal conductivity is radiative transfer characterized by the radiative flux $\mathbf{F}_r$, given within the diffusion approximation by

$$ \mathbf{F}_r = -\frac{4acT^3}{3\kappa\rho}\nabla T = -\chi\nabla T, \qquad (4) $$

with $\kappa$ is the Rosseland mean opacity of the gas and $\chi$ the radiative conductivity. Realistic stellar opacities and equation of states appropriate for the description of stellar interiors are implemented in MUSIC. Opacities are interpolated from the OPAL tables (Iglesias & Rogers 1996) for solar metallicity, and the equation of state is based on the OPAL EOS tables of Rogers & Nayfonov (2002).

## 2.1 Initial stellar model

The initial structures for the 2D simulations come from a 5 M$_{\odot}$ stellar model computed with the one-dimensional (1D) Lyon stellar evolution code (Baraffe & El Eid 1991; Baraffe et al. 1998), using the same opacities and equation of state as implemented in MUSIC. The initial 1D model has an initial helium abundance in mass fraction $Y = 0.28$ and solar metallicity $Z = 0.02$ and is evolved through the pre-main-sequence (PMS) and the early main-sequence. The initial structure of the simulations have burnt 1 per cent of their hydrogen since the ZAMS. There is no overshooting or diffusion considered during the computation of the 1D model. In the MUSIC simulations, the energy generated by nuclear reactions is taken into account in the energy equation as a source term. The nuclear energy profile is the one from the 1D model, and it is assumed to remain constant during the run of the MUSIC simulations. The properties of the initial 1D model are summarized in Table 1.

## 2.2 Spherical-shell geometry and boundary conditions

A detailed description of the set-up of the simulations can be found in Baraffe et al. (2023). In this section, we only summarize some characteristics relevant to this work. 2D simulations are performed in a spherical shell using spherical coordinates, namely $r$ the radius and the polar angle $\theta$, and assuming azimuthal symmetry in the $\phi$-direction. For all simulations, the inner radius $r_{\text{in}}$ is set at 0.02 $R_{\text{star}}$ and the outer radius $r_{\text{out}}$ at 0.91 $R_{\text{star}}$. Extension of the numerical domain to the photosphere ($r = R_{\text{star}}$) is an open challenge for stellar hydrodynamical simulations given the sharp decrease of the pressure scale height with increasing radius. The domain in the co-latitudinal direction ranges from 0 to $\pi$. We use a uniform grid resolution of $N_r \times N_\theta = 1322 \times 668$ cells. This provides a good resolution of the pressure scale height at the convective boundary $H_{\text{P, conv}}/\Delta r \sim 147$, with $\Delta r = 1.242 \times 10^8$ cm, the size of a radial grid cell.

Concerning the radial boundary conditions, we impose a constant radial derivative on the density on the inner and outer radial boundaries, as discussed in Pratt et al. (2016). For the velocity, we impose reflective conditions at the radial boundaries. For our
waves generated at the edge of the convective core propagate towards the surface through a medium of decreasing density, which tend to increase their amplitude. These waves are also damped by radiative diffusion as they travel, and are suggested to transport energy and angular momentum and mix chemical elements through this mechanism (Schatzman 1993). Thus, the evolution of IGWs amplitude will depend on the interplay between growth due to decreasing density and decay due to radiative damping. It remains unclear if these waves should be able to propagate up to the surface. More than a decade ago, Blomme et al. (2011) observed a low frequency power excess in the spectra of O type stars observed by the CoRoT (Auvergne et al. 2009) mission. They concluded that the physical origin of this power excess was unclear. In recent studies, (Bowman et al. 2019, 2020) claim that this low frequency power excess in the spectrum of O and B type stars observed by the CoRoT and TESS (Ricker et al. 2015) satellites is due to IGWs excited by turbulent core convection. This hypothesis is supported by hydrodynamical simulations in 3D (Edelmann et al. 2019) and 2D (Ratnasingam, Edelmann & Rogers 2020). However, theoretical work by Lecoanet et al. (2019) and numerical simulations by Lecoanet et al. (2021) do not agree with this conclusion and state that the origin of this power excess is more probably due to a near surface convection zone (Cantiello et al. 2021). This question will remain difficult to answer while the properties of IGWs in such stars remain poorly known.

Hydrodynamical simulations offer a great opportunity to test theoretical models and guide observations. Particularly numerical modelling of internal waves, is a good way to get constraints on the spectrum generated by convection, their amplitude and damping rate within the star. This has already proven to be efficient in simulations of solar-like stars (see e.g. Rogers & Glatzmaier 2005; Alvan, Brun & Mathis 2014) as well as stars with convective cores (see e.g. Rogers et al. 2013; Horst et al. 2020). In this work, we present two-dimensional (2D) simulations of a 5 solar mass star with different enhancement factors for the luminosity and the radiative diffusivity, as it is widely used for hydrodynamical simulations. This is usually referred to as *boosting* a numerical model. In this study, we focus mainly on the damping of IGWs and the impact of boosting on these waves, and we neglect the effect of rotation, even though most OB stars may be moderate to fast rotators. In Section 2, we present the numerical set-up of our simulations. Then Section 3 and 4 focus on amplitude related properties of IGWs, whereas Section 5 studies the non-linearity of waves and Section 6 studies an increase in the temperature observed in the upper layers. Finally, in Section 7 we summarize and discuss our results.

## 2 NUMERICAL SIMULATIONS

We performed 2D simulations of the interior of a 5 M$_{\odot}$ star model at zero-age-main-sequence (ZAMS) with the fully compressible time-implicit code MUSIC (Viallet, Baraffe & Walder 2011, 2013; Geroux et al. 2016; Pratt et al. 2016; Viallet et al. 2016; Goffrey et al. 2017; Pratt et al. 2017, 2020). The code solves the inviscid hydrodynamical equations in a fully compressible medium

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

$$ \frac{\partial \rho e}{\partial t} = -\nabla \cdot (\rho e \mathbf{v}) - p \nabla \cdot \mathbf{v} - \nabla \cdot \mathbf{F}_r + Q_{\text{nuc}}, \qquad (2) $$

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

where $\rho$ is the density, $e$ is the specific internal energy, $\mathbf{v}$ is the velocity field, $p$ the gas pressure, $Q_{\text{nuc}}$ is the nuclear energy rate,

<table><caption>**Table 1.** Properties of the initial stellar model used for the 2D hydrodynamical simulations.</caption><thead><tr><th>M/M<sub>☉</sub></th><th>L<sub>star</sub>/L<sub>☉</sub><sup>a</sup></th><th>R<sub>star</sub> (cm)</th><th>r<sub>conv</sub>/R<sub>tot</sub></th><th>H<sub>P, conv</sub> (cm)</th></tr></thead><tbody><tr><td>5</td><td>523</td><td>1.8424 × 10<sup>11</sup></td><td>0.1814</td><td>1.828 × 10<sup>10</sup></td></tr></tbody></table>

Notes. With mass $M$, luminosity $L_{\text{star}}$, radius $R_{\text{star}}$, size of the convective core $r_{\text{conv}}$, and pressure scale height at the convective boundary $H_{\text{P, conv}}$. $^{a}$We use $L_{\odot} = 3.839 \times 10^{33}$ erg s$^{-1}$.

and $g$ is the gravitational acceleration, which we assume is radial. For the stellar simulations considered in this work, the major heat transport that contributes to thermal conductivity is radiative transfer characterized by the radiative flux $\mathbf{F}_r$, given within the diffusion approximation by

$$ \mathbf{F}_r = -\frac{4acT^3}{3\kappa\rho}\nabla T = -\chi\nabla T, \qquad (4) $$

with $\kappa$ is the Rosseland mean opacity of the gas and $\chi$ the radiative conductivity. Realistic stellar opacities and equation of states appropriate for the description of stellar interiors are implemented in MUSIC. Opacities are interpolated from the OPAL tables (Iglesias & Rogers 1996) for solar metallicity, and the equation of state is based on the OPAL EOS tables of Rogers & Nayfonov (2002).

## 2.1 Initial stellar model

The initial structures for the 2D simulations come from a 5 M$_{\odot}$ stellar model computed with the one-dimensional (1D) Lyon stellar evolution code (Baraffe & El Eid 1991; Baraffe et al. 1998), using the same opacities and equation of state as implemented in MUSIC. The initial 1D model has an initial helium abundance in mass fraction $Y = 0.28$ and solar metallicity $Z = 0.02$ and is evolved through the pre-main-sequence (PMS) and the early main-sequence. The initial structure of the simulations have burnt 1 per cent of their hydrogen since the ZAMS. There is no overshooting or diffusion considered during the computation of the 1D model. In the MUSIC simulations, the energy generated by nuclear reactions is taken into account in the energy equation as a source term. The nuclear energy profile is the one from the 1D model, and it is assumed to remain constant during the run of the MUSIC simulations. The properties of the initial 1D model are summarized in Table 1.

## 2.2 Spherical-shell geometry and boundary conditions

A detailed description of the set-up of the simulations can be found in Baraffe et al. (2023). In this section, we only summarize some characteristics relevant to this work. 2D simulations are performed in a spherical shell using spherical coordinates, namely $r$ the radius and the polar angle $\theta$, and assuming azimuthal symmetry in the $\phi$-direction. For all simulations, the inner radius $r_{\text{in}}$ is set at 0.02 $R_{\text{star}}$ and the outer radius $r_{\text{out}}$ at 0.91 $R_{\text{star}}$. Extension of the numerical domain to the photosphere ($r = R_{\text{star}}$) is an open challenge for stellar hydrodynamical simulations given the sharp decrease of the pressure scale height with increasing radius. The domain in the co-latitudinal direction ranges from 0 to $\pi$. We use a uniform grid resolution of $N_r \times N_\theta = 1322 \times 668$ cells. This provides a good resolution of the pressure scale height at the convective boundary $H_{\text{P, conv}}/\Delta r \sim 147$, with $\Delta r = 1.242 \times 10^8$ cm, the size of a radial grid cell.

Concerning the radial boundary conditions, we impose a constant radial derivative on the density on the inner and outer radial boundaries, as discussed in Pratt et al. (2016). For the velocity, we impose reflective conditions at the radial boundaries. For our
<table><caption>**Table 2.** Summary of the 2D simulations.</caption><thead><tr><th>Simulation</th><th>L/L<sub>star</sub></th><th>τ<sub>conv</sub><sup>a</sup> (s)</th><th>N<sub>conv</sub><sup>b</sup></th><th>ω<sub>conv</sub><sup>c</sup> (μHz)</th></tr></thead><tbody><tr><td>ref</td><td>1</td><td>1.5 × 10<sup>6</sup></td><td>52</td><td>0.7</td></tr><tr><td>boost1d1</td><td>10<sup>1</sup></td><td>6.0 × 10<sup>5</sup></td><td>35</td><td>1.7</td></tr><tr><td>boost1d2</td><td>10<sup>2</sup></td><td>3.0 × 10<sup>5</sup></td><td>58</td><td>3.4</td></tr><tr><td>boost1d4</td><td>10<sup>4</sup></td><td>6.5 × 10<sup>4</sup></td><td>54</td><td>15.4</td></tr></tbody></table>

<sup>a</sup>Convective turnover time [see equation (5) for its definition], measured from our simulations.

<sup>b</sup>Number of convective turnover times used for this work.

<sup>c</sup>Convective turnover frequency associated with τ<sub>conv</sub>.

reference model *ref*, the energy flux at the inner and outer radial boundaries are set to the value of the energy flux at that radius in the initial 1D model. In this work, we analyse the impact of enhancing the luminosity and thermal diffusivity by factors of 10, 10<sup>2</sup>, and 10<sup>4</sup>. As in Baraffe et al. (2021) and Le Saux et al. (2022), for the artificially boosted simulations, the energy flux, and equivalently the luminosity, at the boundaries is multiplied by an enhancement factor, and the Rosseland mean opacities κ in MUSIC are decreased by the same factor. At the boundaries in the θ-direction we use reflective boundary conditions for the density, velocity and energy.

The characteristics of the four numerical models used in this work are presented in Table 2. We define the convective turnover time τ<sub>conv</sub> by

$$ \tau_{\text{conv}} := \left\langle \int_{r_{\text{in}}}^{r_{\text{conv}}} \frac{dr}{v_{\text{rms}}(r, t)} \right\rangle_t = \omega_{\text{conv}}^{-1}, \quad (5) $$

where r<sub>in</sub> = 0.02R<sub>star</sub> is the inner boundary of the 2D simulations and r<sub>conv</sub> = 0.1814R<sub>star</sub> is the location of the convective boundary of the stellar model as defined by the Schwarzschild criterion. The root-mean-square of the velocity v<sub>rms</sub> is defined as

$$ v_{\text{rms}}(r, t) := \sqrt{\langle v^2(r, \theta, t) \rangle_S}, \quad (6) $$

where v is the velocity vector. The operators 〈·〉<sub>t</sub> and 〈·〉<sub>S</sub> are time and angular averages. The temporal average is defined as

$$ \langle f(t) \rangle_t := \frac{1}{T} \int_0^T f(t) dt, \quad (7) $$

with T the time of integration considered for the average. The lower limit of the integral t = 0 corresponds to the time from which convection is in steady state (see Baraffe et al. 2023, for details on the definition of this steady state.). The angular average is defined as

$$ \langle g(\theta) \rangle_S := \frac{1}{4\pi} \int_S g(\theta) 2\pi \sin\theta d\theta. \quad (8) $$

In equation (5), we also define the convective turnover frequency ω<sub>conv</sub> which is the frequency associated with the characteristic time-scale τ<sub>conv</sub>. Aerts et al. (2021) provide observed values of the convective turnover frequency for Slowly Pulsating B (SPB) stars with mass approximately between 3M☉ and 9M☉. They found typical values of ω<sub>conv</sub> in the range 0.2 to 0.5 μHz, which is close to the value of 0.7 μHz in our reference model. We consider the stellar model studied in this paper to be a template of intermediate-mass stars with a convective core, with mass between 3 and 20 M☉. Therefore, our conclusions regarding wave damping and propagation could be applied, at least qualitatively, to other stellar masses.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Radial profile of the Brunt–Väisälä frequency (left y-axis, blue curve) and the helium mass fraction Y (right y-axis, brown dashed curve) in the initial 1D model. The hatched regions are not considered in the 2D simulations.</figcaption></figure>

## 2.3 Stratification and radiative diffusivity

A medium is stably stratified when it is stable against convection. A particle of fluid moved radially in a stratified medium will oscillate around its equilibrium position with a finite amplitude (see e.g. Lighthill 1978). Therefore, one of the main characteristic of such medium is its ability to support internal waves, for which the restoring force is buoyancy. The maximum frequency of oscillation of this particle is set by the Brunt–Väisälä, or buoyancy, frequency N. Expressed in Hertz, it is defined as

$$ N = \frac{1}{2\pi} \sqrt{g \left( \frac{1}{\Gamma_1} \frac{d \ln p}{dr} - \frac{d \ln \rho}{dr} \right)}, \quad (9) $$

where Γ<sub>1</sub> is the first adiabatic exponent,

$$ \Gamma_1 = \left( \frac{d \ln \rho}{d \ln p} \right)_{\text{ad}}. \quad (10) $$

A fluid is stably stratified if N<sup>2</sup> > 0 and the convective instability arises if N<sup>2</sup> ≤ 0.

Fig. 1 shows the radial profile of the Brunt–Väisälä frequency (left y-axis, blue curve) and the helium mass fraction Y (right y-axis, brown dashed curve) in the initial 1D model. The peak of the Brunt–Väisälä frequency just above the convective core is a result of the small gradient of Y in this region. In our simulations, the hatched regions are excluded. Internal gravity waves are able to propagate in radiative zones and are evanescent in convective zones. This means they can travel in region where N<sup>2</sup> > 0, with the condition on their frequency ω < N. From Fig. 1, we can see that an IGW of frequency ω might not propagate in the whole radiative envelope. For instance, a wave with a frequency of 150 μHz may only propagate between r ~ 0.25R<sub>star</sub> and r ~ 0.6R<sub>star</sub>.

The radiative diffusivity is of major importance as it is the main mechanism that damp waves. It is defined as

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

with χ the radiative conductivity defined by equation (4) and c<sub>p</sub> the specific heat capacity at constant pressure. In our reference simulation, model *ref*, we use a realistic profile of radiative diffusivity for a 5M☉ model. This allows modelling of a realistic damping of IGWs. This profile is displayed in Fig. 2 (left y-axis). This figure also shows the radial profile of the density in our model (right y-axis). Note that the profiles displayed in Fig. 2 are the ones used in our simulation *ref*. This is important as these two quantities define the propagation properties of IGWs.
<table><caption>**Table 2.** Summary of the 2D simulations.</caption><thead><tr><th>Simulation</th><th>L/L<sub>star</sub></th><th>τ<sub>conv</sub><sup>a</sup> (s)</th><th>N<sub>conv</sub><sup>b</sup></th><th>ω<sub>conv</sub><sup>c</sup> (μHz)</th></tr></thead><tbody><tr><td>ref</td><td>1</td><td>1.5 × 10<sup>6</sup></td><td>52</td><td>0.7</td></tr><tr><td>boost1d1</td><td>10<sup>1</sup></td><td>6.0 × 10<sup>5</sup></td><td>35</td><td>1.7</td></tr><tr><td>boost1d2</td><td>10<sup>2</sup></td><td>3.0 × 10<sup>5</sup></td><td>58</td><td>3.4</td></tr><tr><td>boost1d4</td><td>10<sup>4</sup></td><td>6.5 × 10<sup>4</sup></td><td>54</td><td>15.4</td></tr></tbody></table>

<sup>a</sup>Convective turnover time [see equation (5) for its definition], measured from our simulations.

<sup>b</sup>Number of convective turnover times used for this work.

<sup>c</sup>Convective turnover frequency associated with τ<sub>conv</sub>.

reference model ref, the energy flux at the inner and outer radial boundaries are set to the value of the energy flux at that radius in the initial 1D model. In this work, we analyse the impact of enhancing the luminosity and thermal diffusivity by factors of 10, 10<sup>2</sup>, and 10<sup>4</sup>. As in Baraffe et al. (2021) and Le Saux et al. (2022), for the artificially boosted simulations, the energy flux, and equivalently the luminosity, at the boundaries is multiplied by an enhancement factor, and the Rosseland mean opacities κ in MUSIC are decreased by the same factor. At the boundaries in the θ-direction we use reflective boundary conditions for the density, velocity and energy.

The characteristics of the four numerical models used in this work are presented in Table 2. We define the convective turnover time τ<sub>conv</sub> by

$$ \tau_{\text{conv}} := \left\langle \int_{r_{\text{in}}}^{r_{\text{conv}}} \frac{dr}{v_{\text{rms}}(r, t)} \right\rangle_t = \omega_{\text{conv}}^{-1}, \quad (5) $$

where r<sub>in</sub> = 0.02R<sub>star</sub> is the inner boundary of the 2D simulations and r<sub>conv</sub> = 0.1814R<sub>star</sub> is the location of the convective boundary of the stellar model as defined by the Schwarzschild criterion. The root-mean-square of the velocity v<sub>rms</sub> is defined as

$$ v_{\text{rms}}(r, t) := \sqrt{\langle v^2(r, \theta, t) \rangle_S}, \quad (6) $$

where v is the velocity vector. The operators 〈·〉<sub>t</sub> and 〈·〉<sub>S</sub> are time and angular averages. The temporal average is defined as

$$ \langle f(t) \rangle_t := \frac{1}{T} \int_0^T f(t) dt, \quad (7) $$

with T the time of integration considered for the average. The lower limit of the integral t = 0 corresponds to the time from which convection is in steady state (see Baraffe et al. 2023, for details on the definition of this steady state.). The angular average is defined as

$$ \langle g(\theta) \rangle_S := \frac{1}{4\pi} \int_S g(\theta) 2\pi \sin\theta d\theta. \quad (8) $$

In equation (5), we also define the convective turnover frequency ω<sub>conv</sub> which is the frequency associated with the characteristic time-scale τ<sub>conv</sub>. Aerts et al. (2021) provide observed values of the convective turnover frequency for Slowly Pulsating B (SPB) stars with mass approximately between 3M☉ and 9M☉. They found typical values of ω<sub>conv</sub> in the range 0.2 to 0.5 μHz, which is close to the value of 0.7 μHz in our reference model. We consider the stellar model studied in this paper to be a template of intermediate-mass stars with a convective core, with mass between 3 and 20 M☉. Therefore, our conclusions regarding wave damping and propagation could be applied, at least qualitatively, to other stellar masses.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Radial profile of the Brunt–Väisälä frequency (left y-axis, blue curve) and the helium mass fraction Y (right y-axis, brown dashed curve) in the initial 1D model. The hatched regions are not considered in the 2D simulations.</figcaption></figure>

## 2.3 Stratification and radiative diffusivity

A medium is stably stratified when it is stable against convection. A particle of fluid moved radially in a stratified medium will oscillate around its equilibrium position with a finite amplitude (see e.g. Lighthill 1978). Therefore, one of the main characteristic of such medium is its ability to support internal waves, for which the restoring force is buoyancy. The maximum frequency of oscillation of this particle is set by the Brunt–Väisälä, or buoyancy, frequency N. Expressed in Hertz, it is defined as

$$ N = \frac{1}{2\pi} \sqrt{g \left( \frac{1}{\Gamma_1} \frac{d \ln p}{dr} - \frac{d \ln \rho}{dr} \right)}, \quad (9) $$

where Γ<sub>1</sub> is the first adiabatic exponent,

$$ \Gamma_1 = \left( \frac{d \ln \rho}{d \ln p} \right)_{\text{ad}}. \quad (10) $$

A fluid is stably stratified if N<sup>2</sup> > 0 and the convective instability arises if N<sup>2</sup> ≤ 0.

Fig.1 shows the radial profile of the Brunt–Väisälä frequency (left y-axis, blue curve) and the helium mass fraction Y (right y-axis, brown dashed curve) in the initial 1D model. The peak of the Brunt–Väisälä frequency just above the convective core is a result of the small gradient of Y in this region. In our simulations, the hatched regions are excluded. Internal gravity waves are able to propagate in radiative zones and are evanescent in convective zones. This means they can travel in region where N<sup>2</sup> > 0, with the condition on their frequency ω < N. From Fig. 1, we can see that an IGW of frequency ω might not propagate in the whole radiative envelope. For instance, a wave with a frequency of 150 μHz may only propagate between r ~ 0.25R<sub>star</sub> and r ~ 0.6R<sub>star</sub>.

The radiative diffusivity is of major importance as it is the main mechanism that damp waves. It is defined as

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

with χ the radiative conductivity defined by equation (4) and c<sub>p</sub> the specific heat capacity at constant pressure. In our reference simulation, model ref, we use a realistic profile of radiative diffusivity for a 5M☉ model. This allows modelling of a realistic damping of IGWs. This profile is displayed in Fig. 2 (left y-axis). This figure also shows the radial profile of the density in our model (right y-axis). Note that the profiles displayed in Fig. 2 are the ones used in our simulation ref. This is important as these two quantities define the propagation properties of IGWs.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Radial profile of radiative diffusivity $\kappa_{\text{rad}}$ (red, left y-axis) and density $\rho$ (green right x-axis) for the initial 1D model. The blue vertical dotted line indicates the convective boundary. The hatched region is not considered in the 2D simulations.</figcaption></figure>

Indeed, as they propagate towards the surface, IGWs amplitudes grow due to decreasing density and decay due to radiative damping. The variation of both the radiative diffusivity and the density is very large in a star, by 9 orders of magnitude between the centre and the surface. This is one of the main challenges that hydrodynamical simulations have to deal with. In our simulations, that exclude the 9 per cent outer layers, this variation is restricted to approximately 5 orders of magnitude. These variations have to be taken into account in hydrodynamical simulations because they significantly impact propagation of IGWs as we will see in the next sections.

# 3 VELOCITIES

Internal gravity waves manifest themselves as perturbations in density, temperature, luminosity or velocity. In this section, we focus on the radial velocity.

## 3.1 Radial velocity pattern

Fig. 3 shows snapshots of the radial velocity $v_r$ for the four simulations. For better visualization, the radial velocity is normalized by the root-mean-square value of the radial velocity $v_{r,rms}$. The convective core extends from the centre up to $r_{\text{conv}} = 0.1814 R_{\text{star}}$. In this region, the structure of the flows is similar in the four snapshots, with large coherent upflows (red) and downflows (blue).

The three simulations *ref*, *boost1d1* and *boost1d2* present similar patterns in the radiative envelope with the classical spiral structure corresponding to wavefronts of IGWs. The inclination of these wavefronts with respect to the convective boundary, i.e. wrt the horizontal (or angular) direction, defines an angle $\alpha$, to which correspond a particular frequency. This is known as the St. Andrews cross (Sutherland 2010) and results from the dispersion relation for IGWs that set the direction of propagation of the waves (see Eq. (7.61) from Vallis 2017)

$$ \frac{\omega^2}{N^2} = \frac{k_h^2}{k^2} := \cos^2 \alpha, \qquad (12) $$

where $\omega$ is the wave frequency and $k$ the total wavenumber defined by $k = \sqrt{k_r^2 + k_h^2}$, with $k_r$ its radial part and $k_h$ its horizontal part defined at a given radius $r$ by:

$$ k_h^2 := \frac{\ell(\ell + 1)}{r^2}, \qquad (13) $$

with $\ell \ge 0$ the spherical harmonic, or angular, degree. Fig. 3 suggests that the angle $\alpha$ increases with the boosting factor. Therefore, waves of higher frequencies dominate in the radiative zone when the luminosity is increased, an effect already suggested by Stein (1967) in their analytical work. This was also observed in the solar-like star simulations from Le Saux et al. (2022). The structure of the radiative envelope of *boost1d4* does not present the classical spiral IGWs wavefronts pattern. The structure in this most boosted case is larger scale and not periodic any more. Similar patterns are observed in other simulations of intermediate-mass stars using the same or different artificial boosting factors for the luminosity and the radiative diffusivity (see e.g. Edelmann et al. 2019; Horst et al. 2020).

## 3.2 Radial evolution of the velocity amplitude

According to linear theory, the evolution of the amplitude of an IGW is expected to depend on the stratification of the supporting medium and on spatial damping due to radiative effects as well as on the frequency and wavelength of the wave. The analytical formula that expresses these dependences is (Press 1981; Zahn, Talon & Matias 1997)

$$ v_r(r, \ell, \omega) = C \rho^{-1/2} k_h^{3/2} \left( \frac{N^2 - \omega^2}{\omega^2} \right)^{-1/4} e^{-\tau/2}, \qquad (14) $$

where $\omega$ and $N$ are expressed in hertz and $k_h$ is the local horizontal wavenumber, already defined in equation (13). The exponential term in equation (14) represents the effect of radiative damping, and the parameter $\tau$ is defined as

$$ \tau(r, \ell, \omega) = [\ell(\ell + 1)]^{3/2} \int_{r_e}^{r} \kappa_{\text{rad}} \frac{N^3}{\omega^4} \left( \frac{N^2}{N^2 - \omega^2} \right)^{1/2} \frac{dr}{r^3}, \qquad (15) $$

where $r_e$ is the radius at which the waves are excited. It is important to keep in mind that in boosted simulations, the radiative diffusivity is enhanced by the same amount as the luminosity. Consequently, wave damping by radiative diffusion is enhanced in a boosted simulation. In equation (14), $C$ is a constant fixing the amplitude. In this work, we chose to fix it such as the analytical velocity amplitude of a wave at $r = r_e$ matches the amplitude of the velocity in the simulations, i.e. $v_0(\ell, \omega) = v_r(r_e, \ell, \omega)$. Therefore, we obtain

$$ C = v_0(\ell, \omega) \rho_0^{1/2} k_{h,0}^{-3/2} \left( \frac{N_0^2 - \omega^2}{\omega^2} \right)^{1/4}, \qquad (16) $$

with $k_{h,0} = \sqrt{\ell(\ell + 1)}/r_e$, $\rho_0 = \rho(r_e)$, $N_0 = N(r_e)$ and by definition $\tau(r_e, \ell, \omega) = 0$. Finally, we can write the analytical expression of the radial velocity as

$$ v_r(r, \ell, \omega) = v_0(\ell, \omega) \left( \frac{\rho}{\rho_0} \right)^{-1/2} \left( \frac{k_h}{k_{h,0}} \right)^{3/2} \left( \frac{N^2 - \omega^2}{N_0^2 - \omega^2} \right)^{-1/4} e^{-\tau/2}. \qquad (17) $$

This formula is similar to the one obtained by Ratnasingam, Edelmann & Rogers (2019). In our simulations, we determine the amplitude of the radial velocity at a given frequency $\omega$ and angular degree $\ell$ using a temporal Fourier transform and a decomposition on the spherical harmonic basis of the velocities computed by MUSIC. The definition we are using for spherical harmonics and Fourier transform are the same as in Le Saux et al. (2022) (see their appendix A). We obtain the power spectrum of the radial velocity $P[\hat{v}_r^2](r, \ell, \omega)$ which scales as the amplitude squared of a given mode $(\ell, \omega)$.

Fig. 4 compares the analytical expression from equation (17) (dotted lines) to the corresponding wave velocity amplitude from MUSIC
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Radial profile of radiative diffusivity $\kappa_{\text{rad}}$ (red, left y-axis) and density $\rho$ (green right x-axis) for the initial 1D model. The blue vertical dotted line indicates the convective boundary. The hatched region is not considered in the 2D simulations.</figcaption></figure>

Indeed, as they propagate towards the surface, IGWs amplitudes grow due to decreasing density and decay due to radiative damping. The variation of both the radiative diffusivity and the density is very large in a star, by 9 orders of magnitude between the centre and the surface. This is one of the main challenges that hydrodynamical simulations have to deal with. In our simulations, that exclude the 9 per cent outer layers, this variation is restricted to approximately 5 orders of magnitude. These variations have to be taken into account in hydrodynamical simulations because they significantly impact propagation of IGWs as we will see in the next sections.

# 3 VELOCITIES

Internal gravity waves manifest themselves as perturbations in density, temperature, luminosity or velocity. In this section, we focus on the radial velocity.

## 3.1 Radial velocity pattern

Fig. 3 shows snapshots of the radial velocity $v_r$ for the four simulations. For better visualization, the radial velocity is normalized by the root-mean-square value of the radial velocity $v_{r,rms}$. The convective core extends from the centre up to $r_{\text{conv}} = 0.1814 R_{\text{star}}$. In this region, the structure of the flows is similar in the four snapshots, with large coherent upflows (red) and downflows (blue).

The three simulations *ref*, *boost1d1* and *boost1d2* present similar patterns in the radiative envelope with the classical spiral structure corresponding to wavefronts of IGWs. The inclination of these wavefronts with respect to the convective boundary, i.e. wrt the horizontal (or angular) direction, defines an angle $\alpha$, to which correspond a particular frequency. This is known as the St. Andrews cross (Sutherland 2010) and results from the dispersion relation for IGWs that set the direction of propagation of the waves (see Eq. (7.61) from Vallis 2017)

$$ \frac{\omega^2}{N^2} = \frac{k_h^2}{k^2} := \cos^2 \alpha, \qquad (12) $$

where $\omega$ is the wave frequency and $k$ the total wavenumber defined by $k = \sqrt{k_r^2 + k_h^2}$, with $k_r$ its radial part and $k_h$ its horizontal part defined at a given radius $r$ by:

$$ k_h^2 := \frac{\ell(\ell + 1)}{r^2}, \qquad (13) $$

with $\ell \ge 0$ the spherical harmonic, or angular, degree. Fig. 3 suggests that the angle $\alpha$ increases with the boosting factor. Therefore, waves of higher frequencies dominate in the radiative zone when the luminosity is increased, an effect already suggested by Stein (1967) in their analytical work. This was also observed in the solar-like star simulations from Le Saux et al. (2022). The structure of the radiative envelope of *boost1d4* does not present the classical spiral IGWs wavefronts pattern. The structure in this most boosted case is larger scale and not periodic any more. Similar patterns are observed in other simulations of intermediate-mass stars using the same or different artificial boosting factors for the luminosity and the radiative diffusivity (see e.g. Edelmann et al. 2019; Horst et al. 2020).

## 3.2 Radial evolution of the velocity amplitude

According to linear theory, the evolution of the amplitude of an IGW is expected to depend on the stratification of the supporting medium and on spatial damping due to radiative effects as well as on the frequency and wavelength of the wave. The analytical formula that expresses these dependences is (Press 1981; Zahn, Talon & Matias 1997)

$$ v_r(r, \ell, \omega) = C \rho^{-1/2} k_h^{3/2} \left( \frac{N^2 - \omega^2}{\omega^2} \right)^{-1/4} e^{-\tau/2}, \qquad (14) $$

where $\omega$ and $N$ are expressed in hertz and $k_h$ is the local horizontal wavenumber, already defined in equation (13). The exponential term in equation (14) represents the effect of radiative damping, and the parameter $\tau$ is defined as

$$ \tau(r, \ell, \omega) = [\ell(\ell + 1)]^{3/2} \int_{r_e}^{r} \kappa_{\text{rad}} \frac{N^3}{\omega^4} \left( \frac{N^2}{N^2 - \omega^2} \right)^{1/2} \frac{dr}{r^3}, \qquad (15) $$

where $r_e$ is the radius at which the waves are excited. It is important to keep in mind that in boosted simulations, the radiative diffusivity is enhanced by the same amount as the luminosity. Consequently, wave damping by radiative diffusion is enhanced in a boosted simulation. In equation (14), $C$ is a constant fixing the amplitude. In this work, we chose to fix it such as the analytical velocity amplitude of a wave at $r = r_e$ matches the amplitude of the velocity in the simulations, i.e. $v_0(\ell, \omega) = v_r(r_e, \ell, \omega)$. Therefore, we obtain

$$ C = v_0(\ell, \omega) \rho_0^{1/2} k_{h,0}^{-3/2} \left( \frac{N_0^2 - \omega^2}{\omega^2} \right)^{1/4}, \qquad (16) $$

with $k_{h,0} = \sqrt{\ell(\ell + 1)}/r_e$, $\rho_0 = \rho(r_e)$, $N_0 = N(r_e)$ and by definition $\tau(r_e, \ell, \omega) = 0$. Finally, we can write the analytical expression of the radial velocity as

$$ v_r(r, \ell, \omega) = v_0(\ell, \omega) \left( \frac{\rho}{\rho_0} \right)^{-1/2} \left( \frac{k_h}{k_{h,0}} \right)^{3/2} \left( \frac{N^2 - \omega^2}{N_0^2 - \omega^2} \right)^{-1/4} e^{-\tau/2}. \qquad (17) $$

This formula is similar to the one obtained by Ratnasingam, Edelmann & Rogers (2019). In our simulations, we determine the amplitude of the radial velocity at a given frequency $\omega$ and angular degree $\ell$ using a temporal Fourier transform and a decomposition on the spherical harmonic basis of the velocities computed by MUSIC. The definition we are using for spherical harmonics and Fourier transform are the same as in Le Saux et al. (2022) (see their appendix A). We obtain the power spectrum of the radial velocity $P[\hat{v}_r^2](r, \ell, \omega)$ which scales as the amplitude squared of a given mode $(\ell, \omega)$.

Fig. 4 compares the analytical expression from equation (17) (dotted lines) to the corresponding wave velocity amplitude from MUSIC
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Visualization of the radial velocity for the four stellar simulations *ref*, *boost1d1*, *boost1d2*, and *boost1d4* as a function of radius, *r*, and co-latitude, *θ*. The radial velocity is normalized by the rms radial velocity. Positive values of the radial velocity (ref) are outward and negative (blue) are inward.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Wave amplitude as a function of normalized radius for the four simulations *ref* (orange), *boost1d1* (indigo), *boost1d2* (yellow), and *boost1d4* (green) for angular degree $\ell = 5$ and frequencies $\omega = 45.2$ µHz. The vertical black lines indicate the boundaries of the simulation domain and the vertical blue line indicate the convective boundary. The solid lines are the velocity measured in the simulations. Dotted and dashed lines are the theoretical velocity amplitudes computed with equation (17) with and without the damping term, respectively.</figcaption></figure>

simulations, $\sqrt{P[\hat{v}_r^2]}$ (solid lines), as a function of normalized radius for the four numerical models. We have also included the analytical velocity amplitude without the damping term (dashed lines), i.e. we are neglecting the term $e^{-\tau/2}$ in equation (17). The spatial boundaries of the simulation domain at $r_{in} = 0.02R_{star}$ and $r_{out} = 0.91R_{star}$ are specified by the vertical black dashed lines, and the convective boundary $r_{conv} = 0.1814R_{star}$ is indicated by the vertical blue dashed line. For the analytical expressions, we need to set the initial velocity amplitude of the waves in the excitation region, i.e. close to the convective boundary. We arbitrarily chose for $v_0$, the value of the radial velocity in the simulations at $r_e = 0.183R_{star}$, just above the convective core. Fig. 4 presents the velocity amplitude of a wave with angular degree $\ell = 5$ and frequency $\omega = 45.2$ µHz. In the four simulations, the wave amplitude shows a similar oscillatory pattern, where the troughs are the radial nodes of the corresponding standing wave. To form these standing waves, propagating IGWs generated by

core convection travel towards the surface of the star until they reflect at their outer turning point or at the top of the simulation domain. A turning point is defined as the radius where $k_r^2 = 0$, or equivalently $\omega = N$ from equation (12). Travelling back towards the centre they again reflect at their inner turning point. Travelling back and forth, propagating waves of a given frequency and angular degree interfere with themselves and form standing waves. These are called gravity modes, or g modes, and have high amplitudes that can be seen as high narrow peaks in the power spectra of the velocity. For the mode observed in Fig. 4, the oscillation code GYRE (Townsend & Teitler 2013; Townsend, Goldstein & Zweibel 2018; Goldstein & Townsend 2020) predicts a radial order<fnref n="1" /> $n = -9$, which is also the number of nodes present in the simulations. This confirms that we see g modes in Fig. 4.

Compared with the analytical predictions from equation (17) (dotted lines in Fig. 4), the simulation velocities present a very similar global evolution from the convective boundary up to the top of the simulation domain. Except for the oscillations, but these are not taken into account in the linear analytical expression of the velocity amplitude. However, there are notable differences between the four numerical models. First, as expected the velocity amplitude of the waves increases with the enhancement of the luminosity. Secondly, the analytical amplitudes predict a sharp drop at a given radius close to the surface, which corresponds to the location where the wave is totally damped out, i.e. the location where the wave will deposit most of its energy. This abrupt drop is indeed due to radiative damping as it is not present in the case with no damping (dashed lines). For models *ref*, *boost1d1*, and *boost1d2*, this abrupt drop is approximately at $r \simeq R_{star}$, but for model *boost1d4* it is located around $r \simeq 0.91R_{star}$. By looking at equation (17) we can see the radiative diffusivity, which is enhanced by the same amount as the luminosity, is included in the expression of $\tau$. As a result, damping of waves by radiative diffusion increases in a boosted model. This is why the sharp drop of the dotted lines is not located at the same radius for the four cases. However, for the waves considered, with $\ell = 5$ and $\omega = 45.2$ µHz, this sharp drop is located at radii $r \ge 0.91R_{star}$, thus they are able to reach the top of the simulation domain.

<footnote n="1">¹By convention, the radial order *n* is negative for a g mode and positive for a p mode (standing acoustic waves).</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Visualization of the radial velocity for the four stellar simulations *ref*, *boost1d1*, *boost1d2*, and *boost1d4* as a function of radius, *r*, and co-latitude, *θ*. The radial velocity is normalized by the rms radial velocity. Positive values of the radial velocity (ref) are outward and negative (blue) are inward.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Wave amplitude as a function of normalized radius for the four simulations *ref* (orange), *boost1d1* (indigo), *boost1d2* (yellow), and *boost1d4* (green) for angular degree $ℓ = 5$ and frequencies $ω = 45.2$ µHz. The vertical black lines indicate the boundaries of the simulation domain and the vertical blue line indicate the convective boundary. The solid lines are the velocity measured in the simulations. Dotted and dashed lines are the theoretical velocity amplitudes computed with equation (17) with and without the damping term, respectively.</figcaption></figure>

simulations, $\sqrt{P[\hat{v}_r^2]}$ (solid lines), as a function of normalized radius for the four numerical models. We have also included the analytical velocity amplitude without the damping term (dashed lines), i.e. we are neglecting the term $e^{-\tau/2}$ in equation (17). The spatial boundaries of the simulation domain at $r_{in} = 0.02R_{star}$ and $r_{out} = 0.91R_{star}$ are specified by the vertical black dashed lines, and the convective boundary $r_{conv} = 0.1814R_{star}$ is indicated by the vertical blue dashed line. For the analytical expressions, we need to set the initial velocity amplitude of the waves in the excitation region, i.e. close to the convective boundary. We arbitrarily chose for $v_0$, the value of the radial velocity in the simulations at $r_e = 0.183R_{star}$, just above the convective core. Fig. 4 presents the velocity amplitude of a wave with angular degree $ℓ = 5$ and frequency $ω = 45.2$ µHz. In the four simulations, the wave amplitude shows a similar oscillatory pattern, where the troughs are the radial nodes of the corresponding standing wave. To form these standing waves, propagating IGWs generated by

core convection travel towards the surface of the star until they reflect at their outer turning point or at the top of the simulation domain. A turning point is defined as the radius where $k_r^2 = 0$, or equivalently $ω = N$ from equation (12). Travelling back towards the centre they again reflect at their inner turning point. Travelling back and forth, propagating waves of a given frequency and angular degree interfere with themselves and form standing waves. These are called gravity modes, or g modes, and have high amplitudes that can be seen as high narrow peaks in the power spectra of the velocity. For the mode observed in Fig. 4, the oscillation code GYRE (Townsend & Teitler 2013; Townsend, Goldstein & Zweibel 2018; Goldstein & Townsend 2020) predicts a radial order<fnref n="1" /> $n = -9$, which is also the number of nodes present in the simulations. This confirms that we see g modes in Fig. 4.

Compared with the analytical predictions from equation (17) (dotted lines in Fig. 4), the simulation velocities present a very similar global evolution from the convective boundary up to the top of the simulation domain. Except for the oscillations, but these are not taken into account in the linear analytical expression of the velocity amplitude. However, there are notable differences between the four numerical models. First, as expected the velocity amplitude of the waves increases with the enhancement of the luminosity. Secondly, the analytical amplitudes predict a sharp drop at a given radius close to the surface, which corresponds to the location where the wave is totally damped out, i.e. the location where the wave will deposit most of its energy. This abrupt drop is indeed due to radiative damping as it is not present in the case with no damping (dashed lines). For models *ref*, *boost1d1*, and *boost1d2*, this abrupt drop is approximately at $r \simeq R_{star}$, but for model *boost1d4* it is located around $r \simeq 0.91R_{star}$. By looking at equation (17) we can see the radiative diffusivity, which is enhanced by the same amount as the luminosity, is included in the expression of $τ$. As a result, damping of waves by radiative diffusion increases in a boosted model. This is why the sharp drop of the dotted lines is not located at the same radius for the four cases. However, for the waves considered, with $ℓ = 5$ and $ω = 45.2$ µHz, this sharp drop is located at radii $r \ge 0.91R_{star}$, thus they are able to reach the top of the simulation domain.

<footnote n="1">¹By convention, the radial order *n* is negative for a g mode and positive for a p mode (standing acoustic waves).</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Same as Fig. 4 but for waves with $ℓ = 10$ and $ω = 31.0$ µHz (left-hand panel) and with $ℓ = 10$ and $ω = 8.3$ µHz (right-hand panel).</figcaption></figure>

### 3.3 Influence of the boost

From equations (15) and (17), it is clear that when the angular degree $ℓ$ and the frequency $ω$ are changed, the amplitude and the damping of the corresponding wave are impacted. Indeed, equation (17) shows that wave amplitude depends on $ℓ$ and $ω$ and equation (15) shows that waves with higher $ℓ$ and smaller $ω$ will be damped more efficiently. The radial velocity amplitude for a wave with $ℓ = 10$ and $ω = 31.0$ µHz is plotted on the left-hand panel of Fig. 5. As expected, the damping is more important than for the wave with $ℓ = 5$ and $ω = 45.2$ µHz (see Fig. 4). This is highlighted by the abrupt drop of the theoretical amplitude (dotted lines) which is shifted towards smaller radii in the three boosted models. Indeed, for this wave, the drop is located at $r \approx 0.95R_{\text{star}}$ in model *boost1d1*, $r \approx 0.91R_{\text{star}}$ in model *boost1d2* and $r \approx 0.8R_{\text{star}}$ in model *boost1d4*. However, for model *ref* the location of the drop is not changed, meaning that the effect of radiative damping remains weak for this wave with $ℓ = 10$ and $ω = 31.0$ µHz. Note that for the four simulations the agreement with theory is still relatively good. The case with no damping (dashed lines) is now clearly different from the other two (solid and dotted lines). In that case, the shape of the curves remains the same as in Fig. 4. Neglecting the damping for these waves would imply that they would be able to propagate up to the surface for all models, yielding an erroneous prediction.

On the right-hand panel of Fig. 5, the angular degree and frequency are set to $ℓ = 10$ and $ω = 8.3$ µHz respectively, implying more efficient damping. Model *boost1d4* is not plotted here as we do not expect such low frequency IGWs to be excited in this simulation. Indeed, the frequency $ω = 8.3$ µHz is smaller than the convective frequency, $ω_{\text{conv}} = 15.4$ µHz, for this simulation (see Table 2). A convective region with associated frequency $ω_{\text{conv}}$ is expected to generate waves with frequencies $ω \ge ω_{\text{conv}}$ (Lecoanet & Quataert 2013). On this plot, the abrupt drop of the analytically predicted velocity amplitude (dotted lines) for models *boost1d1* and *boost1d2* is located at $r \approx 0.8R_{\text{star}}$ and $r \approx 0.85R_{\text{star}}$ respectively. We thus observe a similar phenomenon as for model *boost1d4* in the left-hand panel of Fig. 5. Namely, that the waves in the simulations are damped before the top of the numerical domain and do not form g modes.

Our results highlight that the artificial enhancement of the luminosity and the radiative diffusivity of a numerical model impacts not only the amplitude of the waves but also their spatial damping. This is particularly important at low frequencies, as expected from equation

(15). This enhanced damping of waves in the low frequency regime in boosted simulations was already suggested by Horst et al. (2020). When the luminosity is boosted, the increased damping implies that waves over a smaller range of frequencies reach the top of the cavity. Particularly, low frequency g modes are fully damped in boosted simulations compared to the one with realistic luminosity. In addition, waves of given frequency and angular degree will be damped out in different locations in boosted simulations, therefore depositing their energy in different regions. In a case where the luminosity is artificially enhanced but not the radiative diffusivity, as in the simulation of Horst et al. (2020), the location where IGWs are damped out is not modified. However, because of the higher luminosity, the dominant frequency range of excited waves will be different from in a non-boosted case (see Le Saux et al. 2022, and Section 4). Therefore, in simulations with different enhancement factors for the luminosity and radiative diffusivity, we do not expect the propagation of IGWs to be identical as in a model with realistic luminosity.

Finally, comparison between the two analytical cases with and without damping also highlights that running simulations with unrealistic radiative diffusivity may not capture the proper propagation properties of IGWs. Particularly for low frequency waves, which are the ones more impacted by radiative diffusion. In these simulations, transport by IGWs should be studied with caution as well as their ability to propagate up to the surface or not.

### 3.4 Radial kinetic energy density

Fig. 6 presents the power spectrum of the radial velocity $P[\hat{v}_r^2]$, as a function of normalized radius and frequency for an angular degree $ℓ = 20$. This quantity $P[\hat{v}_r^2]$ is equivalent to the radial kinetic energy density. This plot offers a general overview on the dependence of radiative damping of IGWs on frequency and radius. In the convective core, between the bottom of the plot and $r_{\text{conv}}$, the spectrum is relatively homogeneous at all frequencies and characteristic of a convective zone. The bright ridges observed in the radiative zone, at $r \ge 0.1814R_{\text{star}}$, are high-amplitude g modes. The dark knots observed in the bright ridges are the radial nodes of the considered mode. The number of nodes for a given mode defines the radial order $n$ of the mode, and it increases as the frequency decreases. This is an important characteristic of g modes (Aerts, Christensen-Dalsgaard & Kurtz 2010).
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Same as Fig. 4 but for waves with $ℓ = 10$ and $ω = 31.0$ µHz (left-hand panel) and with $ℓ = 10$ and $ω = 8.3$ µHz (right-hand panel).</figcaption></figure>

### 3.3 Influence of the boost

From equations (15) and (17), it is clear that when the angular degree $ℓ$ and the frequency $ω$ are changed, the amplitude and the damping of the corresponding wave are impacted. Indeed, equation (17) shows that wave amplitude depends on $ℓ$ and $ω$ and equation (15) shows that waves with higher $ℓ$ and smaller $ω$ will be damped more efficiently. The radial velocity amplitude for a wave with $ℓ = 10$ and $ω = 31.0$ µHz is plotted on the left-hand panel of Fig. 5. As expected, the damping is more important than for the wave with $ℓ = 5$ and $ω = 45.2$ µHz (see Fig. 4). This is highlighted by the abrupt drop of the theoretical amplitude (dotted lines) which is shifted towards smaller radii in the three boosted models. Indeed, for this wave, the drop is located at $r \simeq 0.95R_{\text{star}}$ in model *boost1d1*, $r \simeq 0.91R_{\text{star}}$ in model *boost1d2* and $r \simeq 0.8R_{\text{star}}$ in model *boost1d4*. However, for model *ref* the location of the drop is not changed, meaning that the effect of radiative damping remains weak for this wave with $ℓ = 10$ and $ω = 31.0$ µHz. Note that for the four simulations the agreement with theory is still relatively good. The case with no damping (dashed lines) is now clearly different from the other two (solid and dotted lines). In that case, the shape of the curves remains the same as in Fig. 4. Neglecting the damping for these waves would imply that they would be able to propagate up to the surface for all models, yielding an erroneous prediction.

On the right-hand panel of Fig. 5, the angular degree and frequency are set to $ℓ = 10$ and $ω = 8.3$ µHz respectively, implying more efficient damping. Model *boost1d4* is not plotted here as we do not expect such low frequency IGWs to be excited in this simulation. Indeed, the frequency $ω = 8.3$ µHz is smaller than the convective frequency, $ω_{\text{conv}} = 15.4$ µHz, for this simulation (see Table 2). A convective region with associated frequency $ω_{\text{conv}}$ is expected to generate waves with frequencies $ω \ge ω_{\text{conv}}$ (Lecoanet & Quataert 2013). On this plot, the abrupt drop of the analytically predicted velocity amplitude (dotted lines) for models *boost1d1* and *boost1d2* is located at $r \simeq 0.8R_{\text{star}}$ and $r \simeq 0.85R_{\text{star}}$ respectively. We thus observe a similar phenomenon as for model *boost1d4* in the left-hand panel of Fig. 5. Namely, that the waves in the simulations are damped before the top of the numerical domain and do not form g modes.

Our results highlight that the artificial enhancement of the luminosity and the radiative diffusivity of a numerical model impacts not only the amplitude of the waves but also their spatial damping. This is particularly important at low frequencies, as expected from equation

(15). This enhanced damping of waves in the low frequency regime in boosted simulations was already suggested by Horst et al. (2020). When the luminosity is boosted, the increased damping implies that waves over a smaller range of frequencies reach the top of the cavity. Particularly, low frequency g modes are fully damped in boosted simulations compared to the one with realistic luminosity. In addition, waves of given frequency and angular degree will be damped out in different locations in boosted simulations, therefore depositing their energy in different regions. In a case where the luminosity is artificially enhanced but not the radiative diffusivity, as in the simulation of Horst et al. (2020), the location where IGWs are damped out is not modified. However, because of the higher luminosity, the dominant frequency range of excited waves will be different from in a non-boosted case (see Le Saux et al. 2022, and Section 4). Therefore, in simulations with different enhancement factors for the luminosity and radiative diffusivity, we do not expect the propagation of IGWs to be identical as in a model with realistic luminosity.

Finally, comparison between the two analytical cases with and without damping also highlights that running simulations with unrealistic radiative diffusivity may not capture the proper propagation properties of IGWs. Particularly for low frequency waves, which are the ones more impacted by radiative diffusion. In these simulations, transport by IGWs should be studied with caution as well as their ability to propagate up to the surface or not.

### 3.4 Radial kinetic energy density

Fig. 6 presents the power spectrum of the radial velocity $P[\hat{v}_r^2]$, as a function of normalized radius and frequency for an angular degree $ℓ = 20$. This quantity $P[\hat{v}_r^2]$ is equivalent to the radial kinetic energy density. This plot offers a general overview on the dependence of radiative damping of IGWs on frequency and radius. In the convective core, between the bottom of the plot and $r_{\text{conv}}$, the spectrum is relatively homogeneous at all frequencies and characteristic of a convective zone. The bright ridges observed in the radiative zone, at $r \ge 0.1814R_{\text{star}}$, are high-amplitude g modes. The dark knots observed in the bright ridges are the radial nodes of the considered mode. The number of nodes for a given mode defines the radial order $n$ of the mode, and it increases as the frequency decreases. This is an important characteristic of g modes (Aerts, Christensen-Dalsgaard & Kurtz 2010).
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Power spectrum of the radial velocity for the simulation *ref*. The angular degree is fixed at $ℓ = 20$. The power spectrum is obtained via mode projection on the spherical harmonics basis and a temporal Fourier transform of the radial velocity.</figcaption></figure>

The properties displayed in Fig. 6 for model *ref* are common to the three boosted simulations. We find that IGWs are dampened by radiative effects and more importantly that the damping strength depends on the location in radius, as can be expected from the radiative diffusivity profile showed in Fig. 2. The damping rate seems relatively constant between the convective boundary and $r \simeq 0.8 R_{\text{star}}$, but above this radius, wave damping appears to be strengthened. This means that waves deposit their energy in the radiative cavity but there are regions where this deposition will be more important, particularly close to the top of the simulation domain. This will be further investigated in Section 6.

## 4 WAVES ENERGY FLUX

As in Le Saux et al. (2022), we define the flux $F_{\text{wave}}$ for an individual IGW mode ($ℓ, ω$) at radius $r$ by

$$F_{\text{wave}}(r, ℓ, ω) \sim \frac{1}{2} \rho \frac{N}{k_h} P[\hat{v}_r^2](r, ℓ, ω). \quad (18)$$

Then we express it in a differential form in order to compare with theoretical predictions (see Le Saux et al. 2022, for details):

$$\frac{dF_{\text{wave}}}{d \ln \omega d \ln k_h} \sim \frac{1}{2} \rho T_s r N \omega P[\hat{v}_r^2], \quad (19)$$

with $T_s$ the sampling time.

Fig. 7 presents the dependence of the differential wave energy flux given by equation (19) with frequency at two radii, $r = r_{\text{conv}} + 0.5 H_{\text{p,conv}} \simeq 0.231 R_{\text{star}}$ (purple curves) and $r = r_{\text{conv}} + 6 H_{\text{p,conv}} \simeq 0.777 R_{\text{star}}$ (green curves). The two vertical dashed lines indicate the Brunt–Väisälä frequencies at the radii with the same colour code. The wave fluxes are plotted for an angular degree $ℓ = 10$. Note that for each simulation Fig. 7 shows a different frequency range as we focus on the range that bears most of the energy. The spectra present a relatively flat structure at low frequencies, then peak around a given frequency $\omega_{\text{peak}}$ and finally decrease towards higher frequencies. As done in Le Saux et al. (2022), these fluxes are compared with theoretical predictions for waves generation by Reynolds stress (blue, red and orange straight lines) and penetrative convection (cyan and black dashed curves). The analytical model used for Reynolds stress

excitation is from Lecoanet & Quataert (2013) and for penetrative convection, or plumes, excitation the model from Pinçon, Belkacem & Goupil (2016) is used. Lecoanet & Quataert (2013) predict the impact of the transition between the radiative and the convective regions on the wave flux. The dependence on frequency of the wave flux is modified when considering a different profile for the temperature gradient at the interface. In the case of discontinuous transition, the wave flux scales as $\omega^{-1/2}$ (dark blue solid line). If the transition is abrupt but continuous (piecewise linear) the wave flux scales with $\omega^{-41/6}$ (red solid line). Finally, if the transition is smooth (tanh profile), the wave flux scales as $\omega^{-15/2}$ (orange solid line). Pinçon et al. (2016) predict a Gaussian energy flux spectrum that scales as $e^{-\omega^2/4\nu_p^2}$, with $\nu_p$ the characteristic frequency associated with the plumes' lifetime. This frequency $\nu_p$ is difficult to estimate analytically or from simulations, but Pinçon et al. (2016) suggest that $\nu_p \sim \omega_{\text{conv}}$ (cyan dashed curve) is a good approximation. In our comparison in Fig. 7 we consider a second value for characteristic frequency associated with the plumes' lifetime $\nu_p = 2\omega_{\text{conv}}$ (black dashed curve) in order to test the influence of this parameter. For both models, the excitation mechanism is expected to generate waves with frequencies larger than the convective frequencies, $\omega \ge \omega_{\text{conv}}$. This is why we compare the wave flux from MUSIC to these analytical models in the range of frequencies $\omega \ge \omega_{\text{peak}}$. For a better visual comparison with theoretical predictions, we also compute and plot a running median on the 25th percentile of the fluxes with a window of 100 frequency bins (dark purple and green lines superimposed on the corresponding flux).

After IGWs are excited at the boundary by convection in the core at frequencies $\omega \ge \omega_{\text{conv}}$, they propagate away towards the surface. During this propagation, low frequency waves are damped much more rapidly than their higher-frequency counterparts. Indeed, as explained in Section 3.2 radiative damping of IGWs is modelled by equation (15) which scales as $\omega^{-4}$. At a given radius $r$, IGWs with frequencies in the range [$\omega_{\text{conv}}$, $\omega_{\text{peak}}(r)$] are already damped. Consequently, the peak of the flux is shifted towards higher frequencies at larger radii. This is illustrated in Fig. 7, when comparing the fluxes at two locations. For model *ref*, the peak is located at $\omega_{\text{peak}} \sim 2.5$ µHz at $r = 0.231 R_{\text{star}}$ (purple curve) and at $\omega_{\text{peak}} \sim 5$ µHz at $r = 0.777 R_{\text{star}}$ (green curve). Waves with frequencies between 2.5 and 5 µHz have been damped before being able to reach $r = 0.777 R_{\text{star}}$. For models *boost1d1*, *boost1d2*, and *boost1d4*, the shift is from 4 to 7 µHz, 7 to 10 µHz, and 11 to 13 µHz respectively.

For frequencies larger than $\omega_{\text{peak}}$, the fluxes of the four simulations are decreasing towards higher frequencies up to the Brunt–Väisälä frequency. Because of the high-amplitude narrow peaks, which are g modes, it is not possible to precisely measure the slope of the spectra. However, for models *ref*, *boost1d1*, and *boost1d2*, the wave flux measured in MUSIC is broadly consistent with analytical prediction for Reynolds stress excitation from Lecoanet & Quataert (2013) and for plumes excitation from Pinçon et al. (2016), but in different frequency ranges. Close to the peak the excitation seems to be dominated by penetrative convection, whereas the Reynolds stress takes over at larger frequencies. This result was already suggested by Pinçon et al. (2016) (see their fig. 3). In addition, it seems that the fit is better when using $\nu_p = 2\omega_{\text{conv}}$. However, it is difficult to disentangle the two mechanisms since Reynolds stress and penetrative convection act simultaneously to excite waves in the same frequency range, namely between $\omega_{\text{conv}}$ and $N$.

Model *boost1d4* shows a different behaviour. The MUSIC flux has a Gaussian shape that is broadly consistent with the plumes' excitation, but the proximity of the peak of the flux to the Brunt–
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Power spectrum of the radial velocity for the simulation *ref*. The angular degree is fixed at $ℓ = 20$. The power spectrum is obtained via mode projection on the spherical harmonics basis and a temporal Fourier transform of the radial velocity.</figcaption></figure>

The properties displayed in Fig. 6 for model *ref* are common to the three boosted simulations. We find that IGWs are dampened by radiative effects and more importantly that the damping strength depends on the location in radius, as can be expected from the radiative diffusivity profile showed in Fig. 2. The damping rate seems relatively constant between the convective boundary and $r \simeq 0.8 R_{\text{star}}$, but above this radius, wave damping appears to be strengthened. This means that waves deposit their energy in the radiative cavity but there are regions where this deposition will be more important, particularly close to the top of the simulation domain. This will be further investigated in Section 6.

## 4 WAVES ENERGY FLUX

As in Le Saux et al. (2022), we define the flux $F_{\text{wave}}$ for an individual IGW mode ($ℓ, ω$) at radius $r$ by

$$F_{\text{wave}}(r, \ell, \omega) \sim \frac{1}{2} \rho \frac{N}{k_h} P[\hat{\nu}_r^2](r, \ell, \omega). \qquad (18)$$

Then we express it in a differential form in order to compare with theoretical predictions (see Le Saux et al. 2022, for details):

$$\frac{dF_{\text{wave}}}{d \ln \omega d \ln k_h} \sim \frac{1}{2} \rho T_s r N \omega P[\hat{\nu}_r^2], \qquad (19)$$

with $T_s$ the sampling time.

Fig. 7 presents the dependence of the differential wave energy flux given by equation (19) with frequency at two radii, $r = r_{\text{conv}} + 0.5 H_{\text{p,conv}} \simeq 0.231 R_{\text{star}}$ (purple curves) and $r = r_{\text{conv}} + 6 H_{\text{p,conv}} \simeq 0.777 R_{\text{star}}$ (green curves). The two vertical dashed lines indicate the Brunt–Väisälä frequencies at the radii with the same colour code. The wave fluxes are plotted for an angular degree $ℓ = 10$. Note that for each simulation Fig. 7 shows a different frequency range as we focus on the range that bears most of the energy. The spectra present a relatively flat structure at low frequencies, then peak around a given frequency $\omega_{\text{peak}}$ and finally decrease towards higher frequencies. As done in Le Saux et al. (2022), these fluxes are compared with theoretical predictions for waves generation by Reynolds stress (blue, red and orange straight lines) and penetrative convection (cyan and black dashed curves). The analytical model used for Reynolds stress

excitation is from Lecoanet & Quataert (2013) and for penetrative convection, or plumes, excitation the model from Pinçon, Belkacem & Goupil (2016) is used. Lecoanet & Quataert (2013) predict the impact of the transition between the radiative and the convective regions on the wave flux. The dependence on frequency of the wave flux is modified when considering a different profile for the temperature gradient at the interface. In the case of discontinuous transition, the wave flux scales as $\omega^{-1/2}$ (dark blue solid line). If the transition is abrupt but continuous (piecewise linear) the wave flux scales with $\omega^{-4/6}$ (red solid line). Finally, if the transition is smooth (tanh profile), the wave flux scales as $\omega^{-1/2}$ (orange solid line). Pinçon et al. (2016) predict a Gaussian energy flux spectrum that scales as $e^{-\omega^2/4\nu_p^2}$, with $\nu_p$ the characteristic frequency associated with the plumes' lifetime. This frequency $\nu_p$ is difficult to estimate analytically or from simulations, but Pinçon et al. (2016) suggest that $\nu_p \sim \omega_{\text{conv}}$ (cyan dashed curve) is a good approximation. In our comparison in Fig. 7 we consider a second value for characteristic frequency associated with the plumes' lifetime $\nu_p = 2\omega_{\text{conv}}$ (black dashed curve) in order to test the influence of this parameter. For both models, the excitation mechanism is expected to generate waves with frequencies larger than the convective frequencies, $\omega \ge \omega_{\text{conv}}$. This is why we compare the wave flux from MUSIC to these analytical models in the range of frequencies $\omega \ge \omega_{\text{peak}}$. For a better visual comparison with theoretical predictions, we also compute and plot a running median on the 25th percentile of the fluxes with a window of 100 frequency bins (dark purple and green lines superimposed on the corresponding flux).

After IGWs are excited at the boundary by convection in the core at frequencies $\omega \ge \omega_{\text{conv}}$, they propagate away towards the surface. During this propagation, low frequency waves are damped much more rapidly than their higher-frequency counterparts. Indeed, as explained in Section 3.2 radiative damping of IGWs is modelled by equation (15) which scales as $\omega^{-4}$. At a given radius $r$, IGWs with frequencies in the range [$\omega_{\text{conv}}$, $\omega_{\text{peak}}(r)$] are already damped. Consequently, the peak of the flux is shifted towards higher frequencies at larger radii. This is illustrated in Fig. 7, when comparing the fluxes at two locations. For model *ref*, the peak is located at $\omega_{\text{peak}} \sim 2.5$ µHz at $r = 0.231 R_{\text{star}}$ (purple curve) and at $\omega_{\text{peak}} \sim 5$ µHz at $r = 0.777 R_{\text{star}}$ (green curve). Waves with frequencies between 2.5 and 5 µHz have been damped before being able to reach $r = 0.777 R_{\text{star}}$. For models *boost1d1*, *boost1d2*, and *boost1d4*, the shift is from 4 to 7 µHz, 7 to 10 µHz, and 11 to 13 µHz respectively.

For frequencies larger than $\omega_{\text{peak}}$, the fluxes of the four simulations are decreasing towards higher frequencies up to the Brunt–Väisälä frequency. Because of the high-amplitude narrow peaks, which are g modes, it is not possible to precisely measure the slope of the spectra. However, for models *ref*, *boost1d1*, and *boost1d2*, the wave flux measured in MUSIC is broadly consistent with analytical prediction for Reynolds stress excitation from Lecoanet & Quataert (2013) and for plumes excitation from Pinçon et al. (2016), but in different frequency ranges. Close to the peak the excitation seems to be dominated by penetrative convection, whereas the Reynolds stress takes over at larger frequencies. This result was already suggested by Pinçon et al. (2016) (see their fig. 3). In addition, it seems that the fit is better when using $\nu_p = 2\omega_{\text{conv}}$. However, it is difficult to disentangle the two mechanisms since Reynolds stress and penetrative convection act simultaneously to excite waves in the same frequency range, namely between $\omega_{\text{conv}}$ and $N$.

Model *boost1d4* shows a different behaviour. The MUSIC flux has a Gaussian shape that is broadly consistent with the plumes' excitation, but the proximity of the peak of the flux to the Brunt–
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Wave energy flux as a function of frequency for an angular degree $\ell = 10$ for the four simulations *ref* (top left), *boost1d1* (top right), *boost1d2* (bottom left), and *boost1d4* (bottom right) at two radii, $r = r_{CB} + 0.5H_{p,conv} \approx 0.231 R_{star}$ (purple curves) and $r = r_{CB} + 6H_{p,conv} \approx 0.777 R_{star}$ (green curves). The vertical dashed lines indicate the values of the Brunt–Väisälä at $r \approx 0.231 R_{star}$ (purple curves) and $r \approx 0.777 R_{star}$ (green curves). The spectra are obtained via mode projection on the spherical harmonics basis and a temporal Fourier transform of the radial velocity.</figcaption></figure>

Väisälä frequency makes the comparison difficult. Nevertheless, it seems that for this simulation the excitation by penetrative convection is more efficient on a larger frequency range. The shape of the curve is no longer consistent with any of the three predictions from Lecoanet & Quataert (2013) in any frequency range. This suggests that the excitation of IGWs by penetrative flows is strengthened when the luminosity is increased, and seems to dominate over Reynolds stress excitation in the most boosted simulation. Interestingly, this shape of spectrum is similar to results from other multidimensional simulations, such as the ones from Rogers et al. (2013) or Edelmann et al. (2019). In both these studies, it is suggested that the excitations of IGWs is dominated by penetrative convection. Note that the authors artificially increase the luminosity of their models by factors up to $10^7$. This suggests that enhancing the luminosity by large factors tend to increase the efficiency of plume excitation.

Finally, for the three models *ref*, *boost1d1* and *boost1d2* the slope of the wave flux remains similar close to the convective boundary (purple curve) and at the top of the domain (green curve). However, it

is difficult to draw conclusions about the slope up to the stellar surface because of the wave damping that is significantly increasing from $r \approx 0.8 R_{tot}$ as shown in Fig. 6 and also because of the complexity of the near surface layers that will affect the waves. Moreover, when analysing hydrodynamical simulations with different values for the luminosity enhancement factor, waves of different frequencies would be able to reach the stellar surface. As we have seen earlier, higher frequencies waves will be excited with a larger amplitude and low frequencies waves are more strongly damped, when a larger enhancement factor is used. Comparing results from luminosity enhanced models with observations requires caution, since the former predict the wrong range of frequencies for waves that could reach the stellar surface. Extrapolation of spectra measured in the interior of the stellar model to the surface of the star may not be straightforward due to the strong impact of the near surface layers on waves propagation. Indeed, there is a very important increase of the radiative diffusivity and the Brunt–Väisälä frequency in this region. Consequently, we suggest that quantitative direct comparison between observations
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Wave energy flux as a function of frequency for an angular degree $\ell = 10$ for the four simulations *ref* (top left), *boost1d1* (top right), *boost1d2* (bottom left), and *boost1d4* (bottom right) at two radii, $r = r_{CB} + 0.5H_{p,conv} \approx 0.231 R_{star}$ (purple curves) and $r = r_{CB} + 6H_{p,conv} \approx 0.777 R_{star}$ (green curves). The vertical dashed lines indicate the values of the Brunt–Väisälä at $r \approx 0.231 R_{star}$ (purple curves) and $r \approx 0.777 R_{star}$ (green curves). The spectra are obtained via mode projection on the spherical harmonics basis and a temporal Fourier transform of the radial velocity.</figcaption></figure>

Väisälä frequency makes the comparison difficult. Nevertheless, it seems that for this simulation the excitation by penetrative convection is more efficient on a larger frequency range. The shape of the curve is no longer consistent with any of the three predictions from Lecoanet & Quataert (2013) in any frequency range. This suggests that the excitation of IGWs by penetrative flows is strengthened when the luminosity is increased, and seems to dominate over Reynolds stress excitation in the most boosted simulation. Interestingly, this shape of spectrum is similar to results from other multidimensional simulations, such as the ones from Rogers et al. (2013) or Edelmann et al. (2019). In both these studies, it is suggested that the excitations of IGWs is dominated by penetrative convection. Note that the authors artificially increase the luminosity of their models by factors up to $10^7$. This suggests that enhancing the luminosity by large factors tend to increase the efficiency of plume excitation.

Finally, for the three models *ref*, *boost1d1* and *boost1d2* the slope of the wave flux remains similar close to the convective boundary (purple curve) and at the top of the domain (green curve). However, it

is difficult to draw conclusions about the slope up to the stellar surface because of the wave damping that is significantly increasing from $r \approx 0.8 R_{tot}$ as shown in Fig. 6 and also because of the complexity of the near surface layers that will affect the waves. Moreover, when analysing hydrodynamical simulations with different values for the luminosity enhancement factor, waves of different frequencies would be able to reach the stellar surface. As we have seen earlier, higher frequencies waves will be excited with a larger amplitude and low frequencies waves are more strongly damped, when a larger enhancement factor is used. Comparing results from luminosity enhanced models with observations requires caution, since the former predict the wrong range of frequencies for waves that could reach the stellar surface. Extrapolation of spectra measured in the interior of the stellar model to the surface of the star may not be straightforward due to the strong impact of the near surface layers on waves propagation. Indeed, there is a very important increase of the radiative diffusivity and the Brunt–Väisälä frequency in this region. Consequently, we suggest that quantitative direct comparison between observations
and simulations would require to run simulations as close as possible to the stellar surface. Moreover, observations of stellar oscillations do not resolve the surface of stars and consequently the signal is averaged in a way that makes it difficult to compare to simulations.

In our work, we use radial velocity spectra. Lecoanet et al. (2021) suggest that it should be equivalent to use any local wave perturbation variable. We have checked this suggestion by computing spectra using temperature perturbations and co-latitudinal (horizontal) velocity and find a very good agreement with the ones computed using the radial velocity. From an observational perspective, it can be interesting to estimate the ratio of horizontal to radial velocities. At a location $r = 0.85R_{\text{star}}$ in model ref, we calculate $v_h/v_r \sim 270$ for $\omega = 10 \, \mu\text{Hz}$ and $v_h/v_r \sim 13$ for $\omega = 40 \, \mu\text{Hz}$. This decrease with frequency is expected from the dispersion relation of IGWs, which predicts that the ratio $v_h/v_r$ varies as $N/\omega$, with the Brunt-Väisälä frequency $N$ which is fixed at a given radius. The values we obtain for the ratio of the velocities are in agreement with the ones determined in the simulations of Horst et al. (2020) in the same frequency range. At lower frequency, the value of the ratio keeps increasing, and for $\omega = 2 \, \mu\text{Hz}$ we obtain $v_h/v_r \sim 7800$. This value is larger than the ones calculated by Horst et al. (2020) at similar frequency. We suspect that this discrepancy is the result of the lack of independent data points from their simulations for the temporal Fourier transform, as suggested by the authors. Such high value of $v_h/v_r$ is in agreement with the value calculated from the 2D simulations of Aerts & Rogers (2015). In their simulations of rotating intermediate-mass stars, Aerts & Rogers (2015) obtain ratio up to $10^4$ at $r = 0.99R_{\text{star}}$. In the observational community, this ratio is known as the K value which is approximated by $K \approx GM/4\pi^2 v^2 R^3$, with $v$ the observed intrinsic frequency of a given star and $G$ is the gravitational constant. In their study, De Cat & Aerts (2002) measure this K value for SPB stars, which are mid-B type stars pulsating in high-order g modes. They obtain typical values between $\sim 10$ and $\sim 100$. However, note that this approximated K value is only defined at the stellar surface. Consequently, direct comparison with simulated velocity amplitudes determined deeper in the stellar interior should be taken very cautiously.

In their study, Bowman et al. (2019) are using luminosity perturbations, which is equivalent to look at the perturbations of the effective temperature $T_{\text{eff}}^4$ as we have the relation $L \propto T_{\text{eff}}^4$. First, in both spectra, observed and modelled, we note the presence of g modes, appearing as high-amplitude narrow peaks. In the observed spectra, these peaks are present at frequencies larger than the so-called low-frequency power excess. We suggest that if this low-frequency power excess results from IGWs excited by core convection, we should expect to see g modes (i.e. narrow peaks in the spectra) in this low-frequency range. Indeed, in the spectrum of model ref, there are g modes starting to appear from $10 \, \mu\text{Hz}$. As we will see in Section 5, in this simulation, waves with frequencies lower than $10 \, \mu\text{Hz}$ are damped before being able to reach the stellar surface. However, as suggested by Edelmann et al. (2019) and Horst et al. (2020) it is also possible that the small frequency spacing between modes of different radial order and different angular degree ‘hides’ these individual narrow peaks. To confirm this suggestion would require simulations with a radial domain extending to layers close to the stellar surface. Note that the radial extent in Edelmann et al. (2019) and Horst et al. (2020) is at 90 per cent and 91 per cent of the stellar radius, respectively. Second, the simulation spectrum shape peaks around $\omega_{\text{peak}}$ and decreases towards lower frequencies. This feature is not observed, and the spectra inferred from observations remain mostly flat at low frequencies. This difference was already reported by Edelmann et al. (2019) and Lecoanet et al. (2021). In their study,

Horst et al. (2020) state that this drop towards low frequencies is not present in the spectra measured in their simulations. They suggest that this could be attributed to the low viscosity and thermal diffusivity used in their simulations. However, our results show that this feature should be present even for non-boosted simulations. We suggest that this drop is still present in the simulation of Horst et al. (2020) but it has a very low amplitude (see their fig. 16). This is a consequence of boosting only the luminosity and not the radiative diffusivity, which will results in less damping in the frequency range close to $\omega_{\text{peak}}$. Concerning the observed spectra from Bowman et al. (2019, 2020), this drop at low frequency is not observed. This could be the result of rotational effects that could shift wave frequency to lower values. In this case, rotation could help reconcile our simulations and observations. Rogers et al. (2013) indeed suggest an important impact of differential rotation for frequencies below $10 \, \mu\text{Hz}$, based on simulations which initially impose some differential rotation. Further work including rotation and an appropriate modelling of the radiative damping of waves as close as possible to the surface layers would be required to confirm this effect of rotation.

## 5 NON-LINEAR EFFECTS OF IGWS

Linear theory for wave propagation is based on the assumption that the displacement amplitude of a wave is small compared to its wavelength. As the amplitude increases, non-linear effects can arise, such that the advection terms in the material derivative of the momentum equation start to play an important role (Sutherland 2010). These non-linear effects can modify wave properties and evolution. To quantify the impact of non-linear effects, we use a non-linear parameter $\epsilon$, defined as (Press 1981)

$$ \epsilon := \xi_h k_h = \xi_r k_r, \qquad (20) $$

where $\xi_h$ (resp. $\xi_r$) is the horizontal (resp. radial) displacement of a fluid particle associated with an IGW. This definition is basically the ratio of the horizontal (resp. radial) amplitude of a wave to its horizontal (resp. radial) wavelength. According to Press (1981), the displacement may be approximated as $\xi_i = v_i/\omega$ where $i = r, h$. We thus have for the non-linearity parameter

$$ \epsilon = \frac{v_h}{\omega} k_h = \frac{v_r}{\omega} k_r. \qquad (21) $$

From this definition, an IGW is non-linear if $\epsilon \ge 1$. In the following, we compute $\epsilon$ using the definition based on the radial component of velocity and wavenumber. Using equation (17), an analytical expression for $\epsilon$ can be inferred

$$ \epsilon = v_0(\ell, \omega) \omega^{-1} \left(\frac{\rho}{\rho_0}\right)^{-1/2} \left(\frac{k_h}{k_{h,0}}\right)^{3/2} \left(\frac{N^2 - \omega^2}{N_0^2 - \omega^2}\right)^{-1/4} e^{-\tau/2} k_r. \qquad (22) $$

Using the dispersion relation equation (12) for IGWs, we can calculate the radial wavenumber

$$ k_r = k_h \left( \frac{N^2}{\omega^2} - 1 \right)^{1/2} . \qquad (23) $$

Fig. 8 presents the dependence of this non-linear parameter $\epsilon$ on frequency and radius for an angular degree $\ell = 4$, calculated from the MUSIC simulations (left-hand column) and as predicted from theory (right-hand column) using equation (22). We chose to analyse the angular degree $\ell = 4$ as it is for this degree that Horst et al. (2020) predict that the non-linear effects should be most important. We have performed the same analysis for angular degrees $\ell = 1, 3, 5, 10$, and 20, and we obtained similar results. Note that the radial range
and simulations would require to run simulations as close as possible to the stellar surface. Moreover, observations of stellar oscillations do not resolve the surface of stars and consequently the signal is averaged in a way that makes it difficult to compare to simulations.

In our work, we use radial velocity spectra. Lecoanet et al. (2021) suggest that it should be equivalent to use any local wave perturbation variable. We have checked this suggestion by computing spectra using temperature perturbations and co-latitudinal (horizontal) velocity and find a very good agreement with the ones computed using the radial velocity. From an observational perspective, it can be interesting to estimate the ratio of horizontal to radial velocities. At a location $r = 0.85R_{\text{star}}$ in model ref, we calculate $v_h/v_r \sim 270$ for $\omega = 10 \, \mu\text{Hz}$ and $v_h/v_r \sim 13$ for $\omega = 40 \, \mu\text{Hz}$. This decrease with frequency is expected from the dispersion relation of IGWs, which predicts that the ratio $v_h/v_r$ varies as $N/\omega$, with the Brunt-Väisälä frequency $N$ which is fixed at a given radius. The values we obtain for the ratio of the velocities are in agreement with the ones determined in the simulations of Horst et al. (2020) in the same frequency range. At lower frequency, the value of the ratio keeps increasing, and for $\omega = 2 \, \mu\text{Hz}$ we obtain $v_h/v_r \sim 7800$. This value is larger than the ones calculated by Horst et al. (2020) at similar frequency. We suspect that this discrepancy is the result of the lack of independent data points from their simulations for the temporal Fourier transform, as suggested by the authors. Such high value of $v_h/v_r$ is in agreement with the value calculated from the 2D simulations of Aerts & Rogers (2015). In their simulations of rotating intermediate-mass stars, Aerts & Rogers (2015) obtain ratio up to $10^4$ at $r = 0.99R_{\text{star}}$. In the observational community, this ratio is known as the K value which is approximated by $K \approx GM/4\pi^2 v^2 R^3$, with $v$ the observed intrinsic frequency of a given star and $G$ is the gravitational constant. In their study, De Cat & Aerts (2002) measure this K value for SPB stars, which are mid-B type stars pulsating in high-order g modes. They obtain typical values between $\sim 10$ and $\sim 100$. However, note that this approximated K value is only defined at the stellar surface. Consequently, direct comparison with simulated velocity amplitudes determined deeper in the stellar interior should be taken very cautiously.

In their study, Bowman et al. (2019) are using luminosity perturbations, which is equivalent to look at the perturbations of the effective temperature $T_{\text{eff}}^4$ as we have the relation $L \propto T_{\text{eff}}^4$. First, in both spectra, observed and modelled, we note the presence of g modes, appearing as high-amplitude narrow peaks. In the observed spectra, these peaks are present at frequencies larger than the so-called low-frequency power excess. We suggest that if this low-frequency power excess results from IGWs excited by core convection, we should expect to see g modes (i.e. narrow peaks in the spectra) in this low-frequency range. Indeed, in the spectrum of model ref, there are g modes starting to appear from $10 \, \mu\text{Hz}$. As we will see in Section 5, in this simulation, waves with frequencies lower than $10 \, \mu\text{Hz}$ are damped before being able to reach the stellar surface. However, as suggested by Edelmann et al. (2019) and Horst et al. (2020) it is also possible that the small frequency spacing between modes of different radial order and different angular degree 'hides' these individual narrow peaks. To confirm this suggestion would require simulations with a radial domain extending to layers close to the stellar surface. Note that the radial extent in Edelmann et al. (2019) and Horst et al. (2020) is at 90 per cent and 91 per cent of the stellar radius, respectively. Second, the simulation spectrum shape peaks around $\omega_{\text{peak}}$ and decreases towards lower frequencies. This feature is not observed, and the spectra inferred from observations remain mostly flat at low frequencies. This difference was already reported by Edelmann et al. (2019) and Lecoanet et al. (2021). In their study,

Horst et al. (2020) state that this drop towards low frequencies is not present in the spectra measured in their simulations. They suggest that this could be attributed to the low viscosity and thermal diffusivity used in their simulations. However, our results show that this feature should be present even for non-boosted simulations. We suggest that this drop is still present in the simulation of Horst et al. (2020) but it has a very low amplitude (see their fig. 16). This is a consequence of boosting only the luminosity and not the radiative diffusivity, which will results in less damping in the frequency range close to $\omega_{\text{peak}}$. Concerning the observed spectra from Bowman et al. (2019, 2020), this drop at low frequency is not observed. This could be the result of rotational effects that could shift wave frequency to lower values. In this case, rotation could help reconcile our simulations and observations. Rogers et al. (2013) indeed suggest an important impact of differential rotation for frequencies below $10 \, \mu\text{Hz}$, based on simulations which initially impose some differential rotation. Further work including rotation and an appropriate modelling of the radiative damping of waves as close as possible to the surface layers would be required to confirm this effect of rotation.

## 5 NON-LINEAR EFFECTS OF IGWS

Linear theory for wave propagation is based on the assumption that the displacement amplitude of a wave is small compared to its wavelength. As the amplitude increases, non-linear effects can arise, such that the advection terms in the material derivative of the momentum equation start to play an important role (Sutherland 2010). These non-linear effects can modify wave properties and evolution. To quantify the impact of non-linear effects, we use a non-linear parameter $\epsilon$, defined as (Press 1981)

$$ \epsilon := \xi_h k_h = \xi_r k_r, \qquad (20) $$

where $\xi_h$ (resp. $\xi_r$) is the horizontal (resp. radial) displacement of a fluid particle associated with an IGW. This definition is basically the ratio of the horizontal (resp. radial) amplitude of a wave to its horizontal (resp. radial) wavelength. According to Press (1981), the displacement may be approximated as $\xi_i = v_i/\omega$ where $i = r, h$. We thus have for the non-linearity parameter

$$ \epsilon = \frac{v_h}{\omega} k_h = \frac{v_r}{\omega} k_r. \qquad (21) $$

From this definition, an IGW is non-linear if $\epsilon \ge 1$. In the following, we compute $\epsilon$ using the definition based on the radial component of velocity and wavenumber. Using equation (17), an analytical expression for $\epsilon$ can be inferred

$$ \epsilon = v_0(\ell, \omega) \omega^{-1} \left(\frac{\rho}{\rho_0}\right)^{-1/2} \left(\frac{k_h}{k_{h,0}}\right)^{3/2} \left(\frac{N^2 - \omega^2}{N_0^2 - \omega^2}\right)^{-1/4} e^{-\tau/2} k_r. \qquad (22) $$

Using the dispersion relation equation (12) for IGWs, we can calculate the radial wavenumber

$$ k_r = k_h \left( \frac{N^2}{\omega^2} - 1 \right)^{1/2} . \qquad (23) $$

Fig. 8 presents the dependence of this non-linear parameter $\epsilon$ on frequency and radius for an angular degree $\ell = 4$, calculated from the MUSIC simulations (left-hand column) and as predicted from theory (right-hand column) using equation (22). We chose to analyse the angular degree $\ell = 4$ as it is for this degree that Horst et al. (2020) predict that the non-linear effects should be most important. We have performed the same analysis for angular degrees $\ell = 1, 3, 5, 10$, and 20, and we obtained similar results. Note that the radial range
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Dependence of the non-linear parameter $\epsilon$ on frequency and normalized radius for the four simulations *ref* (top row), *boost1d1* (second row), *boost1d2* (third row), and *boost1d4* (bottom row). The left-hand column presents $\epsilon$ as measured in MUSIC and the right-hand column comes from the analytical expression given by equation (22), normalized with the value of the velocity in the simulations at $r = r_N$, where $r_N$ is the smallest radius such as $\omega = N(r_N)$. The vertical white line in the analytical plots indicates the top of the numerical domain. The horizontal grey dashed line indicates the convective frequency $\omega_{conv}$ for each simulation.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Dependence of the non-linear parameter $\epsilon$ on frequency and normalized radius for the four simulations *ref* (top row), *boost1d1* (second row), *boost1d2* (third row), and *boost1d4* (bottom row). The left-hand column presents $\epsilon$ as measured in MUSIC and the right-hand column comes from the analytical expression given by equation (22), normalized with the value of the velocity in the simulations at $r = r_N$, where $r_N$ is the smallest radius such as $\omega = N(r_N)$. The vertical white line in the analytical plots indicates the top of the numerical domain. The horizontal grey dashed line indicates the convective frequency $\omega_{conv}$ for each simulation.</figcaption></figure>
(x-axis) is different in the two columns, in the left one it extends to the top of the numerical domain $r_{\text{out}} = 0.91 R_{\text{star}}$ and in the right one it extends up to the surface $r = R_{\text{star}}$. For the theoretical plots, the value of $v_0$ needed to normalize the amplitude of the velocity is the value of the radial velocity $v(r_N, \ell, \omega)$ where $r_N$ is the smallest radius, such as $\omega = N(r_N)$ for a given frequency $\omega$. For frequencies up to ~50 µHz, $r_N \simeq 0.183 R_{\text{star}}$, as in Section 3.2. For frequency between ~50 µHz and the maximal value of the Brunt-Väisälä, $N_{\text{max}} \sim 220$ µHz, this radius $r_N$ increases with frequency (see Fig. 1). We introduce $r_N$ because of the condition for IGWs propagation $\omega < N$. Note also that the frequency range (x-axis) available for each simulation is different. The lower frequency is set to 0.6 µHz for the four simulations whereas the maximal frequencies are 50, 100, 200, and 200 µHz for ref, boost1d1, boost1d2, and boost1d4, respectively. This is because we focus on the frequency range that bears most of the energy in each simulation. We have also set a minimum threshold for $\epsilon$ at $10^{-10}$, as the actual value can become very small due to the exponential damping term.

All plots present a similar general aspect, with horizontal ridges corresponding to waves of a given frequency and for which the value of $\epsilon$ varies with radius. For each simulation, there is a clearly defined range of frequencies above the convective frequency (horizontal dashed grey line) with $\epsilon \ge 10^{-3}$. This occurs in different frequency ranges for each simulation. These ranges are approximately [1.0, 5.0], [1.2, 9.0], [2.0, 20.0], and [7.0, 50.0] µHz for models ref, boost1d1, boost1d2, and boost1d4, respectively. For most of these waves, we can see that they seem to conserve their structure in the whole envelope, suggesting that no non-linear effects occur (no energy transfer, mode coupling, etc.).

This is not the case for the lowest frequencies in the bottom region of the radiative envelope, just above the convective core. Indeed, in this region for frequencies close to the convective frequency, we can see that $\epsilon \ge 10^{-2}$, and even $\epsilon \ge 10^{-1}$ for some frequencies, both in the simulations and the theoretical plots. Ratnasingam et al. (2019) suggested that we can observe non-linear effects from $\epsilon \sim 0.1$. Therefore, we may expect non-linear effects to be relevant just above the convective core. It is possible that some IGWs generated at these frequencies have too large amplitudes and break close to the boundary.

Now, if we compare the simulations and theoretical plots, there is a major difference at very low frequencies $\omega \le \omega_{\text{conv}}$. The simulations plots show high values of $\epsilon$ in the whole radiative zone for this frequency range, which do not appear on the theoretical plots and do not present the structure of horizontal ridges as at higher frequencies. This signal with large values of $\epsilon$ is difficult to analyse as it is localized at very low frequencies, for which the wavelengths of IGWs approach the spatial resolution of our grid. Indeed, as the radial wavelength approaches twice the length of a grid cell, it is not possible any more to represent the wave on the grid (see also section 6.2 of Le Saux et al. 2022).

If we compare the four simulations, $\epsilon$ is larger as the luminosity enhancement factor is increased. This is expected as the amplitude of the waves increase with the boost. Consequently, we should expect non-linear effects to be more relevant in boosted simulations. This could result in more mixing (see e.g. Jermyn 2022), angular momentum transport (see e.g. Gervais et al. 2018) and/or wave-wave interactions.

Finally, by looking at the analytical predictions for $\epsilon$ up to $r = R_{\text{star}}$, we can see that non-linear effects are not expected close to the surface of the star. Most of the waves in the range of frequencies with values of $\epsilon \ge 10^{-3}$ are damped before even reaching $r = 0.91 R_{\text{star}}$, the top of the MUSIC radial domain. For waves that can propagate

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Spatially average temperature profile in the four simulations ref (orange), boost1d1 (indigo), boost1d2 (yellow), and boost1d4 (green). The dashed curves are the temperature profiles obtained at time $t_0$, and the plain curves at $t_0 + \Delta t$, where $\Delta t$ is different for each simulation (see the text for values). The temperature profile from the 1D initial model is represented by the plain blue curve. The outer boundary of the MUSIC numerical domain is indicated by the vertical dashed black line. The right axis corresponds to the term $k_{\text{rad}}N^3$ (blue dotted curve) computed from the 1D model.</figcaption></figure>

further, the maximal values of the non-linear parameter close to the surface are smaller than $10^{-4}$. This is consistent with the results of Ratnasingam et al. (2019) who computed a non-linearity parameter from the analytical spectra of Kumar, Talon & Zahn (1999) and Lecoanet & Quataert (2013) for Reynolds stress excitation of IGWs and found that $\epsilon < 10^{-2}$ in a 3M$_{\odot}$ star for the angular degree $\ell = 10$. However, this does not agree with the results of Horst et al. (2020) who conclude that non-linear effects may be expected at the surface of the star. We suggest that the high values of the non-linear parameter they measure in their simulation is a result of the underestimation of the radiative diffusion. Indeed, as only the luminosity of their model is enhanced, the resulting larger amplitude of IGWs is not compensated by an increase in damping.

Finally, according to the theoretical plots, some waves may reach the surface of the star, but these are relatively high frequency waves compared to the excitation frequency $\omega_{\text{conv}}$. The predictions, however, are based on the linear theory which neglects the interaction with a possible subsurface convection zone (see e.g. Cantiello & Braithwaite 2019) and the near-surface layers that are difficult to accurately model (see e.g. Basu & Chaplin 2017).

# 6 WAVE HEATING IN THE UPPER LAYERS

Often overlooked, thermal effects of IGWs can be significant. In the Earth atmosphere, they are known to irreversibly convert kinetic energy into internal energy (see e.g. Medvedev & Klaassen 2003). Through this mechanism, they can heat up the thermosphere (Yiğit & Medvedev 2009). Almost always neglected in main-sequence stars, these thermal effects of acoustic and internal gravity waves have recently been accounted for in evolved stars in order to explain outbursts in supernova progenitors (Fuller 2017; Wu & Fuller 2022) and in an attempt to explain the lithium enhancement of clump stars (Jermyn & Fuller 2022). In these studies, waves deposit heat through radiative diffusion.

## 6.1 Temperature increase in the simulations

In our simulations we observe a significant increase of the temperature in the upper layers of the model. Fig. 9 presents the evolution of the radial profile of the average temperature in our four
(x-axis) is different in the two columns, in the left one it extends to the top of the numerical domain $r_{\text{out}} = 0.91 R_{\text{star}}$ and in the right one it extends up to the surface $r = R_{\text{star}}$. For the theoretical plots, the value of $v_0$ needed to normalize the amplitude of the velocity is the value of the radial velocity $v(r_N, \ell, \omega)$ where $r_N$ is the smallest radius, such as $\omega = N(r_N)$ for a given frequency $\omega$. For frequencies up to ~50 µHz, $r_N \simeq 0.183 R_{\text{star}}$, as in Section 3.2. For frequency between ~50 µHz and the maximal value of the Brunt-Väisälä, $N_{\text{max}} \sim 220$ µHz, this radius $r_N$ increases with frequency (see Fig. 1). We introduce $r_N$ because of the condition for IGWs propagation $\omega < N$. Note also that the frequency range (x-axis) available for each simulation is different. The lower frequency is set to 0.6 µHz for the four simulations whereas the maximal frequencies are 50, 100, 200, and 200 µHz for ref, boost1d1, boost1d2, and boost1d4, respectively. This is because we focus on the frequency range that bears most of the energy in each simulation. We have also set a minimum threshold for $\epsilon$ at $10^{-10}$, as the actual value can become very small due to the exponential damping term.

All plots present a similar general aspect, with horizontal ridges corresponding to waves of a given frequency and for which the value of $\epsilon$ varies with radius. For each simulation, there is a clearly defined range of frequencies above the convective frequency (horizontal dashed grey line) with $\epsilon \ge 10^{-3}$. This occurs in different frequency ranges for each simulation. These ranges are approximately [1.0, 5.0], [1.2, 9.0], [2.0, 20.0], and [7.0, 50.0] µHz for models ref, boost1d1, boost1d2, and boost1d4, respectively. For most of these waves, we can see that they seem to conserve their structure in the whole envelope, suggesting that no non-linear effects occur (no energy transfer, mode coupling, etc.).

This is not the case for the lowest frequencies in the bottom region of the radiative envelope, just above the convective core. Indeed, in this region for frequencies close to the convective frequency, we can see that $\epsilon \ge 10^{-2}$, and even $\epsilon \ge 10^{-1}$ for some frequencies, both in the simulations and the theoretical plots. Ratnasingam et al. (2019) suggested that we can observe non-linear effects from $\epsilon \sim 0.1$. Therefore, we may expect non-linear effects to be relevant just above the convective core. It is possible that some IGWs generated at these frequencies have too large amplitudes and break close to the boundary.

Now, if we compare the simulations and theoretical plots, there is a major difference at very low frequencies $\omega \le \omega_{\text{conv}}$. The simulations plots show high values of $\epsilon$ in the whole radiative zone for this frequency range, which do not appear on the theoretical plots and do not present the structure of horizontal ridges as at higher frequencies. This signal with large values of $\epsilon$ is difficult to analyse as it is localized at very low frequencies, for which the wavelengths of IGWs approach the spatial resolution of our grid. Indeed, as the radial wavelength approaches twice the length of a grid cell, it is not possible any more to represent the wave on the grid (see also section 6.2 of Le Saux et al. 2022).

If we compare the four simulations, $\epsilon$ is larger as the luminosity enhancement factor is increased. This is expected as the amplitude of the waves increase with the boost. Consequently, we should expect non-linear effects to be more relevant in boosted simulations. This could result in more mixing (see e.g. Jermyn 2022), angular momentum transport (see e.g. Gervais et al. 2018) and/or wave-wave interactions.

Finally, by looking at the analytical predictions for $\epsilon$ up to $r = R_{\text{star}}$, we can see that non-linear effects are not expected close to the surface of the star. Most of the waves in the range of frequencies with values of $\epsilon \ge 10^{-3}$ are damped before even reaching $r = 0.91 R_{\text{star}}$, the top of the MUSIC radial domain. For waves that can propagate

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Spatially average temperature profile in the four simulations ref (orange), boost1d1 (indigo), boost1d2 (yellow), and boost1d4 (green). The dashed curves are the temperature profiles obtained at time $t_0$, and the plain curves at $t_0 + \Delta t$, where $\Delta t$ is different for each simulation (see the text for values). The temperature profile from the 1D initial model is represented by the plain blue curve. The outer boundary of the MUSIC numerical domain is indicated by the vertical dashed black line. The right axis corresponds to the term $k_{\text{rad}}N^3$ (blue dotted curve) computed from the 1D model.</figcaption></figure>

further, the maximal values of the non-linear parameter close to the surface are smaller than $10^{-4}$. This is consistent with the results of Ratnasingam et al. (2019) who computed a non-linearity parameter from the analytical spectra of Kumar, Talon & Zahn (1999) and Lecoanet & Quataert (2013) for Reynolds stress excitation of IGWs and found that $\epsilon < 10^{-2}$ in a 3M$_{\odot}$ star for the angular degree $\ell = 10$. However, this does not agree with the results of Horst et al. (2020) who conclude that non-linear effects may be expected at the surface of the star. We suggest that the high values of the non-linear parameter they measure in their simulation is a result of the underestimation of the radiative diffusion. Indeed, as only the luminosity of their model is enhanced, the resulting larger amplitude of IGWs is not compensated by an increase in damping.

Finally, according to the theoretical plots, some waves may reach the surface of the star, but these are relatively high frequency waves compared to the excitation frequency $\omega_{\text{conv}}$. The predictions, however, are based on the linear theory which neglects the interaction with a possible subsurface convection zone (see e.g. Cantiello & Braithwaite 2019) and the near-surface layers that are difficult to accurately model (see e.g. Basu & Chaplin 2017).

# 6 WAVE HEATING IN THE UPPER LAYERS

Often overlooked, thermal effects of IGWs can be significant. In the Earth atmosphere, they are known to irreversibly convert kinetic energy into internal energy (see e.g. Medvedev & Klaassen 2003). Through this mechanism, they can heat up the thermosphere (Yiğit & Medvedev 2009). Almost always neglected in main-sequence stars, these thermal effects of acoustic and internal gravity waves have recently been accounted for in evolved stars in order to explain outbursts in supernova progenitors (Fuller 2017; Wu & Fuller 2022) and in an attempt to explain the lithium enhancement of clump stars (Jermyn & Fuller 2022). In these studies, waves deposit heat through radiative diffusion.

## 6.1 Temperature increase in the simulations

In our simulations we observe a significant increase of the temperature in the upper layers of the model. Fig. 9 presents the evolution of the radial profile of the average temperature in our four
simulations ref (orange), boost1d1 (indigo), boost1d2 (yellow), and boost1d4 (green). The average is performed horizontally, as defined by equation (8). A first average is performed at time $t_0$ (dashed curves) and a second one on $t_0 + \Delta t$ (solid curves), with $\Delta t = 7.0 \times 10^7$, $1.6 \times 10^7$, $1.2 \times 10^7$, and $1.8 \times 10^6$ s for models ref, boost1d1, boost1d2, and boost1d4, respectively. We can clearly see an increase of the temperature close to the top of the domain in the four simulations, from $r \simeq 0.75R_{tot}$ up to $r = 0.91R_{tot}$. The heated region corresponds to the one with strong wave damping observed in Figs 6 and 8. The temperature increase is more important when the enhancement factor of the luminosity is larger. The term $k_{rad}N^3$ (blue dotted curve) computed from the 1D model is also plotted in Fig. 9. In the heated region, this term increases sharply. By looking at equation (15) it is clear that this term drives the radiative damping of a wave at fixed $\omega$ and $\ell$. Therefore, IGWs are strongly damped in this region and will deposit a significant amount of energy. We suggest that the observed heating results from this damping by radiative diffusion.

## 6.2 Theoretical estimate of heat added by waves

In order to test this hypothesis, we analytically estimate the amount of heat $\epsilon_{heat}$ added by waves in model ref through radiative damping. Following Fuller (2017), an IGW is damped by radiative diffusion at a rate $\gamma$ given by

$$ \gamma \approx k_r^2 \kappa_{rad}. \qquad (24) $$

Using equation (24), we can define a damping length for a given wave

$$ l_{damp} = \frac{u_{g,r}}{\gamma}, \qquad (25) $$

with $u_{g,r}$ the radial component of the group velocity of the wave, which can be expressed as (Unno et al. 1989)

$$ u_{g,r} \approx \frac{\omega^2}{N^2 k_h} (N^2 - \omega^2)^{1/2}. \qquad (26) $$

Using equations (24) and (26) and the dispersion relation of IGWs given by equation (12), the damping length expression becomes

$$ l_{damp} = \frac{\omega^3}{N^2 k_h^3 \kappa_{rad}} \left( \frac{N^2}{\omega^2} - 1 \right)^{-1/2}. \qquad (27) $$

Then, the amount of heat deposited per unit mass per unit time by a single wave ($\ell, \omega$) is given by

$$ \epsilon_{heat} = - \frac{dL_{wave}(\ell, \omega)}{dM} = \frac{L_{wave}}{M_{damp}}, \qquad (28) $$

with $L_{wave}$ the wave luminosity and $M_{damp} = 4\pi r^2 \rho l_{damp}$ the mass through which the waves pass before being damped. The wave luminosity in equation (28) corresponds to the luminosity initially injected by convection in a given wave, i.e. it corresponds to the initial amplitude of the wave when it is excited. Finally, radiative damping of IGWs produces a heating by unit time and unit mass for a single wave ($\ell, \omega$) estimated by

$$ \epsilon_{heat}(r, \ell, \omega) = \frac{N^2 k_h^3 \kappa_{rad}}{4\pi r^2 \rho \omega^3} \left( \frac{N^2}{\omega^2} - 1 \right)^{1/2} L_{wave}(r, \ell, \omega) \qquad (29) $$

Using equation (29) we can thus estimate an order of magnitude for the amount of heat theoretically added to the region between $r = 0.75R_{tot}$ and $r = 0.91R_{tot}$. To do this, we measure the wave luminosity $L_{wave}$ for a mode ($\ell, \omega$) in our simulation ref using

$$ L_{wave}(r, \ell, \omega) = L_{wave}(r_e, \ell, \omega) e^{-\tau(r, \ell, \omega)} \qquad (30) $$

with $\tau$ the parameter introduced in equation (5) to take into account radiative damping and

$$ L_{wave}(r_e, \ell, \omega) = 4\pi r^2 F_{wave}(r_e, \ell, \omega) \qquad (31) $$

where $r_e$ is the radius at which waves are excited and $F_{wave}$ has been defined in equation (18), except here we do not assume the low frequency limit. In addition, as we are interested quantitatively in the amplitude of the flux, we multiply $F_{wave}$ by a $3.28/4\pi$ to compensate for the loss of power due to FFT windowing and normalization of the spherical harmonics. We chose $r_e \simeq 0.183R_{star}$ as in Section 3.2 and 5. We consider contribution of waves with frequency $\omega$ in the range [$\omega_{conv}, \omega_{max}$], with $\omega_{max} = 50$ µHz being the maximal frequency available for the simulation. The minimal frequency is set to $\omega_{conv}$ as we do not expect waves with lower frequencies to be excited (Lecanet & Quataert 2013). Similarly, we consider in the estimation the contribution of waves with angular degree $\ell$ from 1 to 200.

Finally, the total amount of heat $Q_{theory}$ added in the region between $r = 0.75R_{tot}$ and $r = 0.91R_{tot}$ during the time $\Delta t$ is

$$ Q_{theory} = \sum_{N_{r_1, r_2}} q_{th}(r), \qquad (32) $$

where $N_{r_1, r_2}$ is the number of radial grid cells in our simulation between $r_1 = 0.75R_{tot}$ and $r_2 = 0.91R_{tot}$ and $q_{th}(r)$ is the amount of heat added in each radial grid cell, it is defined as

$$ q_{th}(r) = \sum_{\ell=0}^{200} \sum_{\omega=\omega_{conv}}^{N_{max}} \rho \epsilon_{heat}(r, \ell, \omega) V_{\text{shell}} \Delta t, \qquad (33) $$

with $V_{shell}$ the volume of the shell between radii $r$ and $r + \Delta r$ where $\Delta r$ is the size of a numerical grid cell in the radial direction. In equation (32), $\Delta t$ is the same time interval as the one used in Fig. 9, which is $\Delta t = 7.0 \times 10^7$ s for model ref. We obtain $Q_{theory} = 3.3 \times 10^{41}$ erg.

Now, in order to compare with the results from the simulations, we measure the total added heat $Q_{MUSIC}$ at the top of the numerical domain using the expression

$$ Q_{MUSIC} = \sum_{N_{r_1, r_2}} \rho c_p \Delta T V_{\text{shell}} \qquad (34) $$

with $c_p$ the specific heat capacity at constant pressure and $\Delta T$ the temperature difference after the same interval of time $\Delta t$. The total added heat in model ref is $Q_{MUSIC} = 1.0 \times 10^{43}$ erg. This value is larger than $Q_{theory}$, but this can be expected. In the simulations, the top boundary $r_{out} = 0.91R_{star}$ is reflecting IGWs, therefore waves that would be damped in the region between $r_{out}$ and the surface in an actual star are damped in the region between $r_1 = 0.75R_{tot}$ and $r_{out}$ in our simulation. As previously mentioned, this heating could be strengthened due to boundary conditions. The issue of the impact of boundary conditions on hydrodynamical simulations is an open challenge that affects all stellar simulations (see e.g. Vlakov et al. 2022). We thus consider that the value of $Q_{MUSIC}$ is in relatively good agreement with the value estimated from linear theory $Q_{theory}$. This strengthens our confidence that the heating observed in the simulations results from the damping of IGWs by radiative diffusion.

The estimation of $Q_{theory}$ is performed between $r_1 = 0.75R_{star}$ and $r_2 = 0.91R_{star}$. However, the location of the inner boundary of this domain has a limited influence on the computation of $Q_{theory}$. This is highlighted in Fig. 10. This figure presents the radial cumulative sum of added heat by IGWs, $S_{heat}$, estimated with the analytical expression
simulations ref (orange), boost1d1 (indigo), boost1d2 (yellow), and boost1d4 (green). The average is performed horizontally, as defined by equation (8). A first average is performed at time $t_0$ (dashed curves) and a second one on $t_0 + \Delta t$ (solid curves), with $\Delta t = 7.0 \times 10^7$, $1.6 \times 10^7$, $1.2 \times 10^7$, and $1.8 \times 10^6$ s for models ref, boost1d1, boost1d2, and boost1d4, respectively. We can clearly see an increase of the temperature close to the top of the domain in the four simulations, from $r \simeq 0.75R_{tot}$ up to $r = 0.91R_{tot}$. The heated region corresponds to the one with strong wave damping observed in Figs 6 and 8. The temperature increase is more important when the enhancement factor of the luminosity is larger. The term $k_{rad}N^3$ (blue dotted curve) computed from the 1D model is also plotted in Fig. 9. In the heated region, this term increases sharply. By looking at equation (15) it is clear that this term drives the radiative damping of a wave at fixed $\omega$ and $\ell$. Therefore, IGWs are strongly damped in this region and will deposit a significant amount of energy. We suggest that the observed heating results from this damping by radiative diffusion.

## 6.2 Theoretical estimate of heat added by waves

In order to test this hypothesis, we analytically estimate the amount of heat $\epsilon_{heat}$ added by waves in model ref through radiative damping. Following Fuller (2017), an IGW is damped by radiative diffusion at a rate $\gamma$ given by

$$ \gamma \approx k_r^2 \kappa_{rad}. \quad (24) $$

Using equation (24), we can define a damping length for a given wave

$$ l_{damp} = \frac{u_{g,r}}{\gamma}, \quad (25) $$

with $u_{g,r}$ the radial component of the group velocity of the wave, which can be expressed as (Unno et al. 1989)

$$ u_{g,r} \approx \frac{\omega^2}{N^2 k_h} (N^2 - \omega^2)^{1/2}. \quad (26) $$

Using equations (24) and (26) and the dispersion relation of IGWs given by equation (12), the damping length expression becomes

$$ l_{damp} = \frac{\omega^3}{N^2 k_h^3 \kappa_{rad}} \left( \frac{N^2}{\omega^2} - 1 \right)^{-1/2}. \quad (27) $$

Then, the amount of heat deposited per unit mass per unit time by a single wave $(\ell, \omega)$ is given by

$$ \epsilon_{heat} = -\frac{dL_{wave}(\ell, \omega)}{dM} = \frac{L_{wave}}{M_{damp}}, \quad (28) $$

with $L_{wave}$ the wave luminosity and $M_{damp} = 4\pi r^2 \rho l_{damp}$ the mass through which the waves pass before being damped. The wave luminosity in equation (28) corresponds to the luminosity initially injected by convection in a given wave, i.e. it corresponds to the initial amplitude of the wave when it is excited. Finally, radiative damping of IGWs produces a heating by unit time and unit mass for a single wave $(\ell, \omega)$ estimated by

$$ \epsilon_{heat}(r, \ell, \omega) = \frac{N^2 k_h^3 \kappa_{rad}}{4\pi r^2 \rho \omega^3} \left( \frac{N^2}{\omega^2} - 1 \right)^{1/2} L_{wave}(r, \ell, \omega) \quad (29) $$

Using equation (29) we can thus estimate an order of magnitude for the amount of heat theoretically added to the region between $r = 0.75R_{tot}$ and $r = 0.91R_{tot}$. To do this, we measure the wave luminosity $L_{wave}$ for a mode $(\ell, \omega)$ in our simulation ref using

$$ L_{wave}(r, \ell, \omega) = L_{wave}(r_e, \ell, \omega) e^{-\tau(r, \ell, \omega)} \quad (30) $$

with $\tau$ the parameter introduced in equation (5) to take into account radiative damping and

$$ L_{wave}(r_e, \ell, \omega) = 4\pi r^2 F_{wave}(r_e, \ell, \omega) \quad (31) $$

where $r_e$ is the radius at which waves are excited and $F_{wave}$ has been defined in equation (18), except here we do not assume the low frequency limit. In addition, as we are interested quantitatively in the amplitude of the flux, we multiply $F_{wave}$ by a $3.28/4\pi$ to compensate for the loss of power due to FFT windowing and normalization of the spherical harmonics. We chose $r_e \simeq 0.183R_{star}$ as in Section 3.2 and 5. We consider contribution of waves with frequency $\omega$ in the range [$\omega_{conv}$; $\omega_{max}$], with $\omega_{max} = 50$ µHz being the maximal frequency available for the simulation. The minimal frequency is set to $\omega_{conv}$ as we do not expect waves with lower frequencies to be excited (Lecanet & Quataert 2013). Similarly, we consider in the estimation the contribution of waves with angular degree $\ell$ from 1 to 200.

Finally, the total amount of heat $Q_{theory}$ added in the region between $r = 0.75R_{tot}$ and $r = 0.91R_{tot}$ during the time $\Delta t$ is

$$ Q_{theory} = \sum_{N_{r_1, r_2}} q_{th}(r), \quad (32) $$

where $N_{r_1, r_2}$ is the number of radial grid cells in our simulation between $r_1 = 0.75R_{tot}$ and $r_2 = 0.91R_{tot}$ and $q_{th}(r)$ is the amount of heat added in each radial grid cell, it is defined as

$$ q_{th}(r) = \sum_{\ell=0}^{200} \sum_{\omega=\omega_{conv}}^{N_{max}} \rho \epsilon_{heat}(r, \ell, \omega) V_{\text{shell}} \Delta t, \quad (33) $$

with $V_{shell}$ the volume of the shell between radii $r$ and $r + \Delta r$ where $\Delta r$ is the size of a numerical grid cell in the radial direction. In equation (32), $\Delta t$ is the same time interval as the one used in Fig. 9, which is $\Delta t = 7.0 \times 10^7$ s for model ref. We obtain $Q_{theory} = 3.3 \times 10^{41}$ erg.

Now, in order to compare with the results from the simulations, we measure the total added heat $Q_{MUSIC}$ at the top of the numerical domain using the expression

$$ Q_{MUSIC} = \sum_{N_{r_1, r_2}} \rho c_p \Delta T V_{\text{shell}} \quad (34) $$

with $c_p$ the specific heat capacity at constant pressure and $\Delta T$ the temperature difference after the same interval of time $\Delta t$. The total added heat in model ref is $Q_{MUSIC} = 1.0 \times 10^{43}$ erg. This value is larger than $Q_{theory}$, but this can be expected. In the simulations, the top boundary $r_{out} = 0.91R_{star}$ is reflecting IGWs, therefore waves that would be damped in the region between $r_{out}$ and the surface in an actual star are damped in the region between $r_1 = 0.75R_{tot}$ and $r_{out}$ in our simulation. As previously mentioned, this heating could be strengthened due to boundary conditions. The issue of the impact of boundary conditions on hydrodynamical simulations is an open challenge that affects all stellar simulations (see e.g. Vlakov et al. 2022). We thus consider that the value of $Q_{MUSIC}$ is in relatively good agreement with the value estimated from linear theory $Q_{theory}$. This strengthens our confidence that the heating observed in the simulations results from the damping of IGWs by radiative diffusion.

The estimation of $Q_{theory}$ is performed between $r_1 = 0.75R_{star}$ and $r_2 = 0.91R_{star}$. However, the location of the inner boundary of this domain has a limited influence on the computation of $Q_{theory}$. This is highlighted in Fig. 10. This figure presents the radial cumulative sum of added heat by IGWs, $S_{heat}$, estimated with the analytical expression
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Cumulative sum of the theoretical estimate of the heat deposited by waves, $S_{\text{heat}}$, as a function of normalized radius. The horizontal dashed line indicates the value $S_{\text{heat}} = 1.0$.</figcaption></figure>

given by equation (32). It is defined as

$$S_{\text{heat}}(r) = \frac{\sum_{N_{r_{0,}}} q_{\text{th}}(r)}{Q_{\text{theory}}} \qquad (35)$$

where $N_{r_{0,}}$ is the number of radial grid cells between $r_0 = 0.2R_{\text{star}}$ and a given radius $r$. We can see that most of the heat that waves can deposit is at radii $r \ge 0.7R_{\text{star}}$. This corresponds to the region where we observe an increase of the temperature in Fig. 9 for model ref.

In summary, our results show an increase of the temperature in the upper layers of the star. The amount of heating is in agreement with theoretical predictions of heat added by IGWs damped by radiative diffusion. This thermal impact of IGWs may thus be relevant in main-sequence stars and is currently under investigation.

## 7 DISCUSSION AND CONCLUSION

This study presents an analysis of the properties of IGWs excited by convection in a 2D fully compressible simulation of a 5 solar mass star model at ZAMS. Our reference simulation is run with a luminosity that is not artificially enhanced and with a realistic radiative diffusivity profile. The simulation radial domain extends from $r_{\text{in}} = 0.02R_{\text{star}}$ to $r_{\text{out}} = 0.91R_{\text{star}}$. We have highlighted that waves propagating in the radiative envelope of such stars are strongly damped by radiative diffusion. In our truncated model, the radiative diffusivity varies radially by 5 orders of magnitude. In the region not modelled in this work, between $r_{\text{out}}$ and $R_{\text{star}}$, we expect the effect of radiative damping to be even more important, since the radiative diffusivity varies by more than 4 orders of magnitude in this region. Our analysis highlights the importance of including radiative diffusion in stellar hydrodynamical simulations. This is particularly important if the goal is to analyse the waves that can reach the stellar surface and to establish a link with observations or to study their transport properties. A limitation of our work is that we do not include rotation in our simulations, since the primary goal is to analyse the effect of radiative damping with realistic radiative diffusivity profiles. A relevant comparison of simulated spectra with observations is a challenge as it would require a proper description of the surface layers and including rotation. Indeed, many OB stars observed showing photometric variability are rotating (e.g. Pedersen et al. 2021; Szewczuk, Walczak & Daszyńska-Daszkiewicz 2021). In addition to the dynamical effects of rotation on convection

(e.g. impacting $\omega_{\text{conv}}$ and the plume dynamics) and thus on wave excitation, rotation may produce a shift in the wave frequencies. As already mentioned, Rogers et al. (2013) suggest a large impact of differential rotation on the low frequency power spectrum. But these conclusions are based on an imposed differential rotation profile and would need confirmation with further numerical simulations. Interestingly, the shape and structure of the power excess observed in these stars (e.g. Bowman et al. 2019; Szewczuk et al. 2021) is very similar for most stars in their sample. Since the rotation rates of B dwarfs vary by more than 1 order of magnitude, it is compelling to observe such similarity, since the expected frequency shift would depend on the rotation rate.

In the simulations presented in this work, we observe an increase of the temperature close to the top of the numerical domain. In this region the damping of IGWs is strong, due to the simultaneous increase of stratification and radiative diffusivity. In the geophysics community, it is well known that IGWs can heat up the atmosphere. Using the linear theory of IGWs, we have estimated an order of magnitude for the amount of heat that could be added by waves in these upper layers. This value is comparable to the value inferred from the MUSIC simulations, suggesting that IGWs may be at the origin of the observed heating. Despite the fact that this heating could be slightly enhanced by the reflection of IGWs on the top boundary of our numerical domain, our results show that it is not an artefact. Since the radiative diffusivity can vary steeply with radius, wave damping is not uniform throughout the star and waves will deposit different amounts of energy at different radii. Once again, this highlights the importance of using a realistic radiative diffusivity profile. However, this heating induced by waves is, to our knowledge, always neglected in main-sequence stars. This may not be always justified and thermal effects of IGWs may be relevant in some cases, particularly as this occurs in the outer part of the star. In a future work, we will study in which context this has to be taken into account and how it could impact stellar structures.

We have also studied non-linear effects linked to IGWs propagating in the radiative envelope of intermediate-mass stars. Waves can be efficiently excited by turbulent convection, with a large amplitude, and as they travel towards the surface IGWs see their amplitude growing as the density of the medium is decreasing. Their amplitude is also decaying due to damping by radiative diffusion, and it turns out that this effect is dominant for large amplitude IGWs, i.e. those generated with frequency close to the convective frequency. Our results suggest that non-linear effects may be relevant only above the convective core, close to where waves are excited. In this region, waves with frequencies close to the convective frequency appear to be strongly damped or to break in MUSIC simulations. More precisely, these waves have a non-linear parameter $\epsilon$ close to 1, meaning that they may be highly non-linear. A possible explanation could be that these waves are generated with a very large amplitude and that they likely break almost immediately due to non-linear effects, just above the convective core. This could result in mode coupling and/or generation of lower frequency waves. However, we cannot exclude that the signal in this very low frequency range results from aliasing as the associated wavelengths may not be resolved properly on the numerical grid of our model. This is a common issue in hydrodynamical simulations that is difficult to quantify.

The results obtained for our reference model ref are compared with three simulations for which the luminosity has been increased by factors 10, $10^2$ and $10^4$. This comparison highlights the impact of this artefact on the generation and propagation of IGWs. In stars, the main IGWs excitation mechanism is linked to convection and has two
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Cumulative sum of the theoretical estimate of the heat deposited by waves, $S_{\text{heat}}$, as a function of normalized radius. The horizontal dashed line indicates the value $S_{\text{heat}} = 1.0$.</figcaption></figure>

given by equation (32). It is defined as

$$S_{\text{heat}}(r) = \frac{\sum_{N_{r_{0,r}}} q_{\text{th}}(r)}{Q_{\text{theory}}} \qquad (35)$$

where $N_{r_{0,r}}$ is the number of radial grid cells between $r_0 = 0.2R_{\text{star}}$ and a given radius $r$. We can see that most of the heat that waves can deposit is at radii $r \ge 0.7R_{\text{star}}$. This corresponds to the region where we observe an increase of the temperature in Fig. 9 for model ref.

In summary, our results show an increase of the temperature in the upper layers of the star. The amount of heating is in agreement with theoretical predictions of heat added by IGWs damped by radiative diffusion. This thermal impact of IGWs may thus be relevant in main-sequence stars and is currently under investigation.

## 7 DISCUSSION AND CONCLUSION

This study presents an analysis of the properties of IGWs excited by convection in a 2D fully compressible simulation of a 5 solar mass star model at ZAMS. Our reference simulation is run with a luminosity that is not artificially enhanced and with a realistic radiative diffusivity profile. The simulation radial domain extends from $r_{\text{in}} = 0.02R_{\text{star}}$ to $r_{\text{out}} = 0.91R_{\text{star}}$. We have highlighted that waves propagating in the radiative envelope of such stars are strongly damped by radiative diffusion. In our truncated model, the radiative diffusivity varies radially by 5 orders of magnitude. In the region not modelled in this work, between $r_{\text{out}}$ and $R_{\text{star}}$, we expect the effect of radiative damping to be even more important, since the radiative diffusivity varies by more than 4 orders of magnitude in this region. Our analysis highlights the importance of including radiative diffusion in stellar hydrodynamical simulations. This is particularly important if the goal is to analyse the waves that can reach the stellar surface and to establish a link with observations or to study their transport properties. A limitation of our work is that we do not include rotation in our simulations, since the primary goal is to analyse the effect of radiative damping with realistic radiative diffusivity profiles. A relevant comparison of simulated spectra with observations is a challenge as it would require a proper description of the surface layers and including rotation. Indeed, many OB stars observed showing photometric variability are rotating (e.g. Pedersen et al. 2021; Szewczuk, Walczak & Daszyńska-Daszkiewicz 2021). In addition to the dynamical effects of rotation on convection

(e.g. impacting $\omega_{\text{conv}}$ and the plume dynamics) and thus on wave excitation, rotation may produce a shift in the wave frequencies. As already mentioned, Rogers et al. (2013) suggest a large impact of differential rotation on the low frequency power spectrum. But these conclusions are based on an imposed differential rotation profile and would need confirmation with further numerical simulations. Interestingly, the shape and structure of the power excess observed in these stars (e.g. Bowman et al. 2019; Szewczuk et al. 2021) is very similar for most stars in their sample. Since the rotation rates of B dwarfs vary by more than 1 order of magnitude, it is compelling to observe such similarity, since the expected frequency shift would depend on the rotation rate.

In the simulations presented in this work, we observe an increase of the temperature close to the top of the numerical domain. In this region the damping of IGWs is strong, due to the simultaneous increase of stratification and radiative diffusivity. In the geophysics community, it is well known that IGWs can heat up the atmosphere. Using the linear theory of IGWs, we have estimated an order of magnitude for the amount of heat that could be added by waves in these upper layers. This value is comparable to the value inferred from the MUSIC simulations, suggesting that IGWs may be at the origin of the observed heating. Despite the fact that this heating could be slightly enhanced by the reflection of IGWs on the top boundary of our numerical domain, our results show that it is not an artefact. Since the radiative diffusivity can vary steeply with radius, wave damping is not uniform throughout the star and waves will deposit different amounts of energy at different radii. Once again, this highlights the importance of using a realistic radiative diffusivity profile. However, this heating induced by waves is, to our knowledge, always neglected in main-sequence stars. This may not be always justified and thermal effects of IGWs may be relevant in some cases, particularly as this occurs in the outer part of the star. In a future work, we will study in which context this has to be taken into account and how it could impact stellar structures.

We have also studied non-linear effects linked to IGWs propagating in the radiative envelope of intermediate-mass stars. Waves can be efficiently excited by turbulent convection, with a large amplitude, and as they travel towards the surface IGWs see their amplitude growing as the density of the medium is decreasing. Their amplitude is also decaying due to damping by radiative diffusion, and it turns out that this effect is dominant for large amplitude IGWs, i.e. those generated with frequency close to the convective frequency. Our results suggest that non-linear effects may be relevant only above the convective core, close to where waves are excited. In this region, waves with frequencies close to the convective frequency appear to be strongly damped or to break in MUSIC simulations. More precisely, these waves have a non-linear parameter $\epsilon$ close to 1, meaning that they may be highly non-linear. A possible explanation could be that these waves are generated with a very large amplitude and that they likely break almost immediately due to non-linear effects, just above the convective core. This could result in mode coupling and/or generation of lower frequency waves. However, we cannot exclude that the signal in this very low frequency range results from aliasing as the associated wavelengths may not be resolved properly on the numerical grid of our model. This is a common issue in hydrodynamical simulations that is difficult to quantify.

The results obtained for our reference model ref are compared with three simulations for which the luminosity has been increased by factors 10, $10^2$ and $10^4$. This comparison highlights the impact of this artefact on the generation and propagation of IGWs. In stars, the main IGWs excitation mechanism is linked to convection and has two
components that are Reynolds stress and penetrative convection, both acting simultaneously. Both these mechanisms generate waves in the same frequency range. It is however difficult to determine whether one mechanism is more efficient than the other. In our simulation with realistic luminosity and radiative diffusivity, the wave flux measured is broadly consistent with excitation by penetrative convection at low frequency and with Reynolds stress at higher frequencies. This is similar in the two less boosted models but not in the most boosted one. In the latter, the wave energy flux is approximately consistent with an excitation by penetrative convection but not by Reynolds stress. This suggests that the efficiency of these two excitation mechanisms is also impacted by an enhancement of the luminosity and radiative diffusivity. We have also shown that the heating close to the top of the domain and non-linear effects are more important when the luminosity is artificially increased. In addition, this should also impact the transport of angular momentum, which depends on the radiative damping of IGWs and on their amplitude. Results are thus quantitatively modified by an enhancement of the luminosity and any prediction based on this artefact should be taken with caution.

Concerning the surface manifestation of IGWs excited by core convection, our main conclusion is that extrapolating simulated spectra determined at an internal radius to the surface in order to compare to observations is likely meaningless. Such a comparison would require numerical simulations extending up to the surface layers, to properly describe the radiative damping in these layers, which is a formidable challenge for stellar hydrodynamics simulations. Finally, we can identify two additional challenges regarding the comparison between observations and numerical simulations. First, simulations require physical simplifications and are thus far from realistic stellar conditions. As shown in this work, these assumptions can impact the physics of waves. Second, observations do not resolve stellar surface and thus only consider global variations of luminosity. Further efforts are thus needed in order to improve the reliability of comparisons between observations and the predictions of numerical simulations.

## ACKNOWLEDGEMENTS

We thank Jim Fuller and Adam Jermyn for interesting discussion about wave-induced heating. This work is partly supported by the consolidated STFC grant ST/R000395/1 and the ERC grant No. 787361-COBOM. 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 DE-AC52-07NA27344.

## DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

## REFERENCES

* Aerts C., Rogers T. M., 2015, ApJ, 806, L33

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

* Aerts C., Augustson K., Mathis S., Pedersen M. G., Mombarg J. S. G., Vanlaer V., Van Beeck J., Van Reeth T., 2021, A&A, 656, A121

* Alvan L., Brun A. S., Mathis S., 2014, A&A, 565, A42

* Appourchaux T. et al., 2010, A&AR, 18, 197

* Auvergne M. et al., 2009, A&A, 506, 411

* 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

* Basu S., Chaplin W. J., 2017, Asteroseismic Data Analysis: Foundations and Techniques Princeton University Press Princeton

* Blomme R. et al., 2011, A&A, 533, A4

* Borucki W. J. et al., 2010, Science, 327, 977

* Bowman D. M. et al., 2019, Nat. Astron., 3, 760

* 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

* Cantiello M., Braithwaite J., 2019, ApJ, 883, 106

* Cantiello M., Lecoanet D., Jermyn A. S., Grassitelli L., 2021, ApJ, 915, 112

* Charbonnel C., Talon S., 2005, Science, 309, 2189

* De Cat P., Aerts C., 2002, A&A, 393, 965

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

* Fuller J., 2017, MNRAS, 470, 1642

* García R. A., Turck-Chièze S., Jiménez-Reyes S. J., Ballot J., Pallé P. L., Eff-Darwich A., Mathur S., Provost J., 2007, Science, 316, 1591

* Geroux C. et al., 2016, A&A, 588, A85

* Gervais A. D., Swaters G. E., van den Bremer T. S., Sutherland B. R., 2018, J. Atmos. Sci., 75, 3703

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

* Goldstein J., Townsend R. H. D., 2020, ApJ, 899, 116

* Horst L., Edelmann P. V. F., Andrássy R., Röpke F. K., Bowman D. M., Aerts C., Ratnasingam R. P., 2020, A&A, 641, A18

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

* Jermyn A. S., 2022, preprint (arXiv:2209.08344)

* Jermyn A. S., Fuller J., 2022, preprint (arXiv:2206.13479)

* Kumar P., Talon S., Zahn J.-P., 1999, ApJ, 520, 859

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

* 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

* Lighthill M. J., 1952, Proc. R. Soc. A, 211, 564

* Lighthill J., 1978, Waves in Fluids. Cambridge Univ. Press, Cambridge

* Medvedev A. S., Klaassen G. P., 2003, J. Geophys. Res.: Atmos., 108, 4

* Pedersen M. G. et al., 2021, Nat. Astron., 5, 715

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

* Pinçon C., Belkacem K., Goupil M. J., Marques J. P., 2017, A&A, 605, A31

* Pratt J. et al., 2016, A&A, 593, A121

* 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., Edelmann P. V. F., Rogers T. M., 2019, MNRAS, 482, 5500

* Ratnasingam R. P., Edelmann P. V. F., Rogers T. M., 2020, MNRAS, 497, 4231

* Renzini A., 1987, A&A, 188, 49

* Ricker G. R. et al., 2015, J. Astron. Telescopes Instrum. Syst., 1, 014003

* Rogers T. M., Glatzmaier G. A., 2005, MNRAS, 364, 1135

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

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

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

* Stein R. F., 1967, Sol. Phys., 2, 385
components that are Reynolds stress and penetrative convection, both acting simultaneously. Both these mechanisms generate waves in the same frequency range. It is however difficult to determine whether one mechanism is more efficient than the other. In our simulation with realistic luminosity and radiative diffusivity, the wave flux measured is broadly consistent with excitation by penetrative convection at low frequency and with Reynolds stress at higher frequencies. This is similar in the two less boosted models but not in the most boosted one. In the latter, the wave energy flux is approximately consistent with an excitation by penetrative convection but not by Reynolds stress. This suggests that the efficiency of these two excitation mechanisms is also impacted by an enhancement of the luminosity and radiative diffusivity. We have also shown that the heating close to the top of the domain and non-linear effects are more important when the luminosity is artificially increased. In addition, this should also impact the transport of angular momentum, which depends on the radiative damping of IGWs and on their amplitude. Results are thus quantitatively modified by an enhancement of the luminosity and any prediction based on this artefact should be taken with caution.

Concerning the surface manifestation of IGWs excited by core convection, our main conclusion is that extrapolating simulated spectra determined at an internal radius to the surface in order to compare to observations is likely meaningless. Such a comparison would require numerical simulations extending up to the surface layers, to properly describe the radiative damping in these layers, which is a formidable challenge for stellar hydrodynamics simulations. Finally, we can identify two additional challenges regarding the comparison between observations and numerical simulations. First, simulations require physical simplifications and are thus far from realistic stellar conditions. As shown in this work, these assumptions can impact the physics of waves. Second, observations do not resolve stellar surface and thus only consider global variations of luminosity. Further efforts are thus needed in order to improve the reliability of comparisons between observations and the predictions of numerical simulations.

## ACKNOWLEDGEMENTS

We thank Jim Fuller and Adam Jermyn for interesting discussion about wave-induced heating. This work is partly supported by the consolidated STFC grant ST/R000395/1 and the ERC grant No. 787361-COBOM. 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 DE-AC52-07NA27344.

## DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

## REFERENCES

* Aerts C., Rogers T. M., 2015, ApJ, 806, L33

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

* Aerts C., Augustson K., Mathis S., Pedersen M. G., Mombarg J. S. G., Vanlaer V., Van Beeck J., Van Reeth T., 2021, A&A, 656, A121

* Alvan L., Brun A. S., Mathis S., 2014, A&A, 565, A42

* Appourchaux T. et al., 2010, A&AR, 18, 197

* Auvergne M. et al., 2009, A&A, 506, 411

* 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

* Basu S., Chaplin W. J., 2017, Asteroseismic Data Analysis: Foundations and Techniques Princeton University Press Princeton

* Blomme R. et al., 2011, A&A, 533, A4

* Borucki W. J. et al., 2010, Science, 327, 977

* Bowman D. M. et al., 2019, Nat. Astron., 3, 760

* 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

* Cantiello M., Braithwaite J., 2019, ApJ, 883, 106

* Cantiello M., Lecoanet D., Jermyn A. S., Grassitelli L., 2021, ApJ, 915, 112

* Charbonnel C., Talon S., 2005, Science, 309, 2189

* De Cat P., Aerts C., 2002, A&A, 393, 965

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

* Fuller J., 2017, MNRAS, 470, 1642

* García R. A., Turck-Chièze S., Jiménez-Reyes S. J., Ballot J., Pallé P. L., Eff-Darwich A., Mathur S., Provost J., 2007, Science, 316, 1591

* Geroux C. et al., 2016, A&A, 588, A85

* Gervais A. D., Swaters G. E., van den Bremer T. S., Sutherland B. R., 2018, J. Atmos. Sci., 75, 3703

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

* Goldstein J., Townsend R. H. D., 2020, ApJ, 899, 116

* Horst L., Edelmann P. V. F., Andrássy R., Röpke F. K., Bowman D. M., Aerts C., Ratnasingam R. P., 2020, A&A, 641, A18

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

* Jermyn A. S., 2022, preprint (arXiv:2209.08344)

* Jermyn A. S., Fuller J., 2022, preprint (arXiv:2206.13479)

* Kumar P., Talon S., Zahn J.-P., 1999, ApJ, 520, 859

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

* 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

* Lighthill M. J., 1952, Proc. R. Soc. A, 211, 564

* Lighthill J., 1978, Waves in Fluids. Cambridge Univ. Press, Cambridge

* Medvedev A. S., Klaassen G. P., 2003, J. Geophys. Res.: Atmos., 108, 4

* Pedersen M. G. et al., 2021, Nat. Astron., 5, 715

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

* Pinçon C., Belkacem K., Goupil M. J., Marques J. P., 2017, A&A, 605, A31

* Pratt J. et al., 2016, A&A, 593, A121

* 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., Edelmann P. V. F., Rogers T. M., 2019, MNRAS, 482, 5500

* Ratnasingam R. P., Edelmann P. V. F., Rogers T. M., 2020, MNRAS, 497, 4231

* Renzini A., 1987, A&A, 188, 49

* Ricker G. R. et al., 2015, J. Astron. Telescopes Instrum. Syst., 1, 014003

* Rogers T. M., Glatzmaier G. A., 2005, MNRAS, 364, 1135

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

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

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

* Stein R. F., 1967, Sol. Phys., 2, 385
Sutherland B. R., 2010, Internal Gravity Waves. Cambridge Univ. Press, Cambridge

Szewczuk W., Walczak P., Daszyńska-Daszkiewicz J., 2021, MNRAS, 503, 5894

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

Townsend R. H. D., Goldstein J., Zweibel E. G., 2018, MNRAS, 475, 879

Unno W., Osaki Y., Ando H., Saio H., Shibahashi H., 1989, Nonradial Oscillations of Stars, 2nd edn. University of Tokyo Press, Tokyo

Vallis G. K., 2017, Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation. Cambridge Univ. Press, Cambridge

Viallet M., Baraffe I., Walder R., 2011, A&A, 531, A86

Viallet M., Baraffe I., Walder R., 2013, A&A, 555, A81

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

Vlaykov D. G., Baraffe I., Constantino T., Goffrey T., Guillet T., Le Saux A., Morison A., Pratt J., 2022, MNRAS, 514, 715

Wu S. C., Fuller J., 2022, ApJ, 930, 119

Yiğit E., Medvedev A. S., 2009, Geophys. Res. Lett., 36

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.
Sutherland B. R., 2010, Internal Gravity Waves. Cambridge Univ. Press, Cambridge

Szewczuk W., Walczak P., Daszyńska-Daszkiewicz J., 2021, MNRAS, 503, 5894

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

Townsend R. H. D., Goldstein J., Zweibel E. G., 2018, MNRAS, 475, 879

Unno W., Osaki Y., Ando H., Saio H., Shibahashi H., 1989, Nonradial Oscillations of Stars, 2nd edn. University of Tokyo Press, Tokyo

Vallis G. K., 2017, Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation. Cambridge Univ. Press, Cambridge

Viallet M., Baraffe I., Walder R., 2011, A&A, 531, A86

Viallet M., Baraffe I., Walder R., 2013, A&A, 555, A81

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

Vlaykov D. G., Baraffe I., Constantino T., Goffrey T., Guillet T., Le Saux A., Morison A., Pratt J., 2022, MNRAS, 514, 715

Wu S. C., Fuller J., 2022, ApJ, 930, 119

Yiğit E., Medvedev A. S., 2009, Geophys. Res. Lett., 36

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.