# On the Origin of Stochastic, Low-Frequency Photometric Variability in
Massive Stars

Matteo Cantiello Center for Computational Astrophysics, Flatiron Institute,
New York, NY 10010, USA  Department of Astrophysical Sciences, Princeton
University, Princeton, NJ 08544, USA  Daniel Lecoanet Department of
Engineering Sciences & Applied Mathematics, Northwestern University, Evanston,
IL 60201, USA  CIERA, Northwestern University, Evanston, IL 60201, USA
Department of Astrophysical Sciences, Princeton University, Princeton, NJ
08544, USA  Adam S. Jermyn Center for Computational Astrophysics, Flatiron
Institute, New York, NY 10010, USA  Luca Grassitelli Argelander-Institut für
Astronomie der Universität Bonn, Auf dem Hügel 71, D-53121 Bonn, Germany
Matteo Cantiello<EMAIL_ADDRESS>
###### Abstract

High-precision photometric observations have revealed ubiquitous stochastic
low-frequency photometric variability in early type stars. It has been
suggested that this variability arises due to either subsurface convection or
internal gravity waves launched by the convective core. Here we show that
relevant properties of convection in subsurface convective layers correlate
very well with the timescale and amplitude of stochastic low-frequency
photometric variability, as well as with the amplitude of macroturbulence. We
suggest that low-frequency, stochastic photometric variability and surface
turbulence in massive stars are caused by the the presence of subsurface
convection. We show that an explanation for the observed surface photometric
variability and macroturbulence relying on convective core driven internal
gravity waves encounters a number of difficulties and seems unlikely to be
able to explain the observed trends.

convection, stars: magnetic field, stars: variability

††journal: ApJ

##  1 Introduction

Massive stars largely drive the dynamical and chemical evolution of gas in
galaxies (e.g. Hopkins et al., 2014). They accomplish this via their stellar
winds, eruptions, and explosive deaths, ultimately producing neutron stars and
black holes (Langer, 2012). These compact remnants can merge and generate the
gravitational waves observed by LIGO/Virgo (Abbott et al., 2016). The
evolutionary trajectory starting with a massive star burning hydrogen in its
core and ending with a compact remnant is understood only qualitatively. We
still do not know how to map initial properties of the star, like mass,
rotation rate, and metallicity, to e.g. the final mass and spin of the compact
remnant it leaves behind. The picture is further complicated by the fact that
the majority of massive stars are found in multiple systems (Sana et al.,
2012), with a large fraction expected to interact with their companions (de
Mink et al., 2014).

The path towards a detailed understanding of massive stars begins with a
quantitative study of their internal structure on the main sequence. In recent
years asteroseismology has opened a new window on these challenging
astrophysical environments, with high precision photometry from space
delivering many new exciting results (e.g., MOST, CoRoT, BRITE, Kepler/K2 and
TESS, see Bowman, 2020). The latest discovery is the detection of a new
ubiquitous phenomenon in massive stars: stochastic low-frequency photometric
variability (SLF variability; Blomme et al., 2011; Bowman et al., 2019a, b;
Pedersen et al., 2019; Bowman et al., 2020; Rauw & Nazé, 2021). This joins a
number of other surface and wind phenomena that are routinely observed in
early-type stars, including surface velocity fluctuations (Macroturbulence;
Simón-Díaz & Herrero, 2014), line profile variability (Fullerton et al.,
1996), and discrete absorption components in UV spectra (Howarth & Prinja,
1989; Cranmer & Owocki, 1996; Fullerton et al., 1997; Kaper et al., 1997).
Surface magnetism and bright spots are harder to observe but could still be
common in these stars (e.g. Ramiaramanantsoa et al., 2014).

The origin of this SLF variability is currently debated. It could be caused by
sub-surface convection zones (Cantiello et al., 2009; Blomme et al., 2011;
Lecoanet et al., 2019) or by internal gravity waves (IGWs) launched by
turbulent core convection (Edelmann et al., 2019; Ratnasingam et al.,
2020)111Classical heat-mechanism pulsations could be responsible for
spectroscopic and photometric variability in specific parts of the HRD, but
they can hardly justify the apparent ubiquity of macroturbulence and SLF in
massive stars (Godart et al., 2017; Simón-Díaz et al., 2017a).. Instabilities
in the stellar wind could also play a role (Krticka & Feldmeier, 2021).
Regardless of its origin, this photometric signal likely carries important
information about stellar structure, complementing asteroseismic studies that
use well-identified oscillation modes (e.g. Aerts, 2019; Aerts et al., 2019;
Bowman, 2020).

Recently the use of high resolution ground-based spectroscopy for the targets
observed by K2 and TESS (Bowman et al., 2020) has allowed precise
determination of stellar parameters, including spectroscopic mass, luminosity,
and macroturbulence (Burssens et al., 2020). The latter is particularly
important if, as seems likely, the mechanism exciting surface turbulent
velocities is the same as that which produces the observed SLF variability
(Grassitelli et al., 2016; Bowman et al., 2020).

Here we combine spectroscopic and photometric data to compare the observed
properties of the stochastic photometric variability and macroturbulence with
predictions from non-rotating 1D stellar models. We make simple predictions
for the amplitude and frequency of the variability that subsurface convection
induces at the stellar surface and examine how these vary with stellar
temperature and luminosity. We find that the predicted trends coming from a
suburface convection zone driven by the iron opacity peak (FeCZ) at $\approx
150$kK match the observations well.

We next show that one way to differentiate between the two proposed mechanisms
is by examining macroturbulence and SLF in massive stars with surface magnetic
fields, since magnetic effects have a larger impact on core IGWs than on the
FeCZ. Macroturbulence is observed in stars with fairly strong magnetic fields,
sufficient to suppress IGWs from the core, favoring a model based on
subsurface convection. Furthermore, the only stars with no observed
macroturbulence are ones where the magnetic field is strong enough to shut off
the FeCZ, and so far as we know all stars with such strong magnetic fields
lack macroturbulence, consistent with a subsurface origin of surface
perturbations (Jermyn & Cantiello, 2020).

Based on these considerations we suggest that subsurface convection represents
a possible unifying mechanism causing SLF variability, surface turbulence, and
magnetic spots in massive stars.

Figure 1:  Evolution of the normalized radial location of the FeCZ in a
20${\rm M}_{\odot}$ model during the main sequence (first 5 Myr of evolution).
The HeCZ is also visible very close to the surface. Locations of
P(r)/P($R_{*}$) = $e^{n}$ for $n=[1,3,5,7]$ are also shown and labelled as
$n{\rm H}_{\rm P}$.

##  2 Methods

We calculated stellar evolution models using the Modules for Experiments in
Stellar Astrophysics (MESA Paxton et al., 2011, 2013, 2015, 2018, 2019)
software instrument. Details on the microphysics inputs to this software
instrument are given in Appendix A. Our models have initial mass ranging from
5 to 120${\rm M}_{\odot}$ and are non-rotating. Since OB stars are known to be
rapidly rotating (e.g. Maeder & Meynet, 2000; Dufton et al., 2013; Ramírez-
Agudelo et al., 2013), we discuss the potential impact of rotation on our
results in Section 5. We also neglect the effect of wind mass-loss. While this
process is important for the evolution of massive stars (Smith, 2014), the
conditions in subsurface convective zones and the properties of convection
depend almost exclusively from the location in the Hertzsprung-Russell
diagram. This is true as long as the outer layers composition is not
substantially altered. We focused on the structure and convective properties
of our models in order to calculate typical frequencies and convective fluxes.
Convection is calculated in the framework of the Mixing Length Theory (Böhm-
Vitense, 1958, MLT hereafter), and we adopted $\alpha_{\rm MLT}=1.6$. While
the properties of efficient convection zones (e.g. stellar cores) are
insensitive to the choice of this parameter, those of inefficient convective
regions close to the stellar surface _do_ depend on $\alpha_{\rm MLT}$ (Jiang
et al., 2015; Cantiello & Braithwaite, 2019). In Section 6 we discuss how
uncertainties in the choice of the $\alpha_{\rm MLT}$ parameter affect our
results.

Since we are interested in the excitation of surface phenomena we focus on the
upper part of convective zones. Following Cantiello et al. (2009) we define
the average of a generic quantity $q$ as

$\overline{q}\equiv\frac{1}{\alpha_{\rm
MLT}\mathrm{H}_{\mathrm{P}}}\int_{\mathrm{R}_{\mathrm{c}}-\alpha_{\rm
MLT}\mathrm{H}_{\mathrm{P}}}^{\mathrm{R}_{\mathrm{c}}}q(r)\;dr,$ (1)

where $\mathrm{H}_{\mathrm{P}}$ is the pressure scale height calculated at the
upper boundary ($\mathrm{R}_{\mathrm{c}}$) of the convective zone of interest.
We tested a variety of different average prescriptions and found that our
results do not depend much on the specific choice of prescription.

Using equation (1) we extracted the average Mach number $\overline{\rm M}_{\rm
c}$, convective velocity $\overline{v}_{\rm c}$ and density
$\overline{\rho}_{\rm c}$ in the convective core and in the subsurface
convection zones of our models. The most important subsurface convection zone
for the massive stars we focus on is the FeCZ (see Fig. 1), although in the
low luminosity regime He-driven convection zones could play a role as well
(Cantiello et al., 2009; Cantiello & Braithwaite, 2019).

We compare relevant properties of our theoretical models with the observed
characteristic frequency $\nu_{\rm char}$ and amplitude $\alpha_{0}$ of SLF
variability. These quantities are derived by fitting the stochastically-
excited, low-frequency power excess in a power density spectrum using a
Lorentzian function (e.g. Bowman et al., 2019b):

$\alpha(\nu)=\frac{\alpha_{0}}{1+\big{(}\frac{\nu}{\nu_{\rm
char}}\big{)}^{\gamma}}\,+P_{\rm W}.$ (2)

This shows that $\alpha_{0}$ represents the amplitude at zero frequency.
$\nu_{\rm char}$ is defined as in eqn. 3 with $\tau$ the typical timescale of
the SLF variability. Finally, $\gamma$ is the gradient of the linear part of
the profile and $P_{\rm W}$ is a white noise term. Due to its stochastic, low-
frequency manifestation in the power density spectrum, in the literature the
SLF variability is also referred to as “red noise” (e.g. Blomme et al., 2011).

##  3 Results

We want to test a possible correlation between the properties of the FeCZ and
observed photospheric phenomena, in particular SLF variability and turbulent
velocities at the stellar surface (macroturbulence). We proceed by calculating
quantities that measure the amplitude of perturbations in the FeCZ. We then
check if some of these properties correlate with the amplitude of observed
photospheric phenomena, including turbulent velocity fluctuations and SLF
variability.

We define the characteristic frequency as

$\nu_{\rm char}\equiv\frac{1}{2\pi\tau},$ (3)

where $\tau$ is a characteristic timescale. For comparing with observations we
set $\tau$ to be the average convective turnover time

$\tau_{\rm c}=\alpha_{\rm
MLT}\overline{\mathrm{H}}_{\mathrm{P}}/\overline{v}_{\rm c},$ (4)

calculated either in the FeCZ or in the convective core.

We calculated the convective flux $F_{\rm c}=\overline{\rho}_{\rm
c}\overline{v}_{\rm c}^{3}$, where $\overline{\rho}_{\rm c}$ and
$\overline{v}_{\rm c}$ are the average density and convective velocity
calculated according to eqn. 1. We did this as a function of both mass and
evolutionary history for stars with initial mass ranging from 5 to
120$~{}M_{\odot}$. Here we present results for an initial metallicity of
Z=0.02, but in Appendix C we report results for model grids with Z=0.006 and
0.002 as well.

###  3.1 FeCZ and Macroturbulence

(a)  (b)

Figure 2:  Left: Average convective velocities in the FeCZ as a function of
the location of stellar models in the spectroscopic H-R Diagram
($\mathscr{L}\equiv{\rm T}^{4}_{\rm eff}/g$). We also show observed stars with
detected macroturbulent velocity as grey circles. The area of the circle is
proportional to the observed macroturbulence (data from Burssens et al., 2020;
Bowman et al., 2020). The FeCZ is absent in models with
$\log\mathscr{L}/\mathscr{L}_{\odot}\lesssim 2.5$. Right: Same as left, but
showing predictions for the maximum of the ratio between turbulent pressure
and total pressure in any subsurface convection zone. The FeCZ largely
dominates, except for stars at low luminosities where turbulent pressure is
provided by a helium convection zone (HeCZ).

Macroturbulence is a spectroscopic measure of velocity fields with a scale
larger than the photons mean free path in the stellar atmosphere. The shape of
spectral lines can be used not only to measure the amplitude of the velocity
field, but also to infer some of its directionality. Note that Simón-Díaz &
Herrero (2014) showed that the line profiles are fitted better by a radial-
tangential velocity function than a Gaussian one, but the observations do not
tell if the dominant velocity component is radial or horizontal222The claim is
that the line profiles are fitted better by the radial-tangential velocity
function, compared to an isotropic one, with either $v_{\rm r}\gg v_{\rm t}$
or $v_{\rm r}\ll v_{\rm t}$, where $v_{\rm r}$ and $v_{\rm t}$ are the radial
and tangential components of the velocity, respectively (Simón-Díaz, private
comm.). Below we compare the observed amplitude of macroturbulent velocities
with predicted properties of the FeCZ.

We show the average convective velocity in the FeCZ in Fig. 2, left panel.
Velocities of the order $10\dots 100\ \rm kms^{-1}$ are found across the upper
spectroscopic H-R diagram ($\mathscr{L}\equiv{\rm T}^{4}_{\rm eff}/g$; Langer
& Kudritzki, 2014), with a trend of increasing $v_{\rm c}$ for higher
luminosities. Note that, contrary to core convection, the FeCZ is just mildly
subsonic, with Mach numbers ranging from 0.01 to 0.3, though these are
uncertain by a factor of $\approx 8$ due to a dependence on the uncertain
$\alpha_{\rm MLT}$ (Cantiello & Braithwaite, 2019).

In the right panel of the same figure we also show the maximum of the ratio of
turbulent pressure to total pressure (P${}_{\rm turb}\propto v_{\rm c}^{2}$),
which, in agreement with the results of Grassitelli et al. (2015b), shows a
strong correlation with the spectroscopically-derived macroturbulent
velocities in massive stars.

The two quantities in Fig. 2 measure the strength of the inefficient
convection. We do not yet know the exact mechanism connecting the FeCZ with
the surface velocity perturbations. If one assumes convective elements
conserve their inertia as they reach layers stable against convection (e.g.,
because they are thermally diffusive, see Jiang et al., 2015), then the
surface velocities should be proportional to the convective velocity; see the
left panel of Fig. 2. Alternatively, one can assume IGWs are excited with
pressure perturbations $\delta p\sim P_{\rm turb}\propto v_{\rm c}^{2}$, the
turbulent pressure of the convection (Press, 1981). This stochastic excitation
can lead to both running waves and standing modes. Using the polarization
relations of adiabatic IGWs (e.g., Sutherland, 2010), a wave with pressure
perturbation $\delta p$ has an associated horizontal velocity
$u_{h}\sim(\delta p/\rho_{0})k_{h}/\omega$, where $k_{h}$ is the horizontal
wavenumber of the wave and $\omega$ is its frequency. The dominant waves will
have $k_{h}\sim 1/{\mathrm{H}}_{\mathrm{P}}$ and $\omega\sim 1/\tau_{c}$
(Cantiello et al., 2009). So running IGWs would have surface velocities which
also scale like $u_{h}\sim v_{\rm c}$ at the surface. Concerning mode
excitation, Grassitelli et al. (2015b) argue that the ratio of turbulent
pressure to total pressure in the FeCZ traces the stochastic Lagrangian
pressure perturbation associated with the convective motions. This is
responsible for local deviation from hydrostatic equilibrium and the
excitation of high-order pulsations with frequencies close to the spectrum of
the fluctuations. The ratio of the turbulent to total pressure in the FeCZ is
reported in the right panel of Fig. 2.

We confirm that convective velocities and the ratio of turbulent pressure to
total pressure in the FeCZ correlate very well with the amplitude of
macroturbulence.

###  3.2 FeCZ and Stochastic, Low-Frequency Variability

(a)  (b)

Figure 3:  Left panel: Characteristic frequency $\nu_{\rm char}$ in the FeCZ
as a function of the location of stellar models in the spectroscopic H-R
Diagram (black contour lines). Evolutionary tracks are shown as grey solid
lines. We also show the observed stars with detected stochastic photometric
variability as grey circles. The area of the circle is proportional to the
observed $\nu_{\rm char}$, derived from eqn. (2) fitting the data in the range
$0.1\leq\nu\leq 360\,$d-1 (Burssens et al., 2020; Bowman et al., 2020). Right
panel: Ratio of FeCZ convective flux to the total stellar flux in the
spectroscopic H-R Diagram ($\mathscr{L}\equiv{\rm T}^{4}_{\rm eff}/g$). We
also show the observed stars with detected stochastic, low-frequency
photometric variability as grey circles. The area of the circle is
proportional to $\alpha_{0}$.

If the SLF variability is caused by the FeCZ, a natural choice of proxy for
the typical timescale of this variability is the convective turnover
timescale. We find typical values of the convective turnover timescale to be
about $\sim 0.1\dots 2$ d in the FeCZ, with a tendency for shorter values in
models with high effective temperature and surface luminosity (Cantiello et
al., 2009). We computed $\nu_{\rm char}$ using eqn. 3 and plot this alongside
observed characteristic frequencies333These were inferred via eqn. (2). of SLF
variability on the spectroscopic H-R diagram, see left panel in Fig. 3.

We see good agreement, with our models reproducing both the typical values of
the observed $\nu_{\rm char}$ and the trend with $\log\mathscr{L}$ and
$\log{\rm T}_{\rm eff}$. While the characteristic frequencies found by Bowman
et al. (2020) seem to be larger by a factor of $\approx$ 3, our predictions
for the turnover timescale are affected by uncertainty in the convective
velocities arising from the MLT treatment ($v_{\rm c}\propto\alpha_{\rm
MLT}^{3}$, so $\nu_{\rm char}\propto\alpha_{\rm MLT}^{2}$), as well as our
limited knowledge of the frequency spectrum generated by turbulent convection.
Since $\alpha_{\rm MLT}$ is uncertain by a factor of 2 or so, our estimates of
the characteristic frequency in the FeCZ are uncertain by a factor of
$\approx$ 4 and so are consistent with the observations.

We also compare the amplitude of SLF variability with the relative convective
flux in the FeCZ (right panel in Fig. 3). Values of Fc/F∗ tend to increase
with increasing $\log\mathscr{L}$ and decreasing $\log{\rm T}_{\rm eff}$, a
trend that is also found for the amplitude of observed SLF variability.

Overall the turnover timescale and relative flux in the FeCZ correlate very
well with the observed timescale and amplitude of SLF variability.

###  3.3 Core Convection and Stochastic, Low-Frequency Variability

(a)  (b)

Figure 4:  Left panel: Square root of the ratio between the flux of gravity
waves launched by the convective core and
$F_{0}=1/2\rho(r)\,r^{3}\omega^{2}\sqrt{N(r)^{2}-\omega^{2}}$, evaluated at
the stellar surface. In the absence of damping, this quantity is expected to
be proportional to the relative radial displacement (Appendix B) and hence to
the relative luminosity fluctuations at the surface. The IGW flux is
calculated multiplying the core convective flux by the average convective Mach
number in the top pressure scale height of the convective region. We also show
the observed stars with detected stochastic, low-frequency photometric
variability as grey circles. The area of the circle is proportional to
$\alpha_{0}$, derived from eqn. (2) fitting the data in the range
$0.1\leq\nu\leq 360\,$d-1 (Burssens et al., 2020; Bowman et al., 2020). Right
panel: Characteristic frequency $\nu_{\rm char}$ in the convective core as
function of the location of stellar models in the spectroscopic H-R diagram.
We also show the observed stars with detected photometric variability as grey
circles. The area of the circle is proportional to $\nu_{\rm char}$.

Early-type stars posses convective cores during their main sequence. Turbulent
convection in these convective cores excites internal gravity waves (Goldreich
& Kumar, 1990; Rogers et al., 2013; Lecoanet & Quataert, 2013; Shiode et al.,
2013; Edelmann et al., 2019; Horst et al., 2020), that can propagate through
the stellar envelope and reach the stellar surface (e.g. Ratnasingam et al.,
2019; Lecoanet et al., 2019; Ratnasingam et al., 2020). Some groups have
argued that such waves could be responsible for both the observed
macroturbulence and SLF variability in early-type stars (Aerts et al., 2009;
Simón-Díaz et al., 2017b; Bowman et al., 2019b, a).

There is substantial uncertainty in the surface brightness fluctuations from
IGWs generated by core convection. Shiode et al. (2013) predicted the typical
amplitude of g-modes excited by convection to be $\approx
10^{-2}-10^{2}\,\mu{\rm mag}$, seemingly at odds with the relatively large
amplitudes $\approx 10-10^{4}\mu{\rm mag}$ observed by e.g. Bowman et al.
(2020). However, Lecoanet et al. (2021) recently found an error in the Shiode
et al. (2013) prediction. The g-mode amplitude should be larger by a factor of
$\approx\sqrt{\nu/\gamma}$, where $\gamma$ is the mode’s damping rate. This
missing factor could increase the predicted g-mode amplitude by a factor of
$10^{4}$ or larger for high-frequency modes near the Brunt-Väisälä frequency
(e.g., $\sim 10\,d^{-1}$ for a $10M_{\odot}$ ZAMS star), but does not change
the predicted mode amplitudes for lower frequency waves (e.g., $\sim
0.3\,d^{-1}$ for a $10M_{\odot}$ ZAMS star). Lecoanet et al. (2019) argued
that there should be very low wave power at low frequencies due to radiative
damping, while the wave signal at frequencies above $0.5\,d^{-1}$ should be
dominated by g-modes, as predicted by Shiode et al. (2013). These features do
not seem to be present in the observed spectra.

Recent numerical simulations of wave generation by convection in a
$3M_{\odot}$ star (Edelmann et al., 2019; Horst et al., 2020) produce wave
fluctuation spectra which are qualitatively similar to those observed.
However, those simulations artificially boost the stellar luminosity by
factors ranging from $10^{3}$ to $10^{7}$. Boosting the luminosity should both
enhance the wave amplitude and increase the typical frequency of excited
waves, making it difficult to quantitatively compare to observations.

Although the detailed physics of wave generation by convection is uncertain,
we can still analyze the properties of core convective of our models. If the
surface variability is due to core convection, one would expect the
characteristic frequency and amplitude of the SLF variability to correlate
with the properties of the core convection. The flux of IGWs excited by
turbulent convection at the core boundary is of order F${}_{\rm IGW}={\rm
F}_{*}\overline{\rm M}_{\rm c}$ (Goldreich & Kumar, 1990), where we evaluate
the average Mach number $\overline{\rm M}_{\rm c}$ using eqn. 1 and the local
adiabatic sound speed. We expect the luminosity fluctuations to be
proportional to the relative surface radial displacement $\xi_{r}/R$ produced
by these waves at the stellar surface (e.g. Dziembowski, 1977; Aerts et al.,
2010). It can be shown that $\xi_{r}/R\propto\sqrt{{\rm F}_{\rm IGW}/{\rm
F}_{0}}$ (see Appendix B). This quantity is presented in Fig. 4 along with the
observed amplitudes of SFL variability.

One important caveat is that this estimate neglects the important role of
radiative damping, which is essential in shaping both the amplitude and the
shape of the spectrum of waves at the surface (e.g. Lecoanet et al., 2019).
Rotation might also be key in setting the amplitude and shape of the surface
fluctuations spectrum (See Section 5). Despite neglecting the important effect
of radiative damping, the trend for the relative radial displacement at the
surface do show some interesting correlations with the observed trend in SLF
variability, though they generally proceed the wrong way, with increasing
amplitude towards higher luminosities where we expect IGW to show the smallest
effects (left panel in Fig. 4).

Next, we focus on the characteristic timescale of waves excited by core
convection. We expect the maximum of the IGW flux to be launched at
frequencies close to $\nu_{\rm char}=(2\pi\tau_{\rm c})^{-1}$, with $\tau_{\rm
c}$ defined in eqn. 4. Fig. 4 shows that the predicted values of $\nu_{\rm
char}$ are in the range $0.02\dots 0.008\,{\rm d}^{-1}$ ($0.2\dots
0.08\,\mu\textrm{Hz}$), in agreement with the results of Shiode et al. (2013).
These values are about 2 orders of magnitude smaller than the typical
characteristic frequencies of SLF variability observed by (Bowman et al.,
2019b, a).

As mentioned earlier, a significant caveat in correlating core and surface
quantities is that wave propagation through the stellar envelope affects the
spectrum, changing the frequency of maximum power of waves reaching the
surface (Lecoanet et al., 2019). Radiative damping efficiently suppresses low-
frequency IGWs (because these waves have high radial wavenumbers), and the
peak of the spectrum observed at the surface is expected to move to higher
frequencies. The amplitude of this effect depends on the envelope properties,
which are a function of $\log{\rm T}_{\rm eff}$ and $\log\mathscr{L}$.
Nevertheless we will assume that trends of $\nu_{\rm char}$, as function of
$\log{\rm T}_{\rm eff}$ and $\log\mathscr{L}$, are still set by the core
convective properties. Under this assumption, Fig. 4 shows that $\nu_{\rm
char}$ should increase with both $\log{\rm T}_{\rm eff}$ and
$\log\mathscr{L}$, so that the characteristic frequencies of peak IGWs should
increase as stars evolve on the main sequence. This is exactly the opposite of
what is observed: the characteristic frequencies of SLF variability is largest
for stars early on on the main sequence, and seem to decreases as stars
evolve. Therefore either radiative damping is able to revert this trend or
else the observed variability is unlikely to be caused by IGWs launched by the
core.

##  4 Macroturbulence in Magnetic Stars

The hypothesis that subsurface convection is responsible for surface
turbulence is corroborated by the match between trends in the FeCZ properties
(i.e. convective velocities and turbulent pressure) and the observed micro and
macroturbulence (Cantiello et al., 2009; Grassitelli et al., 2015b). One
important test to this hypothesis is provided by the disappearance of
macroturbulence in stars with surface magnetic fields above a critical
strength (Sundqvist et al., 2013), closely corresponding to the critical field
needed to stabilize the FeCZ (MacDonald & Petit, 2019; Jermyn & Cantiello,
2020). Note that in OB stars, the FeCZ is deeper and more vigorous than the H
and He convection zones, so a magnetic field stabilizing the FeCZ will
necessarily also stabilize the other subsurface convection zones.

We can analogously define a critical magnetic field strength $B_{\rm crit}$
which suffices to reflect IGWs before they reach the photosphere. Using the
dispersion relation for IGWs in a magnetic medium, the radial component of
this field is (Fuller et al., 2015)

$\displaystyle B_{r,\rm crit}=\frac{\omega}{2k_{r}}\sqrt{4\pi\rho},$ (5)

where $\omega=2\pi\nu$ is the angular frequency of the waves, $\rho$ is the
density, and $k_{r}$ is the radial wave-number. This critical field may be
thought of as the field strength at which the Alfvén frequency computed with
the length-scale $1/k_{r}$ is comparable to the wave frequency and is
therefore analogous to the effect of rotation, which enters in when the
rotation angular velocity is faster than the wave frequency. For IGWs the
radial wave-number is related to the spherical harmonic degree $\ell$ by

$\displaystyle
k_{r}\approx\frac{\sqrt{\ell(\ell+1)}}{r}\left(\frac{N}{\omega}\right),$ (6)

where $N$ is the Brünt-Väisälä frequency. So

$\displaystyle B_{r,\rm crit}$
$\displaystyle\approx\frac{\omega^{2}r}{N}\sqrt{\frac{4\pi\rho}{2\ell(\ell+1)}}$
(7) $\displaystyle\approx\frac{\omega^{2}r}{\ell N}\sqrt{4\pi\rho}.$ (8)

Note that because this decreases with increasing $\ell$, all waves of a given
frequency are reflected if the $\ell=1$ waves are reflected.

We computed $B_{r,\rm crit}$ for several values of $\ell$ and $\nu$ for a
main-sequence model of a $30M_{\odot}$ star as a function of radius, shown in
Fig. 5. For frequencies comparable to those of core convection the critical
magnetic field strength is of order $10^{-2}\,\mathrm{G}$ to
$10^{-1}\,\mathrm{G}$. Strong macroturbulence is observed in similar-mass
O-type stars with magnetic fields up to $2.5\,\mathrm{kG}$ (e.g. HD 191612,
Sundqvist et al., 2013), so macroturbulence in those stars is unlikely to be
due to IGWs coming from their cores if those waves have similar frequencies to
that of core convection.

For frequencies comparable to the observed $\nu_{\rm char}$, on the order of
$3\,\mathrm{d}^{-1}$, the critical magnetic field is much larger, on the order
of $300\,\mathrm{G}$ to $1\,\mathrm{kG}$. This makes an explanation of
macroturbulence based on IGWs marginally inconsistent with observations
showing strong macroturbulence up to field strengths of $2.5\,\mathrm{kG}$.
However, because eqn. 7 is a strong function of $\omega$, and hence of $\nu$,
it is possible that these strongly magnetized stars just have larger $\nu_{\rm
char}$. The full range of observed characteristic frequencies spans
$0.2-10\,\mathrm{d^{-1}}$, corresponding to critical field strengths at the
surface of $3\,\mathrm{G}-10\mathrm{kG}$ for the $\ell=1$ mode. Thus while
IGWs with lower frequencies are inconsistent with observations of strongly-
magnetized stars with substantial macroturbulence, those at higher frequencies
likely make it to the surface and could contribute to the observed
macroturbulence.

A further prediction of this calculation is that, if IGWs are the cause of
macroturbulence in these stars, we should expect the strength of
macroturbulence to decline with increasing magnetic field strength as more and
more modes are reflected before they reach the surface. Sundqvist et al.
(2013) find no such trend, though they do that macroturbulence vanishes when
the magnetic field exceeds the FeCZ shutoff strength (MacDonald & Petit,
2019), and this is consistent with observations of other strongly-magnetized
O/B stars such as HD 215441 (Landstreet et al., 1989) and HD 54879 (Castro et
al., 2015), both of which show little or no macroturbulence and magnetic
fields stronger than the theoretical shutoff field strength. This again points
against an explanation of macroturbulence based on core-generated IGWs.

Figure 5:  The critical radial magnetic field strength needed to reflect
internal gravity waves of different $\ell$ and $\nu$ is show as a function of
fractional depth for a $30M_{\odot}$ stellar model at an age of
$4\,\mathrm{Myr}$. The frequency $\nu_{\rm char,MLT}$ is that given by eqn. 3
but evaluated for the core convection zone instead of the subsurface FeCZ. The
frequency $\nu_{\rm char,Obs}$ is a typical value of $3\,\mathrm{d}^{-1}$ for
the observations. We further show frequencies of $0.2\,\mathrm{d^{-1}}$ and
$10\,\mathrm{d^{-1}}$ as these span the full range of observed characteristic
frequencies.

##  5 Rotation

###  5.1 Subsurface Convection

The sample of stars discussed here have projected rotational velocities
$v_{\rm eq}\sin i$ in the range 7 to 320 $\rm kms^{-1}$, suggesting rotation
could have an impact on the properties of subsurface convection. The effect of
rotation is usually measured by the convective Rossby number
$R_{0}=1/(2\,\Omega\,\tau_{c})$, where $\Omega$ is the stellar rotational
frequency. For $R_{0}>1$ we expect rotation to have a moderate to negligible
impact. On the other hand when the rotational period becomes comparable or
shorter than the convective turnover time ($R_{0}\leq 1$), the properties of
convection can be altered substantially (e.g. Stevenson, 1979; Augustson &
Mathis, 2019). The typical rotational period of OB stars is about 3 days
(assuming a typical equatorial rotational velocity $v_{\rm eq}\approx 150\rm
kms^{-1}$), while the convective turnover timescale in the FeCZ is a few hours
(See e.g. Fig. 3). Therefore the $R_{0}$ in the observed sample is likely in
the range $1\dots 10$. For these values, the convective velocities as
calculated in the 1D MLT approximation are affected only at the $\sim 10\%$
level (see e.g. Fig.4 in Cantiello et al., 2009). This said, the latitudinal
structure of the FeCZ zone is substantially altered at the highest rotation
rates (Maeder et al., 2008), which could have an impact on the way these
regions affect the stellar surface.

###  5.2 Core Convection

In the stellar cores of intermediate and massive stars $R_{0}$ is very likely
$<1$, so rotation is expected to change the properties of convection (e.g.
Stevenson, 1979; Augustson & Mathis, 2019). At the same time, in the presence
of rotation, gravity waves can be perturbed by the Coriolis acceleration and
combine with inertial waves (gravito-inertial waves, GIWs). The stochastic
excitation of gravity and GIWs by rotating convective zones was studied by
Mathis et al. (2014) and Augustson et al. (2020). The main result is that
rotation can enhance the amplitude of stochastically excited waves (Mathis et
al., 2014). The work of Neiner et al. (2020) shows that in some rapidly
rotating stars, stochastically excited GIWs from the core could explain part
of the observed low-frequency variability.

However, we note that it is unlikely that core-generated GIWs are responsible
for the ubiquitous SLF variability. This is because the visibility of GIWs
depends on the inclination of the stellar rotation axis respect to the
observer. Internal waves cannot propagate at the poles for $\omega<2\Omega$,
where $\Omega$ is the stellar rotational frequency and $\omega$ is the wave
frequency. The largest wave flux is expected at low latitudes, with the degree
of equatorial confinement proportional to $\Omega$. This means that the
propagation domain of subinertial ($\omega<2\Omega$) GIWs excludes the pole
and it becomes increasingly concentrated toward the equator for faster
rotation rate (Dintrans & Rieutord, 2000; Prat et al., 2016; Augustson et al.,
2020). Assuming that the stars in Bowman et al. (2020) have spin vectors
randomly oriented, some should be observed nearly pole-on. Then if the SLF
variability was due to GIWs, these objects would show very little power at
frequencies less than $2\Omega$. Such a sharp decline in variability at low
frequencies is not observed in any of the stars, suggesting GIWs are not the
culprit (Lecoanet et al., 2019). A detailed study of the surface amplitude of
waves excited by core convection in rotating early-type stars is beyond the
scope of this paper.

##  6 Discussion

The presence of a subsurface convection zone can result in variability of
photospheric properties via a number of processes, including wave excitation
(Cantiello et al., 2009; Grassitelli et al., 2015b) and magnetic buoyancy
(Cantiello & Braithwaite, 2011). A linear perturbative analysis is limited,
since the turbulent fluctuations in these convective regions can be large
(Grassitelli et al., 2015b). Multi dimensional simulations including radiation
have been performed in a restricted range of the parameter space, and show
that the full turbulent manifestation of these convective regions extends up
to the stellar surface (Jiang et al., 2015, 2017, 2018; Schultz et al., 2020).
In their calculations of OB stars envelopes including the stellar photosphere,
Jiang et al. (2015) observe turbulent velocities reaching the isothermal sound
speed ($\approx 50\rm kms^{-1}$) at the stellar surface, demonstrating that
velocity fields of amplitude comparable to the observed macroturbulence are
naturally explained by the presence of the FeCZ.

Therefore it could be that the observed SLF variability and macroturbulence
simply represent the direct manifestation of turbulent, radiation-dominated
convection at the stellar surface. The simplified MLT treatment in our one
dimensional calculations is unable to capture the complex phenomenology of
these layers, but the fact that it can reproduce both the timescales and the
trends in amplitude of the observed SLF variabilty is compelling. It calls for
extending the radiation hydrodynamics simulations to cover the parameter space
of the TESS observations, in order to unravel the precise mechanism connecting
subsurface convection zones to the observed surface variability.

###  6.1 Perturbation Lengthscale

If the perturbation is due to stochastically excited modes driven by the FeCZ
(Grassitelli et al., 2015b), then we expect the largest fluctuations to be
produced by modes with $\ell\lesssim 20$ (see e.g. Godart et al., 2017). The
situation is different if instead the perturbation is provided by running
waves or by convective motions extending to the stellar surface (Cantiello et
al., 2009; Jiang et al., 2015). In this case a good proxy for the typical
scale of surface perturbations is provided by the size of convective cells in
the FeCZ. This in turn is quantified by the average pressure scale height in
the subsurface convection zone. Note that the pressure scale height only
decreases slightly moving from the FeCZ to the stellar surface. In general one
expects that velocity perturbations induced by the FeCZ should have scales
that are comparable or slightly larger than the line forming region, so
macroturbulence can be explained via this mechanism. Rotation could also be
responsible for organizing the convective flow on slightly larger scales (see
Sec. 5). Since convective turbulence also results in smaller-scale motions,
the FeCZ could also be responsible for the excitation of surface
microturbulence (Cantiello et al., 2009).

We show in Fig. 6 the number of convective cells $N_{\rm
CC}=(R_{\star}/\overline{\rm H}_{\rm P})^{2}$, calculated using the stellar
radius and the average pressure scale height in the FeCZ. We expect
approximately $10^{2}...10^{4}$ convective cells in the FeCZ of OB stars.
Therefore the order of the perturbation $\ell\approx\sqrt{N_{\rm CC}}\approx
10...100$. While it might seem impossible for such high-degree perturbations
to leave a visible signature on the stellar disc integrated properties, we
point out that in this case the surface fluctuations are uncorrelated. So even
high degree $\ell$ fluctuations do not undergo the dramatic cancellation
effects experienced by highly-correlated stellar oscillations. Similarly to
granulation, we expect the amplitude of the integrated surface perturbations
to scale as $1/\sqrt{N_{\rm CC}}$.

It is important to distinguish the signal induced by (sub)surface convection
zones in OBA stars with the granulation pattern expected in cool stars with
convective envelopes. The driving mechanism, properties, and location of these
convection zones change substantially from late- to early-type stars
(Cantiello & Braithwaite, 2019). This is why is not surprising that the
characteristic frequency of the SLF variability in the early-type stars
observed by Bowman et al. (2019b) and Szewczuk et al. (2021) does not follow
the granulation scaling of Kjeldsen & Bedding (1995), which was derived for
late-type stars. On the other hand we have shown here that the observed
frequencies are consistent with the expectation of perturbations arising from
subsurface convection zones (see Fig. 3).

Figure 6:  Number of convective cells in the FeCZ. This was calculated using
the average pressure scale height in the FeCZ.

###  6.2 Metallicity

If the FeCZ is responsible for surface macroturbulence and SLF variability,
then these phenomena should be affected by the stellar metallicity. On the
other hand we do not expect a significant metallicity dependence in the
context of a core-convection origin. The results of (Bowman et al., 2019b)
show that SLF variability is also observed in low-metallicity, LMC stars. The
presence of the FeCZ depends on the luminosity and metallicity of the star,
and in 1D stellar evolution calculations the FeCZ occurs above a luminosity of
L$\approx 10^{3.2}L_{\odot}$ for Z=0.02. This corresponds to a zero age main
sequence star of about 7${\rm M}_{\odot}$ (see e.g. Fig. 2). The LMC has a
metallicity about half solar, and in models at Z=0.008 the FeCZ only appears
at L$\approx 10^{3.9}L_{\odot}$, corresponding to a zero age main sequence
star of about 11${\rm M}_{\odot}$ (Cantiello et al., 2009). In our models with
metallicity Z=0.006 (See Appendix C) the FeCZ appears at
$\log\mathscr{L}/\mathscr{L}_{\odot}\approx 3.2$.

Note that these limits could move downward in the presence of atomic diffusion
and radiative acceleration (Richer et al., 2000), or due to an upward revision
of the uncertain values of Fe opacity (e.g. Bailey et al., 2015). We notice
that the TESS observations of LMC stars reported by Bowman et al. (2019a)
reveal a trend of lower $\nu_{\rm char}$ compared to the galactic sample (See
their Fig. 4). Interestingly, this trend is reproduced by the characteristic
frequency of convection in the FeCZ in our models (compare Fig. 3 with Fig.
8). It will be interesting to see if the TESS observations can probe the
transition region between stars with and without a FeCZ, and determine a
possible change in surface properties. At the same time, X-shooter within the
Ulysses program could more firmly establish the metallicity-dependence of the
macroturbulent line-broadening (Penny & Gies, 2009).

It is important to note that some low-luminosity, main sequence A stars in
Bowman et al. (2019b) show SLF variability. For these stars models do not
predict the presence of the FeCZ. However, these stars still show (sub)surface
convection zones triggered by ionization of H and He (Cantiello & Braithwaite,
2019). The characteristic frequencies of convection in these regions are also
in the range $\sim$tens of $\mu$Hz, although the amplitude of the velocity
fluctuations they can induce is much smaller than for the FeCZ, so it is not
clear if they could be linked to observed SLF variability and macroturbulence.
On the other hand, since the relative perturbation induced by subsurface
convection is weakest in the regime of A and late-B type (Cantiello &
Braithwaite, 2019; Jermyn & Cantiello, 2020), these stars are the best targets
for detecting core-generated IGWs. Depending on the amplitude of core
generated IGWs, the impact of subsurface convection could well be subdominant
in these objects, allowing for a detection.

###  6.3 Stochastic, Low-Frequency Variability in Evolved Stars

SLF variability with similar properties to the massive main sequence stars
discussed by (Bowman et al., 2019b, 2020) was recently observed in evolved
massive stars. Naze et al. (2021) detected SLF variability in both luminous
blue variables (LBV) and Wolf-Rayet (WR) stars, and Dorn-Wallenstein et al.
(2020) found the same photometric signature in yellow supergiants (YSG).
Compared to OB stars, LBV, WR and YSG correspond to later stages of evolution.
In particular, WR stars and YSG are likely burning helium in their cores.

The internal structure of OB, LBV, WR and YSG stars changes dramatically, and
this can affect substantially the generation and propagation of internal
gravity waves. The amount of radiative damping is expected to change due to
the large differences in envelope temperature and densities. We recall here
that the radiative damping rate $\gamma_{\rm rad}$ for a traveling g-mode is
given by

$\displaystyle\rm\gamma_{rad}(\omega,\ell,r)=K_{\rm rad}(r)\,k_{r}^{2},$
$\displaystyle K_{\rm
rad}(r)=\frac{16\,\sigma\,T(r)^{3}}{3\,\rho(r)^{2}\kappa(r)\,c_{p}(r)}.$

where $k_{r}$ is the radial wavenumber (eqn. 6) and $\kappa$ is the opacity.
Compact WR stars have surface temperatures that can exceed $\approx 10^{5}$K,
while OB stars and LBVs have effective temperatures $\approx 10^{4\dots
4.7}$K, with LBVs found at the cooler end of this range. The surface of
extended, low-density, YSGs is cooler than $10^{4}$K. With different types of
core convection (H-burning vs He-burning) and radiative damping rates, it
would be surprising if core-generated IGWs in e.g. OB and WR stars showed up
at the surface with similar properties. On the other hand, the FeCZ is present
below/at the surface in both OB and WR stars, and with similar velocities and
convective turnover times. This seems to strengthen the main thesis of this
work, adding support to a (sub)surface origin of the observed SLF variability.

The discussion is more complicated for LBVs and YSGs, since at lower
temperatures other convective regions can become prominent (driven by H and He
recombination, see e.g. Jiang et al., 2018). Dorn-Wallenstein et al. (2020)
disfavor a near-surface origin for the observed SLF variability in YSG, on the
ground that the observed timescale do not follow the predicted scaling for
granulation (Kallinger et al., 2014). However, we believe that the scaling of
Kallinger et al. (2014) is not applicable in the regime explored by Dorn-
Wallenstein et al. (2020). This scaling was derived for a sample of red giants
and solar-like stars, which all have temperatures well-below the temperature
for the recombination of hydrogen, and it is not directly applicable to
earlier spectral subtypes.

Interestingly, the characteristic timescales of 0.1-1 days found by Dorn-
Wallenstein et al. (2020) at $\log{\rm T}_{\rm eff}\approx$4 are consistent
with the convective turnover timescale in the FeCZ (e.g. Fig. 3). Moreover, in
their lower temperature sample ($\log{\rm T}_{\rm eff}<3.75$) the rapid
increase in amplitudes and characteristic frequencies of the variability is in
agreement with the development of near-surface convection induced by the large
opacities associated with the recombination of hydrogen (Grassitelli et al.,
2015a). We suggest that surface and near-surface convection could be indeed
responsible for the variability observed by Dorn-Wallenstein et al. (2020),
and we plan to systematically study the properties of (sub)surface convection
in these evolved, cool stars in future work.

###  6.4 Towards a unified model for surface phenomena in massive stars

The presence of subsurface convection can simultaneously account for a large
variety of puzzling phenomena that appear ubiquitous at the surface of early-
type stars.

  * •

Microturbulence and Macroturbulence can be accounted for by velocity fields
excited by the underlying FeCZ subsurface convection zone (Cantiello et al.,
2009; Grassitelli et al., 2015b; Jiang et al., 2015). The only stars with no
macroturbulence appear to have magnetic fields strong enough to shut off
convection in the FeCZ, while stars with slightly lower surface magnetic
fields show normal values of macroturbulence (Sundqvist et al., 2013).

  * •

Bright spots in early-type stars have been observed (e.g. Ramiaramanantsoa et
al., 2014), and can be explained with the presence of magnetic spots rising
from subsurface convective layers (Cantiello & Braithwaite, 2011, 2019).

  * •

Discrete absorption components (DACs) in UV spectra (e.g. Howarth & Prinja,
1989; Cranmer & Owocki, 1996; Fullerton et al., 1997; Kaper et al., 1997) can
then be caused by the aforementioned bright spots (Cantiello & Braithwaite,
2011) and associated prominences (Sudnik & Henrichs, 2016).

  * •

Line profile variability is another ubiquitous phenomena in hot stars
(Fullerton et al., 1996), and can be explained by surface velocity and density
perturbations seeded by subsurface convection (e.g. Cantiello et al., 2009;
Jiang et al., 2015).

  * •

Wind clumping can also be seeded by these surface density and velocity
perturbations, which are amplified by the development of instabilities in the
stellar wind (Owocki et al., 1988, but see also Sundqvist & Owocki (2013)).

  * •

SLF variability can also be caused by the presence of subsurface convection
zones, as discussed in this work.

An economical hypothesis emerges: the presence of subsurface convection, and
in particular of the FeCZ, could be the common underlying physical cause for
the appearance of turbulence, magnetic spots, SLF variability, as well
spectroscopic variability in early-type stars.

##  7 Conclusions

We used one-dimensional, non-rotating stellar evolution calculations to study
the predicted trends in the properties of subsurface convection in the
spectroscopic H-R diagram. We found that the trends of relative convective
flux and convective turnover timescale in the FeCZ of our models match very
well the trends in timescale and amplitude of stochastic, low-frequency
photometric variability in OB stars observed by TESS and K2.

Similar to previous works, we show that the observed trends in stellar
macroturbulence are also well reproduced assuming the FeCZ is its driver. This
connection is also supported by the observations of strongly magnetized early-
type stars, showing no macroturbulence only for magnetic fields above the
critical value required to shut-off turbulent convection in the FeCZ. We find
that IGWs coming from the stellar core would be reflected or damped for values
of the magnetic field well below this critical value. The fact that stars with
strong but subcritical magnetic fields show typical values of macroturbulence
points against a convective core origin of this surface perturbation. In the
presence of rotation, GIWs are also expected to propagate and reach the
stellar surface. These waves are increasingly confined to stellar equatorial
regions in rapidly-rotating stars, and for stars seen close to pole-on a sharp
decline in their variability is expected below twice their rotational
frequency. This feature is not detected in any of the early-type stars
observed, suggesting GIWs are also unlikely to explain the observed surface
variability.

Overall the observations support a picture in which subsurface convection, and
in particular the FeCZ, is responsible for the ubiquitous low-frequency,
stochastic photometric variability and macroturbulence detected in OB stars.
These surface manifestations join a number of phenomena observed in early-type
stars and attributed to the presence of subsurface convection, including the
observations of (magnetic) bright spots as well as wind and spectroscopic
variabilty. Radiation (magneto)hydrodynamics simulations of the outer envelope
regions of early-type stars are required to understand the details of how
subsurface convection zones cause the observed surface perturbations.

We thank the anonymous referee for a constructive report which helped improve
the manuscript. We also thank Evan Anders for useful discussions on stellar
convection. The Center for Computational Astrophysics at the Flatiron
Institute is supported by the Simons Foundation. This research was supported
in part by the National Science Foundation under Grant No. NSF PHY-1748958 and
by the Gordon and Betty Moore Foundation through Grant GBMF7392. DL is
supported by a Lyman Spitzer Jr. Fellowship.

##  Appendix A Software Details

Calculations were done with MESA version 15140. The MESA EOS is a blend of the
OPAL (Rogers & Nayfonov, 2002), SCVH (Saumon et al., 1995), FreeEOS (Irwin,
2004), HELM (Timmes & Swesty, 2000), and PC (Potekhin & Chabrier, 2010) EOSes.
Radiative opacities are primarily from OPAL (Iglesias & Rogers, 1993, 1996),
with low-temperature data from Ferguson et al. (2005) and the high-
temperature, Compton-scattering dominated regime by Buchler & Yueh (1976).
Electron conduction opacities are from Cassisi et al. (2007). Nuclear reaction
rates are from JINA REACLIB (Cyburt et al., 2010) plus additional tabulated
weak reaction rates Fuller et al. (1985); Oda et al. (1994); Langanke &
Martínez-Pinedo (2000). Screening is included via the prescription of Chugunov
et al. (2007). Thermal neutrino loss rates are from Itoh et al. (1996). The
inlists, processing scripts, and model output are available as a ZENODO
repository: 10.5281/zenodo.4676427 (catalog
https://zenodo.org/record/4676427).

##  Appendix B Radial Displacement from IGWs

Let’s consider traveling internal gravity waves in the WKB approximation. The
wave flux is $F_{\rm IGW}=\frac{1}{2}\rho\,|u|^{2}u_{g,r}$ with the group
velocity $u_{g,r}=\partial\omega/\partial\kappa$. Since
$u_{g,r}=(\omega/N)^{2}\sqrt{N^{2}-\omega^{2}}/k_{h}$,
$|u|^{2}=(N/\omega)^{2}u_{r}^{2}$, and $u_{r}=\omega\xi_{r}$ we obtain:

$F_{\rm
IGW}=\frac{1}{2}\,\rho\,\xi_{r}^{2}\omega^{2}\frac{\sqrt{N^{2}-\omega^{2}}}{k_{h}}.$
(B1)

Using $k_{h}=\sqrt{\ell(\ell+1)}/r$ we can then write the relative radial
displacement as:

$\frac{\xi_{r}}{r}=\left[\frac{F_{\rm
IGW}\sqrt{\ell(\ell+1)}}{1/2\rho\,r^{3}\omega^{2}\sqrt{N^{2}-\omega^{2}}}\right]^{1/2}\propto\sqrt{\frac{F_{\rm
IGW}}{F_{0}}}$ (B2)

##  Appendix C Grids at Different Metallicities

Here we present results for models at metallicities of Z=0.006 and Z=0.002,
representing early-type stars in the Large Magellanic Cloud (LMC, Fig. 7 and
8) and Small Magellanic Cloud (SMC, Fig. 9 and 10) respectively (e.g. Yusof et
al., 2013). The initial helium content of the grids is Y=0.2559 (LMC) and
Y=0.2508 (SMC). The metallicity is initialized scaling the standard solar
composition of Grevesse & Sauval (1998).

(a)  (b)

Figure 7:  Left: Average convective velocities in the FeCZ as function of the
location of stellar models in the spectroscopic H-R Diagram for models with
Z=0.006 (LMC). The FeCZ is absent in models with
$\log\mathscr{L}/\mathscr{L}_{\odot}\lesssim 3.2$. Right: Same as left, but
showing predictions for the maximum of the ratio between turbulent pressure
and total pressure in any subsurface convection zone. The FeCZ largely
dominates, except for stars at low luminosity where turbulent pressure is
provided by a helium convection zone (HeCZ).

(a)  (b)

Figure 8:  Left panel: Ratio of FeCZ convective flux to the total stellar flux
in the spectroscopic H-R Diagram for models with Z=0.006 (LMC). Right panel:
Characteristic frequency $\nu_{\rm char}$ in the FeCZ as a function of the
location of stellar models in the spectroscopic H-R Diagram (black contour
lines). Evolutionary tracks are shown as grey solid lines.

(a)  (b)

Figure 9:  Left: Average convective velocities in the FeCZ as function of the
location of stellar models in the spectroscopic H-R Diagram for models with
Z=0.002 (SMC). The FeCZ is absent in models with
$\log\mathscr{L}/\mathscr{L}_{\odot}\lesssim 3.2$. Right: Same as left, but
showing predictions for the maximum of the ratio between turbulent pressure
and total pressure in any subsurface convection zone. The FeCZ largely
dominates, except for stars at low luminosity where turbulent pressure is
provided by a helium convection zone.

(a)  (b)

Figure 10:  Left panel: Ratio of FeCZ convective flux to the total stellar
flux in the spectroscopic H-R Diagram for models with Z=0.002 (SMC). Right
panel: Characteristic frequency $\nu_{\rm char}$ in the FeCZ as a function of
the location of stellar models in the spectroscopic H-R Diagram (black contour
lines). Evolutionary tracks are shown as grey solid lines.

## References

  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 061102, doi: 10.1103/PhysRevLett.116.061102
  * Aerts (2019) Aerts, C. 2019, arXiv e-prints, arXiv:1912.12300.  https://arxiv.org/abs/1912.12300
  * Aerts et al. (2010) Aerts, C., Christensen-Dalsgaard, J., & Kurtz, D. W. 2010, Asteroseismology
  * Aerts et al. (2019) Aerts, C., Mathis, S., & Rogers, T. M. 2019, ARA&A, 57, 35, doi: 10.1146/annurev-astro-091918-104359
  * Aerts et al. (2009) Aerts, C., Puls, J., Godart, M., & Dupret, M. A. 2009, A&A, 508, 409, doi: 10.1051/0004-6361/200810471
  * Augustson & Mathis (2019) Augustson, K. C., & Mathis, S. 2019, ApJ, 874, 83, doi: 10.3847/1538-4357/ab0b3d
  * Augustson et al. (2020) Augustson, K. C., Mathis, S., & Astoul, A. 2020, ApJ, 903, 90, doi: 10.3847/1538-4357/abba1c
  * Bailey et al. (2015) Bailey, J. E., Nagayama, T., Loisel, G. P., et al. 2015, Nature, 517, 56, doi: 10.1038/nature14048
  * Blomme et al. (2011) Blomme, R., Mahy, L., Catala, C., et al. 2011, A&A, 533, A4, doi: 10.1051/0004-6361/201116949
  * Böhm-Vitense (1958) Böhm-Vitense, E. 1958, ZAp, 46, 108
  * Bowman (2020) Bowman, D. M. 2020, arXiv e-prints, arXiv:2008.11162.  https://arxiv.org/abs/2008.11162
  * Bowman et al. (2020) Bowman, D. M., Burssens, S., Simón-Díaz, S., et al. 2020, A&A, 640, A36, doi: 10.1051/0004-6361/202038224
  * Bowman et al. (2019a) Bowman, D. M., Burssens, S., Pedersen, M. G., et al. 2019a, Nature Astronomy, 3, 760, doi: 10.1038/s41550-019-0768-1
  * Bowman et al. (2019b) Bowman, D. M., Aerts, C., Johnston, C., et al. 2019b, A&A, 621, A135, doi: 10.1051/0004-6361/201833662
  * Buchler & Yueh (1976) Buchler, J. R., & Yueh, W. R. 1976, ApJ, 210, 440, doi: 10.1086/154847
  * Burssens et al. (2020) Burssens, S., Simón-Díaz, S., Bowman, D. M., et al. 2020, A&A, 639, A81, doi: 10.1051/0004-6361/202037700
  * Cantiello & Braithwaite (2011) Cantiello, M., & Braithwaite, J. 2011, A&A, 534, A140, doi: 10.1051/0004-6361/201117512
  * Cantiello & Braithwaite (2019) —. 2019, ApJ, 883, 106, doi: 10.3847/1538-4357/ab3924
  * Cantiello et al. (2009) Cantiello, M., Langer, N., Brott, I., et al. 2009, A&A, 499, 279, doi: 10.1051/0004-6361/200911643
  * Cassisi et al. (2007) Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M. 2007, ApJ, 661, 1094, doi: 10.1086/516819
  * Castro et al. (2015) Castro, N., Fossati, L., Hubrig, S., et al. 2015, A&A, 581, A81, doi: 10.1051/0004-6361/201425354
  * Chugunov et al. (2007) Chugunov, A. I., Dewitt, H. E., & Yakovlev, D. G. 2007, Phys. Rev. D, 76, 025028, doi: 10.1103/PhysRevD.76.025028
  * Cranmer & Owocki (1996) Cranmer, S. R., & Owocki, S. P. 1996, ApJ, 462, 469, doi: 10.1086/177166
  * Cyburt et al. (2010) Cyburt, R. H., Amthor, A. M., Ferguson, R., et al. 2010, ApJS, 189, 240, doi: 10.1088/0067-0049/189/1/240
  * de Mink et al. (2014) de Mink, S. E., Sana, H., Langer, N., Izzard, R. G., & Schneider, F. R. N. 2014, ApJ, 782, 7, doi: 10.1088/0004-637X/782/1/7
  * Dintrans & Rieutord (2000) Dintrans, B., & Rieutord, M. 2000, A&A, 354, 86
  * Dorn-Wallenstein et al. (2020) Dorn-Wallenstein, T. Z., Levesque, E. M., Neugent, K. F., et al. 2020, ApJ, 902, 24, doi: 10.3847/1538-4357/abb318
  * Dufton et al. (2013) Dufton, P. L., Langer, N., Dunstall, P. R., et al. 2013, A&A, 550, A109, doi: 10.1051/0004-6361/201220273
  * Dziembowski (1977) Dziembowski, W. 1977, Acta Astron., 27, 203
  * Edelmann et al. (2019) Edelmann, P. V. F., Ratnasingam, R. P., Pedersen, M. G., et al. 2019, ApJ, 876, 4, doi: 10.3847/1538-4357/ab12df
  * Ferguson et al. (2005) Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585, doi: 10.1086/428642
  * Fuller et al. (1985) Fuller, G. M., Fowler, W. A., & Newman, M. J. 1985, ApJ, 293, 1, doi: 10.1086/163208
  * Fuller et al. (2015) Fuller, J., Cantiello, M., Stello, D., Garcia, R. A., & Bildsten, L. 2015, Science, 350, 423, doi: 10.1126/science.aac6933
  * Fullerton et al. (1996) Fullerton, A. W., Gies, D. R., & Bolton, C. T. 1996, ApJS, 103, 475, doi: 10.1086/192285
  * Fullerton et al. (1997) Fullerton, A. W., Massa, D. L., Prinja, R. K., Owocki, S. P., & Cranmer, S. R. 1997, A&A, 327, 699
  * Godart et al. (2017) Godart, M., Simón-Díaz, S., Herrero, A., et al. 2017, A&A, 597, A23, doi: 10.1051/0004-6361/201628856
  * Goldreich & Kumar (1990) Goldreich, P., & Kumar, P. 1990, ApJ, 363, 694, doi: 10.1086/169376
  * Grassitelli et al. (2015a) Grassitelli, L., Fossati, L., Langer, N., et al. 2015a, A&A, 584, L2, doi: 10.1051/0004-6361/201527289
  * Grassitelli et al. (2016) —. 2016, A&A, 593, A14, doi: 10.1051/0004-6361/201628912
  * Grassitelli et al. (2015b) Grassitelli, L., Fossati, L., Simón-Diáz, S., et al. 2015b, ApJ, 808, L31, doi: 10.1088/2041-8205/808/1/L31
  * Grevesse & Sauval (1998) Grevesse, N., & Sauval, A. J. 1998, Space Sci. Rev., 85, 161, doi: 10.1023/A:1005161325181
  * Hopkins et al. (2014) Hopkins, P. F., Kereš, D., Oñorbe, J., et al. 2014, MNRAS, 445, 581, doi: 10.1093/mnras/stu1738
  * Horst et al. (2020) Horst, L., Edelmann, P. V. F., Andrássy, R., et al. 2020, A&A, 641, A18, doi: 10.1051/0004-6361/202037531
  * Howarth & Prinja (1989) Howarth, I. D., & Prinja, R. K. 1989, ApJS, 69, 527, doi: 10.1086/191321
  * Iglesias & Rogers (1993) Iglesias, C. A., & Rogers, F. J. 1993, ApJ, 412, 752, doi: 10.1086/172958
  * Iglesias & Rogers (1996) —. 1996, ApJ, 464, 943, doi: 10.1086/177381
  * Irwin (2004) Irwin, A. W. 2004, The FreeEOS Code for Calculating the Equation of State for Stellar Interiors.  http://freeeos.sourceforge.net/
  * Itoh et al. (1996) Itoh, N., Hayashi, H., Nishikawa, A., & Kohyama, Y. 1996, ApJS, 102, 411, doi: 10.1086/192264
  * Jermyn & Cantiello (2020) Jermyn, A. S., & Cantiello, M. 2020, ApJ, 900, 113, doi: 10.3847/1538-4357/ab9e70
  * Jiang et al. (2015) Jiang, Y.-F., Cantiello, M., Bildsten, L., Quataert, E., & Blaes, O. 2015, ApJ, 813, 74, doi: 10.1088/0004-637X/813/1/74
  * Jiang et al. (2017) —. 2017, ApJ, 843, 68, doi: 10.3847/1538-4357/aa77b0
  * Jiang et al. (2018) Jiang, Y.-F., Cantiello, M., Bildsten, L., et al. 2018, Nature, 561, 498, doi: 10.1038/s41586-018-0525-0
  * Kallinger et al. (2014) Kallinger, T., De Ridder, J., Hekker, S., et al. 2014, A&A, 570, A41, doi: 10.1051/0004-6361/201424313
  * Kaper et al. (1997) Kaper, L., Henrichs, H. F., Fullerton, A. W., et al. 1997, A&A, 327, 281
  * Kjeldsen & Bedding (1995) Kjeldsen, H., & Bedding, T. R. 1995, A&A, 293, 87.  https://arxiv.org/abs/astro-ph/9403015
  * Krticka & Feldmeier (2021) Krticka, J., & Feldmeier, A. 2021, arXiv e-prints, arXiv:2103.08755.  https://arxiv.org/abs/2103.08755
  * Landstreet et al. (1989) Landstreet, J. D., Barker, P. K., Bohlender, D. A., & Jewison, M. S. 1989, ApJ, 344, 876, doi: 10.1086/167855
  * Langanke & Martínez-Pinedo (2000) Langanke, K., & Martínez-Pinedo, G. 2000, Nuclear Physics A, 673, 481, doi: 10.1016/S0375-9474(00)00131-7
  * Langer (2012) Langer, N. 2012, ARA&A, 50, 107, doi: 10.1146/annurev-astro-081811-125534
  * Langer & Kudritzki (2014) Langer, N., & Kudritzki, R. P. 2014, A&A, 564, A52, doi: 10.1051/0004-6361/201423374
  * Lecoanet et al. (2021) Lecoanet, D., Cantiello, M., Anders, E. H., et al. 2021, arXiv e-prints, arXiv:2105.04558.  https://arxiv.org/abs/2105.04558
  * Lecoanet & Quataert (2013) Lecoanet, D., & Quataert, E. 2013, MNRAS, 430, 2363, doi: 10.1093/mnras/stt055
  * Lecoanet et al. (2019) Lecoanet, D., Cantiello, M., Quataert, E., et al. 2019, ApJ, 886, L15, doi: 10.3847/2041-8213/ab5446
  * MacDonald & Petit (2019) MacDonald, J., & Petit, V. 2019, MNRAS, 487, 3904, doi: 10.1093/mnras/stz1545
  * Maeder et al. (2008) Maeder, A., Georgy, C., & Meynet, G. 2008, A&A, 479, L37, doi: 10.1051/0004-6361:20079007
  * Maeder & Meynet (2000) Maeder, A., & Meynet, G. 2000, ARA&A, 38, 143, doi: 10.1146/annurev.astro.38.1.143
  * Mathis et al. (2014) Mathis, S., Neiner, C., & Tran Minh, N. 2014, A&A, 565, A47, doi: 10.1051/0004-6361/201321830
  * Naze et al. (2021) Naze, Y., Rauw, G., & Gosset, E. 2021, arXiv e-prints, arXiv:2101.05521.  https://arxiv.org/abs/2101.05521
  * Neiner et al. (2020) Neiner, C., Lee, U., Mathis, S., et al. 2020, A&A, 644, A9, doi: 10.1051/0004-6361/201935858
  * Oda et al. (1994) Oda, T., Hino, M., Muto, K., Takahara, M., & Sato, K. 1994, Atomic Data and Nuclear Data Tables, 56, 231, doi: 10.1006/adnd.1994.1007
  * Owocki et al. (1988) Owocki, S. P., Castor, J. I., & Rybicki, G. B. 1988, ApJ, 335, 914, doi: 10.1086/166977
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4, doi: 10.1088/0067-0049/208/1/4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15, doi: 10.1088/0067-0049/220/1/15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34, doi: 10.3847/1538-4365/aaa5a8
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10, doi: 10.3847/1538-4365/ab2241
  * Pedersen et al. (2019) Pedersen, M. G., Chowdhury, S., Johnston, C., et al. 2019, ApJ, 872, L9, doi: 10.3847/2041-8213/ab01e1
  * Penny & Gies (2009) Penny, L. R., & Gies, D. R. 2009, ApJ, 700, 844, doi: 10.1088/0004-637X/700/1/844
  * Potekhin & Chabrier (2010) Potekhin, A. Y., & Chabrier, G. 2010, Contributions to Plasma Physics, 50, 82, doi: 10.1002/ctpp.201010017
  * Prat et al. (2016) Prat, V., Lignières, F., & Ballot, J. 2016, A&A, 587, A110, doi: 10.1051/0004-6361/201527737
  * Press (1981) Press, W. H. 1981, ApJ, 245, 286, doi: 10.1086/158809
  * Ramiaramanantsoa et al. (2014) Ramiaramanantsoa, T., Moffat, A. F. J., Chené, A.-N., et al. 2014, MNRAS, 441, 910, doi: 10.1093/mnras/stu619
  * Ramírez-Agudelo et al. (2013) Ramírez-Agudelo, O. H., Simón-Díaz, S., Sana, H., et al. 2013, A&A, 560, A29, doi: 10.1051/0004-6361/201321986
  * Ratnasingam et al. (2019) Ratnasingam, R. P., Edelmann, P. V. F., & Rogers, T. M. 2019, MNRAS, 482, 5500, doi: 10.1093/mnras/sty3086
  * Ratnasingam et al. (2020) —. 2020, MNRAS, doi: 10.1093/mnras/staa2296
  * Rauw & Nazé (2021) Rauw, G., & Nazé, Y. 2021, MNRAS, 500, 2096, doi: 10.1093/mnras/staa3310
  * Richer et al. (2000) Richer, J., Michaud, G., & Turcotte, S. 2000, ApJ, 529, 338, doi: 10.1086/308274
  * Rogers & Nayfonov (2002) Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064, doi: 10.1086/341894
  * Rogers et al. (2013) Rogers, T. M., Lin, D. N. C., McElwaine, J. N., & Lau, H. H. B. 2013, ApJ, 772, 21, doi: 10.1088/0004-637X/772/1/21
  * Sana et al. (2012) Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444, doi: 10.1126/science.1223344
  * Saumon et al. (1995) Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713, doi: 10.1086/192204
  * Schultz et al. (2020) Schultz, W. C., Bildsten, L., & Jiang, Y.-F. 2020, ApJ, 902, 67, doi: 10.3847/1538-4357/abb405
  * Shiode et al. (2013) Shiode, J. H., Quataert, E., Cantiello, M., & Bildsten, L. 2013, MNRAS, 430, 1736, doi: 10.1093/mnras/sts719
  * Simón-Díaz et al. (2017a) Simón-Díaz, S., Godart, M., Castro, N., et al. 2017a, A&A, 597, A22, doi: 10.1051/0004-6361/201628541
  * Simón-Díaz et al. (2017b) —. 2017b, A&A, 597, A22, doi: 10.1051/0004-6361/201628541
  * Simón-Díaz & Herrero (2014) Simón-Díaz, S., & Herrero, A. 2014, A&A, 562, A135, doi: 10.1051/0004-6361/201322758
  * Smith (2014) Smith, N. 2014, ARA&A, 52, 487, doi: 10.1146/annurev-astro-081913-040025
  * Stevenson (1979) Stevenson, D. J. 1979, Geophysical and Astrophysical Fluid Dynamics, 12, 139, doi: 10.1080/03091927908242681
  * Sudnik & Henrichs (2016) Sudnik, N. P., & Henrichs, H. F. 2016, A&A, 594, A56, doi: 10.1051/0004-6361/201628529
  * Sundqvist & Owocki (2013) Sundqvist, J. O., & Owocki, S. P. 2013, MNRAS, 428, 1837, doi: 10.1093/mnras/sts165
  * Sundqvist et al. (2013) Sundqvist, J. O., Petit, V., Owocki, S. P., et al. 2013, Monthly Notices of the Royal Astronomical Society, 433, 2497, doi: 10.1093/mnras/stt921
  * Sutherland (2010) Sutherland, B. R. 2010, Internal Gravity Waves (Cambridge University Press), doi: 10.1017/CBO9780511780318
  * Szewczuk et al. (2021) Szewczuk, W., Walczak, P., & Daszyńska-Daszkiewicz, J. 2021, MNRAS, 503, 5894, doi: 10.1093/mnras/stab683
  * Timmes & Swesty (2000) Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501, doi: 10.1086/313304
  * Yusof et al. (2013) Yusof, N., Hirschi, R., Meynet, G., et al. 2013, MNRAS, 433, 1114, doi: 10.1093/mnras/stt794