Photometric variability of blue supergiants: Constraints on stellar masses and evolutionary stages

Introduction

There are still important open questions about massive stars that profoundly impact our understanding of their future evolution. These questions involve the formation process of massive stars [85], the amount of mass lost via stellar winds [93], the interplay between mass loss and stellar pulsations [52, 24], the long-standing mass-discrepancy problem [44, 26, 80, 4], the fraction of binary mergers [55], among others. Massive stars of more than 8 $M_{\odot}$ remain as early B- or O-type stars in the main sequence (MS) for most of their lifetimes. Beyond this phase, the transition to red supergiants (RSGs) occurs gradually over several $10^{4}$ years at an almost constant stellar luminosity while hydrogen-shell burning becomes the dominant energy source. During this evolution, stars pass through the blue supergiant (BSG) phase. In general, stars with initial masses between 20 – 45 $M_{\odot}$ undergo a "blue-loop" phase [29], evolving back toward the BSG region during core helium burning.

In contrast, stars with initial masses below 20 $M_{\odot}$ may undergo blue loops before returning to the red region of the Hertzsprung-Russell (HR) diagram [78]. However, the mass limit separating these evolutionary paths depends on the mixing processes and mass-loss rate, two critical ingredients in stellar evolution, which are often uncertainly modeled. This makes it observationally challenging to identify whether BSGs are in a pre- or post-RSG phase since it cannot be understood from the location of the star on the HR diagram alone. Distinguishing between the two pathways might be plausible by studying the pulsation modes of BSG stars (Saio et al. 2013), through photometric and spectroscopic analyses [77, 53, 52, 78]. Specifically, this distinction relies on the fact that significant mass loss after going through the RSG stage significantly increases the luminosity-to-mass ratio ($L/M$) than during the first crossing. Saio et al. (2013) showed that radial pulsations are excited in the BSG region only in models that were previously RSGs, while, at a given effective temperature, these models also

Notes. The luminosity, $L$, and stellar mass, $M$, are in solar units.

exhibit a much richer spectrum of non-radial pulsation modes than their pre-RSG counterparts. These pulsations cause either slight variations or slight changes in the broadening of line profiles [82, 2, 1, 15, 27].

Typical oscillation modes observed in BSGs are p and g modes excited by the $\kappa$ mechanism at the Fe bump of opacity (see, e.g., Pamyatnykh 1999) and stochastic excitations caused by different mechanisms [64, 3]. For stars with $L/M$ (in solar units) exceeding $10^{4}$, the stellar envelopes are characterized by short thermal timescales and are dominated by radiation pressure [34, 40]. In this regime, the standard spectrum of acoustic modes is crossed by additional modes, the so-called strange modes, that cannot be explained by classical $\kappa$ and $\epsilon$ mechanisms [34, 37, 96]. The adiabatic approximation, which assumes an infinite thermal timescale, is no longer valid for describing the pulsational properties of such stars, and therefore a non-adiabatic stability analysis is needed [37, 38, 40] (Saio et al. 2013). These strange mode instabilities persist over the entire range of effective temperatures from red to BSGs and have been associated with variable mass-loss phenomena [39, 52, 43, 24, 40], with direct consequences for their evolutionary trajectory and final fate.

In this study, we aim to identify the oscillation modes and analyze the frequencies of selected southern BSGs as a means for constraining their properties. Nowadays, satellite data enable studies of light curves of massive stars over long time intervals with a temporal resolution of a few minutes, allowing for the detection of "micro" oscillations. A good example is the Transiting Exoplanet Survey Satellite (TESS) space mission [74]. The availability of high-precision TESS photometric data provides an excellent opportunity to investigate six BSGs (HD 52382, HD 74371, HD 79186, HD 80077, HD 99953, and HD 115842) selected from the sample of Haucke et al. [43]. All exhibit $L/M > 10^{4}$, in solar units, making them promising candidates for strange-mode pulsations. By analyzing their light curves, we aim to characterize their variability and derive constraints on their physical properties while minimizing the risk of identifying spurious periods or aliases [32, 11]. We also aim to compare the observed frequencies with the predictions from linear, nonadiabatic pulsation models. This comparison allows us to estimate the asteroseismic masses of the stars, which can then be contrasted with those derived from evolutionary models or spectroscopic techniques, thereby providing insights into the mass discrepancy problem.

The paper is structured as follows: the observations are detailed in Sect. 2. The time-series analysis methods and the models used in the linear stability analysis are presented in Sect. 3. Our main results and discussion are presented in Sects. 4 and 5, respectively. Additional data and detailed results are compiled in Appendix § 7. Section 6 summarizes our conclusions.

Target selection and TESS data

The objects listed in Table 1 were selected from Haucke et al. [43] based on their high $\log(L/M)$. The final sample comprises stars with effective temperatures between approximately 13 000 K and 26 000 K. These stars were selected to provide approximately uniform coverage across this temperature range. According to the variability classification of these objects, four – (HD 74371, HD 80077, HD 99953, and HD 115842) – are $\alpha$ Cygni variables. For these stars, the variations might not have a single explanation or even be caused by pulsations. Rather, they may result from some other surface-layer instability [42] (Kritička & Feldmeier 2021). The stellar parameters of our sample, such as the effective temperature ($T_{\mathrm{eff}}$), stellar radius ($R_{\star}$), surface gravity corrected by the centrifugal acceleration ($\log g_{\mathrm{cor}}$), and projected rotational velocity $v\sin i$ (where $v$ is the equatorial rotational velocity), were taken from Haucke et al. [43], who performed a quantitative spectroscopic analysis of the photospheric lines. All values are summarized in Table 1.

HD Number	Name	Sp. type	$T_{\mathrm{eff}}$ [K]	$R_{\star}$ [$R_{\odot}$]	$\log(L)$ [dex]	$\log g_{\mathrm{cor}}$ [dex]	$\log(L/M)$ [dex]	$v\sin i$ [km s$^{-1}$]
52382		B1 Ia	21500 $\pm$ 1000	21 $\pm$ 2.5	4.94 $\pm$ 0.13	2.48	4.29	55
74371	LN Vel	B6 Iab/b	13700 $\pm$ 500	73 $\pm$ 10	5.23 $\pm$ 0.14	1.81	4.17	30
79186	GX Vel	B5 Ia	15800 $\pm$ 500	61 $\pm$ 7	5.33 $\pm$ 0.11	2.02	4.21	40
80077	PV Vel	B2 Ia+e	17700 $\pm$ 1000	195 $\pm$ 47	6.53 $\pm$ 0.23	2.20	4.23	10
99953	V808 Cen	B1/2 Iab/b	19000 $\pm$ 1000	25 $\pm$ 2	4.87 $\pm$ 0.12	2.33	4.22	50
115842	ALS 3038	B0.5 Ia/ab	25500 $\pm$ 1000	35 $\pm$ 4	5.67 $\pm$ 0.12	2.76	4.30	50

For the photometric analysis, we retrieved high-cadence (120 s) TESS light curves from the Mikulski Archive for Space Telescopes (MAST) via the Python package Astroquery.Mast [36]. Targets were identified using their TESS Input Catalog (TIC) designation (v8.2, Paegert et al. [68], selecting the entry closest to the coordinates of each target. The time series correspond to PDCSAP_FLUX values processed by the Science Processing Operations Center (SPOC) pipeline, which applies co-trending basis vector corrections to remove systematic effects. The observations span three sectors for HD 52382 and HD 115842, four sectors for HD 99953, and six sectors for HD 74371, HD 79186, and HD 80077. Each sector was detrended independently by fitting a second-order polynomial to the time series. The normalized flux ($F_{\mathrm{n}}$) was then converted to differential magnitude via $\Delta m=-2.5\log(F_{\mathrm{n}})$, with uncertainties propagated from the flux errors.

Methodology

Time series analysis methods

For a detailed study of the light curves of massive stars, it is important to have signal analysis tools that complement each other to provide a complete analysis [5]. Therefore, in this work, we computed the Lomb-Scargle periodogram [60, 79], a well-known method for analyzing stationary signals. This method decomposes unequally spaced time-series data into a sum of sine and cosine functions, aiming to minimize the difference between observed and modeled data. On the other hand, the BSGs can exhibit erratic and quasiperiodic behavior leading to strong photometric and spectroscopic variations (see, e.g., [56, 25, 52, 51]). Therefore, we also employed the weighted wavelet Z-transform (WWZ), a method proposed by Foster [31], based on the Morlet wavelet function. Decomposing a time series into time-frequency space allows us to determine the dominant modes and how those modes vary in time [89]. The wavelet transform provides information on the time evolution of the parameters (frequency, period, amplitude, and power) used to describe periodic and quasiperiodic signals (see Figure 2). The WWZ scalogram is essential for interpreting the results, helping us distinguish between coherent and quasiperiodic frequencies. Coherent pulsation modes are expected to appear as constant frequencies in the scalogram, whereas fluctuations in power and frequency may arise from beating between close frequencies, non-sinusoidal modulations [5, 19], rotational effects [23] (Li et al. 2019), the intrinsic resolution of the WWZ technique itself [31], among others. A detailed description of the implementation of the methods and the criteria applied for the selection of frequencies is provided in Appendix § 7.

Stochastic low-frequency variability

For massive OB-type stars, a gradual increase in power toward low frequencies is ubiquitously observed in the frequency domain, an effect known as stochastic low-frequency variability (SLFV; [13, 72, 66, 41]). This type of variability is intrinsic to a star, but the origin of the signal is yet unknown, with several possible explanations proposed, such as inertia-gravity waves excited at the interface of the convective core [87, 51] (Bowman, D. M. et al. 2019; Bowman et al. 2019), stochastic waves generated by subsurface convection or granulation [13, 21], damped gravity waves (Ma et al. 2024), and line-driven wind instabilities [54].

The SLFV of the stars in our sample, with the exception of HD 52382 and HD 80077, was recently characterized by Kourniotis et al. [51] by fitting the amplitude spectrum with a $\gamma$-index power law:

where $f$ is the frequency, $A_{0}$ is the red-noise amplitude as $f \to0$, $\tau$ is the characteristic timescale of the structures producing the variability, $\gamma$ is the slope of the power-law decay, and $W$ is the white-noise level. The presence of SLFV in the amplitude spectra of these objects reflects stochastic, turbulent variability that must be accounted for in our frequency analysis. To better understand the variability of the sample and compare the characteristic parameters of the SLFV of the stars, they must be measured using the same methodology. Therefore, we calculated the SLFV parameters for HD 52382 and HD 80077 using the Kourniotis et al. [51] method. We refer the readers to that paper for a detailed explanation.

Our fittings are shown in Figure 1 and the resulting parameters for all the stars are listed in Table 2. The reported parameters reveal a clear bimodal structure in the power-law slope $\gamma$, which separates the sample into two physically distinct groups. The hotter, more compact supergiants (HD 52382, HD 99953, and HD 115842) exhibit steep slopes $\gamma$ and short characteristic timescales $\tau$. On the contrary, the cooler, more extended objects (HD 74371, HD 79186, and HD 80077) show flatter slopes and longer characteristic timescales. A further discussion of these results in a broader context is presented in Sect. 5.

ID	Var. type	$(P_{\mathrm{known}})$ [d]	$P_{\mathrm{TESS}}$ [d]	$\log(W)$	$\log(A_0)$	$\tau$	$\gamma$
HD 52382	IA	-	1.1, 2.3, 2.9<br>3.4, 8.7	−5.40	−3.80	0.05	4.04
HD 74371	ACYG, IA	2.9a, 5.2a, 7.7a<br>8.3b,f	2.5, 3.1, 4.1<br>5.4, 7.3, 8.8	−5.51†	−3.72†	0.13†	1.87†
HD 79186	IA	78.6c	3.2, 3.9, 5.0<br>6.1, 8.7	−5.30†	−3.11†	0.25†	1.58†
HD 80077	GCAS, SDOR<br>ACYG	3.1b, 21.2e, 41.5g<br>55.5g, 66.5g, 76.0g	3.5, 4.0, 5.4<br>8.6	−5.63	−3.60	0.20	1.59
HD 99953	ACYG, IA	17.7d	2.5, 3.1, 4.0<br>8.8	−5.26†	−3.75†	0.07†	3.80†
HD 115842	ACYG, IA	10.3b, 13.4f	1.4, 2.1, 4.0<br>6.0, 9.2	−5.51†	−3.58†	0.05†	3.58†

The linear nonadiabatic model

In this subsection, we present the linear nonadiabatic pulsation framework used to investigate the stability properties of stellar models. The first step in such a theoretical analysis is the construction of appropriate hydrostatic stellar models, which serve as the equilibrium configurations for the subsequent pulsation calculations (Ober et al. 1983; Gautschy & Glatzel 1990a,b; Kiriakidis et al. 1993). As the chemical composition, luminosity, effective temperature, and, at least, a guess for the mass of the objects are provided by the observations, it seems reasonable to calculate hydrostatic stellar models by an envelope integration from the photosphere into the interior of the star up to some conveniently chosen cut-off - temperature, typically on the order of $10^{7}$ K [97]. An envelope integration (solving the equations of mass conservation, hydrostatic equilibrium, and energy transport) poses an initial value problem with unique initial conditions specified at the photosphere, i.e., using the prescribed parameters: stellar luminosity, effective temperature, and mass, together with an estimate of the photospheric pressure. For the opacity, we adopted the latest OPAL tables from Iglesias & Rogers [46]. Convection was treated according to the stanNotes. Columns 2, 3, and 4 list the variability type, periods reported in the literature ($P_{\mathrm{known}}$), and photometric periods ($P_{\mathrm{TESS}}$) determined in this work, respectively. The SLFV parameters defined by Eq. (1) are listed in the last four columns. Parameters marked with a dagger (†) were taken from Kourniotis et al. [51] and are shown for comparison. Those for HD 52382 and HD 80077 were derived in this work. More details on the detected frequencies and periods are given in Table A.1. The variability classification was taken from Lefèvre et al. [58]: IA (poorly studied variables of early spectral type), GCAS ($\gamma$ Cassiopeiae type), ACYG ($\alpha$ Cygni type), and S Dor (S Doradus type). The periods reported in the literature were compiled from: van Genderen et al. [91], Lefèvre et al. [58], Underhill [90], Sterken [84], Knoechel & Moffat [49], Koen & Eyer [50], and van Leeuwen et al. [92].

dard mixing length theory (see, e.g., Böhm-Vitense [14]) using 1.5 pressure scale heights for the mixing length. With respect to the mass uncertainty for each object, a sequence of models was generated with a fixed luminosity and effective temperature (as provided by the observations), with the stellar mass being a free parameter within an interval of values comprising the observational guess (see Sect. 4.2). Subsequently, the envelope models were tested for linear nonadiabatic stability in response to radial perturbations. The associated perturbation equations and appropriate boundary conditions were accounted for as explained, for example, in Baker & Kippenhahn [9]. Convection was treated within the frozen-in-approximation (see also, e.g., Baker & Kippenhahn [10]).

As a result of the stability analysis, we obtained a discrete set of complex eigenfrequencies $\sigma$, whose real parts correspond to oscillation frequencies. Positive imaginary parts indicate damping; negative imaginary parts indicate growth and instability of a mode. For convenience, $\sigma$ is usually made dimensionless using the global free-fall timescale. The numerical solution of the boundary eigenvalue problem associated with the linear nonadiabatic stability analysis faces severe difficulties for stellar envelopes, which are characterized by small thermal timescales and a significant fraction of radiation pressure. Such models, to which the sequences considered here belong, exhibit strange modes together with violent instabilities with growth rates in the dynamical range. Therefore, we adopted the Riccati method (see Gautschy & Glatzel [33]) for a numerical treatment of the problem, which has been proven to provide reliable and highly precise results even in highly nonadiabatic environments.

Results

Frequency analysis

We analyzed the TESS light curves of the target stars to identify their independent frequencies. Significant frequencies were extracted through an iterative pre-whitening procedure. To complement the Fourier analysis and assess the temporal stability of the detected signals, we also applied the WWZ method. A detailed presentation of the corresponding periodograms and WWZ scalograms is provided in Appendix § 7.

HD 52382 (B1Ia) is a luminous supergiant located on the border of the constellation of Monoceros, the Unicorn. It is listed as an O-B2 runaway candidate (Tetzlaff et al. 2011). Pulsational activity was suggested based on the high observed macroturbulence (Haucke et al. [43]; Lefever et al. [57]); however, no pulsation periods have previously been documented for this star.

Our analysis of TESS sectors 7, 33, and 87 reveals five independent frequencies (whose associated periods are listed in Table 2), with most showing amplitude variations from sector-to-sector. Two ($f_1$ and $f_2$) also exhibit frequency variability (see Table A.1). Several harmonic and beating frequencies are additionally detected in the WWZ scalogram (Figure 2). The first independent frequency, $f_1 = 0.115\ \mathrm{d}^{-1}$ ($P \approx8.7$ days), is detected in sectors 7 and 87 but is absent in sector 33. The observed frequency difference between both detections remains within the Rayleigh frequency resolution and is interpreted as the same mode undergoing temporal amplitude and frequency modulation. The second frequency, $f_2 = 0.293\ \mathrm{d}^{-1}$ ($P \approx3.4$ days), shows a progressive decrease in frequency across the three sectors. In contrast, $f_3 = 0.343\ \mathrm{d}^{-1}$ ($P \approx2.9$ days), $f_4 = 0.432\ \mathrm{d}^{-1}$ ($P \approx2.3$ days), and $f_5 = 0.912\ \mathrm{d}^{-1}$ ($P \approx1.1$ days) show only marginal sector-to-sector scatter with no monotonic trend, suggesting more stable oscillations.

HD 74371 (B6 Iab/b) and HD 79186 (B5 Ia) are both located in the constellation of Vela. For HD 74371, a previous photometric study using HIPPARCOS data suggested a period of 8.291 days (Lefèvre et al. [58]), while van Genderen et al. [91] reported variations with periods of 2.9, 5.20, and 7.7 days, although these results were inconclusive. For HD 79186, Underhill [90] reported ultraviolet light variations and estimated a rotational period of approximately 78.6 days, proposing that the brightness changes could be due to the rotation of a mottled disk.

Our analysis of the combined TESS sectors 8 + 9, 35 + 36, 61 + 62, and 62 + 63 reveals a complex low-frequency variability pattern for both stars, with several independent frequencies identified alongside numerous transient peaks, beating frequencies, and sector-dependent signals. This is also evident in the WWZ scalograms, where several structures vary in amplitude and coherence over time. Given their similar spectral classifications, location, and comparable frequency distributions, both stars appear to share a similar variability behavior.

For HD 74371, six independent frequencies were identified between 0.11 and 0.40 d$^{-1}$, corresponding to periods between approximately 2.5 and 8.8 days, with the dominant frequencies concentrated between 0.136 and 0.246 d$^{-1}$. Several frequencies display moderate sector-to-sector variations in amplitude and frequency. Notably, $f_1$ ($P \approx8.8$ days), $f_2$ ($P \approx7.3$ days), $f_3$ ($P \approx5.4$ days), and $f_5$ ($P \approx3.1$ days) aligns the periods previously reported by van Genderen et al. [91]. Similarly, HD 79186 exhibits five independent frequencies distributed between 0.11 and 0.31 d$^{-1}$, corresponding to periods from 3.2 to 8.7 days. The frequency pattern is remarkably similar to that of HD 74371, especially for the dominant low-frequency component near 0.11 d$^{-1}$ and the sequence of peaks around 0.20–0.31 d$^{-1}$ (Figs. A.2 and A.3).

Both stars additionally display subtle visual microoscillations in their light curves over two consecutive sectors: HD 74371 in sectors 61 + 62 and HD 79186 in sectors 62 + 63. These oscillations fall below the adopted S/N cutoff and were therefore excluded from the frequency analysis. They are more likely associated with stochastic surface variability rather than coherent stellar p modes, since they do not persist across different sectors.

HD 80077 (B2 Ia+e) has been identified as an extreme supergiant, potentially a member of the open cluster Pismis 11, located at a distance of 3.2 kpc. With an absolute bolometric magnitude of $-10.4 \pm0.5$, this star ranks among the brightest B-type supergiants in our Galaxy (Knoechel & Moffat [49]).

If this classification is accurate, the star should exhibit substantial instability, as it lies above the Humphreys-Davidson limit (Humphreys & Davidson [45]). However, Carpay [22] reported that this star appears relatively stable, and the contrast between the stability of this extreme hypergiant and its neighboring luminous blue variables (LBVs) on the HR diagram warrants further investigation. Using HIPPARCOS and V-band photometric data, van Leeuwen et al. [92] identified significant periodicities of 66.5 days and 55.5 days, along with additional, less prominent peaks at 76.0 days and 41.4 days. Knoechel & Moffat [49] observed a period of 21.2 days through polarimetric analysis. Although no binary system has been confirmed, these authors suggested the presence of a compact companion that was originally more massive than the studied star. Additionally, Lefèvre et al. [58] reported light variations of 0.151 mag over a period of 3.11 days.

Our analysis of three TESS sector groups (8 + 9, 35 + 36, and 62 + 63) yields four independent frequencies. Similar to the behavior observed in the previous stars, the frequency spectrum is strongly affected by transient signals and sector-dependent amplitude variations (Fig. A.4). The dominant independent frequency, $f_1 = 0.116$ d$^{-1}$ ($P \approx8.6$ days), displays a noticeable frequency drift between sectors. The most stable signal corresponds to $f_2 = 0.185$ d$^{-1}$ ($P \approx5.4$ days), with the highest S/N value in sectors 8 + 9. Two additional frequencies were identified at $f_3 = 0.248$ d$^{-1}$ ($P \approx4.0$ days) and $f_4 = 0.287$ d$^{-1}$ ($P \approx3.5$ days), the latter also showing temporal variability. Although none of the TESS frequencies match directly with the long periods reported for this object, the frequency $f_4$ is close to the period reported by Lefèvre et al. [58].

HD 99953 (B1/2 Ia/b) is part of the Centaurus constellation. Haucke et al. [43] reported that the star exhibits a variable P Cygni H$\alpha$ profile, indicating that both the mass-loss rate ($\dot{M}$) and terminal wind velocity ($v_\infty$) vary by a factor of $\sim2.8$, which suggests significant wind variability.

Our analysis of TESS sectors 10, 37 + 38, and 64 reveals four independent signals between 0.11 and 0.41 d$^{-1}$, corresponding to periods between approximately 2.5 and 8.8 days. The frequency spectrum is strongly affected by transient and low-amplitude variability. The independent frequencies interact with each other and with harmonics, generating strong beating patterns, visible in the WWZ power spectra (Fig. A.5).

The dominant frequency, $f_1 = 0.113$ d$^{-1}$ ($P \approx8.8$ days), remains remarkably stable across all sectors. Similarly, $f_2 = 0.250$ d$^{-1}$ ($P \approx4.0$ days) and $f_3 = 0.323$ d$^{-1}$ ($P \approx3.1$ days) show only small sector-to-sector variations. In contrast, $f_4 = 0.398$ d$^{-1}$ ($P \approx2.5$ days) exhibits a frequency drift, decreasing from $\sim0.408$ to 0.385 d$^{-1}$ in almost four years (Table A.1).

The improved frequency resolution of the combined TESS sectors 37 + 38 reveals possible rotational splitting signatures. A triplet centered on $f_2$ is identified, with nearly symmetric separations of $\Delta f \approx0.048$ d$^{-1}$, corresponding to components at $f_- = 0.202$ d$^{-1}$ and $f_+ = 0.298$ d$^{-1}$. A similar separation is observed around $f_3$, with additional peaks at 0.2758 and 0.3726 d$^{-1}$. However, these side components are not recovered in other sectors. In addition, $f_2$, $f_3$, and $f_4$ also share a similar frequency separation between them of $\Delta f \approx0.07$ d$^{-1}$. Interestingly, all the frequencies found are very similar to those detected in HD 74371, which also exhibit the same separation of $\Delta f \approx0.07$ d$^{-1}$.

HD 115842 (B0.5 Ia/b) is a supergiant star located in the Centaurus constellation. Photometric variations with periods of 13.38 days and 10.31 days have been reported by Koen & Eyer [50] and Lefèvre et al. [58], respectively. The star also exhibits variations in the H$\alpha$ line profile [43] and present a large macroturbulent velocity $V_{\mathrm{mac}} = 225\ \mathrm{km\ s}^{-1}$ [30].

Our analysis of TESS sectors 11, 38, and 64 yields five independent frequencies; all of them display amplitude or frequency variations from one sector to another. Numerous beating structures are also visible in the WWZ scalograms, as shown in Fig. A.6. The first independent frequency is $f_1 = 0.109\ \mathrm{d}^{-1}$ ($P \approx9.2$ days), which varies from $0.094$ to $0.118\ \mathrm{d}^{-1}$. The period detected in TESS sector 38 ($P \approx10.64$ d) is in good agreement with the 10.31 day period reported by Lefèvre et al. [58]. Since the observed frequency variations remain within the Rayleigh resolution limit of a single TESS sector, as for HD 52382, $f_1$ is interpreted as the same oscillation mode undergoing temporal modulation. The frequencies $f_2 = 0.168\ \mathrm{d}^{-1}$ ($P \approx6.0$ days) and $f_3 = 0.252\ \mathrm{d}^{-1}$ ($P \approx4.0$ days) both exhibit systematic, monotonically increasing frequency shifts between consecutive sectors. This behavior differs from the more irregular pattern expected from SLFV alone and may instead indicate long-term modulation. In contrast, the higher-frequency signals $f_4 = 0.472\ \mathrm{d}^{-1}$ ($P \approx2.1$ days) and $f_5 = 0.710\ \mathrm{d}^{-1}$ ($P \approx1.4$ days) remain more stable across all sectors, making them the strongest candidates for coherent pulsation modes. To further assess the stability of the detected variability, we constructed phase-folded light curves, shown in Figure 3, for each independent frequency after pre-whitening all remaining independent signals.

Linear nonadiabatic analysis

The resulting periods for all the objects and their associated nonadiabatic model sequences are presented in Figure 4 as a period-mass ($P - M$) diagram, in terms of the lowest-order pulsation periods. At high masses, the eigenfrequency spectrum is well ordered and resembles that of a classical acoustic resonator (thin red lines). As the stellar mass decreases, or equivalently, as the $L/M$ increases, mode crossings become more frequent and give rise to the appearance of strange modes (thick blue lines). These unstable modes are characteristic of highly nonadiabatic stellar envelopes and are associated with increasingly violent instabilities, particularly when $L/M > 10^4$ [37].

As shown in Figure 4, all stars in our sample exhibit unstable strange modes with periods consistent with those observed (magenta circles). This diagram is used to constrain the current stellar mass (i.e., the asteroseismic mass). However, the solutions are not unique, and the possible asteroseismic mass intervals can be estimated when the observed periods (horizontal lines) intersect the loci of unstable modes within a narrow mass interval. Conversely, for a stellar mass derived from independent methods, the predicted periods of unstable modes can be compared with the observed periods, allowing the identification and classification of the detected pulsations.

In this work, we used the observed periods to estimate the stellar mass. We adopted uncertainties of about 0.5 $M_\odot$, which represent the spread of acceptable solutions (vertical bands) defined by the uncertainty in the periods (horizontal bands). In each case, the intersections cluster around a common mass value, producing an approximately vertical alignment in the diagram. This approach assumes that the periods predicted by the linear pulsation analysis are equal to the observed periods. Although this approximation is generally valid for moderately unstable modes, nonlinear effects may become important for modes with dynamical growth rates, particularly at the low-mass end of the model sequences. In such cases, the instabilities tend to inflate the envelopes when evolving into the nonlinear regime, (see, e.g., [97]) and, therefore, the final pulsation periods are longer than their linear counterparts. Consequently, the mass inferred from the observed periods in a linear analysis should be considered as a first approximation to the asteroseismic mass of the star.

The resulting mass estimates corresponding to the different solutions are listed in Table 3. For each mass interval (vertical bands), we adopted the mean mass. When different periods yielded overlapping mass intervals, they were assumed to correspond to the same solution and were grouped into a single mass range. The subsequent column indicates the number of observed periods, relative to the total number of detected periods used in the analysis, successfully matched by the unstable modes.

	HD 52382		HD 74371		HD 79186		HD 80077		HD 99953		HD 115842
$M \pm\sigma$ [$M_{\odot}$]	$N_P/T$		$M \pm\sigma$ [$M_{\odot}$]	$N_P/T$		$M \pm\sigma$ [$M_{\odot}$]	$N_P/T$		$M \pm\sigma$ [$M_{\odot}$]	$N_P/T$
3.90 ± 0.38	3/5		7.30 ± 0.52	4/6		9.37 ± 0.30	3/5		122.58 ± 0.25	1/4
4.75 ± 0.25	2/5		8.70 ± 0.42	2/6		10.06 ± 0.27	3/5		128.45 ± 0.25	1/4
6.08 ± 0.29	3/5		10.35 ± 0.37	2/6		10.67 ± 0.29	3/5		132.66 ± 0.35	3/4
7.20 ± 0.25	1/5		12.09 ± 0.60	3/5					138.37 ± 0.25	1/4
8.38 ± 0.25	1/5		13.85 ± 0.25	1/5					148.37 ± 0.25	1/4
3.80 ± 0.37	3/4
4.43 ± 0.25	2/4
5.54 ± 0.25	1/4
6.59 ± 0.25	2/4
7.67 ± 0.25	1/4
18.49 ± 0.25	1/5
20.08 ± 0.30	3/5
22.49 ± 0.25	2/5
24.86 ± 0.25	1/5
27.33 ± 0.25	1/5
34.90 ± 0.25	1/5
$M_{\mathrm{spec}}$	4.86 ± 0.98		12.56 ± 2.54		14.22 ± 2.86		219.92 ± 44.23		4.88 ± 0.98		25.72 ± 5.17
$M_{\mathrm{evol}}$	16.62 ± 0.96		21.12 ± 1.64		22.77 ± 1.63		87.26 ± 27.22		15.66 ± 0.97		29.17 ± 1.44
					11.23 ± 1.78						15.73 ± 1.87

For none of the stars does a single mass range fully explain all detected periods, indicating either that the asteroseismic solution is not unique or that multiple physical processes may contribute to the observed variability. For HD 52382 and HD 99953, the asteroseismic mass ranges are similar, yielding several discrete mass solutions between $3.8$ to $8.4\,M_{\odot}$. The masses that cluster around the largest number of detected periods correspond to $3.9 \pm0.4\,M_{\odot}$ and $6.1 \pm0.3\,M_{\odot}$ for HD 52382, and $3.8 \pm0.4\,M_{\odot}$ for HD 99953.

HD 74371 and HD 79186 display higher asteroseismic masses, between $7.3$ to $10.4\,M_{\odot}$ and $9.4$ to $13.9\,M_{\odot}$, respectively. For HD 74371, five of the six observed periods are associated with different sequences of unstable modes, whereas the remaining period does not intersect any unstable mode. A mass of $7.3 \pm0.5\,M_{\odot}$ can account for the majority of observed periods. For HD 79186, four discrete mass values are equally probable between $9.4$ and $12.1\,M_{\odot}$, leading to a mean value of $10.7 \pm1.4\,M_{\odot}$. Masses between $122\,M_{\odot}$ to $148\,M_{\odot}$ are predicted for HD 80077, being the most probable value of $132.7 \pm0.4\,M_{\odot}$, since it reproduces three of the four observed periods.

For the last object in the sample, HD 115842, the radial nonadiabatic models predict several possible masses between $18$ and $35\,M_{\odot}$, depending on the observed periods considered. Among all the solutions, the mass of $20.1 \pm0.3\,M_{\odot}$ reproduces three of the five observed periods.

Note: Asteroseismic mass ranges ($M_{\mathrm{ast}}$) are listed together with the fraction of matched periods over the total stable periods observed ($N_P/T$). The spectroscopic mass ($M_{\mathrm{spec}}$) is given for comparison. The evolutionary mass ($M_{\mathrm{evol}}$) corresponds to the position of the star on the HR diagram. Two values are provided when the star lies in a position where the blue loop can occur. Dispersions for $M_{\mathrm{evol}}$ were propagated from luminosity errors. For $M_{\mathrm{ast}}$, the uncertainty corresponds to half the adopted mass band (±0.25 $M_{\odot}$) used to identify period matches, extended by the spread of crossings within each range.

Discussion

Comparison of mass determinations and evolutionary implications

The mass discrepancy problem in massive stars was first identified by Herrero et al. [44], who performed a non-local thermodynamic equilibrium (NLTE) analysis of 25 Galactic O stars and found that spectroscopic masses were systematically lower than evolutionary masses. These masses are derived by comparing the position of the star on the HR diagram with theoretical evolutionary tracks. The discrepancies in both mass and helium abundances are significant, pointing to a failure in either the atmosphere or the stellar evolution models, suggesting that a higher mass loss rate at the beginning of the evolution should be adopted in the calculations to reduce the discrepancy. In subsequent decades, the inclusion of line-blanketed atmospheres, spherical stellar winds, and rotating evolutionary models substantially reduced the magnitude of the problem, though without fully resolving it (Crowther et al. [26]; Searle et al. [80]), particularly for stars below 30 $M_{\odot}$ (Markova et al. [63]). Moreover, Sánchez Arias et al. [78] suggest that the mass discrepancy problem may be part of a broader tension between spectroscopic diagnosis and evolutionary predictions. They show that even when stellar masses and luminosities can be quantitatively matched by evolutionary models, the observed CNO surface abundances may remain inconsistent with evolutionary predictions and highlight the need for large, homogeneously analyzed samples combining multi-epoch spectroscopy and long-baseline photometry. Additional recent studies have also linked the origin of this discrepancy to convective core boundary mixing, suggesting that evolutionary models require more massive helium cores than previously assumed to reproduce spectroscopic masses (Tkachenko, A. et al. [88]; Pedersen et al. [70]).

To evaluate a possible mass discrepancy in this sample of BSGs and establish a comparison with the asteroseismic mass estimates, we derived the spectroscopic masses from the stellar radii and surface gravities reported in Table 1:

The resulting spectroscopic masses are listed in Table 3.

To estimate the evolutionary masses, we interpolated the position of each star on the HR diagram into the nearest evolutionary tracks, using the rotating stellar evolution models computed with the Geneva stellar evolution code (Ekström et al. [29]) (see Figure 5). These models cover a range of initial masses from 16 to 250 $M_{\odot}$ at solar metallicity ($Z = 0.014$), with a ratio between the equatorial surface rotational velocity and critical rotational velocity of $\Omega= 0.4$. The obtained values are listed in Table 3.

For HD 52382 and HD 99953, the evolutionary tracks on the HR diagram suggest initial masses of 17 and 16 $M_{\odot}$, respectively (see Figure 5), leading to a current stellar mass of the same order (16.6 and 15.7 $M_{\odot}$). In contrast, the derived asteroseismic (4 $M_{\odot}$ and 6 $M_{\odot}$) and spectroscopic masses (5 $M_{\odot}$) are in good agreement with each other.

Haucke et al. [43] mention that these objects show signs of interaction with the interstellar medium (see their Figs. 2 and 4). This observation, along with several models that cross each other in the modal diagram as well as low spectroscopic mass estimates, may indicate that both stars are returning from the RSG phase, where they experienced shedding of material into their surroundings. As the adopted models cease to show a blueward evolution at that low-mass regime, possible explanations to support this scenario include enhanced mass-loss rates during the RSG stage (Georgy [35]; Saio et al. 2013) and erroneous distance estimates.

For HD 74371 and HD 79186, their positions on the HR diagram correspond to evolutionary tracks with initial masses of approximately 22 $M_{\odot}$ and 24 $M_{\odot}$, respectively (see Figure 5). At these locations, the stars may either be evolving toward the RSG phase or returning to the BSG phase.

When comparing the most probable asteroseismic mass of 7.3 ± 0.5 $M_{\odot}$ derived for HD 74371 with the spectroscopic mass (12.6±2.5 $M_{\odot}$) and the current predictions of evolutionary tracks (listed in Table 3), we find agreement with a post-RSG scenario, although some discrepancies remain. For HD 79186, the average asteroseismological mass of 10.7±1.4 $M_{\odot}$ closely aligns with the evolutionary mass of 11.23 $M_{\odot}$ expected in a post-RSG phase. This is also consistent with the spectroscopic estimate of 14.2 ± 2.8 $M_{\odot}$, within the uncertainties. Although HD 79186 was not previously classified as an $\alpha$ Cygni variable, such a classification is plausible given the observed variability and the mass ranges inferred from this analysis.

For HD 115842, evolutionary models indicate that the track with an initial mass of approximately $34\,M_{\odot}$ best matches the position of the star on the HR diagram (see Figure 5). The current predicted mass is either $29.2\,M_{\odot}$ or $15.7\,M_{\odot}$, depending on whether the star is in a pre- or post-RSG stage, respectively. The spectroscopic mass of $25.7 \pm5.2\,M_{\odot}$ is consistent, within the uncertainties, with the pre-RSG evolutionary mass. The most probable asteroseismic estimate of $20.1\,M_{\odot}$ is closer to the post-RSG evolutionary mass and is also consistent, within the uncertainties, with the spectroscopic mass. Therefore, both the pre- and post-RSG scenarios remain compatible within the available constraints.

In addition, according to Haucke et al. [43], combined WISE W3 and W4 images reveal a prominent bow-shock structure (see their Fig. 6), indicating a strong interaction between the stellar wind and the interstellar medium. Furthermore, Crowther et al. [26]; Searle et al. [80] reported that this star exhibits a nearly normal CNO abundance pattern for its spectral type, with a slight nitrogen overabundance. Taken together, these observational constraints suggest that the star has undergone significant mass loss through stellar winds after leaving the MS.

The peculiar case of the hypergiant HD 80077

This star exhibits the largest mass discrepancies of the entire sample. Its position on the HR diagram identifies it as a very massive post-MS star near the turning point of its evolutionary track before evolving back toward the hot side of the HR diagram. At that point, the star reaches its minimum effective temperature, placing it on the cool side of the bistability jump associated with the Fe recombination, where a denser and slower wind is expected. The evolutionary track most consistent with its position corresponds to an initial mass of $185\,M_{\odot}$ (see Figure 5), leading to a current predicted mass of approximately $87\,M_{\odot}$.

Instead, the asteroseismological model yields a stellar mass of about $133\,M_{\odot}$ from the simultaneous fit of three observed periods. The spectroscopic value is even larger, reaching $220\,M_{\odot}$, which is almost 2.5 times the evolutionary mass. While part of the discrepancy may be attributed to the large uncertainty in the luminosity arising from the poorly constrained distance, the adopted mass-loss prescription is another significant source of uncertainty, particularly for such highly evolved massive stars. For example, Vink et al. [94, 95] predict an increase in the mass-loss rate across the bistability jump, while Haucke et al. [43] observe a decrease and estimated a mass-loss rate of $5.4 \times10^{-6}\,M_{\odot}\,\mathrm{yr}^{-1}$ for this object.

As the disagreement among the independent mass estimates remains significant, such large discrepancies may indicate that the underlying physics of this object is more complex than that of a single evolving star. Despite the absence of a previously reported binary companion and the lack of evidence for binarity in the TESS data analyzed in this work, this object warrants a more detailed analysis. Particularly, determinations of helium abundances, and $N/C$ and $N/O$ surface ratios would provide strong constraints on the degree of internal mixing and mass-loss history.

SLFV parameters and their relation to stellar properties

Recent studies have revealed correlations between the SLFV parameters and the stellar properties of massive stars. For early-type stars spanning the upper MS, a negative correlation between the red-noise amplitude ($A_0$) and $T_{\mathrm{eff}}$ has been reported (Bowman et al. [17]; Shen et al. [81]). An inverse association between these two parameters was found by Kourniotis et al. [51] in a sample of 41 Galactic post-MS BSGs. Taken together, these results suggest that the SLFV may be amplified at or shortly after the end of the MS phase. Moreover, for BSGs, the characteristic timescale $\tau$ is found to be negatively correlated with $\log T_{\mathrm{eff}}$, while the slope $\gamma$ correlates most strongly with $\log(L/L_{\odot})$, after removing the $\alpha$ Cygni variables (Kourniotis et al. [51]). In addition, Spiecher et al. (2025) report no significant correlations between $A_0$ and the stellar properties for a sample of LBVs and high-luminosity BSGs, suggesting that these stars may be an extension of the $\alpha$ Cygni variables.

To test whether these trends persist in our sample, we compared the SLFV parameters of our six B supergiants spanning both pre- and post-RSG evolutionary stages, all of which exhibit strange modes, as a function of their stellar properties. In Figure 6, we show the SLFV parameters $\log A_0$, $\tau$, and $\gamma$ as functions of the stellar parameters $T_{\mathrm{eff}}$, $R$, $\log(L)$, and $\log g_{\mathrm{cor}}$. The two stars analyzed in this work (HD 52382 and HD 80077) are highlighted in blue, whereas the remaining four stars, whose SLFV properties were adopted from Kourniotis et al. [51], are shown in black diamonds.

No relation is seen between $\log A_0$ and any of the stellar parameters, with HD 79186 appearing as a persistent outlier with respect to $\log A_0$. For $\tau$, a tentative positive trend with $R$ and a negative trend with $T_{\mathrm{eff}}$ are observed when HD 79186 is excluded, suggesting that more extended and cooler supergiants tend to exhibit longer variability timescales. Regarding $\gamma$, the sample reveals a clear bimodal structure: the hotter, more compact supergiants (HD 52382, HD 99953, and HD 115842) exhibit steep power-law slopes, while the cooler, more extended objects (HD 74371, HD 79186, and HD 80077) show flatter slopes.

For none of the stars do the SLFV parameters show a clear segmentation between the post and pre-RSG candidates, suggesting that suggesting that the SLFV properties alone are not sufficient to distinguish evolutionary stages within this sample. As seen in Kourniotis et al. [51], no trends with the SLFV properties are observed in the α Cygni class. These results support the idea that α Cygni variables do not constitute a homogeneous class, neither in evolutionary state nor in the physical origin of their variability [51] (Saio et al. 2013; Spejcher et al. 2025). We note that the derived SLFV parameters are sensitive to the methodology employed [54, 71] (Ma et al. 2024). Model-independent tests, as proposed by Pedersen & Bildsten [71], may provide complementary constraints.

Conclusions

We analyzed TESS light curves of six BSG stars with a large $L/M$ to identify periods potentially associated with strange modes. Using hydrostatic stellar models, we tested their envelope structures for linear nonadiabatic stability under radial perturbations and compared the resulting theoretical pulsation periods with the observational data. This approach allowed us to estimate the most likely asteroseismic masses of either $3.9\,M_\odot$ or $6.1\,M_\odot$ for HD 52382, $7.3\,M_\odot$ for HD 74371, $10.7\,M_\odot$ for HD 79186, $133\,M_\odot$ for HD 80077, $3.8\,M_\odot$ for HD 99953, and $20.1\,M_\odot$ for HD 115842. These values were then compared with current evolutionary masses derived from the Geneva stellar models [29] and with spectroscopic masses obtained according to the stellar radii and $\log g_{\rm corr}$ estimates from

Haucke et al. [43]. From this simultaneous analysis, we find that HD 74371 and HD 79186 are evolving blueward following the RSG phase, while the evolutionary stage of HD 115842 cannot be uniquely determined, with both pre- and post-RSG scenarios compatible with the available mass estimates. We also find that HD 52382 and HD 99953 are likely in the post-RSG phase; however, their low masses are difficult to reconcile with the available evolutionary models. On the other hand, HD 80077 lies close to the turning point of its evolutionary track, exhibiting the largest mass discrepancy in our sample. This indicates that the underlying physics of this object is more complex than that of a single evolving star and warrants a more detailed analysis.

Overall, the SLFV parameters do not provide a clear diagnostic of the evolutionary stage of the stars of our sample, as no systematic trends were found. Nevertheless, the observed bimodal distribution in the characteristic timescales and power-law slopes suggests that the physical mechanisms driving the low-frequency variability may differ among α Cygni variables. Expanding this analysis to a larger sample will be essential to establish whether these trends are intrinsic to the class or reflect the limited size of the current sample.

In summary, we highlighted the potential of asteroseismology as a complementary tool for constraining the evolutionary status of BSGs. However, the mass discrepancies found for some objects indicate that further observational and theoretical efforts are still required. Moreover, since not all of the photometric periods identified in the TESS light curves find a counterpart in our theoretical radial pulsation models, an analysis of non-radial pulsations could help elucidate their physical origin and further characterize the rich variability observed in $\alpha$ Cyg stars. Such studies constitute important next steps toward a more complete understanding of the variability and evolution of BSGs.

Acknowledgements. This research work has been co-funded by the European Union, Project 101183150 - OCEANS and by the Program for Research and Innovation Horizon 2020 (2014-2020) under the Marie Skłodowska-Curie Grant Agreement No. 823734. AA and LC thank the financial support from the Universidad Nacional de La Plata (Programa de Incentivos 11/G160), Argentina. MK and JSA acknowledges financial support from the Czech Science Foundation (project GAČR 25-17532S). The Astronomical Institute Ondřejov is supported by RVO:67985815. AC acknowledge support from Centro de Estudios Atmosféricos y Astroestadística (CEAAS), Universidad de Valparaíso, Chile.

Frequency analysis

Lomb-Scargle periodogram

We performed a Lomb-Scargle periodogram analysis [60, 79] using the timeseries.LombScargle package from the Python module Astropy [8, 7, 6]. The resulting power spectral density (PSD) was converted into amplitude ($A$) using the relation

where $f$ is the frequency and $N$ is the number of data points in the light curve [28, 79, 66]. Significant frequencies were extracted through an iterative pre-whitening process. At each iteration, we selected the frequency with the highest amplitude, fitted the light curve with the corresponding sinusoidal model, and subtracted it. The residuals were then subjected to a new iteration to identify the next frequency. This process was repeated until no frequency exceeded the significance criterion. Baran & Koen [12] recommended a detection threshold of S/N $\geq5$ for TESS data to minimize the risk of spurious peaks. However, as these objects are strongly affected by SLFV, we adopted the traditional S/N $\geq4$ here. The signal-to-noise ratio was computed as the ratio between the amplitude of the extracted frequency and the mean amplitude of the residuals within a symmetric frequency window centered on the peak, excluding a narrow region around the peak itself. Burssens et al. [20] proposed using a $1\ \mathrm{d}^{-1}$ spectral window when significant SLFV is present in TESS data, to prevent the enhanced low-frequency noise from biasing the S/N estimate. When coherent pulsations coexist with SLFV, they suggested a broader $5\ \mathrm{d}^{-1}$ window to obtain a more representative noise level, although there is no strict consensus on this choice. We therefore adopted an intermediate value of $2.5\ \mathrm{d}^{-1}$ for each target.

Frequency uncertainties ($s_f$) were estimated following Kallinger et al. [47] as

where $T$ is the total time span of the observations and $\mathrm{sig}$ is the spectral significance approximated from the signal to noise ratio [73]. This formulation provides an upper bound on the frequency uncertainty, in contrast to the commonly used expression of Montgomery & O’Donoghue [65], which is derived under the assumption of white, uncorrelated noise and therefore represents a formal lower limit. In practice, Kallinger et al. [47] showed through numerical simulations that approximately 30% of cases exceed the Montgomery & O’Donoghue [65] estimate, and that the true uncertainty is further inflated in the presence of correlated noise. Given that the amplitude spectra of BSGs are characterized by SLFV, we adopted the Kallinger et al. [47] expression as a more conservative and physically motivated choice. Period uncertainties were propagated as $s_P=s_f/f^2$. The detected frequencies, derived periods, S/N values, and uncertainties for each object are listed in Table A.1.

Weighted wavelet Z-transform (WWZ)

For the wavelet analysis, we employed the LIBWWZ package available in Python. The wavelet transform provides information on the time evolution of the parameters (frequency, period, amplitude, and power) used to describe periodic and quasiperiodic signals. However, time-averaged values extracted from a wavelet analysis can be misleading due to the high sensitivity of the method to phenomena such as beating between close frequencies or multi-signal interference [89, 19, 5, 23]. For this reason, we relied exclusively on the WWZ scalogram to discuss the temporal behavior of the detected frequencies.

Selection criteria

Following the methodology employed by Kourniotis et al. [51], the extracted frequencies were categorized into four types using the Rayleigh resolution $1/\Delta T$, where $\Delta T$ is the time span of the data:

1. Independent frequencies, $f_i$. These are significant peaks physically interpreted as real signals, which could be either oscillation modes (g- or p-modes), surface spots, or orbital periods.

2. Harmonic frequencies, identified as integer multiples of the independent ones. Harmonics in light curves arise from their non-sinusoidal nature and reflect the nonlinear response of the stellar medium to the propagation of pulsation waves. Consequently, they are not classified as intrinsic pulsation modes and are presumed to replicate the characteristics of their independent parent modes. Strong harmonics are therefore a consequence of nonlinear pulsations, as expected in strange modes, but they also appear in high-amplitude non-radial pulsators. Also, harmonics can appear as a result of rotational modulation, as seen for massive stars when r-modes are present (Li et al. 2019).

3. Combination frequencies, defined as algebraic sums of two previously extracted frequencies. Separately, a frequency was flagged as a linear combination of others if the difference between the observed and predicted frequencies fell within the same $1/\Delta T$ criterion.

4. Unresolved frequencies: Two frequencies were considered unresolved if the difference between them did not exceed the $1/\Delta T$ spectral resolution [61].

It is important to note that after identifying two or three independent frequencies, most of the remaining frequencies can be explained as linear combinations. To avoid discarding potentially real signals that may accidentally satisfy the Rayleigh criterion and, therefore, be misclassified as combination frequencies or unresolved components, we searched for additional evidence by constructing phase diagrams. In this work, we define independent frequencies as physically meaningful oscillations. Therefore, frequencies were classified as independent only if they were detected in all sectors or, at least, in two of three non-consecutive sectors, and they exhibit an unambiguous structure in their phase diagrams. The independent frequencies under this criterion are listed in Table A.1.

Star $f_{\mathrm{id}}$	TESS sector			Avg. freq. $\pm\sigma$ [$\mathrm{d}^{-1}$]	Avg. period $\pm\sigma$ [$\mathrm{d}$]	S/N
HD 52382	7	33	87
$f_1$	0.1212(16)	$\cdots$	0.1093(19)	0.11525 $\pm$ 0.00841	8.676 $\pm$ 0.630	[9.2, 7.1]
$f_2$	0.3026(24)	0.2936(23)	0.2821(22)	0.29277 $\pm$ 0.01025	3.415 $\pm$ 0.121	[6.0, 5.7, 6.2]
$f_3$	0.3469(16)	0.3405(28)	0.3418(32)	0.34307 $\pm$ 0.00342	2.915 $\pm$ 0.029	[9.1, 4.7, 4.2]
$f_4$	0.4325(23)	0.4345(18)	0.4292(32)	0.43207 $\pm$ 0.00267	2.314 $\pm$ 0.014	[6.34, 7.52, 4.27]
$f_5$	0.9149(33)	0.9120(26)	0.9098(29)	0.91223 $\pm$ 0.00255	1.096 $\pm$ 0.003	[4.5, 4.6, 5.1]
HD 74371	8+9	35+36	61+62
$f_1$	0.1135(22)	$\cdots$	0.1127(16)	0.11310 $\pm$ 0.00057	8.842 $\pm$ 0.045	[4.6, 7.4]
$f_2$	0.1374(12)	0.1333(11)	0.1380(17)	0.13623 $\pm$ 0.00255	7.340 $\pm$ 0.138	[18.3, 12.3, 6.9]
$f_3$	0.1812(18)	0.1890(19)	0.1838(18)	0.18467 $\pm$ 0.00404	5.415 $\pm$ 0.118	[8.6, 5.2, 6.3]
$f_4$	0.2409(17)	0.2487(12)	0.2471(15)	0.24557 $\pm$ 0.00415	4.072 $\pm$ 0.069	[8.9, 14.3, 8.4]
$f_5$	0.3205(20)	0.3262(19)	0.3262(17)	0.32430 $\pm$ 0.00329	3.083 $\pm$ 0.031	[5.6, 5.7, 7.3]
$f_6$	0.3923(16)	0.3953(16)	0.3981(18)	0.39523 $\pm$ 0.00290	2.530 $\pm$ 0.018	[7.4, 5.3, 6.8]
HD 79186	8+9	35+36	62+63
$f_1$	0.1121(15)	0.1193(23)	0.1147(16)	0.11536 $\pm$ 0.00364	8.669 $\pm$ 0.273	[10.6, 6.4, 4.3]
$f_2$	0.1633(12)	0.1638(13)	$\cdots$	0.16353 $\pm$ 0.00034	6.115 $\pm$ 0.013	[13.6, 12.6]
$f_3$	0.2026(14)	0.2023(12)	0.2002(13)	0.20170 $\pm$ 0.00133	4.958 $\pm$ 0.033	[9.7, 13.3, 8.9]
$f_4$	0.2538(15)	0.2568(15)	0.2568(17)	0.25579 $\pm$ 0.00174	3.910 $\pm$ 0.027	[9.5, 8.8, 6.5]
$f_5$	0.3107(14)	0.3093(15)	0.3090(14)	0.30967 $\pm$ 0.00088	3.229 $\pm$ 0.009	[6.4, 9.4, 9.8]
HD 80077	8+9	35+36	62+63
$f_1$	0.1219(10)	0.1102(4)	0.1155(5)	0.11584 $\pm$ 0.00588	8.633 $\pm$ 0.444	[11.8, 16.23, 13.3]
$f_2$	0.1829(3)	0.1863(9)	0.1859(15)	0.18503 $\pm$ 0.00187	5.405 $\pm$ 0.055	[20.4, 7.7, 4.4]
$f_3$	0.2459(10)	0.2504(8)	0.2486(11)	0.24830 $\pm$ 0.00225	4.028 $\pm$ 0.036	[6.8, 8.8, 6.2]
$f_4$	0.2853(8)	0.2824(13)	0.2956(11)	0.28777 $\pm$ 0.00698	3.475 $\pm$ 0.084	[8.8, 5.3, 5.9]
HD 99953	10	37+38	64
$f_1$	0.1117(34)	0.1107(12)	0.1176(30)	0.11331 $\pm$ 0.00377	8.825 $\pm$ 0.293	[7.8, 11.8, 8.7]
$f_2$	0.2513(34)	0.2500(16)	0.2500(35)	0.25040 $\pm$ 0.00076	3.994 $\pm$ 0.012	[8.7, 7.3, 5.6]
$f_3$	0.3185(37)	0.3246(14)	0.3246(33)	0.32254 $\pm$ 0.00356	3.101 $\pm$ 0.034	[7.5, 7.8, 7.6]
$f_4$	0.4077(34)	0.3999(17)	0.3853(31)	0.39763 $\pm$ 0.01144	2.515 $\pm$ 0.072	[8.5, 6.4, 8.7]
HD 115842	11	38	64
$f_1$	0.1184(14)	0.09398(39)	0.1154(13)	0.10926 $\pm$ 0.01332	9.152 $\pm$ 1.116	[5.5, 5.0, 6.7]
$f_2$	0.1587(13)	0.16917(34)	0.1750(13)	0.16762 $\pm$ 0.00820	5.966 $\pm$ 0.292	[8.6, 6.7, 6.6]
$f_3$	0.2319(14)	0.25187(28)	0.2718(13)	0.25186 $\pm$ 0.01995	3.970 $\pm$ 0.315	[5.6, 8.7, 7.3]
$f_4$	0.4678(14)	0.4697(14)	0.4782(13)	0.47190 $\pm$ 0.00457	2.119 $\pm$ 0.021	[8.3, 7.9, 7.8]
$f_5$	0.7119(10)	0.7153(36)	0.7088(15)	0.71200 $\pm$ 0.00325	1.405 $\pm$ 0.006	[6.5, 6.3, 6.2]

Notes. The first column lists the name of the star and frequency ID ($f_{i a}$), ordered from lowest to highest. The listed values were obtained with the Lomb-Scargle periodogram for each TESS sector; uncertainties are given in parentheses for the last significant digits. The dispersion $\sigma$ is indicated for the average frequencies and periods. The S/N ratio corresponding to each sector is given in the last column.