1 Introduction

Massive stars play a pivotal role in the chemical, radiative and kinetic evolution of their host galaxies (Maeder 2009; Langer 2012). This is because massive stars, typically defined as having birth masses larger than \(M_{\mathrm{ini}} \gtrsim 8\) M⊙, end their lives as a violent supernova explosion before forming a neutron star or black hole. Hence massive stars are an important driver of feedback to their host galaxy (Bromm et al. 2009; de Rossi et al. 2010). In addition, during their lives massive stars have strong radiatively driven winds that contribute to interstellar feedback and new generations of star formation. Another important aspect of massive stars is that they are commonly found in binary or higher-order multiple systems (e.g. Sana et al. 2012), with dynamical interaction and mass transfer between component stars of a multiple system having a significant impact on a star’s evolution (e.g. de Mink et al. 2013; Schneider et al. 2015; Marchant and Bodensteiner 2023). This makes massive stars, in general, complex but fruitful laboratories for testing a variety of physics.

It is clear that understanding massive star evolution is critical for a broad range of topics within stellar astrophysics. However, in addition to the physics associated to stellar winds and binarity, all evolution models of massive stars currently suffer from large theoretical uncertainties associated with the interior physical processes of rotation and mixing (Langer 2012; Martins and Palacios 2013). In particular, the amount and shape of mixing at the boundary of convective and radiative regions, commonly referred to as convective boundary mixing (CBM), is currently controlled by uncalibrated parameterisations in evolution models (Bowman 2020; Johnston 2021; Anders and Pedersen 2023). Such uncertainties are common to all massive stars regardless of their metallicity, rotation rate or presence of a binary companion, since all massive stars have a hydrogen-burning convective core during the main sequence. Because of the unknown amount and shape of CBM at the boundary of convective cores, the uncertainties in core masses and ages of massive stars can easily reach 100% (Bowman 2020; Johnston 2021). Furthermore, the rotation rate of a massive star influences the amount of interior mixing and interstellar feedback (see e.g. Hirschi et al. 2005; Ekström et al. 2012; Chieffi and Limongi 2013; Georgy et al. 2013), but also the efficiency of angular momentum transport (Aerts et al. 2019a).

Today, evolution models include numerical prescriptions for interior rotation, mixing and angular momentum transport physics based on important theoretical work (see e.g. Zahn 1991). However, testing the validity and applicability of such prescriptions using observational constraints has some catching up to do. A powerful method for such a task is asteroseismology (Aerts et al. 2010; Kurtz 2022), which is the method of quantitative comparison of observed pulsation mode frequencies to a grid of evolution models and their theoretically predicted pulsation mode frequencies to determine the best-fitting model (Aerts 2021). Asteroseismology has been highly successful when applied to low-mass stars, such as red giants (Chaplin and Miglio 2013; Hekker and Christensen-Dalsgaard 2017; García and Ballot 2019). In particularly, asteroseismology has concretely demonstrated large discrepancies in the observed and theoretically predicted radial rotation profiles of low-mass stars, thus calling for improved theory and numerical prescriptions for angular momentum transport across stellar evolution (Aerts et al. 2019a). On the other hand, asteroseismology of massive stars is lagging behind low-mass stars predominantly because such stars are intrinsically rare in the Universe, but also because high-precision, long-duration time-series data for a large sample of stars was not available until recently (Bowman 2020). It is primarily thanks to the ongoing NASA TESS mission (Ricker et al. 2015) that the situation has changed, with the TESS mission providing high-precision light curves of thousands of pulsating massive stars across the entire sky.

In this invited review, I discuss recent advances in our understanding of massive star interiors thanks to asteroseismology. As this review is part of the Springer Nature 2023 Astronomy Prize Awardees Collection, I pay particular attention to work that I have personally been involved with. Indeed, as mentioned previously, it is thanks to the ongoing NASA TESS mission that has revolutionised the observational data-driven paradigm shift of asteroseismology to massive stars. The review by Bowman (2020) discussed the constraints on interior rotation and mixing based on massive star asteroseismology determined up until the launch of the TESS mission. In this complementary review, I focus on the advances made using TESS data and the new questions and puzzles that have been revealed since circa 2020.

2 Diverse variability in massive stars revealed by space photometry

It was the launch of space telescopes in recent decades, such as the MOST (Walker et al. 2003), CoRoT (Auvergne et al. 2009), Kepler (Koch et al. 2010; Borucki et al. 2010), BRITE (Weiss et al. 2014), and TESS (Ricker et al. 2015) missions, that opened pandora’s box for variability studies of massive stars. This is because the typical variability time scales, such as from rotational modulation, pulsations or binarity, can be of the order of a few days to several hours for massive stars making them difficult to probe using ground-based telescopes. Herculean efforts that assembled many nights of observations of individual massive stars with some of the world’s largest ground-based telescopes have yielded significant asteroseismic results (e.g. Aerts et al. 2003; Handler et al. 2004), but only for a handful of massive stars. For an in-depth discussion and examples of asteroseismology of massive stars using ground-based telescopes, I refer the reader to the review by Bowman (2020).

The MOST nanosatellite performed valuable initial studies of massive star variability, including the characterisation of pulsations (Walker et al. 2005b,a; Saio et al. 2006), rotation (Ramiaramanantsoa et al. 2014), and rotation period spin-down caused by large-scale magnetic fields (Townsend et al. 2013). However, it was space photometry assembled by CoRoT that demonstrated the feasibility of massive star asteroseismology at a large scale, although the CoRoT mission was modest in terms of sample size by today’s standards. One of the important contributions from the CoRoT mission was the first detection of a gravity-mode period spacing in a massive star (Degroote et al. 2010). Later, the nominal Kepler mission data indirectly included a few massive stars (see e.g. Aerts et al. 2017; Pope et al. 2019) but generally avoided them to focus instead on low-mass stars with potential exoplanets. As an example, Kepler data allowed pulsations to be characterised in high-mass eclipsing binaries with far better photometric precision than possible from the ground (Tkachenko et al. 2012, 2016). In the revised K2 mission (Howell et al. 2014), a few hundred massive stars were included (Bowman et al. 2019b), but K2 light curves were maximally 80 d in length and of lesser quality than nominal Kepler mission light curves. On the other hand, the BRITE-constellation of cube-satellites, launched over a decade ago and continues to assemble long-term light curves of selected bright massive stars. For a review of the BRITE mission and an ensemble analysis of its first data release, I refer the reader to Weiss et al. (2021) and Zwintz et al. (2023), respectively.

In this review, I focus on work motivated primarily by newly available light curves from the ongoing TESS mission (Ricker et al. 2015), which has a fundamentally different observational strategy compared to previous time-series photometry space missions. The CoRoT and Kepler missions focussed on assembling light curves within fixed fields in the sky, pre-chosen to contain targets with high scientific potential. Whereas the TESS mission is operating as a fully community-driven survey in which essentially all stars brighter than approximately \(V \lesssim 14\) mag in the sky have a light curve spanning at least 27 d with a cadence of 30 min in the nominal mission. As of the extended and second extended TESS missions, this cadence was significantly improved to be 10 min and 200 sec, respectively, in the so-called full-frame images (FFIs).

The community are able to propose a limited number of high-priority targets for even shorter cadences of 2 min and even 20 sec in the form of TESS Guest Investigator (GI) proposals.Footnote 1 Another difference and advantage is that all TESS data products are made immediately publicly available with no proprietary period. Such an approach has facilitated many serendipitous discoveries and has the major advantage of providing a large sample of massive stars across the whole sky, albeit at the expense of providing typically shorter continuous light curves compared to CoRoT and Kepler.

3 Discovery and exploitation of stochastic low-frequency variability in massive stars

A major step forward in our understanding of massive stars is in the ongoing work to comprehend stochastic low-frequency (SLF) variability in their light curves. The first discussion of SLF variability, although it was simply referred to as ‘astrophysical red noise’, in time-series photometry was based on the only three O stars observed by CoRoT (Blomme et al. 2011). However, the conclusions were speculative as to the physical mechanism giving rise to such a phenomenon.

Later, the presence of SLF variability at the surface of these three massive O stars was postulated to be connected to convectively driven gravity waves by Aerts and Rogers (2015) based on a qualitative comparison of the observed and theoretical frequency spectra of gravity waves driven by core convection predicted by 2D hydrodynamical simulations (Rogers et al. 2013; Rogers 2015). At the time, these were the only numerical simulations with the necessary global spatial scale, resolution and input physics that permitted a physical comparison to the inferred gravity wave spectra of observed stars. It was this hypothesis that sparked much of the subsequent research into finding and characterising the signatures of gravity waves in many more massive stars using both photometry and spectroscopy.

3.1 Constraints from time series photometry

The first large-scale search for SLF variability using dozens of early-type stars observed by CoRoT demonstrated that a common morphology exists in the frequency spectra of such stars (Bowman et al. 2019a). The shape of the SLF variability in the amplitude spectra of early-type stars is well described by a semi-Lorentzian function:

$$ \alpha \left ( \nu \right ) = \frac{ \alpha _{0} }{ 1 + \left ( \frac{\nu}{\nu _{\mathrm{char}}} \right )^{\gamma}} + C_{\mathrm{w}} ~ , $$
(1)

where \(\alpha _{0}\) represents the amplitude at a frequency of zero, \(\gamma \) is the logarithmic amplitude gradient, \(\nu _{\mathrm{char}}\) is the characteristic frequency, and \(C_{\mathrm{w}}\) is a frequency-independent (i.e. white) noise term (Bowman et al. 2019a). The characteristic frequency, which represents the inverse of the typical timescale of the variability, is of order several d−1 for massive dwarf stars (Bowman et al. 2019a). Moreover, this was the first time that a quantitative approach to characterising SLF variability in observations was used, with results shown to be consistent with the theoretical frequency spectra of gravity waves from 2D hydrodynamical simulations (Rogers et al. 2013; Rogers 2015) and predictions of standing wave gravity eigenfrequencies of such stars (Bowman et al. 2019a). Examples of the SLF variability morphology in the frequency spectra of two O stars, HD 46223 and HD 46150, observed by CoRoT and TESS are shown in Fig. 1, which demonstrate the broad frequency range of SLF variability, with significant periods ranging from minutes to several days and includes the reported \(\nu _{\mathrm{char}}\) values from Bowman and Dorn-Wallenstein (2022).

Fig. 1
figure 1

Example of SLF variability profiles in the frequency spectra of two O stars, HD 46223 and HD 46150, observed by CoRoT and later re-observed by the TESS mission in 2-min cadence in its sectors 6 and 33 (if available). These two examples demonstrate the broad frequency range over which SLF variability is significant (i.e. periods of tens of days to minutes). The reported \(\nu _{\mathrm{char}}\) values and their 2\(\sigma \) confidence intervals from Bowman and Dorn-Wallenstein (2022) are also shown as the yellow shaded regions

Full size image

With a larger sample of massive stars available from the K2 and TESS missions, it was later shown that SLF variability is essentially ubiquitous in massive stars when sufficiently high-quality light curves are available (Bowman et al. 2019b). A large sample of over 160 massive stars observed by K2 and TESS, which included metal-poor extra-galactic and metal-rich galactic stars across the entire sky, again demonstrated that all stars had a similar morphology in their frequency spectra. Two important conclusions from this large sample are: (i) the morphology of SLF variability and specifically the inferred \(\nu _{\mathrm{char}}\) parameter is seemingly insensitive to metallicity; and (ii) the amplitude and characteristic frequency, \(\nu _{\mathrm{char}}\), of SLF variability depends on the brightness and luminosity of the star (Bowman et al. 2019b).

Using a combination of TESS light curves and high-resolution spectroscopy, including from the HERMES spectrograph mounted on the 1.2-m KU Leuven telescope on La Palma (Raskin et al. 2011), it was demonstrated that a strong correlation exists between the morphology of observed SLF variability and the location of a star in the Hertzsprung–Russell (HR) diagram (Bowman et al. 2020). More specifically, this means that the morphology of SLF variability directly probes the mass and age of a massive star (Bowman et al. 2020). Such a discovery has opened up a new sub-branch of asteroseismology in which constraints on the fundamental properties of massive stars can be achieved without the prerequisite of identifying the spherical geometry of individual pulsation mode frequencies. A demonstration of the change in the measured \(\nu _{\mathrm{char}}\) parameter as stars evolve during the main sequence based on the results of Bowman et al. (2020) is shown in Fig. 2. The conclusion from this transition is that the physical mechanism causing SLF variability tends to produce longer periods for older stars.

Fig. 2
figure 2

Relationship between the measured characteristic frequency, \(\nu _{\mathrm{char}}\), of SLF variability and location in the spectroscopic HR diagram for 30 galactic O stars (Bowman et al. 2020; Bowman and Dorn-Wallenstein 2022), demonstrating that as stars evolve the value of \(\nu _{\mathrm{char}}\) decreases. Evolutionary tracks starting at the ZAMS are shown as solid grey lines with initial masses given in units of M⊙ (Burssens et al. 2020). A typical spectroscopic error bar for the sample is shown in the top-left corner

Full size image

More recently, the correlation between SLF variability morphology and location of massive stars in the HR diagram was confirmed using the novel light-curve fitting method of Gaussian process (GP) regression (Bowman and Dorn-Wallenstein 2022). A GP regression model describes stochastic variability in a light curve by fitting hyperparameters and defining a covariance matrix. This provides posterior distributions that represent structures in their covariance matrices, which are conditioned on the input data to evaluate the optimum representation of the input time series data (Rasmussen and Williams 2006). For a detailed discussion of GP regression methods and applications to astronomical time series data, I refer the reader to Aigrain and Foreman-Mackey (2023). In the study by Bowman and Dorn-Wallenstein (2022), a GP kernel of a damped simple harmonic oscillator was used within the celerite2Footnote 2 software package (Foreman-Mackey et al. 2017), which allowed \(\nu _{\mathrm{char}}\) as well as the quality factor, \(Q\), of the input light curves to be efficiently and accurately inferred. Interestingly, the GP regression methodology revealed that stars are more stochastic (i.e. lower \(Q\) values) near the zero-age main sequence and are more quasi-periodic (i.e. higher \(Q\) values) near the terminal-age main sequence (Bowman and Dorn-Wallenstein 2022). This transition further demonstrates that the physical mechanism responsible for SLF variability not only depends on age, but also that SLF variability has an age-dependent stochasticity and coherency.

3.2 Constraints from complementary spectroscopy

Complementary to photometric time series observations, spectroscopic variability also provides valuable constraints on the physical processes within massive stars. Spectral line profile variability (LPV) and the time-dependence of non-rotational spectral line profile broadening mechanisms, for example macroturbulence, has historically been linked to pulsations (Lucy 1976; Fullerton et al. 1996; Howarth et al. 1997; Aerts et al. 2009). The triangular shapes of spectral lines in massive stars cannot exclusively be explained by rotational broadening (Simón-Díaz and Herrero 2014), but despite the necessity of including a macroturbulent broadening term, the physics explaining this ‘nuisance’ term in spectroscopy has remained illusive.

In a detailed study of macroturbulence in galactic massive stars, it was found that early- to late-B main-sequence stars have a diverse range in macroturbulent broadening values, whereas main-sequence O and B supergiants typically have macroturbulence rather than rotation as the dominant broadening mechanism in their spectral lines (Simón-Díaz et al. 2017; de Burgos et al. 2023a,b). Large amounts of macroturbulence are also necessary to fit the spectral lines of supergiants in the Large Magellanic Cloud (LMC) galaxy (Serebriakova et al. 2023), demonstrating its importance across different metallicity regimes. In order to explain such large values of macroturbulence, which can exceed 100 km s−1 for O stars, non-radial gravity-mode pulsations have been shown to provide the necessary predominantly tangential velocity field near the stellar surface (Aerts et al. 2009). Such pulsations can be excited by an opacity heat-engine mechanism and/or turbulent pressure fluctuations within massive star envelopes (Aerts et al. 2009; Grassitelli et al. 2015). Alternatively, the dynamics of turbulent massive star envelopes can give rise to large velocity fields (Schultz et al. 2022; Jiang 2023), but the interplay of these non-mutually exclusive mechanisms is yet to be investigated.

In addition to coherent gravity-mode pulsations described above, it is also plausible for an ensemble of stochastic gravity waves to produce the necessary predominantly tangential velocity field at the stellar surface to explain LPV and macroturbulence. This has been demonstrated using numerical simulations (Aerts and Rogers 2015). In support of this conclusion, the amplitude of SLF variability in time-series photometry was also found to strongly correlate with the amount of macroturbulence in a large sample of massive stars (Bowman et al. 2020). Therefore, the physics that causes SLF variability in photometry is likely also responsible for macroturbulence in spectroscopy and connected to pulsations.

3.3 Plausible physical mechanisms

Inspired by the new observational discovery that SLF variability is seemingly ubiquitous across a wide range in stellar mass and age, there have been several theoretical interpretations offered in the literature each with their own advantages and disadvantages. As mentioned previously, gravity waves excited by turbulent core convection predicted by 2D and 3D hydrodynamical simulations produce a frequency spectrum with similar morphology to that of SLF variability in time-series observations (Rogers et al. 2013; Rogers 2015; Rogers and McElwaine 2017; Edelmann et al. 2019; Ratnasingam et al. 2019, 2020, 2023; Horst et al. 2020; Varghese et al. 2023; Vanon et al. 2023; Herwig et al. 2023; Thompson et al. 2023).

Of course, all hydrodynamical studies have specific numerical setups and assumptions in an attempt to most accurately reproduce reality. It is an advisable and valuable exercise to compare different numerical methods, especially if different results are obtained from codes attempting to solve the same physical problem (see e.g. Andrassy et al. 2022; Lecoanet and Edelmann 2023). The differences in some setups and apparent disagreement in numerical predictions as outputs of specific simulations, such as the frequency spectrum of core-excited gravity waves, has led to the physical mechanism behind SLF variability being hotly debated in the literature. For example, some hydrodynamical simulations with different numerical setups conclude that gravity waves excited by turbulent core convection do not reach the stellar surface with observable amplitudes because of radiative damping (Lecoanet et al. 2019, 2021; Le Saux et al. 2023; Anders et al. 2023). However, almost all simulations are performed using zero or very small rotation rates compared to observed stars, with rotation being both an important restoring force but also a strong contributory factor for setting the (non-linear) amplitudes of waves within stellar interiors. In their 3D hydrodynamical simulations, Anders et al. (2023) predict around an order of magnitude larger wave amplitudes for moderately rotating simulations compared to slowly rotating simulations, with predicted wave amplitudes that approach the measured amplitude of SLF variability in observed rotating stars. Therefore, for an accurate comparison of numerical predictions from simulations to measured parameters of observed stars, realistic rotation rates and large-scale spherical simulations are necessary before concluding if the cause of SLF variability is (or is not) gravity waves excited by turbulent core convection.

On the other hand, other studies hypothesise that if the physical mechanism of SLF variability is presumably not gravity waves excited by core convection, then it may be instead gravity waves and/or turbulence induced from sub-surface convection zones (Cantiello et al. 2009, 2021). Some hydrodynamical simulations of subsurface convection zones in massive stars also show them to be a plausible mechanism for SLF variability and macroturbulence in massive stars (Schultz et al. 2022). However, it is not clear whether all massive stars have well-defined subsurface convection zones, as they are strongly dependent on a star’s mass, rotation profile, age, opacity profile and metallicity. For example, main-sequence stars with masses between approximately 8 and 20 M⊙ at metallicities similar to the Small Magellanic Cloud (SMC) galaxy do not have significant subsurface convection zones according to the Rayleigh criterion for convection (Jermyn et al. 2022). Yet stars expected to not have subsurface convection zones show similar SLF variability in time series photometry to those that do (Bowman et al. 2024).

For completeness, the modulation of a star’s photosphere and atmosphere by an optically thick wind has also been shown to be a plausible mechanism for SLF variability in time-series photometry of post-main sequence massive stars (Krtička and Feldmeier 2018, 2021). However, such a mechanism is likely negligible for the optically thin winds of main-sequence massive stars. This is especially true for late-O and early-B main-sequence stars which only have a modest radiatively driven wind.

It is still an open question of which physical mechanisms that produce SLF variability (and macroturbulence) dominate in particular parameter regimes of the HR diagram. This is because the properties of (sub-surface) convection, stellar winds, but also the excitation, propagation and damping of gravity waves depend on stellar structure properties, such as mass, age, rotation rate and metallicity. From the observations, it is clear that any mechanism aiming to explain SLF variability in main-sequence massive stars must explain a broad range in mass (i.e. \(3 \lesssim M \lesssim 100\) M⊙), ages that span from the zero-age main sequence (ZAMS) to the terminal-age main sequence (TAMS) and beyond, rotation rates between zero and approaching critical, and metal-poor (i.e. \(Z \leq 0.002\)) and metal-rich (i.e. \(Z \geq 0.014\)) stars. More specifically, any mechanism must additionally explain the broad range in frequency of SLF variability and the correlation between its morphology and mass and age, and apparent insensitivity to metallicity (Bowman et al. 2019a,b, 2020; Bowman and Dorn-Wallenstein 2022; Bowman et al. 2024).

4 Physics of stellar interiors derived from forward asteroseismic modelling

The TESS mission data have and continue to provide a vast treasure trove for constraining the interior physical processes of stars using asteroseismology. Recent reviews focussing on low-to-intermediate-mass stars are available (Aerts 2021; Kurtz 2022), but in this work I focus on massive stars in particular. As mentioned previously, massive stars with coherent pulsations studied from the ground were limited to relatively high-amplitude \(\beta \) Cephei stars, which have masses larger than approximately 6 M⊙ and low-radial order pressure and gravity modes with periods of order several hours (Stankov and Handler 2005). A handful of asteroseismic studies of \(\beta \) Cep stars using ground-based data were able to constrain the interior rotation profiles and to a lesser extent the interior mixing properties (e.g. Aerts et al. 2003; Dupret et al. 2004). Such constraints remain extremely valuable to this day for improving our understanding of angular momentum transport within stellar interiors (e.g. Salmon et al. 2022b).

The constraints on interior rotation and mixing of massive stars from asteroseismology prior to the TESS mission were summarised by Bowman (2020). In this section, I provide an update by highlighting major advances in the field of massive star asteroseismology since circa 2020 and were hence not included in the review by Bowman (2020). It goes without saying that many of the discussed studies would not have been possible without open-source and excellently documented software tools, such as the MESAFootnote 3 evolution code (Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023) and the GYREFootnote 4 pulsation code (Townsend and Teitler 2013; Townsend et al. 2018; Townsend 2020; Goldstein and Townsend 2020; Sun et al. 2023).

Early ensemble papers using the first few sectors of TESS data demonstrated that there is a diverse range of variability mechanisms at work in massive stars, which includes coherent pulsations but also binarity and rotational modulation (Bowman et al. 2019b; Pedersen et al. 2019; David-Uraz et al. 2019; Southworth and Bowman 2022). Moreover, space photometry has demonstrated that \(\beta \) Cep stars are relatively common among massive stars and have a large range in pulsation periods and amplitudes (e.g. Burssens et al. 2019, 2020). With the addition of high-resolution spectroscopy, the pulsations seen in photometry can be placed in the context of a star’s location in the HR diagram, but this also allows some initial constraints on a star’s mass and evolutionary stage (Burssens et al. 2020; Gebruers et al. 2022; Serebriakova et al. 2023). All of the newly discovered pulsating massive stars within cycle 1 (i.e. TESS sectors 1-13) identified by Burssens et al. (2020) and for which high-resolution optical spectroscopy was available are plotted in the spectroscopic HR diagram in Fig. 3, which demonstrates the broad range in mass and age that pulsating massive stars occupy.

Fig. 3
figure 3

The spectroscopic HR diagram of pulsating massive stars identified by Burssens et al. (2020) in TESS sectors 1-13 that have high-resolution spectroscopy are shown as green circles. The ordinate axis shows spectroscopic luminosity such that \(\mathscr{L} := T_{{\mathrm{eff}}}^{4} / g\) (Langer and Kudritzki 2014). The blue and red hatched areas are the theoretical instability regions of pressure and gravity modes, respectively, calculated for non-rotating main-sequence stars, and non-rotating evolutionary tracks with solar metallicity (in units of \(M_{\odot}\)) are given as grey lines from Burssens et al. (2020). A typical error bar for a star’s location is also shown in the bottom-left corner. Figure adapted from Bowman (2020), his Fig. 4

Full size image

The first-light ensemble papers of massive star variability as seen by TESS have demonstrated two important conclusions: (i) the variability fraction is extremely high (\(> 95\%\)) in photometry of massive stars; and (ii) pulsating stars are often found outside of the theoretical pulsation instability regions, but also vice versa. There are two likely explanations for how non-pulsating stars can be found inside theoretical instability regions, and vice versa. The first is our incomplete understanding of stellar opacities and how accurate pulsation excitation calculations are based on opacity tables (Seaton et al. 1994; Iglesias and Rogers 1996). The excitation of coherent pulsations in massive stars is expected because of the opacity (\(\kappa \)) heat-engine mechanism operating at the iron bump in opacity at approximately 200,000 K (Dziembowski and Pamyatnykh 1993; Dziembowski et al. 1993). However, the modification of opacity table data in stellar models can drastically alter the parameter space in which coherent pulsations are (not) expected in the HR diagram (Daszyńska-Daszkiewicz and Walczak 2009; Moravveji 2016; Paxton et al. 2015; Walczak et al. 2015, 2019; Szewczuk and Daszyńska-Daszkiewicz 2018). On the other hand, the opacity within a star is indeed lower for low-metallicity stars such as those in the SMC and LMC galaxies. This has historically been the explanation of why a dearth of coherent pulsators are found outside of the metal-rich disk of our Milky Way galaxy (see e.g. Salmon et al. 2012).

The second important physical influence on the pulsational instability regions of massive stars is rotation. The importance of rotation for pulsations in massive stars is perhaps best illustrated using the sub-group of Be stars (e.g. Baade 1988; Huat et al. 2009; Kurtz et al. 2015; Neiner et al. 2020). These make up around 20% of B dwarfs and are defined by having emission lines in spectroscopy, although this is known to be a transient phenomenon. They are also generally rapidly rotating stars that seem to lack binary companions making them thought to be the products of binary evolution (see e.g. Abt and Cardona 1984; Bodensteiner et al. 2020b). However, Be stars in binary or multiple systems do exist (see e.g. Rivinius et al. 2013), so their origin remains unclear. Despite their interesting properties for binary evolution theory (e.g. Abdul-Masih et al. 2020; Shenar et al. 2020; Bodensteiner et al. 2020a; Marchant and Bodensteiner 2023), which is not the subject of this review, the rapid rotation of Be stars has been found to correlate with the incidence of non-radial pulsation using time-series photometry (Nazé et al. 2020; Labadie-Bartz et al. 2021, 2022). Most importantly, rapid rotation breaks the spherical symmetry of a star and spherical geometry of pulsation modes. Such an impact is significant enough to modify the location of instability regions of pulsating massive stars in the HR diagram (Townsend 2005a,b; Szewczuk and Daszyńska-Daszkiewicz 2017). Moreover, the high incidence of non-linear pulsations in Be stars (e.g. Kurtz et al. 2015; Neiner et al. 2020) also demonstrates the importance of including rotation in future hydrodynamical simulations studying wave generation, propagation and damping.

4.1 Interior rotation and mixing from asteroseismology

Generally speaking, forward asteroseismic modelling of massive stars requires sophisticated modelling techniques, with large grids of evolutionary models containing different prescriptions of physical processes to be constrained, and complex statistical treatment of random and systematic uncertainties (e.g. Daszyńska-Daszkiewicz and Walczak 2010; Daszyńska-Daszkiewicz et al. 2013; Aerts et al. 2019b; Handler et al. 2019; Salmon et al. 2022a). This is because the theoretical uncertainties – i.e. correlations and degeneracies among theoretical parameters to be constrained – are considerably larger than the observed pulsation mode frequencies in the era of space photometry (Bowman and Michielsen 2021).

In principle, forward asteroseismic modelling in a linear framework makes use of a set of observed pulsation mode frequencies with identified mode geometry, \(n\), \(\ell \), and \(m\), in terms of spherical harmonics. These are then compared using a statistical likelihood, or merit function, to predicted pulsation mode frequencies with the same spherical harmonic geometry from a grid of stellar evolution models. Such a theoretical grid contains various prescriptions for physical processes, such as mixing, which are controlled by free parameters. A minimal set of free parameters to be determined in this way includes a star’s mass (\(M\)), initial hydrogen and metal mass fractions (\(X\), \(Z\)), age (\(t\)) or, for example, the central hydrogen mass fraction, \(X_{\mathrm{c}}\) by proxy, and parameterisations for interior mixing profiles such as CBM and radiative envelope mixing, \(f_{\mathrm{CBM}}(r)\) and \(D_{\mathrm{env}}(r)\), respectively. The model grid is calculated to cover the required parameter space in the HR diagram inferred from spectroscopic estimates of effective temperature, \(T_{\mathrm{eff}}\), and luminosity, \(L\), or surface gravity by proxy, \(\log g\). Moreover, the determination of the interior radial rotation profile, \(\Omega (r)\), a priori using a mostly model-independent methodology, such as from rotational splitting within pulsation mode multiplets and/or gravity-mode period spacing patterns, is a major advantage for any subsequent forward asteroseismic modelling (see e.g. Aerts et al. 2018; Burssens et al. 2023).

Historically, a simplistic \(\chi ^{2}\) merit function has been used to ascertain the best fitting model in forward asteroseismic modelling of early-type stars (e.g. Moravveji et al. 2015, 2016). However, such an approach neglects to include the correlation structure of the theoretical model grid, with a more suitable, robust and superior statistical merit function for forward asteroseismic modelling being the Mahalanobis Distance (Aerts et al. 2018). Comparisons of using the Mahalanobis Distance and \(\chi ^{2}\) as merit functions in forward asteroseismic modelling of early-type stars demonstrated how the former is a significant improvement for probing the posterior parameter distribution (Michielsen et al. 2021). This is particularly true for forward asteroseismic modelling using gravity modes that are extremely sensitive to the mass and size of the convective core in main-sequence early-type stars (Michielsen et al. 2021). Building on this framework, detailed forward asteroseismic modelling that robustly handles random and systematic observational and theoretical uncertainties is now possible and facilitates high-precision constraints for pulsation massive stars (Bowman and Michielsen 2021). In some cases, a 1% precision on the core mass has been achieved (Michielsen et al. 2021; Bowman and Michielsen 2021). Nevertheless, asteroseismology across a wide range of masses and ages has concretely shown that core masses are underestimated in current evolution models (Johnston 2021), meaning that our understanding of CBM is incomplete (Anders and Pedersen 2023).

One of the biggest success stories of asteroseismology has been the determination of interior rotation rates of stars across all stages of stellar evolution. In cases of stars with rotationally split pulsation multiplets, demonstrating non-rigid radial rotation rates can be achieved using a largely model-independent method (Aerts et al. 2003; Kurtz et al. 2014; Saio et al. 2015). However the exact rotation profile is model dependent (Salmon et al. 2022b; Burssens et al. 2023). Given that asteroseismology is an ever advancing field of stellar astrophysics, reviews of interior rotation rates become out of date rather soon after publication. Nevertheless, the literature rotation rates for 1210 dwarfs and giant stars (Aerts et al. 2019a; Aerts 2021), and additionally 26 slowly pulsating B (SPB) stars (Pedersen et al. 2021), 611 \(\gamma \) Dor stars (Li et al. 2020), the few available massive stars from the review by Bowman (2020), and the new \(\beta \) Cep star HD 192575 (Burssens et al. 2023) are summarised in Fig. 4.

Fig. 4
figure 4

Top Panel: low- and intermediate-mass stars with asteroseismic estimates of near-core rotation rates (\(\Omega _{\mathrm{near}-\mathrm{core}}\)) versus asteroseismic \(\log g\) as a proxy for age. Middle Panel: sub-sample of stars at different evolutionary stages with an asteroseismic estimate of envelope rotation (\(\Omega _{\mathrm{env}}\)) as a function of asteroseismic \(\log g\). Uncertainties in \(\log g\) range between 0.01 and 0.30 dex, and uncertainties in rotation rates are smaller than the symbol size. Bottom Panel: same as middle panel but exclusively for massive stars, of which there are very few known in the literature. These figures were produced using asteroseismic near-core rotation measurements for 1210 low- and intermediate-mass stars from the review by Aerts et al. (2019a), 26 SPB stars from Pedersen et al. (2021), 611 \(\gamma \) Dor stars from Li et al. (2020), the few available massive stars from the review by Bowman (2020), and the new \(\beta \) Cep star HD 192575 (Burssens et al. 2023). The colours of all symbols are consistent across all panels and denote the evolutionary stage of the stars. Note that individual error bars for the massive stars are not plotted for clarity because they vary considerably within the subsample

Full size image

In the top panel of Fig. 4, the near-core rotation rates for low- and intermediate-mass (i.e. \(0.7 \leq M \leq 9\) M⊙) stars are shown as a function of their so-called asteroseismic \(\log g\), which is calculated from the asteroseismic masses and radii from forward asteroseismic modelling, and acts as a proxy for a star’s age or evolutionary stage. In the middle panel, the stars that also have an envelope rotation rate measured from asteroseismology or rotational modulation are shown as triangle symbols in addition to their corresponding circle symbols denoting the near-core rotation rate. The sample in Fig. 4 demonstrates the difference in differential rotation as a function of evolutionary stage. Most interesting is that core-hydrogen burning main-sequence stars and core-helium burning red clump stars have the smallest amount of differential rotation on average. This implies that the presence of a convective core has a direct impact on the efficiency of angular momentum transport within stars across their evolution (Aerts et al. 2019a; Aerts 2021).

In the bottom panel of Fig. 4, the massive stars (i.e. \(M \geq 9\) M⊙) with asteroseismic near-core and envelope rotation rates have been purposefully separated to demonstrate that not many exist in the literature. It is important to emphasise that, with the exception of HD 192575, there is a large amount of uncertainty for the rotation rates for the sub-sample of massive stars. This is because the other massive stars reported in the literature typically have only an upper or lower limit on the core and envelope rotation rates, which in some cases was derived from interpreting only a single or a few fiducial structure models (see e.g. Dupret et al. 2004). The exception is HD 192575 for which full error propagation was performed to derive its rotation profile. In the next section, I focus on the interesting and novel case of the \(\beta \) Cep star HD 192575, for which its rotation and mixing profiles have been measured to high precision thanks to excellent TESS mission data (Burssens et al. 2023).

4.2 The case of HD 192575

A valuable new addition to the sparsely populated bottom panel of Fig. 4 is the \(\beta \) Cep star HD 192575 (Burssens et al. 2023). This massive star was first discovered to pulsate in many low-radial order pressure and gravity modes thanks to TESS cycle 2 data. Even though it is a relatively bright star, it was not subject to an in-depth photometric study prior to the TESS mission. A quantitative spectroscopic analysis using multi-epoch high-resolution HERMES spectroscopy (Raskin et al. 2011) showed no evidence for a binary companion, and allowed high-precision atmospheric parameters of \(T_{\mathrm{eff}} = 23{,}900 \pm 900\text{ K}\), \(\log g = 3.65 \pm 0.15\), \(v\sin i = 27^{+6}_{-8}\) km s−1, and solar metallicity to be determined. Moreover, based on the Gaia DR3 parallax (Gaia Collaboration et al. 2016), a luminosity of \(L/L_{\odot} = 4.30 \pm 0.07\) was derived. These complementary constraints were used to hugely delimit the parameter space within the HR diagram for calculating structure models and corresponding pulsation mode frequencies with the MESA and GYRE codes (Burssens et al. 2023).

Frequency analysis of the 1-yr cycle 2 TESS light curve revealed several rotational multiplets, which are pulsation modes of the same radial order and angular degree but different azimuthal order, and have \(2\ell + 1\) \(m\)-components that are separated by the rotation rate within each mode-specific pulsation cavity. An illustration of the rotational multiplets in HD 192575 from a frequency analysis of the TESS cycle 2 data is shown in Fig. 5. The differences in the frequency spacings of different rotational multiplets allowed Burssens et al. (2023) to extract a core-to-surface radial rotation profile and conclusively exclude rigid rotation. Based on different assumed prescriptions for the shape of the shear layer between the core and surface, the amount of differential rotation was found to differ, such that core was inferred to be rotating between 1.4 and 6.3 times faster than the envelope (Burssens et al. 2023). This is the first time that such an analysis was performed. However, the most likely scenario is that the chemical composition gradient left behind by the receding convective core during the main sequence is the true shear layer, which yielded a core-to-surface rotation ratio of \(1.49^{+0.56}_{-0.33}\) (Burssens et al. 2023). This result has been subsequently shown to be consistent with 2D stellar structure models using the ESTER code (Espinosa Lara and Rieutord 2011, 2013; Rieutord et al. 2016) by Mombarg et al. (2023). Therefore, the asteroseismic modelling of HD 192575 serves as an extremely useful anchor point for testing angular momentum transport within the massive star regime (see e.g. Salmon et al. 2022b; Mombarg et al. 2023).

Fig. 5
figure 5

Top panel: 10-d section of the cycle 2 TESS light curve for the \(\beta \) Cep star HD 192575. Middle panel: frequency spectrum of the full 1-year cycle 2 TESS light curve, with the average rotational splitting of the utilised rotational multiplets labelled and colour-coded. Bottom panels: zoom-in of the same rotational multiplets and the corresponding azimuthal order, \(m\), of each component pulsation mode as assigned by Burssens et al. (2023) in the best-fitting solution from the forward asteroseismic modelling. Figure adapted from Burssens et al. (2023), their Fig. 1

Full size image

In addition to its interior rotation profile, Burssens et al. (2023) performed forward asteroseismic modelling to constrain the stellar mass, radius, age, and core mass to unprecedentedly high precision for a massive star. Ultimately the fractional uncertainties for mass, age and core mass were better than about 15%, which is remarkable given the complexity of the physics that was included in the input grid of evolution models. The derived parameters from only 1 yr of TESS data and their \(2\sigma \) confidence intervals obtained through forward asteroseismic modelling were a mass of \(M = 12.0 \pm 1.5\) M⊙, a radius of \(R = 9.1^{+0.8}_{-1.7}\) R⊙, an age of \(t = 17.0^{+4.7}_{-5.4}\) Myr, and a core mass of \(M_{{\mathrm{cc}}} = 2.9^{+0.5}_{-0.8}\) M⊙ (Burssens et al. 2023).

An important detail is that the work of Burssens et al. (2023) represents the first forward asteroseismic modelling study of a massive star to incorporate sources of observational and theoretical uncertainties, and implement using the Mahalanobis distance as a statistical merit function, thus appropriately handling the stellar model grid parameter correlations and degeneracies. As discussed briefly by Aerts and Tkachenko (2023), additional TESS data of HD 192575 during cycle 4 have recently become available, which when combined with the previous cycle 2 data provide improved frequency resolution and precision for forward asteroseismic modelling. With these higher quality data a more complete assignment of \(m\) values for components of the \(\ell =2\) rotational multiplets is possible (cf. Fig. 5). With these updated \(m\) values, a slightly more evolved star is found from forward asteroseismic modelling (Vanlaer et al. 2024), but the derived mass, age, radius and core mass remain within the \(2\sigma \) confidence intervals reported by Burssens et al. (2023). Interestingly, there are very slight asymmetries detectable in the \(\ell =2\) multiplets, which can also be exploited to constrain the interior magnetic field and geometry within the deep interior of HD 192575 (Vanlaer et al. 2024).

4.3 Magnetoasteroseismology

All of the asteroseismic modelling discussed in this review up until this point has been based on the use of non-magnetic stellar evolution models. Of course, magnetism is an important aspect of stellar structure and evolution since it can drastically change the interior rotation, mixing and angular momentum transport (Maeder and Meynet 2005; Keszthelyi et al. 2019, 2020, 2021). However, until recently all constraints on the strength and geometry of magnetic fields for massive stars were limited to the stellar surface. Even though the detection of rotational modulation in a massive star’s light curve suggests the presence of a magnetic field, this remains an indirect and limited detection method (e.g. Buysschaert et al. 2018b; David-Uraz et al. 2019, 2021). Using a more targeted and suitable strategy, large campaigns using spectropolarimetry have demonstrated that approximately 10% of massive dwarf stars have strong large-scale magnetic fields at their surfaces (Wade et al. 2016; Grunhut et al. 2017; Shultz et al. 2019). Such fields are presumed to extend much deeper than the stellar photosphere, and thought to be created during star formation (Mestel 1999; Neiner et al. 2015), or formed during binary mergers (Schneider et al. 2019, 2020).

There has been a great deal of theoretical work studying the impact of strong magnetic fields for gravity-mode pulsators. For example, Prat et al. (2019) and Prat et al. (2020) investigated the topology and obliquity of a strong magnetic field in the deep interior of a fiducial intermediate-mass main-sequence star, and specifically how it changes the morphology of gravity-mode period spacing patterns to take on a “saw-tooth” pattern. Later, Van Beeck et al. (2020) applied the same theoretical framework to various structure models for a range of masses and ages on the main-sequence and demonstrated that more evolved stars have the greater potential for detecting deep magnetic fields using pulsations. However, including a strong large-scale interior magnetic field using an a posteriori perturbative approach has limitations when comparing theoretical predictions to observations.

One particular magnetic star that has received much attention in the literature is the B3.5 V star HD 43317 which has a predominantly dipolar field of strength \(B_{\mathrm{p}} = 1312 \pm 332\) G detected and characterised using spectropolarimetry (Briquet et al. 2013; Buysschaert et al. 2017). HD 43317 is also a known pulsator with dozens of high-frequency gravity mode frequencies detected in its 150-d CoRoT light curve (Pápics et al. 2012; Buysschaert et al. 2018a). Interestingly, unlike other magnetic early-type dwarfs, HD 43317 is a relatively rapid rotator with a rotation period of \(0.897673 \pm 0.000004\) d (Buysschaert et al. 2017). Also, it has no signatures of being part of a binary system making HD 43317 a single, rapidly rotating, pulsating magnetic star. Although these characteristics do not make it unique, they do make it rare among the wider population of B dwarfs.

Thanks to space photometry and dedicated spectropolarimetry campaigns to charactertise the geometry of magnetic fields, there is a growing number of magnetic pulsating stars (Shibahashi and Aerts 2000; Briquet et al. 2012; Handler et al. 2012; Aerts et al. 2014; Wade et al. 2020). However, only two confirmed pulsating magnetic early-type stars have undergone forward asteroseismic modelling: the \(\beta \) Cep star V2052 Ophiuchi (Briquet et al. 2012) and the SPB star HD 43317 (Buysschaert et al. 2018a). In both cases, a negligible amount of CBM was required to best fit the observed pulsation mode frequencies. This led to the tentative conclusion, although it is based on a very small sample size, that the presence of a magnetic field in the near core region is inhibiting mixing at the interface of the convective core and the radiative envelope (Briquet et al. 2012; Buysschaert et al. 2018a). Indeed, a magnetic field in the near-core region could be supported by either a global magnetic field as detected at the surface that extends much deeper within a star and/or a convective core dynamo with magnetic field lines that penetrate outside the convective core (see e.g. Augustson et al. 2016; Augustson and Mathis 2019; Augustson et al. 2020).

The study of pulsating magnetic stars naturally leads to defining the field of magneto-asteroseismology, with an important proof-of-concept of its application being the recent analysis of HD 43317 by Lecoanet et al. (2022). In this work, the best-fitting non-magnetic structure models of the pulsating magnetic star HD 43317 from the forward asteroseismic modelling by Buysschaert et al. (2018a) were used as input into 3D rotating magnetohydrodynamical calculations with the DedalusFootnote 5 code (Burns et al. 2020). Some assumptions inevitably had to be made, which included: (i) the utilised structure models based on non-magnetic evolution models were accurate; (ii) the measured rotation period at the surface from spectropolarimetry was used to assume a rigid interior rotation profile; and (iii) the measured surface magnetic field strength was used as an anchor point to assume a purely dipolar magnetic field geometry within the star’s interior. Given all these assumptions, the 3D magnetohydrodynamical calculations allowed wave-wave interactions to be studied, and specifically the critical magnetic field strength required to suppress a standing gravity wave (i.e. gravity mode) and convert it into an Alfven wave. Based on a comparison of the observed pulsation modes identified by Buysschaert et al. (2018a) in the CoRoT light curve of HD 43317, Lecoanet et al. (2022) were able to place an upper limit of the magnetic field strength of approximately \(5 \times 10^{5}\) G in the near-core region. Such a value is consistent with the expected magnetic field strength in a star with a convective core dynamo supported by a large-scale global magnetic field of fossil origin (Augustson et al. 2016). Therefore, not only pulsation frequencies but also their observed amplitudes provide invaluable insight of the physics of stellar interiors and facilitate magneto-asteroseismology.

5 Upcoming instrumentation and projects

It is not only TESS that is currently driving advances in our understanding of massive star interiors. For example, the European Space Agency (ESA) Gaia mission (Gaia Collaboration et al. 2016) is providing astrometry and light curves of a vast number of variable stars, including potentially hundreds of thousands of massive stars (Gaia Collaboration et al. 2019, 2023). The all-sky Gaia survey provides a truly huge data set for identifying new pulsating massive stars to follow-up. There remains much work to do to fully exploit the natural synergies of Gaia with asteroseismology (e.g. Aerts et al. 2023).

Spectroscopic campaigns unprecedented in terms of size and scope for studying massive stars are also currently underway, which are necessary to supplement forward asteroseismic modelling. One such example is the XShootU collaboration,Footnote 6 which is assembling ground-based European Southern Observatory (ESO) XShooter optical and near-infrared spectroscopy and NASA Hubble Space telescope (HST) UV spectroscopy to provide the largest sample of massive stars at low metallicity to date (Vink et al. 2023). Large programs with ESO instrumentation targeting massive stars in the galaxy and beyond are also assembling extremely useful spectroscopy for placing stars in the HR diagram and extracting their abundances, rotational and macroturbulent broadening velocities (Gebruers et al. 2022; Serebriakova et al. 2023), with the number of targets expected to grow much larger in the coming years. An example of a large sample of high-resolution optical spectra of galactic massive stars is the IACOBFootnote 7 project, which continues to act as a robust data set for studying massive stars (e.g. de Burgos et al. 2023a,b).

5.1 The CubeSpec space mission

A particularly exciting project specifically tailored to empower massive star asteroseismology in the coming decade is the CubeSpec space mission (Bowman et al. 2022), which is being designed and built by KU Leuven, Belgium.Footnote 8 The CubeSpec mission is an in-orbit demonstration of a 6-unit cubesat within the European Space Agency’s GSTP technology programme, and is funded by the Belgian federal science policy office (BELSPO). A schematic of the platform layout of CubeSpec is shown in Fig. 6. The primary engineering goal of the CubeSpec mission is to validate as a proof-of-concept the feasibility of high spectral resolution spectroscopy from a cubesat (Raskin et al. 2018).

Fig. 6
figure 6

Schematic of the platform layout of the CubeSpec mission with important components labelled

Full size image

On the other hand, the primary science objective of CubeSpec is an in-orbit demonstration of massive star asteroseismology using long-term, high-cadence, and high-resolution time-series optical spectroscopy. These data will facilitate the study of spectral LPV caused by pulsations, such that the geometry of pulsation modes in terms of spherical harmonics can be identified (Bowman et al. 2022). Such a method is highly complementary to the methods of pulsation mode identification based on time-series photometry, and had huge success using ground-based telescopes in the past (see e.g. Aerts et al. 1992; Fullerton et al. 1996; Telting and Schrijvers 1998; Uytterhoeven et al. 2004, 2005). However, ground-based telescopes are inevitably hampered by daily aliasing and multi-site coordination to achieve a maximal duty cycle is expensive, which is a negligible issue for the space-based CubeSpec mission.

The original design concept of CubeSpec is discussed in Raskin et al. (2018). Updated engineering and scientific specifications as well as a discussion of expected performance and identification of an initial target list of pulsating massive stars are provided by Bowman et al. (2022). The prospect of combining exquisite time-series spectroscopy from CubeSpec with concomitant time-series photometry from TESS is truly a novel endeavour and will lead to exciting new results for massive stars.

5.2 PLATO

In the not too distant future, the ESA PLATO mission (Rauer et al. 2014) will facilitate major advancements across a range of topics in astrophysics. Selected as ESA’s M3 mission in the 2015-2023 Cosmic Vision program and an expected launch date at the end of 2026, PLATO has the primary science goal to detect and accurately characterise large numbers of exoplanets using the transit method in time-series photometry. To achieve this goal, the PLATO platform consists of 26 white-light cameras each with a field of view of 1037 deg2, which partially overlap to create a combined field of view of 2132 deg2 (Rauer et al. 2014). Such a large area to be observed with high photometric precision means that light curves spanning years will be assembled for hundreds of thousands of pulsating stars. Unlike the Kepler mission, PLATO is not purposefully avoiding massive stars, nor is it avoiding (very) bright stars allowing it to maximise the synergy with ground-based spectrographs, which are typically brightness-limited for the purposes of radial velocity follow-up.

The location of PLATO’s long pointing fields of view on sky will be decided no later than two years before launch. Nevertheless, with such a large area many more massive stars will be included compared to previous missions, such as CoRoT, Kepler and TESS. Moreover, the duration of PLATO’s high-precision light curves in its long pointing fields of view will far exceed the 150-d length of a CoRoT campaign and an individual 28-d TESS sector, and be comparable to the 4-yr length of the Kepler mission. Asteroseismic results based on Kepler mission data of early-type stars have demonstrated that years-long light curves are necessary to resolve individual pulsation mode frequencies of coherent pulsators (Bowman et al. 2016; Bowman and Michielsen 2021). In this regard, PLATO data will have a huge advantage compared to the individual 28-d TESS sectors for many pulsating massive stars in the sky.

With a deluge of high-quality light curves for so many massive stars already available from TESS and expected soon from PLATO, which include not only pulsators but also eclipsing binaries and stars with rotational modulation (etc), the asteroseismic community are already moving from a star-by-star approach to more ensemble-driven analysis techniques. Even the first task of classification of different types of variability in light curves is a computationally expensive and time-consuming task to perform manually. Hence machine-learning techniques are becoming increasingly popular because of their efficiency and overall high success rates (see e.g. Armstrong et al. 2016; Audenaert et al. 2021; Audenaert and Tkachenko 2022; Barbara et al. 2022). Similar efforts to maximise the asteroseismic output of massive stars from PLATO are being prepared as part of an international asteroseismic consortiumFootnote 9 and will be needed once the data begin to arrive.

6 Concluding remarks

In this invited review, I have discussed the exciting results that have been made possible because of space-based time-series photometry missions, paying particular attention to the new results since the launch of the NASA TESS mission (Ricker et al. 2015). In this way, the current review is complementary to the pre-TESS review of massive star asteroseismology by Bowman (2020). The all-sky approach taken by the TESS mission has demonstrated that diverse variability caused by pulsations, binarity, rotation and magnetic fields are widespread across the upper part of the HR diagram (see e.g. Burssens et al. 2020), which has opened up many new avenues for stellar astrophysics research.

An important new discovery continuing to fuel debate in the literature is the near-ubiquitous detection of SLF variability in massive star light curves (Bowman et al. 2019b, 2020). Compelling mechanisms proposed in the literature include gravity waves excited by turbulent core convection (Rogers et al. 2013; Rogers 2015; Rogers and McElwaine 2017; Edelmann et al. 2019; Ratnasingam et al. 2019, 2020, 2023; Horst et al. 2020; Varghese et al. 2023; Herwig et al. 2023; Thompson et al. 2023), the dynamics of turbulent massive star envelopes (Schultz et al. 2022; Jiang 2023), and optically thick line-driven winds in evolved massive stars (Krtička and Feldmeier 2018, 2021). The plausibility of some mechanisms is debated since they depend on the specific numerical setup of 2D or 3D hydrodynamical simulations (Lecoanet et al. 2019, 2021; Le Saux et al. 2023; Anders et al. 2023; Herwig et al. 2023; Thompson et al. 2023). However, it is critical to emphasise that these mechanisms are not mutually exclusive in explaining SLF variability, and that there is growing body of observational evidence for how the morphology of SLF variability in both the time and frequency domains depends on a star’s mass, age, rotation rate and metallicity (Bowman et al. 2020; Bowman and Dorn-Wallenstein 2022).

With high-precision time-series photometry from space missions, asteroseismology has been able to place tight constraints on physical processes that are currently controlled by free parameters in state-of-the-art evolution models. For example, precise masses, radii, ages and interior rotation profiles are available for thousands of low-mass stars (Aerts 2021). However, insight of massive stars is lacking behind in numbers because of their intrinsic rarity but also because of their complexity, which includes their rapid rotation and possible presence of a companion star to name a few. The first forward asteroseismic modelling that appropriately takes sources of observational and theoretical uncertainties into account was recently performed for the \(\beta \) Cep star HD 192575 by Burssens et al. (2023). This work yielded fractional uncertainties of less than 15% for mass, age and core mass of a supernovae progenitor (Burssens et al. 2023), which is unprecedented and thanks to the combination of TESS mission photometry, Gaia astrometry and high-resolution spectroscopy. Moreover, a rigid interior rotation profile is conclusively ruled out for HD 192575, which has already proven an invaluable result for improving 2D stellar evolution models and prescriptions for angular momentum transport (Mombarg et al. 2023).

In addition to constraints on interior rotation and mixing, the emergence of magneto-asteroseismology has allowed the strength and geometry of interior magnetic fields to be probed. For example, Lecoanet et al. (2022) infer a magnetic field strength of \(5 \times 10^{5}\) G in the near-core region of the rapidly rotating SPB star HD 43317 (Buysschaert et al. 2018a). With longer and longer light curves becoming available for pulsating massive stars and larger samples of stars targeted with spectropolarimetry, additional magneto-asteroseismic studies are expected soon.

In the near future, building on the successes of CoRoT, Kepler and now TESS, the field of asteroseismology is moving towards a more holistic approach to understanding the physics of stellar interiors rather than focusing on a niche part of it for individual stars. These efforts will be maximised by ESA’s upcoming PLATO mission providing high-precision light curves for potentially tens of thousands of pulsating massive stars (Rauer et al. 2014). By including complementary observables within a forward asteroseismic modelling framework, such as parallaxes to derive distances from Gaia (Gaia Collaboration et al. 2016), magnetic field strengths from spectropolarimetry, and constraints on the spherical harmonic geometry of additional pulsation modes from spectral LPV with the CubeSpec mission (Bowman et al. 2022), asteroseismology is expected to reach new heights in the decades to come. This is ultimately the goal of the ongoing SYMPHONYFootnote 10 project (P.I. Bowman), which has a holistic approach to studying different aspects of massive stars with asteroseismology as its main tool.