Three-dimensional Simulations of Massive Stars. II. Age Dependence

As mentioned in Section 1, theoretical studies on IGW propagation within stars assume that their generation spectra are driven exclusively by either overshooting plume incursions (Townsend 1966; Montalbán & Schatzman 2000; Pinçon et al. 2016) or convective eddies (Lighthill 1952; Goldreich & Kumar 1990; Kumar et al. 1999; Lecoanet & Quataert 2013). The more recent analysis by Le Saux et al. (2023) notes that the dominance of plumes or eddies may be related to the convective forcing. Of course, it is highly likely that both plumes and convective eddies play a part in generating IGWs. Given that the spectra of IGWs throughout the star are sensitive to their generation mechanism, it is important to understand the relevance of both of these mechanisms within our simulations. For this reason, we analyze frequency spectra for the kinetic energy:

$$\begin{eqnarray}&&{\hat{E}}_{\mathrm{kin}}=\displaystyle \frac{1}{2}\overline{\rho }\left({\hat{v}}_{r}^{2}+{\hat{v}}_{\theta }^{2}+{\hat{v}}_{\phi }^{2}\right),\end{eqnarray} \tag{ 12 }$$

where ${\hat{E}}_{\mathrm{kin}}$ represents the Fourier transform of the real quantity E_kin.

Figure 4 shows the frequency spectra for the kinetic energy ${\hat{E}}_{\mathrm{kin}}$ taken 0.07H_p below (Figure 4(a)) and 0.07H_p above (Figure 4(b)) the CB ⁶ (the latter radial location is chosen to match the value used for the generation spectra analysis by Edelmann et al. 2019, which has the added advantage of being located within the BVF spike in midMS and TAMS, allowing us to sample all the waves generated) for all three stellar ages in panel (i). In panel (ii) the same kinetic energy spectrum has been multiplied by the frequency f to account for the integration over $d\mathrm{log}f;$ this helps to emphasize the regions where most of the kinetic energy is contained. In both panels we also plot color-coded vertical lines corresponding to the estimates of the turnover frequencies of the models according to

$$\begin{eqnarray}&&{f}_{\mathrm{TO}}=\displaystyle \frac{{v}_{\mathrm{rms}}^{\mathrm{CZ}}}{\pi {r}_{\mathrm{CZ}}},\end{eqnarray} \tag{ 13 }$$

where the expression assumes that the largest eddy spans the entire radial extent of the convective zone r_CZ and turns at ${v}_{\mathrm{rms}}^{\mathrm{CZ}};$ for convenience we also note here the actual values of the turnover frequencies used, which are ${f}_{\mathrm{TO}}^{{\rm{Z}}}=2.51\,\mu \mathrm{Hz}$, ${f}_{\mathrm{TO}}^{{\rm{m}}}=3.45\,\mu \mathrm{Hz}$, and ${f}_{\mathrm{TO}}^{{\rm{T}}}=3.93\,\mu \mathrm{Hz}$. In Figure 4 the peaks of the spectra in panel (ii) are located at frequencies roughly one order of magnitude higher than the estimated turnover frequencies; furthermore, the spectra do not appear dominated by one single frequency, with most of the kinetic energy of the motions being stored in frequencies within the range 10¹ ≲ f ≲ 10² μHz. This is in contrast with the theoretical idea of a generation spectrum that is mostly influenced by convective eddies, where most of the energy is concentrated at the convective turnover frequency with steep drops away from it. While the lower end of the spectrum within the radiative zone could be affected by viscous and thermal diffusivities, the fact that we do not see a prominent convective frequency within the convection zone indicates that this concept may not be relevant in this geometry.

Zoom In Zoom Out Reset image size

The spectra shown in Figure 4(b), as mentioned previously, are taken just above (0.07H_p ) the CB rather than below it; in the midMS and TAMS models this approximately coincides with the maximum value reached by N in their BVF spikes. Here the spectra, particularly once they are multiplied by f in panel (ii), appear much flatter and it is impossible to pinpoint a dominant frequency (particularly for the ZAMS model) or the location of the turnover frequencies.

While both of these sets of spectra represent the motions generated at the CB, rather than specifically waves, these motions are responsible for driving the waves and the properties of the two are therefore closely linked; given the spectra of the motions on either side of the CB are not dominated by the turnover frequency, it is safe to assume that the IGW energy will likewise not be concentrated at f_TO.

Convective plumes overshooting into the radiative zone play a significant role in the dynamical evolution of stars. Apart from generating IGWs at the CB, as previously mentioned, they also lead to the formation of a layer between the convective and radiative zones where turbulent motions intermittently mix with the radiative flow. The properties of this layer can influence chemical mixing, the generation of the magnetic field, and the evolutionary track in the Hertzsprung–Russell diagram of a star, therefore playing a role in its evolution (Marik & Petrovay 2002; Marques et al. 2006; Tian et al. 2009).

While the physics of this overshooting layer lie beyond the scope of this paper, we aim to understand how plume overshooting varies with age across our models. For this purpose, we carry out an analysis on their incursion into the radiation zone; we scan every output file (which are separated by an interval of 1000 s) in every (ϕ, θ) direction to look for possible plumes overshooting beyond the CB. This is accomplished by monitoring the sign of the radial velocity v_r or, equivalently, of the vertical kinetic energy flux ${E}_{K}=\tfrac{1}{2}{v}_{r}\rho {{\boldsymbol{v}}}^{2}$ (Hurlburt et al. 1986, 1994), whose sign changes when plumes are decelerated. In particular, we define the overshooting length Δr_ov of a plume at a given (ϕ, θ) combination as the distance from the CB at which v_r turns from positive to negative. The results are subsequently refined to remove false positives by filtering all the obtained overshooting lengths by a minimum depth (i.e., the radial extent of a few grid cells, meaning any smaller depths are ignored) and by finding structures which are coherent both spatially and temporally among the filtered data. Each coherent structure is labeled as a plume, and its maximum incursion depth and velocity are recorded.

A similar analysis has been carried out by Pratt et al. (2016, 2017), who employed 2D simulations of spherical shells of various thicknesses to study overshooting over several tens of convective turnover times or more. Apart from the differences in the dimensional and geometrical setup, the main difference is that the analysis by Pratt et al. (2016, 2017) is focused on a prototype of a young Sun, where plumes are launched inward, along an increasing density gradient, from a subsurface convection zone. Of course the opposite is true in our models, where plumes are launched outwards from a convective core.

Figure 5 shows probability density functions (PDFs) of the plume overshooting length in units of the pressure scale height at the CB ${H}_{p}^{\mathrm{CB}}$ (panel (a)) for each stellar age (z: ZAMS, purple; m: midMS, cyan; and T: TAMS, orange) and an illustration of the overshoot of plumes in our simulations as a function of the latitude θ (panel (b)). For reference, the pressure scale heights at the CB for the three models are ${H}_{p,{\rm{Z}}}^{\mathrm{CB}}=0.095{R}_{\star }$, ${H}_{p,{\rm{m}}}^{\mathrm{CB}}=0.079{R}_{\star }$, and ${H}_{p,{\rm{T}}}^{\mathrm{CB}}=0.052{R}_{\star }$.

Zoom In Zoom Out Reset image size

The PDFs presented in panel (a), which have been normalized such that the heights of their bars sum up to 1, show that for all stellar ages a large fraction of plumes only reach an approximate maximum overshooting length of ${\rm{\Delta }}{r}_{\mathrm{ov}}/{H}_{p}^{\mathrm{CB}}=0.035$ (vertical dotted lines). There also appears to be a second overshoot layer, directly above the first, which plumes reach more infrequently. This deeper overshoot layer, which is also observed in the 2D simulations by Pratt et al. (2017), appears fairly accessible in the ZAMS case but becomes harder to reach with increasing stellar age; few plumes are able to overshoot beyond ${\rm{\Delta }}{r}_{\mathrm{ov}}/{H}_{p}^{\mathrm{CB}}=0.035$ in the TAMS case. This can be quantified by calculating the percentage of plumes overshooting beyond ${\rm{\Delta }}{r}_{\mathrm{ov}}/{H}_{p}^{\mathrm{CB}}=0.035$ for the three ages: while 67% of plumes overshoot beyond this threshold for ZAMS, this figure drops to 43% and 31% for midMS and TAMS, respectively. While the ZAMS value is probably slightly inflated by some potential false positives found beyond ${\rm{\Delta }}{r}_{\mathrm{ov}}/{H}_{p}^{\mathrm{CB}}\sim 0.2$, the downward trend with stellar age is clear. The reason for this decline in the number of plumes overshooting beyond this threshold, given that the three models possess comparable generation spectra (Figure 4), can be safely attributed to the emergence of the spike in the BVF in older stars, caused by a steep compositional gradient.

Panel (b) of Figure 5 shows the latitudinal distribution of overshooting plumes for each model as a function of the normalized overshooting length. Here each vertical bar represents a plume (the plume size is not taken into consideration, therefore all bars have equal width), with overlapping plumes resulting in darker regions. There appears to be no preferential latitude for plume overshoot, although fewer plumes are observed near the stellar poles. The mean overshooting length is also indicated for each model by means of a horizontal dashed line; this appears to decrease with stellar age, with a significant drop from ZAMS to midMS (from ${\rm{\Delta }}r/{H}_{p}^{\mathrm{CB}}\approx 0.09$ to ≈0.05), and a minor one from midMS to TAMS (from ≈0.05 to ≈0.04). These results are comparable with the results from the 2D simulations by Pratt et al. (2017), as well as the 3D simulations of a 2 M_⊙ A star by Browning et al. (2004); in the latter, the authors find overshooting depths of roughly 0.2H_p , and notice the overshooting length shrinking with increasing Reynolds number. This would potentially reconcile the results even more, as they consider Reynolds numbers that are 2.5–50 times smaller than the one considered in our ZAMS model. While the overshooting lengths across all models could potentially be inflated by the boosted stellar luminosities used, the observed trend is consistent with the percentage of plumes overshooting to the second overshooting layer mentioned above.

In order to understand the relative importance of overshooting plumes and convective eddies in our generation spectra, we carry out a direct comparison between the spectra shown in Figure 4 and the analytical prescriptions for both formalisms. To begin with, we consider two examples of eddy-dominated generation spectra from Kumar et al. (1999; labeled K) and Lecoanet & Quataert (2013; labeled LD for "Lecoanet discontinuous," to indicate it is the profile obtained with a discontinuous stratification profile at the CB) according to which the wave energy spectrum is described by the following power laws:

$$\begin{eqnarray}&&{E}_{\mathrm{kin}}^{{\rm{K}}}\propto {f}^{-4.34},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{E}_{\mathrm{kin}}^{\mathrm{LD}}\propto {f}^{-6.5}.\end{eqnarray} \tag{ 15 }$$

We also consider a third power-law spectrum, obtained by Rogers et al. (2013; labeled R) from 2D simulations of a 3 M_⊙ star. Their spectra are actually best described by two distinct power laws; for simplicity we pick the power-law exponent stretching the largest portion of the spectrum, encompassing from low frequencies up to f ∼ 100 μHz. Since Rogers et al. (2013) do not analyze the same rotation rates used here (the models in both Figure 4(a) and Figure 4(b) have Ω = 10⁻⁵ rad s⁻¹), we take an average of their results obtained with the two closest rotation rates, resulting in

$$\begin{eqnarray}&&{E}_{\mathrm{kin}}^{{\rm{R}}}\propto {f}^{-0.85}.\end{eqnarray} \tag{ 16 }$$

With regards to IGWs generated by plume overshoot, we adopt the formalisms developed by Montalbán & Schatzman (2000) and Pinçon et al. (2016), which are based on the plume models by Rieutord & Zahn (1995) and Zahn (1991). These studies are only valid at the high Péclet numbers (i.e., Pe >> 1), with the Péclet number being defined as

$$\begin{eqnarray*}&&\mathrm{Pe}=\displaystyle \frac{{uL}}{\kappa },\end{eqnarray*}$$

where u and L are the characteristic velocity and length scale of the convective motions, respectively. At the radii of the spectra shown in Figure 4 all our models show Pe values of order ∼10⁴, allowing us to apply their formalism with confidence. The plume energy spectrum depends on the plume timescale τ_pl, which is approximated as

$$\begin{eqnarray}&&{\tau }_{\mathrm{pl}}\sim \displaystyle \frac{{\rm{\Delta }}{r}_{\mathrm{ov}}}{{v}_{\mathrm{pl}}},\end{eqnarray} \tag{ 17 }$$

where v_pl is the velocity of a plume and Δr_ov is its overshooting length. Its reciprocal gives the plume frequency,

$$\begin{eqnarray}&&{f}_{{\rm{pl}}}\sim \displaystyle \frac{{v}_{{\rm{pl}}}}{{\rm{\Delta }}{r}_{{\rm{ov}}}},\end{eqnarray} \tag{ 18 }$$

where the energy spectrum takes an exponential form in both Montalbán & Schatzman (2000) and Pinçon et al. (2016) cases, with the general form of

$$\begin{eqnarray}&&{E}_{\mathrm{kin}}\left(f\right)\propto {e}^{-{\left(f/{f}_{\mathrm{pl}}\right)}^{2}}.\end{eqnarray} \tag{ 19 }$$

While this simplified expression yields the kinetic energy generation spectrum for a specific plume frequency, it does not provide any information about the frequency distribution of such plumes; each of our simulations of course possesses a unique f_pl distribution, obtained from the plume analysis described in Section 4.2, with the normalized plume frequency PDF for the midMS model shown in Figure 6. The distribution peaks at ≈6 μHz (black, dotted–dashed vertical line) and appears well described by two exponential fits (dashed magenta lines) having gradients of −1/β₁ and −1/β₂, with β₁ = 6.2 μHz and β₂ = 40.9 μHz, with a crossover point at roughly f_pl ≈ 20 μHz. The f_pl distributions of the other models are similarly well described by two exponential fits. We therefore combine the double exponential nature of the f_pl distributions from our models with the original plume prescriptions by Montalbán & Schatzman (2000; labeled M) and Pinçon et al. (2016; labeled P) to ensure the appropriate plume frequency contributions are considered. This is achieved through an extra ${e}^{-{f}_{\mathrm{pl}}/{\beta }_{n}}$ factor (with n = 1, 2); the adapted plume prescriptions become

$$\begin{eqnarray}&&{E}_{\mathrm{kin}}^{{\rm{M}}}\left(f\right)\propto \displaystyle \sum _{n}{\int }_{0}^{\infty }{e}^{-{f}_{\mathrm{pl}}/{\beta }_{n}}\,{e}^{-{\left(f/{f}_{\mathrm{pl}}\right)}^{2}}{{df}}_{\mathrm{pl}},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{E}_{\mathrm{kin}}^{{\rm{P}}}\left(f\right)\propto \displaystyle \sum _{n}{\int }_{0}^{\infty }{e}^{-{f}_{\mathrm{pl}}/{\beta }_{n}}\,{f}_{\mathrm{pl}}^{-1}{e}^{-{\left(f/4{f}_{\mathrm{pl}}\right)}^{2}}{{df}}_{\mathrm{pl}},\end{eqnarray} \tag{ 21 }$$

where −1/β_n is the slope of each exponential fit to the unique f_pl PDF of each model, as shown in Figure 6 by the dashed magenta lines.

Zoom In Zoom Out Reset image size

Figure 7 showcases the results of the comparison between our simulation spectra (black, full lines) and Equations (14)–(16) and (20)–(21) for all three stellar ages, both just below (top panel) and above (bottom panel) the CB. Below the CB (top panel), the adapted plume prescriptions appear to describe the generation spectra exceedingly well throughout most of the frequency range and for all three stellar models. Out of the two plume generation spectra analyzed, the Pinçon et al. (2016) prescription ${E}_{\mathrm{kin}}^{{\rm{P}}}$ (orange dashed line) appears the best fit to the simulation data, while ${E}_{\mathrm{kin}}^{{\rm{M}}}$ (blue, dashed) appears to produce a profile whose gradient changes too sharply around f ≈ 10 μHz compared to the simulation spectra, resulting in either an excessive drop in amplitude at high frequencies, or a profile that is too flat in the low-f range. This result is in line with the findings of Edelmann et al. (2019), who also found their generation spectra to be well described by plume-driven motions.

Zoom In Zoom Out Reset image size

It is worth noting that all three spectra are also sufficiently well described by a broken power law with the best fits among the power laws used being ${E}_{\mathrm{kin}}^{{\rm{R}}}$ (green, dotted–dashed) for low f and ${E}_{\mathrm{kin}}^{{\rm{K}}}$ (magenta, dotted–dashed) for high frequencies (although both of these power laws look somewhat too steep for our data). This is well aligned with the numerical results of Rogers et al. (2013) who also found a broken power law with exponents of α ∼ 1 for low and intermediate frequencies and α ∼ 4 (therefore very similar to the power-law exponent in ${E}_{\mathrm{kin}}^{{\rm{K}}}$) for high f for stellar rotation rates close to those of the three models shown here. However, it is clear that this broken power-law fit does not represent as accurate a description as the plume prescriptions, particularly as it appears to overestimate substantially the spectrum at low frequencies.

At the location just above the CB (bottom panel), the situation is somewhat different. All three spectra have been affected in the mid- and high-frequency ranges, and their profiles now appear better described by a single power law with an index similar to that of ${E}_{\mathrm{kin}}^{{\rm{R}}}$. While this is true for all models, the ZAMS spectrum still appears adequately well described by a Pinçon et al. (2016) plume-driven profile; however, the agreement between the spectra and their respective ${E}_{\mathrm{kin}}^{{\rm{P}}}$ fits deteriorates with stellar age. This may be expected given those theoretical models do not account for BVF spikes as the star ages.

We now move on to consider the propagation of waves throughout the entirety of the radiation zone. Figure 8 shows equatorial slices illustrating the propagation of IGWs through the RZ across the three stellar ages: ZAMS (labeled Z, top row), midMS (m, middle row), and TAMS (T, bottom row). While in the ZAMS model the generated IGWs are seen to propagate freely throughout the whole radiation zone, in the midMS and TAMS models the propagation of waves is significantly affected beyond their BVF spikes (the extent of which is marked by red dashed circles), leaving the rest of their RZs more sparsely populated by IGWs. This is further highlighted by the velocity amplitudes of the slices, which clearly decrease with age (see the associated color bar).

Zoom In Zoom Out Reset image size

The depletion of IGWs in the top part of the RZ of older stars is due to a combination of three separate effects. The one we believe to be the most important, both in these simulations and in real stars, is IGWs becoming evanescent as they propagate toward the stellar surface. As explained by Ratnasingam et al. (2020) in their 2D cylindrical model, the 2D linearized anelastic equations can be reduced to a single second-order differential equation describing the evolution of IGWs:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\partial }^{2}a}{\partial {r}^{2}} & & +\left(\displaystyle \frac{{N}^{2}}{{\omega }^{2}}-1\right){k}_{h}^{2}a\\ & & +\displaystyle \frac{1}{2}\left[\displaystyle \frac{\partial {h}_{\rho }}{\partial r}-\displaystyle \frac{{h}_{\rho }^{2}}{2}+\displaystyle \frac{{h}_{\rho }}{r}\right]a+\displaystyle \frac{1}{4{r}^{2}}a\,=\,0,\end{array}\end{eqnarray} \tag{ 22 }$$

where $a={v}_{r}{\overline{\rho }}^{1/2}{r}^{3/2}$. IGWs maintain their wave-like properties as long as the ratio between the oscillatory term (OT; second term on the left-hand side of the equation) and the density term (DT; third term on the left-hand side of the equation) remains above unity, i.e.,

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{OT}}}{{\rm{DT}}}=\displaystyle \frac{\left(\tfrac{{N}^{2}}{{\omega }^{2}}-1\right)\tfrac{{m}^{2}}{{r}^{2}}}{\tfrac{1}{2}\left[\tfrac{\partial {h}_{\rho }}{\partial r}-\tfrac{{h}_{\rho }^{2}}{2}+\tfrac{{h}_{\rho }}{r}\right]}\gt 1.\end{eqnarray} \tag{ 23 }$$

Within the radiation zone the value of this ratio decreases with radius, potentially reaching values below one depending on the frequency chosen and the value of m. The radius at which the ratio reaches unity is called the turning point, and IGWs lose their wave-like behavior beyond it. Ratnasingam et al. (2020) also observed the location of the turning point receding to smaller r/R_⋆ values with increasing age for waves of the same frequency, with IGWs in older stars therefore becoming evanescent earlier in their propagation, hence contributing to the scarcity of IGWs observed in our midMS and TAMS models. The reason older stars are more susceptible to this phenomenon is the overall decrease in the BVF as a function of stellar age, as seen in Figure 1, which results in a decreased OT for all frequencies. To put the magnitude of this effect into perspective, a wave of frequency f = 15 μHz (m = 1) loses its wave-like behavior past r/R_⋆ ≈ 0.72 for ZAMS, but its turning point recedes significantly to r/R_⋆ ≈ 0.56 and r/R_⋆ ≈ 0.38 for midMS and TAMS, respectively. This effect, which has not been considered before when analyzing the propagation of IGWs in stars and their impact on stellar dynamics and observational signatures, is believed to also be responsible for the change in wave structure seen across the RZ in Figure 8 as a progressively increasing fraction of IGWs become evanescent as a function of radius. This is supported by the magenta and green dashed circles in the right panels of Figure 8, which represent the locations at which m = 1 waves with frequencies of 15 μHz (green, inner circle) and 10 μHz (magenta, outer) ⁷ become evanescent in each model as predicted by Equation (23), roughly coinciding with the location where the wave structure changes. While we see this effect in our simulations, it is likely to be much more relevant in actual stars, where viscous and thermal damping are far less relevant. Moreover being a linear effect, it is unaffected by most of the simulation assumptions and shortcomings.

On top of IGW propagation in older stars being more restricted due to turning-point limitations, a substantial portion of the generated IGW frequencies in the midMS and TAMS models becomes trapped within the respective BVF spikes; this is caused by the drop in N at the end of the BVF spike, which causes high-frequency IGWs to no longer be permitted to propagate as ω < N would be violated. This phenomenon is visualized in the left panels of Figure 8 for midMS and TAMS, where it is possible to see IGWs trapped within the BVF spike (the extent of which is marked by the dashed red circles). There we see the angle of the phase lines change dramatically, indicating reflection. As shown in the BVF profile from Figure 1, for TAMS this jump in N² is of a factor of ≈30, meaning frequencies in the 70 ≲ f ≲ 400 μHz range are trapped within the BVF spike.

Lastly, the enhanced radiative damping within the BVF spike due to the large value of N within the spike itself also contributes to the scarcity of IGWs in the top part of the radiative zones in midMS and TAMS. Linear theory suggests the velocity amplitudes of waves propagating through the radiation zone will evolve according to several effects, such as density stratification, geometrical effects, and radiative damping. Ratnasingam et al. (2019), building on the work by Press (1981) and Kumar et al. (1999), provide the following expression for the linear wave velocity amplitude:

$$\begin{eqnarray}&&{v}_{r}\propto {\left(\displaystyle \frac{{r}_{0}}{r}\right)}^{3/2}\sqrt{\displaystyle \frac{{\rho }_{0}}{\rho }}{\left(\displaystyle \frac{{N}^{2}-{\omega }^{2}}{{N}_{0}^{2}-{\omega }^{2}}\right)}^{1/4}{e}^{-\tau /2},\end{eqnarray} \tag{ 24 }$$

with ω = 2π f, where r₀ is the radius at which waves begin propagating, and N₀ and ρ₀ are the BVF and density at r₀; τ is the damping opacity, whose full expression is

$$\begin{eqnarray}\begin{array}{rcl}\tau & = & {\displaystyle \int }_{{r}_{0}}^{r}\displaystyle \frac{{\gamma }_{\mathrm{rad}}}{| {v}_{{\rm{g}}}| }\,{dr},\\ & = & {\displaystyle \int }_{{r}_{0}}^{r}\displaystyle \frac{\kappa {\left[l\left(l+1\right)\right]}^{3/2}}{{r}^{3}{\omega }^{4}}{N}^{3}\sqrt{1-\displaystyle \frac{{\omega }^{2}}{{N}^{2}}}\,{dr},\end{array}\end{eqnarray} \tag{ 25 }$$

where γ_rad is the radiative damping rate and v_g is the wave group velocity. The amount of radiative damping of waves is therefore very sensitive to the values of N, implying that a specific frequency ω of a given l value would be significantly more damped within the BVF spike of an older star than in the corresponding radial location of a younger one due to the higher value of N. To quantify this increased radiative damping within the BVF spike of older stars, we use Equation (25) to calculate the damping factor $\exp (-\tau /2)$ for an f = 10 μHz wave with l = 2 at r = 0.5 R_⋆ for the ZAMS and TAMS models. We find that at the given location, the TAMS wave is subject to a damping that is ∼5 times stronger than its ZAMS counterpart. Radiative damping for a fixed l value is at its strongest in low-frequency waves (ω/N ≪ 1) since γ_rad/∣v_g ∣ ∝ ω⁻⁴, meaning these waves would be particularly affected. With Talon et al. (2002) treating viscosity and thermal diffusion as additive when considering wave damping, it is possible viscosity dampens IGWs in our simulations at least as much as κ; in any case, comparing radial velocity amplitudes from our simulations to those obtained using stellar κ values, we find waves with f ≳ 5 μHz and m ≲ 5 are consistent with MESA. However, it is important to note that while changes in radiative damping as a function of age must be considered in our simulations due to the excessive thermal diffusivities used for numerical reasons, we do not expect this phenomenon to be significant in real stars. As real stars feature much smaller values for κ, the changes in radiative damping across different stellar ages is negligible for most wavenumber and frequency permutations unless considered near the stellar surface. Our focus is, therefore, centered on the turning point and the IGW trapping mechanisms, both nondamping in nature, which also play an important role in IGW propagation in real stars.

The two effects linked to the BVF spikes, IGW trapping, and enhanced radiative damping (in simulations) can be seen at play in the snapshots of Figure 9, which shows a comparison for the radial velocity v_r at the same r/R_⋆ locations for ZAMS (left) and midMS (right); specifically, we look at the flows at the CB (top) as well as 0.02 R_⋆ (middle) and 0.15 R_⋆ (bottom) above it. The two stars have qualitatively similar flows at the CB (top row), with both showing plumes hitting the convective interface and generating IGWs. At 0.02 R_⋆ above the CB (middle row), which for the midMS model roughly corresponds to the midway point within its BVF peak, both stars show significant IGW generation by the plumes observed at the CB (these can be seen as concentric circles in the figure). Looking past the BVF peak, at 0.15 R_⋆ above the convective interface (bottom row) the ZAMS star still shows clear signs of both small- and large-scale IGW propagation throughout the (ϕ, θ) domain; the midMS model, on the other hand, shows lower amplitude fluctuations, as both low- and high-frequency waves have been heavily affected by the presence of the BVF spike, with the former being particularly susceptible to the enhanced radiative damping and the latter remaining trapped in the BVF spike. This overall trend is likely to be reproduced in real stars where low-frequency waves will be radiatively damped (though at a lower rate than these simulations) and high-frequency waves will be trapped within the spike.

Zoom In Zoom Out Reset image size

In order to visualize the frequencies excited in a particular model and how its spectrum evolves as a function of stellar radius, we create spectral heat maps integrated over all l by sampling multiple longitudinal points across the stellar equator at all simulation time steps; Figure 10 shows such spectral heat maps for the radial velocity v_r as a function of both frequency f and fractional stellar radius r/R_⋆ for models 7ZR10 (labeled "z," left) and 7mR10 ("m," right). The convection zones, which stretch up to r/R_⋆ ≈ 0.22 and r/R_⋆ ≈ 0.13, respectively, are characterized by the presence of all possible frequencies, thanks to convective motions being distributed over a large range of timescales, with f ≲ 30 μHz (corresponding to timescales longer than 0.39 days) being the most dominant range.

Zoom In Zoom Out Reset image size

Within the radiation zone, the signal is instead constrained by the BVF (white dashed lines), with a sharp drop for f > N/2π as expected for IGWs. This behavior can also be observed in the inset axes of the spectral heat map for model 7mR10 (right), which shows a zoomed-in view around the BVF spike, where the signal is observed to mirror the N/2π profile closely. Exempt from the restriction of the BVF are the fundamental modes and gravity modes, the strong vertical features which are seen spanning significant radial portions of the RZ. Low-frequency waves appear strongly damped, with frequencies below approximately 10 μHz lacking throughout most of the RZ in both models. In our simulations this is largely due to enhanced thermal and viscous diffusion, however, one would expect a lack of low frequencies in real stars as well due to thermal diffusion.

Figure 11 shows integrated frequency spectrum profiles for v_r for all three stellar ages, obtained by considering slices at fixed r/R_⋆ from their respective spectral heat maps (such as those shown in Figure 10). Each panel (top: ZAMS, middle: midMS, and bottom: TAMS) shows spectrum profiles at five distinct fractional radial locations: in each case, the first profile is taken at approximately 0.02 R_⋆ below the CB (maroon lines), while the remaining four are spread throughout the radiation zone. The profiles within the radiation zones are taken at roughly 0.03 R_⋆ above the CB, which for the midMS and TAMS models falls within their BVF spikes (green lines); at the location where the BVF reaches its maximum value, excluding BVF spikes (orange lines); at r/R_⋆ = 0.60 (magenta lines) and at approximately 0.03 R_⋆ below the outer edge of the domain (light blue lines). The minimum BVF value throughout the radiation zone, N_min, for each model is indicated by means of a vertical gray dotted line; additionally, the value of the BVF at the chosen r/R_⋆ closest to the surface, N_surf, is indicated by light blue dashed lines, matching the color of the profile obtained at the same location.

Zoom In Zoom Out Reset image size

The most obvious difference between the models is the relative lack of evolution of the ZAMS spectrum as a function of r/R_⋆ when compared to its older counterparts; the midMS and TAMS profiles are seen to vary considerably from within the BVF spike (green profiles) to the outer domain boundary, due to several effects such as filtering by the BVF spike, radiative damping, and a more restrictive BVF profile. While the profiles taken just above the CZ (green) are qualitatively similar across age, the phenomena mentioned above mean the surface signatures (light blue) of the three stellar ages differ greatly. The steepness of the spectra taken near the surface appears to increase from ZAMS, which presents a flatter spectrum with a power-law exponent of roughly −0.6, to the older midMS and TAMS models whose exponents are found to be roughly −2. This effect appears to be a consequence of IGWs with frequencies approaching N_surf not reaching the outer boundary of the domain in the midMS and TAMS models, as evidenced by their gradual drop-offs in signal around N_surf (light blue dashed vertical lines) as opposed to the sharp decline observed for ZAMS. This is believed to be caused by large quantities of IGWs becoming evanescent in older stars due to the position of the turning point receding inwards with age; IGWs with frequencies approaching the local N² value at any particular r/R_⋆ would be the most heavily affected by the turning point receding inward and would become evanescent further away from the stellar surface than waves with lower frequencies, as suggested by Equation (23). The steepness of all models considered is however consistent with other 2D and 3D simulations, such as those of Rogers & MacGregor (2010), Rogers et al. (2013), Alvan et al. (2014), Augustson et al. (2016), and Edelmann et al. (2019), who observed spectra with power-law exponents roughly between −1 and −3.

Both Figures 10 and 11 exhibit gravity and fundamental modes in the form of vertical features; a more thorough analysis on their nature and the accuracy of their frequencies is carried out in Figure 12, which shows the frequency spectra of the radial velocity v_r for all radii for the ZAMS (labeled "Z," left column) and midMS (labeled "m," right column) models for individual angular degrees: l = 2 (top panels), l = 3 (middle panels), and l = 4 (bottom panels). This allows us to disentangle the contributions from different l values and obtain a clearer view of the resonant coherent modes within the stellar cavity. The nature of the modes can be established by the number of radial nodes they present; the modes without any nodes are f-modes, or fundamental modes, while modes exhibiting an increasing amount of radial nodes with decreasing frequency are g-modes, or gravity modes. The modes observed within each panel are compared to their predicted frequencies obtained from the stellar oscillation code GYRE (Townsend & Teitler 2013; Townsend et al. 2018), which are represented by means of vertical white lines and labeled according to the Eckart–Osaki–Scuflaire–Takata scheme (see, e.g., Aerts et al. 2010); this means the f-mode is indicated by 0, g-modes are represented by negative numbers, and p-modes by positive numbers. The latter are however not picked up in our simulations due to the anelastic nature of the code. ⁸

Zoom In Zoom Out Reset image size

For the ZAMS model there is good agreement between the simulation and GYRE concerning the g-modes with g₂, g₃, and g₄ matching particularly well, especially for l = 2 and l = 3. The level of agreement for the f-modes varies more noticeably across l; for l = 2 the observed f-mode sits roughly halfway between the frequencies predicted by GYRE for the g₁ mode and the f-mode, but for l = 4 the simulation and the data from GYRE match closely. This improved agreement for higher values of l is actually in contrast with the analogous analysis from Edelmann et al. (2019), where the level of agreement between GYRE and the simulation data for the fundamental modes was found to instead decrease with increasing l.

The agreement with GYRE is however more erratic for the midMS model. The high-frequency f-mode (clearly located within the BVF spike at f > 150 μHz) appears completely mismatched with the frequencies expected by GYRE and although the g-modes match well for l = 2 and l = 3, their agreement worsens for l = 4 where the simulation overestimates the frequencies of modes g₁ and g₂. The reason for this inconsistent agreement between our midMS model and GYRE is not clear, but two likely explanations have been singled out (as suggested by R. H. D. Townsend in a private conversation). One option is that it could be caused by an eigenfrequency distortion in GYRE due to the software not properly dealing with mass conservation within the BVF spike. Alternatively, it could be a natural limitation of our anelastic code when dealing with older stars: the emergence of the BVF spike in the midMS model (and TAMS models) increases the maximum frequency allowed for the propagation of gravity waves; this has the consequence of overlapping the frequency ranges for gravity waves and sound waves (which is demonstrated by our midMS simulation finding the f-modes at frequencies higher than the GYRE p₁ frequency). Any modes falling within this overlap region would take a "mixed mode" character (Unno et al. 1989), and the anelastic nature of our code might struggle to interpret and reproduce these frequencies properly.

IGWs play an important role in the transfer of angular momentum within stellar interiors thanks to their ability to propagate and dissipate within radiation zones; this is particularly true for massive stars, such as those studied here, where the stellar density stratification causes some IGW amplitudes to grow during their propagation toward the stellar surface.

Transport of angular momentum via IGWs is however only possible when the waves are attenuated, as demonstrated by Eliassen & Palm (1960). There exist three main mechanisms that can dissipate IGWs within stellar interiors: radiative damping, the formation of critical layers, and nonlinear wave breaking due to large amplitude. Radiative damping (Equations (24)–(25)) affects wave amplitudes nonlinearly and is very sensitive to the wave frequency, being at its most effective for low frequencies. Critical layers are the product of a corotation resonance, which is achieved when the Doppler-shifted frequency of a wave,

$$\begin{eqnarray}&&\omega (r)={\omega }_{\mathrm{gen}}-m\left[{\rm{\Omega }}(r)-{{\rm{\Omega }}}_{\mathrm{gen}}\right],\end{eqnarray} \tag{ 26 }$$

becomes zero, where ω_gen is the frequency of such a wave at the generation radius (i.e., at the CB), and Ω_gen and Ω(r) are the angular velocities of the flow at the CB and at radius r. Interactions between IGWs and critical layers can either result in wave attenuation or in the IGWs being reflected/transmitted, depending on the value of the Richardson number at the critical layer (Alvan et al. 2013). While our models are started with a solid body rotation rate of either 10⁻⁶ or 10⁻⁵ rad s⁻¹ (see Table 1), they quickly develop a weak differential rotation between the convective and radiative regions and near the simulated surface; given the strength of such differential rotation, we would expect critical layers to occur in these regions for waves possessing frequencies ω_gen ∼ m × 10⁻⁷ rad s⁻¹. Lastly, nonlinear wave breaking can occur either through convective instability or a Kelvin–Helmholtz instability; the former is restricted to regions where N ∼ 0, such as the interfaces between convective and radiative zones, meaning a Kelvin–Helmholtz instability is the most likely source of wave breaking within the bulk of stellar radiative zones. Kelvin–Helmholtz instabilities develop when the destabilizing shearing motions between adjacent layers of a stratified shear flow prevail over the buoyancy; numerically, this is controlled by the Richardson number:

$$\begin{eqnarray}&&\mathrm{Ri}=\displaystyle \frac{{N}^{2}}{{k}_{r}^{2}{u}_{h}^{2}},\end{eqnarray} \tag{ 27 }$$

with Ri ≲ 1 yielding a Kelvin–Helmholtz instability. For IGWs, this instability criterion can be expressed as = k_h u_h /ω ≳ 1 (Press 1981). It is unclear whether this instability occurs due to IGWs in stellar radiative interiors. While these 3D simulations do not demonstrate this instability, it was seen in previous 2D simulations, with the difference being due to slightly larger forcing and lower damping in 2D compared to 3D. Given it is unclear what combination of forcing/damping might accurately represent the radiative interior of a star (if any), the jury is still out. However, the work by Le Saux et al. (2023) indicates that nonlinearity of the waves is likely only in limited circumstances.

Whenever an IGW is attenuated, it transfers angular momentum to the mean flow (Eliassen & Palm 1960). The dynamical interaction between the waves and the mean flow are complex and highly nonlinear, being dependent on the properties of both. Neglecting meridional circulation, the evolution of angular momentum transport via IGWs is described by

$$\begin{eqnarray}\begin{array}{rcl}\rho & & {\partial }_{t}\left({r}^{2}\langle {\rm{\Omega }}\rangle \right)\\ & & =\,-\displaystyle \frac{1}{{r}^{2}}{\partial }_{r}\left({r}^{2}\rho r\langle \sin \theta {v}_{r}{v}_{\phi }\rangle \right)+\displaystyle \frac{1}{{r}^{2}}{\partial }_{r}\left(\rho \nu {r}^{4}{\partial }_{r}\langle {\rm{\Omega }}\rangle \right),\end{array}\end{eqnarray} \tag{ 28 }$$

where angle brackets 〈·〉 represent horizontal averages and $\langle {\rm{\Omega }}\rangle ={\int }_{0}^{\pi }{\rm{\Omega }}{\sin }^{3}\theta d\theta /{\int }_{0}^{\pi }{\sin }^{3}\theta d\theta$ (Zahn 1992; Maeder & Zahn 1998; Mathis & Zahn 2004). As seen from Equation (28), the evolution of the mean flow is controlled by the divergence of the horizontally averaged Reynolds stress, ${F}_{\mathrm{Rey}}=\rho r\langle \sin \theta {v}_{r}{v}_{\phi }\rangle$, and by the divergence of the viscous flux. Comparing the magnitudes of these two terms allows us to establish which term is dominant and, consequently, understand the dynamics of the stellar interiors. With that aim in mind, we define the ratio between the divergences of the Reynolds stress and the viscous flux as

$$\begin{eqnarray}&&\xi =\displaystyle \frac{{\partial }_{r}\left({r}^{2}\rho r\langle \sin \theta {v}_{r}{v}_{\phi }\rangle \right)}{{\partial }_{r}\left(\rho \nu {r}^{4}{\partial }_{r}\langle {\rm{\Omega }}\rangle \right)}.\end{eqnarray} \tag{ 29 }$$

The evolution of ξ over time as a function of radius is shown in Figure 13 for both ZAMS (z, left) and midMS (m, right) models. The results from the TAMS model are not shown because the low amplitudes of IGWs, due to the effects described above, mean this model is viscously dominated everywhere within the radiative region. In both models, the convection zones are dominated by Reynolds stresses (blue); immediately above the CB, both models present an area dominated by viscous flux (red), due to the strong shear flow which develops between the convective and radiative regions. This shear flow partly overlap with the overshooting regions of the two models, where the value of the BVF is low, causing the overshoot regions to be potentially susceptible to shear instability or convective instability, the latter of which could result in nonlinear wave breaking as mentioned earlier. Additionally, both stars exhibit another shear flow at the outer boundary of their domain, with the midMS model featuring a shear flow that is seemingly more radially extended. Above the first shear flow some differences between the two models can be observed: the bulk of the radiation zone in the ZAMS model is dominated by Reynolds stresses, indicating the prominence of IGWs in the dynamics of this region at this age; the same region for the midMS model, on the other hand, exhibits both Reynolds stresses and viscous fluxes alternating in a seemingly stochastic nature. This disparity in the importance of IGWs within the bulk of the radiation zones is due to the same effects described in Section 4.4, i.e., the trapping of waves within the BVF spike and waves becoming evanescent at a lower fractional radius within the radiative zone. While we expect these effects to be relevant in stars, they are exacerbated in these simulations due to the enhanced viscosity. Hence, while we expect older stars to have lower Reynolds stress amplitudes, we also expect them to have significantly lower viscous flux. Which is dominant in actual stars is unclear.

Zoom In Zoom Out Reset image size

The conclusions drawn from the ξ plots are reinforced by Figure 14, which shows the time-averaged radial profiles of the angular velocity fluctuation Ω for all three stellar ages (models 7ZR10, 7mR10, and 7TR10). As mentioned previously, these three models are all set up with a uniform stellar rotation of $\overline{{\rm{\Omega }}}={10}^{-5}\,\mathrm{rad}\,{{\rm{s}}}^{-1}$ (as indicated in Table 1), but are seen to develop differentially rotating radial profiles. In all models the convective cores exhibit a sizeable positive angular velocity fluctuation, while the shear flows above them are rotating slower than the background value. The depths of these shear flows are highlighted by the color-coded shaded areas; all models present similarly deep shear layers, with their radial extents ranging from 0.116 R_⋆ (midMS) to 0.174 R_⋆ (TAMS). This results in the central part of the radiative zone rotating slower than the background ration for the ZAMS model, while in the midMS and TAMS models the same region rotates slightly faster than the initial uniform rotation. Near the outer boundary we observe another band of negative Ω in both models, representing the surface shear flows observed in the ξ plot, which are likely maintained by IGWs.

Zoom In Zoom Out Reset image size

Finally, Figure 14 shows the radial profile for the rate of specific angular momentum transport, ${\partial }_{t}\left({r}^{2}\langle {\rm{\Omega }}\rangle \right)$, as obtained from Equation (28) for the ZAMS model. This profile shows strong similarities to the Ω profile in Figure 14, with values of $| {{\rm{\partial }}}_{t}\left({r}^{2}\langle {\rm{\Omega }}\rangle \right)| \sim {10}^{10}\,{{\rm{c}}{\rm{m}}}^{2}\,{{\rm{s}}}^{-2}$ recorded in the core and ∼10⁶ cm² s⁻² in the bulk of the RZ. This latter value implies that IGWs would be able to alter the rotation rate of the radiative zone in a timescale,

$$\begin{eqnarray}&&{t}_{{\rm{IGW}}}=\displaystyle \frac{{R}_{{\rm{RZ}}}^{2}{\rm{\Omega }}}{{{\rm{\partial }}}_{t}\left({r}^{2}\langle {\rm{\Omega }}\rangle \right)}\sim 8\times {10}^{8}\,{\rm{yr}},\end{eqnarray} \tag{ 30 }$$

where R_RZ is the radial extent of the radiative zone. This timescale is significantly longer than the one estimated by Fuller et al. (2014) for solar-type stars, giants, and subgiants. It is however important to highlight that the models presented in this work self-consistently account for both prograde and retrograde IGWs and autonomously develop differential rotation as a result of their IGW propagation properties.