Chemical Mixing Induced by Internal Gravity Waves in Intermediate-mass Stars

We generated the stellar models for three masses (3, 7, and 20 M_⊙) at three ages (zero age main sequence (ZAMS), core hydrogen mass fraction, X_c = 0.7), mid-main sequence (midMS), X_c = 0.35), and terminal age main sequence (TAMS), X_c = 0.01) from MESA. We set the mixing length parameter to 1.8 and the convective overshoot profile to exponential, implemented the mass loss through van Loon stellar wind scheme and considered zero rotation in our models. The inlists used to generate the models are given in Zenodo at doi:10.5281/zenodo.2596370

We calculated the Brunt–Väisälä frequency, N based on the work of Rogers et al. (2013),

$$\begin{eqnarray}&&{N}^{2}=\displaystyle \frac{\bar{g}}{\bar{T}}\left(\displaystyle \frac{d\bar{T}}{d\bar{r}}-(\gamma -1)\bar{T}{h}_{\rho }\right).\end{eqnarray} \tag{ 1 }$$

where $\bar{g}$, $\bar{T}$, γ, and h_ρ are the reference state gravity, temperature, adiabatic index, and negative inverse density scale height. The location of the convective–radiative interface is determined by the sign of the Brunt–Väisälä frequency, where

$$\begin{eqnarray}&&{N}^{2}\gt 0\end{eqnarray} \tag{ 2 }$$

is the radiation zone and

$$\begin{eqnarray}&&{N}^{2}\lt 0\end{eqnarray} \tag{ 3 }$$

is the convection zone. The background density and the Brunt–Väisälä frequency profiles of all the models considered in this work are shown in Figure 1.

Zoom In Zoom Out Reset image size

To study IGWs in stellar interiors, we solved the Navier–Stokes equations in the anelastic approximations using 2D hydrodynamic simulations considering an equatorial slice of the star with stress-free, isothermal, and impermeable boundary conditions (Rogers et al. 2013). The numerical model followed in this work is similar to that of Rogers et al. (2013) and the 2D simulations are run by R.P Ratnasingam (more details on the simulations can be found in Ratnasingam 2020, R. P. Ratnasingam et al. 2022 in preparation).

Figure 2 shows the time snapshots of vorticity for all the models considered in this work.

Zoom In Zoom Out Reset image size

The thermal diffusivity (κ) and viscosity (ν) are set to constant values shown in Table. 1. Our objective was to attain the highest Reynolds number,

$$\begin{eqnarray}&&\mathrm{Re}\,=\,\displaystyle \frac{{v}_{\mathrm{rms}}L}{\nu },\end{eqnarray} \tag{ 4 }$$

for each model studied. Here v_rms is the root mean squared velocity averaged over the convection zone and L is the characteristic length scale (the radial extent of the convection zone in this case). We also aimed to attain convective velocity ratios for different masses from our simulations which are comparable to those predicted by the mixing-length theory (MLT) velocities in MESA for any given age. (We chose to compare across masses as the uncertainties in the MLT velocities in MESA are larger across ages.) We, therefore, chose the values such that the above are satisfied while maintaining numerical stability. The simulation is run for an average time of approximately, 2.7 × 10⁷ s for the models studied.

The simulation domain extends up to 90% of the total stellar radius for all the 3 and 7 M_⊙ models and up to 80% of the total stellar radius for all the 20 M_⊙ models. The cutoff radius is determined such that the density does not vary more than ∼6 orders of magnitude as we move from the inner to the outer boundary to maintain numerical stability. As evident from Figure 1, this is satisfied only until 0.9R_⋆ for the 3 and 7 M_⊙ models and up to 0.8R_⋆ for the 20 M_⊙ models. We, therefore, chose these respective domains for our simulations. For our analysis, we chose the velocity data from a time, 4 × 10⁶ s for all the 3 and 7 M_⊙ models and from 6 × 10⁶ s ⁵ for all the 20 M_⊙ models. The total time window considered is given in terms of convective turnover times in Table. 1. We then used this velocity data which is saved at a regular time interval to study the mixing by IGWs in stellar interiors.

To measure the mixing by IGWs in stellar interiors, we introduce ${ \mathcal N }$ tracer particles into our simulation and track them over time T. This procedure is carried out in post-processing and we use the relations given by RM17 to calculate the diffusion coefficient D(r, τ) at radius r, for a given time difference, τ.

$$\begin{eqnarray}&&D(r,\tau )=\displaystyle \frac{Q(r,\tau )}{2\tau n(r,\tau )}-\displaystyle \frac{P{\left(r,\tau \right)}^{2}}{2\tau n{\left(r,\tau \right)}^{2}}.\end{eqnarray} \tag{ 5 }$$

where n(r, τ), is the number of sub-trajectories starting at r with a duration of τ, P(r, τ) is the sum of the lengths of the sub-trajectories and Q(r, τ) is the sum of the square of these sub-trajectories. More details on the calculation of n, P, Q, and D(r, τ) can be found in RM17.

We then plotted the diffusion coefficients calculated using the above equation as a function of radius. We varied different parameters such as the number of particles, radial grid size, and velocity time steps in our simulation to check for numerical convergence and obtained the diffusion coefficient for each set of parameters. Figure 3 shows the radial diffusion profiles with (a) different velocity time steps Δt, (b) different radial grid sizes, and (c) different numbers of particles ${ \mathcal N }$, for a 7 M_⊙ ZAMS model. In each case, the profiles appear converged with respect to the variables changed. We extended the same treatment to all the other stellar models and found the radial diffusion profile to be already converged for a time resolution of 1000 s, radial resolution of 500, and particle resolution of 10,000. This is consistent with RM17, who found the mixing profile robust to parameter changes. Henceforth, we choose these values for all our further analysis and obtain the radial diffusion profile at a time difference ⁶ of 1.3 × 10⁷s for all three masses.

Zoom In Zoom Out Reset image size

Figure 4 shows the radial diffusion profiles at ZAMS (blue), midMS (green), and TAMS (pink) for the 3, 7, and 20 M_⊙ models. In general, we find that the mixing profile varies significantly across the ages for all the masses studied. The amplitude of the diffusion profile depends on several factors such as the wave driving, density stratification, thermal diffusivity, the Brunt–Väisälä frequency, and the geometric effect. Among these, density stratification and thermal diffusivity dominate in both ZAMS and midMS stars; hence, the diffusion coefficient increases with stellar radius following the density stratification modulated by the thermal diffusion.

Zoom In Zoom Out Reset image size

However, the mixing profile at TAMS shows an increasing trend from the convective–radiative interface followed by a decreasing trend toward the surface, which is considerably different from that of the younger stars. To understand this behavior, we looked at the following second-order differential equation for wave propagation,

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

which is obtained by reducing the linearized 2D hydrodynamic equations considering no thermal or viscous diffusion with $\alpha ={v}_{r}{\bar{\rho }}^{\tfrac{1}{2}}{r}^{\tfrac{3}{2}}$ and m is the 2D Fourier basis wavenumber. In younger stars, the first two terms dominate over most of their radii. However, in older stars, the density term,

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}\left[\displaystyle \frac{\partial {h}_{\rho }}{\partial r}-{h}_{\rho }^{2}+\displaystyle \frac{{h}_{\rho }}{r}\right],\end{eqnarray} \tag{ 7 }$$

becomes dominant at larger radii.

We define turning point as the radius at which the ratio of the oscillatory term,

$$\begin{eqnarray}&&\left[\displaystyle \frac{{N}^{2}}{{\omega }^{2}}-1\right]\displaystyle \frac{{m}^{2}}{{r}^{2}},\end{eqnarray} \tag{ 8 }$$

to that of the density term $\left(\tfrac{1}{2}\left[\tfrac{\partial {h}_{\rho }}{\partial r}-{h}_{\rho }^{2}+\tfrac{{h}_{\rho }}{r}\right]\right)$ is equal to 1 (Ratnasingam et al. 2020) and it is located at a lower fraction of the total radius in older stars. Waves lose their wave-like behavior when this ratio is less than 1, and as a result, they experience extra damping. This leads to waves becoming evanescent toward the surface and leads to the diffusion profile we see in Figure 4. The shaded region in the plot indicates the location of turning points for the TAMS models using the range of frequencies (4–14 μ Hz) determined in Section 3.4.

Another factor contributing to the diffusion profiles of TAMS models is the Brunt–Väisälä frequency spike near the convective–radiative interface left behind by a steep composition gradient as the star evolves from ZAMS (Figure 1). This peak damps many waves near the convective–radiative interface and traps the higher frequencies resulting in fewer waves contributing to mixing, thereby leading to a decrease in the diffusion coefficient. In general, TAMS models have a significantly reduced diffusion coefficient than the ZAMS and midMS models.

Figure 5 shows the diffusion profiles as a function of radius for 3 M_⊙ (green), 7M_⊙ (pink), and 20 M_⊙ (blue) at different ages. As evident from the figure, the mixing is stronger in massive stars, and this can be accounted for by the stronger convective driving in massive stars. Another factor contributing to the increase in the diffusion coefficient with mass is the Brunt–Väisälä frequency. As noted in Figure 1, the Brunt–Väisälä frequency decreases with increasing mass. This leads to lower damping of waves in stars with higher mass and hence results in stronger mixing than lower mass stars.

Zoom In Zoom Out Reset image size

RM17 determined the diffusion profile in the radiation zone to be set by the square of the wave amplitude given by

$$\begin{eqnarray}&&D={{Av}}_{\mathrm{wave}}^{2}(\omega ,l,r),\end{eqnarray} \tag{ 9 }$$

where the coefficient A was speculated to depend on the diffusivities, rotation, and dimensionality and was found to be ∼1 s by RM17.

To test this theory against the new data, we calculated the wave amplitudes using the linear theory given by Ratnasingam et al. (2018) where they solved the anelastic, linearized hydrodynamic equations within the WKB approximation without considering the rotational, thermal, and viscous effects. The equation is given by

$$\begin{eqnarray}\begin{array}{rcl}{v}_{{\rm{wave}}}(\omega ,l,r) & = & {v}_{0}(\omega ,l,r){\left(\displaystyle \frac{\rho }{{\rho }_{0}}\right)}^{-\tfrac{1}{2}}{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{-1}\\ & & {\left(\displaystyle \frac{{N}^{2}-{\omega }^{2}}{{N}_{0}^{2}-{\omega }^{2}}\right)}^{-\tfrac{1}{4}}{e}^{-\frac{{ \mathcal T }}{2}},\end{array}\end{eqnarray} \tag{ 10 }$$

where ρ₀, r₀, and N₀ are the density, radius, and the Brunt–Väisälä frequency at the initial reference point. v₀(ω,l,r) is the initial wave amplitude for a given frequency and wavenumber. The radial velocity is related to the tangential velocity by the following relation:

$$\begin{eqnarray}&&\displaystyle \frac{{v}_{\theta }}{{v}_{r}}={\left(\displaystyle \frac{{N}^{2}}{{\omega }^{2}}-1\right)}^{\tfrac{1}{2}}.\end{eqnarray} \tag{ 11 }$$

The above can be approximated to $\tfrac{N}{\omega }$ for frequencies much smaller than the Brunt–Väisälä frequency. We, therefore, assume the wave amplitude to be equal to the rms velocity at the reference point averaged over the time window considered (see Section 2.1) multiplied by a factor of $\tfrac{\omega }{N}$ to obtain the theoretical radial diffusion profiles. The term $\exp (-\frac{{ \mathcal T }}{2})$ in the above equation is the attenuation factor multiplied to the wave amplitude as the radiative damping is considered within the quasi-adiabatic limit (Zahn et al. 1997). The damping coefficient ${ \mathcal T }$ according to Kumar et al. (1999) is given as

$$\begin{eqnarray}&&{ \mathcal T }(\omega ,l,r)={\int }_{{r}_{0}}^{r}\displaystyle \frac{16\sigma {T}^{3}}{3{\rho }^{2}\kappa {c}_{p}}\left(\displaystyle \frac{{\left(l\left(l+1\right)\right)}^{\tfrac{3}{2}}{N}^{3}}{{r}^{3}{\omega }^{4}}\right)\left(1-\displaystyle \frac{{\omega }^{2}}{{N}^{2}}\right){dr},\end{eqnarray} \tag{ 12 }$$

where σ, T, c_p , κ, and l are the Stefan–Boltzmann constant, temperature, specific heat capacity at constant pressure, opacity, and wavenumber. RM17 studied a 3 M_⊙ midMS model considering up to 70% of the total stellar radius without taking into account the effect of the Brunt–Väisälä frequency on the wave amplitude. We extend their analysis to stellar models presented in this work by considering a modified relation to calculate the wave amplitude (Equation (10)) based on the linear theory given by Ratnasingam et al. (2018).

In order to incorporate this mixing into 1D stellar evolution codes, we needed to determine the dominant waves contributing to the mixing profiles, so we launched waves from a radius outside the convective–radiative interface and calculated the theoretical diffusion coefficient using Equations (9)–(12) for a set of frequencies (1–40 μHz) and wavenumbers (1–10) for each stellar model. The initial reference points where the waves are launched are highly influenced by the overshooting motions. We, therefore, initially determined the overshooting depth for each model by fitting the following Gaussian function to the mixing profile near the convective–radiative interface,

$$\begin{eqnarray}&&D={D}_{0}\exp \left[\displaystyle \frac{{\left(x-{cr}+{f}_{0}{H}_{p}\right)}^{2}}{2{\left({f}_{1}{H}_{p}\right)}^{2}}\right],\end{eqnarray} \tag{ 13 }$$

where D₀ is the amplitude, cr is the radius of the convection zone, H_p is the pressure scale height, and f₀ and f₁ are the overshooting parameters. We measured the overshooting depth as the radius up to which we could get a good fit of the Gaussian profile. Figure 6 shows the Gaussian fit along with the mixing profile near the convective–radiative interface for a 7M_⊙ ZAMS model as an example. Here, we obtained the mixing profile (blue) near the convective–radiative interface from the simulation and then fitted this profile with the Gaussian function (pink) given by Equation (13). The overshooting depth determined is indicated by the pink vertical dashed line. We repeated this procedure for all the models and the overshooting depth for the different models found are given in Table 2. We then launched the waves at a radius equal to 1.3 ⁷ times the overshooting depth in all the models studied. More details on the overshooting profiles along with the comparison of the overshooting depth determined from the theory and simulation will be discussed in an upcoming paper.

Zoom In Zoom Out Reset image size

In order to implement a prescription for wave mixing into a 1D code, it would be helpful to narrow down the range of frequencies and wavenumbers contributing to the mixing profile. We, therefore, compared the theoretical profiles calculated using Equations (9)–(12) at different frequencies and wavenumbers with our simulation diffusion profiles for all the models to determine if a simpler prescription, with limited waves, could be found. We observed that the waves with higher wavenumbers provide a smaller contribution toward the diffusion coefficient since they experience higher thermal damping as expected from Equation (12), and therefore, we discard their contributions and consider only l = 1 waves in our further calculations. We show the theoretical diffusion profiles corresponding to the dominant frequencies that provide the best fit to the profiles from the simulations in Figure 7. The figure also shows the theoretical profiles at different frequencies as an example (dotted profiles in Figure 7). In general, the mixing profiles obtained from the linear theory corresponding to the dominant frequencies are in good agreement with the simulations diffusion profiles for all ZAMS (Figures 7(a)–(c)) and midMS (Figures 7(d)–(f)) models.

Zoom In Zoom Out Reset image size

We found the dominant frequencies for the different models studied to be in the range of 4–9 μHz at wavenumber l = 1; with the lower frequencies being dominant for massive stars. This can be attributed to the thermal damping given by Equation (12). Even though the lower frequencies are generated with larger amplitude near the convective–radiative interface, these waves are highly affected by thermal damping. As noted in Figure 1, the average Brunt–Väisälä frequency decreases with mass, causing the effect of damping to be less dominant in massive stars and thereby resulting in lower dominant frequencies. We also noted that the waves experience extra damping near the convective–radiative interface in all the midMS models because of the Brunt–Väisälä frequency peak as noted in Figure 1. The theoretical profiles do not match the numerical simulations well at the surface, nor the interface. At the surface, this is because the radial velocity is forced to zero in our hydrodynamic simulations, causing the amplitude of our numerical profile to drop near the surface, which is not necessarily physical. At the interface, convective overshooting means the region is not totally dominated by waves (Rogers et al. 2006). In our models, the deviation of the theoretical diffusion profile from that of the simulation near the interface, particularly in the case of ZAMS can be explained by these overshooting motions, as they influence the determination of the initial reference point at which the waves are to be launched.

We followed the same procedure to calculate the theoretical diffusion profiles for the TAMS models. We found that none of the theoretical profiles could explain the profile from the numerical simulation accurately. Figure 8 shows the best match that we could get from our analysis for the different TAMS models (pink dotted–dashed profiles). The pink vertical dashed line represents the turning point for the frequencies shown. We can clearly see that the profiles start to diverge as they approach the turning point beyond which the waves become evanescent. As stated previously, the waves lose their wave-like behavior at a lower fraction of the total stellar radius because of the density term in Equation (6) being dominant for a major fraction of the radius in older stars. Therefore, the WKB approximation applied to the equation for wave propagation (Equation (6)) is not valid in the TAMS models and hence the theoretical description outlined in RM17 and Section 3.3 is not valid in older stars. More analyses need to be carried out in the future for a better understanding of the diffusion profiles of the TAMS models. However, at this stage, it appears that a substantially lower diffusion coefficient, which is virtually flat throughout the radiation zone is a decent approximation of wave mixing at the late stages of stellar evolution.

Zoom In Zoom Out Reset image size

We also computed the kinetic energy spectrum, ${\hat{E}}_{\mathrm{kin}}$,

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

where $\hat{{v}_{r}}$ and $\hat{{v}_{\theta }}$ are the Fourier transforms of radial and tangential velocities. Figure 9 shows the frequency distribution of the kinetic energy for l = 1 at three different radii where the top panels (a)–(c) present the spectra at 1.3 times the overshoot depth (radius at which the waves are launched), the middle panels (d)–(f) at 0.3R_⋆, and the bottom panels (g)–(i) at 0.6R_⋆ for all the models studied. The vertical dashed lines indicate the dominant frequencies shown in Figures 7 and 8. We note that the dominant frequencies lie close to where the kinetic energy density peaks, thereby justifying our choice of frequencies. This is more evident from the middle and bottom panels of Figure 9. We are not concerned with the TAMS models as we have already shown that the linear theory does not explain their mixing profiles accurately.

Zoom In Zoom Out Reset image size

We conducted a number of tests to determine the value of the parameter "A" from our simulations. We found that the simulation should be run for longer convective turnover times for the correct estimate of this parameter. Regardless of how long the simulations were run, the mixing profiles maintained the same trends discussed in Sections 3.1 and 3.2. But, we noted that the amplitudes of the profiles decrease further and converge at higher time differences. This is shown in Figure 10(a) for the 7 M_⊙ ZAMS model run for 87 convective turnover times. We see that the amplitude of the mixing profile converges at ∼τ = 3 × 10⁷ s. We then compared this profile to the square of wave amplitude (here the amplitude is summed over a range of frequencies (1–40 μHz) and wavenumbers (1–20) ⁸ ) obtained from the simulation (pink in Figure 10(b)). As evident from the figure, both the profiles match to a good extent, giving us the value of the parameter "A" ∼1. We expect all the models to show similar behavior once the simulation is run for longer convective turnover times. Our finding is in agreement with that of RM17. A better estimate of this parameter can only be obtained by the comparison of our theoretical prescription with that of the observations.

Zoom In Zoom Out Reset image size