Asteroseismology of high-mass stars: new insights of stellar interiors with space telescopes

Massive stars are important metal factories in the Universe. They have short and energetic lives, and many of them inevitably explode as a supernova and become a neutron star or black hole. In turn, the formation, evolution and explosive deaths of massive stars impact the surrounding interstellar medium and shape the evolution of their host galaxies. Yet the chemical and dynamical evolution of a massive star, including the chemical yield of the ultimate supernova and the remnant mass of the compact object, strongly depend on the interior physics of the progenitor star. We currently lack empirically calibrated prescriptions for various physical processes at work within massive stars, but this is now being remedied by asteroseismology. The study of stellar structure and evolution using stellar oscillations - asteroseismology - has undergone a revolution in the last two decades thanks to high-precision time series photometry from space telescopes. In particular, the long-term light curves provided by the MOST, CoRoT, BRITE, Kepler/K2 and TESS missions provided invaluable data sets in terms of photometric precision, duration and frequency resolution to successfully apply asteroseismology to massive stars and probe their interior physics. The observation and subsequent modelling of stellar pulsations in massive stars has revealed key missing ingredients in stellar structure and evolution models of these stars. Thus asteroseismology has opened a new window into calibrating stellar physics within a highly degenerate part of the Hertzsprung-Russell diagram. In this review, I provide a historical overview of the progress made using ground-based and early space missions, and discuss more recent advances and breakthroughs in our understanding of massive star interiors by means of asteroseismology with modern space telescopes.

stars allow them to be observed at large distances (see e.g. Stark 2016), hence allow us to study the early epochs of the Universe including the re-ionisation of the Universe and the formation of the first galaxies (Bromm and Larson, 2004;Robertson et al., 2010).
Massive stars typically form in dense, cold and large molecular clouds with one of the important signatures of massive star formation being giant filament structures and powerful bi-polar outflows when they are embedded in such dense clouds -see the recent review by Rosen et al. (2020). After the formation phase, massive stars enter the so-called main sequence phase of stellar evolution, which is defined by the onset of hydrogen fusion in the core via the CNO cycle whilst maintaining hydrostatic equilibrium. The length of the main sequence is governed by the nuclear time scale, with more massive stars having shorter main sequence life times. During their lives, massive stars produce intense radiation fields from their high luminosities and experience line-driven winds, which together play major roles in the shaping of their environment (Langer, 2012;Kippenhahn et al., 2012). The majority of massive stars experience an explosive death as a supernova, which provides mechanical and chemical feedback to the interstellar environment (Mac Low and Klessen, 2004;de Rossi et al., 2010;Hopkins et al., 2014;Stark, 2016;Crowther et al., 2016) and can trigger a new generation of stars and planets. The exotic compact remnants of massive star evolution are neutron stars and black holes, with the latter being end products for many of the higher-mass progenitor stars. These remnants facilitate important tests of Einstein's theory of General Relativity and the study of the Universe using gravitational waves when they coalesce (Abbott et al., 2016(Abbott et al., , 2019. Hence, understanding the evolution of massive stars and their roles as supernovae and gravitational-wave progenitors represent fundamental questions in astronomy (Smartt, 2009;Langer, 2012). It is especially important to understand massive star evolution since there is a strong dependence of a supernova's chemical yield and the mass of the remnant on the interior physics of the progenitor star (Nomoto et al., 2006;Langer, 2012;Stark, 2016), and because of the large diversity in observed supernovae light curves (Dessart and Hillier, 2019).
Despite the importance of massive stars in our Universe, their physics is not yet fully understood. There are still many questions spanning all evolutionary phases, and specifically how their formation, evolution and inevitable deaths differ to those of the more common intermediate-and low-mass stars (Langer, 2012). A major shortcoming of current stellar evolution models is that they contain large theoretical uncertainties for massive stars, which is evident already during the earliest phases of stellar evolution including the main sequence. Consequently these uncertainties propagate and strongly impact the post-main sequence stage of stellar evolution (Maeder and Meynet, 2000;Ekström et al., 2012;Chieffi and Limongi, 2013). Hence the power of models for predicting if a massive star will explode as a supernova, the corresponding chemical yield and the mass of the compact remnant are limited by the accuracy of models reproducing the observed properties of stars prior to them exploding as supernovae.
Since massive stars have convective cores and radiative envelopes during the main sequence, the physics and numerical implementation of convection and convective-boundary mixing is crucial in determining their core masses and subsequent evolution (Gabriel et al., 2014;Georgy et al., 2014;Paxton et al., 2018Paxton et al., , 2019. The mixing profile at the interface of convective and radiative regions, and the mixing profile within the envelope directly impact the amount of hydrogen available for nuclear burning. With more internal mixing, a massive star experiences a longer main sequence and produces a larger helium core mass at the end of the main sequence since fresh hydrogen from the envelope is readily supplied to the convective core (Miglio et al., 2008b;Kippenhahn et al., 2012;Pedersen et al., 2018;Michielsen et al., 2019). In the case of single, slowly rotating and non-magnetic massive stars, it is the internal mixing profile and helium core mass that dictates the evolution beyond the main sequence and determines the ultimate end state. However, Rotation plays a major role among massive stars, specifically because rotationally induced mixing is expected within their interiors (Zahn, 1992;Maeder and Meynet, 2000). Yet prescriptions for such mixing profiles are currently assumed in models and controlled by numerous free parameters that have yet to be empirically calibrated. This represents a large source of uncertainty within theoretical evolution models when estimating the masses and ages of high-mass stars (see e.g. Aerts 2020 andSerenelli et al. 2020 for recent detailed reviews). The combined influence of various rotation rates, different metallicity regimes and mass loss through stellar winds also introduce strong degeneracies within evolutionary models of massive stars (Maeder and Meynet, 2000;Georgy et al., 2011;Ekström et al., 2012;Georgy et al., 2013;Chieffi and Limongi, 2013;Groh et al., 2019). Furthermore, based on dedicated studies such as MiMeS, , the BOB campaign (Morel et al., 2015), and the BRITE spectropolarimetric survey (Neiner et al., 2017), approximately 10 % of massive stars are inferred to host a large-scale magnetic field with polar field strengths that range from approximately 100 G up to a few kG. The degeneracies within evolutionary models are also more complex when dealing with the presence of magnetic fields in massive stars (Alecian et al., 2014;Shultz et al., 2018;Keszthelyi et al., 2019Keszthelyi et al., , 2020. It is known that many massive stars are members of multiple systems, which interact over the course of their lifetimes, so theoretical uncertainties are further compounded by the effects of binarity and mass transfer (Podsiadlowski et al., 1992;Sana et al., 2012;de Mink et al., 2013;Duchêne and Kraus, 2013;Moe and Di Stefano, 2017). However, despite the complexities associated with rotation, metallicity, mass loss and magnetic fields, massive stars in multiple systems have proven extremely useful in probing and mitigating model parameter uncertainties (Guinan et al., 2000;Hilditch et al., 2005;Torres et al., 2010;Tkachenko et al., 2014Tkachenko et al., , 2016Almeida et al., 2015;Abdul-Masih et al., 2019;Johnston et al., 2019;Mahy et al., 2020b,a). This is primarily because binary studies have the potential to provide masses and radii from the stars' relative orbital motion around a common centre of mass. Moreover, accurate, absolute and modelindependent masses and radii are achievable in the case of eclipsing binary systems (see e.g. Tkachenko et al. 2020 andSouthworth et al. 2020) owing to the ability to simultaneously model spectroscopic radial velocities together with the eclipse depths in the light curve of a binary system.
To truly maximise the predictive power of evolutionary models for massive stars, it is essential to calibrate their physical prescriptions and parameters using stringent observational constraints on stellar interiors. One of the most successful and novel methodologies for this is called asteroseismology, which uses the resonant oscillation frequencies of stars to probe their structure -see the research monograph by Aerts et al. (2010). Until recently, most observational studies of massive stars have focussed on the determination of global and/or average properties of these stars, such as the effective temperature and surface gravity derived from spectroscopy being used to estimate masses and ages. On the other hand, the recent space photometry revolution has truly brought asteroseismology to the forefront of astronomy as a means to calibrate stellar evolution theory across the Hertzsprung-Russell (HR) diagram (Chaplin and Miglio, 2013;Hekker and Christensen-Dalsgaard, 2017;García and Ballot, 2019;Aerts, 2020).
In this review, I discuss the progress that has been made in constraining the interiors of massive stars by means of asteroseismology and the space photometry revolution made possible thanks to modern space missions. In section 2, I provide a brief overview of asteroseismology, its methodology, and the types of pulsating massive stars it can be applied to. In section 3, an overview of the space telescopes that have led to the space photometry revolution is described. In section 4, I discuss the recent advances in massive star interiors made by means of asteroseismology, and section 5 describes the state-of-the-art of variability studies in some of the most massive stars. I finish by discussing the current challenges and future prospects for asteroseismology of high-mass stars in section 6.

ASTEROSEISMOLOGY
A powerful method for probing and constraining the physics of stellar interiors is asteroseismology, which uses stellar oscillations to probe the physics of stellar structure . The pulsation modes of stars are standing waves exhibiting nodes and anti-nodes and are described by spherical harmonics. In the case of non-rotating and non-magnetic stars, the wavefunctions of stellar pulsations are separable into the radial and angular directions. The radial parts of the wavefunction solutions are characterised by the radial order n. Whereas the angular dependence is characterised by the angular degree (number of surface nodes), and the azimuthal order m (where |m| is the number of surface nodes that are lines of longitude). The simplest example of a pulsation mode is a radial mode for which { , m} = 0 such that the surface of a star expands and contracts during a pulsation cycle. More complex examples of pulsation modes include non-radial modes (i.e. > 0), for which the indices and m define the surface geometry of the oscillation. As an example, the axisymmetric dipole mode (i.e. { , m} = {1, 0}) has the stellar equator as a node. Thus, the northern and southern hemispheres of a star expand and contract in anti-phase with one another.
Although they have a common structure comprising a convective core and a radiative envelope during the main sequence, there are different types of pulsations that can be excited in massive stars. In general, however, the excitation mechanism of pulsation modes has been shown to be the heat-engine mechanism operating in the local maximum of the Rosseland mean opacity caused by iron-group elements -the so-called Z-bump Gautschy and Saio, 1993;Pamyatnykh, 1999;Miglio et al., 2007). This κ-mechanism gives rise to pulsation modes with properties and excitation physics which depend on the host star's mass, age and chemical composition. There are two main types of pulsation modes excited by the κ mechanism in massive stars, which are defined based on their respective restoring force: pressure (p) modes and gravity (g) modes.

Pressure modes
Pressure (p) modes are standing waves for which the pressure force acts as a restoring force . Typically, p modes have high frequencies (i.e. pulsation periods of order several hours in massive stars), can be radial or non-radial and are mostly sensitive to the radiative envelopes of massive stars. For radial p modes, the entire interior of the star acts as a pulsation cavity, with the centre of a star being a node and its surface an anti-node. In the cases of non-radial p modes, the depth of the pulsation cavity from the surface is determined by the local adiabatic sound speed c(r). A non-radial pulsation mode encounters an increasing c(r) when travelling inward from the stellar surface, which causes it to travel faster and be refracted. The depth a non-radial p mode can reach is called its turning radius, r t , which is proportional to ( + 1) and defined outwards from the centre of the star . Thus, higher degree p modes have smaller pulsation cavities that are more sensitive to the stellar surface. The power of asteroseismology is that each pulsation has a cavity defined by a star's structure, such that each pulsation mode can be used as a direct probe of the physical processes at work within its cavity.
If the radial orders of p modes are sufficiently large such that the modes satisfy n , which is called the asymptotic regime, the pulsations are approximately equally-spaced in frequency (Tassoul, 1980). Deviations from a constant frequency spacing are possible, and become more prevalent for more evolved stars. During the main sequence phase of evolution the radius of a massive star increases and the core contracts, which drives the g-and p-mode pulsation cavities closer to one another as a result of an increasing Brunt-Väisälä frequency . Consequently, this can cause a form of pulsation mode interaction called avoided crossings, in which p and g modes can exchange character whilst retaining their identities if their frequencies approach one another (Osaki, 1975). In more evolved cases, such as post main-sequence stars, the evanescent region between the p-and g-mode cavities inside massive stars decreases. This can allow p and g modes to couple with each other and form mixed modes, which are modes with the character of a p mode in the envelope and the character of a g mode in the deep interior . The regularities of asymptotic p modes in the amplitude spectra of low-and intermediate-mass stars has greatly simplified the issue of mode identification and facilitated asteroseismology for low-mass stars (see e.g. Chaplin and Miglio 2013;Hekker and Christensen-Dalsgaard 2017;García and Ballot 2019), but are rarely observed in massive stars (see e.g. Belkacem et al. 2010;Degroote et al. 2010b). Such high-radial order p modes are generally not expected for massive stars owing to the excitation physics of the κ-mechanism being inefficient in driving such modes in massive stars Gautschy and Saio, 1993;Pamyatnykh, 1999;Miglio et al., 2007).
In the presence of rotation the frequency degeneracy of non-radial pulsation modes with respect to m is lifted, which serves as a unique method of mode identification in certain pulsating stars. The simplest case is for stars that rotate (very) slowly and rigidly -i.e. with a uniform interior rotation angular frequency Ωsuch that the splitting of non-radial pulsation frequency, ω n m , is given by where C n is the Ledoux constant which sets the size of the splitting due to the Coriolis force. In this idealised example, the result of Equation 1 produces a multiplet of pulsation frequencies separated by the stellar rotation frequency in the amplitude spectrum for p modes of high radial order or high-angular degree since C n 0 in such cases . An example of rotationally-split quadrupole p modes is shown in Fig. 1, using the example of KIC 11145123 originally discovered by Kurtz et al. (2014). The amplitude spectrum of the resultant quintuplet split by rotation shown in Fig. 1 uses both 1-yr and 4-yr light curves to emphasise the significant improvement in the resolving power of longer light curves for asteroseismic studies of rotation. Therefore, if the rotation rate is sufficiently slow, p-mode multiplets serve as a means of determining the interior rotation rates of stellar envelopes using an almost model-independent methodology.
Beyond the first-order perturbative approach for including the Coriolis force in slow and rigid rotators given in Equation 1, second-and third-order perturbative formalisms have been discussed by, for example, Dziembowski and Goode (1992), Daszyńska-Daszkiewicz et al. (2002) and Suárez et al. (2010). As described by Suárez et al. (2010), it is important to note that the first-order perturbative treatment of the Coriolis force applied to p modes is only applicable for stars with rotation velocities below approximately 15 % of their critical breakup velocity, with faster rotating stars requiring more complex formalisms.

Gravity modes
Gravity (g) modes are standing waves for which buoyancy (i.e. gravity) acts as a restoring force . Typically, g modes have low frequencies, can only be non-radial and are mostly sensitive to the deep interiors of massive stars near their convective cores. In the asymptotic regime, g modes are equally (1 C nl ) env 4 yr 1 yr Figure 1. Example of rotational splitting of quadrupole p modes into a quintuplet using both 1-yr and 4-yr light curves of the star KIC 11145123 (Kurtz et al., 2014). Horizontal red lines correspond to the rotational splitting value of the modes. spaced in period (Tassoul, 1980), and exhibit a characteristic period Π 0 . In the case of a non-rotating and chemically-homogeneous star, Π 0 can be calculated from the individual g-mode periods, P n, , given by in which α is a phase term independent of the mode degree, , and where r 1 and r 2 are the inner and outer boundaries of the g-mode pulsation cavity, and N (r) is its Brunt-Väisälä frequency. Thus, Equation 2 defines a constant spacing in period for g modes of the same angular degree, , and consecutive radial order, n. Equation 3 demonstrates that the characteristic period, Π 0 , is largely determined by the Brunt-Väisälä frequency, N (r), which has a strong dependence on the mass of the convective core, and hence the mass and age of a star (Miglio et al., 2008a).

Interior rotation
Since all massive stars rotate to some extent, the Coriolis force is also a dominant restoring force for gravity modes. Therefore, it is more appropriate to describe massive stars having gravito-inertial modes, for which both the Coriolis force and buoyancy are important. This is particularly true for pulsation modes with frequencies in the co-rotating frame below twice the rotation frequency (see Aerts et al. 2019a). As discussed in detail by Bouabid et al. (2013), the period spacing increases with period in the co-rotating frame for prograde modes, and decreases in the inertial frame. This is because in the co-rotating frame the effective ( + 1) (cf. Equation 2) for prograde sectoral modes decreases with rotation due to the effect of Coriolis force. Whereas in the inertial frame an increasing period spacing is caused by the frequency increase (i.e. period decrease) due to the effect of advection |m|Ω. Therefore, for rotating stars as viewed in the inertial frame by an observer, one expects a decreasing period spacing for prograde g modes and an increasing period spacing for retrograde g modes. Consequently, a powerful diagnostic in interpreting the oscillation spectrum of a rotating star pulsating in g modes is its period spacing pattern, which is defined as the period differences, ∆P , of consecutive radial order (n) gravity modes of the same angular degree ( ) and azimuthal order (m) as a function of the pulsation mode period, P . An example of an observed period spacing pattern for a series of prograde dipole g modes in the star KIC 3459297  is shown in Fig. 2, in which a fit to the g-mode period spacing pattern reveals the near-core rotation rate -see Van Reeth et al. (2016), Ouazzani et al. (2017) and Pápics et al. (2017) for the application of this technique. Under the asymptotic approximation, gravity modes in a non-rotating, chemically homogenous star are equally spaced in period (cf. Equation 2), yet rotation and a chemical gradient left behind from nuclear burning within a receding convective core introduce perturbations in the form of a "tilt" and "dips", respectively (Miglio et al., 2008a;Bouabid et al., 2013). Higher rotation rates induce a larger "tilt" with the gradient being negative for prograde modes and positive for retrograde modes in the inertial frame.
The commonly-used and mathematically appropriate approach to including rotation in the numerical computation of pulsation mode frequencies is the use of the Traditional Approximation for Rotation (TAR; Eckart 1960;Lee and Saio 1987a,b;Bildsten et al. 1996;Townsend 2003b). Such a treatment of the Coriolis force is necessary for g-mode pulsators if 2Ω/ω 1 . This is because high-radial order gravity modes in moderately and rapidly-rotating stars are in the gravito-inertial regime . Within the formalism of the TAR, the horizontal component of the rotation vector is ignored, which is a reasonable assumption for gravito-inertial modes in main sequence stars given that the Lagrangian displacement vector is predominantly horizontal. The differential equations for non-radial pulsations in rotating stars are almost equivalent to those of non-rotating stars (using the Cowling approximation) if ( + 1) is replaced by the eigenvalue of the Laplace tidal equation, λ. Hence, the mathematical framework of the TAR allows the asymptotic approximation to be used for high-radial g modes in rotating stars (see Lee and Saio 1997, Townsend 2003band Townsend 2003a. Today, the TAR has been used to probe the impact of rotation on g-mode period spacing patterns both theoretically (e.g. Bouabid et al. 2013) and observationally (e.g. Van Reeth et al. 2016), and has been extended by Mathis (2009) and Mathis and Prat (2019) to take into account differential rotation and the slight deformation of stars, respectively. The TAR is implemented within the state-of-the-art pulsation code GYRE (Townsend and Teitler, 2013;Townsend et al., 2018) and has been used by various observational studies to probe (differential) rotation inside g-mode pulsators (Van Reeth et al., 2016. We refer the reader to Aerts (2020) for a detailed discussion of the TAR and its application to pulsating stars.

Interior mixing
Since massive stars have convective cores and radiative envelopes during the main sequence, the physics of convection and convective-boundary mixing is crucial in determining their core masses and evolution (Kippenhahn et al., 2012). The mixing profile at the interface of convective and radiative regions, and the mixing profile within the envelope directly impact the amount of hydrogen available for nuclear burning. Mixing at the boundary of convective regions, such as near the convective core in a main sequence star, is typically implemented as "overshooting" in numerical codes and expressed in terms of the local pressure scale height (Freytag et al., 1996;Herwig, 2000). This is predicated on the non-zero inertia of convective bubbles at a convective boundary causing them to overshoot into a radiative layer. In massive stars, the overshooting of the convective core (also known as convective-boundary mixing) entrains hydrogen from the envelope into the core resulting in a longer main sequence lifetime and a larger helium core mass Michielsen et al., 2019). This has a direct impact on the characteristic g-mode period, Π 0 , of main-sequence stars with convective cores (Mombarg et al., 2019).
A non-zero amount of convective core overshooting is necessary when interpreting pulsations in massive stars using 1D stellar evolution codes (Dupret et al., 2004;Briquet et al., 2007;Daszyńska-Daszkiewicz et al., 2013b). Yet, the amount and shape of convective-boundary mixing remains largely unconstrained for such stars. Two examples of typical shapes of convective-boundary mixing profiles currently implemented in evolution codes include a "step overshoot" and an "exponential overshoot" (see e.g. Herwig 2000; Paxton et al. 2015). Typically, these two prescriptions in the shape of convective core overshooting are referred to as α ov and f ov , respectively, in the literature and differ approximately by a factor of 10-12 . However, it is only recently that asteroseismology has demonstrated the potential to discriminate them in observations using g-mode period spacing patterns Pedersen et al., 2018). Moreover, there is considerable ongoing work using 3D hydrodynamical simulations  and g-mode pulsations to probe the temperature gradient within an overshooting layer and ascertain if it is adiabatic, radiative, or intermediate between the two (Michielsen et al., 2019).
In addition to the need for convective-boundary mixing in massive stars, the origin of mixing within their radiative envelopes is also unconstrained within evolutionary models. Direct evidence for needing increased envelope mixing comes from enhanced surface nitrogen abundances in massive stars (Hunter et al., 2009;Brott et al., 2011). Since nitrogen is a by-product of the CNO cycle of nuclear fusion in a massive star (Kippenhahn et al., 2012), an efficient mixing mechanism in the stellar envelope must bring it to the surface. Rotationally-induced mixing has been proposed as a possible mechanism (Maeder and Meynet, 2000), but it is currently unable to explain observed surface nitrogen abundances in slowly-rotating massive stars in the Milky Way and low-metallicity Large Magellanic Cloud (LMC) galaxies (Hunter et al., 2008;Brott et al., 2011). Nor can rotational mixing fully explain surface abundances in massive overcontact systems (Abdul-Masih et al., 2019. Furthermore, there was no statistically-significant relationship between the observed rotation and surface nitrogen abundance in a sample of galactic massive stars studied by Aerts et al. (2014). In fact, the only robust correlation with surface nitrogen abundance in the sample was the dominant pulsation frequency , which suggests that pulsations play a significant role in determining the mixing properties within the interiors of massive stars.

Period spacing patterns as probes of interior rotation and mixing
As previously illustrated in Figure 2, an observed g-mode period spacing pattern provides direct insight of the interior rotation rate of a star. Such patterns also allow the amount of interior mixing in terms of both convective core overshooting and envelope mixing to be determined. Since massive stars have a receding core whilst on the main sequence, they develop a chemical gradient in the near-core region as they evolve. Gravity modes are particularly sensitive to this molecular weight (µ) gradient and in the absence of large amounts of internal mixing this leads to mode trapping (Aerts, 2020). Thus the g modes get trapped which leads to "dips" in the g-mode period spacing pattern on top of the overall "tilt" caused by rotation (cf. Figure 2).
An illustration of the effect of different amounts of interior mixing and rotation on the g-mode period spacing patterns of prograde dipole modes in a 12 M star about halfway through the main sequence is shown in Figure 3. In the left column of Figure 3, the effect of increasing the amount of envelope mixing (denoted by D mix ) is shown going from top to bottom. Whereas in the right column, the effect of increasing the amount of convective core overshooting (denoted by f ov ) is shown from top to bottom. Such values of D mix and f ov represent the range of values typically found by asteroseismic studies of stars with convective cores. Clearly, even for moderate values of envelope mixing (i.e. D mix 10 4 cm 2 s −1 ), the presence of dips in the g-mode period spacing pattern are strongly diminished, as shown in the bottom-left panel of Figure 3. Thus, the observed presence of strong dips in g-mode period spacing patterns already places an upper limit on the amount of envelope mixing possible in such stars. In all panels of Figure 3, three different rotation rates are shown in green, blue and red, demonstrating the significant affect of rotation for g modes, which correspond approximately to rotation frequencies of 0.0, 0.1 and 0.2 d −1 , respectively. The effect of rotation was calculated assuming rigid rotation using the TAR implemented in the GYRE pulsation code (Townsend and Teitler, 2013;Townsend et al., 2018).
Yet, it is possible that only some, or even none, of the pulsation modes shown in Fig. 3 are observed in massive stars, since the excitation of a given pulsation mode depends on stellar parameters, including mass, age, metallicity. The example using a 12 M shown in Fig. 3 does not include any predictions of mode excitation. Nevertheless, Fig. 3 serves as a schematic example of similar behaviour for any main sequence star born with a convective core. In summary, the morphology of an observed gravity-mode period spacing pattern facilitates mode identification and offers a direct measurement of the near-core rotation and chemical mixing within a star. In practice, asteroseismology of g modes requires long-term and high-precision time series (space) photometry to extract g-mode period spacing patterns. These patterns consequently allows one to measure the interior rotation and constrain the envelope mixing of the host star, and the corresponding characteristic g-mode period, Π 0 , which places constraints on its mass and age. From a large and multi-dimensional grid of stellar structure models covering the possible values of mass, age, and interior mixing parameterised by f ov and D mix , a quantitative comparison of observed  Figure 3. Theoretical gravity-mode period spacing patterns for prograde dipole modes of a 12 M star about halfway through the main sequence (i.e. X c = 0.4). The left and right columns are for different envelope mixing (D mix ) in cm 2 s −1 , and exponential convective-boundary mixing (f ov ) expressed in local pressure scale heights, respectively, calculated using the MESA stellar evolution code (Paxton et al., 2019). For each panel, three different rotation rates expressed as a fraction of the critical rotation rate, Ω crit , calculated using the Traditional Approximation for Rotation (TAR) using the GYRE pulsation code (Townsend and Teitler, 2013) are shown. and theoretical gravity-mode period spacing patterns facilitates asteroseismology to derive the interior properties of massive stars .

Instability domains of high-mass stars
The common interior structures of massive stars includes a convective core and radiative envelope, with pulsations in these stars being driven by the κ-mechanism operating within the local opacity enhancements caused by the Z-bump associated with iron-peak elements in their near-surface layers Gautschy and Saio, 1993;Pamyatnykh, 1999;Miglio et al., 2007). The depth of the Z-bump at T 200 000 K depends on the effective temperature of the star and in turn defines an upper and lower temperature boundary for the instability region of the κ-mechanism in the HR diagram, which are sometimes referred to the blue and red edges of an instability region, respectively. In addition to the mass, radius and effective temperature of a star, which define its thermal structure, the metallicity and choice of opacity table in models are also important parameters (Dziembowski and Pamyatnykh, 2008;Paxton et al., 2015;Walczak et al., 2015;Daszyńska-Daszkiewicz et al., 2017). Since the κ-mechanism operates in the Z-bump, it requires a sufficiently-large opacity enhancement to block radiation and excite coherent pulsation modes. This is supported by theoretical models of pulsation excitation, and observations which indicate a dearth of massive pulsators in low-metallicity environments such as the LMC galaxy with Z 0.5 Z (see e.g. Salmon et al. 2012).
Furthermore, rotation plays an important role in defining the instability regions of massive stars. From an observational perspective, moderate and fast rotation distorts the spherical symmetry of a star. This has significant implications for the spectroscopic determination of atmospheric parameters such as the effective temperature and surface gravity, as these parameters are significantly affected by gravity darkening (von Zeipel, 1924;Townsend et al., 2004;Espinosa Lara and Rieutord, 2011). From a more theoretical perspective, the distorted spherical symmetry of rapidly-rotating stars impacts the applicability of using 1D models, and because phenomena associated with rapid rotation such as rotationally-induced mixing, can significantly influence evolutionary tracks in the HR diagram (Maeder and Meynet, 2000;Maeder, 2009;Lovekin, 2020). Moreover, as discussed in section 2.2.1, the Coriolis force perturbs the pulsation frequencies of a rotating star and consequently also the expected parameter range of instability regions in the HR diagram (Townsend, 2005;Bouabid et al., 2013;Szewczuk and Daszyńska-Daszkiewicz, 2017).
The calculation of instability regions for stars requires non-adiabatic calculations, and specifically the calculation of a pulsation mode's growth rate (Unno et al., 1989;Aerts et al., 2010). Following the laws of thermodynamics, heat-driven pulsation modes require that heat is gained in phase with compression during a pulsation cycle. A non-adiabatic pulsation calculation yields a mode's eigenfrequency, and its imaginary component yields the growth rate. For heat-driven modes excited by the κ-mechanism, the growth rate is a positive quantity for modes that are effectively excited and negative for modes that are damped (Unno et al., 1989;Aerts et al., 2010). The instability regions of pulsations in massive stars for different masses, ages, metallicities and rotation rates can be readily calculated for early-type stars by means of stellar structure and evolution codes, such as MESA (Paxton et al., 2011(Paxton et al., , 2013(Paxton et al., , 2015(Paxton et al., , 2018(Paxton et al., , 2019, when coupled to non-adiabatic stellar pulsation codes, such as GYRE (Townsend and Teitler, 2013;Townsend et al., 2018). We refer the reader to Moravveji (2016), Godart et al. (2017) and Szewczuk and Daszyńska-Daszkiewicz (2017) for instability regions of main-sequence massive stars calculated with and without rotation, and to Daszyńska-Daszkiewicz et al. (2013a) and Ostrowski and Daszyńska-Daszkiewicz (2015) for instability calculations for post-main sequence massive stars.
Amongst the early-type stars, there are two main groups of stars that pulsate in coherent pulsation modes excited by the κ-mechanism: the β Cephei (β Cep) stars and the Slowly Pulsating B (SPB) stars, which together span the approximate mass range from 3 to 25 M . Although not traditionally classified as a distinct pulsator group amongst massive stars, the pulsating Be stars are also discussed in this section for completeness. Pulsations in more evolved and/or more massive stars have also been detected, such as those in periodically variable supergiant (PVSG) stars . However, they are not included here since there are currently very few asteroseismic studies of these objects. We refer the reader to Saio et al. (2006) and Ostrowski et al. (2017) for insightful work on the variability of such stars.

β Cephei stars
The β Cephei (β Cep) stars are Population I stars with spectral types ranging from late O to early B on the main sequence, and have birth masses larger than some 8 M and up to 25 M . They pulsate in low-radial order g and p modes excited by the κ-mechanism operating in the Z-bump , and have pulsation periods that range between about 2 and 8 hr (Stankov and Handler, 2005;Pigulski and Pojmański, 2008a,b;Aerts et al., 2010). Most β Cep stars are dwarf stars, making them likely main sequence stars, although a significant fraction are giants or supergiants. Our understanding of the driving of low-radial order g modes in high-mass β Cep stars remains somewhat elusive (Handler et al., , 2017Aerts et al., 2004a), since a substantial overabundance of iron and nickel in the Z-bump is typically needed for the κ-mechanism to be efficient at exciting g modes in such stars (Pamyatnykh et al., 2004;Moravveji, 2016;Daszyńska-Daszkiewicz et al., 2017). Moreover, in addition to the mode excitation by the κ-mechanism, β Cep stars have been shown to exhibit non-linear mode excitation (Degroote et al., 2009) and stochastically-excited pulsation modes (Belkacem et al., 2009(Belkacem et al., , 2010Degroote et al., 2010b), which demonstrate that multiple pulsation driving mechanisms exist in these stars.

Slowly Pulsating B stars
The Slowly Pulsating B (SPB) stars are the lower-mass counterparts of the β Cep stars, with the original definition of this type of variable star made by Waelkens (1991) although individual examples of SPB stars were known previously (e.g. 53 Persei; Smith and McCall 1978). The SPB stars are Population I stars with spectral types that range from B3 to B9 on the main sequence, thus they have birth masses between approximately 3 and 8 M . They pulsate in high-radial order and predominantly prograde dipole g modes excited by the κ-mechanism operating in the Z-bump Gautschy and Saio, 1993) and have pulsation periods that range between a few days and several hours (Waelkens et al., 1998b;Aerts et al., 1999b;De Cat and Aerts, 2002).

Pulsating Be stars
A subset of approximately 20% of non-supergiant massive stars are classified as Be stars (Porter and Rivinius, 2003). This group comprises stars that have shown Balmer lines in emission on at least one occasion, since such emission lines are known to be transient (Zorec and Briot, 1997;Neiner et al., 2011;Rivinius et al., 2013). The Be stars are rapid rotators with circumstellar decretion disks, and their nearcritical rotation rates are thought to be related to their evolutionary history. More specifically, Be stars may have accreted mass from a companion or are the result of a stellar merger, which is supported by the relative rarity of Be stars with main-sequence companions (Bodensteiner et al., 2020). Many Be stars show evidence of pulsations and experience outbursts of material thought to be driven by pulsations Huat et al., 2009;Kurtz et al., 2015). The pulsational behaviour of Be stars is quite diverse, with such stars showing coherent g modes and stochastically excited gravito-interial waves, with the amplitude of their pulsations being connected to whether the star is mid-outburst or in quiescence (Kambe et al., 1993;Porter and Rivinius, 2003;Neiner et al., 2009Neiner et al., , 2012bBaade et al., 2016). Given the diverse variability seen in Be stars, it remains unclear if a single excitation mechanism is unanimously responsible for exciting pulsations in Be stars, with their fast rotation being their main common characteristic (Porter and Rivinius, 2003;Townsend et al., 2004).

THE SPACE PHOTOMETRY REVOLUTION
The successful application of asteroseismology requires long-term, continuous and high-precision time series data to resolve individual pulsation mode frequencies, and perform unambiguous mode identification. As discussed in sections 2.1 and 2.2.3, mode identification using the continuous photometry provided by space telescopes can be readily achieved using rotationally-split modes and period spacing patterns, respectively. However, prior to space telescope missions different techniques were more common. Nevertheless, once mode identification has been achieved, a quantitative comparison of observed pulsation mode frequencies and those predicted by theoretical models in addition to atmospheric constraints such as the effective temperature and metallicity from spectroscopy reveals the physics that best represents the observed star -a fitting process known as forward seismic modelling .

Prior to space telescopes
Prior to the near-continuous photometry provided by space telescopes, pulsation modes in massive stars were identified by means high-resolution and high-cadence spectroscopic and/or photometric time series data assembled using ground-based telescopes (e.g. Smith 1977;Waelkens 1991;De Cat et al. 2005). In the early days of massive star asteroseismology, mode identification via spectroscopic line profile variations (LPVs) targeted the silicon triplet at 4560Å in slowly-rotating β Cep stars (Gies and Kullavanijaya, 1988;Aerts et al., 1992;Aerts and Waelkens, 1993;Aerts et al., 1994a,b;Telting and Schrijvers, 1997;Uytterhoeven et al., 2004Uytterhoeven et al., , 2005, whereas for SPB stars the silicon doublet at 4130Å was also useful (Aerts et al., 1999a;Aerts and De Cat, 2003). A spectral resolving power of at least 50 000 and a very high signal-to-noise being preferable. In fast rotating B stars, such as Be stars, these silicon multiplets can be blended and one must target isolated lines such as the helium I 6678Å line (e.g. Balona and Kambe 1999;Štefl et al. 2003;Maintz et al. 2003). Balmer lines are not suitable for LPV studies as they are less sensitive to the radial and non-radial velocities of pulsations since they are dominated by intrinsic (i.e. Stark) broadening. For a complete discussion of mode identification using spectroscopic time series, we refer the reader to Aerts et al. (2010).
An alternative and complementary methodology to using LPVs in the identification of pulsation modes in massive stars is to use amplitude ratios from multi-colour photometry. Originally devised by Watson (1988) and Heynderickx et al. (1994), this approach was particularly effective when applied to β Cep stars owing to their relatively high amplitude and short-period pulsations. We refer the reader to Shobbrook et al. (2006), Handler et al. (2012) and Handler et al. (2017) for examples of this technique. The analysis of ground-based photometry and/or spectroscopy has been successful in demonstrating the importance of constraining parameters such as convective core overshooting, rotation and metallicity in massive star evolution Handler et al., 2004Handler et al., , 2006Briquet et al., 2007;Daszyńska-Daszkiewicz et al., 2013b;Szewczuk and Daszyńska-Daszkiewicz, 2015).

Early space missions
When the era of space telescopes dawned, the revolution of asteroseismology for stars across the HR diagram truly began. Space telescopes not only provide long-term and near-continuous observations of multiple stars simultaneously, but the typical precision of space-based photometry is at least two orders of magnitude better than what is possible from the ground. Such early space missions, which were not designed for asteroseismic surveys but nonetheless were extremely useful, included the Hipparcos mission (van Leeuwen et al., 1997), which discovered hundreds of pulsating B stars (Waelkens et al., 1998a;Aerts et al., 2006c;Lefèvre et al., 2009). The drastically improved time series data allowed low-amplitude pulsation modes to be detected in stars previously believed to be constant, and yielded precise pulsation frequencies in multi-periodic pulsators that are dominated by complex beating patterns (Koen and Eyer, 2002).
The first space telescope dedicated to asteroseismology was the MOST mission, which was launched in 2003 (Walker et al., 2003). The MOST spacecraft may have been a small telescope, but it detected pulsations in a wide variety of different types of variable stars. Among the massive stars, MOST studied SPB, β Cep and Be stars (Walker et al., 2005;Aerts et al., 2006b,a;Saio et al., 2007b,a;Cameron et al., 2008), including rare instances of β Cep stars in eclipsing binary systems (Desmet et al., 2009c,a). Simultaneous MOST photometry and ground-based spectroscopy proved vital in understanding and modelling the first massive star exhibiting hybrid p-and g-mode behaviour, γ Peg (Handler et al., 2009;Walczak et al., 2013). Furthermore, MOST detected pulsations in blue supergiants, such as Rigel (Saio et al., 2006;Moravveji et al., 2012), and interestingly it also detected pulsations in some Wolf-Rayet stars (Lefèvre et al., 2005;Moffat et al., 2008a) but the absence of pulsations in others (Moffat et al., 2008b). Although the time series photometry assembled by the MOST mission was limited in length, it served as a valuable proof-of-concept exercise demonstrating the power of massive star asteroseismology using space telescopes compared to facilities on the ground.
The next major milestone in the space photometry revolution was the French-led CoRoT mission launched in 2006 ). The CoRoT spacecraft was a combined planet hunting and asteroseismology focussed mission, which delivered short-cadence (i.e. 32 sec) time series photometry for different fields of view across the sky for up to 150 d (Baglin et al., 2009). Over the mission duration, CoRoT contributed a wealth of information concerning variability in massive stars, with some aspects remaining unchallenged in quality and impact more than a decade later. The CoRoT fields of view were optimised to contain hundreds of candidate pulsating B stars (Degroote et al., 2009), which included β Cep, SPB and Be stars (Degroote et al., 2009;Degroote, 2013;Briquet et al., 2011;Pápics et al., 2011;Neiner et al., 2012b;Aerts et al., 2011Aerts et al., , 2019b. The mission was a great success for massive star variability studies, with CoRoT detecting regular frequency spacings in the O8.5 V star HD 46149 (Degroote et al., 2010b), and providing firm confirmation that massive stars can pulsate in both p-and g-mode frequencies . Furthermore, CoRoT led to the first discovery of deviations from a constant period spacing in the B3 V star HD 50230 (Degroote et al., 2010a). In addition to high-precision asteroseismology of targeted β Cep stars , CoRoT also discovered stochastic variability caused by pulsations in massive star photospheres. This includes stochastic non-radial pulsations in B stars (Belkacem et al., 2009;Neiner et al., 2012b) and stochastic low-frequency variability in the three O stars HD 46223, HD 46150, and HD 46966 (Blomme et al., 2011;Bowman et al., 2019a).
Of all the space photometry missions available for asteroseismology, the BRITE-constellation of nanosatellites is ranked amongst the highest for providing excellent asteroseismic returns considering its budget and etendue (Weiss et al., 2014;Pablo et al., 2016). The constellation of nanosatellites from the collaboration of Austria, Poland and Canada, was originally launched in 2013 and have provided long-term time series photometry of some of the brightest stars in the sky. This includes β Cep, SPB and Be stars Baade et al., 2016;Handler et al., 2017;Daszyńska-Daszkiewicz et al., 2017;Walczak et al., 2019), but also pulsating stars in multiple systems (Kallinger et al., 2017;Pablo et al., 2017Pablo et al., , 2019 and stochastic variability in O supergiants and Wolf-Rayet stars (Buysschaert et al., 2017b;Ramiaramanantsoa et al., 2018aRamiaramanantsoa et al., ,b, 2019. The multi-colour and long-term photometry of BRITE and its observing strategy are naturally complementary to combing BRITE data with simultaneous ground-based photometry and/or spectroscopy (e.g. Handler et al. 2017).

The Kepler and K2 missions
Perhaps the most famous of all space telescopes providing time series photometry, the Kepler space telescope was launched in 2009 and had a primary goal of finding Earth-like planets orbiting Sun-like stars Koch et al., 2010). Although by the end of the nominal Kepler mission, the 4-yr light curves of more than 200 000 stars proved invaluable for asteroseismology as well as exoplanet studies. Massive stars (> 8 M ) were purposefully avoided by the Kepler mission since they are bright and typically saturated the CCDs, although dozens of SPB stars and rotationally-variable B stars were discovered (McNamara et al., 2012;Balona et al., 2015). Since the end of the nominal 4-yr Kepler mission, the extremely high precision light curves have been used to discover and analyse rotationally-split modes and period spacing patterns of prograde dipole g modes in dozens of SPB stars . Such asteroseismic studies have revealed the interior rotation rates, convective core overshooting and mixing in main-sequence B stars Szewczuk and Daszyńska-Daszkiewicz, 2018). Despite massive stars not being included in the Kepler field of view, ingenious techniques to extract the light curves of nearby massive stars have been successful (Pope et al., 2016(Pope et al., , 2019a, which includes the scattered-light variability of the O9.5 Iab star HD 188209 . After approximately 4 yr the Kepler spacecraft lost a second vital reaction wheel, which meant that the field of view could not be maintained without a significant expenditure of fuel. A new mission, K2: Kepler's second light, was devised by NASA, which consisted of 80-d campaigns pointing in the direction of the ecliptic (Howell et al., 2014). Since the K2 campaign fields included young star-forming regions, a plethora of massive stars were available to study using K2 light curves. The first proof-of-concept of O-star asteroseismology was demonstrated by Buysschaert et al. (2015), which acted as an important proof-of-concept for massive star variability. The K2 mission allowed the rotation, pulsation and binary properties of the seven sisters of the Pleiades to be studied in detail for the first time thanks to its long time base and high precision (White et al., 2017). Since the end of the K2 mission in 2018, these high-precision data have revealed dozens of previously unknown β Cep and SPB stars (Pope et al., 2016(Pope et al., , 2019bBurssens et al., 2019) and ubiquitous stochastic low-frequency variability in the photospheres of massive stars (Bowman et al., 2019b).

The TESS mission
The ongoing Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) is currently providing high-precision and short cadence (i.e. 2 min) observations for hundreds of thousands of stars across the sky. Each ecliptic hemisphere (|b| > 6 degrees) is divided into 13 sectors which are each observed for up to 28 d. However, there is overlap of the observational sectors near the ecliptic poles, such that stars within the TESS continuous viewing zones (CVZ) have uninterrupted light curves spanning 1 yr. Such an observing strategy is optimised to find transiting exoplanets orbiting bright stars across the sky, but TESS data are also extremely valuable for massive star asteroseismology. Amongst the earliest studies of massive star variability as viewed by TESS are those by Handler et al. (2019), Bowman et al. (2019b) and Pedersen et al. (2019) in which the diverse variability of massive stars is demonstrated using a total sample of more than 200 massive stars. The study by Handler et al. (2019) considers one of the first β Cep stars observed by the TESS mission, which revealed it to be multi-periodic and subsequent asteroseismic modelling indicated it was a runaway star because of its inferred age.
Since the first 13 sectors of TESS data in the southern ecliptic hemisphere have become available, Burssens et al. (2020) have completed a census of massive star variability using TESS data, and coupled this to high-resolution spectroscopy from the IACOB (Simón-Díaz et al., 2011;Simón-Díaz and Herrero, 2014;Simón-Díaz et al., 2015) and OWN (Barbá et al., 2010(Barbá et al., , 2014(Barbá et al., , 2017 surveys. This allowed them to place more than 100 pulsating massive stars in the spectroscopic HR diagram. A summary of the  . Spectroscopic HR diagram of the pulsating massive stars observed by TESS in its sectors 1-13, which have spectroscopic parameters derived by high-resolution spectroscopy, are shown as green circles. Note that the ordinate axis shows spectroscopic luminosity such that L := T 4 eff /g (Langer and Kudritzki, 2014). Red and blue hatched regions denote the theoretical instability regions of g and p modes, respectively, for main sequence stars. Non-rotating evolutionary tracks at solar metallicity (in units of M ) are shown as solid grey lines, and a typical spectroscopic error bar is shown in the bottom-left corner. variability catalogue by Burssens et al. (2020) is shown in the spectroscopic HR diagram in Figure 4, in which the observed stars are shown as green circles overlaid on top of evolutionary tracks (grey lines) and theoretical instability regions for p and g modes (blue and red hatched regions, respectively), which were calculated using the MESA evolution code (Paxton et al., 2019) and the GYRE pulsation code (Townsend and Teitler, 2013). Although there is general agreement between the observed location of pulsating stars and the predicted instability regions, the study by Burssens et al. (2020) clearly demonstrates that our understanding of the variability mechanisms in massive stars is far from complete given the overall distribution of stars within the HR diagram.

EMPIRICAL CONSTRAINTS ON MASSIVE STAR INTERIORS
In the previous section, it was described how early asteroseismic studies of massive star interiors were predominantly based on ground-based campaigns. The acquisition of these necessary data was painstakingly complex and required substantial dedication from all those involved, since it is non-trivial to collect and analyse such fragmented time series. However, despite the complications some aspects of these groundbased studies remain unrivalled to this day owing to the carefully-selected sample of stars and the relative exclusion of massive stars in later space missions. In particular, there remain only a handful of truly massive stars that have been undergone forward seismic modelling. In this section, some important case studies of pulsating massive stars that have studied using asteroseismology are presented. In particular the range of parameters that have been obtained via forward seismic modelling, where available, are summarised in Table 1, but we refer the reader to the individual studies for full details.

Ground-based studies
The first detailed study of the interior of a massive star using asteroseismology is arguably that of the slowly-rotating B3 V star V836 Cen (HD 129929) (Aerts et al., , 2004bDupret et al., 2004). The analysis of the extensive ground-based data set consisting of multicolour Geneva photometry and high-resolution spectroscopy revealed HD 129929 to be a multiperiodic β Cep star pulsating in low-radial order g and p modes (Aerts et al., , 2004b. The rotationally-split multiplets in HD 129929 yielded a non-rigid interior rotation rate where the near-core region was determined to be rotating a factor of ∼4 times faster than the envelope. This was the first measurement of the interior rotation rate of a main-sequence star. Subsequent forward seismic modelling allowed the best-fitting mass, age, metallicity, core hydrogen mass fraction and convective core overshooting to be determined (Dupret et al., 2004), which are provided in Table 1.
An important example of extensive ground-based multicolour photometry and spectroscopy being used to probe the interior of a β Cep star is that of B2 III star ν Eri (HD 29248) (Handler and Aerts, 2002;Handler et al., 2004;Aerts et al., 2004a;Pamyatnykh et al., 2004;De Ridder et al., 2004;Ausseloos et al., 2004;Jerzykiewicz et al., 2005;Suárez et al., 2009;Daszyńska-Daszkiewicz and Walczak, 2010). At the time, and arguably still true today, ν Eri represents one of the best studied β Cep stars owing to the substantial data set having been assembled for the star. The spectroscopic and frequency analyses of this slowly-rotating (v sin i 20 km s −1 ) star revealed rotationally-split low-radial order dipole g and p modes, and possibly high-radial order g modes Aerts et al., 2004a). Forward seismic modelling of ν Eri was unable to find a satisfactory theoretical model unless an iron enhancement throughout the star, and in particular in the Z-bump, was included in the models to explain the mode excitation (Pamyatnykh et al., 2004;Ausseloos et al., 2004). The best-fitting parameters of ν Eri yielded a non-rigid interior rotation rate of approximately 3-5 times faster in the near-core region compared to the envelope, although this is quite uncertain given the asymmetry of the small number of rotationally-split modes available (Pamyatnykh et al., 2004;Ausseloos et al., 2004). Later, Dziembowski and Pamyatnykh (2008), Suárez et al. (2009), andWalczak (2010) revisited the modelling of ν Eri and investigated its interior rotation, overshooting and mode excitation properties based on different opacity tables, specifically focussing on the relative abundance of iron in the Z-bump. In particular, these asteroseismic modelling studies were able to produce a better fit to the observed pulsation mode frequencies using larger overshooting values and an enhancement of iron in the Z-bump, albeit at the expense of not necessarily re-producing the inferred location of the star in the HR diagram. The best-fitting parameters from these studies are provided in Table 1, although the variance in the parameters can be understood as arising from the different metallicities and opacity tables being used.
A similarly famous β Cep star is 12 Lac (B1.5 III; HD 214993), which is comparable in mass to ν Eri although it is somewhat more rapidly rotating with a projected rotational velocity of 30 v sin i 40 km s −1 (Gies and Lambert, 1992;Abt et al., 2002). However, producing a satisfactory asteroseismic model for 12 Lac has been difficult owing to inconsistencies between observed and theoretically-predicted mode identifications, the consequential impact on the inferred interior rotation profile, and stellar opacity data being unable to explain the excitation of all observed pulsation mode frequencies (Aerts, 1996;Dziembowski and Jerzykiewicz, 1999;Handler et al., 2006;Dziembowski and Pamyatnykh, 2008;Desmet et al., 2009b). Despite these difficulties, Dziembowski and Pamyatnykh (2008) conclude based on the then-available observations of 12 Lac that it also exhibits an interior rotation rate of approximately 4-5 times faster in the near-core region compared to the envelope. Desmet et al. (2009b) were able to constrain the parameters of 12 Lac, but the most advanced study of 12 Lac to date was performed by Daszyńska-Daszkiewicz et al. (2013b), who performed mode identification, asteroseismic modelling and studied how rotation, metallicity and opacity data significantly affect the modelling results.
The slowly-rotating (v sin i 30 km s −1 ) β Cep star θ Oph (B2 IV; HD 157056) was also subject to intense photometric campaigns to detect, extract and identify its variability Briquet et al., 2005;Daszyńska-Daszkiewicz and Walczak, 2009). The frequency spectrum of θ Oph resembled that of V836 Cen in that it also contained a single radial mode and a handful of low-radial order rotationally-split multiplets, although one difference is that θ Oph is known to be a multiple system . However, contrary to previous examples of β Cep stars discussed in this section, forward seismic modelling of θ Oph by Briquet et al. (2007) revealed an approximately rigid interior rotation rate and a significantly larger amount of convective core overshooting with the best-fitting model yielding α ov = 0.44 ± 0.07. Such a high value of convective core overshooting is rare in the context of modern asteroseismology of stars with convective cores , but does highlight the trend that stellar models typically underestimate the mass of the convective cores in high-mass main-sequence stars.
A final example case study of forward seismic modelling of a slowly-rotating β Cep star based on groundbased data is the B2 IV/V star V2052 Oph (HD 163472;Briquet et al. 2012). The importance of V2052 Oph in this context is that it is a known magnetic star (Neiner et al., 2003(Neiner et al., , 2012a. A large-scale magnetic field is thought to suppress the near-core mixing caused by convective core overshooting in massive stars. The data set of V2052 Oph assembled and analysed by Briquet et al. (2012) included more than 1300 spectra from 10 different telescopes around the world, and the resultant frequency spectrum of included a radial mode and two prograde non-radial modes identified by means of spectroscopy. Forward seismic modelling of these identified pulsation modes yielded a relatively fast rotation rate with v eq 75 km s −1 , and only a small amount of convective core overshooting with α ov ∈ [0.00, 0.15], which was concluded to be low because of the star's magnetic field .

Space-based studies
Early space missions, such as MOST (Walker et al., 2003) have targeted massive stars for the purposes of asteroseismology. However, in this review we focus on selected case studies from the BRITE, CoRoT and Kepler missions. This is because of their high photometric precision and their long observational base lines Baglin et al., 2009;Borucki et al., 2010;Koch et al., 2010;Weiss et al., 2014), which are necessary for successful forward seismic modelling.
The CoRoT mission provided near-continuous light curves up to 150 d in length at a cadence of 32 s, which at the time was a major revolution for asteroseismology. Early important results from the mission included the seismic modelling of the β Cep star HD 180642 , and detection of the g-mode period spacing pattern in the B3 III SPB star HD 50230 (Degroote et al., 2010a). Such a gmode period spacing pattern was a huge step forward in asteroseismology, as such patterns enable mode identification and allow interior rotation and mixing to be measured directly. Subsequent modelling of HD 50230 revealed it to be a mid-main sequence star with a mass of approximately 7 M , and convective core overshooting of approximately 0.2 ≤ α ov ≤ 0.3 (Degroote et al., 2010a). Later, Degroote et al. (2012) discovered that HD 50230 is a wide-binary system, and confirmed its ultra-slow rotation rate of v eq 10 km s −1 using spectroscopy, although they concluded that this does not impact the aforementioned g-mode period spacing pattern. More recently, Wu and Li (2019) have revisited the analysis of HD 50230 and conclude it to be a metal-rich hybrid pulsator with a modest amount of convective core overshooting. Not long after the early studies by Degroote et al. (2010a) and Aerts et al. (2011), another important result for massive star asteroseismology from the CoRoT mission was made based on the O9 V star HD 46202 . Spectroscopy of HD 46202 confirmed its literature spectral type and yielded a projected surface rotation rate of v sin i 25 km s −1 and an effective temperature of T eff = 34 100±600 K, which makes it amongst the hottest β Cep pulsators known . Several radial and nonradial pulsation modes were identified using the CoRoT photometry of HD 46202. Forward seismic modelling of these identified modes yielded a precise mass and age for this massive pulsator, confirming it as the most massive modelled β Cep star to date. A significant amount of convective core overshooting was required to fit the observed pulsation frequencies and the best-fitting asteroseismic parameters from Briquet et al. (2011) are provided in Table 1.
A third pivotal study of a massive star using asteroseismology from CoRoT data is the magnetic and fast-rotating B3.5 V star HD 43317 Briquet et al., 2013;Buysschaert et al., 2017aBuysschaert et al., , 2018. This star was originally studied by Pápics et al. (2012) using the 150-d CoRoT light curve and it was concluded to exhibit g and p modes based on its frequency spectrum, since it appeared to have independent pulsation modes in both high and low frequency regimes. High-resolution spectroscopy confirmed HD 43317 as an early-type star with an effective temperature of T eff = 17 350 ± 750 K, and as a fast rotator with a projected surface rotational velocity of v sin i = 115 ± 9 km s −1 . Later, HD 43317 was detected to host a large-scale magnetic field with a polar field strength of approximately 1.3 kG (Briquet et al., 2013;Buysschaert et al., 2017a). After revisiting the light curve with the knowledge of how such fast rotation perturbs the frequencies of g-mode pulsations in B stars (see e.g. Kurtz et al. 2015), Buysschaert et al. (2018) extracted all significant pulsation mode frequencies and performed forward seismic modelling of HD 43317 and determined its fundamental parameters, which are provided in Table 1. The small amount of convective core overshooting in the rapidly-rotating and magnetic star HD 43317 led Buysschaert et al. (2018) to a similar conclusion to that of Briquet et al. (2011) for the β Cep star V2052 Oph (cf. 4.1): the presence of a large-scale magnetic field is a plausible cause for the suppression of additional mixing (i.e. overshooting) in the near-core region.
As discussed in section 3.3, during the life time of the Kepler space mission, massive stars were avoided as targets in the field of view. Hence, whereas much of the asteroseismic inference of massive stars based on ground-based data were focussed on β Cep stars, more recent asteroseismic results based on the 4-yr light curves of the Kepler telescope are typically derived from SPB stars . Although not all SPB stars are "massive" stars, they do share a common interior structure whilst on the main sequence, so their discussion is relevant as part of this review.
Amongst the early days of Kepler asteroseismology, two SPB stars were identified as high priority targets: KIC 10526294 and KIC 7760680. A long series of rotationally-split g modes in the B8.3 V star KIC 10526294 allowed Pápics et al. (2014) to determine a near-core rotation period of approximately 188 d. Later, Moravveji et al. (2015) performed the first in-depth asteroseismic modelling of this main-sequence B star, from which it was concluded that a modest amount of exponential diffusive overshooting (i.e. f ov ) fit the Kepler data significantly better than when using the step overshooting prescription (i.e. α ov ). Furthermore, a small but non-negligible amount of additional envelope mixing (D mix 100 cm 2 s −1 ) was needed despite KIC 10526294 being an ultra-slow rotator . The best-fitting asteroseismic modelling parameters, including the measured overshooting of 0.017 ≤ f ov ≤ 0.018, are included in Table 1. The exquisite Kepler data allowed Triana et al. (2015) to compute a near-rigid interior rotation profile through an asteroseismic inversion. Interestingly, at the 1-σ confidence level, the inversion by Triana et al. (2015) confirmed the ultra-slow rotation rate, but revealed that KIC 10526294 had a counter-rotating envelope compared to that of its core. Similar conclusions in terms of a non-negligible amount of convective core overshooting and envelope mixing were also reached for the B8 V SPB star KIC 7760680 based on the period spacing pattern of prograde dipole g modes by Moravveji et al. (2016). However, it should be noted that the uncertainties obtained by Moravveji et al. (2015) and Moravveji et al. (2016) are typically quoted as the step size within the computed grid of evolution models and may be unrealistically small, because they ignore degeneracies amongst the model parameters (see Aerts et al. 2018). There has been substantial development in the correct way to treat parameter correlations and degeneracies for asteroseismology of g-mode pulsators by moving away from the often-used χ 2 merit function and towards the use of the Mahalanobis distance which includes heteroscedasticity. We refer the reader to Aerts et al. (2018) for further details.
At somewhat higher masses, Szewczuk and Daszyńska-Daszkiewicz (2018) performed forward seismic modelling of the SPB star KIC 3240411 using its axisymmetric (m = 0) period spacing pattern extracted from Kepler data. KIC 3240411 was found to be a relatively young star near the zero-age main sequence with the upper limit on its convective core overshooting being f ov ≤ 0.030. The best-fitting model parameters (i.e. model #2 in Table 2 of Szewczuk and Daszyńska-Daszkiewicz 2018) are included in Table 1. Importantly, Szewczuk and Daszyńska-Daszkiewicz (2018) also studied the effect of rotation and opacity data within their modelling for the purposes of explaining mode excitation of high-radial order g modes in such stars. This is because the excitation of such a large number and radial order range for high-radial order g modes remains somewhat difficult to explain from a theoretical perspective.
Recently, the nanosatelittes of the BRITE-constellation mission have been targeting pulsating massive stars with the aim of performing forward seismic modelling of the detected pulsation mode frequencies.
Notable examples of the importance of the contribution of the BRITE mission in this aspect include ν Eri  and θ Oph (Walczak et al., 2019). The BRITE data marked the first time that these famous β Cep stars were observed from space, thus the first time an asteroseismic analysis of them benefitted from having near-continuous, short-cadence and long-term monitoring. Handler et al. (2017) discovered new g-mode pulsations in ν Eri using BRITE data, and Daszyńska-Daszkiewicz et al. (2017) used the extracted pulsation mode frequencies to perform forward seismic modelling. Daszyńska-Daszkiewicz et al. (2017) determined optimum parameters for ν Eri and emphasised the improvement of using modified opacity tables to explain all the observed pulsation mode frequencies. Moreover, the modified opacity data increases the efficiency of convection within the Z-bump, an effect also predicted when increasing the metallicity of the star (Cantiello et al., 2009) or increasing the rotation (Maeder et al., 2008). In the analysis of θ Oph using BRITE data, Walczak et al. (2019) discovered several g modes. From their forward seismic modelling, similar conclusions made by Daszyńska-Daszkiewicz et al. (2017) for ν Eri were made. Specifically that an increase in the mean opacity within the Z-bump is needed to explain all pulsation mode frequencies in θ Oph (Walczak et al., 2019). However, since θ Oph is a triple system with ∼ 8.5 and ∼ 5.5 M primary and tertiary components, respectively, it was difficult to ascertain the origin of the g-mode frequencies (Walczak et al., 2019).
Recently, a new method of performing forward seismic modelling by means of a machine learning was developed by Hendriks and Aerts (2019). The authors trained a deep neural network using more than 62 million pulsation mode frequencies, which were calculated from a vast grid of stellar models covering main sequence intermediate-and high-mass stars spanning from 2 to 20 M . Optimum models were selected using a genetic algorithm making such an automated pipeline extremely quick compared to more conventional forward seismic modelling techniques. Hendriks and Aerts (2019) test their methodology using several well-studied massive stars to benchmark the accuracy of their technique and good agreement is found overall. However, as discussed in detail by Hendriks and Aerts (2019), the modelling results based on their deep neural network depend on the choice of hyperparameters, which affect the ability to find the global minimum in the solution space. In Table 1, the model parameters resulting from the optimised tuning of these hyperparameters are provided (cf. gray points), but such values are claimed to be an "optimal starting point" for further complex modelling. We refer the reader to Hendriks and Aerts (2019) for full details.
Asteroseismology using Kepler data has also been applied to hundreds of intermediate-mass stars covering masses between approximately 1 and 8 M , rotation rates up to 80% of critical, and evolutionary stages spanning from the main sequence through to the red giant branch. More specifically, g-mode period spacings have also been used to probe interior rotation in hundreds of intermediate-mass main-sequence stars, which are more commonly known as γ Dor stars (Van Reeth et al., 2015aOuazzani et al., 2017Ouazzani et al., , 2019Li et al., 2019aLi et al., ,b, 2020). An important conclusion from these works is that current angular momentum transport theory is erroneous by more than an order of magnitude . The situation is less clear for massive stars owing to the much smaller sample size currently available, but significant progress has already been made in recent years because of ground-and space-based data sets and asteroseismology.
The pioneering asteroseismic studies of massive stars using ground-and space-based data sets clearly demonstrate the importance of constraining the interior properties of massive stars and the effectiveness of using g-mode period spacing patterns in SPB stars and low-radial order g and p modes for β Cep stars. Thanks to the ongoing TESS mission (Ricker et al., 2015), the asteroseismic sample of pulsating massive stars has increased in size by at least two orders of magnitude Burssens et al., 2020). Thus detailed forward seismic modelling of many high-mass TESS targets are expected in the not-so-distant future. We refer the reader to Handler et al. (2019) for the first modelling study of a β Cep star using TESS data. Despite the relatively small sample size so far compared to intermediate-mass stars, asteroseismic studies of massive stars have demonstrated the need to include improved prescriptions for convective core overshooting and envelope mixing given that current evolution models underestimate the core masses of massive stars and consequently also their main-sequence lifetimes.

Implications of mixing for the post-main sequence
The majority of asteroseismic studies of massive stars have been of main sequence stars. Of course, based on evolutionary time scales, stars spend more than 90% of their lives in this phase (Kippenhahn et al., 2012). However, pulsating post-main sequence massive stars do exist (Saio et al., 2006;Aerts et al., 2010), and to truly understand their interior physics requires us to first understand the interiors of main-sequence stars. For example, the interplay of different mixing processes can cause some post-main sequence massive stars to undergo blue loops in the HR diagram. However, the exact nature of why some massive stars undergo blue loops and others do not remains unknown (Langer, 2012), especially given the differences in the physics and numerics in various stellar evolution codes (e.g. Martins and Palacios 2013). It is known, however, that convection, mixing and mass loss all play significant roles in determining if stars undergo blue loops in the HR diagram and their pulsational properties (Georgy, 2012;Georgy et al., 2014;Saio et al., 2013;Wagle et al., 2019). To demonstrate how the difference in the amount of the mixing at the boundary of a convective core within a massive main-sequence star drastically impacts its post-main sequence evolution, evolutionary tracks calculated with the MESA stellar evolution code (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019r11554) for initial masses of 12-and 13.5-M and two moderate values for their convective core-overshooting value using the diffusive exponential prescription (i.e. f ov ) are shown in Fig. 5. The evolutionary tracks shown in Fig. 5 use initial hydrogen and metal mass fractions of X = 0.71 and Z = 0.02, respectively, OP opacity tables (Seaton, 2005), and the chemical mixture of Nieva and Przybilla (2012) and Przybilla et al. (2013)  . Hertzsprung-Russell diagram containing non-rotating evolutionary tracks starting from the zero-age main-sequence for stars of initial masses of 12-and 13.5-M calculated using the MESA stellar evolution code (Paxton et al., 2019), using two different values for the diffusive exponential convective core-overshooting value (f ov ), which are expressed in pressure scale heights. With increased mixing in the near-core region, a massive star has a longer main-sequence lifetime and produces a larger convective core mass at the terminal age main-sequence. The increased core mass and interior mixing during the main sequence phase of evolution also impacts the post-main sequence evolution of a massive star, including whether it performs a "blue loop", which is demonstrated by the evolutionary track of the 12-M star with f ov = 0.02.
Currently, there exists a knowledge gap concerning the interiors of blue supergiants yet to be filled by asteroseismology. This is primarily because until recently few pulsating blue supergiants had been found within the Hertzsprung gap in the HR diagram (Bowman et al., 2019b), despite theoretical models predicting such stars should pulsate (Saio et al., 2013;Daszyńska-Daszkiewicz et al., 2013a;Ostrowski and Daszyńska-Daszkiewicz, 2015). Asteroseismology of post-main sequence massive stars represents a major future goal for the community, especially since g-mode period spacing patterns and the presence of radial p modes have the capability to distinguish shell-hydrogen burning and core-helium burning massive stars (Saio et al., 2013;Bowman et al., 2019b). Furthermore, when coupled to spectroscopic surface abundances of core-processed material (e.g. C, N and O), convection and mixing in massive stars can be significantly improved within stellar evolution codes, which has important consequences for stars that eventually explode as a supernova (e.g. Saio et al. 1988;Georgy et al. 2014;Wagle et al. 2019). In particular, the helium core masses of post-main sequence massive stars represent a critical deliverable to the wider astronomical community as it is a fundamental parameter for predicting the chemical properties of supernovae (Smartt, 2009;Langer, 2012).

DIVERSE PHOTOMETRIC VARIABILITY IN MASSIVE STARS
A recent discovery made using the CoRoT, K2 and TESS missions was that the vast majority of early-type stars have significant low-frequency variability in photometry (Bowman et al., 2019a,b). Such stochastic variability is not predicted from the κ-mechanism, but is expected from convectively-driven internal gravity waves (IGWs) excited by core convection (Rogers et al., 2013;Rogers, 2015;Edelmann et al., 2019;Horst et al., 2020). 2D and 3D hydrodynamical simulations predict that IGWs reach the surface with significant amplitudes and provide detectable perturbations in temperature and velocity (Edelmann et al., 2019), and that IGWs are also extremely efficient at transporting angular momentum and chemical elements (Rogers, 2015;Rogers and McElwaine, 2017). A snapshot of a 3D simulation of IGWs propagating within a 3-M main-sequence star from Edelmann et al. (2019) is shown in Figure 6. The morphology of the low-frequency variability in more than 160 massive stars was found to be similar across a large range of masses and ages for both metal-rich galactic and metal-poor LMC stars (Bowman et al., 2019b). The insensitivity of the stochastic variability to a star's metallicity together with the fact that evolutionary timescales predicted most stars were likely on the main sequence was concluded as strong evidence that the observed stochastic variability is likely caused by IGWs excited by core convection (Bowman et al., 2019b). More recently, Bowman et al. (2020) demonstrated that the morphology of the IGWs in the frequency spectrum probes the evolutionary properties of the star, such as mass and radius. Furthermore, the amplitudes of IGWs in photometry were found to correlate with the macroturbulent broadening in the spectral lines of dozens of massive stars observed by the TESS mission , with macroturbulence also having been associated with pulsations (Aerts et al., 2009;Simón-Díaz et al., 2010, 2017. Thus, mixing and angular momentum transport caused by IGWs are important to take into account when studying massive star evolution , and are currently not implemented in most stellar evolution codes. There are currently four known excitation mechanisms that can trigger variability in early-type stars: (i) coherent p and g modes excited by the κ-mechanism Gautschy and Saio, 1993;Szewczuk and Daszyńska-Daszkiewicz, 2017); (ii) IGWs generated by turbulent core convection (Edelmann et al., 2019); (iii) IGWs generated by thin sub-surface convection zones (Cantiello et al., 2009); and (iv) waves generated from tides in binary systems (Fuller, 2017). Whereas essentially all of the asteroseismic results discussed in section 4 are based on heat-driven modes, the other variability mechanisms in massive stars remain somewhat under-utilised despite having great probing power of stellar interiors. In the future, it is expected that binary asteroseismology in particular is going to play a more substantial role in providing excellent constraints of massive star interiors because recent space telescopes such as BRITE and TESS are providing the necessary high-quality time-series observations (see e.g. Jerzykiewicz et al. 2020 andSouthworth et al. 2020).

CONCLUSIONS AND FUTURE PROSPECTS
Our knowledge of the interior physics of massive stars has historically been based on a handful of asteroseismic studies using time series data obtained with ground-based telescopes. Such data sets are notoriously limited by their poor duty cycles which complicates mode identification Dupret et al., 2004;Handler et al., 2004;Pamyatnykh et al., 2004;Ausseloos et al., 2004;Briquet et al., 2007Briquet et al., , 2012Desmet et al., 2009b;Daszyńska-Daszkiewicz et al., 2013b, 2017. These initial asteroseismic studies, however, provided the first evidence of the interior rotation profiles of massive stars and that evolutionary models of such stars lacked calibrated prescriptions for interior processes such as convective core overshooting. Consequently standard evolutionary models typically underestimate the mass of convective cores and hence the main sequence lifetimes of massive stars. Furthermore, ground-based studies yielded the first measurements of the interior rotation rates of main sequence stars , with both rigid and non-rigid rotation rates being inferred for a handful of β Cep stars Pamyatnykh et al., 2004;Ausseloos et al., 2004;Briquet et al., 2007). More recently, the MOST (Walker et al., 2003), CoRoT , Kepler/K2 Howell et al., 2014), and the ongoing BRITE-Constellation (Weiss et al., 2014) missions have advanced asteroseismic studies beyond measuring only the masses, ages and metallicities of massive stars. More advanced forward seismic modelling techniques and improved observations have constrained both macroscopic and microscopic physical processes and mechanisms inside stars, such as the amount and shape of mixing in their radiative envelopes Buysschaert et al., 2018;Szewczuk and Daszyńska-Daszkiewicz, 2018;Walczak et al., 2019) and angular momentum transport as a function of stellar evolution .
Today, thanks to the ongoing TESS mission (Ricker et al., 2015), there is huge asteroseismic potential for massive stars as we are no longer limited by the biases of having a small sample of pulsating massive stars. The long-term and high-photometric precision provided by space telescopes is unrivalled by ground-based telescopes, and the sample of massive stars has expanded to hundreds of stars because of the all-sky TESS mission Burssens et al., 2020). Crucially, TESS is observing massive stars which span a large range in mass and age, but also massive stars in different metallicity regimes. This is because the southern CVZ of TESS includes the LMC galaxy, which allows pulsation excitation models to be tested for metal-rich and metal-poor stars (Bowman et al., 2019b). TESS also offers the opportunity to revisit "old friends" in terms of previously studied massive stars with high-precision photometry, which is particularly useful to probe the long-term stability in pulsation mode amplitudes and frequencies in relatively short-lived stars (e.g. Neilson and Ignace 2015). The diverse variability of massive stars, which includes both coherent pulsation modes excited by the κ-mechanism and IGWs excited by core convection Bowman et al., 2019b;Burssens et al., 2020), enables asteroseismology for a sample of massive stars larger by two orders of magnitude compared to any that came before.
The important future goals of asteroseismic studies based on the large and diverse TESS data set include constraining the helium core masses, near-core and envelope mixing processes, interior rotation profiles and angular momentum transport mechanisms inside massive stars. Insight of the physics in the near-core region of stars above approximately 8 M is currently lacking compared to intermediate-and low-mass stars . Moreover, HD 46202 remains the only massive star above 15 M to have undergone forward seismic modelling . In turn TESS data combined with high-resolution spectroscopy and asteroseismology will mitigate the currently large uncertainties in stellar evolution theory and lead to constraining why only some massive stars undergo blue loops during the post-main sequence phase of evolution (e.g. Bowman et al. 2019b), and improved predictions of supernovae chemical yields and remnant masses. The future is very bright for massive stars, and the goal to calibrate stellar structure and evolution models of massive stars using asteroseismology is now within reach.

CONFLICT OF INTEREST STATEMENT
The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

AUTHOR CONTRIBUTIONS
The author confirms being the sole contributor of this work and has approved it for publication.

ACKNOWLEDGMENTS
DMB is grateful to S. Burssens who graciously re-produced Fig. 4 and to P. V. F. Edelmann who made Fig. 6. DMB is grateful to the anonymous reviewers and Prof. C. Aerts for constructive feedback. DMB is also grateful to the MESA and GYRE developers, in particular B. Paxton and R. H. D. Townsend, for continually supporting the development of state-of-the-art and open-source tools for modelling massive stars. The Kepler, K2 and TESS data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute (STScI), which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support to MAST for these data is provided by the NASA Office of Space Science via grant NAG5-7584 and by other grants and contracts. Funding for the Kepler/K2 mission was provided by NASA's Science Mission Directorate. Funding for the TESS mission is provided by the NASA Explorer Program. Funding for the TESS Asteroseismic Science Operations Centre is provided by the Danish National Research Foundation (Grant agreement no.: DNRF106), ESA PRODEX (PEA 4000119301) and Stellar Astrophysics Centre (SAC) at Aarhus University, Denmark. This research has made use of the SIMBAD database and the VizieR catalog access tool operated at CDS, Strasbourg, France, and the SAO/NASA Astrophysics Data System. The research leading to these results has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 670519: MAMSIE), and a senior postdoctoral fellowship from the Research Foundation Flanders (FWO) with grant agreement No. 1286521N. Table 1. Best-fitting parameters of high-mass stars derived from forward seismic modelling of pulsations discussed in this review. Note that not all studies listed here use the same evolution codes, opacity tables (e.g. standard versus modified), nor necessarily the same numerical setup for selecting the statistically best-fitting asteroseismic model, so we refer the reader to each individual study for full details.