© The Authors 2023

1. Introduction

Given the progress made in the field of asteroseismology (the study of pulsations and oscillations in stars apart from the Sun), it is clearly an auspicious time for theoretical research on stellar interiors to take the centre stage. Direct observables from stars, such as brightness variations from photometric data and stellar light spectrum variations from spectroscopic data, can be analysed further with asteroseismological techniques to quantify fluid oscillations in stellar interiors. These oscillations take the form of p-modes (pressure modes) and g-modes (gravity modes). Fortunately, hydrodynamical simulations can provide information on how these g-modes and p-modes are formed and how they propagate through a star. With a rapidly increasing sample of stars available for asteroseismology from space missions, there is a strong call for studies of the effect of stellar mass and age on surface observables.

Thanks to helioseismology (the study of pulsations and oscillations in the Sun) providing a huge amount of information on solar-like stars, we turn to more massive stars, which have not been studied as much. One particular characteristic of more massive stars is the presence of a convective core below a layer of radiative envelope, instead of a radiative core encompassed by a convective envelope in less massive stars (e.g. the Sun). The radiative layer is always stably stratified and subadiabatic, which allows for the propagation of outwardly travelling internal gravity waves (IGWs), generated mostly at the radiative-convective boundary.

It is important to study these waves in massive stars, as they provide information about the top of the convective cores. Currently, one of the more important unanswered questions in this field is the role of convection in generating these waves. Kumar et al. (1999) suggests that the bulk material inside the convection zone generates waves, which then tunnel through and manifest themselves at the top of the zone. This notion was developed further in Lecoanet & Quataert (2013) and a theoretical spectrum with a cut-off point at the convective turnover frequency was introduced. A more numerical approach was done in Rogers et al. (2013) and the results from the two-dimensional (2D) equatorial annulus simulation of a three-solar mass (M_⊙) star showed a shallower generation spectrum. This work was improved upon with three-dimensional (3D) simulations in Edelmann et al. (2019), where the authors arrived at very similar results. More recently, Pinçon et al. (2016) theoretically showed how plume models can lead to observations of more than one slope in the frequency spectrum of IGWs.

In addition to the issue of uncertainty in the generation spectrum, IGW propagation in the radiation zone is also a subject of widespread study. One of the very first studies on this was carried out by Press (1981), who showed that using a locally Bousinesseq and globally anelastic assumption, the linear propagation of IGWs can be represented with an analytical model. This work took into consideration the effect of pseudo-momentum conservation, geometry, and radiative damping on IGW propagation. Zahn et al. (1997) revisited this work and used the same expression for both viscous and radiative damping interchangeably, depending on which effect was more dominant for a given case. Ratnasingam et al. (2020) showed that the 2D and 3D linear models differ by a factor of the square root of the radial coordinate. Apart from linear analyses, non-linear simulations (Edelmann et al. 2019; Horst et al. 2020) have also been done to investigate the propagation of IGWs.

Until now, numerical simulations of a convective core underlying a radiative envelope were done for mainly a 3 M_⊙ star, as done in Rogers et al. (2012) and Edelmann et al. (2019). With the assumption that the IGW spectrum can be scaled to represent that from stars of different masses and ages, Aerts & Rogers (2015) showed that there was a match between stellar brightness variations and the IGW signatures predicted by numerical simulations. However, the generation spectrum and, hence, the surface IGW spectrum will be different for stars of different ages and masses for a range of reasons, such as stellar size, where heavier and older stars tend to be larger, and convective velocities, where heavier and older stars have larger convective velocities.

In Sect. 2, we present our numerical setup, which describes the numerical methods, equations, and the code used in this study. Section 3 explains the properties of the stellar interiors and how they were modified from models produced by the 1D stellar evolution code named Modules for Stellar Astrophysics or MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019) to ensure a stable running of the hydrodynamical simulations. Section 4 describes our results from simulations of different stellar models. Section 5 compares these results with some photometric data from a recent TESS survey. The final section of the paper presents a summary and our conclusions.

2. Numerical setup

We used a pseudo-spectral method to solve the hydrodynamical momentum and energy equations, namely, Eqs. (1)–(3) in the anelastic approximation (Rogers & Glatzmaier 2005) within the geometry of a 2D equatorial slice. An analysis of how these equations differ from the fully compressible Navier–Stokes equations and variations of the anelastic approximation was done in Brown et al. (2012), which referred to the following equations as the Rogers–Glatzmaier (RG) equations. These are:

$$\begin{array}{cc}& \mathrm{\nabla }·\overline{\rho }\mathit{v}=0,\hfill \end{array}$$(1)

$$\begin{array}{cc}& \frac{\partial \mathit{v}}{\partial t}+(\mathit{v}·\mathrm{\nabla })\mathit{v}=-\mathrm{\nabla }\left(\frac{P}{\overline{\rho }}\right)-C\overline{g}\mathit{r}+2(\mathit{v}\times \widehat{\mathit{z}}\mathrm{\Omega })\hfill \\ \hfill & +\overline{\nu }({\mathrm{\nabla }}^2\mathit{v}+\frac{1}{3}\mathrm{\nabla }(\mathrm{\nabla }·\mathit{v})),\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\partial T}{\partial t}+(\mathit{v}·\mathrm{\nabla })T=(\gamma -1)T{h}_{\rho }{v}_{\mathrm{r}}\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill & -v_{\mathrm{r}}(\frac{\mathrm{d}\overline{T}}{\mathrm{d}r}-(\gamma -1)\overline{T}{h}_{\rho })\hfill \\ \hfill & +\frac{1}{c_v\overline{\rho }}\mathrm{\nabla }·(c_{\mathrm{p}}\overline{K}\overline{\rho }\mathrm{\nabla }T)\hfill \\ \hfill & +\frac{1}{c_v\overline{\rho }}\mathrm{\nabla }·(c_{\mathrm{p}}\overline{K}\overline{\rho }\mathrm{\nabla }\overline{T})\hfill \end{array}$$(4)

where ρ, T, and P are the perturbation density, temperature, and pressure, respectively. Any quantity with an over-line or bar represents the background state variables, which vary only with radius. The radial velocity, v_r, and the tangential velocity, v_θ, together form the velocity vector, v. The negative inverse density scale height is represented by h_ρ = d$(ln\overline{\rho })$/dr. The rotational angular velocity and viscosity are represented by Ω and $\overline{\nu }$. The adiabatic index, γ is set to be 5/3. The specific heat capacity at constant volume and specific heat capacity at constant pressure are represented by c_v and, c_p respectively. The reference state gravity is represented by $\overline{g}$. The co-density (Braginsky & Roberts 1995; Rogers & Glatzmaier 2005), labelled as C, is

$$\begin{array}{c}\hfill C=-\frac{1}{\overline{T}}(T+\frac{1}{\overline{g}\overline{\rho }}\frac{\partial \overline{T}}{\partial r}P)·\end{array}$$(5)

In the radial direction, the RG equations are solved using a uniform finite difference scheme, specifically: the central difference method (e.g. Grossmann et al. 2007). In the horizontal direction, the spectral solution method was applied, where stellar variables are expanded in sines and cosines. The non-linear terms of the RG equations are integrated in time using the Adams-Bashforth method. The linear terms are integrated in time using the Crank-Nicolson method.

The temperature perturbation, T, and radial velocity, v_r, at the top boundary of the simulation domain were set to zero. For the horizontal velocities, v_θ, at the top boundary, the following stress-free boundary condition was imposed as ∂v_θ/∂r = 0.

2.1. Setup of the convection and radiation zones

The 1D stellar evolution code, MESA, was used to generate models of various stars in the intermediate-mass regime at the middle of the main sequence (midMs), when the central hydrogen mass fraction becomes approximately 0.35 These models were then set as the reference state for these simulations. The stellar metallicity, Z, was set to be equal to the solar value of Z = 0.02. The mixing-length parameter was set to 1.8 and the convective overshoot profile was set to exponential (see Eq. (2) in Pedersen et al. 2018 for details). Due to limitations in the numerical solver that we utilised, the maximum density variation between the bottom and top of the simulation domain had to be capped at approximately six orders of magnitude (see Fig. 1a). This was found to be within 90% of the total stellar radius (as defined in MESA), for all of our models, the MESA models have been cropped. We note that even with trimming, stars with larger mass have larger radii, which causes the radial resolution to decrease with increasing mass as we are limited by the total number of radial grids in our simulations. Table 1 summarises this information for all the models we simulated.

Fig. 1. Density (a) and Brunt–Väisälä frequency (b) profiles as functions of stellar radius, in units of total solar radius. All the models are cut at 90% the total stellar radius to ensure numerical stability.

We study intermediate-mass stars, specifically stars with a convective core and radiative envelope. In our simulations, the convection and radiation zones are set using the subadiabaticity term in Eq. (3), as given below:

$$\begin{array}{c}\hfill -(\frac{\mathrm{d}\overline{T}}{\mathrm{d}r}-(\gamma -1)\overline{T}{h}_{\rho })·\end{array}$$

This term represents the difference between the stellar temperature gradient and the adiabatic temperature gradient, where a positive value is necessary for convection and conversely, a negative value indicates a convectively stable region. It is derived with the ideal gas assumption and can be used to calculate the Brunt–Väisälä frequency using:

$$\begin{array}{c}\hfill N^2=\frac{\overline{g}}{\overline{T}}(\frac{\mathrm{d}\overline{T}}{\mathrm{d}\overline{r}}-(\gamma -1)\overline{T}{h}_{\rho })·\end{array}$$(6)

However, stellar models from MESA are produced with a more realistic equation of state, which disagree with our expression for the subadiabaticity. Thus, to better represent the subadiabaticity in the radiation zone, we calculated this term from the Brunt–Väisälä frequency profile for all of our MESA models, by rearranging Eq. (6). However, in the convection zone, MESA provides negative Brunt–Väisälä frequencies which are physically inaccurate. Thus, we find a value between 1 × 10⁻⁶ and 1 × 10⁻⁷ that leads to the desired convective velocities once the convection in this region reaches a steady-state evolution. The realistic values of subadiabaticity in stellar convection zones have yet to be determined by the larger stellar research community.

Figure 1b shows the Brunt–Väisälä frequencies for MESA stellar models with different masses. Generally, we see a decrease in the average Brunt–Väisälä frequency with increasing stellar mass. We also see a significant peak at the beginning of each profile, which is an artefact of the retracting convective core as a star ages.

3. Properties of the stellar interior

3.1. Core convection

The 1D stellar evolution code, MESA, simulates the convection zones in stars using a formulation known as mixing-length theory (MLT). Although this method is computationally inexpensive, there is much debate on the accuracy of this technique among the stellar research community. In our simulations, the average convective velocity for any particular stellar model was forced to be higher than their corresponding MLT velocity given by MESA (refer to Table 2). This was done to counterbalance the effect of the unrealistically large viscosities and thermal diffusivities used in our simulations, which cause stronger damping on waves than expected in a real star. Given that we are generally more concerned with wave amplitudes at the surface, the higher convective velocities allow IGWs to be better represented, as they are generated with higher amplitudes and less affected by thermal diffusion and viscosity. Additionally, it is not certain whether the MLT velocities accurately represent convective velocities in stars (Jones et al. 2017).

One of the key differences in core convection between different stellar masses is the magnitude of convective velocities. With increasing stellar mass, convective velocities are larger due to higher temperatures, as shown in the last column of Table 2, which summarises different properties of the models used in our investigation. The error values indicate the general variation of MLT velocities in the convection zone¹. Convective velocities are also expected to increase with stellar age, which will be investigated in Rogers et al. (in prep). Mainly due to the larger convective velocities in more massive stars, several convective parameters such as the Reynolds and Rayleigh numbers were found to be larger too, as seen in the 4th and 6th columns in Table 2. The Reynolds number, Re, is defined as

$$\begin{array}{c}\hfill \mathrm{Re}=\frac{\mathit{uL}}{\nu },\end{array}$$(7)

where ν is the kinematic viscosity, u is the characteristic velocity in the convection zone, and L is the characteristic length scale. It is usually defined as the radial extent of the convection zone. This quantity describes the ratio of the inertial forces to the viscous forces in the system. Low values signify laminar flow, while high values show turbulent flow. For more massive stars, larger Re indicates stronger and more turbulent convection.

For the real intermediate-mass stars studied in this paper, we expect the average Prandtl number $\mathrm{Pr}=\nu /\overline{K}$, to be somewhere close to 1 × 10⁻⁶, which is much smaller than the values we use in our simulation (Garaud et al. 2015). The work done in Garaud et al. (2015) further shows that Pr is expected to decrease with increasing stellar mass and with increasing radial distance from the stellar core within the star (Fig. 7 in Garaud et al. 2015). Generally, with increasing age, the stellar interior profiles of an intermediate-mass star follow closely those of higher mass stars, so we expect Pr to decrease as a star ages and also with radial distance from the stellar centre. In our models, we show in Fig. 2 that Pr varies approximately from 10 at the bottom boundary to 0.01 at the top boundary. Unfortunately, this mismatch between realistic and simulated Pr is a computational limitation faced by many such studies, especially when turbulent convection is involved.

Fig. 2. Prantl number profiles used in the simulations of all the stellar models used in this work.

3.2. The radiation zone

The radiation zone is a convectively stable region, where radiative diffusion represents the main form of heat transport. The stratified density in this region allows the propagation of IGWs, generated mainly in the convection-radiative boundaries. As these waves travel through the radiation zone, radiative diffusion causes these waves to be damped, which primarily depends on the thermal diffusivity and the Brunt–Väisälä frequency profiles. For linear IGWs, the effect of thermal diffusion (Press 1981; Kumar et al. 1999) on wave amplitudes can be written as:

$$\begin{array}{c}\hfill A_{\mathrm{IGW}}\propto A_0{e}^{-\tau }\end{array}$$(8)

where

$$\begin{array}{c}\hfill \tau ={\displaystyle {\int }_{r_{\mathrm{generation}}}^r}\overline{K}\left(\frac{m^3{N}^3}{r^3{\omega }^4}\right){(1-\frac{{\omega }^2}{N^2})}^{1/2}\mathrm{d}{r}^{\prime }.\end{array}$$(9)

The damping opacity, τ, is a function of radius, r, the wave frequency, ω, and wavenumber, m.

The thermal diffusivity is a function of stellar radius and increases by at least four orders of magnitude in the radiation zone (see solid lines in Fig. 3). However, to ensure numerical stability, we have to boost the thermal diffusivities in our simulations. Although we cannot reproduce the realistic thermal diffusivity profile in stars, we mimic the increase using the profile shown below:

$$\begin{array}{c}\hfill \overline{K}=K{\rho }^{-0.5},\end{array}$$(10)

Fig. 3. Thermal diffusivities (top panel) as a function of stellar radius and the ratio of linear IGW amplitudes obtained with real stellar thermal diffusivities to those obtained with $\overline{K}$ profiles used in our main simulations (bottom panel) as a function of wave frequencies at the top of the radiation zone. The dotted lines in the top panel indicate the varying thermal diffusivity profiles used for our main simulations, whilst the solid lines indicate the realistic thermal diffusivities from MESA. In the bottom panel, the wavenumbers for all the waves have been set to 1. The initial radius for these IGWs was chosen at random from the middle of the radiation zone and the initial amplitude was set to unity.

where K is a constant, chosen to ensure numerical stability. Our $\overline{K}$ values are generally larger than in real stars except at the surface, but this is compensated by our larger convective velocities (see previous section for details). This thermal diffusivity profile was applied on five different stellar models and the kinematic viscosity was chosen to be a constant (see Table 2 for details on the thermal diffusivity and viscosity profiles). The results of this investigation are elaborated further in Sect. 4.

4. Main results

To investigate the effect of our artificial thermal diffusion profile compared to the more realistic ones from MESA, we started by conducting a small study using the IGW linear theory. In Fig. 3, we show two plots, where the thermal diffusivities used in these simulations (top panel) are higher than the realistic values from MESA. From the bottom panel, which shows the ratio of linear IGW amplitudes for the thermal diffusivity from MESA to those for the artificial thermal diffusivity, we see that at lower frequencies, the artificial $\overline{K}$ predicts smaller IGW amplitudes than those from MESA’s $\overline{K}$. We find that for any wave frequency above 8 μHz, simulations of all the stellar models should predict wave amplitudes with higher accuracy. However, below this frequency, we see large variations in the ratios which increase with increasing stellar mass.

Moving on to the simulations, we show the average v_rms in the convection zone for all of our stellar models as a function of time in Fig. 4. Generally, we see that increasing the stellar mass leads to increasing convective velocities. Looking at individual convective velocity profiles, we see the following trend: a rapid increase in fluid velocity, indicating the onset of convection, followed by a gradual decrease and settling into a steady-state evolution. We chose a time interval of 6 × 10⁶ s (only when the convective core is in a steady-state evolution) to perform all of our time-series analysis. Figure 5 shows the log of radial velocities from our simulations as a function of stellar radius and wave frequency for a randomly chosen ray or angle. All five spectra exhibit stationary modes, which are more distinctly recognisable at higher frequencies. The lower masses exhibit more modes, mainly because of the generally larger Brunt–Väisälä frequency (white, solid lines) for models with lower stellar masses. In the convection zone, velocities are larger at smaller frequencies, which shows that most of the energy in the convection zone is concentrated on larger timescale motions (e.g. 10 to 30 μHz, which corresponds to approximately motions on timescales between 0.8 and 3 days).

Fig. 4. Average vrms in the convection zones as a function of time of all the models in Table 2.

Taking slices of each plot in Fig. 5 for all five stellar models at r = 0.5 r_CZ, where r_CZ is the extent of the convection zone, we obtain the spectra of motions as shown in Fig. 6 (radial velocities in the top panel, tangential velocities in the middle panel, and temperatures in the bottom panel). We find two main observations, with the first one being that for all the models, there is a trend of high amplitudes at very low frequencies, which decreases rapidly at mid-range frequencies. Second, the homogenous mixing of convection has allowed the tangential and radial velocities of each stellar model to become similar and, thus, their slopes (indicated by the dashed lines) end up becoming similar as well. To confirm this, we show the ratio of the tangential velocity to the radial velocity in the top panel of Fig. 7, which shows oscillatory behaviour around 1.

Fig. 5. Radial velocities as functions of stellar radius and frequency for different stellar masses, all at midMs. The solid, white line shows the Brunt–Väisälä frequency profile for each stellar model, respectively.

Fig. 6. Radial velocities (top panels), tangential velocities (middle panels), and temperatures (bottom panels) as a function of frequency for different stellar masses in the middle of the convection zone for a randomly chosen radial ray. The dashed lines show the line fits to each profile at frequencies between 8 μHz and 30 μHz.

The panels in row 2 and 3 of Fig. 7 show the ratio, v_θ/v_r at 0.1 H_P and 0.25 H_P, where H_P is the pressure scale height, from the top of the convection zone, respectively. For linear IGWs, we can use the dispersion relationship for IGWs to get

$$\begin{array}{c}\hfill \frac{v_{\theta }}{v_{\mathrm{r}}}={(\frac{N^2}{{\omega }^2}-1)}^{1/2}.\end{array}$$(11)

Fig. 7. Ratios of the tangential velocity (vθ) to the radial velocity (vr) as a function of frequency, at different radii throughout the stellar interior for different stellar masses. The black line represents the ratio of the Brunt–Väisälä frequency to the linear IGW frequency, at that particular radius. The first row shows vθ/vr inside the convection zone, the second and third rows show the ratios at 0.1 HP and 0.25 HP, respectively, and the final two rows show the ratios at a randomly chosen radius inside the radiation zone and close to the top of the simulation domain, respectively.

The term on the right in Eq. (11) is plotted as the solid, black line in Fig. 7. In general, we find that the simulations agree well with the linear theory. There are a few exceptions, such as the disagreements seen between the black line and the velocities at very low frequencies for some radii and masses. At such low frequencies, our simulations do not resolve IGWs completely, so we did not include these waves in our analysis.

We show the velocity and temperature perturbations as a function of frequency at a randomly chosen angle and 0.25 H_P from the top of the convection zone in Fig. 8. The radius of 0.25 H_P from the top of the convection zone chosen here are all inside the Brunt–Väisälä frequency peaks seen in Fig. 1b, allowing for stationary modes at higher wavenumbers to form within this peak. Thus, we observe stationary wave peaks at higher frequencies, f > 50 μHz, in almost all the spectra. As mentioned previously, wave amplitudes from our simulations match those calculated using linear wave theory with realistic thermal diffusivities up to 50% accuracy above 8 μHz. Thus, we chose a frequency range of 8 μHz to 30 μHz for the rest of the analysis. We fit these spectra with a simple straight line form, i.e. Y = mX + C, where X is the logged frequency and Y is the logged velocity or temperature perturbation, represented by the dashed lines in Fig. 8. The value m, is defined as the slope of these fits. Generally, we observe flatter slopes (>  − 1) for the radial velocities and steeper slopes for the tangential velocities and temperature perturbations (<  − 1). We also observe the temperature perturbation slopes to be more similar to the tangential velocity slopes, compared to the radial velocity slopes.

Fig. 8. Radial velocities (top panels), tangential velocities (middle panels), and temperatures (bottom panels) as a function of frequency for different stellar masses at 1/4th the pressure scale height from the top of the convection zone for a randomly chosen radial ray. The dashed lines show the line fits to each profile at frequencies between 8 μHz and 30 μHz.

To further check whether these motions are waves, we can manipulate the following linear, asymptotic relation for IGWs,

$$\begin{array}{c}\hfill \frac{v_{\theta }}{v_{\mathrm{r}}}=\frac{k_{\mathrm{r}}}{k_{\theta }},\end{array}$$(12)

where k_r and k_θ are the IGW radial and tangential wavenumbers, respectively. Clearly, this expression does not prove the presence of non-linear waves, so satisfying this relationship indicates the presence of only linear waves. Taking the log of both sides of Eq. (12), replacing k_r and k_θ with the appropriate expressions for IGW wavenumbers, and differentiating both sides in terms of ∂log₁₀ω, we obtain

$$\begin{array}{c}\hfill \frac{\partial {log}_{10}{v}_{\theta }}{\partial {log}_{10}\omega }=\frac{\partial {log}_{10}{v}_r}{\partial {log}_{10}\omega }-1.\end{array}$$(13)

The above expression indicates that for linear waves, the slopes of the radial velocity and tangential velocity should differ by 1, and this is what we observe approximately in Fig. 8. For example, the v_r slope for 3 M_⊙ is −0.94 and that for v_θ is −2.06. It is obvious that we should not expect the difference to be perfectly equal to 1, which is the case for all the stellar models, as we are fitting stochastic perturbations. However, the fact that the difference is close to 1 is reassuring.

Returning to Fig. 7, the bottom two panels show the ratio, v_θ/v_r, above the middle of the radiation zone and close to the top of the simulation domain. We can see that except at very low frequencies, the linear wave properties are still observed². These waves should, then, propagate through the radiation zone and manifest themselves close to the top of the domain, with properties that can be explained with linear wave theory. In Fig. 9, we show the velocity and temperature spectra close to the top of the simulation domain, where we see a general trend of smaller variations and lower amplitudes at lower frequencies. As the frequency increases, we observe larger variations and higher amplitudes at mid-range frequencies, but with a negative slope for all the stellar models.

Fig. 9. Radial velocities (top panels), tangential velocities (middle panels), and temperatures (bottom panels) as a function of frequency for different stellar masses at 0.82 Rstar for a randomly chosen radial ray. The dashed lines are line-fits to each profile at frequencies between 8 μHz and 30 μHz.

The spectra that we have presented so far are related to a radial ray, chosen at a random horizontal angle. Due to the stochastic behaviour of these perturbations, variations are expected for radial rays at different angles. Additionally, we chose a fixed frequency range (between 8 μHz and 30 μHz) to calculate the slopes for all the models, so we expect further variations for different frequency intervals. Thus, in Table 3, we present the average slopes with standard deviations for the radial velocities, tangential velocities, and temperatures close to the surface of each model (see the panels at the bottom of Fig. 7 for the exact radii). These slope averages and standard deviations are calculated for radial rays at different angles and different frequency intervals, all starting at 8 μHz. Equation (13) is approximately fulfilled for the slopes of the radial and tangential velocities. Finally, as seen in the figures shown in this section, the temperature slopes are still similar to the tangential velocity slopes. Comparing these slopes with the slopes at the convective-radiative interface (as shown in Fig. 8), s similar behaviour is observed for the averaged slopes, as shown in Table 4. This table suggests that during propagation in the stellar interior, IGWs do not experience any frequency dependent changes within the chosen frequency range here. However, the amplitude of the perturbations in Figs. 8 and 9 do show that there is a constant increase at all frequencies as IGWs propagate, which is consistent with the predictions of linear IGW theory (see A₀ in Eq. (8), where A₀ can be dependent on stellar radius).

5. Comparison with photometric data

The velocity and temperature perturbations seen at the top of our simulation domain are likely to be related to brightness variations seen in photometric data from stars. One of the first direct comparisons for intermediate-mass stars was done in Aerts & Rogers (2015). This work noted the similarity between the velocity perturbations from the simulations (Rogers et al. 2013) of a 3 M_⊙ star, with the brightness variations from young massive O-type stars observed with high-precision CoRoT space photometry. Many works from both the numerical and observation fronts followed, but a direct comparison between a set of refined photometric data and simulation results of different masses remains unexplored. To address this, we utilised some results of a recent survey of β Cephei stars in the Southern Hemisphere, conducted by a few co-authors of this paper, and whose results are yet to be published. In the survey, photometric data from the TESS satellite were used along with spectroscopy to classify and carry out statistical studies of β Cephei stars. In total, 486 stars were observed in low resolution (421 stars) and high resolution (65 stars) spectroscopic mode, and their spectroscopic parameters were derived, to complement the photometric aspect of the survey. Both photometrically and spectroscopically, the target selection was made through various means, which will be described in detail in Handler et al. (in prep.) and Chowdhury et al. (in prep.). To derive luminosities for these 486 spectroscopically observed stars, Gaia EDR3 parallaxes and TIC V magnitudes were utilized. All bolometric corrections were applied using interpolated values derived from Flower (1996). The extinction was calculated based on MWDUST python package (Bovy et al. 2016), using the combined (Marshall et al. 2006; Green 2015; Drimmel et al. 2003) reddening map provided by Bovy et al. (2016).

Having their T_eff and luminosity values derived, these 486 stars were then plotted in the HR diagram with the evolutionary tracks calculated using the Warsaw-New Jersey stellar evolution code (for details see, e.g. Pamyatnykh 1999, and references therein). To calculate the tracks, they used OPAL opacity tables, hydrogen abundance (X) as 0.74, heavy element abundance (Z) as 0.02, and initial rotational velocity as 150 km s⁻¹ (sample average). The masses were then extracted by comparing the position of stars in the HR diagram, with the evolutionary tracks for 61 different masses between 2.5 M_⊙ and 60 M_⊙ with varying step size, and nine different rotational velocities between 0 and 350 km s⁻¹ with a step size of 50 km s⁻¹. For this work, we utilize these derived masses.

Before comparing our simulation data with the observational data, we have to consider a few issues in choosing the sample. Our output data are 2D in space: radial and tangential. Thus, for a time-series analysis at any particular radius, there is freedom with respect to choosing the angle at which this is done. For the purposes of the analysis of a single radial ray (inside the convection zone, close to the interface, and at the top of the simulation domain), observations at a few angles and obtaining an error estimate would suffice (as shown in the previous section). However, for the purpose of comparing our simulation results with the photometric data of stars, a few approaches can be considered, but we chose to perform a weighted average of spectra over half the domain circumference. This is done to account for the fact that telescopes receive light from only half of any star that is observed. Additionally, we receive fewer photons from the edges of a star compared to the centre. To account for all of these effects, we performed a weighted average as shown in the equation below:

$$\begin{array}{c}\hfill f_{\mathrm{avg}}=\frac{1}{N}{\displaystyle \sum _i^N{f}_i-f_{i}b(1-sin\theta ),}\end{array}$$(14)

where N is the total number of angular grid points over half the domain circumference, b is a decimal number between 0 and 1, and θ is the angle of observation from the edge of half the simulation domain. We chose a simple linear expression to represent the photon count drop due to the geometrical effect mentioned above and a value of 0.6 for the coefficient b.

As mentioned earlier, the data from the survey done using the TESS satellite data cover a variety of stars at different masses. Thus, we compared our simulation results with the observational data in the following manner: 2 M_⊙–4 M_⊙ with 3 M_⊙ midMs, 4 M_⊙–6 M_⊙ with 5 M_⊙ midMs, 6 M_⊙–9 M_⊙ with 7 M_⊙ midMs, 9 M_⊙–11 M_⊙ with 10 M_⊙ midMs, and 11 M_⊙–15 M_⊙ with 13 M_⊙ midMs. The next step in the process was to decide which of the three simulation stellar parameters was best to compare with the brightness variations seen in the TESS photometric data. For this work, we make a choice of comparing the temperature perturbations, in units of Kelvin, with the brightness variations, in units of μmag, seen in the TESS photometric data. The temperature perturbations were chosen mainly because of direct relations between stellar luminosity and temperature (i.e. Stefan–Boltzmann law).

The main objective of this part of the study was to search for similarities and differences between the slopes of the brightness variations from the TESS photometric data and the slopes of temperature perturbations in our simulation results. To achieve this, we selected a fixed frequency window in logarithmic space and calculated the slopes of both the temperature perturbations and brightness variations with straight line fits within this frequency range. We then moved the window along the frequency axis in small increments and repeated the slope calculation. This scanning process was done from 8 μHz to 100 μHz, for different frequency intervals. The frequency range at which the difference between the temperature perturbation slope and the brightness variation slope is less than 5% (chosen arbitrarily) was selected as the best match. Due to our strict tolerance value, the number of stars that we found suitable for comparison reduced from the initial 486 in the survey to 22, with more higher mass stars fitting our criteria than lower mass stars.

The periodograms of the TESS data are over-plotted on the temperature perturbations from the simulations and presented in Appendix A. An important point to note here is that due to the difference in the units of the brightness variations and temperature perturbations, we adjusted the temperature perturbation amplitudes by using the y-intercepts of the straight-line fits done to these data. We also present the best matches for frequency slopes, the frequency range for these slopes, and the stellar masses for the simulation models and TESS stars in Table 5. The most negative frequency slope is found to be −1.87 and the least negative frequency slope is −0.74, with most of the slopes being more negative than −1.0. We found that for most of the TESS data, the frequency range at which the simulated temperatures agree best with the TESS surface brightness variations, lies above 10 μHz. For higher mass stars, however, this frequency range shifts to lower values (≥8 μHz), which, with the assumption that the TESS brightness variations are due to IGWs from core convection, could be caused by the decrease in the average Brunt–Väisälä frequency when the stellar mass is increased. In some cases (e.g. model TIC 412543692, 461021821, 466592679, 078636551 and 065516748), the periodograms even show that the amplitudes start to rise³ (orange lines), when simulated temperatures (blue lines) start to fall rapidly due to the Brunt–Väisälä frequency limit.

6. Discussion and conclusions

We ran 2D simulations to study the internal fluid evolution of five stellar models: 3 M_⊙, 5 M_⊙, 7 M_⊙, 10 M_⊙, and 13 M_⊙, which are all chosen at the middle of the main sequence (midMs). These stars all possess a convective core and a radiative envelope, which is cut off at 90% of the total stellar radius, to ensure numerical stability. Compared to previous studies with such anelastic hydrodynamical simulations, where constant thermal diffusivities (as a function of radius) were implemented, this study utilises varying thermal diffusivity profiles that are more realistic and allows for any additional effect due to the variations in thermal diffusivity to exist.

Before analysing our simulation results, we used a simple linear theory model for internal gravity waves (IGWs) to calculate the effect of thermal diffusivity profiles that we chose for our models. We found that in the case of linear IGWs, any perturbations above 8 μHz, will be approximately represented up to 50% accuracy, compared to the same perturbations with the real stellar thermal diffusivities. Thus, we chose 8 μHz as the lower limit for all of our analysis and we proceeded to study the radial velocity, tangential velocity, and temperature perturbations in frequency space. One of our main conclusions is that the perturbations generally show linear IGW-like behaviour throughout the radiation zone. Second, the frequency slopes at the domain surface are, on average, <  − 1 for the tangential velocity and temperature perturbations and >  − 1 for the radial velocity perturbation. This is also true at the convective-radiative interface. Third, there is no clear trend for these frequency slopes across the different stellar masses, although there is a clear trend of increasing amplitudes with increasing mass, mainly due to the higher convective velocities for higher stellar masses.

Although we did not delve extensively into the topic of IGW excitation in this work, we would like to offer a remark on this here. Generally, we observe that plumes are dominant in the convection zone (similar to Rogers et al. 2013). This is mainly seen from the IGW frequency slopes at generation, which are all shallower than −2 (between 0 and −2). These flatter slopes agree well with the theoretical work on plume forcing done in Pinçon et al. (2016), and the numerical work done in Edelmann et al. (2019) and Horst et al. (2020) on IGW generation. Earlier works such as that of Kumar et al. (1999), where the authors carried out a theoretical study of IGWs generated through Reynolds’ stresses inside the convection zone and then tunnel through predicted slopes steeper than −3, whilst Lecoanet & Quataert (2013) predicted slopes steeper than −4. Both these works on IGW generation through Reynolds’ stress predict much steeper slopes than what we have observed, which leads us to make the claim that the IGWs in our simulations are generated primarily through plumes.

In the previous section, we attempted to utilise the photometric data from a recent survey by several co-authors of this paper to directly compare the frequency slopes from these observations with the frequency slopes from our simulation. We chose the slopes from the simulated temperature perturbations and found the best slope match for multiple radial IGW rays (which vary with horizontal angles), and frequency intervals. For these 26 stars in particular, which are mainly a mix of β-Cepheid variables, slowly pulsating B-type stars, and stars which just typically show low-frequency variability, the best-match frequency slopes range from −1.87 to −0.74, which we found for different frequency intervals, all starting from 8 μHz.

This work mainly shows the low-frequency variabilities, seen in the photometric data of some stars, could be explained by the IGWs from core convection. Our results indicate that the surface frequency slopes are consistent with the observational results to-date (Bowman et al. 2019, 2020). We found a more recent observational work, Bowman et al. (2019) – where it was found that the frequency slope magnitudes of IGW signatures from a variety of blue supergiants (BSGs) are less than 3.5 – to be the most comparable to the findings of the present work. Unfortunately, we were not able to extract information about the stellar mass from our simulated data. The largest differences between the models that we used lie in the stellar radius, density, Brunt–Väisälä frequency, and thermal diffusivity profiles (which introduces 4 degrees of freedom, calling for a study on how to reduce this). In addition, investigating the effect of varying just one parameter while keeping the others constant would be beneficial, albeit tedious.

In this work, we used a constant solid-body rotation rate for all the simulations. However, O- and B-type stars are known to have angular velocities that vary with radius and polar angles. Thus, the next natural step for this work would be to study the effect of differential rotation for different masses. Studies on the propagation of IGWs through a radial differential rotation have already been done on one stellar mass (Rogers et al. 2013), but repeating this for a broader range of masses is crucial as we do not expect the rotation profile to be the same for all stellar masses. This study on differential rotation will also greatly help us understand the transport of angular momentum and chemical composition, which are two of the more poorly understood fundamental problems in stellar evolution.

Acknowledgments

P.V.F.E. was supported by the US Department of Energy through the Los Alamos National Laboratory (LANL). LANL is operated by Triad National Security, LLC, for the National Nuclear Security Administration of the US Department of Energy (contract no. 89233218CNA000001). This work has been assigned a document release number LA-UR-22-26062.

References

Appendix A: Periodograms of TESS stars versus simulated temperature perturbations

We present the simulated surface temperature perturbations from multiple stellar models, 3 M_⊙ to 13 M_⊙ (see Table 5) in blue, overplotted on the periodograms of selected stars, and in orange, from the TESS catalogue, specifically: from the study in Handler et al. (in prep) and Chowdhury et al. (in prep). The titles of each plot provide the TIC number and mass calculated from HR diagram fitting (see Section 5). The legends show the mass of the stars used in our simulations. The y-axes represent brightness perturbations, where we have converted the temperature perturbations from simulations using the multiplication factor shown in brackets in the legend to align both datasets.

Fig. A.1. Normalised temperature perturbation (solid blue line) from a 3 M⊙ star simulation and brightness variation (orange, solid line) from a 2.51 M⊙ star from the TESS catalogue as a function of wave frequencies.

Fig. A.2. Normalised temperature perturbation (solid blue line) from a 3 M⊙ star simulation and brightness variation (orange, solid line) from a 2.68 M⊙ star from the TESS catalogue as a function of wave frequencies.

Fig. A.3. Normalised temperature perturbation (solid blue line) from a 5 M⊙ star simulation and brightness variation (orange, solid line) from a 4.02 M⊙ star from the TESS catalogue as a function of wave frequencies.

Fig. A.4. Normalised temperature perturbation (solid blue line) from a 5 M⊙ star simulation and brightness variation (orange, solid line) from a 4.58 M⊙ star from the TESS catalogue as a function of wave frequencies.

Fig. A.5. Normalised temperature perturbation (solid blue line) from a 5 M⊙ star simulation and brightness variation (solid orange line) from a 5.50 M⊙ star from the TESS catalogue as a function of wave frequencies.

Fig. A.6. Normalised temperature perturbation (solid blue line) from a 5 M⊙ star simulation and brightness variation (solid orange line) from a 5.99 M⊙ star from the TESS catalogue as a function of wave frequencies.

Fig. A.7. Normalised temperature perturbation (solid blue line) from a 5 M⊙ star simulation and brightness variation (solid orange line) from a 5.99 M⊙ star from the TESS catalogue as a function of wave frequencies.

Fig. A.8. Normalised temperature perturbation (solid blue line) from a 7 M⊙ star simulation and brightness variation (solid orange line) from a 6.18 M⊙ star from the TESS catalogue as a function of wave frequencies.

Fig. A.9. Normalised temperature perturbation (solid blue line) from a 7 M⊙ star simulation and brightness variation (solid orange line) from a 7.55 M⊙ star from the TESS catalogue as a function of wave frequencies.

Fig. A.10. Normalised temperature perturbation (solid blue line) from a 10 M⊙ star simulation and brightness variation (solid orange line) from a 9.55 M⊙ star from the TESS catalogue as a function of wave frequencies.

Fig. A.11. The normalised temperature perturbation (solid blue line) from a 10 M⊙ star simulation and brightness variation (solid orange line) from a 9.82 M⊙ star from the TESS catalogue as a function of wave frequencies.

Fig. A.12. Normalised temperature perturbation (solid blue line) from a 10 M⊙ star simulation and brightness variation (solid orange line) from a 10.10 M⊙ star from the TESS catalogue as a function of wave frequencies.

Fig. A.13. Normalised temperature perturbation (solid blue line) from a 10 M⊙ star simulation and brightness variation (solid orange line) from a 10.47 M⊙ star from the TESS catalogue as a function of wave frequencies.

Fig. A.14. Normalised temperature perturbation (solid blue line) from a 10 M⊙ star simulation and brightness variation (solid orange line) from a 10.50 M⊙ star from the TESS catalogue as a function of wave frequencies.

Fig. A.15. Normalised temperature perturbation (solid blue line) from a 10 M⊙ star simulation and brightness variation (solid orange line) from a 10.79 M⊙ star from the TESS catalogue as a function of wave frequencies.

Fig. A.16. Normalised temperature perturbation (solid blue line) from a 10 M⊙ star simulation and brightness variation (solid orange line) from a 10.96 M⊙ star from the TESS catalogue as a function of wave frequencies.

Fig. A.17. Normalised temperature perturbation (solid blue line) from a 13 M⊙ star simulation and brightness variation (solid orange line) from a 11.26 M⊙ star from the TESS catalogue as a function of wave frequencies.

Fig. A.18. Normalised temperature perturbation (solid blue line) from a 13 M⊙ star simulation and brightness variation (solid orange line) from a 11.31 M⊙ star from the TESS catalogue as a function of wave frequencies.

Fig. A.19. Normalised temperature perturbation (solid blue line) from a 13 M⊙ star simulation and brightness variation (solid orange line) from a 12.60 M⊙ star from the TESS catalogue as a function of wave frequencies.

Fig. A.20. Normalised temperature perturbation (solid blue line) from a 13 M⊙ star simulation and brightness variation (solid orange line) from a 12.96 M⊙ star from the TESS catalogue as a function of wave frequencies.

Fig. A.21. Normalised temperature perturbation (solid blue line) from a 13 M⊙ star simulation and brightness variation (solid orange line) from a 13.20 M⊙ star from the TESS catalogue as a function of wave frequencies.

All Tables

All Figures

	Fig. 1. Density (a) and Brunt–Väisälä frequency (b) profiles as functions of stellar radius, in units of total solar radius. All the models are cut at 90% the total stellar radius to ensure numerical stability.
In the text

	Fig. 2. Prantl number profiles used in the simulations of all the stellar models used in this work.
In the text

Fig. 3. Thermal diffusivities (top panel) as a function of stellar radius and the ratio of linear IGW amplitudes obtained with real stellar thermal diffusivities to those obtained with $\overline{K}$ profiles used in our main simulations (bottom panel) as a function of wave frequencies at the top of the radiation zone. The dotted lines in the top panel indicate the varying thermal diffusivity profiles used for our main simulations, whilst the solid lines indicate the realistic thermal diffusivities from MESA. In the bottom panel, the wavenumbers for all the waves have been set to 1. The initial radius for these IGWs was chosen at random from the middle of the radiation zone and the initial amplitude was set to unity.

In the text

	Fig. 4. Average vrms in the convection zones as a function of time of all the models in Table 2.
In the text

	Fig. 5. Radial velocities as functions of stellar radius and frequency for different stellar masses, all at midMs. The solid, white line shows the Brunt–Väisälä frequency profile for each stellar model, respectively.
In the text

	Fig. 6. Radial velocities (top panels), tangential velocities (middle panels), and temperatures (bottom panels) as a function of frequency for different stellar masses in the middle of the convection zone for a randomly chosen radial ray. The dashed lines show the line fits to each profile at frequencies between 8 μHz and 30 μHz.
In the text

Fig. 7. Ratios of the tangential velocity (vθ) to the radial velocity (vr) as a function of frequency, at different radii throughout the stellar interior for different stellar masses. The black line represents the ratio of the Brunt–Väisälä frequency to the linear IGW frequency, at that particular radius. The first row shows vθ/vr inside the convection zone, the second and third rows show the ratios at 0.1 HP and 0.25 HP, respectively, and the final two rows show the ratios at a randomly chosen radius inside the radiation zone and close to the top of the simulation domain, respectively.

In the text

	Fig. 8. Radial velocities (top panels), tangential velocities (middle panels), and temperatures (bottom panels) as a function of frequency for different stellar masses at 1/4th the pressure scale height from the top of the convection zone for a randomly chosen radial ray. The dashed lines show the line fits to each profile at frequencies between 8 μHz and 30 μHz.
In the text

	Fig. 9. Radial velocities (top panels), tangential velocities (middle panels), and temperatures (bottom panels) as a function of frequency for different stellar masses at 0.82 Rstar for a randomly chosen radial ray. The dashed lines are line-fits to each profile at frequencies between 8 μHz and 30 μHz.
In the text

	Fig. A.1. Normalised temperature perturbation (solid blue line) from a 3 M⊙ star simulation and brightness variation (orange, solid line) from a 2.51 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text

	Fig. A.2. Normalised temperature perturbation (solid blue line) from a 3 M⊙ star simulation and brightness variation (orange, solid line) from a 2.68 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text

	Fig. A.3. Normalised temperature perturbation (solid blue line) from a 5 M⊙ star simulation and brightness variation (orange, solid line) from a 4.02 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text

	Fig. A.4. Normalised temperature perturbation (solid blue line) from a 5 M⊙ star simulation and brightness variation (orange, solid line) from a 4.58 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text

	Fig. A.5. Normalised temperature perturbation (solid blue line) from a 5 M⊙ star simulation and brightness variation (solid orange line) from a 5.50 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text

	Fig. A.6. Normalised temperature perturbation (solid blue line) from a 5 M⊙ star simulation and brightness variation (solid orange line) from a 5.99 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text

	Fig. A.7. Normalised temperature perturbation (solid blue line) from a 5 M⊙ star simulation and brightness variation (solid orange line) from a 5.99 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text

	Fig. A.8. Normalised temperature perturbation (solid blue line) from a 7 M⊙ star simulation and brightness variation (solid orange line) from a 6.18 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text

	Fig. A.9. Normalised temperature perturbation (solid blue line) from a 7 M⊙ star simulation and brightness variation (solid orange line) from a 7.55 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text

	Fig. A.10. Normalised temperature perturbation (solid blue line) from a 10 M⊙ star simulation and brightness variation (solid orange line) from a 9.55 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text

	Fig. A.11. The normalised temperature perturbation (solid blue line) from a 10 M⊙ star simulation and brightness variation (solid orange line) from a 9.82 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text

	Fig. A.12. Normalised temperature perturbation (solid blue line) from a 10 M⊙ star simulation and brightness variation (solid orange line) from a 10.10 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text

	Fig. A.13. Normalised temperature perturbation (solid blue line) from a 10 M⊙ star simulation and brightness variation (solid orange line) from a 10.47 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text

	Fig. A.14. Normalised temperature perturbation (solid blue line) from a 10 M⊙ star simulation and brightness variation (solid orange line) from a 10.50 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text

	Fig. A.15. Normalised temperature perturbation (solid blue line) from a 10 M⊙ star simulation and brightness variation (solid orange line) from a 10.79 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text

	Fig. A.16. Normalised temperature perturbation (solid blue line) from a 10 M⊙ star simulation and brightness variation (solid orange line) from a 10.96 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text

	Fig. A.17. Normalised temperature perturbation (solid blue line) from a 13 M⊙ star simulation and brightness variation (solid orange line) from a 11.26 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text

	Fig. A.18. Normalised temperature perturbation (solid blue line) from a 13 M⊙ star simulation and brightness variation (solid orange line) from a 11.31 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text

	Fig. A.19. Normalised temperature perturbation (solid blue line) from a 13 M⊙ star simulation and brightness variation (solid orange line) from a 12.60 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text

	Fig. A.20. Normalised temperature perturbation (solid blue line) from a 13 M⊙ star simulation and brightness variation (solid orange line) from a 12.96 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text

	Fig. A.21. Normalised temperature perturbation (solid blue line) from a 13 M⊙ star simulation and brightness variation (solid orange line) from a 13.20 M⊙ star from the TESS catalogue as a function of wave frequencies.
In the text