© The Authors 2024

1. Introduction

Macroturbulent line broadening in stellar spectra refers to the broadening of spectral lines caused by large-scale turbulent motions in the stellar atmospheres. Historically, macroturbulent line broadening was identified as a key mechanism decades ago, when it was noticed that rotational broadening alone could not account for the observed widths and shapes of spectral lines. This led to the development of models incorporating large-scale atmospheric turbulence. A comprehensive review on the matter exists thanks to Gray (1978, 2008) and, concerning massive stars, to Howarth (2004).

Particularly strong macroturbulent broadening was found in spectral lines of massive stars. It was characterised by macroturbulent velocities Θ ten times larger than those found in solar-type stars, often dominating over the rotational broadening. Conti & Ebbets (1977) analysed spectra of 205 O-type stars and compared rotational line broadening in main-sequence (MS) and evolved populations. The authors proposed the need for additional turbulent broadening to explain observed differences in distributions. The excess broadening was further quantified by Howarth et al. (1997), who measured the total line broadening in 373 O- and B-type supergiants in a self-consistent way by cross-correlating all observed spectra with one observed template spectrum. They reported the need for additional broadening mechanisms on top of the rotation to be able to explain the observed total line broadening of over 100km s⁻¹ for such a large sample, with presumably randomly distributed inclination angles.

With the development of models for macroturbulent line profiles, it became possible to fit the total shapes of individual line profiles by convolution of intrinsic line profiles with both rotational and macroturbulent profiles. Ryans et al. (2002) analysed high-resolution spectra of 12 B-type supergiants and found that rotational broadening alone was not sufficient to describe observed line profile shapes, and extra macroturbulent broadening of 20–60 km s⁻¹ had to be applied. Even higher macroturbulent velocities of up to 100 km s⁻¹ were found by Simón-Díaz et al. (2010) in a sample of 15 supergiants of late O to B spectral types. The same results were later reported in larger samples, such as the sample of 430 O- and B-type Galactic stars analysed by Simón-Díaz et al. (2017) and the sample of 126 B-type supergiants from the LMC analysed by Serebriakova et al. (2023).

The physical origin of such large velocities present in the photospheres of hot massive stars is still debated. In the Sun, photospheric granulation is directly observed and leads to the macroturbulent broadening of lines in the solar spectrum (and spectra of other solar-type stars). Massive hot stars have radiative envelopes; hence, the macroturbulent broadening in these stars requires another physical explanation. Often it is attributed to the excitation of multiple non-radial pulsations modes Lucy (1976), Kaufer et al. (1997), Ryans et al. (2002), Howarth (2004)). Kaufer et al. (1997) performed a long-term monitoring of six BA supergiants and reported variability in radial velocities and line strengths with multiple-mode cyclic behaviour on timescales of five to 50 days. The authors attributed the detected variability to the effect of non-radial oscillations on the stellar surface. Aerts et al. (2009a,b) modelled spectral line profiles through the collective effect of hundreds of low-amplitude non-radial gravity modes on top of rotation. The authors reported that the resulting macroturbulent profile mimics these collective pulsational effects and that ignoring time-dependent waves leads to single snapshot line profiles with global velocity fields of unrealistic (sometimes supersonic) scale.

While there seems to be agreement that macroturbulent velocities originate from pulsation modes and internal gravity waves (IGWs), the debates are focused on where they get excited: at the interface of the convective core and the radiative envelope (Rogers et al. 2013) or subsurface convective layers. The subsurface convective layers originate in radiative envelopes of massive stars with an initial mass larger than some 10 M_⊙ due to the iron (FeCZ) or helium (HeCZ) opacity bumps. Cantiello et al. (2009) computed models of interiors of stars with masses from 5 to 100 M_⊙ and found correlations between the presence of FeCZ and observed microturbulent velocities ξ and wind clumping. Furthermore, Grassitelli et al. (2015b, 2016) investigated FeCZ and found that areas of high turbulent pressure in these subsurface convective layers, when mapped to the positions in a Hertzsprung-Russell diagram (HRD), aligned nicely with regions where observations showed an excess of stars with high macroturbulent velocity values. This result has been reinforced by the conclusions of Simón-Díaz et al. (2017), who analysed macroturbulent velocities in a large sample of 430 massive stars.

The debates were boosted with the discovery of stochastic low-frequency variability (SLFV) in space photometry of O- and B-type stars. The SLFV manifests itself as a power excess in the low-frequency part of the Fourier spectrum. While the first observational detections of SLFV date back to CoRoT and early Kepler times (e.g. Blomme et al. 2011; Tkachenko et al. 2014; Aerts et al. 2017), detailed observational studies and interpretation were first presented in Bowman et al. (2019a,b). Furthermore, Burssens et al. (2020) and Bowman et al. (2020) find a correlation between SLFV properties and macroturbulent broadening, a finding suggestive of a connection between the two observed phenomena, which are both caused by the modes excited near the boundary of the convective core.

Lecoanet et al. (2019, 2021) provide a theoretical model for IGWs and conclude that the calculated surface manifestation does not correlate with the observed SLFV nor macroturbulence. At the same time, Edelmann et al. (2019) and Vanon et al. (2023) perform 3D simulations of core convection for 3 and 7 M_⊙ stellar models, respectively, and report theoretical frequency spectra to have properties similar to those of the observed SLFV. In particular, the authors note a pronounced low-frequency excess present in the models, with the frequency range decreasing with stellar age. Furthermore, Ratnasingam et al. (2020) and Ratnasingam et al. (2023) perform 2D simulations and compute frequency spectra of IGWs from 3, 5, 7, 10, and 13 M_⊙ stellar models. The authors report the slopes of the frequency spectra to be in agreement with observed properties of SLFV. Anders et al. (2023) challenged higher masses and performed comprehensive simulations of core convection for 3, 15, and 40 M_⊙ stellar models at the zero-age MS (ZAMS). Based on their simulations, the authors predict the frequency spectra of the gravity waves propagate to the surface. Contrary to the previous findings based on 2D and 3D simulations, Anders et al. (2023) conclude that all three models predict variability amplitudes several orders of magnitude lower than both observed SLFV amplitudes and current observational detection capabilities.

Apart from IGWs, the subsurface convective layers were also simulated and tested as a plausible explanation of SLFV and macroturbulent broadening. Cantiello et al. (2021) find that model properties of FeCZ from 1D MESA simulations in a large range of masses and ages correlate well with both SLFV and macroturbulence. A similar conclusion is presented in Ma et al. (2024). Three-dimensional simulations of exteriors are computed by Schultz et al. (2022) and Schultz et al. (2023a) for a 13 M_⊙ star at the terminal-age MS (TAMS), a 35 M_⊙ star at ZAMS and at mid-MS, and 56 and 80 M_⊙ stars in the Hertzsprung gap. The authors conclude that FeCZ and HeCZ extend significantly further in the radial direction compared to 1D models and reach the stellar photosphere. The resulting large turbulent velocities at the stellar surface are found to be in agreement with similarly large observed macroturbulent velocities. Moreover, the same models were used by Schultz et al. (2023b) to synthesise spectral line profiles broadened by the clumpy turbulent surface. The authors conclude on the similarity between the shapes of the macroturbulent broadening-dominated profile and their simulated line profiles.

In this study, we present a statistical analysis of a sample of 594 stars for which observational inference of the macroturbulent velocity is available. The sample covers masses from 2.5 to 80 M_⊙ and two metallicity regimes: the Large Magellanic Cloud (LMC) and the Galaxy. For 86 stars, we present macroturbulent velocities that are measured for the first time. In addition to the observed sample, we compute a dense grid of MESA stellar structure and evolution models and investigate the properties of FeCZ and HeCZ across the HRD. We compare the statistical properties of the observed sample and model properties of FeCZ in order to examine their possible connection. Moreover, for the full grid of models, we estimate the theoretical possibility of IGWs tunnelling through convective layers and assess its possible connection to the observations.

A full description of the observations and analyses involved in the sample is given in Sect. 2. Section 3 summarises input physics included in the computation of the MESA models. Section 4 is dedicated to qualitative and quantitative analyses of the observed macroturbulent velocities and model convective velocities, delving into multivariate linear regression in Sect. 4.1 and principal component analysis in Sect. 4.2. Finally, Sect. 5 explores the propagation properties of IGWs. The conclusions of our study are presented in Sect. 6.

2. Observed sample, data reduction, and spectroscopic parameter determination

Our study focuses on a comprehensive sample of stars covering a mass range from 2.5 to 80 M_⊙ with measured macroturbulent velocities Θ. The sample consists of data from two ESO Large Programs using the UVES and FEROS spectrographs: 126 supergiants in the LMC from Serebriakova et al. (2023), 86 Galactic MS stars from Gebruers et al. (2022), and 382 Galactic MS and supergiant stars from Simón-Díaz et al. (2017). This results in a total of 594 stars, covering metallicity regimes of the Milky Way (Solar metallicity of Z = 0.014 is adopted for this regime) and LMC (Z = 0.006) galaxies. The full sample is shown on the HRD in Fig. 1. We note that we plot the so-called spectroscopic luminosity defined by ℒ ≡ T_eff⁴/g for the observed stars in all figures, following Simón-Díaz et al. (2017). This spectroscopic luminosity can be inferred directly from observed T_eff and log g values without the requirement of knowing the distance and extinction, while it behaves similarly to the actual luminosity of the star (Langer & Kudritzki 2014) and hence can be used as a proxy for it. The use of the observed spectroscopic luminosity in the spectroscopic HRD may introduce shifts of individual objects compared to their position in the classical HRD, thus affecting evolutionary mass determination. This is particularly the case when uncertainties in log g are high (as seen in supergiants; see Markova et al. 2018; Holgado et al. 2020). However, it remains a valid tool for statistical analyses of large stellar samples as long as no conclusions are drawn about the individual stellar masses.

Fig. 1. Full observed sample of high-mass stars with measured macroturbulent velocities. Left panel: LMC supergiants from Serebriakova et al. (2023) over evolutionary tracks for Z = 0.006. Middle and right panels: Galactic objects from Gebruers et al. (2022) and Simón-Díaz et al. (2017), respectively, over evolutionary tracks for Z = 0.014.

The data reduction, atmospheric parameter determination, and line broadening measurement for the three subsamples were handled primarily within their respective studies. Our work involved using the processed data and applying additional steps where necessary, specifically for the Galactic MS stars from Gebruers et al. (2022) where macroturbulent velocity measurements were not included in the analysis. The data reduction and processing varied slightly across the different subsamples due to the distinct methodologies employed in each of the referenced studies. Here, we provide a brief summary of the analyses employed for each subsample.

2.1. Subsample of LMC supergiants from Serebriakova et al. (2023)

The observations were conducted as part of the large programs using the Ultraviolet and Visual Echelle Spectrograph (UVES; Dekker et al. 2000) mounted on UT2 at the VLT and the Fiber-fed Extended Range Optical Spectrograph (FEROS; Kaufer et al. 1999) mounted on the MPG/ESO 2.2-m telescope. These programs aimed to deliver high-resolution, high signal-to-noise ratio spectroscopy of a sample of apparently single OB stars in the LMC. Each target was observed at two epochs separated by at least two to three weeks to detect potential spectral variability. The data reduction process employed the ESOREFLEX data reduction pipeline for the UVES spectra and the CERES (Brahm et al. 2017) pipeline for the FEROS spectra. Custom order-by-order normalisation was applied to correct for artefacts resulting from the order merging of the original pipelines.

The atmospheric parameters were determined using a grid search method involving the TLUSTY (Hubeny 1988), CMFGEN (Hillier & Lanz 2001), and GSSP (Tkachenko 2015) model spectra, depending on the parameter range. The projected rotational velocities v sin i were inferred using the Fourier transform (FT) method. This technique identifies the position of the first zero in the FT of a rotationally broadened profile, directly translating to the v sin i of the star. Specific spectral lines such as Si II, Si III, Mg II, and He I were selected based on their suitability for different types of stars in the sample. Macroturbulent velocities Θ were deduced by fitting synthetic spectra to observed line profiles of the Mg II 4481 Å line. The synthetic intrinsic profiles in the form of specific intensities I(μ) were computed for every object’s given atmospheric parameters and convolved with a radial-tangential macroturbulent kernel. With the integration of I(μ) over the sphere, rotational broadening naturally included realistic limb darkening for the exact atmosphere model of each object. The projected rotational velocities were allowed to vary within 1σ uncertainties of the values inferred with the FT method.

2.2. Subsample of Galactic main-sequence stars from Gebruers et al. (2022)

The observations were conducted as part of the same FEROS large program as in Serebriakova et al. (2023) albeit focused on Galactic MS OB-type stars. The modified CERES pipeline was employed for the extraction of stellar spectra. In the original study, the normalisation function was optimised together with atmospheric parameters. For the sake of consistency with Serebriakova et al. (2023), we employed the same order-by-order approach to normalise the spectra. Atmospheric parameters were derived using the ZETA-PAYNE (Straumit et al. 2022) pipeline, which employs a pre-trained neural network to predict normalised synthetic spectra for a given set of stellar labels (T_eff, log g, v sin i, [M/H]).

Since v sin i derived by ZETA-PAYNE represents the total broadening due to the collective effect of rotation and macroturbulence (as the latter was not included in the fitting), we additionally employed the same method as in Serebriakova et al. (2023) to infer v sin i values, namely, using the FT method to separate the rotational broadening. The original study did not include measurements of Θ; hence we repeated the same procedure as in Serebriakova et al. (2023) to infer Θ, which implied fitting the Mg II 4481 Å line profile with both rotationally and macroturbulent broadened synthetic intrinsic profiles.

2.3. Subsample of Galactic main-sequence and supergiant stars from Simón-Díaz et al. (2017)

The sample was observed as part of the IACOB project (Simón-Díaz et al. 2011), utilising high-resolution spectra from the Fiber-fed Echelle Spectrograph (FIES; Telting et al. 2014) and the High-Efficiency and high-Resolution Mercator Echelle Spectrograph (HERMES; Raskin et al. 2011) attached to the 2.56-m Nordic Optical Telescope and 1.2-m Mercator Telescope, respectively. The spectra were reduced using the FIESTOOL and HERMESDRS pipelines for FIES and HERMES, respectively, followed by normalisation using custom procedures implemented in IDL.

The atmospheric parameters were obtained using the IACOB-GBAT tool with a grid of FASTWIND (Puls et al. 2005) synthetic spectra. Projected rotational velocities were determined using the IACOB-BROAD tool, which combines the FT and goodness-of-fit (GOF) methodologies. The same IACOB-BROAD tool was used to fit radial-tangential macroturbulent and rotational profiles. The FT method provided initial estimates, while the GOF method was used to refine the final values for both broadenings.

2.4. A look into the full combined sample of the observed macroturbulent velocities Θ

The macroturbulent velocities Θ of the full sample are shown as histograms and kernel density estimates (KDE) in Fig. 2. The left panel shows the entire sample. The middle panel shows two subsamples separated by metallicity (dashed light grey for the LMC and dark grey for the Galaxy). The right panel shows two subsamples separated by the method of measurement (dashed light grey for the method we applied and dark grey for the one used in Simón-Díaz et al. 2017), which we treat separately to analyse potential systematic shifts due to methodological differences. One can see three peaks in the distributions around 20 km s⁻¹, 40 km s⁻¹, and 60–80 km s⁻¹. Their meaning becomes clear when comparing Θ to the other observables such as temperature and luminosity in Fig. 3, where colours show the temperature, the symbol size shows the spectroscopic luminosity, and the symbol itself (circles and triangles) shows whether the measurement originates from us or is taken from Simón-Díaz et al. (2017), respectively. The first and the last peaks of the histograms are only visible in the Galactic subsample and belong to the MS and evolved populations, respectively. This is seen in Fig. 3, with the bulk of low-Θ objects corresponding to low-luminosity and high-log g at all temperatures. The highest Galactic peak of 60–80 km s⁻¹ corresponds to the evolved population with low-log g and high luminosity at all temperatures. The middle 40 km s⁻¹ peak originates clearly from the LMC subsample and is analogous to the peak of 60–80 km s⁻¹ shifted to lower values due to the difference in the metallicity of the subsamples. There is no second low peak in the LMC subsample since all the observed targets from the LMC fell in the supergiant regime (see left panel in Fig. 1). The peak at 80 km s⁻¹, while seen in the KDE as a separate feature, belongs to the same group of galactic supergiants as can be seen in Fig. 3.

Fig. 2. Histograms and KDE (solid line) of the measured Θ. From left to right: Full sample of 594 stars, Galactic (dark grey) versus LMC (light grey) stars, and stars for which macroturbulent velocities were inferred with the methods of Serebriakova et al. (2023, light grey) and Simón-Díaz et al. (2017, dark grey).

Fig. 3. Macroturbulent velocity as a function of stellar parameters. From top left to bottom right: Θ versus v sin i, Θ versus Teff, Θ versus log g, and Θ versus log(ℒ/ℒ⊙). The marker size is proportional to ℒ, and the marker colour shows temperature Teff. The error bars in the lower right panel are only shown for our measurements as they are not available for the sample of Simón-Díaz et al. (2017).

Another notable feature is visible in the KDE in the right panel of Fig. 2. A small shift of the low-Θ peak occurs between our and Simón-Díaz et al. (2017) measurements of the Galactic MS stars. The same difference may be noticed in the lower right panel of Fig. 3 (showing a pairwise plot of Θ and spectroscopic luminosity) as a slight kink present around 2–2.5 in log(ℒ/ℒ_⊙). It is unclear whether the shift is caused by methodological differences or is due to differences in the parameter region covered by our samples. Our Galactic sample covers masses from 2.5 to 5 M_⊙, while the sample of Simón-Díaz et al. (2017) starts at 4 M_⊙; hence the overlap is limited to a few stars only.

The upper-left panel of Fig. 3 shows relations between the projected rotational velocity v sin i and macroturbulent velocity Θ. One can see qualitatively three populations showing a linear growth of Θ with v sin i, as noticed by Simón-Díaz et al. (2017), with different slopes for each population. The clear linear trends are likely caused by the ability of the FT method to separate rotational and macroturbulent broadenings, which strongly depends on whether they operate in different frequency ranges of the FT or not. The limitations of the FT approach for measuring rotational velocities are extensively discussed by Aerts et al. (2014) and Serebriakova et al. (2023); thus we restrain from the detailed discussion of this matter. We note though that the FT method relies on several assumptions. First, it assumes that the total broadening may be represented by a convolution of the rotational profile with all other mechanisms as independent functions. Second, the rotational profile itself is thought to be a simplified function with a fixed linear limb darkening coefficient, which allows for the FT of the rotational profile of a given v sin i to be characterised by a single parameter – the position of the first zero. Moreover, the position of the first zero in the FT should ideally not overlap features caused by other broadening functions, which implies rotation dominating over all other mechanisms. Lastly, the analysed line profiles are assumed to be symmetric. In reality, these assumptions are rarely met and the FT of the line profiles shows several zeros that are smeared out and are thus difficult to interpret. This causes, in our case, high degeneracy between Θ and v sin i, and leads to the large errors for Θ in the lower right panel of Fig. 3.

In all pairwise plots of Fig. 3, the macroturbulent velocities Θ show the same trends of increasing with temperature and luminosity and decreasing with surface gravity. The Θ of Galactic stars are systematically higher than those of the LMC stars, although it is not clear whether the difference is caused by the metallicity and not the bias of the LMC sample that is overpopulated with evolved stars.

3. MESA models

A grid of stellar structure and evolution models was computed with the MESA code (Paxton et al. 2011, 2013, 2015; Paxton et al. 2018, 2019; Jermyn et al. 2023) of version r15140. The chosen input physics is taken from Michielsen et al. (2023), where a detailed description is provided. For convection, we used the mixing length theory (MLT) with a mixing length parameter α_MLT = 2.0. The Ledoux criterion for convection was adopted with no semi-convection included. Michielsen et al. (2023) includes extra core-boundary mixing, implemented as both a step-function increase of the core size by a factor α_ov (often referred to as step overshooting parameter) and by the exponentially decreasing mixing characterised by the f_ov parameter. We set these parameters to 0.05 and 0.005, respectively, to keep the extra mixing minimal.

The grid was computed for the full mass range in our sample and uniform coverage in the HRD (2, 2.5, 3, 4, 5, 7, 9, 12, 15, 20, 25, 32, 40, 60, and 80 M_⊙). The evolution was stopped at log T_eff = 3.8 for masses up to 9 M_⊙ and log T_eff = 4.0 for higher masses, which is dictated by the span of our observed parameter ranges.

Figure 4 shows Kippenhahn diagrams for initial masses 60, 20, 12, 7, and 3 M_⊙ with Z = 0.014 (while Z = 0.006 are available in the Appendix A). The radiative transfer zones are hatched, and Ledoux-unstable convective zones are colour-coded based on convective velocity. The odd panels show the upper 10% of the stellar radius, where one can see the development of subsurface convective layers with age and initial mass, as well as ranges of the convective velocity. This velocity range is near-zero for lower masses and only reaches some 10 km s⁻¹ late in the post-MS stage, while high-mass stars have up to 80 km s⁻¹ already at the early MS.

Fig. 4. Kippenhahn diagrams for 60, 20, 12, 7, and 3 M⊙ (rows) stellar models at metallicity Z = 0.014. Odd columns show the full radius (in mass units), and even columns show the zoom into the outer 10% of the radius (in radius units). Hatched areas denote radiative zones. The Ledoux-unstable convective zones are mapped by colour based on convective velocity. The two left columns are plotted with the age in megayears, while the two right columns are plotted with the age in the model number, which allows for a better view of the short-lived post-MS stages.

These models are shown in more detail in complementary figures at Zenodo (see the Data availability section for links) as Appendices C and D, where full internal profiles are plotted for three snapshots at three evolutionary stages: mid-MS, TAMS, and HG. While Fig. C.1 shows the convective velocity and opacity versus temperature, Figs. C.2–C.4 show the same quantities as a function of the stellar radius. In Fig. C.1, one can clearly see the iron and helium opacity bumps depending on the mass and evolutionary stage, triggering the formation of a convective zone. The iron opacity bump at T ∼ 2 × 10⁵ K is responsible for the high convective velocities in the high-mass regime (7 M_⊙ and higher). The He II opacity bump at T ∼ 4 × 10⁴ K causes convective velocities to occur at lower masses. This opacity bump also becomes significant at the late stages of evolution for all mass regimes in the expanded envelope, where it forms the second peak in the convective velocities with a much smaller amplitude than the velocities due to the iron peak. Despite the absolute values of the FeCZ convective velocities being large in some cases, we show that they are subsonic by plotting velocities as a fraction of the sound speed in Figs. C.4 and C.8.

4. Comparing convective velocity and macroturbulent velocity

Fig. 5 shows the distribution of the convective velocities as predicted by MESA models in the HRD, separately for the Fe and He opacity bumps. Measurements of the macroturbulent velocity are presented with the coloured symbols (see figure caption for details). It is well visible that the convective velocities in the Fe opacity bump have very similar trends as the observed macroturbulent velocities in the high-mass regimes (M > 10 M_⊙ or log(ℒ/ℒ_⊙) > 2.5) with the same ranges of variation. For the lower masses, the observed macroturbulent velocities have much larger scatter and amplitudes than the convective velocities in both the Fe and He opacity bumps. The same is seen in Fig. 6 where the same macroturbulent and interpolated convective velocities are plotted along the x-axis vs spectroscopic luminosity log(ℒ/ℒ_⊙) > 2.5, while colour maps the effective temperature. The He-bump velocity is too small to explain any observed macroturbulent velocities that are two orders of magnitudes larger. The Fe-bump only appears for stars with log(ℒ/ℒ_⊙) > 2.5 or M > 10 M_⊙. Therefore, we hypothesise that the Fe-bump and macroturbulent velocities are connected in the high-mass regime while there is a lack of such relation for lower stellar masses.

Fig. 5. Convective velocities in the HRD interpolated from MESA models for continuous representation. Top to bottom: Metallicity Z = 0.006 and 0.014. Left to right: Convective velocities in the regions of iron and helium opacity bumps. Evolutionary tracks that were used for the interpolation are shown in each panel. Coloured circles show observed macroturbulent velocities measured in this work, triangles those measured in Simón-Díaz et al. (2017).

Before embarking upon a comparison between the macroturbulence and the convective velocities, we point out that Θ measured from spectral lines cannot be directly fitted by the convective velocities in the subsurface convection layers of 1D models, since those layers are located much deeper in the stellar interior than the line formation region in the photosphere. This is a limitation of 1D evolution models based on mixing-length theory. Schultz et al. (2022, 2023a) showed that subsurface velocities connected to turbulence do reach the photosphere in their 3D hydrodynamical simulations, but these only cover dynamical timescales and cannot be used for evolutionary studies. Moreover, 2D unified box-in-a-star simulations by Debnath et al. (2024) reveal a highly structured and turbulent subsurface atmosphere causing turbulent velocities of order 30 to 100 km s⁻¹ in the surface layers of galactic O stars. Their simulations give rise to turbulent gas interacting with a line-driven wind above the photosphere and result in an outflow of gas rather than the parcels turning over and returning back into the stellar envelope. This interplay between the subsurface turbulent pressure and the stellar wind implies variations in the flux and creates turbulent velocities in the line-forming regions that scale linearly with the convective velocity obtained from MLT (Moens et al., in prep.). For this reason, it is a good strategy to analyse the possible connection between Θ and the convective velocities from 1D models based on MLT. We thus study their trends from a statistical point of view for a large sample of stars.

4.1. Multivariate linear regression

As the first step of this analysis, we follow the statistical approach of Aerts et al. (2023) and Serebriakova et al. (2023) and employ multivariate linear regression. Multivariate linear regression is a statistical technique used to model the relationship between two or more predictors (independent variables) and a single predicted (dependent) variable. It extends the concept of simple univariate linear regression, which deals with one predictor and one predicted variable, to accommodate multiple predictors, thereby allowing for a more comprehensive analysis of complex data sets. Multivariate linear regression can be employed to investigate the relationships between stellar parameters and show how various factors simultaneously influence a variable of interest, which in our case is macroturbulent velocity Θ. We investigate trends in the observed Θ quantitatively from multivariate linear regression models in the following form:

$$\begin{array}{c}\hfill V=c_0+{\displaystyle \sum _{i=0}^N}{c}_i{P}_i,\end{array}$$(1)

where V is a dependent variable (in our case, macroturbulent velocity Θ), P_i – predictors (other observables or model parameters), and c_i – coefficients of linear regression for a corresponding i-th predictor.

In Serebriakova et al. (2023), based on adjusted R²-statistics, 96.4% of variability in the observed macroturbulent velocities of 126 LMC stars were found to be explained by only two predictors: T_eff and v sin i, while the set of observables (T_eff, log g, log(ℒ/ℒ_⊙), and v sin i) was tested and the rest observables were found insignificant for predicting Θ. Due to the absence of any analogue of v sin i in our non-rotating MESA models, we need to repeat the test of different combinations of predictors without v sin i to find a set of significant observables that are available in models. In absence of v sin i, we found both log g and log(ℒ/ℒ_⊙), in addition to log T_eff, to be significant, although log(ℒ/ℒ_⊙) is computed from ℒ = T_eff⁴/g. Therefore, we further test all three parameters, log(ℒ/ℒ_⊙), log T_eff, and log g, as predictors on the observed data. The results of the test are presented in Table 1 as p-values and coefficients c_i per each predictor, and adjusted R² per tested dataset.

We repeated the exercise for various subsamples with cut-off by log(ℒ/ℒ_⊙) in order to look at behaviour of the measured microturbulent velocities in different luminosity regimes. These regimes are evident in Fig. 6 where a turn of Θ is seen towards higher values above log(ℒ/ℒ_⊙) = 2.5. For the LMC sample, we only use the full sample as all stars in that sample already have log(ℒ/ℒ_⊙) > 3. In Table 1, one can see how the fraction of explained variance (reflected in adjusted R²) increases for subsamples with higher log(ℒ/ℒ_⊙), and reaches 60%. Moreover, the samples that only include stars with log(ℒ/ℒ_⊙) < 3 show large p-values and near-zero R², suggesting the predictors fail to explain variability in Θ for low-log(ℒ/ℒ_⊙) subsample.

Fig. 6. Dependence of the macroturbulent and convective velocities on spectroscopic luminosity of the star. From left to right: Observed macroturbulent velocities (left) and convective velocities for the iron (middle) and helium (right) opacity bumps predicted with MESA. The convective velocities in the middle and right panels have been interpolated on the observed positions in the HRD from Fig. 5. Symbol size and shape are the same as in Fig. 3.

Furthermore, we can compare the roles of each predictor in explaining observed Θ with the roles of the same predictors in explaining model convective velocities v_conv by looking at their coefficients. We repeated the same multivariate linear regression for the model sample, using each computed internal profile as a point of a sample with model v_conv as the dependent variable. First, we repeated tests of sets of predictors for their predicting power (which are not limited by a modest set of observables that we have to stick to with observed data). While it appeared possible to describe up to 96% of variability using model predictors such as age and mass (as expected), a subset of observable parameters had to be selected for qualitative comparison with observations. Same log(ℒ/ℒ_⊙), log T_eff, and log g were found to all be significant and explain up to 79% of variability in the model subsets. The results of multivariate linear regression of the model data are presented in Tables B.1 and B.2 as the model data provided more possible subsets. The subsets include both metallicities with full coverage of the parameter space, convective velocity in both iron (v_conv, Fe) and helium (v_conv, He) convective layers, as well as maximum ($v_{\mathrm{conv},\mathrm{x}}^{\max }$) or integrated velocity ($v_{\mathrm{conv},\mathrm{x}}^{\mathrm{integr}}$) chosen as dependent variable.

The higher adjusted R² values in the model data compared to observations (for example of $v_{\mathrm{conv},\mathrm{Fe}}^{\max }$ subsets with log(ℒ/ℒ_⊙) > 3, R² is 0.794 against 0.602 for the Solar metallicity, or 0.722 against 0.420 for the LMC metallicity) suggest additional complexities (rotation) or observational uncertainties (in both the position in HRD and measured Θ) in the observed datasets. The overall consistency in coefficient trends between the model and observed datasets for each predictor suggests that the fundamental physical relationships between the predictors and the $v_{\mathrm{conv},\mathrm{Fe}}^{\max }$ and Θ are similar. The similar coefficient trends across different luminosities and metallicities also indicate that the predictors log(ℒ/ℒ ⊙ ), T_eff, and log g play comparable roles in influencing both macroturbulent and convective velocities, which was seen in Figs. 5 and 6. The largest discrepancies are found in low-log(ℒ/ℒ_⊙), where model dataset have near-zero $v_{\mathrm{conv},\mathrm{Fe}}^{\max }$ or $v_{\mathrm{conv},\mathrm{He}}^{\max }$ with no scatter, while observations have a large scatter around still noticeable velocities with weak to no trends with other observables.

The comparison of all tests shows that the observed scatter in macroturbulent velocities Θ for the low-mass regime (log(ℒ/ℒ_⊙) < 3 or M < 12 M_⊙) cannot be explained by same predictors that partially describe model convective velocities. For high-mass regime (log(ℒ/ℒ_⊙) > 3 or M > 12 M_⊙) and Solar metallicity, 60.2% of the observed Θ variability is explained by the same predictors and similar coefficients as Fe-bump convective velocities, pointing to the common origin of the two phenomena. Even though it is noticeably less than 79.4% described in the model convective velocities, the observed sample has a scatter due to the measurement errors that smear out the trends. Moreover, the coefficients of all predictors show consistently similar trends with luminosity and metallicity, which illustrates that the predictors play comparable roles in influencing both macroturbulent and convective velocities in the observations and models, respectively.

4.2. Principal component analysis

Principal component analysis (PCA) is a statistical technique widely utilised in the context of multivariate data analysis to investigate the underlying structure and trends within high-dimensional datasets. PCA operates by transforming the original n-dimensional data into a new coordinate system where the axes (principal components) are ordered by the amount of variance they capture from the data. The first principal component captures the maximum variance, the second captures the second most variance, and so on. This transformation is achieved through an eigenvalue decomposition of the covariance matrix of the data, yielding eigenvalues and eigenvectors that define the principal components.

In the context of comparing two n-dimensional point clouds, PCA can be applied separately to each dataset to derive their respective principal components. By examining the principal components, we can compare the direction and magnitude of the variance captured in each dataset. Specifically, if the principal components of both clouds align similarly (i.e. they have similar eigenvectors and eigenvalues), it suggests that the two datasets share similar trends. Conversely, significant discrepancies in the principal components would indicate divergent trends between the two point clouds.

In order to compare the trends in Θ and $v_{\mathrm{conv},\mathrm{Fe}}^{\max }$, we independently apply PCA to four clouds of points: observed LMC, model LMC, observed Galactic, and model Galactic datasets. In each case, we used same variables, with the difference that Θ and $v_{\mathrm{conv},\mathrm{Fe}}^{\max }$ were only used in the observed and model data, respectively, so the full list of dimensions per each point in the cloud consisted of T_eff, log g, log(ℒ/ℒ_⊙), and either Θ or $v_{\mathrm{conv},\mathrm{Fe}}^{\max }$. We included the first two components only since they provided cumulative explained variances of 89 to 95 % in each of the four datasets, as with the ratio per component shown in Table 2.

Each principal component may be represented as a linear combination of the initial set of parameters (T_eff, log g, log(ℒ/ℒ_⊙), and either Θ or $v_{\mathrm{conv},\mathrm{Fe}}^{\max }$). The coefficients for each parameter in each principal component for each dataset are provided in Table 3. Comparison of these coefficients for four datasets allowed us to draw the following conclusions. The primary component (PC1) shows consistent contributions from all parameters across all datasets, indicating that the primary trend is similar in both observations and models. The coefficients for $v_{\mathrm{conv},\mathrm{Fe}}^{\max }$ and Θ are comparable, suggesting these parameters are interchangeable in capturing the primary trend. However, the secondary component (PC2) exhibits significant differences in both the sign and magnitude of the coefficients for most parameters. These discrepancies may be attributed to specifics of the observed sample, such as uncertainties; the presence of the extra parameter v sin i, which is absent in the models; and the limited sample size and parameter span.

The visual representation of all four datasets in PC1-PC2 space is shown in Fig. 7, where circles represent model data points, and stars represent observational data points. The colour indicates the value of Θ for observations and $v_{\mathrm{conv},\mathrm{Fe}}^{\max }$ for models. Although the four datasets were analysed independently, they align very similarly in the PC1-PC2 space. The alignment of points and matching colour gradients strongly indicate that the primary trend is consistent across all four datasets. This suggests that the underlying structure and trends in the full 4-dimensional clouds are captured similarly by the first principal component across all datasets, supporting the hypothesis of Θ and $v_{\mathrm{conv},\mathrm{Fe}}^{\max }$ playing the same role in forming the trends. On the contrary, if we compare how the observed sample is misaligned from model tracks (compared to how the plot looks in the HRD instead of PC1-PC2), we can see obvious stretch and shrink of spans of the observed samples. Looking at the Galactic sample (second panel in Fig. 7), one can see the shift and shrink of the observed sample while matching colours. The shift indicates systematic differences in absolute values of Θ and $v_{\mathrm{conv},\mathrm{Fe}}^{\max }$ with the former being systematically higher. The shrink of the cloud to match lower-masses observations with higher-masses tracks indicates that the observed Θ in lower-mass regime is inconsistent with near-zero convective velocities of the models, while high-mass observed points are kept aligned with high-mass tracks. We interpret this as a support for the hypothesis of having two principally different regimes, with Θ in high- and low-mass regimes caused by different mechanisms. The LMC sample, in contrast to the Galactic sample, seems to be stretched compared to model tracks and misaligned in colour, which we interpret to be caused by the small sample size in a very limited span of masses, while the Galactic sample covers the full range of parameters.

Fig. 7. Projection of datasets on the principal components space. Left panel: Models with Z = 0.006 (circles) and observed LMC sample (stars). Right panel: Models with Z = 0.014 (circles) and observed Galactic sample (stars). The colours show the convective velocity for the model data and the macroturbulent velocities for the observed data. All four datasets were analysed independently.

5. On wave tunnelling through subsurface convective layers

Internal profiles of models computed with MESA allow us to look into the wave propagation by comparing the Brunt-Väisälä N and Lamb S_l frequencies across the stellar radius. The first two panels of Fig. 8 show an example of such a propagation diagram, i.e. internal profiles of N and S_l for a 20 M_⊙ MS star model at the Solar metallicity. The blue line shows N, the dashed brown lines – S_l for l = 1,2,5, and the black dashed line – the acoustic cut-off frequency. The left panel shows the full radius, and the middle one – zoom-in into the outer 10% of the envelope. The blue-shaded area indicates regions where g-modes and IGW can propagate. This region is defined by frequencies ω less than both N and S_l and where N² is positive (in radiative layers, as negative corresponds to convective layers where g-modes cannot propagate). The brown shaded area shows regions where acoustic modes and waves can propagate, which is defined by frequencies ω greater than both the Brunt-Väisälä frequency and the Lamb frequency but less than the acoustic cut-off frequency. Finally, the white regions where N < ω < S_l or S_l < ω < N represent evanescent zones, where waves exponentially decay. In this particular model, one can see present g-mode and p-mode cavities between the core boundary at 2 R_⊙ and iron-bump convective zone FeCZ at some 6.85 R_⊙, as well as some between the upper boundary of the FeCZ and the lower boundary of the HeCZ. However, a combination of the extensive FeCZ between 6.85 and 7 R_⊙ with rapidly decreasing S_l and next HeCZ right below the surface does not leave any potential ω from g-modes or p-modes cavities that could reach the surface. The same propagation plots for models of different masses and evolutionary states are provided in Appendix D at Zenodo (see the Data availability section for links) for both metallicity regimes – Figs. D.1 and D.2 for the Solar metallicity and Figs. D.4 and D.5 for the LMC metallicity.

Fig. 8. Propagation diagram for a 20 M⊙ MS star model. The x-axis represents the radius in solar radii, and the y-axis shows the cyclic frequency in microhertz. The left and right panels show the full radius, while the middle panel shows the external 10% of the radius. The blue line represents the Brunt-Väisälä frequency N, the brown dashed lines correspond to the Lamb frequencies Sl for different spherical harmonic degrees l, and the black dashed line shows the profile of the acoustic cut-off frequency. The blue-shaded region shows the g-mode cavity, and the brown-shaded region – p-mode cavity. The right panel shows profiles of tunnelling ωmin for l = 1,2,5 in dashed red lines and one horizontal value of ωmin taken near the convective layer (see text for details). Profiles of ωmin collapse to an apparent vertical line in the He convection zone, and this particular zone is better resolved in zoomed version of the figure in Appendix D.

As a useful addition to the propagation diagrams, we can roughly estimate the minimal wavelength a wave must have to tunnel through the subsurface convective layers before being decayed. The minimal requirement of the wavelength λ being larger than the span of the region we want the waves to tunnel through may be expressed through wavenumber k as

$$\begin{array}{c}\hfill \lambda =\frac{2\pi }{|k|}.\end{array}$$(2)

The frequency ω is connected to the wavenumber through the dissipation relation of IGW (Aerts et al. 2010):

$$\begin{array}{c}\hfill {\omega }^2=\frac{k_h^2{N}^2}{k_r^2+k_h^2}.\end{array}$$(3)

Here, k_h and k_r are the radial and horizontal components of wavenumber |k|² = k_r² + k_h².

The horizontal wavenumber k_h can be expressed in terms of the spherical harmonic degree l and radius r as

$$\begin{array}{c}\hfill k_h=\sqrt{\frac{l(l+1)}{r^2}}.\end{array}$$(4)

The condition λ > r₂ − r₁, where r₁ and r₂ are the positions of lower and upper boundaries of a convective layer, gives the minimal frequency that allows IGW to tunnel through the evanescent zone.

$$\begin{array}{c}\hfill {\omega }_{min}=\frac{(r_2-r_1)N(r)\sqrt{l(l+1)}}{2\pi r}.\end{array}$$(5)

The last equation shows that ‘tunnelling’ frequency follows the profile of the Brunt-Väisälä frequency scaled by relations of r₂ − r₁ and radius r for degrees l. The right panel of Fig. 8 shows internal profile ω_min for l = 1,2,5 as red dashed lines, where the span of FeCZ was taken as r₂ − r₁. A single minimum may be taken as a value of the ω_min profile for l = 1 at a point before the convective zone and is shown in the same figure as a solid red horizontal line. This gives a minimal value for a frequency below which waves from below the convective layer cannot tunnel through the evanescent zone. This may serve as an additional constraint for the g-mode cavity, extending the condition ω < N < S_l into ω_{min, l = 1, r = r1^*} < ω < N < S_l. Frequencies satisfying this condition occur in the blue-shaded regions in the right panel of Fig. 8. Same modified propagation diagrams with ‘tunnelling’ ω_min for models of different masses and evolutionary states are provided in the Appendix D (Figs. D.3 and D.6 for the Solar and LMC metallicities, respectively). We note that r₂ − r₁ was considered individually for each model based on which subsurface convective layer is the most extended for a given mass and age, i.e. the most restrictive on λ. Such a definition also allowed us to introduce a flag for each model depending on the extent of the blue-shaded region between the core boundary and the lower boundary of the thickest convective layer. In other words, if a certain model does not allow any ω satisfying ω_{min, l = 1, r = r1^*} < ω < N < S_l due to r₂ − r₁ requiring ω_min > N (or > S_l), we can flag such a model as forbidding IGWs from tunnelling through to the surface. We note that such a crude definition of ω_min makes it rather conservative and may only be used to mark models that forbid tunnelling, but that is not to say that the rest of the models allow for it.

Based on the above definition of models with forbidden tunnelling, we may define corresponding regions in the HRD and compare them with the observed sample. Figure 9 shows the full observed sample in the HRD with colour representing observed macroturbulent velocities Θ on top of evolutionary tracks. A thick black dashed line defines a boundary between models that forbid tunnelling of core-excited IGWs through subsurface convective layers (above the line) and the models not forbidding it (below the line). One can see that the vast majority of observed high-macroturbulence stars lie above the line, where tunnelling is not possible. This is a strong evidence against core-generated IGWs being responsible for high macroturbulent velocities observed in blue supergiants. This region occurs above 40 M_⊙ at ZAMS, 20 M_⊙ at TAMS, and 12 M_⊙ in the Hertzsprung gap for the Solar metallicity, and above 60 M_⊙ at ZAMS, 32 M_⊙ at TAMS, and 12 M_⊙ in the Hertzsprung gap for the LMC metallicity. in all these cases, the FeCZ is extended enough to block waves from the inside. These cases also have convective velocities up to 80 km s⁻¹, sufficiently large to describe the high observed macroturbulent velocities.

Fig. 9. Observed sample with panels as in Fig. 1. Left panel: LMC sample and metallicity. Middle and right panels: Galactic samples and solar metallicity. The colours show the observed macroturbulent velocities Θ. The thick dashed black line separates the region of forbidden tunnelling (above the line) through subsurface convective layers (see text for details).

While the region in the lower HRD below the dashed line lacks truly high-macroturbulent stars (except several outliers), the observed microturbulent velocities almost never fall below 20 km s⁻¹. In the same region, models predict (near-)zero convective velocities (see Fig. 10), rendering the FeCZ unlikely to cause the observed macroturbulent broadening, yet allowing for core-generated IGWs to tunnel through the subsurface convective layers. This justifies a regime of core-generated IGWs causing moderate macroturbulence for the stars with initial mass below some 12 M_⊙.

Fig. 10. Level plot of FeCZ maximal convective velocities (in solid red lines) and tunnelling frequency ωmin (in dashed blue lines; see text for details). The thick dashed black line separates the region of forbidden tunnelling (above the line) through subsurface convective layers.

Finally, the region between 12 M_⊙ and the dashed line allows for both regimes to co-exist. This co-aligns with the spectroscopic luminosity log(ℒ/ℒ_⊙) around 3.5, where we can as well see a sharp turn in the distribution of the observed macroturbulent velocities towards the highest values in Fig. 6.

6. Discussion and conclusions

In this study, we have investigated the physical origin of macroturbulent spectroscopic line broadening in hot massive stars by analysing a comprehensive sample of 594 stars with masses ranging from 2.5 to 80 M_⊙ across two different metallicity regimes: the LMC and the Galaxy with Z = 0.006 and Z = 0.014, respectively. Our sample included 86 stars with macroturbulent broadening measured for the first time. Our study was supplied with a grid of MESA stellar structure models in order to examine the properties of subsurface convective zones. These models also allowed us to explore wave propagation and the potential for IGWs to tunnel through subsurface convective layers.

To understand the underlying mechanisms behind the macroturbulent broadening, we performed several analyses. We used multivariate linear regression to identify significant predictors in both macroturbulent velocities and convective velocities predicted by 1D models. The analyses showed that in the low-mass regime, log(ℒ/ℒ_⊙) < 3 or M < 12 M_⊙, macroturbulent velocities Θ cannot be explained by the same predictors that describe the model convective velocities. At higher masses and luminosities, log(ℒ/ℒ_⊙) > 3 or M > 12 M_⊙, the observed variability in Θ may be largely explained by the same predictors and similar coefficients as the velocities occurring in the FeCZ zone. The coefficients of all predictors in this high-mass regime show consistently similar trends with luminosity and metallicity, which points towards a common origin for the observed macroturbulence and convective velocities in the models.

We also conducted a PCA in order to compare trends between the observed and model datasets. We found that the underlying structure and trends are captured similarly by the first principal component across independently analysed observed and model datasets. In the high-mass regime, the clouds of the observed and model points align in PC1-PC2 space. In contrast, in the low-mass regime, they are misplaced, showing that model values fail to predict the observed Θ.

We interpret the consistent results of both PCA and multivariate linear regression as strong statistical evidence that macroturbulent velocities in stars in different mass regimes are caused by different mechanisms. In stars with initial masses larger than 12 M_⊙, or spectroscopic luminosities log(ℒ/ℒ_⊙) > 3, macroturbulent velocities are closely related to the convective velocities in the subsurface convection zones formed due to the iron opacity bump (FeCZ), as originally proposed by Grassitelli et al. (2015a). At the same time, our analysis revealed a clear dependence of macroturbulent velocities on metallicity. In particular, Galactic stars (higher metallicity) exhibited systematically higher macroturbulent velocities compared to stars in the LMC (lower metallicity). This result aligns well with the model predictions: Subsurface convection zones are found to develop at both Galactic and LMC metallicities – but the magnitude of the convective velocities is generally lower in the latter case.

Since in all of our 1D models the regions of FeCZ are situated much deeper than the line formation region in the photosphere, it is likely that the surface macroturbulent velocities are formed by IGWs excited at the boundary of FeCZ. Excitation of IGWs by the FeCZ in massive stars has hardly been studied in multi-D simulations. However, Rogers et al. (in prep.; private communication) find that such a subsurface convection zone can generate low-frequency waves that reach the stellar surface, in conjunction with higher frequency IGWs generated by core convection. In 3D models computed by Schultz et al. (2022, 2023a), the subsurface convection layers expand further than in 1D models and cause a turbulent surface with velocity ranges matching observed macroturbulent velocities. However, such simulations are computationally demanding, and hence they were computed for several snapshots in time and a few values of stellar mass. More simulations would be required to draw definitive conclusions as to whether the subsurface convection could fully explain the observed macroturbulent broadening.

For stars with initial masses below 12 M_⊙, the convective velocities in both the FeCZ and (HeCZ) are not sufficient to explain the observed macroturbulent velocities, which are smaller than in higher-mass stars but still significant. This points to the presence of some other mechanism that is responsible for the observed line broadening in these intermediate-mass stars. This may be seen as the change in the distribution of Θ with luminosity in Fig. 6.

Our analysis of wave propagation through the stellar interiors indicated that core-generated IGWs are highly unlikely to tunnel through extensive subsurface convection layers of stars with initial masses above ∼30 M_⊙¹. The models from the same region in the HRD also have convective velocities up to 80 km s⁻¹, which is sufficient to describe the high observed macroturbulent velocities. We conclude that macroturbulent velocities are linked to the FeCZ properties in this high-mass regime based on statistical properties and trends in 1D models and observations while keeping in mind the results from modern 2D and 3D hydrodynamical simulations (Schultz et al. 2023b; Debnath et al. 2024).

The region in the HRD characterised by the initial masses below some 12 M_⊙ generally lacks stars with high macroturbulent velocity values. At the same time, only a few stars in that region show macroturbulent velocities below 20 km s⁻¹. In the same region, models predict (near-)zero convective velocities, as seen in Fig. 10, rendering the FeCZ unlikely to cause the observed macroturbulent broadening but, at the same time, allowing the core-generated IGWs to tunnel through the subsurface convective layers. This may allow for a regime of core-generated IGWs that cause moderate macroturbulence for the stars with an initial mass below some 12 M_⊙.

Lastly, the region between 12 M_⊙ and 30 M_⊙ allows IGWs to propagate to the surface and at the same time still has large FeCZ convective velocities, indicating a possible region of coexistence for the two mechanisms. More detailed 3D modelling is required to investigate the interplay of subsurface convective layers and core-generated IGWs in order to determine a more physical estimation of the boundary between the two discussed regimes.

The conclusion on the existence of different regimes aligns with the recent research of the stochastic low-frequency variability from space-based photometry in SMC, LMC, and the Galaxy (Bowman et al., in prep.). The authors conclude that FeCZ plays different roles based on age and is only crucial for evolved stars. We draw similar conclusions but emphasise the importance of the mass range apart from the age, as the FeCZ is extensive from the very ZAMS for the higher masses. Moreover, Bowman et al. (in prep.) present a convection stability window based on a Rayleigh number that aligns with our region where IGWs are allowed to tunnel through the FeCZ. Furthermore, the authors demonstrate a lack of dependence of the observed SLFV characteristics on the metallicity of the star. Assuming that the SLFV and macroturbulence phenomena originate from the same mechanism(s), the conclusion by Bowman et al. (in prep.) contradicts our findings of the macroturbulent velocities being dependent on the metallicity of the star. This discrepancy leaves room for the possibility that macroturbulence and SLFV are caused (at least partially) by different physical mechanisms. However, we may see systematically lower macroturbulent velocities in the spectra of LMC stars due to a more evolved sample than the Galactic sample we compare it to.

Finally, we highlight the complex interplay between different physical mechanisms contributing to macroturbulence in different mass regimes as well as limitations in the used stellar structure models. In particular, the models’ 1D nature and inability to account for rotating and time-dependent convection represent major uncertainties in the interpretation of the state-of-the-art (spectroscopic and photometric) observations. Furthermore, the observed sample shows a strong correlation between macroturbulence and projected rotational velocity, which we could not account for in models. Proper treatment of the interplay between these two phenomena would require development of models that incorporate both subsurface and core-generated IGWs with rotational dynamics. Such 3D rotating models could provide deeper insights into the physical processes driving macroturbulence and SLFV and improve our understanding of the atmospheric dynamics in massive stars.

Data availability

Table containing the relevant parameters of 210 stars observed in UVES/FEROS LP in the sample is available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/692/A245

Appendices C and D, which contain additional figures for Sects. 3 and 5, are available at Zenodo (DOI: zenodo.13984157). These figures include internal profiles and propagation diagrams from MESA models. Appendix C focuses on opacity bumps and subsurface convection through evolutionary stages, while Appendix D presents propagation diagrams and wave tunnelling through subsurface convective layers.

Acknowledgments

The authors thank the anonymous referee for encouraging remarks and acknowledge fruitful discussions with Jon Sundqvist, Nico Moens, and Dwaipayan Debnath on the comparison between convective velocities from 1D MESA models and turbulent velocities in the atmospheres and outflows of their unpublished 2D simulations. We also extend our heartfelt thanks to Prof. Dominic Bowman and Prof. Tamara Rogers for their valuable insights and stimulating discussions during our private communications that significantly enhanced quality of this work. The research leading to these results has received funding from the KU Leuven Research Council (grant C16/18/005: PARADISE), from the Research Foundation Flanders (FWO) under grant agreement G089422N, as well as from the BELgian federal Science Policy Office (BELSPO) through PRODEX grant PLATO. CA acknowledges funding from the European Research Council (ERC) under the Horizon Europe programme (Synergy Grant agreement N°101071505: 4D-STAR). While funded by the European Union, views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. The research is based on observations collected at the European Southern Observatory under ESO programs 1104.D-0230 and 0104.A-9001. This research has made use of the SIMBAD database (Wenger et al. 2000) and the VizieR catalogue access tool (Ochsenbein et al. 2000), CDS, Strasbourg Astronomical Observatory, France. We would like to express our sincere gratitude to Sergio Simón-Díaz for generously providing tables with the data from Simón-Díaz et al. (2017), which were crucial for our analysis.

References

Appendix A: Kippenhahn diagrams for LMC

Fig. A.1. Kippenhahn diagrams for 60, 20, 12, 7, and 3 M⊙ (rows) for metallicity Z = 0.006. Odd columns show the full radius (in mass units), and even columns show the zoom into outer 10% of radius (in radius units). Hatched areas denote radiative zones. The Ledoux-unstable convective zones are colour-coded based on convective velocity. The two left columns are plotted with the age in megayears, while the two right columns are plotted with the age in model number, which allows for a better view of the short-lived post-MS stages.

Appendix B: Full tables with results of multivariate linear regression

All Tables

All Figures

	Fig. 1. Full observed sample of high-mass stars with measured macroturbulent velocities. Left panel: LMC supergiants from Serebriakova et al. (2023) over evolutionary tracks for Z = 0.006. Middle and right panels: Galactic objects from Gebruers et al. (2022) and Simón-Díaz et al. (2017), respectively, over evolutionary tracks for Z = 0.014.
In the text

	Fig. 2. Histograms and KDE (solid line) of the measured Θ. From left to right: Full sample of 594 stars, Galactic (dark grey) versus LMC (light grey) stars, and stars for which macroturbulent velocities were inferred with the methods of Serebriakova et al. (2023, light grey) and Simón-Díaz et al. (2017, dark grey).
In the text

Fig. 3. Macroturbulent velocity as a function of stellar parameters. From top left to bottom right: Θ versus v sin i, Θ versus Teff, Θ versus log g, and Θ versus log(ℒ/ℒ⊙). The marker size is proportional to ℒ, and the marker colour shows temperature Teff. The error bars in the lower right panel are only shown for our measurements as they are not available for the sample of Simón-Díaz et al. (2017).

In the text

Fig. 4. Kippenhahn diagrams for 60, 20, 12, 7, and 3 M⊙ (rows) stellar models at metallicity Z = 0.014. Odd columns show the full radius (in mass units), and even columns show the zoom into the outer 10% of the radius (in radius units). Hatched areas denote radiative zones. The Ledoux-unstable convective zones are mapped by colour based on convective velocity. The two left columns are plotted with the age in megayears, while the two right columns are plotted with the age in the model number, which allows for a better view of the short-lived post-MS stages.

In the text

Fig. 5. Convective velocities in the HRD interpolated from MESA models for continuous representation. Top to bottom: Metallicity Z = 0.006 and 0.014. Left to right: Convective velocities in the regions of iron and helium opacity bumps. Evolutionary tracks that were used for the interpolation are shown in each panel. Coloured circles show observed macroturbulent velocities measured in this work, triangles those measured in Simón-Díaz et al. (2017).

In the text

Fig. 6. Dependence of the macroturbulent and convective velocities on spectroscopic luminosity of the star. From left to right: Observed macroturbulent velocities (left) and convective velocities for the iron (middle) and helium (right) opacity bumps predicted with MESA. The convective velocities in the middle and right panels have been interpolated on the observed positions in the HRD from Fig. 5. Symbol size and shape are the same as in Fig. 3.

In the text

	Fig. 7. Projection of datasets on the principal components space. Left panel: Models with Z = 0.006 (circles) and observed LMC sample (stars). Right panel: Models with Z = 0.014 (circles) and observed Galactic sample (stars). The colours show the convective velocity for the model data and the macroturbulent velocities for the observed data. All four datasets were analysed independently.
In the text

Fig. 8. Propagation diagram for a 20 M⊙ MS star model. The x-axis represents the radius in solar radii, and the y-axis shows the cyclic frequency in microhertz. The left and right panels show the full radius, while the middle panel shows the external 10% of the radius. The blue line represents the Brunt-Väisälä frequency N, the brown dashed lines correspond to the Lamb frequencies Sl for different spherical harmonic degrees l, and the black dashed line shows the profile of the acoustic cut-off frequency. The blue-shaded region shows the g-mode cavity, and the brown-shaded region – p-mode cavity. The right panel shows profiles of tunnelling ωmin for l = 1,2,5 in dashed red lines and one horizontal value of ωmin taken near the convective layer (see text for details). Profiles of ωmin collapse to an apparent vertical line in the He convection zone, and this particular zone is better resolved in zoomed version of the figure in Appendix D.

In the text

	Fig. 9. Observed sample with panels as in Fig. 1. Left panel: LMC sample and metallicity. Middle and right panels: Galactic samples and solar metallicity. The colours show the observed macroturbulent velocities Θ. The thick dashed black line separates the region of forbidden tunnelling (above the line) through subsurface convective layers (see text for details).
In the text

	Fig. 10. Level plot of FeCZ maximal convective velocities (in solid red lines) and tunnelling frequency ωmin (in dashed blue lines; see text for details). The thick dashed black line separates the region of forbidden tunnelling (above the line) through subsurface convective layers.
In the text

Fig. A.1. Kippenhahn diagrams for 60, 20, 12, 7, and 3 M⊙ (rows) for metallicity Z = 0.006. Odd columns show the full radius (in mass units), and even columns show the zoom into outer 10% of radius (in radius units). Hatched areas denote radiative zones. The Ledoux-unstable convective zones are colour-coded based on convective velocity. The two left columns are plotted with the age in megayears, while the two right columns are plotted with the age in model number, which allows for a better view of the short-lived post-MS stages.

In the text