© The Authors 2025

1. Introduction

β Cephei (β Cep) stars are core-hydrogen burning massive stars (8≲M/M_⊙≲25) pulsating in lower-order pressure (p-) or gravity (g-) modes. Although they have been known for more than a century (Frost 1902, 1906; Guthnick 1913) their mode excitation mechanism was only explained 30 years ago by Moskalik & Dziembowski (1992), Gautschy & Saio (1993), and Dziembowski & Pamiatnykh (1993). The pulsations are driven by the heat mechanism (also known as the κ mechanism) that is active in the partial ionisation zone of iron-like elements in the outer envelope. However, while mode excitation in β Cep stars is roughly understood, our knowledge is incomplete. More modes are detected in the data than are predicted by theory, as revealed in some of the prototypical class members (Pamyatnykh et al. 2004; Daszyńska-Daszkiewicz & Walczak 2010; Daszyńska-Daszkiewicz et al. 2013). Their observations show g-modes; instead, they are predicted to be stable modes (i.e. not excited). Moravveji (2016) provided a summary of observed modes in β Cep stars and illustrated the impact of higher opacities in the driving zone. Additional physical mechanisms may also be active in the interiors of such stars, while being absent in the current excitation theory. Microscopic atomic diffusion, notably radiative levitation, was recently shown to be one such missing ingredient (Rehm et al. 2024).

In addition to shortcomings in the understanding of the mode excitation, the mechanisms of mode selection are even less understood. It is still unclear whether the mode selection is similar over evolutionary timescales and why certain modes are excited to observable amplitudes while others are not. The first step to understand the selection of modes is the identification of their geometrical configuration, which remains challenging in β Cep stars given that their mode frequencies are not in an asymptotic regime. The geometry of the low-order modes of most β Cep stars is well described by a spherical harmonic, characterised by (l,m) representing the mode degree and azimuthal order. The identification of (l,m) is a prerequisite for asteroseismic modelling (Salmon et al. 2022). Our current work aims to improve this aspect of mode identification for this class of massive pulsators by relying on space photometry. Before diving into this topic, it is useful to recall the history and basic properties of β Cep star pulsations.

The interest in β Cep stars came in waves of discoveries. The earliest work focused on the period determination from radial velocity time series of the very few known stars similar to β Cep¹ and the lack of understanding of the instability mechanism (e.g. review by Struve 1952). In the following years, more β Cep stars were discovered and the efforts to understand the pulsation mechanism were continued, although the physical origin was still elusive (Lesh & Aizenman 1978). The discovery and explanation of the mode excitation mechanism (see Pamyatnykh 1999, for a summary) sparked new interest in β Cep stars. Aerts et al. (1992, 1994) conducted successful mode identification from time series spectroscopy for a few bright class members, while Heynderickx et al. (1994) identified the degree of the dominant pulsation modes for a substantial number of class members from multi-colour photometry. While both methods worked, only a limited number of modes per star could be identified, often just one or two. Recently, Shitrit & Arcavi (2024) showed that ground-based multi-colour photometry is still a viable way to achieve mode identification for the asteroseismology of β Cep stars, complementary to single-band modern space-based photometry.

Based on the numerous discovered β Cep stars, Stankov & Handler (2005) composed a catalogue containing 93 genuine and 77 candidate β Cep stars. However, the majority of them lack a solid identification of their pulsation modes. The β Cep stars discovered from large ground-based surveys (e.g. Pigulski & Pojmański 2008; Labadie-Bartz et al. 2020) or dedicated observations of a few selected young open clusters (Saesen et al. 2013; Moździerski et al. 2019) affected by the same limitation in terms of seismic probing power: too few or no pulsation modes have been identified.

Inferences of a β Cep star interior were achieved thanks to the 21-year photometric monitoring of HD 129929 (Aerts et al. 2003, 2004a; Dupret et al. 2004). These three dedicated studies kick-started β Cep asteroseismology from the identification of six modes: a radial mode, a rotationally split p-mode triplet, and two components of a g-mode quintuplet. It led to the first seismic probing of the internal rotation of a main sequence star other than the Sun, with the bottom of the envelope rotating about four times faster than the outer envelope. Inspired by this result, in coordinated efforts, months-long multi-site campaigns involving tens of astronomers worldwide were organised. Four of the brightest β Cep stars were scrutinised seismically from combined multi-colour photometric and spectroscopic time series data. These stars are still the best modelled β Cep stars today, thanks to the unambiguous identification of their modes: ν Eri (Aerts et al. 2004b; De Ridder et al. 2004; Pamyatnykh et al. 2004), θ Oph (Handler et al. 2005; Briquet et al. 2005, 2007), 12 Lac (Handler et al. 2006; Dziembowski & Pamyatnykh 2008; Desmet et al. 2009), and V2052 Oph (Handler et al. 2012; Briquet et al. 2012). The seismic estimates of the internal rotation, age (in terms of core hydrogen mass fraction), and core boundary mixing are summarised in Bowman (2020).

With the new era of space photometry, asteroseismology has gained a treasure trove of high-quality high-cadence data for many pulsating stars (see Aerts 2021; Kurtz 2022, for recent reviews). For β Cep stars in particular, the earliest space-based observations with the Microvariability and Oscillations of Stars (MOST) satellite revealed the multi-periodic behaviour of δ Ceti, which was thought to be a mono-periodic radial pulsator until then (Aerts et al. 2006). The CoRoT mission advanced our knowledge of only three more β Cep stars (Belkacem et al. 2009; Degroote et al. 2009; Briquet et al. 2011; Aerts et al. 2019), and the Kepler space satellite did not reveal any additional β Cep pulsators with identified modes (Balona et al. 2011; Lehmann et al. 2011). This led to a paucity of research on these stars with space-based time series photometry, although asteroseismology as a whole was revolutionised by the Kepler mission.

Already with the first two sectors of the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015) mission it became obvious that massive star asteroseismology had entered a new era (Bowman et al. 2019; Pedersen et al. 2019; Burssens et al. 2020; Balona & Ozuyar 2020; Sharma et al. 2022). Shi et al. (2024) provide a catalogue of all β Cep stars observed by TESS with 2 min cadence data. Southworth & Bowman (2022) identified eight β Cep pulsators in eclipsing binaries, while Eze & Handler (2024) recently expanded this list to 78 systems. However, to best study these stars, long-baseline time series observations are needed to resolve the anticipated close-by frequencies and to identify the pulsation modes.

Burssens et al. (2023) analysed the 352 d light curve of the bright β Cep star HD 192575 located in the TESS northern continuous viewing zone and were able to model its structure parameters with high precision as the fifth β Cep star with tight seismic constraints. They also estimated its internal rotation profile from forward modelling, but encountered large uncertainty due to the use of only three rotationally split multiplets. Vanlaer et al. (in prep.) meanwhile analysed the five-year TESS light curve (with two year-long gaps) and resolved some of the ambiguity in the mode identification, allowing the authors to perform a rotation inversion.

Next to CoRoT, Kepler, and TESS, the Gaia mission has proven to be a highly successful gateway to discover pulsating stars. Although it was not designed for asteroseismology, the recurring visits to millions of stars in a quasi-random pattern allowed the discovery of >100 000 new main sequence pulsators (Gaia Collaboration 2023). Hey & Aerts (2024) revisited these stars and analysed their first two years of TESS photometry, confirming the main pulsation frequency and classification for more than 80% of these Gaia-identified pulsators. Given these findings, we are encouraged to focus on the 222 β Cep stars in the sample of Gaia Collaboration (2023) and to follow up on the work of Hey & Aerts (2024) by considering all available TESS data.

Here we use TESS light curves with a considerably longer time base than Hey & Aerts (2024), thus enabling improved frequency resolution. Our aim is to search for resolved low (radial) order pressure and gravity modes and rotational splittings in the 222 Gaia-detected β Cep stars, opening pathways to mode identification and seismic modelling from combined TESS photometry and Gaia data, in a homogeneous way for an ensemble of β Cep stars that is as large as possible with publicly available light curves assembled with both satellites.

2. Data reduction

We selected all 222 variable stars classified by Gaia Collaboration (2023) as β Cep stars as our initial sample. Since our work focuses on the synergy between the TESS and Gaia DR3 light curves, we only included these stars. However, we emphasise that the methods developed in this work could potentially be applied to TESS observations of known β Cep stars with Gaia time series photometry. To prepare for the data reduction, we obtained for each star on our target list its position from the Gaia archive and cross-matched it with the TESS input catalogue (Stassun et al. 2019; Paegert et al. 2021) to obtain its TIC ID. Out of our initial sample of 222 stars, 216 have been observed by TESS and we can obtain time series photometry for them.

2.1. TESS data reduction and light curve detrending

We employed the TESS data reduction pipeline tessutils (Garcia et al. 2022) which enables the end-to-end creation of detrended light curves from TESS full-frame images (FFI) and is optimised for asteroseismic analysis. We briefly recall its work-flow here. For every TIC number, tessutils downloads cut-outs from the FFIs using tesscut (Brasseur et al. 2019). In the present work we chose a cut-out size of 25 px by 25 px. Subsequently, tessutils searches for the optimal aperture for each sector while trying to keep the contamination from neighbouring stars below 1%. Afterwards, it extracts the light curve of that TESS sector as the total flux within the chosen aperture. The number of sectors in the light curves range from one to seven with a median of four visits by TESS over the current mission lifetime spanning about five years.

Fifteen of our stars are located in very crowded fields near the Galactic plane or in open clusters. Since the targeted β Cep stars are typically among the brightest stars in the field they provide the main flux contribution. To reduce the influence of neighbouring stars, we set the aperture for these stars to a single pixel centred on the target's position (see Appendix A for an example). With this approach, we were able to obtain light curves for stars that would otherwise have (partially) been missed.

After aperture photometry, tessutils detrends the extracted light curves per sector with an automated principal component analysis with up to seven components to remove long-term trends in the data. This is justified for the asteroseismic analysis of β Cep stars as we are interested in variability with periodicities of typically several hours for the individual modes. The normalised light curves were binned to the longest cadence of the TESS photometry, that is 30 min. This guarantees a homogeneous treatment of all TESS FFI data of the sample stars detected by Gaia, while it leads to excellent agreement among the dominant frequencies detected in the TESS and Gaia DR3 light curves, as we show below. We tested that restricting to the higher-cadence TESS FFI data does not change the results for our β Cep population study. Finally, the data of all sectors were stitched to a single long-term light curve per star. However, despite the initial PCA detrending, we still noticed significant instrumental trends in many light curves. Therefore, we removed those parts of the light curves where instrumental effects were dominant and detrended the remaining light curve. The full description of the process and an example are given in Appendix A.

2.2. Frequency extraction from pre-whitening

For the extraction of frequencies, we used a pre-whitening procedure. This method consists in identifying the dominant signal frequency in the Fourier transform of the light curve, finding the amplitude and phase of it and subtracting the matching sinusoid from the light curve to construct a residual light curve. Subsequently, the next strongest frequency is identified in the residuals, and a sinusoid is built and subtracted. This procedure is followed until a user-specified stop criterion is reached (see, e.g. Chapter 5 in Aerts et al. 2010, for details).

In order to extract all significant frequencies up to the Nyquist frequency of 24 d⁻¹, we used a customised version of the pre-whitening module of the publicly available frequency analysis code STAR SHADOW (IJspeert et al. 2024). While we refer to that paper for details on the methodology, we point out that STAR SHADOW's pre-whitening procedure is computationally efficient because it does not rely on time-consuming non-linear regression for all extracted signals simultaneously. Rather, it only resorts to non-linear regression for those signals that strongly affect each other. On top of that, STAR SHADOW is well suited to analyse blended frequency peaks in Fourier spectra, which are quite common in our set of light curves due to the limited coverage of the baseline.

In each pre-whitening step, STAR SHADOW demands the Bayesian Information criterion (BIC) (Schwarz 1978) decreases by at least 2 to accept a signal as significant. In order to discourage overextraction of frequencies from the high-quality space photometry, we increased STAR SHADOW's BIC-threshold in each pre-whitening step to 10. Additionally, we employed a hybrid method for the selection of the next signal to be pre-whitened. Herein, the pre-whitening initially uses amplitude hinting (e.g. Van Beeck et al. 2021) and swaps to signal-to-noise ratio (S/N) hinting once a frequency is rejected, where the noise level is computed over a range of 1 d⁻¹ in the Fourier spectrum. If another potential signal is rejected, the pre-whitening is terminated. For the final significance assessment of each signal, STAR SHADOW demands the signal's S/N to be greater than the threshold defined in Eq. (6) of Baran & Koen (2021). We increased the S/N threshold by 0.25 for light curves with long gaps as recommended by Baran & Koen (2021).

2.3. Extracted frequencies

We found between four and 36 (median 13) independent signals in each of our TESS light curves. These values already show that β Cep stars pulsate with few modes. It also reveals that our detrending removed most of the instrumental effects, that would contaminate the Fourier domain and lead to more extracted yet spurious frequencies.

To analyse a large number of pulsators, stacking and sorting the periodograms can reveal hidden structures (e.g. Gilliland et al. 2010; Michel et al. 2017; Li et al. 2020a; Hey & Aerts 2024; Read et al. 2024). Figure 1 shows the extracted significant frequencies from the TESS photometry for our Gaia-observed sample sorted by the frequency of highest amplitude. Since β Cep stars show typically few excited low-order modes that do not occur in an asymptotic regime of high or low frequencies, no obvious structure emerges from this diagram. Nevertheless, this visual representation of the sample shows the range of detected significant frequencies. We observe that the distribution flattens towards low frequencies while declining more abruptly towards high frequencies. This is fully in line with instability computations of low-order modes in young and evolved β Cep stars (Pamyatnykh 1999). However, the sharp cut-off near 8 d⁻¹ may also be partly due to the natural transition towards the δ Sct star frequency regime adopted in the classification algorithms (Gaia Collaboration 2023). Further, it becomes obvious from the figure that the observed significant frequencies for a given star depend on the position of its dominant frequency.

Fig. 1. Stacked representation of the extracted independent, significant TESS frequencies for 212 of the Gaia-detected β Cep stars, sorted by the TESS frequency with the largest amplitude. The colour-scale indicates the amplitude of each detected signal on a logarithmic scale. This figure includes only the stars whose main frequency occurs in the plotted β Cep range. Three stars have lower main frequencies and turn out to be SPB stars or β Cep and SPB hybrid pulsators with dominant high-order g-modes, or are dominated by other variability than pulsations.

2.4. Comparison to Gaia frequencies

Figure 2 compares the main frequency from our TESS analysis with the frequencies derived from Gaia DR3 by Gaia Collaboration (2023). As already found by Hey & Aerts (2024), most of the pulsation frequencies agree within a tolerance of Δf<0.006 d⁻¹, highlighting the potential of the Gaia time series photometry to discover non-radial pulsators at a level of 4 mmag or above. For the twelve stars without a matching TESS primary frequency, we searched through other signals in order of decreasing amplitudes in the TESS data and compared them to the Gaia frequency. Once a match within our permitted tolerance was found we stopped the search. For seven stars², we found no matching signal in the TESS data. The median difference between the TESS and Gaia matched frequencies is $\mathrm{\Delta }\tilde{f}=0.0013$ d⁻¹.

Fig. 2. Comparison between the highest amplitude Gaia and TESS frequencies for the 215 stars with significant TESS frequencies. Top: Dominant TESS frequency fTESS against the Gaia derived frequency fGaia. The black symbols indicate stars for which we found a match in frequency, the blue symbols show stars matching a frequency other than the dominant signal, and the red symbols indicate the seven stars for which the Gaia frequency does not match any frequency in the TESS data. The lines connect the offset frequencies to the matched counterpart. For three stars the main frequency is below the plotting range. Bottom: Difference between the Gaia and TESS frequencies. The black dashed line indicates our chosen threshold for matching signals (Δf<0.006 d−1).

We limited the current analysis to the primary frequency detected with Gaia. Although we possess measurements of the secondary frequency from the Gaia DR3 time series photometry, we noticed that these frequencies match our extracted frequencies only for a minority of cases. Often, we find f_Gaia>10 d⁻¹, whereas we rarely extracted any significant frequency in this range from TESS. It was already discussed by Gaia Collaboration (2023) that the Gaia DR3 data suffer from instrumental effects at mmag level for the secondary frequencies, which may be connected to multiples of the satellite frequency of 4 d⁻¹. Hence for the β Cep population study in our work, using the secondary frequency from Gaia is not informative. We note that detailed work on certain stars might be able to derive additional constraints for the pulsation mode of the secondary frequency. Further, longer time series data delivered by future Gaia data releases might provide more robust secondary frequencies in agreement with TESS.

3. Stellar properties

3.1. Gaia derived properties

We aim to analyse the β Cep stars in our sample in a homogeneous way and hence chose to use stellar parameters from Gaia DR3. Of particular interest are the effective temperature (T_eff) and the luminosity (log L/L_⊙). Shi et al. (2023) show that the Gaia T_eff derived by the Extended Stellar Parametrizer for Hot Stars (ESP-HS, Fouesneau et al. 2023) is reliable for pulsating hot main sequence stars. We follow their example and use these values as the basis for our analysis. Further, we extracted positions, parallaxes, and G band magnitudes for every target. Where available, we also added stellar parameters from ESP-HS and the General Stellar Parametrizer from Photometry (GSPPHOT). Unfortunately, the Gaia ESP-HS parameters are only available for 190 of the 222 stars. Hence, we are not able to place all the observed stars in the Hertzsprung-Russell diagrams (HRD). Further, this limits the sample of stars that we can model asteroseismically.

From the Gaia data, we find the majority of our sample stars to be in the magnitude range G = 8.5−13.5 mag (cf. Table 1). Nine stars, including the well studied stars BW Vul and 12 Lac, are brighter than G = 8.5 mag while no star is fainter than G = 13.6 mag. To derive the stellar luminosity, we require a bolometric correction. The Gaia bolometric correction (Creevey et al. 2023) is only valid up to 10 000 K and β Cep stars are hotter with 15 000 K≤T_eff≤25 000 K. Pedersen et al. (2020) derived bolometric correction functions for massive stars including the range of β Cep stars. We used these and calculated the luminosities log L/L_⊙ based on the Gaia DR3 absolute magnitudes taking into account the reddening based on the ESP-HS estimate and the distance as determined by the GSPPHOT pipeline for OB stars, the general GSPPHOT values, or the inverse of the parallax depending on the availability.

The uncertainties of the parameters in the astrophysical parameters table are often underestimated (Fouesneau et al. 2023). In particular the effective temperature has an uncertainty much larger than the typically given 500 K (∼2%). For the purpose of this work, we assumed a conservative uncertainty of 10% (Δlog T_eff = 0.043 dex) for all stars in the sample, which is a reasonable estimate given the comparison in Fig. 15 of Fouesneau et al. (2023). The uncertainty of the luminosity is strongly dependent on the effective temperature through its use in the bolometric correction.

An estimate of the projected surface rotation frequency is useful when searching for rotationally split multiplets and examining differential rotation. The Gaia ESP-HS parameters also include the projected surface rotational velocity v sin i calculated from the line broadening of the RVS spectra under the assumption that rotation determines the line broadening, which is a valid assumption (cf., Aerts et al. 2023). Together with the luminosity and effective temperature, which gives an estimate of the stellar radius, we can calculate the projected surface rotation rate of the stars. This parameter, f_rotsin i, is available for 167 stars. Although it is not very precise with a typical uncertainty ∼40%, it still proves useful in our subsequent analysis.

3.2. Gaia data of known β Cep stars

To compare the properties of the Gaia-detected β Cep stars with these of known β Cep stars, we extracted their Gaia DR3 data as well. As the literature sample, we use the genuine β Cep stars from the catalogues of Stankov & Handler (2005), Pigulski & Pojmański (2008), Labadie-Bartz et al. (2020), Shi et al. (2024), and Eze & Handler (2024). For each of the catalogues, we derived the same Gaia properties as described above. From the four catalogues in the literature, we find 582 stars of which 408 are classified as genuine β Cep stars. Among with the 222 β Cep candidates in our sample – 145 are not listed in any of the literature catalogues – the total number of β Cep stars in the combined sample is 553 stars.

4. Population of β Cep stars

Before analysing our Gaia-detected sample in detail, we present the overall properties of the ensemble. Figure 3 shows all 442 stars (out of 553) with adequate Gaia ESP-HS parameters in an HRD. Most of the stars in our sample fall inside the instability strips for OB-type pulsators calculated by Burssens et al. (2020). These are based on one particular choice of input physics (relying on exponential core overshoot with parameter f_ov = 0.02) at solar metallicity, so we expect some stars to occur outside that region, as also stressed by Aerts & Tkachenko (2024), Rehm et al. (2024), and Hey & Aerts (2024).

Fig. 3. Hertzsprung-Russell diagram based on Gaia data of all 321 β Cep stars from the literature (circles) and the 121 new stars in our Gaia-informed sample (squares, new detections). The typical uncertainty is indicated in the lower left. The black outlines show the instability regions for p-modes (solid, β Cep, l = 0−2), and g-modes (dashed, SPB, l = 1−2) based on the calculations from Burssens et al. (2020). The main sequence evolutionary tracks of Burssens et al. (2020) are shown in grey. The colour-coding gives the main pulsation frequency of the stars cut off at 12 d−1. Stars marked with a red cross do not exhibit β Cep variability in TESS. The top axis marks the spectral type.

Several stars from the literature occur below the β Cep instability region. We checked these stars with publicly available light curves from the TESS Quick Look Pipeline (QLP, Huang et al. 2020) using a classification similar to our main sample and outlined below in Sect. 4.1. We find the majority of stars close to the instability region to be β Cep stars pointing either to imprecisions of the stellar parameters or a larger extend of the instability region as found for other pulsators in this area of the HRD too (Gaia Collaboration 2023; Mombarg et al. 2024). Stars further away from the instability region are either short-period eclipsing binary stars or SPB stars. The stars³ with obvious non-β Cep variability (and available light curves) are marked in Fig. 3 with red crosses. All the new Gaia-derived β Cep stars fall inside or close to the β Cep instability strip and can hence safely be assumed to be genuine β Cep stars from their HRD position and TESS variability characterisation (see below).

At log T_eff≈4.5, we find a gap between the bulk of the β Cep stars and a smaller group mostly containing stars from our new sample. With only the Gaia ESP-HS data, we cannot be certain whether this gap is astrophysical or introduced by the intricate parameter estimation procedure. It is present in the entire ESP-HS data set. We note the position of the gap overlaps with the transition from spectral type O to B (Pecaut & Mamajek 2013). In case the gap is physical the stars on its hotter side would be of great interest for follow-up observations and detailed asteroseismic modelling as these stars might be products of binary evolution (Fabry et al. 2022, 2023).

4.1. Detailed analysis of the Gaia sample

After verifying that all of our sample stars with Gaia ESP-HS parameters indeed fall into the parameter range of β Cep stars, we analyse their TESS light curves in more detail. In a visual inspection of the light curves and their Fourier transforms, three main sub-classes emerge (Fig. 4). The first class comprises of mono-periodic pulsators which show one dominant coherent oscillation in the entire light curve. The periodograms of these stars often reveal they pulsate in multiple frequencies, but their secondary frequencies have small amplitudes compared to the dominant mode. Clearly multi-periodic pulsators constitute the second class. Their light curves exhibit modulations due to mode beating with periods between one and tens of days, as already known from ground-based multisite campaigns prior to the space asteroseismology era. In some cases, the beating pattern is significantly longer than one TESS sector. Hence, multiple TESS sectors are needed to resolve the close frequency peaks. We distinguish between β Cep multi-periodic light curves (due to the presence of multiple low-order p- or g-modes) and hybrid pulsators with SPB components as a third sub-class because their secondary signals have a much smaller frequency, resulting in a differently shaped light curve. (Fig. 4). We manually classified all 216 light curves based on their shape and the content of their periodogram into these three main classes.

Fig. 4. Examples of single sector TESS light curves (left) of β Cep variables of the different light curve sub-classes and their associated periodograms (right). The top light curve is a typical mono-periodic pulsator with a constant pulsation amplitude and a single peak in the frequency spectrum. The second and third panels show multi-periodic pulsators with different lengths of their beating patterns. The bottom panel shows a typical β Cep and SPB hybrid pulsator.

In addition to the manual classification, we quantified the light curve shapes by calculating two measures of the amplitude variability. Firstly, we use the length of the main beating pattern given by

$$P_{\mathrm{beat}}={\left|f_1-f_2\right|}^{-1}.$$(1)

Secondly, we calculate the strength of amplitude variability

$$A_{\mathrm{var}}=A_{90}/A_{95}·0.95/0.90,$$(2)

in which A₉₅ (A₉₀) is the difference between the 95th (90th) and 5th (10th) percentile of the light curve. The percentile differences (A₉₅) are known to robustly estimate stellar variability in space photometry (Basri et al. 2011). We expand on this idea and use the ratios of the percentile differences. By normalising this ratio, a pure sine variability (i.e. a mono-periodic pulsator) reaches unity, while multiperiodic stars can be found at lower values. We used the full light curve itself to compute A_var rather than the amplitudes of the extracted frequencies because it includes all variability beyond the two main frequencies in this measure.

Figure 5 shows that we can clearly separate the above defined groups based on A_var and the beating period. The initial manual classifications are mostly confirmed. Stars with a dominant period in their variable light curves fall towards the right. Multi-periodic pulsators populate the upper left because they are not close to a single sinusoid, although they have a long beating pattern, while hybrids lie near the lower left because their primary and secondary frequencies differ greatly. Some stars fall into regimes in Fig. 5 that are not associated with their manual classification. Among the mono-periodic stars, several stars display additional rotational variability that modulated the light curve without creating a beating pattern, a feature that cannot be accounted for with the current parameters. Some stars have long beating patterns with low amplitudes which makes them appear mono-periodic. The grouping of mono-periodic stars at long beating periods indicate that these stars are indeed not truly mono-periodic but exhibit multi-periodicity with very small frequency differences, a well-known phenomenon of some classical β Cep stars (e.g. HD 129929, β CMa). Given the few differences between our manual classification and the regimes in Fig. 5, we provide A_var and P_beat in Table 1 to enable quantified distinctions of the different groups.

Fig. 5. Beating period of the two main frequencies (Pbeat) against the normalised percentile amplitude ratio (Avar). Light curves manually classified as mono-periodic are shown in purple, multi-periodic in olive, hybrids in red, and other cases in cyan. The dashed lines indicate the approximate dividing lines of the groups of pulsating stars. Uncertainties on the beating period are within the symbol sizes (Avar has no formal uncertainty due to its definition).

Based on the manual classification, we find 49 (almost) mono-periodic, 126 multi-periodic, and 24 hybrid pulsators. The remaining 17 light curves show variability that does not fall into the three classes. Using A_var and P_beat the above introduced quantification for the number of multi-periodic stars stays constant at 124, while we find fewer true mono-periodic stars (33) and more hybrid pulsators (42). Stankov & Handler (2005) classified the stars in their catalogue into similar categories and found a fraction 40% of mono-periodic pulsators. However, they expected this fraction to be lower due to unobserved low-amplitude pulsations in the ground-based light curves, which have much higher noise level than the TESS space photometry. Our fraction of 24% (16.5% when using the classification in Fig. 5) of pulsators with one very dominant oscillation mode is in agreement with this expectation.

The stars of the three sub-classes are found across the whole instability region in the HRD. While most of the hybrid pulsators are close to or inside of the overlap region of the two instability strips, we find some notable exceptions. In particular the distinct group at high T_eff (see above), hosts two hybrid pulsators. A HRD with the stars colour-coded by the classifications is shown in Fig. B.1..

Among the 17 stars we exclude from the β Cep sub-classes, three stars (TIC 22370803, TIC 311328485, and TIC 460977279 (HD 137405)) are likely short-period eclipsing binaries mimicking a pulsation signal in the right frequency domain. We remove these stars in the subsequent analysis. TIC 234644092 might also be of this class but without compelling evidence, we kept the star in our analysis. In contrast, TIC 285370131, TIC 450383484, TIC 329992262 (V1155 Cas), TIC 399436828, and TIC 410100850 show β Cep pulsations but their light curves are dominated by long-term quasi-periodic behaviour. We postulate that they are either Be stars or strongly contaminated. TIC 248275346 (HD 159792) has additional high-frequency components, which create a beating pattern unlike any other star in our sample. The three stars TIC 326815593, TIC 327282466, and TIC 365623418 are noisy and modulated with additional non-periodic variability which stems likely from flux contamination of neighbouring stars. TIC 217602397, TIC 464732844, and TIC 153657168 (HD 114733) exhibit very noisy light curves with hidden periodic variability. Finally, TIC 406696133 is likely contaminated by a nearby, very red (G−J = 4.3 mag) Gaia long-period variable. In the TESS observations this contaminant is brighter than our target star.

The stars TIC 90964721, TIC 395218466, and TIC 1861614014 (HD 339003) are potentially (longer-period) eclipsing binaries with β Cep pulsations, making them interesting targets for detailed modelling. The latter of the three is also included in the sample of Eze & Handler (2024).

4.2. Rotational properties

4.2.1. Surface rotation

Rotation leaves a signature in a star's light curve by bringing chemical spots in and out of view. Some stars in our sample (e.g. TIC 445256790, HD 236535) show this kind of regular sinusoidal modulation in their light curves (see Figures in the online appendix) which is different from the typical multi-periodic pulsational behaviour. Rather, the amplitude of the pulsations is constant but the light curves reveal an added sinusoidal component with a period of a few days. When comparing the typical longer variability to the above derived estimate of the rotation frequency, we find them to match. Hence, the longer-term variability is caused by rotational modulation, a feature commonly observed in high-resolution spectroscopy and high-precision photometry among the brightest β Cep stars (e.g. λ and κ Sco Uytterhoeven et al. 2004, 2005; Handler & Schwarzenberg-Czerny 2013).

In order to select all possible β Cep stars with rotational modulation, we searched for matching signals among the pre-whitened independent frequencies and the Gaia ESP-HS spectroscopic rotation rate. For each star with f_rotsin i, we select the closest signal and accept it as the rotation frequency if it satisfies 0.9<f/f_rot<2.5, allowing for 1.1≲sin i<0.4 to include the unknown projection effect of the rotational axis in the line of sight. Given the large uncertainty on the projected rotation rate the chosen range might exclude some true rotational signals but it ensures a clean selection. Overall, we find some variability within the probed range in 74 out of the 182 (41%) stars with an estimate for f_rotsin i.

Most of the potential photometric rotational signals are close to the estimated f_rotsin i indicating inclinations closer to edge-on. This skewed distribution stems from two selection effects favouring sin i∼1 over sin i∼0. Firstly, the inclination angles of least cancellation of low-degree tesseral and sectoral modes deliver sin i∈[0.8,1.0] (see Appendix 1 in Aerts et al. 2010, for a definition and value of these viewing angles). Moreover, for the fastest rotators the pulsation modes tend to be sectoral modes concentrated to a band around the equator, maximising the visibility of such pulsations from the equatorial region (e.g. Reese et al. 2006). Secondly, stars with low sin i fall below the detection threshold in the Gaia v sin i measurements and hence do not occur in the sub-sample of stars with this quantity available.

A substantial fraction (41%) of our sub-sample features detectable rotational modulation, likely of magnetic origin through chemical spots. This is a higher fraction than the typical ∼10% of magnetic field detections in spectro-polarimetric observations of massive stars (e.g. Grunhut et al. 2017). Hence, the fraction of magnetic β Cep stars is larger than among non-pulsating massive stars (as suggested e.g. by Hubrig et al. 2009), space photometry is a more sensitive detection diagnostic, or the detected signals are not only connected to magnetism.

4.2.2. Influence of internal rotation on pulsations

Due to the shift between a reference frame rotating with the star and the observer's inertial reference frame, the observed pulsation frequencies of non-zonal modes get shifted with respect to their value in a frame of reference corotating with the star. We thus expect f₁ and f_rot to be correlated for sectoral (l=|m|>0) or tesseral (0>|m|>l) modes. Depending on the mode's azimuthal order, a smaller or larger frequency shift of m〈Ω〉 occurs for the mode's frequency between the two reference frames, where 〈Ω〉 is the averaged internal rotation frequency throughout the stellar interior. In general, the internal rotation of stars is not rigid and may hence differ substantially from the surface rotation f_rot. Asteroseismology has shown the level of differentiality between the near-core region, the envelope, and the surface rotation is limited to about 10% for most intermediate-mass main sequence single stars with a convective core (Li et al. 2020b; Aerts 2021). However, the few β Cep stars with such measurement show a large spread in the level of their differential rotation (Burssens et al. 2023).

Figure 6 shows the main frequency observed in the TESS light curves against the estimated projected surface rotation frequency. We find a positive correlation between the two properties, as anticipated. We only have a (rather uncertain) direct estimate of the projected value of the surface rotation frequency in the line of sight, while the pulsation frequency's shift is mainly determined by the rotation in the region of the star where the mode has most of its probing power. The low-order modes in β Cep stars are dominant in the envelope and hence mainly probe the rotation in that region, denoted here as Ω_env. This, and the fact that the shift of any non-radial mode frequency in the line of sight is proportional to mΩ_env, implies that the main pulsation frequency f₁ is not tightly related to the projected surface rotation frequency, explaining the scatter in Fig. 6.

Fig. 6. Main pulsation frequency f1 against the projected surface rotation frequency frotsin i of the 162 stars with available v sin i in the Gaia ESP-HS data. Uncertainties in the pulsation frequency are typically within the symbol size.

5. Mode identification from TESS and Gaia data

Identifying the degree l of at least one pulsation mode is required to reduce degeneracy and perform meaningful asteroseismic modelling of β Cep stars (Hendriks & Aerts 2019). Additionally, an observational distribution of the mode degrees could advance our understanding of the mode selection mechanisms. Hence, we attempt to identify the pulsation modes of the Gaia sample to (1) learn about the distribution of pulsation modes in β Cep stars from a large sample based on space photometry, (2) provide the input for detailed future modelling of the most promising sample stars, (3) open pathways towards ensemble asteroseismology of β Cep pulsators observed across the entire sky from combined TESS and Gaia space data.

5.1. Multi-colour photometry

How a mode's amplitude changes with wavelength depends on its degree, with low-degree modes showing a stronger decrease with wavelength. This has historically been used to identify pulsation modes of β Cep stars using ground-based multi-colour photometry or time-resolved spectroscopy (see Sect. 1). In multi-colour photometry the amplitude ratios are calculated from the amplitudes of the pulsation modes as measured in different photometric filters. The bluest filter is typically used as the reference filter because of its higher S/N. In our work, the highest S/N is provided by the TESS filter. However, we decided to use the bluest filter as reference band to enable comparison with earlier literature work on photometric mode identification (cf. Heynderickx et al. 1994). As already discussed in Hey & Aerts (2024), the Gaia time series photometry offers multi-colour photometry which may help in mode identification.

Figure 7 shows amplitude ratios of the Gaia G_RP and TESS T band observations of stars with the same dominant pulsation frequency in Gaia and TESS. The ratios are computed with respect to the G_BP amplitude. In order to understand the position of different stars in this diagram and assess the feasibility of our approach, we marked stars with unambiguous known mode identifications from the literature (l = 1: TIC 128821888 (12 Lac), TIC 293680998 (IL Vel), TIC 436285033; l = 2: TIC 169786455 (V856 Cen), TIC 288489491 (SY Equ), TIC 327656776 (KZ Mus)). Merely six stars in our sample have an identified mode from long-term or multi-site campaigns in the literature since these observations focused on stars that are either too bright for Gaia or have an amplitude below 4 mmag. Only TIC 419354107 (BW Vul) has an identified dominant radial pulsation in the literature, though its amplitude is so large that non-linear effects occur (Aerts et al. 1995), which may affect its amplitude ratios. Hence no stars with known dominant linear radial mode are marked in Fig. 7.

Fig. 7. Pulsation amplitude ratio AT/ABP against ARP/ABP for the main mode detected in the Gaia photometry of 205 β Cep stars. The colour-coding is according to the amplitude ratio AG/ABP. Crosses mark a few well-studied bright pulsators with secure mode identifications from the literature.

As expected the quadrupole (l = 2) modes fall in the top right corner with amplitude ratios close to unity independent of the filter. Identified dipole (l = 1) modes are found at lower amplitude ratios, while we expect to find the stars with dominant radial (l = 0) mode towards the lower left because their amplitudes decrease strongest with wavelength (see Heynderickx et al. 1994). No clearly separated groups emerge from Fig. 7 in terms of the mode degree. Hence, we cannot achieve unambiguous identification for any given star from Gaia and TESS multi-colour amplitude ratios alone. However, we also use the sample properties and additional information from the light curves. Particularly the trend of having l = 2 in the upper right corner and l = 1 in the centre of the plot does emerge and is encouraging.

5.2. Rotationally split multiplets and standstill stars in the TESS data

Given the high photometric precision of the TESS light curves of these high-amplitude β Cep pulsators, mode identification might also be possible by identifying rotationally split multiplets. Following the prototypical case of HD 129929 (Aerts et al. 2003), we use the pre-whitened pulsation frequencies in the TESS data. A pulsation mode with l>0 can be split into 2l+1 observable mode components of azimuthal order m. In the corotating frame, the offset from the m = 0 central mode is proportional to 〈Ω〉 when treating the Coriolis force perturbatively up to first order (Ledoux 1951).

We exploited the connection between the two measurements f₁ and f_rotsin i and the unknown 〈Ω〉 required to find rotational multiplets, by relying on the correlation seen in Fig. 6. Up to first order in 〈Ω〉 the splittings are symmetrical, which is a good approximation for slow to modest rotators among the β Cep pulsators. This enabled us to distinguish between chance alignments of frequency peaks and true multiplets in the periodograms.

Figure 8 shows two examples of rotationally split multiplets we found in the TESS data. For the example in the top panel, the dominant frequency in the Gaia DR3 light curve matches a significant frequency in the TESS light curve, where it is not the dominant one due to strong beating effects. In total, we identified 19 stars with a rotational multiplet involving the dominant Gaia frequency in the TESS data satisfying at least four of the following five criteria: (1) the multiplet is complete, (2) the splitting is on the same order of magnitude as f_rotsin i, (3) the splittings are symmetric up to |δ_A|<0.1 (see Sect. 6.2), (4) the multiplet and its splittings do not conflict with any other potential multiplet in the Fourier spectrum, and (5) none of the frequencies in the multiplets originated from contaminating sources. The majority of these stars (17, representing 89%) pulsate in a dominant l = 1 mode and only two in a l = 2 mode (representing 11%). As we only seek three evenly spaced frequencies to identify an l = 1 mode compared to five for an l = 2 mode, the preponderance of l = 1 modes may suffer from selection bias. In addition, we identified 14 stars with potential split multiplets that almost satisfied four conditions. We come back to them after the mode identification in Sect. 5.3.

Fig. 8. Examples of rotationally split multiplets found in the amplitude spectra of the TESS data. The top panel shows an l = 1 mode split into three azimuthal components. The bottom panel shows one of the few identified complete l = 2 multiplets. The splittings are marked with grey markers connected by a line. The red markers at the top of each panel indicate the extracted frequencies from pre-whitening. The dashed lines mark the dominant frequency in the Gaia DR3 light curve.

Unlike split multiplets for pulsators with l≥1, radial modes cannot easily be identified from the periodograms. However, the light curve of a high-amplitude radial pulsation in β Cep stars sometimes exhibits the standstill phenomenon⁴, just as for Cepheids (e.g. Bersier et al. 1994). During the rise of the light curve, at the time of maximum radial velocity, the brightness briefly stalls before increasing further and dropping rapidly after the maximum (Walker 1954). The standstill phenomenon stems from the development of shock fronts in the expanding photosphere (Mathias et al. 1998; Fokin et al. 2004). The most prominent example among the known β Cep stars revealing a standstill in its photometric light and spectroscopic radial-velocity curves is the high-amplitude non-linear radial pulsator BW Vul (TIC 419354107, e.g. Sterken et al. 1987; Aerts et al. 1995). Based on phase-folded and binned TESS light curves, we identified five more stars with standstills (TIC 199208037, TIC 274070426, TIC 347520615, TIC 395218466, TIC 444357430). The folded TESS light curve of TIC 274070426 is shown alongside the one of BW Vul in Fig. 10, while those of the other four are displayed in the online appendix. To our knowledge, none of these stars have been reported before to show this phenomenon. TIC 395218466 is of particular interest as it is a member of an eccentric eclipsing binary as detailed in the online appendix.

Three additional stars have been identified to have a radial mode due to the position of its frequency with respect to those of split multiplets. In total, we thus identified nine radial mode pulsators. Figure 9 shows the amplitude ratio distributions of all the mode degrees identified from rotational splitting or a standstill in the TESS data, as well as from securely identified modes from the literature. We show the ratios as a function of the effective wavelength of the considered filters⁵. For each photometric band and pulsation mode, we show the mean value of the amplitude ratio and a violin plot of the distribution. For all three considered degrees l the distributions overlap strongly, which explains why we could not find a clean separation in Fig. 7 from the Gaia data alone and had to work with extra information from the rotational splitting or standstill in the TESS data. We find that the mean amplitude ratio values are distributed as expected, with a larger amplitude ratio for higher l. The radial modes show clearly smaller amplitude ratios than the non-radial modes. The distinction between l = 1 and l = 2 modes is small, following the results in Fig. 7.

Fig. 9. Violin plots of amplitude ratios for identified pulsation modes over the filter's effective wavelengths. The lines connect the mean values for each mode degree. The filter names are given at the top. The GRP and T filters are displaced by ±30 nm to avoid overlapping distributions.

5.3. Probabilistic mode identification from TESS and Gaia

With the secured mode identifications from rotational splittings and standstills in the TESS data, we now return to the Gaia-TESS multi-colour photometry for all pulsators, to identify the dominant mode of the stars without rotational splitting. Given the overlap of the distributions of the amplitude ratios for modes of l = 0,1,2 in Fig. 9 and the spread in Fig. 7, we aimed to estimate the mode identifications for the full sample based on a probabilistic framework including all photometric measurements and their uncertainties. To find the probability of a given mode identification under the measured Gaia and TESS amplitude ratios, we used the distributions shown in Fig. 9 as our training set and modelled them with a 3D Gaussian kernel density estimator (KDE, K_l). As only very few dominant pulsation modes have been identified with a degree l>2 from ground-based observations (e.g. Heynderickx et al. 1994), we assumed l∈{0,1,2}. By cross-validation, we ensured that a Gaussian KDE indeed describes our data best. Similarly, we used cross-validation to choose the optimal value for KDE bandwidth parameter.

Fig. 10. Light curves folded over the dominant frequency of two stars with clear standstills in our sample (stalling of the light curve near phase 0.5). The folded light curve is binned into 35 bins, shown as black circles, to which we fitted a sinusoidal with a fixed frequency for reference (red). The top panel shows the archetypical standstill star BW Vul (Walker 1954).

Since our training set is sparse, we augmented our observations by bootstrapping. In order to account for correlations in the data, we resampled the underlying pulsation amplitudes using a Gaussian perturbation based on the measurement uncertainties. We then calculated the amplitude ratios for each augmented data point. For every securely identified mode, we added 1000 points to capture the uncertainty distribution properly. Based on the KDE calculated over the augmented data set and for every star, we calculated the likelihood L_l for a given mode l in a Naïve Bayesian scheme (Zhang 2004) by integrating K_l weighted by a Gaussian distribution N with the uncertainties as its covariance σ:

$$L_l=\int \int \int K_l(\text{A})N(\text{A},\mathrm{\sigma })\mathrm{d}\text{A}.$$(3)

To estimate the probability for the mode identification, we initially chose a flat prior but noticed that many stars with clear l = 1 splittings were assigned either the correct degree but with a low probability or the wrong degree. This is an effect of the overlap of the distributions which reduces the probability of l = 1 modes. Given the overall clear separation of l = 1 modes, we chose to use a non-flat prior. The best prior is the one that enables the correct recovery for all known mode identifications. From a grid search of possible priors, we found the weights 0.33, 0.47, 0.20, for radial, dipole, and quadrupole mode occurrences, respectively, to provide the best result.

With the weighted probabilities, we identified 143 mode degrees (l = 0: 23, l = 1: 104, l = 2: 16) with probability p>60%. The large number of l = 1 modes follows from the abundance of triplets in the training set and their prior probability. Five stars from the training set do not reach the chosen threshold of p>0.6 but only one of them is assigned the wrong degree. We have thus defined a robust probabilistic mode identification methodology suitable to treat all sample stars. For reference and to inform any future independent modelling efforts of the stars by community members, we provide all mode identifications (including those with p<0.6) and their probabilities in Table 1.

Some stars in Fig. 7 exhibit a low A_T/A_BP amplitude, while near unity in A_RP/A_BP. A detailed analysis revealed that the majority of these stars contain a fraction of contaminating flux in their light curves. Although we tried to minimise the impact of contaminating light, it cannot always be avoided, for example in cases of increased sky background in TESS observations or very nearby bright stars as found in open clusters. Thanks to the two Gaia amplitude ratios included in the degree determination, the impact of contaminating light lowering the TESS amplitude is limited. However, it highlights the vulnerability of TESS amplitudes due to contamination from nearby sources and the sky background.

In summary, we are able to identify the dominant mode degree of 143 Gaia β Cep pulsators with a probability p>60% based on their TESS light curve, Gaia multi-colour amplitude ratios, and the estimated surface rotation rates. For five additional stars, the multicolour photometry alone was inconclusive; even so, we identified their mode degrees from splittings or a standstill. This sample of 148 stars represents the largest sample of identified pulsation modes in β Cep stars to date. Even if each of these pulsators revealed a limited number of identified pulsation modes compared to the handful of β Cep stars modelled in great detail (as summarised by Burssens et al. 2023), our results offer a prime sample for further and more detailed asteroseismic population studies of β Cep stars. Only this type of sample studies can offer suitable input to calibrate stellar evolution theories for galactic massive field stars as a whole. In the following, we present our initial global seismic analysis and grid modelling of this unique and homogeneously composed sample.

6. β Cep ensemble asteroseismology

To gain deeper insights into the population of the Gaia-observed β Cep stars, we used the probabilistic identification of the dominant modes. In the following, we characterise the ensemble and subject all sample stars to one homogenous type of basic modelling by relying on a grid of stellar evolution models. The aim of this work is to understand the properties of the population of field β Cep stars instead of doing detailed studies of only a few individual stars. Our approach is a more powerful way to understand global aspects of stellar evolution theory and pinpoint some of its possible shortcomings for the entire mass regime covered by the sample.

6.1. Distributions of the identified modes

The mean behaviour of the identified modes of the full sample over the wavelength is very similar to the hand-selected sample in Fig. 9. However, we find a slightly clearer separation between the two considered non-radial modes. Notably the distribution of A_T/A_BP for l = 1 shown in Fig. B.3. has a heavy tail in the TESS amplitude but not in the other photometric bands.

Figure 11 shows an amplitude ratio diagram similar to Fig. 7 colour-coded by the determined mode degree. The distribution of mode degrees follows the expected distributions with l = 2 modes in the top right (orange) at the highest amplitude ratios. Stars with l = 1 as their primary mode degree are located in the central part (blue) and radial modes l = 0 in the lower left (green).

Fig. 11. Amplitude ratio diagram based on ARP/ABP and AT/ABP, similar to Fig. 7. The colour-coding shows the mode degree as determined with our probabilistic framework: l = 0 (green), l = 1 (blue), l = 2 (orange). The grey symbols denote stars with p<0.6. The symbol size is proportional to the probability of the mode identification. The stars marked with black outlines constitute the underlying sample from Fig. 9.

When going into details, we find a complicated pattern with several overlapping regions. Due to the three dimensionality of the distributions these overlap regions create fuzzy borders in the two-dimensional projection in Fig. 11. Stars to which our algorithm assigned low probabilities in the mode identification belong to the overlapping regions between the different degrees as also seen from the shown contours. An interesting feature is the group of high-probability l = 1 modes at the bottom of Fig. 11. This region creates the tail in the A_T/A_BP amplitude ratios in Fig. B.3.. From the TESS data, we find that some of these stars have multiplets in which the Gaia observed frequency (and often also the main TESS frequency) is not the central frequency. The multi-colour photometry also carries information on the intrinsic mode excitation amplitudes, aside from the mode geometry and inclination angle of the star. Our findings illustrate that the mode amplitudes are not equal among the prograde and retrograde components of the same |m| within a multiplet for these heat-driven modes. This is in contrast to the case for stochastically excited modes in slowly rotating low-mass stars, where the amplitudes allow for the derivation of the star's inclination angle. Mode visibilities purely based on geometrical arguments are not the only factor that determines the amplitudes (and hence detectability) within a multiplet for β Cep stars, as has been known since long from ground-based photometry.

When placing the stars with identified modes in the HRD, we find that the β Cep stars with a dominant radial mode are mainly of lower mass (≲15 M_⊙). The radial pulsators seem to concentrate in the later half of the main sequence in the overlap region between the β Cep and SPB instability strips. In contrast, non-radial pulsators can be observed throughout the whole β Cep mass range (cf. Fig. B.2.), although two stars with radial modes stand out due to their relatively high luminosity and mass.

6.2. Rotationally split multiplets and asymmetries

As introduced in Sect. 5.2, rotation splits non-radial modes into multiplets with splittings proportional to the averaged envelope rotation rate. We now return to the 14 stars with likely multiplets (not used in the definition set to construct the KDE in Sect. 5.3). We now have their probabilistic mode identification. Adding them to the 19 stars with deterministic identification from the full multiplet structure we had already from Sect. 5.2, they compose a sample of 33 stars with multiplets which we now use to study the stellar rotation properties. The full sample of these 33 stars with multiplets is shown in the online appendix.

In Fig. 12, we show the observed splittings against the estimated (projected) surface rotation rate. As expected, a positive correlation between the rotation rate and the splitting is observed. However, we recall that β Cep stars can have strong radial differential rotation (Aerts et al. 2003; Pamyatnykh et al. 2004; Burssens et al. 2023) and the surface rotation rate may differ from the internal rotation rate at the fractional radius where the mode is most sensitive. Stars at v sin i = 0 km s⁻¹ have undetectable line broadening in Gaia; however, we include them here as they are included as such in Gaia DR3.

Fig. 12. Observed mode splittings against the projected surface rotation rate (frotsin i, encircled in thin grey) or photometric surface rotation rate (frot, encircled in thick solid black). The three lines show the splitting values when equal to the surface rotation rate (solid) as well as one-half (dotted) and twice (dashed) the surface rotation rate. The uncertainties on the splittings indicate the maximum difference of the individual splittings to the mean splitting of the multiplet and are a measure for the asymmetry. The colours show our deterministic identifications as l = 1 (blue) and l = 2 (orange) modes.

We also mark the absolute asymmetries (maximum difference to the mean splitting) in Fig. 12. For most of the stars this asymmetry is small, though it is pronounced for a minority. A purely rotationally split multiplet would consist of evenly and symmetrically spaced components if the Coriolis force due to the mean envelope rotation has only a small influence. This allows for a first-order perturbative approach to compute the oscillation mode frequencies. However, second-order effects are relevant for β Cep stars and include rotational deformation of the star (Saio 1981). This effect shifts the split frequencies asymmetrically with respect to the central frequency (Dziembowski & Goode 1992; Goupil et al. 2000). In addition, magnetic fields (e.g. Bugnet et al. 2021) or binarity (e.g. Guo 2021) also influence the split multiplets, leading to asymmetries in the observed frequency spectra and enabling detailed asteroseismic inferences (Guo et al. 2024).

Guo et al. (2024) (see also Goupil et al. 2000; Deheuvels et al. 2017) define the observed asymmetry as

$${\delta }_A=-\frac{f_{-m}+f_{+m}-2f_0}{f_{+m}-f_{-m}},$$(4)

with f₀ the central frequency and f_±m the frequencies of azimuthal order m=±1. In Fig. 13, we show the 33 most obvious asymmetries occurring in the multiplets detected in our sample stars against the central frequency f₀ of the multiplet. Most of these splittings are positive as expected for rotational deformation (Guo et al. 2024). Negative asymmetries in β Cep stars can hint towards internal magnetic fields (Mathis & Bugnet 2023), especially for dipole p-modes. Our understanding of core magnetism in massive stars is still limited (Lecoanet et al. 2022; Vanlaer et al., in prep.), despite its potential impact on stellar evolution. A large sample of measured internal magnetic fields from asteroseismology would deliver constraints for evolutionary models, similarly to the recent advances in red giants (Li et al. 2023; Deheuvels et al. 2023).

Fig. 13. Asymmetries of 33 rotationally split multiplets δA against their central frequency f0 colour-coded by the projected surface rotation rate. The square symbols mark stars with measured photometric surface rotation rates (frot). The symbols connected by a line belong to the same (l = 2) multiplet. The black outlines indicate stars where the Gaia frequencies do not match the TESS frequencies within our threshold. The measurement uncertainties are within the symbol sizes.

6.3. Modelling the population of β Cep stars

To understand the population of galactic β Cep stars better and gain insights into the distributions of their physical properties with respect to the selected main pulsation mode, we applied grid modelling. In doing so, we used extra seismic information compared to classical grid modelling where one just relies on evolution tracks and the measured effective temperature and luminosity to determine stellar properties. Here, we used the β Cep MESA/GYRE stellar evolution and pulsation grids presented by Burssens et al. (2023). The evolution grid was calculated for masses between 9 M_⊙ and 21.5 M_⊙ at solar metallicity with envelope mixing by internal gravity waves set to D = 1000 cm² s⁻¹ at the base of the envelope and increasing outwards as ρ⁻¹ and an exponentially decaying diffusive core overshoot with the parameter range 0.005≤f_ov≤0.035. For details of the physical implementation, we refer to Burssens et al. (2023).

The required mode identification for the modelling is not limited to the degree l but we also need the radial order n and azimuthal order m. All models will be fit in the azimuthal order m = 0 necessitating us to identify the central peak of split modes denoted here as f₀. For complete multiplets or radial pulsators, f₀ is measured. For the few quadrupole pulsators with incomplete multiplets, we assumed the highest peak of the central two peaks to be the one with m = 0.

The radial order can neither be deduced from the observations nor from a direct fit of the model. Rather, we fitted models with the radial orders n∈{0,1,2} for radial pulsators and n∈{−2,−1,0,1,2} for non-radial modes⁶ using a χ² minimisation over log T_eff, log L, and f₀. In principle, the radial order can be larger but higher-order modes have not been securely established in β Cep stars, hence we did not probe them.

Stars without a model within the 1σ uncertainty ellipse were excluded from the subsequent analysis, reducing the number of modelled stars to 119. For each star, we chose the best fitting model for each n and compared the position of its other modes to the observed frequency spectra. Our focus laid on the relative positions of the main pulsation mode and other visible modes with respect to those predicted by the model. In most but not all cases our choice of n coincides with the lowest χ² value among the set of best fitting models. We show one illustrative example for the identification in Appendix C, including a detailed reasoning for our choice of n.

We find the majority of the radial pulsators to pulsate in their fundamental mode (eight out of twelve modelled radial pulsators). Three radial pulsators pulsate in their first overtone and one stars is likely a second overtone pulsator. A similar picture emerges for the 107 non-radial pulsators. The dominant pulsation mode in 50 stars is the p₁-mode (n = 1) and 24 pulsate in g₁ (n=−1). Only 29 stars pulsate in their second-order mode with 28 identified p₂ and one g₂ pulsators. The remaining four stars were identified to be dominant f-mode pulsators.

For the final model selection and uncertainty estimation, we employed a MCMC approach with a χ² merit function, similar to Fritzewski et al. (2024), to which we refer for details. Since we rely on a grid of stellar models with fixed input physics on some aspects (notably the level of envelope mixing, where the low-order modes have dominant probing power), we allow for an overall mismatch between observed and fitted frequencies (see Aerts et al. 2018, for an extensive discussion). In order to encapsulate this theoretical uncertainty caused by the fixed yet uncertain input physics adopted in the model grid, we follow the predictions for zonal p-modes by Aerts et al. (2018) in the β Cep mass regime and allow for a typical frequency uncertainty of 0.1 d⁻¹. This captures any systematic errors in the predicted frequencies while allowing to recognise the correspondence between the observed and predicted frequency pattern (cf.  the procedure explained in Appendix C).

The histograms in Fig. 14 show distributions of the obtained masses, normalised core hydrogen fraction ($X_c^{\prime }$), age, relative convective core mass, radius, and horizontal versus radial displacement of the dominant mode. As expected from the distribution of β Cep stars in the HRD in Fig. 3, the number of β Cep stars declines with mass, with most β Cep stars found in the range 10≤M_ini/M_⊙≤15. Eight stars could not be modelled because they lie above the upper mass limit of 21.5 M_⊙, as shown by the evolutionary tracks in Fig. 3. In agreement with Fig. B.2., radial pulsators are mostly found among the lower-mass β Cep stars.

Fig. 14. Distribution of β Cep parameters for the modelled stars. Each panel contains the histograms for the full sample (grey) and each degree l = 0 (green), l = 1 (blue), and l = 2 (orange). From top to bottom, these panels show the initial mass distribution, the normalised core hydrogen fraction $X_c^{\prime }$, the stellar age in years, the relative convective core mass, the radius, and the K-value.

The distribution of $X_c^{\prime }$ shows more structure when focusing on the pulsation degree. Overall, the majority of our sample is in their second half of the main sequence evolution. In particular, the small group of radial pulsators are found in the second half of their main sequence evolution, while non-radial pulsators are found over the whole range.

We do not show a distribution of the core overshoot because most of the stars are modelled with f_ov≈0.02, which corresponds to the central value in the grid. However, the uncertainties span the whole range of f_ov covered by the grid, which can be explained by the limited sensitivity of the observed modes to the near-core region and the degeneracies between the mass, metallicity, and core overshoot when relying on just one mode (e.g. Ausseloos et al. 2004). The estimated core overshoot values result in a value for the mass-normalised convective mass (M_core/M_★). This inferred parameter is the more interesting quantity from a stellar evolution point of view and it is a function of both the luminosity and age. The majority of stars has a value in the range 0.2≤M_core/M_★≤0.3, in agreement with the distribution of initial mass shown in the top panel and the HRD. Notably the radial pulsators are concentrated around M_core/M_★≈0.25 due to their common evolutionary status. The few stars with a relatively high convective core mass belong to the higher-mass young O-stars well separated in the HRD (cf. Fig. 3). The radius distribution is peaked at intermediate radii, especially for dominant radial mode pulsators.

As a final distribution, we show in the bottom panel the distribution for the K-value, defined as K≡GM/(2πf₀)²R³, for the non-radial modes. This dimensionless quantity is connected with the boundary conditions for the pulsation equations (Aerts et al. 2010, Chapter 3). In the case of non-radial modes, it represents the ratio between the horizontal and radial displacement of the fluid elements at the stellar surface caused by the mode (as this ratio is zero for radial modes, these are not included in the bottom panel of Fig. 14). This parameter has high values above unity (typically between 10 and 100) for high-order g-modes of observed SPB stars (De Cat & Aerts 2002), as their modes dominantly cause horizontal displacements. Here, as expected for low-order modes, we find the opposite with dominantly radial displacements but with non-negligible contributions from horizontal motions. The values for K range from about 0.03 to about 0.2 for the dominant mode of most of the Gaia β Cep stars in our sample.

To unravel the correlation between the convective core mass and the core hydrogen fraction, we show these two parameters in Fig. 15. The overall distribution shows clearly the distinction between stars with M_core/M_★≤0.3 and these with more massive convective cores. The stars with higher relative core masses are typically more massive and younger as expected from stellar evolution theory. As seen from the evolutionary tracks of different masses, only the lowest mass β Cep stars can be found among full breadth of the main sequence evolution, a consequence of the position of the instability strip lying closer to the zero-age main sequence for these stars.

Fig. 15. Normalised convective core mass Mcore/Mini over the normalised core hydrogen mass fraction $X_c^{\prime }$. The stars are colour-coded by their primary degree: l = 0 (green), l = 1 (blue), and l = 2 (orange). The shape represents the radial order: g2-modes are represented by thin diamonds, g1 by thin triangles, radial fundamental modes by triangles, p1 by circles, and p2 by squares. The grey tracks in the background show the evolution for the indicated masses with fov = 0.005 (solid, labelled at low $X_c^{\prime }$) and fov = 0.035 (dashed, labelled at high $X_c^{\prime }$). The typical uncertainty is indicated in the upper right corner.

In summary, we successfully modelled 119 β Cep stars with a grid modelling approach adding the seismic information on the dominant mode to the two classical observables, effective temperature and luminosity. We provide distributions of their masses, ages, convective core masses, and radii. As such, these distributions provide a first insight into the properties of the population of β Cep stars. Of particular interest for future investigations on the mode excitation of this population is the connection between the stellar age and the degree of the dominant mode, if any. We give the best fitting model properties in Table 2. A discussion of the accuracy of the asteroseismic parameters in comparison with the input values from Gaia can be found in Appendix D.1.

6.4. Period-luminosity relations

Radial and low-order non-radial p-modes probe the internal density of a star. For a given type of pulsator this translates into a connection between a mode's frequency and the luminosity as established for Cepheids by Leavitt & Pickering (1912). For the δ Sct stars, which are situated in the low-luminosity part of the classical instability strip, this also gives rise to a period-luminosity relation (Gaia Collaboration 2023) but it is not so clean as the one found for Cepheids because a fraction of these p-mode pulsators have a dominant non-radial mode (Ziaali et al. 2019).

While the δ Sct stars occupy a small region in mass and radius, this is not the case for the β Cep stars. Given their large spread in mass and radius along the main sequence, one does not expect a clear period-luminosity relation for this class as a whole. Nevertheless, several authors attempted to find such a relation in the past (e.g. Blaauw & Savedoff 1953; Jakate 1980; Waelkens 1981). However, the uncertain mode identifications led to a large scatter in these relations (see review by Lesh & Aizenman 1978).

Here we show the relations and their scatter occurring for our sample of β Cep star when using only those with identified pulsation modes. Figure 16 shows the period-luminosity relations in the log P–log L/L_⊙ plane, colour-coded by the mode degree. Using a linear regression, we arrived at a period luminosity relation for each sub-set of pulsators according to their dominant mode

Fig. 16. Period-luminosity relations for stars with modes of different degrees: l = 0 (green), l = 1 (blue), and l = 2 (orange). The period-luminosity relations include only the radial fundamental mode (n = 0) for radial pulsators and p-modes for the non-radial pulsators. The black line shows the relation of the Shi et al. (2024) sample. For comparison, a similar period-luminosity relation for our sample is plotted in purple (dashed). The green shaded area shows the uncertainty band of the fundamental radial-mode pulsators. The small black data points in the background show dominant g-mode and higher-order radial pulsators excluded from the fits.

$$logL/L_{\odot }=(2.21\pm 0.52)logP+(5.8\pm 0.36)(l=0)$$(5)

$$logL/L_{\odot }=(1.81\pm 0.36)logP+(5.8\pm 0.27)(l=1)$$(6)

$$logL/L_{\odot }=(0.94\pm 1.34)logP+(5.1\pm 1.02)(l=2)$$(7)

with P the mode period in days. Figure 16 shows only the 1σ-uncertainty band for the stars with a dominant fundamental radial mode. Those with a dominant dipole p-mode have a similar uncertainty band while it is larger for the quadrupole mode pulsators.

We find the three relations to be stacked by their mode degree, such that the lowest degree modes (radial) follow a period-luminosity relation with the lowest luminosity for a given period. The relation for the quadrupole modes intersects with the relation for the dipolar modes due to fewer stars in the sample and larger spread in the radial order. Shi et al. (2024) recently provided a similar period-luminosity relation for their sample of β Cep stars. In comparison to a linear regression of all 155 stars from Shi et al. (2024), regardless of their mode degree, our period-luminosity relations are steeper. When including all stars of our sample, we find a similar slope to theirs, highlighting the importance of mode identification.

6.5. Average envelope rotation

Radial differential rotation was observed in the handful of β Cep stars with a measurement of the internal rotation (e.g. Aerts et al. 2003; Pamyatnykh et al. 2004; Dziembowski & Pamyatnykh 2008; Burssens et al. 2023; Vanlaer et al., in prep.). Based on the observed splittings (Sect. 5.2) and our modelling results for the Gaia β Cep stars, we calculated the internal envelope rotation rate Ω_env from the Ledoux constant (Ledoux 1951) for the best grid model for the star. This rotation rate then also allows us to calculate the level of differential rotation between the envelope and the surface for stars with a measured f_rotsin i, as well as the spin parameter of the mode.

In Fig. 17, we show the estimated Ω_env and derived parameters. Unlike the observed projected surface rotation and as expected, the envelope rotation deduced from the modelling does not have any correlation with the pulsation frequency of the zonal mode, as f₀ is unaffected by Ω_env when treating the rotation perturbatively and only up to first order as we do here (top panel). The highest envelope rotation rate is not found for the fastest surface rotators but for stars rotating at an intermediate velocity (second panel), reflecting the effect of the unknown factor sin i and/or differential rotation between the envelope and the surface. This is visually represented in the third panel of Fig. 17, which shows the level of differential rotation assuming an almost equator-on view (sin i≃1). This lower limit of the true surface rotation frequency is denoted here as Ω_surf. It is used to estimate an upper limit for the differential rotation between envelope and surface, which is plotted against the evolutionary stage $X_c^{\prime }$. Aside from the projection effect, we find the most differential rotators to have an internal rate that is below four times the surface rotation. This is completely in line with the first such measurement ever done for a β Cep star (Aerts et al. 2003), which relied on rotational splitting of both a dipole and quadrupole mode. However, the majority of stars reveal Ω_env/Ω_surf≲2. Hence, despite the few cases of strong radial differential rotation as summarised by Burssens et al. (2023), the majority of β Cep stars have a level of envelope-to-surface rotation between 0.3 and 2.

Fig. 17. Envelope rotation Ωenv and derived parameters against observed frequencies and estimated core hydrogen fractions for 25 stars. The top panel shows the envelope rotation rate against the zonal mode frequency of the multiplet (f0) and the second panel against the projected surface rotation rate frotsin i for the stars having this quantity available. The third and fourth panels display derived rotational properties against the evolutionary stage. In the former, we show an upper limit for the level of differential rotation (Ωenv/Ωsurf) and in the bottom panel the spin parameter (s) of the mode in the envelope. The colour-coding corresponds to the mode degree: l = 1 in blue and l = 2 in orange.

We also present the spin-parameter, s, of the dominant mode as a function of $X_c^{\prime }$. It is defined here as s≡2Ω_env/f₀ and characterises the regime in which the pulsation mode occurs (Aerts & Tkachenko 2024). All the Gaia β Cep stars have their dominant mode in the super-inertial regime with s<1. Noticeably, the star with the largest spin parameter (TIC 282881728) also has the highest level of differential rotation and is in the later half of its main sequence phase.

Finally, we show in Fig. 18 histograms of the dimensionless quantities ϵ²≡(2πΩ_env)²R³/GM and K⁻¹=(2πf₀)²R³/GM, following Goupil et al. (2000). The latter quantity is the inverse of K already introduced above. It represents an assessment of the importance of the centrifugal force with respect to the 2nd-order Coriolis force in the treatment of the pulsations, while the former evaluates the effect of the centrifugal force with respect to gravity. Goupil et al. (2000) reported the centrifugal force (and hence the stellar oblateness) to be more important for the computation of the p-modes of rapidly rotating δ Sct stars compared to the 2nd-order rotational effect caused by the Coriolis force. While the mode frequencies of β Cep stars are in general lower than those of δ Sct stars, we still find the second-order Coriolis force to be less important than the centrifugal force. However, the difference in importance between the two forces is smaller than for δ Sct stars.

Fig. 18. Histograms of ϵ2 (top) giving the influence of the centrifugal force with respect to gravity and the inverse of the K-value (bottom).

7. Conclusions

We analysed the TESS light curves of 216 Gaia-detected β Cep candidates (Gaia Collaboration 2023), including 145 new detections. By placing these stars in a Gaia HRD, we confirmed that all of the stars are late O or early B stars near the β Cep instability region. However, a detailed light curve analysis shows that at least three stars are short-period eclipsing binaries. By categorising the TESS light curves, we find the majority of β Cep stars in our sample to be multi-periodic.

The combination of the Gaia and TESS time series photometry, allowed us to determine the mode degrees for 148 stars in our sample from their amplitude ratios. Most observed β Cep stars pulsate in non-radial modes with l = 1. Our identified modes constitute the largest set of identified modes of β Cep stars available to date.

From the TESS light curves, we found 34 rotationally split multiplets accompanying the main pulsation mode. The majority of these multiplets have positive asymmetries in their splittings. However, several stars show negative asymmetries and are prime targets for detailed follow-up analyses. We provide the first distributions of the envelope rotation rates and of upper limits for the level of differential rotation for a meaningful sample of β Cep stars. We find that most stars have envelope-to-surface ratios lower than two and many of them rotate almost rigidly. We do point out that our sample only consists of β Cep stars with a dominant mode in the super-inertial regime, and thus lacks the representation of the fastest rotating class members due to our sample selection from Gaia.

Further, we performed the first ensemble modelling of 119 β Cep stars based on the grid calculated by Burssens et al. (2023). We provide distributions in age, convective core mass, and mass in order to inform population studies of β Cep stars. With the large set of identified modes, we are in a position to provide these distributions for each mode degree. As expected from instability calculations, we find the stars to have completed at least one-third of the main sequence evolution. The radial mode pulsators in particular occur mostly in the intermediate main sequence stage. The period-luminosity relations for the different mode degrees are steeper than that obtained for the whole ensemble.

With the upcoming Gaia data releases and the ongoing TESS mission our approach will be applicable to the majority of known and still-to-be-discovered β Cep stars. With a longer time base and improved instrumental corrections, Gaia will deliver many more new β Cep stars with future data releases. Even with the current data, expanding our work to all known β Cep stars by searching for lower-amplitude pulsations in the Gaia data is possible. Identifying the majority of main modes in β Cep stars will further improve our understanding of β Cep stars from a population point of view while delivering the highly needed input to ensemble forward modelling. Complementary to our photometric analysis of fainter β Cep stars, the upcoming CubeSpec mission (Raskin et al. 2018; Bowman et al. 2022) will target bright β Cep stars with high-resolution time series spectroscopy from space to determine their pulsation modes, while the Global Asteroseismology Project is targeting mode identification from long-term multi-colour continuous ground-based observations (Shitrit & Arcavi 2024). These ongoing and future projects will revive studies of massive pulsators and unlock their full asteroseismic potential.

Data availability

Full Tables 1 and 2 are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/698/A253. Additional figures containing all light curves and frequency spectra, standstill stars, and splittings can be found on Zenodo via https://zenodo.org/records/15094728

Acknowledgments

We thank the referee for the detailed and critical review that helped to improve the paper. We thank Tim Bedding, Dominic Bowman, Siemen Burssens, Luc IJspeert, Kalarickal Sreenivas, and Vincent Vanlaer for useful discussions. The research leading to these results has received funding from the KU Leuven Research Council (grant C16/18/005: PARADISE). CA and MV acknowledge financial support from the European Research Council (ERC) under the Horizon Europe programme (Synergy Grant agreement N°101071505: 4D-STAR). While partially funded by the European Union, views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. JDR and CA acknowledge the Belgian Federal Science Policy Office (BELSPO) for the provision of financial support in the framework of the PRODEX Programme of the European Space Agency (ESA) This research has made use of NASA's Astrophysics Data System Bibliographic Services and of the SIMBAD database and the VizieR catalogue access tool, operated at CDS, Strasbourg, France. This paper includes data collected by the TESS mission. Funding for the TESS mission is provided by the NASA's Science Mission Directorate. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Software: This work made use of Topcat (Taylor 2005). This research made use of Lightkurve, a Python package for Kepler and TESS data analysis (Lightkurve Collaboration 2018). This research made use of the following Python packages: corner (Foreman-Mackey 2016); emcee (Foreman-Mackey et al. 2013); IPython (Pérez & Granger 2007); MatPlotLib (Hunter 2007); NumPy (Harris et al. 2020); Pandas (McKinney 2010); SciPy (Virtanen et al. 2020); scikit-learn (Pedregosa et al. 2011).

References

Appendix A: Detrending and special photometry

In order to use the pulsation modes of the β Cep stars in our sample we ensured 1) minimal contribution by contaminating stars (Sect. 2.1) and 2) a clean light curve containing as few instrumental signals as possible. We removed instrumental and long-term trends after excluding intervals which proved problematic in multiple stars. We preferred excluding these intervals over complicated detrending steps to avoid overcorrection. Further, some weak long-term trends sometimes remained and the rising or falling curves at the sector edges were not always easily corrected. Hence, we detrended blocks of data separated by a gap of at least three hours with second-order polynomials. In the example shown in Fig. A.1., more than ten stars in our sample display similar instrumental effects at the selected times. We note that in the top panel, a downward slope is present in the first half of the sector which has been successfully removed by our detrending. Afterwards, we σ-clipped the light curves with a threshold of 5 σ. Finally, in 18 light curves, hand-picked time intervals or even entire sectors were removed if they displayed unusual behaviour likely caused by instrumental effects or contamination and σ-clipped them again.

Fig. A.1. Example light curve before (top) and after (bottom) the final detrending. Sections highlighted in red in the top panel have been removed.

For 15 stars, we used single-pixel photometry, in order to avoid contamination from nearby sources. One example is given in Fig. A.2.. Although the main β Cep pulsation is clearly visible in the periodograms of nearby pixels, several additional peaks emerge for these stars too, indicating a periodic contamination from other sources. The particular case of TIC 264026227 is a member of the open cluster NGC 5606, which host other massive stars in close proximity of TIC 264026227. We note that the Gaia detected frequency does correspond to the main TESS frequency, although the light curves in nearby pixels have other periodic components that are more prominent. Hence, using a single pixel for the photometry of this star ensures the cleanest possible extraction.

Fig. A.2. Example of a single pixel TESS light curve. The target TIC 264026227, member of the open cluster NGC 5606, is contaminated by several nearby stars. Top left: The TESS FFI cut-out of the field around TIC 264026227, with the pixel used for or final light curve marked in red. The colour scale shows the flux intensity with yellow being brighter pixels. Overplotted in each pixel is its light curve. Top right: Same as left but with the periodogram overplotted. Bottom: Final light curve shows clear β Cep variability after additional cleaning and detrending beyond what is shown in the top left panel.

Appendix B: Additional Hertzsprung-Russell diagrams and mode identification

For convenience, we provide some extra HRDs with colour-coding by the light curve classification (Fig. B.1., Sect. 4.1) and mode degree as identified from multi-colour photometry (Fig. B.2., Sect. 5.3). Both diagrams show no obvious structure in the light of the current study.

Fig. B.1. Hertzsprung-Russell diagram similar to Fig. 3 but including only the 216 stars from our Gaia-detected sample having TESS light curves. Each star is colour-coded according to its light curve shape and content in the Fourier domain with purple indicating mono-periodic stars, olive multi-periodic, and red hybrid pulsators. The stars classified as other are shown in cyan and are discussed in the main text.

Fig. B.2. Hertzsprung-Russell diagram similar to Fig. 3 featuring only the stars with identified modes from multi-colour photometry and colour-coded by their mode degree: l = 0 (green), l = 1 (blue), l = 2 (orange).

Figure B.3. shows the distribution of amplitude ratios for the set of mode degrees as determined by our probabilistic approach. Qualitatively, the diagram resembles Fig. 9. While the Gaia amplitude ratios of the non-radial pulsators overlap significantly, the TESS amplitudes for the non-radial modes are well separated. As expected, the amplitude ratios of the radial modes drop stronger over the whole wavelength range than the amplitude ratios of the non-radial modes.

Fig. B.3. Same as Fig. 9 but based on the probabilistic assignment of mode identifications from multi-colour photometry for all modes with p>0.6.

Appendix C: Identification of the radial order

The grid modelling presented in Sect. 6.3 needs a full mode identification in terms of n, l, and m. Since our model selection is based on only the dominant mode degree, the radial order is still ambiguous and cannot be fitted in our scheme. To find the most likely radial order, we search the best fitting model for n∈{−2,−1,0,1,2} given the mode degree as determined by the multi-colour photometry.

For each star, we produced a diagnostic plot similar to Fig. C.1.. It shows for all relevant radial order out of the selected five radial orders the positions of all twelve considered low-order radial and zonal p- and g-modes overplotted on the frequency spectrum. In these diagrams, we examined the positions of these modes with respect to all identified (and unidentified low S/N) mode frequencies.

Fig. C.1. Diagnostic plot to identify the radial order n. In the background, we show in black the frequency spectrum of TIC 105756939. The thin red vertical lines indicate the positions of extracted frequencies. The short coloured lines show the positions of the model predictions for the modes of different radial orders. From bottom to top, we show we show n=−2,−1,1,2 in orange, blue, purple, brown, respectively.

In the chosen example, the radial order n = 1 is able to approximately explain detected modes at frequencies both higher and lower than the identified dominant mode. However, it fails to reproduce the highest detected frequency in Fig. C.1.. Hence, we can conclude that TIC 105756939 is pulsating with (n,l) = (1,1).

Appendix D: Properties of the modelling results

D.1. Recovered values

By comparing the input and recovered values from the grid modelling, we can assess the accuracy and understand potential biases of the modelling. The three panels of Fig. D.1. show the measured and the differences to the recovered log T_eff, log L, and f₀. There is a slight bias towards higher effective temperatures as the majority of best fitting modes are located below zero. A similar effect can be found for the luminosity, where the lower-luminosity stars tend to be fitted with higher-luminosity models. We also note that stars close to the edges of the model grid are affected by the lack of models and their luminosity is shifted towards the model grid, leading to the structure seen at low and high luminosities. However, all models are within the estimated uncertainty of the measurements. The frequencies are well recovered and deviate less than the uncertainty interval (0.1 d⁻¹) from the measurements.

Fig. D.1. Comparison between the measured values and the properties recovered from the grid modelling. From top to bottom we show the effective temperature, the luminosity, and the pulsation frequency of the dominant mode. The dotted lines mark the uncertainty interval used in our fitting procedure. The colour-coding indicates the main pulsation mode: l = 0 green, l = 1 blue, l = 2 orange.

D.2. Ages

Due to their short lifetimes, massive stars are young in a Galactic context and can often still be found in their birth environment such as OB associations or open clusters. The Gaia mission has also enabled the large scale identification and description of open clusters, as well the attribution of individual stars to these groups. We use the catalogues from Cantat-Gaudin et al. (2020), Kounkel et al. (2020), and Hunt & Reffert (2023) to identify β Cep stars in open clusters and get an independent handle on their ages. We find 34 (15 %) of our sample stars to be hosted by an open cluster or co-moving group.

All the three open cluster catalogues infer the cluster ages based on their colour-magnitude diagrams and isochrones. Fig. D.2. compares our asteroseismic age to these isochronal ages. Overall our asteroseismic ages are younger than the isochronal estimates, pointing to systematic uncertainty caused by the different input physics of the isochrone and asteroseismic models in the adopted grids. However, taking into account the typical ages uncertainty of 0.2 dex for the isochronal measurements (Hunt & Reffert 2023) and 0.1 dex for the asteroseismic age, the ages are mostly in agreement.

Fig. D.2. Comparison between ages for stars in open clusters and our asteroseismic ages. Members of NGC 637 are shown with squares and those of NGC 5606 with triangles. The typical uncertainty is given in the upper left corner.

For the two open cluster NGC 637 and NGC 5606 multiple β Cep members are in our sample. We show their members with distinct symbols in Fig. D.2.. As these stars are part of coeval populations, their ages should agree within the uncertainties, unless they would be binary or multiple products. We find their asteroseismic ages exhibiting a large spread even beyond the uncertainties in the case of NGC 5606. Detailed modelling of these stars should take into account their coevality to better constrain the asteroseismic inferences and investigate if some of them could be multiples or merger products.

All Tables

All Figures

Fig. 1. Stacked representation of the extracted independent, significant TESS frequencies for 212 of the Gaia-detected β Cep stars, sorted by the TESS frequency with the largest amplitude. The colour-scale indicates the amplitude of each detected signal on a logarithmic scale. This figure includes only the stars whose main frequency occurs in the plotted β Cep range. Three stars have lower main frequencies and turn out to be SPB stars or β Cep and SPB hybrid pulsators with dominant high-order g-modes, or are dominated by other variability than pulsations.

In the text

Fig. 2. Comparison between the highest amplitude Gaia and TESS frequencies for the 215 stars with significant TESS frequencies. Top: Dominant TESS frequency fTESS against the Gaia derived frequency fGaia. The black symbols indicate stars for which we found a match in frequency, the blue symbols show stars matching a frequency other than the dominant signal, and the red symbols indicate the seven stars for which the Gaia frequency does not match any frequency in the TESS data. The lines connect the offset frequencies to the matched counterpart. For three stars the main frequency is below the plotting range. Bottom: Difference between the Gaia and TESS frequencies. The black dashed line indicates our chosen threshold for matching signals (Δf<0.006 d−1).

In the text

Fig. 3. Hertzsprung-Russell diagram based on Gaia data of all 321 β Cep stars from the literature (circles) and the 121 new stars in our Gaia-informed sample (squares, new detections). The typical uncertainty is indicated in the lower left. The black outlines show the instability regions for p-modes (solid, β Cep, l = 0−2), and g-modes (dashed, SPB, l = 1−2) based on the calculations from Burssens et al. (2020). The main sequence evolutionary tracks of Burssens et al. (2020) are shown in grey. The colour-coding gives the main pulsation frequency of the stars cut off at 12 d−1. Stars marked with a red cross do not exhibit β Cep variability in TESS. The top axis marks the spectral type.

In the text

Fig. 4. Examples of single sector TESS light curves (left) of β Cep variables of the different light curve sub-classes and their associated periodograms (right). The top light curve is a typical mono-periodic pulsator with a constant pulsation amplitude and a single peak in the frequency spectrum. The second and third panels show multi-periodic pulsators with different lengths of their beating patterns. The bottom panel shows a typical β Cep and SPB hybrid pulsator.

In the text

Fig. 5. Beating period of the two main frequencies (Pbeat) against the normalised percentile amplitude ratio (Avar). Light curves manually classified as mono-periodic are shown in purple, multi-periodic in olive, hybrids in red, and other cases in cyan. The dashed lines indicate the approximate dividing lines of the groups of pulsating stars. Uncertainties on the beating period are within the symbol sizes (Avar has no formal uncertainty due to its definition).

In the text

	Fig. 6. Main pulsation frequency f1 against the projected surface rotation frequency frotsin i of the 162 stars with available v sin i in the Gaia ESP-HS data. Uncertainties in the pulsation frequency are typically within the symbol size.
In the text

	Fig. 7. Pulsation amplitude ratio AT/ABP against ARP/ABP for the main mode detected in the Gaia photometry of 205 β Cep stars. The colour-coding is according to the amplitude ratio AG/ABP. Crosses mark a few well-studied bright pulsators with secure mode identifications from the literature.
In the text

Fig. 8. Examples of rotationally split multiplets found in the amplitude spectra of the TESS data. The top panel shows an l = 1 mode split into three azimuthal components. The bottom panel shows one of the few identified complete l = 2 multiplets. The splittings are marked with grey markers connected by a line. The red markers at the top of each panel indicate the extracted frequencies from pre-whitening. The dashed lines mark the dominant frequency in the Gaia DR3 light curve.

In the text

	Fig. 9. Violin plots of amplitude ratios for identified pulsation modes over the filter's effective wavelengths. The lines connect the mean values for each mode degree. The filter names are given at the top. The GRP and T filters are displaced by ±30 nm to avoid overlapping distributions.
In the text

	Fig. 10. Light curves folded over the dominant frequency of two stars with clear standstills in our sample (stalling of the light curve near phase 0.5). The folded light curve is binned into 35 bins, shown as black circles, to which we fitted a sinusoidal with a fixed frequency for reference (red). The top panel shows the archetypical standstill star BW Vul (Walker 1954).
In the text

Fig. 11. Amplitude ratio diagram based on ARP/ABP and AT/ABP, similar to Fig. 7. The colour-coding shows the mode degree as determined with our probabilistic framework: l = 0 (green), l = 1 (blue), l = 2 (orange). The grey symbols denote stars with p<0.6. The symbol size is proportional to the probability of the mode identification. The stars marked with black outlines constitute the underlying sample from Fig. 9.

In the text

Fig. 12. Observed mode splittings against the projected surface rotation rate (frotsin i, encircled in thin grey) or photometric surface rotation rate (frot, encircled in thick solid black). The three lines show the splitting values when equal to the surface rotation rate (solid) as well as one-half (dotted) and twice (dashed) the surface rotation rate. The uncertainties on the splittings indicate the maximum difference of the individual splittings to the mean splitting of the multiplet and are a measure for the asymmetry. The colours show our deterministic identifications as l = 1 (blue) and l = 2 (orange) modes.

In the text

Fig. 13. Asymmetries of 33 rotationally split multiplets δA against their central frequency f0 colour-coded by the projected surface rotation rate. The square symbols mark stars with measured photometric surface rotation rates (frot). The symbols connected by a line belong to the same (l = 2) multiplet. The black outlines indicate stars where the Gaia frequencies do not match the TESS frequencies within our threshold. The measurement uncertainties are within the symbol sizes.

In the text

Fig. 14. Distribution of β Cep parameters for the modelled stars. Each panel contains the histograms for the full sample (grey) and each degree l = 0 (green), l = 1 (blue), and l = 2 (orange). From top to bottom, these panels show the initial mass distribution, the normalised core hydrogen fraction $X_c^{\prime }$, the stellar age in years, the relative convective core mass, the radius, and the K-value.

In the text

Fig. 15. Normalised convective core mass Mcore/Mini over the normalised core hydrogen mass fraction $X_c^{\prime }$. The stars are colour-coded by their primary degree: l = 0 (green), l = 1 (blue), and l = 2 (orange). The shape represents the radial order: g2-modes are represented by thin diamonds, g1 by thin triangles, radial fundamental modes by triangles, p1 by circles, and p2 by squares. The grey tracks in the background show the evolution for the indicated masses with fov = 0.005 (solid, labelled at low $X_c^{\prime }$) and fov = 0.035 (dashed, labelled at high $X_c^{\prime }$). The typical uncertainty is indicated in the upper right corner.

In the text

Fig. 16. Period-luminosity relations for stars with modes of different degrees: l = 0 (green), l = 1 (blue), and l = 2 (orange). The period-luminosity relations include only the radial fundamental mode (n = 0) for radial pulsators and p-modes for the non-radial pulsators. The black line shows the relation of the Shi et al. (2024) sample. For comparison, a similar period-luminosity relation for our sample is plotted in purple (dashed). The green shaded area shows the uncertainty band of the fundamental radial-mode pulsators. The small black data points in the background show dominant g-mode and higher-order radial pulsators excluded from the fits.

In the text

Fig. 17. Envelope rotation Ωenv and derived parameters against observed frequencies and estimated core hydrogen fractions for 25 stars. The top panel shows the envelope rotation rate against the zonal mode frequency of the multiplet (f0) and the second panel against the projected surface rotation rate frotsin i for the stars having this quantity available. The third and fourth panels display derived rotational properties against the evolutionary stage. In the former, we show an upper limit for the level of differential rotation (Ωenv/Ωsurf) and in the bottom panel the spin parameter (s) of the mode in the envelope. The colour-coding corresponds to the mode degree: l = 1 in blue and l = 2 in orange.

In the text

	Fig. 18. Histograms of ϵ2 (top) giving the influence of the centrifugal force with respect to gravity and the inverse of the K-value (bottom).
In the text

	Fig. A.1. Example light curve before (top) and after (bottom) the final detrending. Sections highlighted in red in the top panel have been removed.
In the text

Fig. A.2. Example of a single pixel TESS light curve. The target TIC 264026227, member of the open cluster NGC 5606, is contaminated by several nearby stars. Top left: The TESS FFI cut-out of the field around TIC 264026227, with the pixel used for or final light curve marked in red. The colour scale shows the flux intensity with yellow being brighter pixels. Overplotted in each pixel is its light curve. Top right: Same as left but with the periodogram overplotted. Bottom: Final light curve shows clear β Cep variability after additional cleaning and detrending beyond what is shown in the top left panel.

In the text

Fig. B.1. Hertzsprung-Russell diagram similar to Fig. 3 but including only the 216 stars from our Gaia-detected sample having TESS light curves. Each star is colour-coded according to its light curve shape and content in the Fourier domain with purple indicating mono-periodic stars, olive multi-periodic, and red hybrid pulsators. The stars classified as other are shown in cyan and are discussed in the main text.

In the text

	Fig. B.2. Hertzsprung-Russell diagram similar to Fig. 3 featuring only the stars with identified modes from multi-colour photometry and colour-coded by their mode degree: l = 0 (green), l = 1 (blue), l = 2 (orange).
In the text

	Fig. B.3. Same as Fig. 9 but based on the probabilistic assignment of mode identifications from multi-colour photometry for all modes with p>0.6.
In the text

Fig. C.1. Diagnostic plot to identify the radial order n. In the background, we show in black the frequency spectrum of TIC 105756939. The thin red vertical lines indicate the positions of extracted frequencies. The short coloured lines show the positions of the model predictions for the modes of different radial orders. From bottom to top, we show we show n=−2,−1,1,2 in orange, blue, purple, brown, respectively.

In the text

	Fig. D.1. Comparison between the measured values and the properties recovered from the grid modelling. From top to bottom we show the effective temperature, the luminosity, and the pulsation frequency of the dominant mode. The dotted lines mark the uncertainty interval used in our fitting procedure. The colour-coding indicates the main pulsation mode: l = 0 green, l = 1 blue, l = 2 orange.
In the text

	Fig. D.2. Comparison between ages for stars in open clusters and our asteroseismic ages. Members of NGC 637 are shown with squares and those of NGC 5606 with triangles. The typical uncertainty is given in the upper left corner.
In the text