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

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.


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 Ratnasingam et al. 2020) 1 . 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 ) has allowed precise determination of stellar parameters, including spectroscopic mass, luminosity, and macroturbulence ). 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 ≈ 150kK 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).
1 Classical 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).
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 20M 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 nHP.

METHODS
We calculated stellar evolution models using the Modules for Experiments in Stellar Astrophysics (MESA Paxton et al. 2011, 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 120M 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 α 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 α MLT (Jiang et al. 2015;Cantiello & Braithwaite 2019). In Section 6 we discuss how uncertainties in the choice of the α 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 where H P is the pressure scale height calculated at the upper boundary (R 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 M c , convective velocity v c and density ρ 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 ν char and amplitude α 0 of SLF variability. These quantities are derived by fitting the stochastically-excited, lowfrequency power excess in a power density spectrum using a Lorentzian function (e.g. Bowman et al. 2019b): This shows that α 0 represents the amplitude at zero frequency. ν char is defined as in eqn. 3 with τ the typical timescale of the SLF variability. Finally, γ is the gradient of the linear part of the profile and P 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).

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 where τ is a characteristic timescale. For comparing with observations we set τ to be the average convective turnover time calculated either in the FeCZ or in the convective core. We calculated the convective flux F c = ρ c v 3 c , where ρ c and v 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 . 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.

FeCZ and Macroturbulence
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 horizontal 2 We show the average convective velocity in the FeCZ in Fig. 2, left panel. Velocities of the order 10 . . . 100 kms −1 are found across the upper spectroscopic H-R diagram (L ≡ T 4 eff /g; Langer & Kudritzki 2014), with a trend of increasing v 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 ≈ 8 due to a dependence on the uncertain α 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 turb ∝ v 2 c ), which, in agreement with the re- 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 L /L 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).
sults 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 δp ∼ P turb ∝ v 2 c , 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 δp has an associated horizontal velocity u h ∼ (δp/ρ 0 )k h /ω, where k h is the horizontal wavenumber of the wave and ω is its frequency. The dominant waves will have k h ∼ 1/H P and ω ∼ 1/τ c (Cantiello et al. 2009). So running IGWs would have surface velocities which also scale like u h ∼ v 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 highorder 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.

FeCZ and Stochastic, Low-Frequency Variability
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 ∼ 0.1 . . . 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 ν char using eqn. 3 and plot this alongside observed char- Left panel: Characteristic frequency ν 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 ν char , derived from eqn. (2) fitting the data in the range 0.1 ≤ ν ≤ 360 d −1 Bowman et al. 2020). Right panel: Ratio of FeCZ convective flux to the total stellar flux in the spectroscopic H-R Diagram (L ≡ T 4 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 α0. acteristic frequencies 3 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 ν char and the trend with log L and log T eff . While the characteristic frequencies found by Bowman et al. (2020) seem to be larger by a factor of ≈ 3, our predictions for the turnover timescale are affected by uncertainty in the convective velocities arising from the MLT treatment (v c ∝ α 3 MLT , so ν char ∝ α 2 MLT ), as well as our limited knowledge of the frequency spectrum generated by turbulent convection. Since α 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 ≈ 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 F c /F * tend to increase with increasing log L and decreasing log T 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.

Core Convection and Stochastic, Low-Frequency Variability
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 ≈ 10 −2 − 10 2 µmag, seemingly at odds with the relatively large amplitudes ≈ 10 − 10 4 µmag observed by e.g. Bowman  Left panel: Square root of the ratio between the flux of gravity waves launched by the convective core and F0 = 1/2ρ(r) r 3 ω 2 N (r) 2 − ω 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 α0, derived from eqn. (2) fitting the data in the range 0.1 ≤ ν ≤ 360 d −1 Bowman et al. 2020). Right panel: Characteristic frequency ν 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 ν char . The g-mode amplitude should be larger by a factor of ≈ ν/γ, where γ 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., ∼ 10 d −1 for a 10M ZAMS star), but does not change the predicted mode amplitudes for lower frequency waves (e.g., ∼ 0.3 d −1 for a 10M 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 star 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 IGW = F * M c (Goldreich & Kumar 1990), where we evaluate the average Mach number M c using eqn. 1 and the local adiabatic sound speed. We expect the luminosity fluctuations to be proportional to the relative surface radial displacement ξ r /R produced by these waves at the stellar surface (e.g. Dziembowski 1977;Aerts et al. 2010). It can be shown that ξ r /R ∝ F IGW /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 essen-tial 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 ν char = (2πτ c ) −1 , with τ c defined in eqn. 4. Fig. 4 shows that the predicted values of ν char are in the range 0.02 . . . 0.008 d −1 (0.2 . . . 0.08 µ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 T eff and log L . Nevertheless we will assume that trends of ν char , as function of log T eff and log L , are still set by the core convective properties. Under this assumption, Fig. 4 shows that ν char should increase with both log T eff and log 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.

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;Gras-sitelli 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 , closely corresponding to the critical field needed to stabilize the FeCZ (Mac-Donald & 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 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) where ω = 2πν is the angular frequency of the waves, ρ 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 lengthscale 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 by where N is the Brünt-Väisälä frequency. So B r,crit ≈ ω 2 r N 4πρ 2 ( + 1) ≈ ω 2 r N 4πρ.
Note that because this decreases with increasing , all waves of a given frequency are reflected if the = 1 waves are reflected.
We computed B r,crit for several values of and ν for a main-sequence model of a 30M 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 G to 10 −1 G. Strong macroturbulence is observed in similar-mass O-type stars with magnetic fields up to 2.5 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 ν char , on the order of 3 d −1 , the critical magnetic field is much larger, on the order of 300 G to 1 kG. This makes an explanation of macroturbulence based on IGWs marginally inconsistent with observations showing strong macroturbulence up to field strengths of 2.5 kG. However, because eqn. 7 is a strong function of ω, and hence of ν, it is possible that these strongly magnetized stars just have larger ν char . The full range of observed characteristic frequencies spans 0.2−10 d −1 , corresponding to critical field strengths at the surface of 3 G − 10kG for the = 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.  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.

Subsurface Convection
The sample of stars discussed here have projected rotational velocities v eq sin i in the range 7 to 320 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 Ω τ c ), where Ω 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 ≤ 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 eq ≈ 150kms −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 . . . 10. For these values, the convective velocities as calculated in the 1D MLT = 1, ν = ν char, MLT = 2, ν = ν char, MLT = 3, ν = ν char, MLT = 1, ν = ν char, Obs = 2, ν = ν char, Obs = 3, ν = ν char, Obs = 1, ν = 0.2d −1 = 1, ν = 10d −1 Figure 5. The critical radial magnetic field strength needed to reflect internal gravity waves of different and ν is show as a function of fractional depth for a 30M stellar model at an age of 4 Myr. The frequency ν char,MLT is that given by eqn. 3 but evaluated for the core convection zone instead of the subsurface FeCZ. The frequency ν char,Obs is a typical value of 3 d −1 for the observations. We further show frequencies of 0.2 d −1 and 10 d −1 as these span the full range of observed characteristic frequencies.
approximation are affected only at the ∼ 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.

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 lowfrequency variability.
However, we note that it is unlikely that coregenerated GIWs are responsible for the ubiquitous SLF variability. This is because the visibility of GIWs depends on the inclination of the stellar rotation axis re-spect to the observer. Internal waves cannot propagate at the poles for ω < 2Ω, where Ω is the stellar rotational frequency and ω is the wave frequency. The largest wave flux is expected at low latitudes, with the degree of equatorial confinement proportional to Ω. This means that the propagation domain of subinertial (ω < 2Ω) 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Ω. 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.

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(Jiang et al. , 2017(Jiang et al. , 2018Schultz 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 (≈ 50kms −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 pre-cise mechanism connecting subsurface convection zones to the observed surface variability.

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 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 CC = (R /H 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 ≈ √ N CC ≈ 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 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/ √ N 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 latetype 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

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≈ 10 3.2 L for Z=0.02. This corresponds to a zero age main sequence star of about 7M (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≈ 10 3.9 L , corresponding to a zero age main sequence star of about 11M (Cantiello et al. 2009). In our models with metallicity Z=0.006 (See Appendix C) the FeCZ appears at log L /L ≈ 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 ν 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 ∼tens of µ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.

Stochastic, Low-Frequency Variability in Evolved Stars
SLF variability with similar properties to the massive main sequence stars discussed by (Bowman et al. 2019b 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 γ rad for a traveling g-mode is given by γ rad (ω, , r) = K rad (r) k 2 r , K rad (r) = 16 σ T (r) 3 3 ρ(r) 2 κ(r) c p (r) .
where k r is the radial wavenumber (eqn. 6) and κ is the opacity. Compact WR stars have surface temperatures that can exceed ≈ 10 5 K, while OB stars and LBVs have effective temperatures ≈ 10 4...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 coregenerated 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 T eff ≈4 are consistent with the convective turnover timescale in the FeCZ (e.g. Fig. 3). Moreover, in their lower temperature sample (log T 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 nearsurface 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.

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 ).
• 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 ).
• 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.

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 ob-served, 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.

ACKNOWLEDGMENTS
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.

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