Stochastic low-frequency variability of 50 massive stars in the Cygnus OB associations and the Small Magellanic Cloud

ABSTRACT

In recent years, high-precision high-cadence space photometry has revealed that stochastic low-frequency (SLF) variability is common in the light curves of massive stars. We use the data from the Transiting Exoplanet Survey Satellite to study and characterize the SLF variability found in a sample of 49 O- and B-type main-sequence stars across six Cygnus OB associations and one low-metallicity Small Magellanic Cloud star AV 232. We compare these results to 53 previously studied SLF variables. We adopt two different methods for characterizing the signal. In the first, we follow earlier work and fit a Lorentzian-like profile to the power density spectrum of the residual light curve to derive the amplitude |$\alpha _0$|⁠, characteristic frequency |$\nu _{\rm char}$|⁠, and slope |$\gamma$| of the variability. In our second model-independent method, we calculate the root mean square (RMS) of the photometric variability as well as the frequency at 50 per cent of the accumulated power spectral density, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, and the width of the cumulative integrated power density, w. For the full sample of 103 SLF variables, we find that |$\alpha _0$|⁠, |$\gamma$|⁠, RMS, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, and w correlate with the spectroscopic luminosity of the stars. Both |$\alpha _0$| and RMS appear to increase for more evolved stars, whereas |$\nu _{\rm char}$| and |$\nu _{50~{{\ \rm per\ cent}}}$| both decrease. Finally, we compare our results to 2D and 3D simulations of subsurface convection, core-generated internal gravity waves, and surface stellar winds, and find good agreement between the observed |$\nu _{\rm char}$| of our sample and predictions from subsurface convection.

1 INTRODUCTION

Stochastic low-frequency (SLF) variability is ubiquitous in massive stars, from main-sequence O- and B-type stars to the more evolved blue, yellow, and red supergiants (e.g. Balona 1992; Kiss, Szabó & Bedding 2006; Blomme et al. 2011; Tkachenko et al. 2014; Aerts & Rogers 2015; Aerts et al. 2018; Ramiaramanantsoa et al. 2018; Dorn-Wallenstein, Levesque & Davenport 2019; Pedersen et al. 2019; Bowman et al. 2019a, b, 2020, 2024; Burssens et al. 2020; Bowman & Dorn-Wallenstein 2022; Elliott et al. 2022; Ma et al. 2024; Shen et al. 2024; Zhang et al. 2024; Kourniotis et al. 2025) as well as the more exotic Wolf–Rayet (WR) stars (e.g. Lamontagne & Moffat 1987; Lépine & Moffat 1999; Moffat et al. 2008; Chené et al. 2011; Ramiaramanantsoa et al. 2019; Nazé, Rauw & Gosset 2021; Lenoir-Craig et al. 2022). In the frequency domain, such variability is seen as an increase in power towards low frequencies as characteristic of red noise, and has been speculated to be caused by surface granulation (e.g. Kiss et al. 2006), subsurface convection in the Fe peak opacity zone (e.g. Cantiello et al. 2021; Schultz, Bildsten & Jiang 2022, 2023b), internal gravity waves (IGWs) excited by core convection (e.g. Aerts & Rogers 2015; Edelmann et al. 2019; Thompson et al. 2024), or stellar winds (e.g. Krtička & Feldmeier 2018, 2021).

SLF variability has been characterized in different ways. In early studies of SLF variability detected in the space photometry of massive stars, the signal was generally characterized by fitting a Lorentzian-like profile to the amplitude or power density spectrum (e.g. Blomme et al. 2011; Tkachenko et al. 2014; Bowman et al. 2018, 2019a, b, 2020; Nazé et al. 2021; Shen et al. 2024; Kourniotis et al. 2025). Some more recent works have adopted Gaussian process regression to characterize the SLF variability in the time domain (Dorn-Wallenstein et al. 2020; Bowman & Dorn-Wallenstein 2022; Bowman et al. 2024; Zhang et al. 2024). Ma et al. (2024) used a modified Lorentzian-like profile to fit the amplitude spectrum of a sample of stochastically variable blue supergiants in the Large Magellanic Cloud (LMC). Yet other studies suggest that a broken power law may be more appropriate to characterize the SLF variability depending on the origin of the signal (e.g. Krtička & Feldmeier 2021).

The goal of this work is twofold. First, using data from the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2014), we aim to identify and characterize the SLF variability of a sample of O- and B-type stars found in six different Cygnus OB associations containing 655 Galactic OB stars (Quintana & Wright 2021, 2022), and compare these to a smaller well-studied sample of 13 O-type stars in the lower metallicity Small Magellanic Cloud (SMC, Bouret et al. 2021). Secondly, we test a new model-independent method for characterizing the SLF variability which aims to both simplify the SLF characterization process as well as circumvent having to make different choices in which model to use to characterize the variability. We discuss the sample selection in Section 2 as well our new method for characterizing the SLF variability. For comparison, we also choose to fit a Lorentzian-like profile to the power density spectrum and leave a similar comparison to the Gaussian process regression method as future work. We repeat the analysis for a sample of 70 SLF variables previously studied by Bowman et al. (2020) in Section 3. We present our results in Section 4 and compare them to predictions for surface granulation, subsurface convection, IGWs excited by core convection, and stellar winds in Section 5. Our final conclusions are provided in Section 6.

2 METHOD

2.1 Target selection

2.1.1 Cygnus OB association sample

Starting with the initial sample of 655 O- and B-type stars in the six Cygnus OB associations identified by Quintana & Wright (2021, 2022) (group A–F), we include all stars with |$\log L/\mathrm{ L}_\odot \ge 4$|⁠. Amongst these we exclude contaminated stars, as determined from a visual inspection of the TESS pixel data, and stars where the dominant variability is due to coherent mode pulsations (⁠|$\beta$| Cep and SPB stars), eclipsing binaries, or rotational variability. Finally, we add three additional stars with |$\log L/\mathrm{ L}_\odot < 4$| where the dominant variability is clearly SLF. The resulting sample of 54 stars is shown in the Hertzsprung-Russell (HR) diagram in Fig. 1. We further reduce the sample to 49 O- and B-type stars after excluding stars where the signal-to-noise ratio (S/N) for the SLF variability is too low (see Section 2.4). The five stars below the detection threshold are shown as triangles in the figure. These are thus the only non-variable stars in the six Cygnus OB associations with |$\log L/{\rm L}_\odot \ge 4$|⁠, demonstrating that the SLF variability is a common feature among massive stars in line with previous observations (see Introduction). For the rest of this paper, we will refer to this sample as the Cyg OB sample.

2.1.2 SMC sample

To compare our results to findings in low-metallicity environments, we consider also an initial sample of 13 O-type stars in the SMC, whose ultraviolet (UV) and optical spectra have recently been analysed in detail by Bouret et al. (2021). Stars observed by TESS in the SMC have a higher risk of contamination due to the higher crowding combined with the large |$21{{\rm arcsec} } \times 21{{\rm arcsec} }$|  TESS pixel sizes. The tess_localizepython package (Higgins & Bell 2023) has recently been used successfully to localize the source of variability in TESS light curves of evolved massive stars in the LMC and SMC (Pedersen & Bell 2023). However, this tool takes as input a list of frequencies of coherent signals to determine the location where the amplitude of these signals is highest and therefore by design cannot be used to localize the source of SLF variability. Instead, we extract, background correct, and normalize the light curves of the nearby pixels surrounding the target following the method outlined in Section 2.2 and study these light curves in combination with their corresponding power density spectra (PDS) to determine if the detected variability originates from the target or a nearby star.

Based on this, we find one of the stars (AV 232) to be an SLF variable with an additional three (AV 43, AV 83 and AV 327) candidate SLF variables. Inspecting the pixel data of AV 83 shows that the star does appear to be an SLF variable, however, due to the high crowding the SLF signal is very likely diluted which would impact its characterization. AV 327 has an adjacent similar brightness star that is always in the same pixel as the target and it is therefore unclear if the signal is coming from AV 327 or its neighbour. AV 43 shows a small amount of SLF variability but is contaminated by a nearby Cepheid variable in some of its sectors. We therefore choose to label these three stars as SLF candidates and only characterize the signal of AV 232. Of the remaining 10 stars, two are found to be eclipsing binaries (AV 25 and AV 75), one is contaminated by a much brighter nearby star and another by a nearby Cepheid variable. The remaining two stars only contain white noise with (some) occasional contamination from nearby Cepheid variables.

2.2 The TESS data

The members of the Cygnus OB associations were observed by TESS in its observing cycle 2, 4, a few in cycle 5, and in cycle 6. Out of the 49 stars discussed in this work, 20 have 2-min cadence available, whereas the rest only have full-frame image (FFI) data available. Over these four TESS cycles, the cadence of the FFI data changed from 30-min, to 10-min and 200 s, corresponding to a change in Nyquist frequency from 277 to 833 and |$2499\, \mu$|Hz. An overview of the sample and the available TESS data is given in Table A1.

The SMC was observed by TESS in cycles 1, 3, and 5, with the FFI data of cycles 1 and 3 sharing the same observing cadence as cycles 2 and 4, respectively. AV 232 was observed in sectors 1, 28, 67, and 68, and only has 2-min cadence data available for the last two sectors. An overview of the SMC O-type sample falling in the category of SLF variable, SLF candidate, and not variable (i.e. white noise only) is provided in Table A2.

The light curves were extracted from the TESS FFI data and 2-min cadence target pixel files using the python packages lightkurve (Lightkurve Collaboration 2018) and tesscut (Brasseur et al. 2019) using custom target pixel masks. To correct for background signals and systematics, a background pixel mask was created from which a time-series of the background signals were extracted from each background pixel and a principal component analysis subsequently carried out to reduce the background signals in the time domain to their first one to seven principal components depending on the target. These principal components were used to remove the background signals from the target light curve using an adapted version of the RegressionCorrector functionality of lightkurve. The target light curves were then normalized by dividing by a low (typically zeroth or first) order polynomial fit to the first and second half of each sector (corresponding to two different TESS orbits) and the flux units changed to parts-per-million (ppm). Finally, the time stamps of the TESS FFI data were corrected following the methodology adapted in the TESS Asteroseismic Science Operations Centre pipeline developed by TESS Data for Asteroseismology group under the TESS Asteroseismic Science Consortium (Handberg et al. 2021). For more information on the light-curve extraction and preparation, please see Pedersen et al. (in preparation).

2.3 Iterative prewhitening

Before a study of the red noise, i.e. the SLF variability, can be carried out we first remove any remaining periodic coherent signals from the light curves, which are unlikely to be related to the dominant SLF variability. We use an iterative prewhitening scheme, where sinusoidal signals are removed one at a time until the S/N ratio of an extracted signal goes below a predefined S/N limit. The optimal choice of this limit depends on the probability that the extracted signal may be a random noise peak instead of a genuine intrinsic signal of the extracted light curve. For one sector of TESS data, Pedersen et al. (in preparation) have shown that an S/N limit of 3.4, 4.0, and 4.6, where the noise is calculated using a |$11.57\, \mu {\rm Hz}$| window around the extracted signal, corresponds to a probability of 50 per cent, 75 per cent, and 90 per cent that a signal of the given S/N value is not a random noise peak. In this work, we adopt |$\mathrm{ S/N} = 4$| as the limit to be consistent with previous studies of SLF variability (Bowman et al. 2018, 2019a, b, 2020). To carry out the iterative prewhitening, we extract signals one by one from the observed light curve by calculating its corresponding Lomb–Scargle periodogram (Lomb 1976; Scargle 1982) to identify the frequency |$\nu _{\rm peak}$| and amplitude |$A_{\rm peak}$| of the highest S/N peak at each step of the iteration. Using the identified |$\nu _{\rm peak}$| and |$A_{\rm peak}$| as initial guesses, we fit the equation

to the light curve using the non-linear least-squares fitting algorithm provided by the lmfit python package (Newville et al. 2024), allowing |$\nu _i$|⁠, |$A_i$|⁠, and |$\phi _i$| of all the extracted |$i=1,\dots , N$| signals to be optimized at each step of the iteration. Here, |$A_i$| is the amplitude of signal i, |$\nu _i$| is the frequency, |$\phi _i$| is the phase, and C is a constant. The derived model in equation (1) is subsequently subtracted from the original light curve and a periodogram calculated for the resulting residual light curve to determine the next |$\nu _{\rm peak}$| and |$A_{\rm peak}$|⁠. This iterative process continues until the S/N of the next signal drops below four.

Following Pápics et al. (2017), we carry out the iterative prewhitening in the order of the highest to lowest significant signals rather than highest to lowest amplitude as for data with a significant amount of red noise the latter carries the risk of reaching the stopping criterion before all significant signals have been extracted, as lower amplitude higher- S/N signals at higher frequencies may exist. By prewhitening according to peak significance, we effectively remove twice as many coherent signals compared to prewhitening according to the amplitudes.

2.4 Fitting the power density spectrum

Following earlier studies of SLF variability in massive stars (Blomme et al. 2011; Tkachenko et al. 2014; Bowman et al. 2018, 2019a, b, 2020), we model the red noise in the power density spectrum of the residual light curve using the Lorentzian-like function

where |$\nu$| is the frequency, |$\alpha _0$| is the PDS at |$0\, \mu {\rm Hz}$|⁠, |$\gamma$| defines the slope of the variability, |$\nu _{\rm char}$| is the characteristic frequency, and |$C_w$| is the white noise contribution to the PDS. The spectrum below |$\nu _0 = 1.157\, \mu {\rm Hz}$| is influenced by how the light curve is detrended, and therefore we do not model the spectrum for frequencies below this value. To convert the power to power density we multiply the power by the effective length of the TESS sector, calculated as the inverse of the area under the spectral window (Kjeldsen et al. 2005).

The first term of equation (2) has been multiplied by the power attenuation (Chaplin et al. 2011; Huber et al. 2022)

to account for the time averaging of high-frequency signals caused by long exposure times of the earlier TESS FFI data (see Fig. 2). If the true power of a signal is |$P_{\rm true}$|⁠, then the observed power is modulated by the power attenuation and corresponds to |$P_{\rm obs} = \eta P_{\rm true}$|⁠. As attenuation will have the highest impact on the results obtained from the 30-min cadence data, we choose to exclude data from TESS cycles 1 and 2 for the remained of this paper. See Appendix  C for further discussion of cadence dependence of the parameter estimates.

To fit equation (2) to the PDS we derive posterior probability distributions and Bayesian evidence with the nested sampling Monte Carlo algorithm MLFriends (Buchner 2016, 2019) using the UltraNest python package (Buchner 2021). We adopt 400 live points to sample the parameter space assuming that |$\alpha _0$|⁠, |$\nu _{\rm char}$|⁠, |$\gamma$|⁠, and |$C_W$| all have uniform priors. At each iteration a model is calculated and evaluated using a log-likelihood function. Under the assumption of Gaussian noise in the time domain, which translates to a |$\chi ^2$| distribution in power with two degrees of freedom in the frequency domain, the corresponding log-likelihood function can be shown to be (Duvall & Harvey 1986; Toutain & Appourchaux 1994)

and has commonly been used when fitting Lorentzian functions (or Harvey profiles; Harvey 1985) to the granulation background in solar-like stars as well as their oscillation frequencies (e.g. Davies et al. 2016; Lund et al. 2017; Li et al. 2020; Nielsen et al. 2021). We adopt the same log-likelihood function in equation (4) to evaluate equation (2) at each iteration.

The final parameter estimates are taken as the 50th percentiles of the samples and the 16th and 84th percentiles their corresponding lower and upper values. For the SLF variability detection to be considered to be significant, we require that |$\alpha _0/C_w > 10$| for more than half of the sectors that a given star was observed in. Five stars in our initial Cyg OB sample of 54 O- and B-type stars failed this criterion, leaving us with the final sample of 49 stars reported in Table A1. An example residual FFI light curve (top panel) and its corresponding fitted PDS (bottom left) is shown in Fig. 3.

2.5 RMS and integrated power

As an alternative to using the Lorentzian-like model in equation (2) to characterize the SLF variability, we also calculate the root-mean-squared (RMS) variability of the residual light curve as a separate measure of the amplitude of the SLF variability, as well as the frequencies at which 20 per cent, 50 per cent, and 80 per cent of the power is located. This approach has the benefit that it does not rely on any model fits, but instead the quantities are calculated directly from the data. The calculated RMS is shown against the TESS magnitude of the 49 Cyg OB stars in Fig. 4 (circles), with the TESS noise models from Sullivan et al. (2015) and Schofield et al. (2019) shown for comparison. The results for the five excluded stars from this sample are shown as open triangles, while the SMC star AV 232 is indicated by a star.

We calculate the cumulative integrated power of the power density spectrum |$S (\nu)$| as

where |$\nu \in [\nu _0, \nu _{\rm norm}]$|⁠. Once again we use |$\nu _0 = 1.157\, \mu {\rm Hz}$| as the lower frequency limit. The upper frequency |$\nu _{\rm norm}$| defines the frequency at which the cumulatively integrated power density spectrum is normalized to one. Just as in Section 2.4 we exclude data from TESS cycle 1 + 2 from this analysis and choose to use the Nyquist frequency of the 10-min cadence data as |$\nu _{\rm norm}$|⁠.

We use the integrated power to calculate the frequencies below which 20 per cent, 50 per cent, and 80 per cent of the power resides (⁠|$\nu _{20~{{\ \rm per\ cent}}}$|⁠, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, and|$\nu _{80~{{\ \rm per\ cent}}}$|⁠), corresponding to |$P_{\rm int} (\nu _{20~{{\ \rm per\ cent}}}) = 0.2$|⁠, |$P_{\rm int} (\nu _{50~{{\ \rm per\ cent}}}) = 0.5$|⁠, and |$P_{\rm int} (\nu _{80~{{\ \rm per\ cent}}}) = 0.8$|⁠. Additionally, we also calculate the ‘width’

of the spectrum, which we expect to depend on |$\gamma$| and |$\nu _{\rm char}$|⁠. The bottom-right panel of Fig. 3 provides an example of the normalized cumulative integrated power density for the 10-min cadence FFI data of the star Gaia EDR3 2059070135632404992 from sector 41. The vertical lines show the value of |$\nu _{20~{{\ \rm per\ cent}}}$| (orange), |$\nu _{50~{{\ \rm per\ cent}}}$| (blue), and |$\nu _{80~{{\ \rm per\ cent}}}$| (green) for this light curve.

As previously mentioned, we choose the same |$\nu _{\rm norm}$| in the derivation of |$\nu _{20~{{\ \rm per\ cent}}}$|⁠, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, and |$\nu _{80~{{\ \rm per\ cent}}}$| for all sectors and observing cadences. This is done for the sake of consistency, allowing for direct comparisons between the variation of these parameters from sector to sector due to the intrinsic variability of the star. This enforces the use of the Nyquist frequency of the longest cadence FFI data as the highest meaningful value of |$\nu _{\rm norm}$|⁠. To demonstrate how the |$\nu _{20~{{\ \rm per\ cent}}}$|⁠, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, |$\nu _{80~{{\ \rm per\ cent}}}$|⁠, and w parameters would change if a different choice of |$\nu _{\rm norm}$| was made, we used the Lorentzian-like model in equation (2) to simulate the SLF variability in PDS for a 2-min cadence data set. The white noise was fixed to |$\log C_W = 0.2$|⁠, while |$\log \alpha _0 \in [1,10]$|⁠, |$\nu _{\rm char} \in [1,100]$|⁠, and |$\gamma \in [1,5]$| were varied within the given ranges. The four parameters |$\nu _{20~{{\ \rm per\ cent}}}$|⁠, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, |$\nu _{80~{{\ \rm per\ cent}}}$|⁠, and w were then calculated for four different values of |$\nu _{\rm norm}$|⁠, corresponding to the Nyquist frequency of a 30-min, 10-min, 200-s, and 2-min cadence light curve, respectively.

We find that the differences between the estimated parameters for the four different cadences are minimal provided that the amplitude of the SLF variability (⁠|$\log \alpha _0$|⁠) is sufficiently high with little dependence on the observing cadence. The differences are less than 1 |$\mu$|Hz for |$\nu _{20~{{\ \rm per\ cent}}}$|⁠, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, and |$\nu _{80~{{\ \rm per\ cent}}}$| if |$\alpha _0$| is larger than |$10^2$|⁠, |$10^3$|⁠, and |$10^4$| ppm|$^2\mu$|Hz|$^{-1}$|⁠, respectively, and less than 0.1 for w if |$\alpha _0 > 10^4$| ppm|$^2\mu$|Hz|$^{-1}$|⁠. In comparison, when considering the dependence of the four parameters |$\nu _{20~{{\ \rm per\ cent}}}$|⁠, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, |$\nu _{80~{{\ \rm per\ cent}}}$|⁠, and w on |$\gamma$| and |$\nu _{\rm char}$| a much stronger cadence dependence is found, especially for |$\nu _{20~{{\ \rm per\ cent}}}$|⁠, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, and |$\nu _{80~{{\ \rm per\ cent}}}$| with the differences being largest for |$\nu _{80~{{\ \rm per\ cent}}}$|⁠. These differences also show a clear correlation between |$\gamma$| and |$\nu _{\rm char}$| and become larger when both |$\gamma$| and |$\nu _{\rm char}$| increase in value. The difference between the estimated |$\nu _{20~{{\ \rm per\ cent}}}$|⁠, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, and |$\nu _{80~{{\ \rm per\ cent}}}$| values of the 30- and 2-min cadence data are larger than 5 |$\mu$|Hz (20 |$\mu$|Hz) for |$\gamma$|-values larger than 2.5 (1.5), 3 (2.2), and 3.5 (2.8), respectively, for all values of |$\nu _{\rm char}$|⁠. When comparing the results for the 10- and 2-min cadence simulated data, the differences between the estimated |$\nu _{20~{{\ \rm per\ cent}}}$|⁠, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, and |$\nu _{80~{{\ \rm per\ cent}}}$| values are less than 5 |$\mu$|Hz for |$\gamma \gtrsim 1.5$|⁠, 2, and 2.5, respectively.

Based on the results of these comparisons, we choose to exclude the FFI data from TESS cycles 1 and 2 from our analysis including both the calculations of |$\nu _{20~{{\ \rm per\ cent}}}$|⁠, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, |$\nu _{80~{{\ \rm per\ cent}}}$|⁠, and w as well as the Lorentzian-like model fit to the PDS described in Section 2.4 and adopt |$\nu _{\rm norm} = 833\, \mu$|Hz, corresponding to the Nyquist frequency of the 10-min cadence data.

3 COMPARISON STARS

To compare our results to prior studies of SLF variability we consider the sample of 70 O- and B-type stars studied by Bowman et al. (2020) who also used a Lorentzian-like model to characterize the SLF variability using TESS 2-min cadence data. One important difference is that Bowman et al. (2020) did the fitting in the amplitude spectrum instead of in PDS, meaning that we cannot draw a direct comparison between their resulting parameter estimates and the ones presented here.

To circumvent this issue, we redo the fitting for the Bowman et al. (2020) sample in the PDS and also derive RMS, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, and w for the sample. To ensure that we are comparing similar data to those used by Bowman et al. (2020), we consider only the 2-min cadence TESS Pre-search Data Conditioning Simple Aperture Photometry (PDCSAP) light curves for this sample. These light curves were downloaded from the Mikulski Archive for Space Telescopes (MAST) using lightkurve. We normalized the light curves by dividing by a first-order polynomial fit to each sector and change the flux units to ppm. To otherwise keep the analysis method similar to the approach adopted in this work, we estimate the |$\alpha _0$|⁠, |$\nu _{\rm char}$|⁠, |$\gamma$|⁠, |$C_W$|⁠, RMS, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, and w parameters separately for each sector. Finally, we exclude stars with a low S/N detection of the SLF variability as well as stars where the dominant variability is due to binarity, rotation, or coherent pulsations. This reduces the sample of comparison stars to 53 O- and B-type stars, see Table B1 in Appendix  B for an overview. For the remainder of this paper, we will refer to this sample as the B20 sample.

4 RESULTS

In Table A3 in Appendix  A, we provide the average values across all sectors and observing cadences of the estimated parameters |$\alpha _0$|⁠, |$\nu _{\rm char}$|⁠, |$\gamma$|⁠, |$C_W$|⁠, RMS, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, and w for the Cyg OB sample and the SMC star AV 232. The corresponding parameter estimates for the B20 sample are listed in Table B2 in Appendix  B. Following previous studies of SLF variability, we now investigate the interparameter dependencies as well as the dependence on stellar evolution stage by placing the stars in an HR diagram. Because the B20 sample only has spectroscopic luminosities (⁠|$\mathcal {L} = T_{\rm eff}^4/g$|⁠) available, we first convert our bolometric luminosities L to spectroscopic luminosities using

derived using |$L = 4\pi \sigma R^2 T_{\rm eff}^4$| and |$g = GM/R^2$|⁠. The derived spectroscopic luminosities are likewise listed in Table A3. One low-luminosity star in our Cyg OB sample is a clear outlier from the rest of our sample, and so we choose to exclude it from the rest of this discussion but nevertheless list its SLF parameter estimates in Table A3.

Fig. 5 shows the six estimated parameters as a function of |$\mathcal {L}$| for the Cyg OB sample (blue circles), SMC star AV 232 (pink star), and the B20 sample (orange diamonds). The grey bars on |$\alpha _0$|⁠, |$\nu _{\rm char}$|⁠, |$\gamma$|⁠, |$C_W$|⁠, RMS, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, and w represent the full retrieved range in the estimated parameter across all sectors and observing cadences for a given star, while the symbols correspond to the average value. Large bars are therefore indicative of large variations in these parameters over the different sectors, which means that the SLF variability changes characteristics over time.

The trends in the individual panels of Fig. 5 indicate that there is a correlation between most of these parameters and the spectroscopic luminosity, independently of if all three stellar samples are considered simultaneously or independently. To quantify this, we calculate the Spearman’s rank correlation coefficients |$r_s$| (Spearman 1904), which are applicable and can be calculated for two variables that share a monotonic relationship. The coefficient takes values between |$-1$| and 1, where |$r_s=1$| and −1 correspond to a perfect positive and negative correlation between the two variables, respectively, and |$r_s=0$| means that the parameters are uncorrelated. How close |$\left| r_s \right|$| is to one determines the strength of the correlation: |$\left| r_s \right| \in ]0, 0.19]$| is a very weak correlation, |$\left| r_s \right| \in [0.2, 0.39]$| is weak, |$\left| r_s \right| \in [0.4, 0.59]$| is moderate, |$\left| r_s \right| \in [0.6, 0.79]$| is strong, and |$\left| r_s \right| \in [0.8, 1.0]$| is a very strong correlation. Finally, each calculated |$r_s$| is associated with a p-value which determines the significance of the correlation. For |$p > 0.05$|⁠, the correlation is considered insignificant, while |$0.001 < p \le 0.05$| corresponds to significant, and |$p < 0.001$| is highly significant.

The first row of Table 1 lists the Spearman’s ranks correlation coefficients between |$\log \mathcal {L}/\mathcal {L}_\odot$| and the six estimated parameters used to characterize the SLF variability for the full sample of stars considered in this paper (Cyg OB  + AV 232  + B20). The colour indicates the p-value for the considered correlation. As seen in the table, significant weak to strong positive correlations are found for |$\alpha _0$|⁠, |$\gamma$|⁠, and RMS, while significant weak to moderate negative correlations are found for |$\nu _{50~{{\ \rm per\ cent}}}$| and w. For |$\nu _{\rm char}$|⁠, the correlation is insignificant and we cannot reject the null hypothesis that |$r_s = 0$|⁠. Similar generally stronger correlations are found if the Cyg OB sample is considered on its own, whereas the B20 sample also shows no correlation between |$\log \mathcal {L}/\mathcal {L}_\odot$| and the two parameters |$\gamma$| and |$\nu _{50~{{\ \rm per\ cent}}}$|⁠.

The full sample of stars including the Cyg OB, AV 232, and B20 samples are shown in the spectroscopic HR diagrams in Fig. 6, where the colours of the symbols in each subplot correspond to the derived average parameter estimate for each of the six parameters |$\alpha _0$|⁠, |$\nu _{\rm char}$|⁠, |$\gamma$|⁠, |$C_W$|⁠, RMS, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, and w as indicated by the colour bars. The shape of the symbols have the same meaning as in Fig. 5. Based on this, figure we can see that stars in the same part of the HR diagram generally have, (i) high RMS values when |$\alpha _0$| is high, (ii) high |$\nu _{50~{{\ \rm per\ cent}}}$| values when |$\nu _{\rm char}$| is low, and (iii) high w values when |$\gamma$| is low, and vice versa. The trends with evolutionary stage are less clear especially when the full sample is considered, but are indicative of the |$\alpha _0$| and RMS both increase when the stars get older, while |$\nu _{\rm char}$| and |$\nu _{50~{{\ \rm per\ cent}}}$| both decrease. No clear correlation with age is seen for either |$\gamma$| or w. Comparing our results in Fig. 6 to the study of the trends in macroturbulence velocities of |$\approx 600$| O- and B-type stars from Serebriakova, Tkachenko & Aerts (2024) shows excellent agreement, where stars with high |$\alpha _0$| and RMS are situation in the same parts of the HR diagram as stars with high macroturbulence velocities. This is consistent with earlier findings of relations between photometric SLF variability and spectroscopic macroturbulence by Bowman et al. (2020).

The correlation coefficients between |$\log T_{\rm eff}$| and each of the six estimated SLF parameters are shown in the second rows of each of the table segments in Table 1. When the full sample is considered, we find that neither |$\alpha _0$|⁠, RMS, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, or w are significantly correlated with |$\log T_{\rm eff}$|⁠, but if we consider the Cyg OB and B20 samples separately a different pattern emerges for these four parameters. When a parameter is significantly positively correlated with |$\log T_{\rm eff}$| for Cyg OB it is negatively correlated with |$\log T_{\rm eff}$| for B20, and vice versa. Because of this, when the two samples are combined no correlation with |$\log T_{\rm eff}$| is found. For |$\nu _{\rm char}$|⁠, a very strong correlation with |$\log T_{\rm eff}$| is found for B20, while it is insignificant for Cyg OB, resulting in an overall significant moderate correlation for the full sample. The opposite is seen for the RMS parameter. The |$\gamma$| parameter is the only one that is simultaneously found to be both significant and positively correlated with |$\log T_{\rm eff}$| for both the Cyg OB and B20 sample, but the correlation is stronger for the B20 sample. For the Cyg OB sample, we find a positive correlation between |$\log T_{\rm eff}$| and both |$\alpha _0$| and RMS, which is consistent with recent findings for a sample of blue supergiants of lower effective temperatures than our sample (Kourniotis et al. 2025).

Fig. 7 plots the six estimated parameters that we use to characterize the SLF variability. The colour of the symbols indicate the spectroscopic luminosity. In combination with the correlation coefficients in Table 1, we find that |$\alpha _0$| is correlated with all other parameters except |$\gamma$| for the B20 sample. The |$\nu _{\rm char}$| parameter is correlated with |$\alpha _0$| and |$\gamma$| for Cyg OB, whereas it is also correlated with |$\log T_{\rm eff}$|⁠, |$\alpha _0$|⁠, RMS, |$\nu _{50~{{\ \rm per\ cent}}}$| and w for the B20 sample. The |$\gamma$| parameter is correlated with all other parameters for the Cyg OB sample, but not with |$\log \mathcal {L}$|⁠, |$\alpha _0$|⁠, RMS or w for B20. The RMS is correlated with all other parameters except for |$\nu _{\rm char}$| for Cyg OB, and correlates with |$\log \mathcal {L}$|⁠, |$\alpha _0$|⁠, |$\nu _{\rm char}$|⁠, |$\nu _{50~{{\ \rm per\ cent}}}$| and w for the B20 sample. The |$\nu _{50~{{\ \rm per\ cent}}}$| parameter is correlated with all parameters except for |$\log T_{\rm eff}$| and |$\nu _{\rm char}$| for Cyg OB and |$\log \mathcal {L}$| for B20. Finally, the width w is correlated with all other parameters except for |$\nu _{\rm char}$| for the Cyg OB sample and for |$\log T_{\rm eff}$| and |$\gamma$| for the B20 sample. When all stars are combined, the only insignificant correlations are found for |$\log T_{\rm eff}$| versus |$\alpha _0$|⁠, |$\log T_{\rm eff}$| versus RMS, |$\log \mathcal {L}$| versus |$\nu _{\rm char}$|⁠, |$\alpha _0$| versus |$\gamma$|⁠, |$\log T_{\rm eff}$| versus |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, |$\log T_{\rm eff}$| versus w, and |$\gamma$| versus |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, i.e. for seven parameter combinations out of 21. Six of the significant correlations are found to be both strong and highly significant.

For the lower metallicity SMC star AV 232 we find that its SLF parameters are generally at the edges of the parameter estimates for both the Cyg OB and B20 sample, but are similar to the values and correlations found for stars with similar spectroscopic luminosities.

5 DISCUSSION

5.1 Surface granulation

For solar-like oscillators, the background signal arising from a combination of surface granulation noise and stellar activity is usually modelled and removed by fitting a sum of power laws (i.e. Harvey models) with 2–5 components to the PDS and including a constant for the white noise component, resulting in

for the background PDS (⁠|${\rm PDS}_{\rm Bkg}$|⁠), which simplifies to equation (2) when only one power-law component is included in the background model. Unlike in this work and previous studies of SLF variability in massive stars, the |$c_i$| parameter is held fixed when fitting equation (8) to the background noise of solar-like oscillators. When initially introducing this model, Harvey (1985) adopted |$c_i = 2$| for the Sun, while later it was shown that |$c_i = 4$| appears to be a more appropriate value (e.g. Aigrain, Favata & Gilmore 2004; Michel et al. 2009). As seen in Fig. 5, most of our estimated |$\gamma$| fall between these two |$c_i$| values.

In a study of the connection between granulation noise and oscillation signal in solar-like, subgiant and red giant stars observed by the Kepler Space Telescope, Kallinger et al. (2014) considered eight different models for the granulation background including a two-component Harvey model of the form of equation (8). Their derived correlation between the two characteristic granulation frequencies |$b_1$| and |$b_2$| and the frequency at maximum power |$\nu _{\rm max}$| of the solar-like oscillations are shown as dashed and full grey lines in Fig. 8, respectively, and compared to our estimates for the Cyg OB, SMC star, and B20 sample. Here, we have adopted

following Kallinger & Matthews (2010) based on Kjeldsen & Bedding (1995) and Huber et al. (2009). The derived |$M R^{-2} T_{\rm eff}^{-1/2}$| values are given with respect to solar values such that |$M R^{-2} T_{\rm eff}^{-1/2} = 1$| for the Sun. If a single component is used to fit the background, as done for the stars in this work, the granulation frequency |$\nu _{\rm gran}$| should fall between the two grey lines.

As seen in Fig. 8, an offset is seen between the OB stars and the derived |$\nu _{\rm char}$| versus |$M R^{-2} T_{\rm eff}^{-1/2}$| relations found by Kallinger et al. (2014) for the lower mass solar-like oscillators, with the massive stars showing consistently lower |$\nu _{\rm char}$| values. This is in line with previous findings by Bowman et al. (2019b), who argued that the observed SLF variability in their sample of (near-) main-sequence OBA stars is unlikely to be caused by surface granulation for this reason. If granulation is the cause, then the properties of the resulting red noise cannot be directly extrapolated from the relations found for cooler stars (see also similar discussion by Blomme et al. 2011). For comparison, we also show the characteristic frequencies derived by Kallinger & Matthews (2010) for two |$\delta$| Sct stars (red triangles), which are oscillating main-sequence A-type stars expected to have very shallow convective envelopes. Contrary to the OB stars, their |$\nu _{\rm char}$| values do align with the relation found for the cooler stars, suggesting that their observed red noise could be caused by surface granulation.

For the more evolved massive red supergiants (RSGs), the situation may be different. These stars are known to have large surface convective cells with convective turnover timescales on the order of hundreds of days (e.g. Goldberg, Jiang & Bildsten 2022), and have also been shown to exhibit SLF variability (Kiss et al. 2006; Zhang et al. 2024). Using ground-based photometric light curves of 48 RSGs from the database of the American Association of Variable Star Observers with a mean time span of 61 yr, Kiss et al. (2006) showed that the variability followed a |$1/\nu$| frequency dependence consistent with expectations for solar granulation background, but could not derive a |$\nu _{\rm char}$| due to the much sparser sampling of the data compared to current photometric space missions. On the other hand, Zhang et al. (2024) find granulation time-scales on the same order as the predictions from Goldberg et al. (2022) by modelling the SLF variability seen in ground-based Optical Gravitational Lensing Experiment (OGLE) and All-Sky Automated Survey for Supernovae (ASAS-SN) light curves of more than 6000 RSGs in the LMC and SMC using Gaussian processes.

5.2 Sub-surface convection

Stothers & Chin (1993) were some of the first to show that the opacity bumps near the stellar surface of hot massive stars caused by the partial ionization of iron and helium (Iglesias, Rogers & Wilson 1992) give rise to dynamical instabilities, also known as subsurface convection zones. Cantiello et al. (2009) showed that while the He opacity bump in the majority of cases does not form a significant subsurface convection zone, convection is always developed around the iron opacity bump for |$L \gtrsim 10^{3.2}\, {\rm L}_\odot$| provided that the metallicity of the star is sufficiently high. In general, this subsurface convection zone from the iron opacity bump (FeCZ) increases in importance and becomes more prominent for lower surface temperatures and gravities and higher luminosities and metallicities. Later studies by Jermyn, Anders & Cantiello (2022) demonstrated that the subsurface convection zones are present over a narrower range in stellar properties than previously estimated and is expected to be absent below |$\approx 16$| and |$\approx 35\, {\rm M}_\odot$| for stars in the LMC and SMC, respectively, while the FeCZ should be present for all masses above |$\approx 8\, {\rm M}_\odot$| for stars of solar metallicity. Bowman et al. (2024) compared the predictions by Jermyn et al. (2022) to the observed positions of SLF variables in the HR diagram, finding SMC stars located within regions predicted to be stable against subsurface convection in the FeCZ at low metallicities.

The estimated surface velocities arising from gravity waves excited by the FeCZ appear to correlate with the occurrence of significant microturbulence in O- and B-type stars, indicating a possible physical connection between the two (Cantiello et al. 2009). Additionally, it has been shown that SLF photometric variability and macroturbulence could be caused by these subsurface convection zones (Cantiello et al. 2021). With the exception of three stars, the SLF variables in our Cyg OB sample all fall in the regime where substantial subsurface FeCZs are expected to be found. This is largely due to our sample selection, in which we focused on stars with |$L/{\rm L}_\odot \ge 4$| and deliberately avoided stars where the photometric variability is dominated by rotational variability, pulsations, or binarity. Of the stars above this threshold, 96  per cent show SLF variability. The SMC star AV 232, however, is found right at the boundary (⁠|$M_{\rm spec} = 35.3 \pm 8.2\, {\rm M}_\odot$|⁠, Bouret et al. 2021) where prominent subsurface FeCZ are expected to be developed according to Jermyn et al. (2022) and shows similar SLF amplitudes as the Galactic OB Cyg sample at similar spectroscopic luminosities.

The theoretical considerations of sub-surface convection zones mentioned above all relied on the use of 1D stellar models. In recent years, detailed 3D numerical simulations of these subsurface convection zones have been carried out, partially in an attempt to explain the SLF variability. Schultz et al. (2022) performed 3D radiation hydrodynamical simulations of the outer |$\approx 15~{{\ \rm per\ cent}}$| of the stellar envelopes of two |$35\, {\rm M}_\odot$| stars at two different evolutionary stages, one on the zero-age main sequence (model T42L5.2) and the other half-way through the main-sequence evolution (model T35L5.0). They find that the amplitudes of the tangential and radial velocities from the simulations are comparable to the predictions from 1D models, however, the velocity profiles from the 3D simulations are much broader and extend all the way to the surface. While the simulated convection from the Fe opacity bump carries minimal flux, the turbulent motions reach the surface and there is therefore no convectively quiet zone in the outer |$\approx 7~{{\ \rm per\ cent}}$| of the simulated stars. Through an analysis of the synthetic light curves produced by the two models, Schultz et al. (2022) find similar slopes as comparison SLF variables from Bowman et al. (2020) found around the same position in the HR diagram, while their derived |$\nu _{\rm char}$| are consistent with the thermal time-scale in the FeCZ. The same authors later added a |$13\, {\rm M}_\odot$| terminal-age main-sequence 3D envelope model simulation to their sample of simulated stars (M13TAMS, Schultz et al. 2023b) and likewise determine the |$\nu _{\rm char}$| of the model. All of these simulations were run at solar metallicity.

The results for all three simulated stars are included in Fig. 8. A dotted line connects the two |$35\, {\rm M}_\odot$| models from Schultz et al. (2022, open black cross and plus symbols) and are indicative of the direction the stars should be moving over time if the SLF variability is caused by subsurface convection. As seen in the figure, this trend aligns with the parameter range of the majority of the stars in the Cyg OB sample with the more evolved T35L5.0 model even crossing the observed trend for surface granulations in solar-like oscillators (dashed grey). Similarly, the results for the simulated |$13\, {\rm M}_\odot$| star matches the observations of a part of the Cyg OB sample.

The same general trends are seen in Fig. 9 where |$\nu _{\rm char}$| is shown as a function of |$\log \mathcal {L}/\mathcal {L}_\odot$|⁠. As seen in the figure, the observed scatter in |$\nu _{\rm char}$| at similar |$\log \mathcal {L}/\mathcal {L}_\odot$| values could be the result of a change in characteristics of subsurface convection zone during the main-sequence evolution. To confirm this, more simulations at similar masses but different evolutionary stages would be required. Additionally, the weakness of the subsurface convection theory lies in the disappearance of such prominent zones at lower metallicities as predicted from 1D models (Cantiello et al. 2009; Jermyn et al. 2022), while SLF variability has previously been identified in stars located in the LMC (e.g. Pedersen et al. 2019; Bowman et al. 2019a). With our |$M_{\rm spec} = 35.3 \pm 8.2\, {\rm M}_\odot$| SMC star being located right at the boundary where substantial FeCZ should be able to develop at SMC metallicity, having a 3D hydrodynamical simulation of such a star would be highly valuable, however, to our knowledge no such simulations currently exist. For an extensive review of simulations of massive star envelopes, we refer the reader to Jiang (2023).

5.3 Internal gravity waves excited by core convection

At interfaces between convective and radiative zones, IGWs are stochastically exited by processes that causes disturbances at the convective boundary such as the overshooting/penetration of convective plumes (Townsend 1966; Montalbán & Schatzman 2000) and Reynolds stresses of turbulent convective eddies (Lighthill 1952; Lecoanet & Quataert 2013). Such excited IGWs propagate through the radiative zones and away from the convective regions either towards higher (low-mass stars) or lower (high-mass stars) density regions. As discussed in Section 5.2, such waves are expected to be generated both above and below subsurface convection zones. For stars with convective cores, IGWs excited by core convection have been demonstrated to be capable of efficient angular momentum transport (Rogers et al. 2013), efficient mixing of chemical elements (Rogers & McElwaine 2017; Varghese et al. 2023), and potentially cause enhanced mass loss in massive stars towards the end of their evolution (Quataert & Shiode 2012). For an extensive review on simulations of convective cores, please see Lecoanet & Edelmann (2023).

While the earliest simulations of convective cores focused on the study of convective boundary mixing (e.g. Deupree 2000; Browning, Brun & Toomre 2004), increased attention has been put on the study of the resulting IGW spectrum and its possible connection to SLF variability. Aerts & Rogers (2015) took the IGW spectrum from the 2D hydrodynamical simulations of the inner 0.5–98  per cent in radius of a |$3\, {\rm M}_\odot$| star (Rogers et al. 2013) and rescaled both the frequency and amplitude of the spectrum to match it to the photometric observations of three young O-type stars showing SLF variability, finding good general agreement between the two except at the very low-frequency end. A similar approach was taken by Ramiaramanantsoa et al. (2018) to study the red noise component of the O-type star V973 Scorpii. They showed that doing the rescaling of the simulated IGW spectrum assuming a |$45\, {\rm M}_\odot$| TAMS star resulted in a better match to the observations than for a similar evolutionary stage |$34\, {\rm M}_\odot$| star.

Edelmann et al. (2019) carried out the first 3D hydrodynamical simulations covering 1 per cent–90  per cent in radius of a rotating |$3\, {\rm M}_\odot$| zero-age main-sequence star allowing them to study the properties of the IGWs generated by the convective core throughout the majority of the radiative envelope. They find the generated spectrum to be well represented by a broken power law, where the frequency position of the break point depends on the angular degree |$\ell$| of the waves. In Figs 8 and 9, we show the range in frequency of these break points found for the Fourier transform of the kinetic energy spectrum of |$\ell \in [1,20]$| as a vertical full black line. The overlapping cross is the break-point frequency estimated for the combined Fourier spectrum of the corresponding temperature fluctuations over all |$\ell$| values. As seen in Fig. 8, the range in break-point frequencies matches well with the range in |$\nu _{\rm char}$| for the B20 sample at the given |$M R^{-2} T_{\rm eff}^{-1/2}$| value. In Fig. 9, the results from the 3D simulations do not overlap with the observations as the stars studied in this work are of much higher masses and therefore higher |$\log \mathcal {L}/\mathcal {L}_\odot$|⁠. This is different for the 3D hydrodynamical simulations of a |$25\, {\rm M}_\odot$| main-sequence star carried out by Herwig et al. (2023) whose simulations cover the inner 54 per cent of the star. Based on these simulations, Thompson et al. (2024) derived the |$\nu _{\rm char}$| resulting from fitting equation (2) to the resulting IGW spectrum. Their result is indicated in Fig. 9 by the black square and is situated at the lower end in the parameter ranges for the observations at similar |$\log \mathcal {L}/\mathcal {L}_\odot$|⁠.

The simulations by Rogers et al. (2013), Herwig et al. (2023), and Edelmann et al. (2019) all required a boosting of the convective core to compensate for the high thermal diffusivity in the simulations. In the case of the 2D simulations by Rogers et al. (2013), Aerts & Rogers (2015) argued that such boosting could mean that the generated IGW spectrum might be more representative of a |$\approx 30\, {\rm M}_\odot$| star. Assuming a similar argument can be made for the 3D simulations by Edelmann et al. (2019), this would move the full vertical black line in Fig. 9 to the position of the black dashed line for a |$34\, {\rm M}_\odot$| star.¹ In this case, the range in break-point frequency approximately spans the observed range in |$\nu _{\rm char}$| at the given |$\log \mathcal {L}/\mathcal {L}_\odot$|⁠, potentially implying that the scatter is caused by a variation in which |$\ell$| values are more dominant.

The required boosting of the convective core is one of the criticisms of the IGW theory as the source of the observed SLF variability as it does not allow for a prediction of the amplitude of the waves. Recently, Anders et al. (2023) carried out 3D simulations of three zero-age main-sequence stars at 3, 15, and |$40\, {\rm M}_\odot$| covering the inner 93  per cent of the star in radius and without boosting the convective core in the simulations. They find the amplitudes of the generated IGW spectrum to be a few orders of magnitudes lower than the observed red noise in O- and B-type stars. The derived |$\nu _{\rm char}$| are likewise lower than the observations that they compare their results to (Bowman et al. 2020) and decrease with mass. Here, we have likewise included the results from Anders et al. (2023) in Figs 8 and 9. As seen in Fig. 8, the |$\nu _{\rm char}$| values from the simulations fall along the lower boundary of the observed values. Contrary to the findings by Anders et al. (2023, see their fig. 3), the |$\nu _{\rm char}$| values from the simulations do fall in regions in Fig. 9 that are covered by the observations, however, the majority of the stars in our combined sample have higher |$\nu _{\rm char}$| values at similar |$\log \mathcal {L}/\mathcal {L}_\odot$|⁠. Furthermore, as seen in both figures the break-point frequencies found for the |$3\, {\rm M}_\odot$| simulation with core luminosity boosting are much higher than the |$\nu _{\rm char}$| value from corresponding |$3\, {\rm M}_\odot$| simulation from Anders et al. (2023). Another point of contingency for the IGW theory is the question if it is possible for gravity waves, that are damped in convective regions, to reach the surface of stars with subsurface convection zones. Based on 1D stellar models, Serebriakova et al. (2024) demonstrated that IGWs can indeed tunnel through the FeCZ if the star is less than |$30\, {\rm M}_\odot$| and thereby contribute to the observed variability. Full tests of this are currently awaiting 2D and 3D simulations as none (to our knowledge) of the available published simulations of core generated IGWs include the subsurface convection zones.

5.4 Stellar winds

Massive stars (⁠|$M \gtrsim 15\, {\rm M}_\odot$|⁠, Abbott 1982) have radiatively driven stellar winds which show direct observational features in the stellar spectra when their luminosities are higher than |$10^4\, {\rm L}_\odot$| (Abbott 1979; Kudritzki & Puls 2000). When the mass-loss rates caused by such stellar winds are variable as due to, e.g. line-drive wind instabilities (Feldmeier, Puls & Pauldrach 1997; Runacres & Owocki 2002) they result in photometric variability due to the effect of wind blanketing (Abbott & Hummer 1985).

Relying on simulations of stellar winds of O-type stars, Krtička & Feldmeier (2018) demonstrated that the resulting photometric variability from stellar winds is stochastic and shows similar characteristics to the SLF variability found in four comparison blue supergiants, with amplitudes decreasing for more structured stellar winds. Krtička & Feldmeier (2021) expanded upon this study by considering the impact of the choice of inner boundary atmospheric perturbations on the stellar wind. For stochastic perturbations expected by, e.g. subsurface convection zones, Krtička & Feldmeier (2021) find that the resulting frequency spectrum follows a uniform power law close to the base of the perturbations but develops into a broken power law at larger heights above the boundary. Additionally, they find that the disc-integrated photometric variability in mmag shows negative skewness calculated as |$\overline{\left(x(t)-\overline{x(t)}\right)^3}/\sigma _{\rm std}^3$|⁠, where |$\sigma _{\rm std}$| is the standard deviation of the photometric variability |$x(t)$| in mmag as a function of time t.

The studies of the simulated SLF variability caused by stellar winds carried out by Krtička & Feldmeier (2018, 2021) do not provide enough information to place the predictions in either Fig. 8 or Fig. 9. Instead we calculate the skewness of the residual light curves of the Cyg OB, B20, and SMC samples after converting the flux to mmag and plot the skewness as a function of spectroscopic luminosity in Fig. 10. While 66  per cent of the stars show a negative skewness when averaged across all TESS sectors, we find no clear dependence on |$\log \mathcal {L}$|⁠. The lower metallicity SMC star shows a positive skewness. As Krtička & Feldmeier (2018, 2021) predict the signal of the SLF variability caused by stellar winds to be larger for higher mass-loss rates |$\dot{M}$|⁠, this could potentially be consistent with the lower |$\dot{M}$| found for lower metallicity stars (Mokiem et al. 2007), however, the sample size is too small to say anything conclusive. In comparison, Kourniotis et al. (2025) also calculated the skewness of their sample of 41 Galactic blue supergiants, finding values in the range of −0.69 to 0.34 with |$\approx 54~{{\ \rm per\ cent}}$| of the stars having negative skewness values.

The hydrodynamical simulations of subsurface convection and convective core generated IGWs discussed above do not provide any estimates for the skewness of the predicted SLF variability. It is therefore unclear if negative skewness of SLF light curves can be uniquely linked to stellar winds. However, given how different choices of boundary perturbations impact the predicted photometric variability (Krtička & Feldmeier 2021) combined with the results from 3D hydrodynamical simulations showing that the velocities generated by the iron subsurface convection zone extend all the way to the surface (Schultz et al. 2022) and that the wind structure starts to form already at the iron opacity bump (Debnath et al. 2024), it is highly likely that these processes are closely related.

6 CONCLUSIONS

In this work, we used two different methods to characterize the SLF variability found in a sample of 49 O- and B-type stars located within six Cygnus OB associations and one low-metallicity star in the SMC, AV 232. An additional three SMC stars out of an initial sample of 13 are found to be candidate SLF variables. For comparison, we also included 53 additional O- and B-type stars with SLF variability which was previously studied by Bowman et al. (2020). For the first method we modelled the SLF variability in the power density spectrum using a Lorentzian-like profile, resulting in estimates of four parameters: |$\alpha _0$|⁠, |$\nu _{\rm char}$|⁠, |$\gamma$|⁠, and |$C_W$|⁠. The second method uses the photometric time-series combined with the cumulative integrated power density spectrum to characterize the SLF variability using three parameters: RMS, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, and w.

We find that both |$\nu _{50~{{\ \rm per\ cent}}}$| and w are much less impacted by the observing cadence of the TESS data than |$\nu _{\rm char}$| and |$\gamma$| and recommend using the RMS, |$\nu _{50~{{\ \rm per\ cent}}}$| and w parameters in future characterization of SLF variability for this reason. Similarly, Bowman & Dorn-Wallenstein (2022) found the use of Gaussian process regression to be less sensitive to the data quality than fitting models to the variability in the Fourier domain. Upon studying the sector-to-sector dependence of the derived SLF parameters, we find that the SLF variability changes characteristics over time. Additionally, we find that all derived parameters (except |$\nu _{\rm char}$|⁠) are significantly correlated with the spectroscopic luminosity of the stars for the combined sample of 103 studied O- and B-type stars, whereas these correlations are no longer significant for |$\gamma$| and |$\nu _{50~{{\ \rm per\ cent}}}$| if the B20 sample is considered on its own. Studying the values of RMS, |$\alpha _0$|⁠, |$\nu _{\rm char}$|⁠, and |$\nu _{50~{{\ \rm per\ cent}}}$| relative to the position of the stars in the spectroscopic HRD indicate that both |$\alpha _0$| and RMS increase for more evolved stars, whereas |$\nu _{\rm char}$| and |$\nu _{50~{{\ \rm per\ cent}}}$| decrease. No similar clear correlation is found for either |$\gamma$| or w.

In line with previous studies (Bowman et al. 2018, 2019b), we find that the SLF variability is unlikely to be caused by surface granulation. When comparing the position of the full sample of stars in the |$\nu _{\rm char}$| versus |$\log \mathcal {L}/{\mathcal {L}_\odot }$| diagram to predictions from hydrodynamical simulations of subsurface convection and IGWs excited by core convection, we find good agreement between the observations and simulations of subsurface convection zones (Schultz et al. 2022, 2023b). Similar agreement can be found for the simulations of IGWs generated by convective cores in some parts of the diagram.

Having more 2D and 3D simulations available covering different masses, ages, and metallicities would be highly valuable to compare against observations and help determine the dominating process causing the SLF variability. Particularly, also publishing and making publicly available the simulated light curves is especially useful for anyone attempting to make comparisons to observations where different methods for characterizing the SLF variability have been applied, as is the case for this work. Here, we were limited to making comparisons to |$\nu _{\rm char}$|⁠. For the same reasons, we could not make similar direct comparisons to simulations of photometric variability caused by variable stellar winds and had to derive a new parameter, the skewness of the photometric variability, to compare to predictions from one simulation. With more than half of the 103 stars studied here showing negative skewness, this might support a connection to stellar winds.

The low-metallicity SMC star AV 232 shows similar SLF characteristics to the 102 studied Galactic OB stars, and is located at the boundary where Jermyn et al. (2022) predict that a substantial subsurface FeCZ should be able to develop according to 1D models. Having additional 3D simulations of the envelope of a low-metallicity SMC star near this boundary would therefore be very valuable. Bowman et al. (2024) studied a larger sample of SMC, LMC, and Galactic OB stars, demonstrating that the lower metallicity SLF variables are found both within and outside the same predicted FeCZ stability windows. Again, having 3D simulations of stars in the same position of the HRD but at different metallicities would provide an interesting comparison for the predicted SLF signals.

Finally, in this work, we focused on characterizing the SLF variability present in the TESS light curves of our sample of stars. Previous studies have also linked such derived SLF parameters with spectroscopic observations of macroturbulence (Bowman et al. 2020; Shen et al. 2024), which may likewise be caused by core-generated IGWs (Aerts & Rogers 2015) and subsurface convection (Cantiello et al. 2009, 2021; Schultz et al. 2023a) with the dominating process likely depending on the position of the star in the HRD (e.g. Serebriakova et al. 2024). Krtička & Feldmeier (2018, 2021) suggested that studying the variability in the UV could also be used to distinguish between different effects causing SLF variability (Krtička 2016). Most recently, the SLF variability observed in TESS light curves of the O4 supergiant |$\zeta$| Puppis has also been linked to variations in the linear polarization of the light (Bailey et al. 2024). While carrying out similar studies are beyond the scope of this work, we look forward to following the continued utilization of synergies between different observations and predictions from simulations to study the origin of the SLF variability in massive stars.

ACKNOWLEDGEMENTS

The authors are grateful for discussions with Tim Bedding, Dennis Stello, and Courtney Crawford at earlier stages of the manuscript, and appreciate the comments by Conny Aerts, Dominic Bowman, and Pieterjan Van Daele to the initial submitted version of the manuscript. We thank the anonymous referee for their comments which helped improve the manuscript. This research was supported in part by the NASA ATP grant ATP80NSSC22K0725, and by the National Science Foundation under grant no. NSF PHY-2309135 at the Kavli Institute for Theoretical Physics (KITP), as well as through the TESS Guest Investigator program cycle 4 under grant no. 80NSSC22K0743 from NASA, and by the Professor Harry Messel Research Fellowship in Physics Endowment, at the University of Sydney. This paper includes data collected with the TESS mission, obtained from the MAST data archive at the Space Telescope Science Institute (STScI). Funding for the TESS mission is provided by the NASA Explorer Program. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.

We acknowledge the use of the following python packages: astropy (Astropy Collaboration 2013, 2018, 2022), astroquery (Ginsburg et al. 2019), lightkurve (Lightkurve Collaboration 2018), tesscut (Brasseur et al. 2019), numpy (Harris et al. 2020), matplotlib (Hunter 2007), pandas (Wes McKinney 2010; pandas development team 2020), lmfit (Newville et al. 2024), scipy (Virtanen et al. 2020), MLFriends (Buchner 2016, 2019), and UltraNest (Buchner 2021).

DATA AVAILABILITY

The original TESS FFI data, target pixel files, and PDCSAP light curves for the studied SLF variables are publicly available on the Michulski Archive for Space Telescopes (MAST). Residual light curves for each star, data type (FFI or 2-min cadence), and corresponding PDS are publicly available on Zenodo (DOI:10.5281/zenodo.15261328). Similarly, Tables A1–B2 are publicly available on Zenodo in machine readable format, as are all necessary data for reconstructing figures (Figs 8–10) when not already included in one of the tables.

Footnotes

REFERENCES

APPENDIX A: OVERVIEW OF AND FITTING RESULTS FOR THE CYG OB AND SMC SAMPLE

An overview of the sample fo 49 SLF variables found across the six Cygnus OB associations including their name (Gaia EDR3 and TIC ID), association group, TESS (T) and Gaia (G) magnitudes, effective temperatures, luminosities and masses, as well as the available TESS sectors with FFI and 2-min cadence data is provided in Table A1. A similar overview of the six SMC stars from Bouret et al. (2021) that are found to either an SLF variable, SLF candidate or not variable is provided in Table A2. The corresponding derived SLF parameters |$\alpha _0$|⁠, |$\nu _{\rm char}$|⁠, |$\gamma$|⁠, |$C_W$|⁠, RMS, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, and w are provided in Table A3 along with the derived spectroscopic luminosity |$\log \mathcal {L}/\mathcal {L}_\odot$|⁠.

APPENDIX B: OVERVIEW OF AND FITTING RESULTS FOR THE B20

Similar to Appendix  A, an overview of the sample of 53 comparison OB stars from Bowman et al. (2020) known to be SLF variables is provided in Table B1 and their corresponding derived SLF parameters are given in Table B2.

APPENDIX C: CADENCE DEPENDENCE OF ESTIMATED PARAMETERS

With a few exceptions (Nazé et al. 2021; Lenoir-Craig et al. 2022), the analysis of SLF variability based on TESS observations has generally relied on the use of the 2-min cadence data. There are a couple of reasons for this. One is that the TESS 2-min cadence light curves are easily accessible as they are automatically made available on MAST by the TESS Science Processing Operations Center (SPOC; Jenkins et al. 2016), whose TESS science pipeline produces both simple aperture photometry and PDCSAP light curves. The same is not the case for light curves constructed from FFI data, which therefore require additional steps to be prepared and ready for analysis such as those described in Section 2.2. The TESS SPOC team has recently started running their pipeline also on FFI data (Caldwell et al. 2020), while other authors have also used their own pipelines to produce TESS FFI light curves that are made available publicly on MAST (e.g. Bouma et al. 2019; Nardiello et al. 2019; Handberg et al. 2021; Hon et al. 2021), none of which are complete.

Another reason for using the 2-min cadence TESS data instead of FFI data is related to amplitude attenuation (e.g. Chaplin et al. 2014) also sometimes called amplitude suppression (e.g. Bowman et al. 2020), where the time averaging of high-frequency signals due to increased exposure times results in an attenuation of their amplitudes. In the case of TESS data, the exposure time equals the observing cadence and the fractional attenuation in power can be written as equation (3).

As most of the stars our OB Cyg sample only have FFI data available, we quantify the impact of the sampling rate on the estimated SLF parameters by taking all of the 2-min cadence residual light curves for both the OB Cyg and B20 sample and resampling them to match an observing cadence of 240 s, 10 min, and 30 min. We use 240 s instead of 200 s as it corresponds to a multiple integer of the 2-min cadence sampling. The resampling is done by binning the normalized residual flux to the new cadence and dividing by the number of data points in each bin to get the averaged flux, while the time stamps are changed to the center values of each bin. Once the light curves have been resampled, we derive the |$\alpha _0$|⁠, |$\nu _{\rm char}$|⁠, |$\gamma$|⁠, |$C_W$|⁠, |$RMS$|⁠, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, and w parameters again and compare them to the estimated values of the 2-min cadence data.

Fig. C1 shows the differences in the estimated |$\log \alpha _0$| (top), |$\nu _{\rm char}$| (middle), and |$\gamma$| (bottom) parameters obtained when the 2-min cadence data are resampled to 240-s (green), 10-min (orange), and 30-min (blue) cadence. As an example, the y-axis on the top panel corresponds to |$\Delta \log \alpha _0 = \log \alpha _{0,{\rm resampled}}-\log \alpha _{0,{\rm 2min}}$|⁠, where |$\log \alpha _{0,{\rm resampled}}$| is the estimated value for, e.g. the resampled 240-s cadence data and the |$\log \alpha _{0,{\rm 2min}}$| is the corresponding estimate for the 2-min cadence data. Instead of showing the individual parameter estimates for all sectors and all stars in the Cyg OB and B20 sample with 2-min cadence data available, we calculate the averages in 20 equally sized bins (full connected lines) as well as the corresponding |$1\sigma _{\rm std}$| indicated by the shaded green (240-s), hatched orange (10-min), and dotted blue (30-min) regions in Fig. C1. For all three parameters, the estimates approach the 2-min cadence values as the cadence increases as expected. In the 10- and 30-min cadence cases, the high-amplitude SLF signals are underestimated at longer cadences. In the case of |$\nu _{\rm char}$|⁠, both the 240-s and 10-min cadence results are close to the 2-min cadence estimates, while a clear trend is seen for the 30-min cadence data with the differences increasing towards higher characteristic frequencies. The bottom panel of Fig. C1 shows that while the curves representing the average differences are generally flat the standard deviations of the differences become smaller when the cadence is increased.

Fig. C2 is the same as Fig. C1, but for the three parameters RMS, |$\nu _{\rm 50~{{\ \rm per\ cent}}}$| and w. Because |$\nu _{\rm norm}$| is set to the Nyquist frequency of the 10-min cadence data, we excluded the 30-min cadence comparison for the two bottom panels of Fig. C2. We limited the comparisons to stars and sectors with |$\nu _{\rm 50~{{\ \rm per\ cent}}} < 50\, \mu {\rm Hz}$| and |$w < 12$| as only few data points are available above these values. An opposite trend with observing cadence is seen from RMS compared to |$\alpha _0$| where the differences in the RMS estimates decrease for increasing RMS, contrary to the increasing differences found for increasing |$\alpha _0$| values in Fig. C1.

A summary of the impacts of changing the observing cadence on the estimated parameters is provided in Table C1. The values in the table are the total 1|$\sigma$| standard deviations of the differences between the estimated parameters for the resampled data and the original 2-min cadence residual light curves, and not just the binned standard deviations shown in Figs C1 and C2. As in Fig. C2, we excluded stars and sectors with |$\nu _{\rm 50~{{\ \rm per\ cent}}} \ge 50\, \mu {\rm Hz}$| and |$w \ge 12$| in the derivation of the total standard deviations of these two parameters and provide the corresponding values for the full sample in the parentheses. A comparison between Figs C1 and C2 as well as of the |$1\sigma _{\rm std}$| values in Table C1 shows that the RMS, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, and w parameters are generally less sensitive to the observing cadence than |$\alpha _0$|⁠, |$\nu _{\rm char}$|⁠, and |$\gamma$| at least for |$\nu _{\rm 50~{{\ \rm per\ cent}}} < 50\, \mu {\rm Hz}$| and |$w < 12$|⁠.

In their study of SLF variability in Galactic WR stars, Lenoir-Craig et al. (2022) also considered the impact of longer observing cadences on the estimated |$\alpha _0$|⁠, |$\nu _{\rm char}$|⁠, |$\gamma$|⁠, and |$C_W$| parameters derived from fitting equation (2) to the amplitude spectrum rather than the PDS. To do so, they looked at the averages and standard deviations of the ‘normalized parameters’, e.g. |$\alpha _0 {\rm [new\ cadence]}/\alpha _0 {\rm [2min]}$|⁠, obtained from resampling the 2-min cadence TESS data to 10, 30, and 100 min. For the sake of comparison, we repeat this exercise here and show in Fig. C3 the change in normalized parameters for the resampled 240-s, 10-min, and 30-min cadence data for the combined Cyg OB and B20 sample. The filled blue points show the averages and associated standard deviation obtained for each cadence; these values are repeated in Table C2. The orange triangles in Fig. C3 and values in parenthesis in Table C2 show the same results after limiting the sample to stars and sectors with |$\nu _{\rm 50~{{\ \rm per\ cent}}} < 50\, \mu {\rm Hz}$| and |$w < 12$|⁠. As seen in the figure, the effects of limiting the sample are negligible for |$\log \alpha _0$|⁠, |$\nu _{\rm char}$|⁠, and |$\gamma$|⁠, whereas the standard deviations significantly decrease for the normalized |$\log {\rm RMS}$|⁠, |$\nu _{50~{{\ \rm per\ cent}}}$|⁠, and w parameters, while their averages generally move closer to one, especially for |$\log {\rm RMS}$|⁠. The open blue circles and open brown triangles show the corresponding result obtained if |$\eta (\nu)$| is excluded from the first term in equation (2) used to fit the PDS, demonstrating that while including |$\eta (\nu)$| improves the results it does not fully solve the issues with the longer cadences. For this reason, we chose to exclude the TESS cycles 1 and 2 data from our analysis to minimize the impact of possible cadence dependencies on our results.

For their sample of WR stars, Lenoir-Craig et al. (2022) found the amplitudes and slopes to be most affected by the change in observing cadence, with |$\alpha _0$| (⁠|$\gamma$|⁠) showing a negative (positive) trend and larger standard deviations for longer observing cadences (cf. their fig. 5). The same albeit smaller trends are observed for our sample of OB stars. However, in contrast to our results they find |$\nu _{\rm char}$| to be least affected, while in our case it is the normalized parameter that changes the most with itsaverage value becoming almost twice as large for the 30-min cadence sampling. We attribute this in part to be due to the WR stars having a smaller range in observed |$\nu _{\rm char} \in [1.2, 36.6]\ \mu$|Hz. As shown in Fig. C1, it is above |$\approx 30\ \mu$|Hz that the largest deviations in |$\nu _{\rm char}$| between the observing cadences occur, with our considered full sample of OB stars reaching |$\nu _{\rm char}$| values more than twice as high as the WR stars studied by Lenoir-Craig et al. (2022).

Finally, we find again that both |$\nu _{50~{{\ \rm per\ cent}}}$| and w are much less affected by the change in observing cadence than |$\nu _{\rm char}$| and |$\gamma$| also when the normalized parameters are considered.