Low-frequency Variability in Massive Stars: Core Generation or Surface Phenomenon?

Massive stars play an important role in astrophysics, driving galactic turbulence that regulates star formation (Hayward & Hopkins 2016) and enriching the chemical composition of the interstellar medium (Nomoto et al. 2013). It is thus important to test and calibrate models of massive star structure and evolution. Asteroseismology is an unparalleled technique to probe the interiors of stars (Aerts et al. 2010), and space-based photometry has revolutionized our understanding of stellar interiors (e.g., Aerts et al. 2019).

The detection of internal gravity waves, which propagate to the bottom of a star's radiative zone, provides precious information about the cores and interiors of massive stars (e.g., Zwintz et al. 2017; Ouazzani et al. 2019). Bowman et al. (2019b, hereafter B19) detected low-frequency variability in a large population of massive stars (following similar detections in fewer stars by, e.g., Blomme et al. 2011). In over 100 massive stars, B19 detected luminosity variations with roughly constant amplitude for frequencies ≲1 days⁻¹, and power laws with slope −2 to −3 at higher frequencies. This "red noise" signal is very different from previous detections of discrete peaks in the power spectrum, which is likely due to linearly unstable g-modes (Pápics et al. 2017). B19 interpret this low-frequency variability as internal gravity waves, which could be excited by core convection.

Much of the early work on wave generation by convection was motivated by the Sun and stars with convective envelopes (e.g., Press 1981; Goldreich & Kumar 1990). These theories were extended to massive stars with convective cores (e.g., Lecoanet & Quataert 2013), and were validated by high-resolution, turbulent, 3D numerical simulations in Cartesian geometry (Couston et al. 2018). They find an excitation spectrum that decreases rapidly for frequencies above the convective turnover frequency. In contrast, work by other authors has found extremely shallow wave spectra at high frequencies (Rogers et al. 2013; Edelmann et al. 2019). There is currently no consensus on the spectrum of internal gravity waves at the radiative–convective boundary.

Simulations of wave excitation by convection do not follow the propagation of internal waves to the surface. To predict surface luminosity fluctuations, Samadi et al. (2010) and Shiode et al. (2013) solved the wave propagation problem using eigenmodes and the Wentzel–Kramers–Brillouin (WKB) approximation, but only considered standing modes. These calculation cannot predict the wave amplitude at frequencies between or below the frequencies of the standing modes. Thus, they do not replicate the "red noise" signal observed in B19. Here we test whether or not the signal can be explained by traveling waves; there is no previous work studying the surface manifestation of traveling waves in massive stars.

The luminosity variations at the surface of a star δL(f; R_⋆) are related to the radial velocity variations at the radiative–convective boundary $\delta {u}_{r}(f;{r}_{\mathrm{RCB}})$ by

$$\begin{eqnarray}&&\delta L(f;{R}_{\star })=T(f)\delta {u}_{r}(f;{r}_{\mathrm{RCB}})\end{eqnarray} \tag{ 1 }$$

for linear waves. Here f is the frequency of the wave, and T(f) is the wave transfer function. We check the linearity assumption in Section 3. Equation (1) separates out the nonlinear physics related to convection, which determines $\delta {u}_{r}(f;{r}_{\mathrm{RCB}})$ and is uncertain, from the linear physics of wave propagation, which gives T(f).

We will discuss the wave transfer function for a 10 M_⊙, solar-metallicity star near the zero-age main sequence (ZAMS). The stellar structure model is calculated using MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019). Appendix G shows the wave transfer function for stars from 3 M_⊙ to 20 M_⊙, and from the ZAMS to the terminal-age main sequence (TAMS). The broad features of the transfer functions are similar in this range of masses and ages.

We calculate the wave transfer function in two ways. First, we use non-adiabatic eigenmodes calculated with GYRE (Townsend et al. 2018). We assume the non-adiabatic eigenmodes form a complete basis, and expand the perturbations in this basis. We force the horizontal displacement with $\exp (-i2\pi {ft})$ at a forcing radius r_f above the radiative–convective boundary. This produces a luminosity variation at the surface L(R_⋆) with frequency f. We then calculate L(R_⋆) for multiple forcing radii in a 0.02R_⋆ range above the radiative–convective boundary. We force at multiple radii to avoid the nodes of the eigenfunctions. T(f) is the average of the amplitude of L(R_⋆) over these forcing radii.

To validate this calculation, we also time-evolve the linearized non-adiabatic oscillation equations in the Cowling approximation (Equations (33)–(39)) using Dedalus (Burns et al. 2019). We solve the equations from the radiative–convective boundary (r_RCB ≈ 0.22R_⋆) to the bottom of the surface convection zone (r_top ≈ 0.98R_⋆). We cannot include the core or surface convection zones in the simulation; otherwise, linearly unstable convective modes dominate. We force the waves with a boundary condition $\delta {u}_{r}({r}_{\mathrm{RCB}})=\sin (2\pi {ft})$. After an initial transient, the luminosity at the top of the domain, $\delta L({r}_{\mathrm{top}})$, also varies sinusoidally, which allows us to measure T(f). We find good agreement between these two methods of calculating the wave transfer function. Appendices A–F include more details of the calculations, as well as the comparison between them.

Figure 1 shows the wave transfer function calculated using GYRE for spherical harmonic degrees ℓ = 1, 2, and 3. The ℓ = 1 mode dominates the surface luminosity fluctuations (Figure 2). For each ℓ, we see the transfer function is smooth at low frequencies, as is observed. However, the transfer function is also very small due to strong wave damping at low frequencies. If the observed low-frequency variability is due to waves from the core, there should be very little power for f ≲ 0.1 days⁻¹, because these waves are strongly damped. The observed luminosity fluctuations are roughly constant even for f ≲ 0.1 days⁻¹ (Figure 2 and Blomme et al. 2011). However, it is likely the low-frequency observations are dominated by systematic errors, so this inconsistency at low frequencies does not rule out the detection of propagating internal gravity waves at higher frequencies.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 1 also shows large-amplitude, regularly spaced peaks at frequencies ≳1 day⁻¹. The peaks are at the frequencies of standing g-modes. These frequencies are the resonant frequencies for propagating internal gravity waves. Peaks with this type of regular spacing are not present in the observed data. Thus, we believe the observed low-frequency variability in high-mass stars is not due to linearly propagating internal gravity waves generated in the core.

Because the observed power spectra do not have regularly spaced peaks, one may be tempted to attribute the low-frequency variability to propagating, rather than standing, internal gravity waves. However, propagating internal gravity waves only produce a smooth spectrum when the damping rate is similar or greater than the mode frequency spacing. This is only true for waves that are strongly damped, which are too low amplitude to be observed. This dichotomy is illustrated by the wave transfer functions in Figure 1: the transfer function is smooth for low frequencies, but is also very small; whereas at higher frequencies it is dominated by large-amplitude peaks.

To illustrate this point, in Figure 2 we compare the observed spectrum of EPIC 234517653 from B19 to two theoretical spectra. We choose this star because it has spectral type B2, roughly corresponding to our 10 M_⊙ stellar model, and has a relatively simple spectrum without strong rotation effects or signature of linearly unstable modes. Although we plot the "residual" spectrum, the "original" spectrum is qualitatively similar.

The theoretical spectrum depends on the wave excitation spectrum $\delta {u}_{r}(f;{r}_{\mathrm{RCB}})$ via Equation (1). To show a range of possibilities, we consider two very different excitation spectra reported in numerical simulations. Rogers et al. (2013) ran a series of 2D cylindrical simulations using the anelastic equations, and found $\delta {u}_{r}\sim {f}^{0.4}{{\rm{\Lambda }}}^{-0.9}$, where ${\rm{\Lambda }}=\sqrt{{\ell }({\ell }+1)}$ (Ratnasingam et al. 2019); we refer to results using this spectrum with R. Couston et al. (2018) ran a series of 3D cartesian simulations using the Boussinesq equations, and found $\delta {u}_{r}\sim {f}^{-3.25}{{\rm{\Lambda }}}^{5/2};$ we refer to results using this spectrum with LQ, as they match with the theoretical prediction in Lecoanet & Quataert (2013). We then sum luminosity variations over ℓ using the Dziembowski (1977) relations as in Shiode et al. (2013). To maintain the same frequency resolution as the observations, we calculate the maximum brightness fluctuation in frequency bins with the same width as the frequency spacing of the observational data from B19.

Figure 2 shows the luminosity variation spectrum for EPIC 234517653, as well as the R and LQ spectra. The overall normalization of the theoretical spectra are scaled to more easily compare to the observations. The predicted normalization of the LQ spectrum is similar what we plot, whereas the predicted normalization of the R spectrum is about 4 orders of magnitude higher than what we plot (Appendix H). The theoretical spectra inherit two important properties from the wave transfer function (Figure 1): (1) at low frequencies (f ≲ 0.1 days⁻¹), the luminosity perturbations become very small due to strong wave damping; and (2) at frequencies ≳0.5 day⁻¹, there are regularly spaced peaks due to the near-resonant effects of standing modes. Neither of these features are in the observations, which show a fairly flat spectrum at low frequencies, without regularly spaced peaks at higher frequencies. Our calculations do not take into account rotation. Rotation would increase the number of resonant peaks at high frequencies, but also decrease the luminosity variation at low frequencies. Appendix I shows that the high-frequency waves would still manifest as a series of resonant peaks, so rotation does not change the main features of Figure 2.

Although here we only compare to a single star, none of the spectra in B19 show a large number of regularly spaced peaks similar to the theoretical spectra in Figure 2. Even stellar power spectra that do have peaks (e.g., EPIC 202061164) only have a couple, and thus do not match the theoretical spectra for convectively excited internal gravity waves. The standing waves in the R spectrum reach very large amplitudes, which would probably lead to nonlinear interactions. However, in the next section we show that the observed variability is too small to be nonlinear.

In the previous section, we showed that linearly propagating waves have a surface luminosity spectrum inconsistent with observations. Here we will investigate possible nonlinear effects. We argue that traveling waves do not become strongly nonlinear, i.e., break. While standing waves could in principle experience weak nonlinearities, we show that the observed luminosity fluctuations are so low that an equivalent standing wave would be stable to nonlinear instabilities. Thus, all the peaks of the predicted spectrum should be observable.

3.1. Traveling Wave Breaking

In the absence of damping, waves conserve the wave luminosity as they propagate outward,

$$\begin{eqnarray}&&{L}_{\mathrm{wave}}=\displaystyle \frac{4\pi {r}^{2}\rho {\omega }^{4}}{{{Nk}}_{h}}| {\xi }_{h}{| }^{2}.\end{eqnarray} \tag{ 2 }$$

Here ξ_h is the horizontal wave displacement, ω is the wave's angular frequency, N is the Brunt–Väisälä frequency, and ${k}_{h}\equiv {\rm{\Lambda }}/r$ and k_r are the horizontal and radial components of the wavenumber. In deriving Equation (2) we have approximated the dispersion relation as $\omega ={{Nk}}_{h}/{k}_{r}$.

As waves propagate outward, ρ decreases, leading ξ_h to increase. When ξ_h k_h ∼ 1, waves break quickly in about a wave period (e.g., Staquet & Sommeria 2002; Liu et al. 2010; Eberly & Sutherland 2014). Thus, the criterion for wave breaking is

$$\begin{eqnarray}&&{L}_{\mathrm{wave}}\gtrsim \displaystyle \frac{4\pi {r}^{2}\rho {\omega }^{4}}{{{Nk}}_{h}^{3}}.\end{eqnarray} \tag{ 3 }$$

The total wave luminosity is expected to be smaller than the convective luminosity by a factor of the convective Mach number (e.g., Goldreich & Kumar 1990; Lecoanet & Quataert 2013); this scaling was verified in the simulations of Rogers et al. (2013) and Couston et al. (2018).

Radiative damping is significant for some waves. Waves are weakly affected by damping if they propagate further than one pressure scaleheight H over a damping time,

$$\begin{eqnarray}&&\displaystyle \frac{\omega }{{k}_{r}H}\gtrsim {K}_{\mathrm{rad}}{k}_{r}^{2},\end{eqnarray} \tag{ 4 }$$

where K_rad is the radiative diffusivity. This is a conservative estimate: waves need to propagate many scaleheights for them to reach the surface without damping. We can rewrite this in terms of the radiative luminosity using ${L}_{\mathrm{rad}}\sim (4\pi {r}^{2}){K}_{\mathrm{rad}}\rho {N}^{2}H$. Then a requirement for waves to not experience significant damping is

$$\begin{eqnarray}&&\displaystyle \frac{4\pi {r}^{2}\rho {\omega }^{4}}{{{Nk}}_{h}^{3}}\gtrsim {L}_{\mathrm{rad}}.\end{eqnarray} \tag{ 5 }$$

If waves are radiatively damped, they will not reach sufficient amplitudes to break. So combining this with Equation (3), we find a requirement for waves to break is

$$\begin{eqnarray}&&{L}_{\mathrm{wave}}\gtrsim {L}_{\mathrm{rad}}.\end{eqnarray} \tag{ 6 }$$

Since the wave luminosity is significantly smaller than the convective luminosity, this relation is not satisfied in almost all stars (some stars in the last year of their life may have ${L}_{\mathrm{wave}}\gt {L}_{\mathrm{rad}}$; see Quataert & Shiode 2012). Thus, traveling waves will not experience wave breaking. Note, however, that many numerical simulations artificially increase the convective luminosity by many orders of magnitude (e.g., Rogers et al. 2013; Jones et al. 2017; Edelmann et al. 2019), which likely explains reports of breaking waves generated by convection in Rogers et al. (2013) and Edelmann et al. (2019). Artificially increasing the convective (and hence wave) luminosity makes it difficult to study wave excitation and propagation, because the wave physics can be very different.

3.2. Nonlinear Stability of Standing Waves

Figure 1 shows that waves excited near the star's g-mode frequencies can experience resonant amplification. Because the resonant waves are coherent for many wave periods, weakly nonlinear effects can build up over time, and lead to significant transfer of energy between frequencies without necessitating wave breaking.

However, the observed luminosity fluctuations are very small. Even if the entire luminosity variation was due to standing modes, the modes would be stable to nonlinear instability for most low-ℓ standing modes in 10 M_⊙ ZAMS stars. As discussed in Weinberg et al. (2012), nonlinear instabilities of a "parent" mode with frequency f and amplitude ξ_h k_h requires the existence of two other modes satisfying

$$\begin{eqnarray}&&{\xi }_{h}{k}_{h}f\kappa \gt \sqrt{{\gamma }_{1}{\gamma }_{2}}{\left[1+\displaystyle \frac{\delta {f}^{2}}{{\left({\gamma }_{1}+{\gamma }_{2}\right)}^{2}}\right]}^{1/2},\end{eqnarray} \tag{ 7 }$$

where κ is a spatial overlap integral, γ₁ and γ₂ are the damping rates of the two "daughter" modes, and δf is the frequency difference between the parent and daughter modes, all measured in units of day⁻¹. The difference between the ℓ of the daughter modes can be no greater than the ℓ of the parent mode.

We calculate the threshold for nonlinear instability by finding the smallest possible ξ_h k_h satisfying Equation (7) for each mode. We search over modes with ℓ ≤ 20 with daughter mode frequencies less than the parent mode frequency so they can efficiently damp. We assume that the spatial overlap integral κ is small if the difference between the radial mode orders is greater than two; otherwise we take κ = 1. In Figure 3, we plot the surface luminosity fluctuations associated with these marginally stable modes. All but one ℓ = 1 mode would be nonlinearly stable, even if we assume the observed fluctuations are entirely due to those modes. We also find nonlinear stability for modes with f ≤ 1 day⁻¹.

Zoom In Zoom Out Reset image size

Less-massive and more evolved stars have more weakly damped modes, so some standing modes experience nonlinear instabilities in these stars. We also note that some stars in B19 show several peaks, likely due to linearly unstable g-modes. In many cases, those peaks are at amplitudes that are similar to the marginally stable amplitude calculated in Figure 3. Many other stars show signatures of nonlinear wave interaction (e.g., Degroote et al. 2009; Bowman et al. 2016). In pulsating white dwarfs (ZZ Ceti stars), the mode amplitudes are also believed to be limited by weak nonlinearities (e.g., Wu & Goldreich 2001; Luan & Goldreich 2018). However, in none of these cases is there any indication that these interactions can produce the "red noise" spectra observed in B19.

We have argued that the observed low-frequency variability is massive stars is unlikely due to internal gravity waves generated by core convection. Another possible source of the variability is surface phenomena.

Stars with spectral types later than F have well-developed surface convection zones driven by hydrogen ionization. In these stars, turbulent convection results in a specific photometric signature: granulation. In the power spectrum, granulation shows up as an excess of power at low frequencies (red noise). 3D radiative hydrodynamical simulations of surface convection can match observations of the disk-integrated intensity of main-sequence stars and red giants (e.g., Samadi et al. 2013a, 2013b, and references therein). Convective surface velocity fluctuations also affect stellar spectra. If the surface velocity fluctuations are correlated over the entire line-forming region, this is referred to as macroturbulence. Otherwise it is referred to as microturbulence.

Both micro- and macro- turbulence have also been observed in OB stars, with amplitudes ranging from a few to tens of km s⁻¹. These stars are not convective at the surface, but possess subsurface convective layers due to partial ionization of helium and iron group elements (Cantiello et al. 2009; Cantiello & Braithwaite 2019). Cantiello et al. (2009) and Grassitelli et al. (2015a, 2015b) show that the presence and vigor of subsurface convection is correlated to spectroscopically derived surface velocity fluctuations in OB stars. This is true for all stars except those with magnetic fields strong enough to substantially affect the iron subsurface convection zone (Sundqvist et al. 2013). This suggests that subsurface convection could play an important role in inducing stellar surface variability, even for OB stars.¹²

It is thus natural to ask if the observed low-frequency variability in massive stars is also caused by subsurface convection. This was discussed in Bowman et al. (2019a), which showed that the variability in massive stars has a characteristic frequency ν_char of a few tens of μHz, with little dependency on the star's position on the HR-diagram. This is inconsistent with the "KB-scaling" of convective frequency with stellar properties derived in Kjeldsen & Bedding (1995) (see Figure 8 in B19).

However, Kjeldsen & Bedding (1995) derived this scaling for surface convection. The KB-scaling matches observations of stars with surface temperatures below approximately 10–11,000 K, which have surface convection driven by hydrogen and helium ionization. For hotter stars, convection occurs below the surface around temperatures corresponding to partial ionization of He ii (≈45 kK) and Fe (≈150 kK). These subsurface convective regions (He iiCZ and FeCZ) occur at the specified temperatures; if a star expands, they will move to deeper mass coordinate. As such, they do not directly trace stellar surface properties, and one does not expect their induced surface photometric fluctuations to follow the KB-scaling. Instead, models predict a fairly constant value of their characteristic frequency as function of the stellar parameters, with typical frequencies of 6–60 μHz (Cantiello et al. 2009; Cantiello & Braithwaite 2011, 2019), which are consistent with the range of characteristic frequencies observed in the Bowman et al. (2019a) sample. Hence, we argue the distribution of characteristic frequencies of the observed low-frequency variability may support a subsurface origin.

The low-frequency variability is also observed in low-metallicity Large Magellanic Cloud (LMC) stars (B19). Since the occurrence of the FeCZ is metallicity dependent, this could rule out subsurface convection as the cause of the variability. In 1D stellar evolution calculations, the FeCZ occurs above a luminosity of $L\approx {10}^{3.2}\,{L}_{\odot }$ for solar-metallicity stars, corresponding to a ZAMS star of about $7\,{M}_{\odot }$. At the lower metallicity of the LMC, the FeCZ appears at higher luminosities (L ≈ 10^3.9 ${L}_{\odot }$), corresponding to a ZAMS star of about $11\,{M}_{\odot }$ (Cantiello et al. 2009). However, it is unclear if the objects in the TESS LMC sample of B19 are hot, main-sequence stars probing this FeCZ transition luminosity. For example, if the sample contains mostly luminous giants and supergiants, then the stars could be either above the FeCZ luminosity threshold for LMC metallicity, or posses surface convection due to H and He ionization. This would also be consistent with the much larger amplitude of the photometric variability of the LMC sample, as expected in more massive, brighter stars, independent of the scenario considered (subsurface convection or core internal gravity waves).

Finally, some of the stars showing low-frequency, stochastic variability are low-mass, main-sequence A stars. These are not expected to show the presence of an FeCZ. However, they are affected by other (sub)surface convection zones due to the ionization of H and He (Cantiello & Braithwaite 2019). While the amplitude of the velocity fluctuations induced by these convective regions tend to be smaller than for the FeCZ, characteristic frequencies are also in the range of tens of μHz.

Therefore, subsurface convection remains a viable explanation for the observed variability. More work is required to test this hypothesis and firmly establish—or rule out—a connection between subsurface turbulent convective motions and the wide-spread, low-frequency photometric variability observed in massive stars.

• 1.  
  We argue that the low-frequency variability presented in B19 is not due to linearly propagating waves from the core.

  • i.  
    At low frequencies, linear waves have a smooth spectrum, but very little surface luminosity variation due to strong radiative damping. This is not observed.
  • ii.  
    Linearly propagating waves at frequencies f ≳ 0.5 day⁻¹ form many standing modes (Figure 2), which are also not observed.
• 2.  
  Wave propagation can be approximated as linear.

  • i.  
    Low-frequency traveling waves do not break.
  • ii.  
    If high-frequency variability was due to standing modes, many would be at such low amplitude they would be stable to nonlinear interactions.
• 3.  
  We conclude that the spectra presented in B19 are not due to IGWs generated by core convection.
• 4.  
  However, almost all massive stars have subsurface convection zones that could produce the observed low-frequency variability.

We are grateful for extremely useful discussions with Conny Aerts, Dominic Bowman, and Siemen Burssens regarding the observations presented in B19. We would like to thank Rich Townsend for valuable discussions about GYRE and stellar oscillation modes. We would also like to thank Tamara Rogers, Rathish Ratnasingam, and Philipp Edelmann for useful conversations about wave nonlinearity. We would also like to thank the anonymous referee who pointed out some additional reasons why the observations of B19 are unlikely to be due to waves generated in the core. D.L. is supported by a PCTS and Lyman Spitzer Jr. fellowship. D.L. would like to thank his fellow members of the 10th New York County Grand Jury #6 of 2019 for their amiable company while he composed this paper during his four weeks of jury service. The Center for Computational Astrophysics at the Flatiron Institute is supported by the Simons Foundation. This work has benefited from many conversations at meetings supported by Grant GBMF5076 from the Gordon and Betty Moore Foundation. E.Q. thanks the Princeton Astrophysical Sciences department and the theoretical astrophysics group and Moore Distinguished Scholar program at Caltech for their hospitality and support. L.A.C. acknowledges funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 793450. L.A.C., B.F., and M.L.B. acknowledge funding by the European Research Council under the European Union's Horizon 2020 research and innovation program through grant No. 681835-FLUDYCO-ERC-2015-CoG. This work was performed in part under contract with the Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute. B.J.S.P. acknowledges being on the traditional territory of the Lenape Nations and recognizes that Manhattan continues to be the home to many Algonkian peoples. We give blessings and thanks to the Lenape people and Lenape Nations in recognition that we are carrying out this work on their indigenous homelands.

We use the non-adiabatic g-mode oscillations of a star as a basis to expand velocity perturbations. For a given spherical harmonic degree ℓ, we have

$$\begin{eqnarray}&&{{\boldsymbol{u}}}_{{\ell }}(r,t)=\displaystyle \sum _{\omega ^{\prime} }{A}_{{\ell }}(t;\omega ^{\prime} ){\tilde{{\boldsymbol{u}}}}_{{\ell }}(r;\omega ^{\prime} )\exp (-i\omega ^{\prime} t),\end{eqnarray} \tag{ 8 }$$

where ${\tilde{{\boldsymbol{u}}}}_{{\ell }}$ corresponds to the velocity eigenfunction, and $\omega ^{\prime}$ its (complex) frequency. For simplicity of notation, here we will write frequencies in terms of the angular frequency $\omega =2\pi f$, and write velocity perturbations as ${\boldsymbol{u}}$, instead of $\delta {\boldsymbol{u}}$ as in the main text. The velocity eigenfunctions are related to the displacement eigenfunctions via ${\tilde{{\boldsymbol{u}}}}_{{\ell }}=-i\omega ^{\prime} {\tilde{{\boldsymbol{\xi }}}}_{{\ell }}$. All vectors are assumed to be composed of two components: a radial component, e.g., u_r; and a horizontal component which is assumed to be proportional to the horizontal gradient of a spherical harmonic function, e.g., ${{\boldsymbol{u}}}_{h}={u}_{h}{{\boldsymbol{\nabla }}}_{h}{Y}_{{\ell },m}(\theta ,\phi )$, where ${{\boldsymbol{\nabla }}}_{h}$ represents a horizontal angular derivative.

The eigenfunctions satisfy

$$\begin{eqnarray}&&-i\omega ^{\prime} {\tilde{{\boldsymbol{u}}}}_{{\ell }}={ \mathcal L }({\tilde{{\boldsymbol{u}}}}_{{\ell }}),\end{eqnarray} \tag{ 9 }$$

where ${ \mathcal L }$ is the linear wave operator. We define the dual basis ${\tilde{{\boldsymbol{u}}}}_{{\ell }}^{\dagger }$ by

$$\begin{eqnarray}&&\langle {\tilde{{\boldsymbol{u}}}}_{{\ell }}^{\dagger }(\omega ^{\prime} ),{\tilde{{\boldsymbol{u}}}}_{{\ell }}(\omega ^{\prime\prime} )\rangle ={\delta }_{\omega ^{\prime} ,\omega ^{\prime\prime} },\end{eqnarray} \tag{ 10 }$$

where δ is the Kronecker delta function, and the inner product is defined as

$$\begin{eqnarray}&&\langle {\boldsymbol{f}},{\boldsymbol{g}}\rangle \equiv {\int }_{0}^{{R}_{\star }}4\pi \rho {r}^{2}{{\boldsymbol{f}}}^{* }\,{\boldsymbol{\cdot }}\,{\boldsymbol{g}}{dr}.\end{eqnarray} \tag{ 11 }$$

We can define the dual basis in this way because the eigenvalues are non-degenerate. The set of all adiabatic modes are a complete basis (e.g., Eisenfeld 1969). Non-adiabatic modes are likely not complete (R. H. D. Townsend 2019, personal communication). Nevertheless, we expect they can represent almost all functions in the regions of the star that do not experience strong damping. The good agreement with calculations using Dedalus indicates the non-adiabatic modes can faithfully represent traveling waves.

We now introduce a forcing to the system,

$$\begin{eqnarray}&&{\partial }_{t}{{\boldsymbol{u}}}_{{\ell }}(r,t)={ \mathcal L }({{\boldsymbol{u}}}_{{\ell }})+{{\boldsymbol{f}}}_{{\ell }}(r,t;\omega ).\end{eqnarray} \tag{ 12 }$$

The forcing has magnitude ${F}_{\omega ,{\ell }}$, is assumed to act in only the horizontal directions, at the forcing radius r_f, and at frequency ω,

$$\begin{eqnarray}&&{{\boldsymbol{f}}}_{{\ell }}(t;\omega )={F}_{\omega ,{\ell }}\,\delta (r-{r}_{f})\exp (-i\omega t){{\boldsymbol{\nabla }}}_{h}{Y}_{{\ell },m}(\theta ,\phi ),\end{eqnarray} \tag{ 13 }$$

where δ is the Dirac delta function. We can solve Equation (12) by using the expansion in normal modes (Equation (8)) and projecting out with ${\tilde{{\boldsymbol{u}}}}_{{\ell }}^{\dagger }(\omega )$. This gives

$$\begin{eqnarray}&&\begin{array}{l}\left[-i\omega ^{\prime} {A}_{{\ell }}(t;\omega ^{\prime} )+{\partial }_{t}{A}_{{\ell }}(t;\omega ^{\prime} )\right]\exp (-i\omega ^{\prime} t)\\ \quad =\,-i\omega ^{\prime} {A}_{{\ell }}(t;\omega ^{\prime} )\exp (-i\omega ^{\prime} t)+\langle {\tilde{{\boldsymbol{u}}}}_{{\ell }}^{\dagger }(\omega ^{\prime} ),{\boldsymbol{f}}(\omega )\rangle .\end{array}\end{eqnarray} \tag{ 14 }$$

Defining ${\tilde{{\boldsymbol{u}}}}_{h}^{\dagger }={\tilde{u}}_{h}^{\dagger }{{\boldsymbol{\nabla }}}_{h}{Y}_{{\ell },m}(\theta ,\phi )$ and assuming ρ is a function of radius, we can simplify this to

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}{A}_{{\ell }}(t;\omega ^{\prime} ) & = & 4\pi \rho ({r}_{f}){r}_{f}^{2}{\tilde{u}}_{h}^{\dagger }({r}_{f};\omega ^{\prime} ){\ell }({\ell }+1){F}_{\omega ,{\ell }}\\ & & \times \,\exp \left[-i\left(\omega -\omega ^{\prime} \right)t\right].\end{array}\end{eqnarray} \tag{ 15 }$$

This can be solved together with the boundary condition A_ℓ = 0 at t = 0,

$$\begin{eqnarray}\begin{array}{rcl}{A}_{{\ell }}(t;\omega ^{\prime} ) & = & i4\pi \rho ({r}_{f}){r}_{f}^{2}{\tilde{u}}_{h}^{\dagger }({r}_{f};\omega ^{\prime} ){\ell }({\ell }+1){F}_{\omega ,{\ell }}\\ & & \times \,\displaystyle \frac{\exp \left[-i\left(\omega -\omega ^{\prime} \right)t\right]-1}{\omega -\omega ^{\prime} }.\end{array}\end{eqnarray} \tag{ 16 }$$

Substituting this into Equation (8), we find

$$\begin{eqnarray}&&\begin{array}{l}{{\boldsymbol{u}}}_{{\ell }}(r,t)=4\pi \rho ({r}_{f}){r}_{f}^{2}{\ell }({\ell }+1){F}_{\omega ,{\ell }}\\ \ \times \,\left[\displaystyle \sum _{\omega ^{\prime} }i{\tilde{u}}_{h}^{\dagger }\left({r}_{f};\omega ^{\prime} \right){\tilde{{\boldsymbol{u}}}}_{{\ell }}\left(r;\omega ^{\prime} \right)\displaystyle \frac{\exp \left(-i\omega t\right)-\exp \left(-i\omega ^{\prime} t\right)}{\omega -\omega ^{\prime} }\right].\end{array}\end{eqnarray} \tag{ 17 }$$

If all the eigenmodes are damped, the imaginary part of $\omega ^{\prime}$ is negative. As $t\to \infty$, the second exponential term drops out, and we find

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{u}}}_{{\ell }}(r,t) & = & 4\pi \rho ({r}_{f}){r}_{f}^{2}{\ell }({\ell }+1){F}_{\omega ,{\ell }}\\ & & \times \,\left[\displaystyle \sum _{\omega ^{\prime} }\displaystyle \frac{i{\tilde{u}}_{h}^{\dagger }({r}_{f};\omega ^{\prime} ){\tilde{{\boldsymbol{u}}}}_{{\ell }}(r;\omega ^{\prime} )}{\omega -\omega ^{\prime} }\right]\exp (-i\omega t).\end{array}\end{eqnarray} \tag{ 18 }$$

The surface luminosity fluctuations are then

$$\begin{eqnarray}\begin{array}{rcl}\delta {L}_{{\ell }}({R}_{\star },t) & = & 4\pi \rho ({r}_{f}){r}_{f}^{2}{\ell }({\ell }+1){F}_{\omega ,{\ell }}\\ & & \times \,\left[\displaystyle \sum _{\omega ^{\prime} }\displaystyle \frac{i{\tilde{u}}_{h}^{\dagger }({r}_{f};\omega ^{\prime} )\delta {\tilde{L}}_{{\ell }}({R}_{\star };\omega ^{\prime} )}{\omega -\omega ^{\prime} }\right]\exp (-i\omega t),\end{array}\end{eqnarray} \tag{ 19 }$$

where $\delta {\tilde{L}}_{{\ell }}({R}_{\star };\omega ^{\prime} )$ is the surface luminosity perturbation of the mode with frequency $\omega ^{\prime}$. In practice, for many stellar models, we find that all modes calculated using GYRE are damped. Some stellar models do have some unstable modes with small positive growth rates. Nevertheless, we treat the time dependence of these modes in the same way as damped modes.

Theories (e.g., Lecoanet & Quataert 2013) and simulations (e.g., Rogers et al. 2013) make predictions for how the radial velocity perturbations at the radiative–convective boundary depends on frequency and spherical harmonic degree, ${u}_{r}({r}_{\mathrm{RCB}};\omega ,{\ell })$. To use Equation (19), we need to establish a link between the radial velocity perturbations and the forcing amplitude ${F}_{\omega ,{\ell }}$.

To illustrate this connection, we will solve the bulk and interfacing forcing problems for an inviscid, adiabatic, Boussinesq fluid in Cartesian geometry. We assume the background has a constant Brunt–Väisälä frequency, N, ranges from z = 0 to L, that gravity is in the z direction, and that the perturbations are periodic in the horizontal directions. For impenetrable boundary conditions (u_z(0) = u_z(L) = 0), the eigenmodes are

$$\begin{eqnarray}&&{u}_{z}(z;{k}_{h})=\sin ({k}_{z}z),\end{eqnarray} \tag{ 20 }$$

where

$$\begin{eqnarray}&&{k}_{z}^{2}(\omega )={k}_{h}^{2}\left(\displaystyle \frac{{N}^{2}}{{\omega }^{2}}-1\right).\end{eqnarray} \tag{ 21 }$$

The horizontal wavenumber k_h is defined by $-{k}_{h}^{2}{u}_{z}={{\rm{\nabla }}}_{h}^{2}{u}_{z}$. The eigenvalues are frequencies ω which satisfy k_z L = nπ, for an integer n.

Now consider the problem with boundary forcing. Our new boundary conditions are ${u}_{z}(0)={F}_{B}\sin (\omega t)$, and u_z(L) = 0. The solution to this new problem is

$$\begin{eqnarray}&&{u}_{z}={F}_{B}\left[\cos ({k}_{z}z)-\displaystyle \frac{\cos ({k}_{z}L)}{\sin ({k}_{z}L)}\sin ({k}_{z}z)\right]\sin (\omega t).\end{eqnarray} \tag{ 22 }$$

Recall that k_z is a function of the forcing frequency (Equation (21)). One can verify that this satisfies the boundary conditions. We can calculate the amplitude of u_z using

$$\begin{eqnarray}&&{\left(\displaystyle \frac{1}{L}{\int }_{0}^{L}{u}_{z}^{2}{dz}\right)}^{1/2}=\sqrt{\displaystyle \frac{1}{2}}\displaystyle \frac{{F}_{B}}{| \sin ({k}_{z}L)| }+{ \mathcal O }\left({\left({k}_{z}L\right)}^{-1/2}\right).\end{eqnarray} \tag{ 23 }$$

The amplitude diverges when k_z L = n π, for n an integer. This corresponds to the frequencies of the eigenmodes of the system with impenetrable boundary conditions. Thus, the solution to the forced problem is related to the unforced eigenvalue problem, even though they have different boundary conditions. Because the problem has no dissipation, the amplitude divergences when forced at an eigenfrequency; if thermal dissipation is included, the divergence is regularized as the eigenfrequencies become complex.

Next consider the bulk forcing problem. Here we solve the problem with impenetrable boundary conditions, but including a forcing term to the horizontal velocity equation, similar to Equation (12),

$$\begin{eqnarray}&&{\partial }_{t}{{\boldsymbol{u}}}_{h}+{{\boldsymbol{\nabla }}}_{h}p={F}_{\omega ,{k}_{h}}\delta (z-{z}_{f})\cos (\omega t){{\boldsymbol{\nabla }}}_{h}g({{\boldsymbol{x}}}_{h}),\end{eqnarray} \tag{ 24 }$$

where $g({{\boldsymbol{x}}}_{h})$ satisfies ${{\rm{\nabla }}}_{h}^{2}g=-{k}_{h}^{2}g$. Away from z_f, u_z is sinusoidal with wavenumber k_z. There are two additional conditions on u_z at z = z_f. The first is that ${\partial }_{z}{u}_{z}$ is continuous, and the second is

$$\begin{eqnarray}&&{\partial }_{t}{\rm{\Delta }}{u}_{z}={k}_{h}^{2}{F}_{\omega ,{k}_{h}}\cos (\omega t)g({{\boldsymbol{x}}}_{h}),\end{eqnarray} \tag{ 25 }$$

where ${\rm{\Delta }}{u}_{z}={u}_{z}({z}_{f}+\epsilon )-{u}_{z}({z}_{f}-\epsilon ),$ in the limit that $\epsilon \to 0$. The solution is

$$\begin{eqnarray}{u}_{z}=\left\{\begin{array}{ll}-\displaystyle \frac{{k}_{h}^{2}}{\omega }{F}_{\omega ,{k}_{h}}\left[\cos ({k}_{z}{z}_{f})\tfrac{\cos ({k}_{z}L)}{\sin ({k}_{z}L)}+\sin ({k}_{z}{z}_{f})\right]\sin ({k}_{z}z)\sin (\omega t) & z\lt {z}_{f}\\ -\displaystyle \frac{{k}_{h}^{2}}{\omega }{F}_{\omega ,{k}_{h}}\tfrac{\cos ({k}_{z}{z}_{f})}{\sin ({k}_{z}L)}\sin ({k}_{z}(z-L))\sin (\omega t) & z\gt {z}_{f}\end{array}\right.\end{eqnarray} \tag{ 26 }$$

where for simplicity we have dropped the horizontal dependence. Averaging over the region above z_f, we find

$$\begin{eqnarray}&&\begin{array}{l}{\left(\displaystyle \frac{1}{L-{z}_{f}}{\int }_{{z}_{f}}^{L}{u}_{z}^{2}{dz}\right)}^{1/2}\equiv {\overline{u}}_{z}({z}_{f})\\ \quad =\,\sqrt{\displaystyle \frac{1}{2}}\displaystyle \frac{{k}_{h}^{2}}{\omega }{F}_{\omega ,{k}_{h}}\displaystyle \frac{| \cos ({k}_{z}{z}_{f})| }{| \sin ({k}_{z}L)| }+{ \mathcal O }\left({\left({k}_{z}L\right)}^{-1/2}\right).\end{array}\end{eqnarray} \tag{ 27 }$$

As in the boundary-forced problem, we find a divergence in the response when ${k}_{z}L=n\pi$, for n an integer (assuming $\cos ({k}_{z}{z}_{f})$ is not also zero). If we also average ${\overline{u}}_{z}$ over a range ${\rm{\Delta }}z$ of z_f, we find

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{1}{{\rm{\Delta }}z}{\displaystyle \int }_{{z}_{0}}^{{z}_{0}+{\rm{\Delta }}z}{\overline{u}}_{z}({z}_{f}){{dz}}_{f} & = & \displaystyle \frac{\sqrt{2}}{\pi }\displaystyle \frac{{k}_{h}^{2}}{\omega }{F}_{\omega ,{k}_{h}}\displaystyle \frac{1}{| \sin ({k}_{z}L)| }\\ & & +\,{ \mathcal O }\left({\left({k}_{z}{\rm{\Delta }}z\right)}^{-1/2}\right),\end{array}\end{eqnarray} \tag{ 28 }$$

provided k_zΔz ≫ 1.

Thus, if we average over a sufficiently large range of forcing locations, we find the boundary-forced and bulk-forced responses are equivalent, with the identification

$$\begin{eqnarray}&&{F}_{\omega ,{k}_{h}}=\displaystyle \frac{\pi }{2}\displaystyle \frac{\omega }{{k}_{h}^{2}}{F}_{B}.\end{eqnarray} \tag{ 29 }$$

One can perform similar calculations including the results of thermal diffusivity (not shown here), and find the same relation between the boundary-forced and bulk-forced responses. We hypothesize that this relation is universal, and also carries over to the more realistic calculations with density variation, spherical geometry, etc. Then we have the following equation:

$$\begin{eqnarray}&&{F}_{\omega ,{\ell }}=\displaystyle \frac{\pi }{2}\displaystyle \frac{\omega {r}_{f}}{{\ell }({\ell }+1)}{u}_{r}({r}_{\mathrm{RCB}};\omega ,{\ell }).\end{eqnarray} \tag{ 30 }$$

Substituting into Equation (19), we find

$$\begin{eqnarray}\begin{array}{rcl}\delta {L}_{{\ell }}({R}_{\star },t) & = & {u}_{r}({r}_{\mathrm{RCB}};\omega ,{\ell })\left(2{\pi }^{2}\rho ({r}_{f}){r}_{f}^{3}\omega \right)\\ & & \times \,\left[\displaystyle \sum _{\omega ^{\prime} }\displaystyle \frac{i{\tilde{u}}_{h}^{\dagger }({r}_{f};\omega ^{\prime} )\delta {\tilde{L}}_{{\ell }}({R}_{\star };\omega ^{\prime} )}{\omega -\omega ^{\prime} }\right]\exp (-i\omega t).\end{array}\end{eqnarray} \tag{ 31 }$$

Finally, we have that the transfer function $T(f)\equiv T(\omega /2\pi )$, defined in Equation (1) of the main text, is given by

$$\begin{eqnarray}\begin{array}{rcl}T(\omega ) & = & \displaystyle \frac{1}{{\rm{\Delta }}r}{\displaystyle \int }_{{r}_{0}}^{{r}_{0}+{\rm{\Delta }}r}2{\pi }^{2}\rho ({r}_{f}){r}_{f}^{3}\omega \\ & & \displaystyle \times \,\left[\displaystyle \sum _{\omega ^{\prime} }\displaystyle \frac{i{\tilde{u}}_{h}^{\dagger }({r}_{f};\omega ^{\prime} )\delta {\tilde{L}}_{{\ell }}({R}_{\star };\omega ^{\prime} )}{\omega -\omega ^{\prime} }\right]{{dr}}_{f}.\end{array}\end{eqnarray} \tag{ 32 }$$

We calculate the transfer function for twelve stellar structure models generated by MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019). We use models with ZAMS masses of 3 M_⊙, 7 M_⊙, 10 M_⊙, and 20 M_⊙, at three different evolutionary stages during the main sequence: when their core H fraction is about 0.66 (near the ZAMS), 0.33, and 0.01 (TAMS). Models were calculated assuming a metallicity Z = 0.02. Convective boundaries were determined using the Schwarzschild criterion, with convective velocities calculated using the mixing length theory in the Henyey formulation with α_MLT = 1.6. Above the convective core, a step overshooting with α_ov = 0.29 was included. The models were computed using gold tolerances and the dedt-form of the energy equation, which provides superior results for energy conservation (Paxton et al. 2019). Details of the calculations, including a typical inlist, can be found at https://github.com/matteocantiello/gyre_igw.

We use GYRE (Townsend et al. 2018) to calculate non-adiabatic normal modes for each stellar model. We use a second-order Gauss–Legendre collocation scheme, which we found to give more robust eigenfunctions than higher-order schemes. The grid is calculated using n_inner=5, alpha_osc=20, and alpha_exp=4. Increasing the grid resolution did not change the transfer function at high frequencies, but leads to more numerical noise at very low frequencies. We search for modes using the grid_type=''INVERSE'' method, and pick the maximum frequency such that we find all of the highest-frequency g-modes, and the minimum frequency such that we can resolve the fast decay of the transfer function for strongly damped modes. Although p-modes are important for the (near) completeness of the eigenmode basis, we find that they have a negligible contribution to the transfer function. This is because ${\tilde{u}}_{h}^{\dagger }({r}_{f})$ is small for p-modes, which are localized near the surface of the star, and there is poor frequency matching, so $\omega -\omega ^{\prime}$ is large. All GYRE inlists are at https://github.com/lecoanet/massive_star_igw.

Once we have the g-mode eigenmodes, we calculate the dual basis according to Equation (10). We then calculate the surface luminosity perturbation for a single r_f. Although convection likely forces waves simultaneously over a range of radii, we use a delta-function forcing as we only study the propagation of the waves after they have been excited. To calculate the transfer function, we average over 50 equally spaced values of r_f with Δr = 0.02R_⋆. All our stellar models have sharp composition gradients just outside the convection zone. We find that if we pick r_f below this composition gradient, the transfer function is dominated by numerical noise. Thus, we pick r₀ = 0.25R_⋆ for all models except the 20 M_⊙ model near the ZAMS, for which we pick r₀ = 0.35R_⋆. Increasing r₀ above these values changes the overall normalization of the transfer function—waves experience less amplification if they are generated at larger radii—but not the shape.

At the lowest frequencies, we sometimes find the transfer function abruptly changes slope. This stems from numerical errors in the eigenfunctions and/or eigenvalues. The exact position of this region changes with any change in numerical method (e.g., changing the grid resolution, number of modes, or r₀). We expect this region would not be present in more accurate calculations.

The transfer function becomes extremely large near the standing-mode frequencies at high frequencies. To find these peaks, we initially calculate the transfer function for a grid of frequencies which are roughly equally spaced logarithmically. Then we recalculate the transfer function on a grid which is refined near the peaks. We recursively refine until the peaks of the transfer function are well resolved. It is important to resolve the peaks of the transfer function to ensure we correct measure the frequency-integrated power in each peak.

We now describe the solution of the forcing problem via numerical integration of the evolution equations using Dedalus. Here we force the waves with a bottom boundary condition, which is equivalent to the bulk forcing used in GYRE (Section B). There are two important reasons for performing these calculations using Dedalus. First, we can test the robustness of the GYRE calculation by comparing with a second calculation using a completely different method. Second, the Dedalus calculations make a different set of simplifying assumptions than the GYRE calculations, so that we can test the sensitivity to assumptions made when using GYRE. The GYRE calculations assume the g-modes form a basis for the traveling wave calculation, and they also force waves at r_f, which is above the radiative–convective boundary. In contrast, our calculations with Dedalus expand perturbations in Chebyshev polynomials, which are known to form a complete basis, and we force the waves at the radiative–convective boundary. On the other hand, some shortcomings of the Dedalus simulations include the following. (1) We solve the equations in the Cowling approximation. (2) Our simulations domain does not extent to r = 0 or r = R_⋆, because if it include a convection zone, simulations become dominated by exponentially growing convection modes. (3) Very long integrations are required to resolve the standing-mode peaks. (4) The Dedalus calculations are much slower than the GYRE calculations.

We use Dedalus to solve the linearized evolution equations in the Cowling approximation,

$$\begin{eqnarray}&&{\partial }_{t}{u}_{r}+{\partial }_{r}p+p\displaystyle \frac{d\mathrm{log}{\rho }_{0}}{{dr}}+\displaystyle \frac{{g}_{c}{\rho }_{c}{R}_{\star }}{{p}_{c}}{g}_{0}\rho =0,\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{\partial }_{t}{u}_{h}+\displaystyle \frac{p}{r}=0,\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{r}^{2}{\rm{\Upsilon }}-{\partial }_{r}({r}^{2}{u}_{r})-{\ell }({\ell }+1){{ru}}_{h}=0,\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{\partial }_{t}\rho +{u}_{r}\displaystyle \frac{d\mathrm{log}{\rho }_{0}}{{dr}}+{\rm{\Upsilon }}=0,\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&{\nu }_{T}T-\left({\nu }_{T}{{\rm{\nabla }}}_{\mathrm{ad}}+{{\rm{\Gamma }}}_{1}^{-1}\right)\displaystyle \frac{{T}_{0}{\rho }_{0}}{{p}_{0}}p+{T}_{0}\rho =0,\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&L-{L}_{0}\left[{\left(\displaystyle \frac{{{dT}}_{0}}{{dr}}\right)}^{-1}{\partial }_{r}T+3\displaystyle \frac{T}{{T}_{0}}-\rho \right]=0,\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}p+\displaystyle \frac{1}{{\rho }_{0}}\displaystyle \frac{{{dp}}_{0}}{{dr}}{u}_{r}+{c}_{s,0}^{2}{\rm{\Upsilon }}\\ \,+\,\displaystyle \frac{{{Cc}}_{s,0}^{2}{\nu }_{T}}{4\pi {r}^{2}{\rho }_{0}{T}_{0}{c}_{p}}\,\left({\partial }_{r}L+{L}_{0}\displaystyle \frac{{{dT}}_{0}}{{dr}}\displaystyle \frac{{\ell }({\ell }+1)}{{r}^{2}}T\right)=0.\end{array}\end{eqnarray} \tag{ 39 }$$

Here u_r is the radial velocity and u_h is the horizontal velocity given by ${{\boldsymbol{u}}}_{h}={u}_{h}{{\boldsymbol{\nabla }}}_{h}{Y}_{{\ell },m}(\theta ,\phi )$, and ϒ is the divergence of the velocity. The background density ρ₀, pressure p₀, gravitational acceleration g₀, temperature T₀, luminosity L₀, and specific heat at constant pressure c_p, are all taken from the MESA model, and are all normalized to their values at the bottom of the domain (at the radiative–convective interface). Their values at the bottom of the domain are denoted with subscript c. The thermodynamic quantities are defined as ${\nu }_{T}=-{(\partial \mathrm{log}\rho /\partial \mathrm{log}T)}_{p}$, ${{\rm{\nabla }}}_{\mathrm{ad}}\,={(\partial \mathrm{log}T/\partial \mathrm{log}p)}_{s}$, and ${{\rm{\Gamma }}}_{1}={(\partial \mathrm{log}p/\partial \mathrm{log}\rho )}_{s}$, where s is the entropy. Radii are normalized to the stellar radius R_⋆. p and ρ are the Eulerian pressure and density perturbation, both normalized to ρ₀. T and L are the Eulerian temperature and luminosity perturbations, normalized to T_c and L_c, respectively. The normalized adiabatic sound speed squared is ${c}_{s,0}^{2}={{\rm{\Gamma }}}_{1}{p}_{0}/{\rho }_{0}$. The constant C is given by

$$\begin{eqnarray}&&C=\displaystyle \frac{{L}_{c}}{{T}_{c}{c}_{p,c}\sqrt{{p}_{c}{\rho }_{c}}{R}_{\star }^{2}}.\end{eqnarray} \tag{ 40 }$$

We normalize time using

$$\begin{eqnarray}&&{t}_{{\rm{norm}}}={R}_{\star }\sqrt{\displaystyle \frac{{\rho }_{c}}{{p}_{c}}}.\end{eqnarray} \tag{ 41 }$$

For boundary conditions, we use

$$\begin{eqnarray}&&{u}_{r}({r}_{\mathrm{bot}})=\sin (\omega t)\displaystyle \frac{1}{2}\left[1+\tanh \left(\displaystyle \frac{t-{t}_{0}}{{\rm{\Delta }}t}\right)\right],\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&{u}_{r}({r}_{\mathrm{top}})=s({r}_{\mathrm{bot}})=s({r}_{\mathrm{top}})=0,\end{eqnarray} \tag{ 43 }$$

where $s=\rho -p/{c}_{s,0}^{2}$ is related to the entropy perturbation. We ramp up the boundary forcing starting at time t₀ ≈ 10P over a timescale Δt ≈ P, where P is the wave period.

We discretize the problem by representing all background and perturbation variables in terms of Chebyshev series. We split up the domain into three segments. The first segment comprises r between ≈0.216R_⋆, the top of the convection zone, and 0.5R_⋆, and is represented using 512 Chebyshev modes. The second and third segments comprise r between 0.5R_⋆ and 0.8R_⋆, and between 0.8R_⋆ and ≈0.976R_⋆, which is just below the surface convection zone; both are represented using 256 Chebyshev modes. Boundary conditions impose continuity of all perturbation variables at the matching points, 0.5R_⋆ and 0.8R_⋆. Background quantities are expanded in Chebyshev polynomials, but coefficients less than 10⁻⁵ of the maximum coefficient in each segment are set to zero to make the problem sparser. For timestepping, we use a four stage, third-order accurate implicit-explicit Runge–Kutta timestepping scheme (Ascher et al. 1997), with timestep ≈P/100. We also ran simulations with double the spatial resolution and smaller timestep size, and found neither led to significant changes in the transfer function.

To calculate the transfer function, we wish to evolve the oscillation equations until they reach a statistically steady state. Unfortunately, this takes a very long time, especially when the forcing frequency is close to the oscillation frequency of one of the higher-frequency standing modes. For low-frequency modes, we evolve the system for ≈100P, which long enough for the luminosity amplitude at the top boundary to saturate. For higher-frequency waves, we evolve for ≈550P. This integration time is sufficient, except for forcing frequencies very close to the oscillation frequencies of a standing mode. We have integrated as long as ≈10⁴P for some frequencies, but the luminosity at the top boundary still continues to grow with time. Thus, we can only provide lower limits on the transfer function near the four highest-frequency standing waves. We also only calculate the transfer function up to ≈1 days⁻¹. To calculate the transfer function, we calculate the root-mean-square luminosity perturbation at the top boundary (recall that the velocity perturbation at the bottom boundary has an amplitude of unity), over a portion of the integration time.

Because the Dedalus calculations are more computationally expensive, we only calculated the transfer function for the 10 M_⊙ stellar model near the ZAMS. In Figure 4 we compare our two methods for calculating the transfer function. We report the transfer function based off the luminosity in two locations: at r = 0.85R_⋆, and also at r ≈ 0.976R_⋆ (at the top of the Dedalus simulations). Note that in GYRE, the transfer calculation evaluated at r ≈ 0.976R_⋆ is nearly identical to the transfer calculation evaluated at the surface.

Zoom In Zoom Out Reset image size

The transfer functions measured with the two numerical techniques are very similar. At low frequencies, the transfer functions are very small, and there are regularly spaced peaks at high frequencies.

However, the peaks are at different frequencies in the two codes. We believe this is due to differences in boundary conditions. At the surface, GYRE assumes no lagrangian pressure fluctuations, and a perturbed Stefan–Boltzmann equations; in Dedalus we set the radial velocity and eulerian entropy perturbation to zero at r ≈ 0.976R_⋆. At the inner boundary, GYRE includes a convection zone, whereas the Dedalus domain does not. We do not expect the standing modes of the two problems to be the same, so the peaks in the transfer function should be at different frequencies. Note that there is always a Dedalus peak between every pair of GYRE peaks; this suggests that the period-spacing between the modes in the two calculations is very similar. Also recall that the amplitude of the four highest-frequency peaks in Dedalus are lower bounds on the transfer function, as we were not able to integrate long enough to establish the height of the peak. For this reason, we also did not calculate the transfer function in Dedalus for frequencies greater than ≈1 days⁻¹.

The two calculations' transfer functions also have different amplitudes. In both panels of Figure 4, we divided the Dedalus transfer function by six. With this normalization, we find excellent agreement between the amplitudes of the two signals at r = 0.85R_⋆. We expect the Dedalus transfer function to be larger because the waves in Dedalus are excited at r ≈ 0.216R_⋆, whereas in GYRE we excited waves between r = 0.25R_⋆ and r = 0.27R_⋆. However, the waves only amplify by a factor of ≈2 between these points. Future work should examine the remaining difference in normalization between the calculations. In this work, we do not use the transfer function in any way that depends on its overall normalization.

We believe the difference between the normalization of the transfer functions calculated at r ≈ 0.976R_⋆ is also due to differences in outer boundary conditions. It appears that the choice of boundary condition in Dedalus leads to luminosity fluctuations which are smaller by a factor of ≈3 at the upper boundary than the choice used in GYRE. We find that the eigenfunctions in both codes match very well until they reach the surface. This explains the agreement at r = 0.85R_⋆.

Although there are some differences between the GYRE and Dedalus calculations, we find them to be minor. In fact, it may even be surprising that the agreement is so good, given the completely different numerical methods and the various inconsistent assumptions made in the two calculations. This good agreement gives us confidence in the main features of the transfer function: the transfer function is small at low frequencies due to damping, and is characterized by regularly spaced peaks due to standing modes at high frequencies.

In Figure 5, we plot the transfer function, calculated using GYRE, for a range of stellar models of different ages and initial masses. We only plot the transfer function for ℓ = 1 so it is easier to identify the trends with mass and age. The ℓ = 1 modes are expected to dominate the disk-averaged luminosity variation, as higher ℓ modes have significant cancellation. The three ages we include correspond to core hydrogen fractions of 0.66 (near the ZAMS), 0.33 (labeled "mid"), and nearly zero (TAMS). In all cases, we find a smooth transfer function corresponding to traveling waves at low frequencies. These waves are strongly damped, so the transfer function is small. At higher frequencies, the transfer function is dominated by regularly spaced peaks corresponding to standing modes.

Zoom In Zoom Out Reset image size

More massive stars have fewer standing-mode peaks. This is because the luminosity of these stars is larger, so radiative damping is more important. Thus, damped traveling waves exist at higher frequencies for higher-mass stars. Nevertheless, there are always several weakly damping standing modes, at least up to 20 M_⊙. Lower-mass stars have an extremely large number of weakly damped standing modes. It is likely that ≈100 days TESS or K2 observations of 3 M_⊙ stars would not be able to resolve most the individual standing-mode peaks. However, lower-mass stars also have a larger Brunt–Väisälä frequency, so g-modes extend to higher frequency. We thus believe that standing modes peaks with frequencies ∼4–5 days⁻¹ should be detectable in observations of ≈3 M_⊙ stars, if such waves were efficiently excited by convection.

There are only minor changes in the transfer function as the stars evolve from the ZAMS to midway between ZAMS and TAMS. The frequencies of the highest modes increase modestly, and radiative damping becomes somewhat weaker, leading to more standing modes. However, we find substantial changes when the stars reach TAMS. Radiative damping becomes significantly weaker at TAMS, leading to many more standing modes. The frequency spacing between the modes also seems smaller. Thus, it may be more difficult to detect individual standing modes for stars at the TAMS.

Several of our transfer functions have an abrupt change of slope at ≈0.05 day⁻¹, e.g., ZAMS stars ≥7 M_⊙. This is a numerical effect; this feature changes magnitude and frequency range when we change the resolution of the modes, the number of modes used, etc. The Dedalus calculation shows the transfer function decreases smoothly past these frequencies (see Figure 4, left panel). We believe the uptick in the transfer function for our 7 M_⊙ TAMS model at the lowest frequencies is also due to this type of numerical error. For our 3 M_⊙ TAMS model, we were not able to calculate the transfer function at low frequencies due to repeated floating-point exception errors thrown by GYRE.

At high frequencies, the transfer functions are dominated by very sharp peaks corresponding to standing modes. For many stellar models, the transfer function between these peaks seems to be nearly constant. We do not understand this behavior, but note that it seems to be robust to changes in numerical parameters, e.g., the resolution of the modes, the number of modes used, etc. It is unlikely that the behavior of the transfer function between the peaks have any observational implications.

Finally, we note that our 10 M_⊙ TAMS model and 20 M_⊙ "mid" model have well-defined peaks at ≈0.04 day⁻¹. These appear to be numerically robust. We do not have a physical understanding of these features. If these features are physical, they may be detectable in long-duration observations of massive stars.

Both LQ and R spectra have the normalization (Lecoanet & Quataert 2013; Rogers et al. 2013)

$$\begin{eqnarray}&&{F}_{w}\sim { \mathcal M }{F}_{c},\end{eqnarray} \tag{ 44 }$$

where F_w is the total wave flux, and F_c is the convective flux. The efficiency of wave generation is given by ${ \mathcal M }$, the convective Mach number, which we calculate using ω_c/N₀, where ω_c is the convective frequency, and N₀ is a typical value of the Brunt–Väisälä frequency in the radiative zone. The wave flux is dominated by waves near the convective frequency ω_c and ΛH/r_RCB ∼ 1, where r_RCB and H are the radius and the scaleheight at the radiative–convective boundary, and ${\rm{\Lambda }}=\sqrt{{\ell }({\ell }+1)}$ as above. Thus, we have

$$\begin{eqnarray}&&{ \mathcal M }{F}_{c}\sim {\rho }_{\mathrm{RCB}}{{HN}}_{0}{u}_{r}^{2},\end{eqnarray} \tag{ 45 }$$

where u_r is the radial velocity of waves with frequencies near ω_c and ${\rm{\Lambda }}H/{r}_{\mathrm{RCB}}\sim 1$, and ρ_RCB is the density at the radiative–convective boundary. We define ${F}_{c}={\rho }_{\mathrm{RCB}}{u}_{c}^{3}$, and ${\omega }_{c}={u}_{c}/H$. Then we have

$$\begin{eqnarray}&&{u}_{r}\sim {u}_{c}{ \mathcal M }.\end{eqnarray} \tag{ 46 }$$

Recall that this estimate is only for the most energetic waves with frequencies near ω_c and ${\rm{\Lambda }}H/{r}_{\mathrm{RCB}}\sim 1$.

We can now combine this with the predicted frequency and wavenumber scalings (Section 2). The predictions for the radial velocity spectra are

$$\begin{eqnarray}&&{u}_{r}\sim {u}_{c}{ \mathcal M }{\left(\displaystyle \frac{{\rm{\Lambda }}H}{{r}_{\mathrm{RCB}}}\right)}^{-0.9}{\left(\displaystyle \frac{f}{{f}_{c}}\right)}^{0.4}\ (R\,\mathrm{spectrum}),\end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray}&&{u}_{r}\sim {u}_{c}{ \mathcal M }{\left(\displaystyle \frac{{\rm{\Lambda }}H}{{r}_{\mathrm{RCB}}}\right)}^{5/2}{\left(\displaystyle \frac{f}{{f}_{c}}\right)}^{-13/4}\ (\mathrm{LQ}\,\mathrm{spectrum}).\end{eqnarray} \tag{ 48 }$$

To evaluate these expressions, we calculate the average convective velocity using

$$\begin{eqnarray}&&{u}_{c}^{3}=\displaystyle \frac{{\int }_{0}^{{r}_{\mathrm{RCB}}}{L}_{c}{dr}}{{\int }_{0}^{{r}_{\mathrm{RCB}}}4\pi {r}^{2}\rho {dr}},\end{eqnarray} \tag{ 49 }$$

where L_c is the convective luminosity. The convective angular frequency is ${\omega }_{c}={u}_{c}/H$, and the convective frequency is ${f}_{c}={\omega }_{c}/(2\pi )$. We use $N({r}_{\mathrm{RCB}}+H)={N}_{0}$ as our estimate of the Brunt–Väisälä frequency.

We calculate each of these terms using our 10 M_⊙ MESA model near the ZAMS. We find ${u}_{c}\approx 0.01{R}_{\star }/d$, ${ \mathcal M }\approx 0.002$, ${r}_{\mathrm{RCB}}/H\approx 2.3$, and f_c ≈ 0.018 days⁻¹. Thus, the theoretical predictions for the radial velocities are

$$\begin{eqnarray}&&{u}_{r}\sim (2\times {10}^{-4}\,{R}_{\star }/d)\,{{\rm{\Lambda }}}^{-0.9}{f}^{0.4}\ (R\,\mathrm{spectrum}),\end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}&&{u}_{r}\sim (6\times {10}^{-12}\,{R}_{\star }/d)\,{{\rm{\Lambda }}}^{5/2}{f}^{-13/4}\ (\mathrm{LQ}\,\mathrm{spectrum}).\end{eqnarray} \tag{ 51 }$$

For comparison, in Figure 2, we use a multiplying prefactor of 5 × 10⁻⁸ R_⋆/d for the R spectrum, and 2 × 10⁻¹¹ R_⋆/d for the LQ spectrum. Thus, the R spectrum predicts surface variability that is many orders of magnitude higher than what is observed. It may seem that the LQ spectrum predicts a similar magnitude to the observed variability, even if the shape of the spectrum is very different. However, we believe this is likely coincidental, as the very steep frequency dependence of the LQ spectrum makes it difficult to predict the overall amplitude of the spectrum. For instance, if we had used a different definition of the convective frequency, leading to a change by 2π, the overall amplitude of the spectrum could be as low as 2 × 10⁻¹⁴ or as high as 3 × 10⁻⁹. However, such a change in the definition of the convective frequency would lead to a very minor change in the R spectrum, because it has a very weak dependence on f_c. This illustrates the need to calibrate the spectrum and convective frequency from 3D simulations of core convection.

We have not included the effects of rotation in any of our calculations. Rotation plays an important role in the convective excitation of internal waves (Mathis et al. 2014). First, rotation can strongly constrain the convection, leading to flows that mostly align with the rotation axis (e.g., Augustson et al. 2016). This will modify the wave excitation process. Second, rotation changes the internal waves themselves. The eigenvalue problem for oscillation modes is 1D without rotation, but becomes 2D with rotation. This is because rotation couples together modes with different spherical harmonic degree ℓ, but does not couple modes with different azimuthal degree m. Oscillation modes with rotation can have non-trivial latitudinal structure.

Because of the 2D nature of the oscillation problem with rotation, it is difficult to calculate the transfer function. However, we will comment on some expected properties of the transfer function with rotation. Internal waves cannot propagate at the poles for f < 2/T_rot, where T_rot is the rotation period (not to be confused with the transfer function, T). The waves concentrate closer and closer to the equator in the limit that fT_rot/2 goes to zero. These waves will cause negligible surface variability in stars observed from the poles. It is likely that at least some of the stars observed in B19 are observed nearly pole-on. These stars would show very little power at frequencies less than 2/T_rot if the low-frequency variability was due to waves. However, a sharp decline in variability at low frequencies due to rotation is not observed in any of the stars, again suggesting that the variability is not due to waves.

If rotation is sufficiently weak (i.e., f ≫ 2/T_rot), then oscillation modes are Doppler-shifted by rotation, $f\to f+m/{T}_{\mathrm{rot}}$, where m is the azimuthal degree of the mode. If the convection is approximately isotropic, energy will be distributed between the different m modes evenly. Thus, each peak in the non-rotating transfer function will be evenly split into 2ℓ + 1 peaks under the influence of weak rotation. This will cause more peaks in the spectrum, and will change the spacing between peaks.

To illustrate these two effects, we plot rotationally modified surface luminosity spectra in Figure 6. This requires choosing nominal rotation periods for massive stars. Dufton et al. (2013) and Ramírez-Agudelo et al. (2013) measured the rotation velocities of B & O type stars as part of the VLT-FLAMES Tarantula Survey. In both cases, they find rotation velocities around 100–300 km s⁻¹. Assuming a typical radius of 10R_⊙, these roughly correspond to rotation rates of 2–5 days. This is in line with Nielsen et al. (2013), who reported a typical photometric rotation period for B stars of around 4 days. With this in mind, we consider two possible rotation periods, 2 days and 5 days, to demonstrate the possible effects of rotation. For frequencies less than 2/T_rot (shaded in gray) we do not modify the spectrum. However, waves at these frequencies would be localized near the equator, and would have negligible surface manifestation for stars viewed pole-on. For frequencies greater than 2/T_rot, we split the power of each peak evenly into 2ℓ + 1 peaks with frequencies f + m/T_rot for $| m| \leqslant {\ell }$, assuming the weak-rotation limit. This leads to more peaks at high frequencies, and a more complicated spectrum. However, most of the individual peaks should still be resolved by the observations, so the theoretical spectra still do not resemble the smooth observed spectra of B19.

Zoom In Zoom Out Reset image size