Variability of Blue Supergiants in the LMC with TESS

The blue supergiant (BSG) problem, namely, the overabundance of BSGs inconsistent with classical stellar evolution theory, remains an open question in stellar astrophysics. Several theoretical explanations have been proposed, which may be tested by their predictions for the characteristic time variability. In this work, we analyze the light curves of a sample of 20 BSGs obtained from the Transiting Exoplanet Survey Satellite (TESS) mission. We report a characteristic signal in the low-frequency (f ≲ 2 day−1) range for all our targets. The amplitude spectrum has a peak frequency of ∼0.2 day−1, and we are able to fit it by a modified Lorentzian profile. The signal itself shows strong stochasticity across different TESS sectors, suggesting its driving mechanism happens on short (≲months) timescales. Our signals resemble those obtained for a limited sample of hotter OB stars and yellow supergiants, suggesting their possible common origins. We discuss three possible physical explanations: stellar winds launched by rotation, convection motions that reach the stellar surface, and waves from the deep stellar interior. The peak frequency of the signal favors processes related to the convective zone caused by the iron opacity peak, and the shape of the spectra might be explained by the propagation of high-order, damped gravity waves excited from that zone. We discuss the uncertainties and limitations of all these scenarios.


INTRODUCTION
The formation pathways of blue supergiants (BSGs) are still not well understood.Despite being a unique class in the Hertzsprung-Russell (H-R) diagram, their observed population disagrees with predictions from the classical theory of massive star evolution (Figure 1).Namely, if BSGs are massive single stars that have triggered hydrogen shell burning after central hydrogen depletion, their lifetime should be characterized by a thermal timescale on which their stellar envelope expands, typically 100 − 1000 times faster than their main-sequence lifetime (Hoyle 1960;Hofmeister et al. 1964).This means the occurrence rate of BSGs should Corresponding author: Linhao Ma lma3@caltech.edube 100 − 1000 times lower than massive main-sequence stars, clearly not consistent with their observed population (Castro et al. 2014(Castro et al. , 2018;;de Burgos et al. 2023).This discrepancy is known as the "blue supergiant problem," and it remains one of the biggest open questions in the theory of stellar evolution (Bellinger et al. 2023).
Several alternative models for single stellar evolution have been proposed to solve the blue supergiant problem.Some authors argue that BSGs are actually centralhelium-burning stars with smaller cores, and typically longer lifetimes.However, such theories usually require little or no mixing near the convective core boundary (Walborn et al. 1987;Weiss 1989;Schootemeijer et al. 2019), a scenario disfavored by our current-day knowledge of massive stars (Kaiser et al. 2020;Martinet et al. 2021).Some other authors argue that BSGs could be core-hydrogen-burning stars that extend their main se-Ma et al.
quences to cooler temperatures and larger luminosities, where most BSGs lie on the H-R diagram (Kaiser et al. 2020).Such solutions typically require an excess amount of mixing beyond the convective core which produces large cores.They are the most satisfactory for hotter BSGs (Brott et al. 2011), but they appear to fail to match the very gradual observed decrease in the number of BSGs toward cooler temperatures (e.g.Bellinger et al. 2023).
In addition to single stellar evolution pathways, it has been proposed that binary interactions may produce long-lived BSGs, motivated by the fact that a large fraction of young massive stars are found in binaries where close interactions or mergers may occur (Podsiadlowski et al. 1992;Kobulnicky & Fryer 2007;Sana et al. 2012;de Mink et al. 2014).Interactions between two unevolved stars may produce rapidly rotating stars (de Mink et al. 2013), which may experience extra mixing, resulting in an extended main sequence.If one of the stars has already evolved beyond the main sequence, it may give rise to a star with a smaller helium core compared to single stellar evolution products, either after accretion of hydrogen from a nearby star (Braun & Langer 1995), or by merging with a hydrogen-rich companion (Vanbeveren et al. 2013;Justham et al. 2014).Such a star could then evolve to a blue supergiant in the core-helium-burning phase (Farrell et al. 2019).Finally, stars that are partially stripped may also spend some time in the BSG (Klencki et al. 2020;Laplace et al. 2020), although this effect seems to work best at low metallicity.
While it is uncertain which of these scenarios contribute to producing BSGs, they all concern different evolution pathways, which result in different internal stellar structures and core properties compared to expectations from classical stellar evolution theory.Thus, asteroseismology, which studies stellar pulsations whose frequencies are finely tuned by the internal structure of stars (Aerts et al. 2010b), is a promising tool to investigate the structure, and hence, the origin of BSGs.In recent decades, there has been much interest in using asteroseismology to probe the structure and origin of BSGs from single-star evolutionary channels (Stothers & Simon 1968;Saio et al. 2006;Saio 2011;Bowman 2020).In particular, Saio et al. (2006) and Daszyńska-Daszkiewicz et al. (2013) identified an instability strip of gravity (g) mode pulsations in BSGs for shell-hydrogen and corehelium-burning models that arise due to trapping by the convective hydrogen burning shell and near core chemical gradient, respectively.Additional work by Godart et al. (2009) and Ostrowski & Daszyńska-Daszkiewicz (2015) demonstrated that varying the physics, such as mass loss or including convective boundary mixing dur-ing the evolution of the main-sequence, can impact the presence of the chemical gradients and intermediate convective zones required to trap and excite modes in BSG models.Saio et al. (2013) and Ostrowski & Daszyńska-Daszkiewicz (2015) further discussed the excitation of pulsations in models before and after core helium ignition, indicating that they produce different pulsation periods and excited mode spectra.While Daszyńska-Daszkiewicz et al. (2013) argue that regular oscillation patterns are not excited in BSGs, Bowman et al. (2019a) demonstrate that, if observed, g mode period spacing patterns can be used to discriminate between BSGs crossing the Hertzsprung-gap and those experiencing blue loops.Recently, Bellinger et al. (2023) made the first predictions of the asteroseismic properties of BSGs that result from binary evolutionary channels.They further demonstrate that asteroseismic signals may help to distinguish different single star an binary evolution BSG formation scenarios.Crucially, several of these studies predict the excitation of a dense spectrum of g mode pulsations with frequencies between 0.1 and 1 day −1 (Saio et al. 2006;Daszyńska-Daszkiewicz et al. 2013;Ostrowski & Daszyńska-Daszkiewicz 2015).
In the past decades, several space missions have delivered high-precision time series measurements for tens of thousands of stars, including CoRoT (Auvergne et al. 2009), Kepler/K2 (Borucki et al. 2010) and the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015).While a handful of samples of BSG variability has been discovered with CoRoT and Kepler/K2 photometry (see, e.g., Aerts et al. 2010bAerts et al. , 2017Aerts et al. , 2018;;Sánchez Arias et al. 2023), the all-sky TESS mission offers a unique opportunity for a systematic study for BSGs due to its large field of view that covers hundreds of potential candidates.This has motivated recent work to investigate the variability of BSGs with TESS data.Specifically, Bowman et al. (2019aBowman et al. ( ,b, 2020) ) have identified a stochastic low-frequency (SLF) excess in both main-sequence OB stars and BSGs.This SLF feature is characterized by an increase in the overall amplitude in the periodogram toward lower frequencies, with a constant amplitude profile as the frequency tends toward zero.Bowman et al. (2019a) argue that this feature is the signature of internal gravity waves (IGWs) generated by core convection (Rogers et al. 2013;Shiode et al. 2013).However, this claim is contested.Specifically, Cantiello et al. (2021) and Schultz et al. (2023a,b) demonstrate that these signals can be caused by subsurface convective zones and rotationally modulated stellar winds launched at the surface.Furthermore, Anders et al. (2023) use 3D hydrodynamical simulations to argue that IGWs generated by core convection should not have observable amplitudes at the stellar surface.Despite the increased attention in recent years, it is clear that more observations are required to identify the causes of the various signals seen in BSGs and then interpret them.
In this work, we analyze the light curves of a sample of 20 BSGs in the Large Magellanic Cloud (LMC).We find a characteristic time variability signal in our sample that universally appears in LMC BSGs.The signal shows similarity to the low-frequency variability of other massive pulsators, suggesting they may have similar origins.The manuscript is organized as follows: We describe our sample selection procedure in Section 2 and present our results in Section 3. We compare our signals to the low-frequency variability of other systems in Section 4, and in Section 5 we discuss the possible physical origins of our signal.We conclude in Section 6.

SAMPLE SELECTION
We selected our sample from the 124 OB pulsators with the spectroscopic measurements in Serebriakova et al. ( 2023), who obtained their data with the Ultraviolet and Visual Échelle Spectrograph (UVES; Dekker et al. 2000) and the Fiber-fed Extended Range Optical Spectrograph (FEROS; Kaufer et al. 1999) instruments on board the UT2@VLT and MPG/ESO 2.2 m telescopes.Their sample is largely based on the massive stars in Bowman et al. (2019a) and Pedersen et al. (2019) with photometric variability detected in the first sector(s) of TESS data, and they carried out high-resolution (R ≳ 40000), high-signal-to-noise-ratio (S/N ∼ 100) two-epoch spectroscopy such that binaries were excluded.The targets all locate in the Southern Continuous Viewing Zone (CVZ-S) of TESS and its vicinity, with high-duty-cycle photometric data of between approximately 200 days and 1 year in duration.
Out of the 148 pulsators, we selected 41 BSG candidates on the spectroscopic H-R diagram, with 3.2 < log (L/L ⊙ ) < 4.0 and T eff < 30000 K, where L := T 4 eff /g is the spectroscopic luminosity (Langer & Kudritzki 2014).The more luminous stars may still be BSGs of higher masses, yet we exclude them from our candidates since their evolution stages are even less clear (see Section 4 for a comparison between our sample and a limited sample of hotter B stars).Bolometric corrections turn out to be difficult for these candidates since many of them do not have reliable measurements of extinction, hence we did not select our sample from their real luminosities.
As each TESS pixel has an angular size of 21 ′′ , stars around our targets could cause contamination to the light curve, and this is especially an issue when the star locates in crowded regions of the sky (Pedersen & Table 1.Names, Tess Input Catalog (TIC) numbers, spectroscopic data and TESS observing sectors of our BSG sample.T eff , log g and v sin i are taken from spectroscopic measurements from Serebriakova et al. (2023) (UVES preferred when measurements from both instruments available), where we quoted their v sin iSP for references.
Bell 2023).To exclude false signals caused by nearby stars, we used the LATTE package (Eisner et al. 2020) to obtain star fields for all our candidates from both TESS and Gaia Data Release 2 (Gaia Collaboration et al. 2018).We compared the star fields and only kept the candidates without any bright (defined by the difference of magnitudes ∆m < 4) neighbors in the same TESS pixel as them.After such treatment, 20 stars are left in our final sample.It turns out that they are all in the LMC (Figure 2), and their spectroscopic properties are listed in Table 1.
We used the Python package Lightkurve (Lightkurve Collaboration et al. 2018) to download the TESS light curves for our sample from the Mikulski Archive for Space Telescopes (MAST1 ), using the data processed by a pipeline developed by the Science Processing Operations Center (SPOC; Jenkins et al. 2016; for its validity in processing low-frequency signals, see discussions in Appendix A).We used Astropy to analyze the data (Astropy Collaboration et al. 2013, 2018, 2022).Our sample has 10-29 observing sectors available from TESS, with typical observation period of ∼ 3 years (sector 27-69, from July 5, 2020 to September 30, 2023).Their details are summarized in Table 1.All the TESS data used in this paper can be found in MAST: https://doi.org/10.17909/df38-ax53.

RESULTS
Here we present our results of identifying and characterizing the variability for the BSGs in our study.

Universal Signal in BSGs
Figure 3 shows the TESS data and periodogram of one typical example system in our sample: SK -67 171 (TIC 425084139).The results for all our sample can be seen from Figures 10-29 in Appendix B. While we do not see any significant individual modes from the periodogram of this particular system, we do identify a low-frequency (f < 2 d −1 ) variability signal from the Fourier transform of the time series, as shown.The signal goes to zero amplitude toward high and low frequencies, with a characteristic peak frequency at ∼ 0.25 d −1 .We see in Appendix B that this signal appears in all of the targets in our sample.We suggest that this signature may be a common, or potentially ubiquitous phenomenon in BSGs.We note that other works have previously identified a variety of both stochastic and periodic variability in hot, main-sequence and evolved stars (e.g., van Leeuwen et al. 1998;Saio et al. 2006;Bowman et al. 2019a).However, this is the first reporting of a signature that tends to zero amplitude toward higher and lower frequencies in a sample of BSGs.The high-frequency tail above the peak frequency in our signal looks very similar to a characteristic signal that has already been widely discussed in hot massive stars, commonly referred as the stochastic lowfrequency (SLF) variability or "red noise" (Bowman et al. 2019a(Bowman et al. ,b, 2020;;Szewczuk et al. 2021;Bowman & Dorn-Wallenstein 2022;Dorn-Wallenstein et al. 2022).We note that, however, the signal pattern in our sample is distinct as it has decaying power toward zero frequency, while the red noise signal mostly concerns the high frequency tail.As the red noise signal is usually fitted with a super Lorentzian profile (e.g., Kallinger et al. 2014), we are motivated to characterize our amplitude spectra with the following modified Lorentzian profile: which is a linear profile times a super-Lorentzian profile with three parameters (ν char , α 0 and γ), plus a background noise term W 0 .This profile captures the shape of the signal with a peak frequency ν peak = (γ−1) −1/γ ν char and a tendency toward zero at both high and low frequencies.It also resembles the high frequency tail of super-Lorentzian profile above the peak frequency.We smoothed our amplitude spectra with a window size of 0.1 d −1 and fitted the smoothed spectra with the above function using a nonlinear least squares method with the curve fit function in the Python package ScipPy (Virtanen et al. 2020).The fitting parameters we found are shown in Table 2 and an example fitting curve is shown in Figure 3 (for all other stars in our sample, see Figures 10-29).We found that our fitting profile matches the smoothed spectrum very well, with the background-noise term W 0 at least one order-of-magnitude lower than the characteristic amplitude term α 0 .The physical peak frequency ν peak is around 0.2 days, and the power-law index γ ranges from 2 to 5, revealing a universal pattern in all our stars with small dispersion.We note that, however, the profile is not physically motivated and we did not find any strong correlations between the fitting parameters and the spectroscopic properties of our sample, which may give insights into the physical origins of the signal.We confirmed that our fitting parameters are robust against tests with different choices of smoothing windows sizes.

Time Variability of the Signal
The above analysis was carried out with the full TESS data obtained for our BSG sample, across multiple sectors.With an average observing time of 27 days per sector, the full TESS light curve gives us a frequency resolution of 0.001 − 0.005 d −1 .We note that, however, the signal actually varies in time as well, as seen by analyzing each observed sector individually.
Figure 4 shows the amplitude spectrum obtained from individual sectors for the example BSG SK -67 171.Each bar represents an amplitude spectra from one sector, with the color mapping the amplitude at different frequencies (such that a lighter color region represents a peak at that frequency).The full amplitude spectra combining all sectors are shown at the bottom for reference.We see that the shape of the amplitude spectrum for individual sectors is different from each other and the full amplitude spectra.Some features in one sector are preserved in the next sector (e.g., in sectors 27 and 28), while some only return after a few sectors (e.g., in sectors 33 and 36), or never repeats.We note that as the frequency resolution of individual sectors is as low  2. Fitting parameters and the physical peak frequency ν peak = (γ − 1) −1/γ ν char for the modified Lorentzian profile (Equation 1) that best matches our signal.
as ∆f ∼ 1/27 d = 0.04 d −1 , we cannot identify any of the features as distinct modes excited.Nevertheless, we found similar variability of the features for all our samples.
As each TESS sector only observes for a period of ∼ 27 days, such time variability in our signal suggests that the physical mechanism related to them must be very short, occurring on timescales shorter than ∼ 30 days.If the signal is stochastically excited from the stellar interior (see discussions in Section 5.3), the driving mechanisms must happen on such short timescales as well, and the characteristic signals they excite must be at least equally short-lived.This put additional constraints on understanding the physics behind them.
Low-frequency photometric variability has been found and discussed in different kinds of massive stars on the H-R diagram, including hot OB stars (Bowman et al. 2019b) and evolved yellow supergiants (Dorn-Wallenstein et al. 2019, 2020).To look for their similarities and possible common physical origins, we analyze a limited sample of other LMC stars and compare them with our BSG sample.(2023).The data are acquired with TESS with similar analysis as described in 2. With this limited sample, we clearly see similarities between the variability of hot OB stars and our BSG sample.Some of the alpha Cygni variables and yellow supergiants resemble the the signal we found in BSGs, yet with much lower power, while some other alpha Cygni variables show a massive peak at very low frequency with no clear turnover.The red supergiants, on the other hand, form the most distinct class on the H-R diagram and typically show some multi-peak pattern in the periodogram, very different from BSGs.
Previous photometric analyses on OB stars and yellow supergiants mostly interpret the low-frequency variability as "red noise", characterized by some function that levels off at zero frequency (Bowman et al. 2019b;Dorn-Wallenstein et al. 2019, 2020), which is clearly not what we see.This is probably because these authors mostly carried out their time series analysis for light curves shorter than ≲ 100 days, such that the characteristic low-frequency turnover cannot be well resolved in Fourier space.Nevertheless, the existing theoretical explanations for their signals are still relevant to our sample, given the similarities between these signals, which is what we probe in the next section.

DISCUSSIONS
In this section, we discuss the possible physical mechanisms that may cause the signal we found in our BSG sample and assess their viability in explaining the signal that we discussed in this work.

Stellar Winds
It is believed that stellar winds launched from the surface may cause low-frequency variability in massive stars.These clumpy and inhomogeneous winds are caused by some radiative driving mechanism whose exact details are still under debate (Owocki & Rybicki 1984;Puls et al. 1996Puls et al. , 2006Puls et al. , 2008)).They might be responsible for the "red noise" signal (see, e.g., Aerts et al. 2018;Ramiaramanantsoa et al. 2018;Krtička & Feldmeier 2018, 2021;Bailey et al. 2024) or even the large macroturbulence observed in spectroscopy of massive stars (e.g., Simón-Díaz et al. 2017).Blomme et al. (2011) argue that if stellar winds cause some low-frequency signal, its frequency should be compatible with the stellar rotation frequency at which the winds are launched, and one can then compare the equatorial velocity v eq,derive := 2πν peak R star derived from the peak frequency of the signal with spectroscopic v sin i measurements of the star.To test this scenario, we derive the bolometric luminosity L bol of 12 BSGs in our sample with reliable G-band magnitude and extinction measurements from Gaia (Gaia Collaboration et al. 2018, 2023), with the bolometric correction table provided by the Mesa Isochrones and Stellar Tracks (MIST) grids (Dotter 2016;Choi et al. 2016), assuming a distance of 49.97 kpc (the distance of LMC).We then derived their radii R bol with the spectroscopic T eff from Serebriakova et al. (2023), and calculated their derived equatorial velocities based on the fitted characteristic frequency in the signal.The results are shown in Table 3.Despite the uncertainties we introduced in this process, we found that v eq,derive for these stars are generally one order of magnitude larger than their measured v sin i, suggesting their inclinations are all less than ∼ 10 degrees.Hence we do not believe the signal is likely caused by stellar winds.
We note, however, whether v sin i measurements from spectroscopy should be interpreted as stellar rotation is actually unclear.For example, Simón-Díaz & Herrero (2014) and Simón-Díaz et al. (2017) point out that macroturbulent velocities may contribute to the total broadening of spectral lines in OB stars at least as equally as rotation, such that the measured v sin i overestimates the real stellar rotation rate (Schultz et al. 2023b).Nevertheless, this picture only disproves the  wind scenario more since it predicts even lower peak frequency compared to what we found from the signal.

Subsurface Convection
Historically, stochastic and nonperiodic variability is mainly discussed in solar-like oscillators and red giants, in the context of its association with surface convection and granulation (Schwarzschild 1975;Michel et al. 2008;Chaplin & Miglio 2013;Kallinger et al. 2014;Hekker & Christensen-Dalsgaard 2017).In massive stars, subsurface convective zones caused by opacity peaks associated with iron and helium ionization may create convective motion in the stellar photosphere, observed as low-frequency photometric variability (Cantiello et al. 2009b(Cantiello et al. , 2021)).
To see if these subsurface convective zones are related to our signal, we looked into a 15 M ⊙ BSG model from Bellinger et al. (2023) with a metallicity of 0.008 (a typical value for LMC stars), made with the Modules for Experiments in Stellar Astrophysics (r23.05.1, Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015Paxton et al. , 2018Paxton et al. , 2019;;Jermyn et al. 2023).The model starts as a post-merger star with a helium core and a hydrogen envelope, and it evolves to a BSG during its central-helium-burning phase, with a BSG lifetime of  6).The upper panels are in linear coordinates and the lower panels are in logarithm to make the core structures easier to see.Left: zone types of the star.We see two subsurface convective zones and a convective core.Before 0.6 Myr there is also a convective zone above the shell-burning region; Right: the opacity and nuclear generation rate of the model.The opacity profile shows two major peaks corresponding to the subsurface convective zones, and nuclear energy indicates that the two additional convective regions are caused by shell and core burning.
1 Myr (see Figure 6 for a comparison between the model and our sample on the H-R diagram).While some authors argue that single stellar evolution could also produce BSGs (see, e.g., Walborn et al. 1987;Weiss 1989;Schootemeijer et al. 2019;Kaiser et al. 2020), the BSG subsurface structures should be similar in those alternative models, since the occurrence of element opacity peaks is mostly determined by temperature.
In Figure 7 we show the structure (zone types, opacity and nuclear burning rates) of our model throughout the BSG phase.We see two subsurface convective zones, corresponding to the opacity peaks caused by helium and iron ionization, at log(T ) = 4.7 and 5.3.We note, however, that the subsurface convective zones in our model are much deeper than those found by Cantiello et al. (2009b), who mostly concerned less evolved main-sequence stars.Our convective zones are typically ≳ 0.5 R ⊙ below the stellar photosphere, one order of magnitude deeper than their findings (see Figure 2 in Cantiello et al. 2009b for comparison).This means subsurface convection is less likely to be observed as photosphere variability for BSGs.Nevertheless, as we used mixing length theory (MLT) to treat convection in our 1D stellar evolution model, we may seriously underestimate the spatial extension of convection, as confirmed in some recent 3D simulations on subsurface convection (Schultz et al. 2022(Schultz et al. , 2023a,b),b).In reality, the subsurface convective zones may extend to the surface to be observed.
If subsurface convection causes low-frequency variability in some signals, Cantiello et al. (2021) points out that the characteristic frequency of the signal should match the convective turnover frequency on top of the subsurface convection zone.In the left panel of Figure 8 we plot these frequencies of the convective zones in our BSG models (including two subsurface zones caused by opacity peaks, a convective zone above the burning shell before 0.6 Myr, and burning the convective core, see Figure 7 for reference), calculated by ν turn = (2πα MLT HP /v c ) −1 (with α MLT = 2.0, see Section 3 in Cantiello et al. 2021 for details), along with the typical peak frequency in our signal ν peak = 0.2 d.The results seem to favor the convective zone related to iron opacity peak, whose turnover frequency on top is only lower than the signal frequency by a factor of a few, and the inconsistency may only be caused by the uncertainty of α MLT (Cantiello et al. 2021).The frequencies of other convective zones are either too low, or the zones themselves are too deep to create any surface convective motion.However, as 3D simulations found that convective zones could extend spatially to larger scales than MLT predictions (Schultz et al. 2022), we should keep in mind that the two subsurface convective zones may actually merge together, and make our MLT calculations unreliable.
In addition, we note that (sub)surface convection and granulation are mostly discussed to explain red noise signals (see, e.g., Kallinger et al. 2014), which is different from the characteristic shape of our signal, as the former usually levels off at zero frequency, where our amplitude spectra vanishes.Physically, one would not expect incoherent convective motion to have some longest cutoff periodicity, creating a decay of power toward lower frequency.
The sector-to-sector time variability of our signal is also difficult to be explained by convection.3D hydrodynamical simulations have predicted that (sub)surface convection usually leads to a steady-state time granulation-like form of SLF variability (Schultz et al. 2022), unlike the long-term modulation timescales on the order of weeks to months in our signal.This means convective motion alone is not satisfactory to explain the shape and variability of our signal in general.

Internal Gravity Waves and Pulsating g-modes
Internal oscillations have also been suggested as origins of low-frequency photometric variability for massive stars.These oscillations could be waves that propagate to the stellar surface (Cantiello et al. 2009b;Bowman et al. 2019a;Ratnasingam et al. 2020), or pulsating modes that are trapped inside the stars (Shiode et al. 2013).They can be driven by the κ-mechanism associated with iron opacity peaks, which is responsible for nonstable pulsations in β Cephei stars and slowly pulsating B stars (see, e.g., Pamyatnykh 1999;Dupret 2001), or excited stochastically by turbulent convective motion (see, e.g., Lecoanet et al. 2021;Thompson et al. 2023).In BSGs with convective zones caused by the iron opacity peak, in principle both mechanisms could be related.Nevertheless, the stochasticity in our signal seems to favor convective excitation, and recent 3D simulations showed that modes excited this way can have lifetimes as short as ∼ 10 days (Thompson et al. 2023), consistent with our findings (see Section 3.2).
It has been proposed that several types of waves may propagate in stars that create low-frequency signals, including Rossby waves that are restored by Coriolis force (Van Reeth et al. 2016;Saio et al. 2018) and gravity waves that are restored by buoyancy (Cantiello et al. 2021;Schultz et al. 2022).Among them, we point out that damped gravity waves may naturally explain the shape of our signal at the low-frequency end: with shorter wavelengths at longer periods, gravity waves are preferentially damped at lower-frequency when they propagate in stellar radiative zones, creating smaller photometric amplitudes when they reach the surface.This is exactly what we found from our signal, and has been verified by recent 3D simulations on wave propagation (see Figure 2 in Anders et al. 2023).Goldreich & Kumar (1990) showed that if gravity waves are excited by convective motion, they should have a peak power around the convective turnover frequency on top of that convective zone.We hence showed these frequencies calculated for the different convective zones from our BSG model in the left panel of Figure 8, compared to the typical peak frequency in our signal ν peak = 0.2 d.We see that the convective turnover frequencies on top of the iron-opacity convective zone and the convective core are both within an order of magnitude lower than the observed peak frequency of the signal (an uncertainty tolerated by the choice of α MLT , see Section 5.2), suggesting they might generate a spectrum of waves that can be interpreted as our signal.We note, however, Anders et al. (2023) point out that waves generated from the core cannot get to the observed amplitude when they reach the surface, though in their simulations, they did not have other convective zones above the core, which may amplify or further damp the wave amplitude and it travels through them.
To see whether gravity waves excited by convection can propagate, we show a propagation diagram (i.e., the parameter space of the angular frequencies and locations inside a star where waves can propagate) in the right panel of Figure 8, for our BSG model at a typical stellar age of 0.9 Myr.Gravity waves can only propagate in regions where their frequencies are below both the Brunt-Väisälä frequency (Väisälä 1925) and the Lamb frequency (Unno et al. 1989).We see that for l = 1 waves, the cutoff frequency above the iron opacity peak is below the peak frequency we find in our signal, which means waves are at our peak frequency evanescent.Nevertheless, Cantiello et al. (2009a) and Shiode et al. (2013) demonstrated that gravity waves excited by subsurface convection can have very high degrees (l ≳ 30), whose cutoff frequencies are above the peak frequency of our signal (see Figure 8 for the propagation of an example l = 40 wave), such that they may propagate to the surface.However, these high-degree oscillations may only produce small amplitudes because of geometric cancellation effects (see, e.g., Aerts et al. 2010a), and it is questionable whether they can produce the observed amplitude of our signal (up to a few percent of stellar luminosity).
The sector-to-sector time variability of our signal is compatible with an origin of stochastically excited waves, as the long-period beating patterns among many superimposed waves of various spatial scales would lead to modulation timescales of months to years.An additional constraint from the signal is that there are no individual modes identified for any of our BSGs around the characteristic peak frequency.Bellinger et al. (2023) proposed that if BSGs are formed by mergers, the l = 1 g-mode period spacing would be ∆Π 1 ≈ 20 − 30 min.To identify individual g modes around ν peak ∼ 0.2 d −1 , this requires the frequency resolution to be less than ∆ν = ν 2 peak ∆Π ≈ 10 −3 d −1 , marginally below the resolution of the full TESS light curves of our sample.For higher-degree g modes, their geometric cancellation effects make it impossible to identify them from photometry.Therefore, the lack of identified modes from the signal still could not exclude the possibility that they are a spectrum of individual g-modes.Nevertheless, we expect further observations with increased photometric precision (such as those with the PLATO mission; Miglio et al. 2017) and observations seeking for spectroscopic time variability will help to distinguish these modes.
We note, however, that our wave analysis is based on timescale arguments on this single 1D stellar evolution model.More detailed investigation including solving oscillation modes/3D modeling of convection should be carried out to verify this explanation in the future.

CONCLUSION
In this manuscript, we analyzed TESS light curves for 20 blue supergiants (BSGs) in the LMC, and found a characteristic signal in the low-frequency (f ≲ 2 d −1 ) range for all our targets.The signals show strong stochasticity across different TESS sectors, yet their full amplitude spectrum show a peak frequency around 0.2 d −1 , below which the amplitude tends to zero.We were able to fit the signal with a modified Lorentzian profile (Equation 1), and the fitting parameters have no strong correlations with spectroscopic parameters measured for these systems.
We compared our signals to those obtained from a limited sample of hotter OB stars, yellow supergiants, alpha Cygni variables and red supergiants.We see similarities between our signal and the low-frequency variability of hot OB stars, yellow supergiants and some alpha Cygni variables.This comparison, while limited, may suggest the origins of these signals are similar.
We discussed three possible mechanisms that may explain this signal: stellar winds launched by rotation, subsurface convective motion, and internal oscillations of the star.The spectroscopically measured v sin i seems to disfavor the wind mechanism, as it would produce too low a characteristic frequency compared to the signal.We created a BSG model to test the latter two scenarios, and we found that the turnover frequency on top of the convective zone caused by the iron opacity peak is consistent with the peak frequency of the signal, despite the uncertainties introduced by the choice of mixing length parameter.While convective motion in this zone may be interpreted as the signal, it might be too deep to be directly observed in the photosphere, and such motions typically do not predict the shape and short-term time variability of our amplitude spectra.High-order, damped gravity waves excited from this convective zone seem to be a better explanation, which explains the peak frequency, the shape of the signal and the short-term time variability, though this picture will require more thorough investigation.
We restricted our sample to 20 BSGs with reliable TESS time series and spectroscopic parameters in the LMC.While this gives a clean sample to draw away the signal, we do not know whether it is still present for higher-metallicity BSGs.In future work, we will extend our analysis to more Galactic BSGs.A greater population of samples may also help us to seek potential correlations between the signal and stellar parameters.
Our theoretical analysis is somewhat limited by the single BSG model we chose, and future work should also focus on more detailed modeling of the BSGs in our sample.
We thank the anonymous referee for the constructive report that helps to improve this work.We thank Dominic Bowman, Jim Fuller, Matteo Cantiello, Stephen Justham, Lars Bildsten, Yanqin Wu and Norbert Langer for useful discussions.This work was initiated at the Kavli Summer Program in Astrophysics 2023, hosted at the Max Planck Institute for Astrophysics.We thank the Kavli Foundation and the MPA for their support.CJ gratefully acknowledges support from the Netherlands Research School of Astronomy (NOVA).
Facilities: Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute.The specific observations analyzed can be accessed via https://doi.org/10.17909/df38-ax53.STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS526555.Support to MAST for these data is provided by the NASA Office of Space Science via grant NAG57584 and by other grants and contracts.This work presents results from the European Space Agency (ESA) space mission Gaia.Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC).Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia Mul-tiLateral Agreement (MLA).The Gaia mission website is https://www.cosmos.esa.int/gaia.The Gaia archive website is https://archives.esac.esa.int/gaia.

A. IMPACT OF LIGHT CURVE REDUCTION ON LOW FREQUENCY SIGNAL
In this work, we made use of the TESS two minute cadence pre-search data conditioning simple aperture photometry, or PDCSAP, data downloaded from MAST that was reduced and treated using the SPOC pipeline (Jenkins et al. 2016).This pipeline removes systematic trends introduced by various instrumental effects.Specifically, the SPOC pipeline does bias, flat-fielding and background corrections and removes trends that are introduced by Earth-shine events, thermal brightening events, loss-of-fine-pointing events, changes in pixel sensitivity, and regular momentum dumps.Generally, these events introduce noise at low frequencies that can hide astrophysical signals at low frequencies.To mitigate these events, the SPOC pipeline utilizes co-trended basis vectors (CBVs) which encapsulate the common trends per CCD channel that are present in a sample of strategically selected photometrically quiet stars (Kinemuchi et al. 2012;Stumpe et al. 2012).These CBVs are fit to each star individually, and the pipeline decides which removal method is best.To investigate the impact that the removal of these trends has on the resulting signals in the periodogram, we investigate five light curves with different reductions.In particular, we investigate the raw SAP flux data, SAP flux data that has been corrected using co-trending basis vectors two different application methods (e.g.multi-scale regularisation and single-scale Gaussian prior), a custom pixel-level principal component analysis (PCA) de-trending, and the SPOC produced PDCSAP flux data.In the cases that use CBVs, we only use the first five CBVs.The different light curves and their accompanying periodograms are shown in the left and right columns of Fig. 9, respectively.The raw SAP and both CBV corrected light curves display clear trends that indicate that not all instrumental effects were properly treated.The impact of the remaining trends is seen as a low frequency noise excess below ∼0.1 d −1 in the top three rows.The PCA corrected light curve shows a clear decrease in the overall noise at the lowest frequencies, however, there are still trends at the beginning and end of observing sectors that are not totally  removed by the PCA method (demonstrated in the insets), resulting in some remaining low frequency noise.We note that only the PDCSAP SPOC light curve with the full calibrated reduction process sufficiently removes instrumental trends to explore low frequency variability.

Figure 1 .
Figure 1.Blue supergiants (blue stars) form a distinct class in the Hertzsprung-Russell diagram, different from Alpha Cygni variables (green stars), yellow supergiants (yellow stars), red supergiants (red stars) or hot OB stars (cyan stars).When we put black dots with equal time intervals on the stellar tracks obtained from single stellar evolution model (from Bellinger et al. 2023), BSGs lie in the region where the number of dots are orders of magnitude less than the dots on the main sequence, suggesting their occurrence rates should be orders of magnitude lower, inconsistent with the observed population.References: BSGs: 12 of the 20 stars in this paper with bolometric luminosity estimate; Alpha Cygnis: van Leeuwen et al. (1998); YSGs: Dorn-Wallenstein et al. (2020); RSGs: Neugent et al. (2020); hot OB Stars: from Bowman et al. (2019a) with bolometric luminosity estimate.

Figure 2 .
Figure 2. Locations of our BSG sample (red circles) in the Large Magellanic Cloud (obtained with the colored Digitized Sky Survey).The figure is made with SIMBAD (Wenger et al. 2000).

Figure 3 .
Figure 3. TESS data and amplitude spectra of SK -67 171 (TIC 425084139).The smoothed and fitted spectrum are shown in dashed and solid red lines.The amplitude spectrum shows a pattern with a peak frequency at ∼ 0.25 day −1 , and can be fitted by a modified Lorentzian profile (see main text).

Figure 4 .
Figure 4.The amplitude spectrum from individual TESS sectors for BSG SK -67 171.Each bar represents an amplitude spectra from one sector, with the color mapping the amplitude at different frequencies (such that a "lighter" color region represents a peak at that frequency).The full amplitude spectra combining all sectors are shown in the bottom for reference.We see clear time variability of the spectra across different sectors.

Figure 5 .
Figure 5.The amplitude spectra of 12 BSGs in the 0 − 2 d −1 frequency range plotted on their locations on the H-R diagram, compared to other types of stars (see Figure 1).We see similar kind of variability in hotter OB stars and around half of the Alpha Cygni variables, suggesting they may share similar origins.References: BSGs: 12 of the 20 BSG in this paper with bolometric luminosity estimate; Alpha Cygnis: van Leeuwen et al. (1998); YSGs: Dorn-Wallenstein et al. (2020); RSGs: Neugent et al. (2020); hot OB Stars: 22 stars from Bowman et al. (2019a) with bolometric luminosity estimate and spectroscopy measurements from Serebriakova et al. (2023).

Figure 6 .
Figure 6.The stellar track of a 15 M⊙ BSG model with a metallicity of 0.008 on the H-R diagram.This is a central helium burning star formed after a stellar merger (Bellinger et al. 2023), with a BSG lifetime of 1 Myr (highlighted with thick blue line).The positions of 12 BSGs with derived bolometric luminosities in our sample are shown as blue stars.

Figure 7 .
Figure7.The stellar structure upon radius of our 15 M⊙ model throughout its evolution in the BSG phase (thick blue line in Figure6).The upper panels are in linear coordinates and the lower panels are in logarithm to make the core structures easier to see.Left: zone types of the star.We see two subsurface convective zones and a convective core.Before 0.6 Myr there is also a convective zone above the shell-burning region; Right: the opacity and nuclear generation rate of the model.The opacity profile shows two major peaks corresponding to the subsurface convective zones, and nuclear energy indicates that the two additional convective regions are caused by shell and core burning.

Figure 8 .
Figure 8. Left: The convective turnover frequencies for the convective zones in our BSG model.The frequencies corresponding to the iron opacity peak convective zone and core convection are the closest to the peak frequency of the signal.Right: The propagation diagram for the stellar model at 0.9 Myr after merger, where the colored/gray shaded regions show the frequencies where l = 1 and l = 40 gravity waves can propagate.The regions where the Brunt frequency vanishes are the two subsurface convective zones.

Figure 9 .
Figure 9. Example light curves (left) and accompanying amplitude spectra (right) for TIC 425084139, with different reduction methods.

Figure 10 .
Figure 10.TESS data and amplitude spectra of SK -68 53 (TIC 31179797).The smoothed and fitted spectrum are shown in dashed and solid red lines.

Figure 11 .
Figure 11.TESS data and amplitude spectra of SK -66 125 (TIC 425086354).The smoothed and fitted spectrum are shown in dashed and solid red lines.
B. TESS LIGHT CURVES AND PERIODOGRAM FOR FULL SAMPLEIn Figures 10 to 29 we present the light curves and periodograms of all our BSG sample.They show a universal pattern in the low-frequency range.

Figure 12 .
Figure 12.TESS data and amplitude spectra of SK -67 283 (TIC 31511729).The smoothed and fitted spectrum are shown in dashed and solid red lines.

Figure 13 .
Figure 13.TESS data and amplitude spectra of SK -69 31 (TIC 30190076).The smoothed and fitted spectrum are shown in dashed and solid red lines.

Figure 14 .
Figure 14.TESS data and amplitude spectra of SK -67 275 (TIC 389864558).The smoothed and fitted spectrum are shown in dashed and solid red lines.

Figure 15 .
Figure 15.TESS data and amplitude spectra of SK -66 142 (TIC 276860494).The smoothed and fitted spectrum are shown in dashed and solid red lines.

Figure 16 .
Figure 16.TESS data and amplitude spectra of HD 269639 (TIC 287400996).The smoothed and fitted spectrum are shown in dashed and solid red lines.

Figure 17 .
Figure 17.TESS data and amplitude spectra of SK -67 7 (TIC 29987961).The smoothed and fitted spectrum are shown in dashed and solid red lines.

Figure 18 .
Figure 18.TESS data and amplitude spectra of SK -67 279 (TIC 31311824).The smoothed and fitted spectrum are shown in dashed and solid red lines.

Figure 19 .
Figure 19.TESS data and amplitude spectra of SK -70 31 (TIC 30534618).The smoothed and fitted spectrum are shown in dashed and solid red lines.

Figure 20 .
Figure 20.TESS data and amplitude spectra of SK -68 152 (TIC 389366376).The smoothed and fitted spectrum are shown in dashed and solid red lines.

Figure 21 .
Figure 21.TESS data and amplitude spectra of SK -67 88 (TIC 179638852).The smoothed and fitted spectrum are shown in dashed and solid red lines.

Figure 22 .
Figure 22.TESS data and amplitude spectra of HD 269721 (TIC 425084965).The smoothed and fitted spectrum are shown in dashed and solid red lines.

Figure 23 .
Figure 23.TESS data and amplitude spectra of HD 269510 (TIC 373682056).The smoothed and fitted spectrum are shown in dashed and solid red lines.

Figure 24 .
Figure 24.TESS data and amplitude spectra of SK -66 92 (TIC 373845622).The smoothed and fitted spectrum are shown in dashed and solid red lines.

Figure 25 .
Figure 25.TESS data and amplitude spectra of SK -67 151 (TIC 391809264).The smoothed and fitted spectrum are shown in dashed and solid red lines.

Figure 26 .
Figure 26.TESS data and amplitude spectra of SK -70 26 (TIC 30403638).The smoothed and fitted spectrum are shown in dashed and solid red lines.

Figure 27 .
Figure 27.TESS data and amplitude spectra of SK -67 171 (TIC 425084139).The smoothed and fitted spectrum are shown in dashed and solid red lines.

Figure 28 .
Figure 28.TESS data and amplitude spectra of SK -67 133 (TIC 287401176).The smoothed and fitted spectrum are shown in dashed and solid red lines.

Figure 29 .
Figure 29.TESS data and amplitude spectra of SK -67 72 (TIC 179038240).The smoothed and fitted spectrum are shown in dashed and solid red lines.

Table 3 .
ObjectT eff,spec (K) log(L bol /L⊙) R bol /R⊙ ν peak (d −1 ) v eq,derive (km/s) v sin i (km/s) (sin i) derive The derived bolometric luminosity, radius and equatorial velocity of 12 BSGs in our sample.The derived equatorial velocities are not compatible with spectroscopic v sin i measurements, typically one order of magnitude too large.