© ESO, 2016

1. Introduction

In 1992, the new stellar opacities from Rogers & Iglesias (computed with the OPAL code) revolutionized asteroseismology of massive stars by providing an explanation to what was called the β Cephei problem: none of the commonly considered excitation mechanisms could account for the presence of pulsations in these B-type stars. As a consequence of these new computations, the mean Rosseland opacity was increased by a factor of three at log T ≃ 5.2, allowing the excitation of modes in these stars through the κ-mechanism in the iron opacity bump (e.g., Kiriakidis et al. 1992). The inclusion of some additional metals in the OPAL computations by Iglesias & Rogers (1996) increased even further the opacity in that region. Meanwhile, similar results – though with a slightly hotter metal-opacity bump – were also obtained in parallel by the OP project (Badnell et al. 2005).

Pamyatnykh (1999) was the first to present an updated overview of the theoretical predictions for instability domains in the upper main sequence (MS) accounting for these new opacities. These domains, which mainly refer to the regions of the Hertzsprung-Russell (HR) diagram where the so-called β Cephei and SPB (slowly pulsating B) stars are located, were subsequently revised by Miglio et al. (2007a,b) considering the updated solar metal mixture¹ of Asplund et al. (2005), and different choices in the assumed opacities (see also Pamyatnykh 2007).

All these computations improved the agreement between models and observations. However, some problems still remained. For example, standard theoretical models were not able to fit all the observed frequencies in some Galactic β Cephei stars (e.g., 12 Lac and ν Eri, Ausseloos et al. 2004; Pamyatnykh et al. 2004; Dziembowski & Pamyatnykh 2008; Desmet et al. 2009). In addition, some Galactic late O-type stars were found to exhibit β Cephei-type pulsations, despite being located outside the predicted instability domains (e.g., Briquet et al. 2011). Moreover, with the progress in detection capabilities, B pulsating stars were also discovered in the Small Magellanic Cloud (e.g., Kourniotis et al. 2014, and references therein) while, adopting the typical metallicity of this galaxy (Z ≈ 0.003, e.g., Buchler 2008), no modes could be excited in a theoretical standard β Cephei model (Salmon et al. 2009). The solution to these discrepancies between theoretical predictions and observations may still be related to the considered opacities (see, e.g., Salmon et al. 2012; Turck-Chièze & Gilles 2013; Moravveji 2016), or to the existence of other types of driving mechanisms of stellar oscillations not accounted for in the computations.

Moving up in the HR diagram to the O star and B supergiant domain (hereafter OB stars), we enter a more uncharted territory in terms of instability predictions, in particular for the more evolved models. Due to the very high contrast in density between the core and the superficial layers in post-MS models (compared to MS models for the same initial mass), the number of nodes of the pulsation modes in the central layers become extremely high and numerical problems appear. However, it has still become clear from both an observational and a theoretical point of view that OB stars present various types of oscillations. Indeed, as an extension to the β Cephei and SPB instability strips in the B star domain, these more massive stars are predicted to present pressure (p-) and gravity (g-) pulsation modes excited by the κ-mechanism. Depending on the stellar effective temperature (T_eff), the iron or the helium opacity bumps are mainly responsible for the occurrence of these modes in stars with masses ≳9 M_⊙ (see, e.g., Dziembowski et al. 1993; Pamyatnykh 1999; Saio 2011). In addition, Saio (2011) showed that oscillatory convective modes are expected to be excited in the HR diagram in a large portion of the MS and post-MS region of the massive-star domain. These are modes punctually excited thanks to the convective motions in the iron and helium opacity bumps. Furthermore, modes excited by the ϵ-mechanism (e.g., Noels & Scuflaire 1986; Scuflaire & Noels 1986; Unno et al. 1989), strange mode instabilities (e.g., Saio et al. 1998; Glatzel 2009, and references therein), and even stochastically excited waves excited by various driving mechanisms (Belkacem et al. 2010; Samadi et al. 2010; Shiode et al. 2013; Mathis et al. 2014; Aerts & Rogers 2015; Grassitelli et al. 2015) are predicted in massive OB stars. Actually, some of these predictions have already found empirical support (e.g., Aerts et al. 2010a; Degroote et al. 2010; Moravveji et al. 2012; Buysschaert et al. 2015). Recent reviews describing the current state of theoretical knowledge about pulsations in massive OB stars can be found in Saio (2011), Godart et al. (2014), Saio et al. (2015) and Samadi et al. (2015).

The slow but sure increase in the amount of observational data coming from space missions (e.g., MOST, CoRoT, Kepler, K2, BRITE), supported by on-ground spectroscopic material, is not only showing us the rich and complex variety of pulsational-type phenomena present in the OB-star domain (see Aerts 2015, and references therein), but these data are also providing us with invaluable empirical information to step forward in the characterization and understanding of the pulsational properties of these massive stellar objects. This situation will improve thanks to future missions such as Transiting Exoplanet Survey Satellite (TESS) and Planetary Transits and Oscillations of stars (PLATO), or coordinated on-ground telescope facility networks (e.g., Stellar Observations Network Group, SONG). It is timely to combine all the promptly-available empirical information with our current theoretical knowledge of stellar oscillations in massive stars. In this regard, it is interesting to note that, despite the remarkable progress made in the last decades, we still lack a comprehensive, homogeneous, in-depth pulsational stability analysis in the massive star domain (from the zero age MS to the more evolved phases, and the complete massive star range above 3 M_⊙). Indeed, results are limited to low degree modes (ℓ ≤ 3), even combining heterogeneous computations by several authors (Miglio et al. 2007a; Pamyatnykh 2007; Godart 2011; Saio 2011; Daszyńska-Daszkiewicz et al. 2013; Saio et al. 2015; Moravveji 2016), except for the study of Balona & Dziembowski (1999) that limited their high-degree computations to MS stars.

Higher-degree modes are usually disregarded because they are difficult to detect photometrically (e.g., Balona & Dziembowski 1999; Aerts et al. 2010b), but they can be important when interpreting spectroscopic observations (see, e.g., Schrijvers et al. 1997). An example of this situation is highlighted in Simón-Díaz et al. (2016), where we present the first comprehensive attempt to evaluate empirically the proposed hypothesis that the so-called macroturbulent broadening in O and B stars is produced by the collective effect of multiple non-radial pulsation modes (e.g., Lucy 1976; Aerts et al. 2009). Using high-spectral-resolution, single-snapshot observations of ≈430 Galactic stars² with spectral types ranging from O4 to B9, we investigated the potential connection between stellar oscillations driven by a heat mechanism and this ubiquitous and dominant non-rotational broadening component shaping the line profiles of O stars and B Supergiants. The result was not very promising when accounting only for dipole modes predictions (ℓ = 1). However, the consideration of a larger number of modes of a higher degree might help to improve the situation. As a continuation to the work by Godart (2011) and motivated by these recent advances in the investigation of macroturbulent broadening in the OB star domain, we present in this paper new theoretical predictions for instability domains of ℓ ≤ 20 non-radial modes for stars with masses between 3 and 70 M_⊙. We aim to provide a global homogeneous overview of the instability domains predicted in the upper part of the HR diagram, together with some indications about the characteristics (in terms of type of modes and frequency spectrum) of the various unstable modes found in the whole region. We also expect that this theoretical work will become a useful tool for the interpretation of data gathered in on-going and future observational campaigns, especially in spectroscopy.

The paper is structured as follows. The model-input physics is presented in Sect. 2, while we introduce in Sect. 3 the diagnostic diagrams used for the determination of the instability boundaries for pulsations. The instability domains resulting from a comprehensive pulsational instability analysis for MS and post-MS models from 3 to 70 M_⊙ for ℓ = 1 and 2 ≤ ℓ ≤ 20 are presented in Sects. 4.1 and 4.2, respectively. This is complemented with a discussion about the behavior of high-degree modes as the stars evolve from the MS to the post-MS in Sect. 4.3. We then investigate, in Sect. 5, possible correlations between the location of stars characterized by having an important contribution of so-called macroturbulent broadening to their line profiles, and the predicted instability domains associated with heat-driven pulsating modes (Sect. 5.1). We then discuss the physical ingredients of the models that can alter the results, such as a change in the metallicity or in the adopted opacities (Sect. 5.2). With this study, we want to assess whether this spectroscopic feature is an observational signature of the existence of a distribution of motions in the line-formation region induced by this specific type of stellar pulsations. Lastly, other possible sources of macroturbulent broadening are briefly discussed in Sect. 5.3, and the main conclusions of our work – along with some future prospects – are summarized in Sect. 6.

2. Stellar models and non-radial modes

We concentrate in this paper on the instability domain predictions for heat-driven pulsation modes produced by the κ-mechanism in the iron opacity bump. To this aim, we benefit from the computations performed in Godart (2011). The stellar models were computed using the ATON evolutionary code (Ventura et al. 2008). In order to cover the whole mass range of observations, we extended the grid of stellar models to 3–70 M_⊙. We briefly mention here the main physical ingredients of these models:

• convective transport is treated following the Mixing LengthTheory of convection (Böhm-Vitense 1958) andby adopting the Schwarzschild criterion for convection(Schwarzschild 1906),

• no overshooting is included: asteroseismology of massive stars shows a rather large dispersion in the determination of the overshooting parameter (e.g., Noels et al.2015, though Castro et al.2014, have argued that a mass-dependent overshooting is necessary to recover the HR diagram features at high masses),

• mass loss is taken into account for masses larger than 7 M_⊙ following the prescription of Vink et al. (2001),

• we used the OPAL opacity tables (Iglesias & Rogers 1996), extended with the opacities of Ferguson et al. (2005) for T< 6000 K,

• the considered metallicity is Z = 0.015, and we used the metal mixture determined by Grevesse & Noels (1993),

• microscopic diffusion and radiative levitation are ignored.

At this point, we remark that our main set of computations (see also below) do not use the most recently accepted values of metallicity (Z = 0.014 and the metal mixture from Asplund et al. 2009). In addition, while we considered the OPAL opacity tables, we are aware of the availability of other computations, such as the OP tables, which favor the excitation of modes in the considered domain of stellar parameters (e.g., Miglio et al. 2007a; Pamyatnykh 2007). The computations performed here are still valid for the main purpose of this work, namely the investigation of the possible connection between heat-driven modes and the so-called macroturbulent broadening in O and B stars. A detailed explanation of the reason for using the set of assumptions quoted above is presented in Sect. 5.2, along with notes on the effect on the computed instability domains (and therefore on the conclusion of our study) of considering a higher metallicity, a different metal mixture, and/or the opacities from the OP project.

Adiabatic and non-adiabatic computations for non-radial pulsations in the whole range of masses described above were performed using the Liège Oscillation code (LOSC, Scuflaire et al. 2008a) and the non-adiabatic code MAD (Dupret et al. 2003) combined with the numerical method described in Godart et al. (2009) for the computations of g-modes in models presenting an intermediate convective shell surrounding the He core. This method allows us to avoid the computation of the full, dense spectrum of non-radial modes by pre-selecting the modes that are reflected at the bottom of the convective shell and which are potentially excited. The adiabatic and non-adiabatic frequencies were computed in the range ω = 0.05 to 20 where ω is the dimensionless frequency³. This range covers periods from minutes to several hundreds of days. The degree ℓ of the mode was chosen to vary from one to 20 considering zonal modes only (m = 0). We highlight the novelty and importance of the ℓ> 3 (and up to 20) computations for the purpose of the study presented here. As already mentioned, higher-degree modes are usually disregarded because of the difficulty in detecting them using photometric observations (the amplitudes are averaged out over the unresolved stellar disc; see Balona & Dziembowski 1999, for an analysis of the visibility of high-degree modes). However, these modes, especially when they are combined together in a dense frequency spectrum, have a more important effect on the line profiles. Therefore, they become interesting in a spectroscopic context and, in particular, for the investigation of line-profile variability and line-broadening effects (see, e.g., Schrijvers et al. 1997; Aerts et al. 2009; Simón-Díaz et al. 2016).

To identify the excited modes, we first evaluate whether the work produced by the star during a pulsation cycle is positive (e.g., Unno et al. 1989). This criterion, used in all non-adiabatic codes, gives a first indication of the stability of the pulsation. However, it may not be sufficiently reliable when studying instabilities of massive stars. Under certain circumstances, some modes, which would be considered excited according to the work integral criterion, present peculiar eigenfunctions. In particular, in the most massive and evolved models, we have found excited modes that are misinterpreted by the non-adiabatic code due to the huge number of nodes compared to the limited number of mesh points we imposed in the code. Indeed, the asymptotic formula for g-modes (Tassoul 1980), ${\mathit{k}}^{\mathrm{2}}\mathrm{\approx }\frac{{\mathit{N}}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}\frac{\mathit{\ell }\mathrm{(}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}}{{\mathit{r}}^{\mathrm{2}}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{(}\mathrm{if}\hspace{0.17em}\hspace{0.17em}\mathit{\sigma }\mathrm{\ll }\mathit{N}\mathit{,}{\mathit{L}}_{\mathit{\ell }}\hspace{0.17em}\mathrm{)}\mathit{,}$(1)shows that for a given frequency (σ), the number of nodes (k) for a mode is very large in the region of the peaks of N² (i.e., the Brunt-Väisälä frequency, see Eq. (3)) and/or for large values of ℓ. When studying the eigenfunctions in detail, these modes are actually not physically excited.

These issues are recurrent over the whole spectra of computed frequencies for the most evolved and/or the most massive models, but the modes with the largest periods are the most problematic to compute due to their very large number of nodes (for ω ≲ 0.2 though, this upper limit increases with increasing ℓ). Therefore, we decided to inspect visually every modal diagram (see Sect. 3) resulting from the non-adiabatic code computations in order to identify the potential unphysical modes. For all these dubious modes, we evaluate the reliability of the results by studying in detail the corresponding eigenfunctions. In particular, whenever necessary, we increased the number of mesh points (usually around 5000 points) in the most critical regions of the star such as in the non-zero mean molecular weight gradient zone or where the work produced by the star reaches high absolute values. As a result, we ended up with what we call cleaned modal diagrams in which only these physically meaningful modes are included.

Fig. 1 Left: modal diagram: range of dimensionless frequencies for a 50 M⊙ MS model (labeled B in Table 1) as a function of the degree ℓ of the mode. Stable modes are represented by gray points while excited modes are shown with solid blue circles. Right: propagation diagram for the same model. The solid red line stands for the Brunt-Väisälä frequency, and the Lamb frequency for three different degrees (ℓ = 1, 3 and 20) is represented by the dotted, dashed, and dot-dashed lines, respectively. The propagation zones for p and g modes with ℓ = 1 are highlighted by the dark gray and light blue shadow regions, respectively. The solid horizontal black line stands for an example of the propagation zones for a mode of ω = 1.0 with ℓ = 3.

3. Diagnostic diagrams

In order to obtain a global picture of the stable and unstable modes for each computed model, we used a version of the so-called modal diagrams⁴ in which the dimensionless frequency ω (or, alternatively, the frequency f) is displayed as a function of the degree ℓ of the mode.

The left panel of Fig. 1 shows an illustrative modal diagram of this type for the particular case of a 50 M_⊙ MS model (see model B in Table 1 for the global properties). Stable and unstable modes are represented by gray and big blue points, respectively. As discussed in the previous section, all the modal diagrams presented along this paper have been cleaned of the unphysical modes encountered in the computations (i.e., only the excited modes that have survived the visual inspection of the eigenfunctions are marked in blue).

The information provided by the modal diagrams is also complemented by analyzing the associated propagation diagrams. This second type of diagnostic diagram allows us to identify the cavities in which p and g modes can propagate in a given model and for a certain value of ℓ (e.g., Unno et al. 1989). This is achieved by plotting together as a function of, for example, log T, the Lamb frequency, ${\mathit{L}}_{\mathit{\ell }}^{\mathrm{2}}\mathrm{=}\mathit{\ell }\mathrm{(}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}\hspace{0.17em}{\mathit{c}}^{\mathrm{2}}\mathit{/}{\mathit{r}}^{\mathrm{2}}\mathit{,}$(2)where c is the sound speed and r is the radius, and the Brunt-Väisälä frequency, ${\mathit{N}}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{g}}_{\mathit{r}}}{\mathit{r}}\frac{\mathrm{d}\ln \mathit{P}}{\mathrm{d}\ln \mathit{r}}\hspace{0.17em}[\hspace{0.17em}(\frac{\mathit{\partial }\ln \mathit{\rho }}{\mathit{\partial }\ln \mathit{T}}{)}_{\mathit{P,\mu }}\mathrm{(}{\mathrm{\nabla }}_{\mathrm{ad}}\mathrm{-}\mathrm{\nabla }\mathrm{)}\mathrm{-}{\mathrm{\nabla }}_{\mathit{\mu }}\hspace{0.17em}(\frac{\mathit{\partial }\ln \mathit{\rho }}{\mathit{\partial }\ln \mathit{\mu }}{)}_{\mathit{T,P}}]\mathit{,}$(3)where g_r is the local gravity, P, ρ, T, and μ are the pressure, density, temperature, and mean molecular weight, respectively, ∇_ad and ∇ are the adiabatic and real temperature gradient, and ∇_μ is the mean molecular weight gradient ∇_μ = dlnμ/dlnP. This is illustrated in the right panel of Fig. 1, which shows the propagation diagram corresponding to the model B mentioned above. The Lamb frequency (for ℓ = 1,3 and 20) and the Brunt-Väisälä frequency are drawn with dashed and solid lines, respectively. The propagation zones for p and g modes for the case ℓ = 1 are indicated by the dark gray and light blue areas, respectively, while the solid horizontal black line stands for an example of the propagation zones for a mode of ω = 1 with ℓ = 3. Generally speaking, those modes with ω²> max$\mathrm{\left(}{\mathit{N}}^{\mathrm{2}}\mathit{,}{\mathit{L}}_{\mathit{\ell }}^{\mathrm{2}}\mathrm{\right)}$ are p modes and propagate in the more superficial layers, while those modes with ω²< min$\mathrm{\left(}{\mathit{N}}^{\mathrm{2}}\mathit{,}{\mathit{L}}_{\mathit{\ell }}^{\mathrm{2}}\mathrm{\right)}$ correspond to g modes and propagate in the deeper layers. However, as illustrated by the horizontal line, there are also mixed modes, which present both behaviors depending on the regions of the stellar interior where they are propagating. For example, the ω = 1.0, ℓ = 3 mode in Fig. 1 will have characteristics of p modes in the outer regions and of g modes in the inner one.

In summary, as described above, the combined information provided by the modal and the propagation diagrams shown in Fig. 1 allows us to identify the frequency spectrum, the associated propagation cavity, and the mode type of excited modes in each considered stellar model.

3.1. Some notes on the propagation diagrams

For a better understanding of the following sections, we describe here in more detail the propagation diagram for massive stars (Fig. 1). We concentrate on the behavior of the Brunt-Väisälä frequency since it plays a major role in the excitation of the instabilities computed in this paper. As shown in the propagation diagram, the Lamb frequency (L_ℓ) is very large in the center and decreases monotonously towards the surface, following the decreasing c and the increasing r (see also Eq. (2)). Its behavior is very similar from one model to another. On the contrary, the Brunt-Väisälä frequency (N) can be very different depending on the specific mass and evolutionary state of the associated stellar model. During the MS, N² = 0 in the convective core (see Eq. (3)). In the layers surrounding the core, the presence of a ∇_μ (due to the receding convective core during the evolution) produces a peak in the Brunt-Väisälä frequency whose intensity depends on the sharpness of the μ gradient (e.g., Noels et al. 2010). As a matter of fact, MS stars with intermediate masses (3–20 M_⊙), usually present a unique peak in N²(r) while the star evolves from the zero age main-sequence (ZAMS) due to the increasing ∇μ region. Therefore, in that case, the star presents mainly two big cavities (p and g), and the frequency spectrum is characterized by avoided crossings between low-order p and g modes during which the modes exchange their behavior.

In more massive stars, the behavior of N²(r), and hence the propagation diagram, is more complex. In addition to the sharp peak due to the ∇μ region on the MS, the Brunt-Väisälä frequency presents other bumps almost equally important, leading to new cavities in which different types of modes may propagate. The propagation diagram (Fig. 1, right panel) illustrates the most important characteristics of the Brunt-Väisälä frequency in massive stars. Starting from the innermost region of the star (A), the actual limit of the convective core is shown by a first discontinuity in N² at log T ~ 7.4 (B). The subsequent decrease in N² (B-C) is due to the decrease of g and μ (see Eq. (3) where P/ρ ~ T/μ) in the region of variable μ. It is followed by another discontinuity at log T ~ 7.1 (C) once ∇μ becomes 0. This discontinuity corresponds to the maximal extension of the convective core during MS. Since above the ∇μ region, ∇_rad becomes smaller following the decreasing opacity, the shift between both temperature gradients in Eq. (3) gets larger and N² increases (C-D). This feature is discussed in more detail below. The next minimum in N² is produced by the deep opacity bump at log T ~ 6.25 due to the L-transition of iron (E). One can also see the convective zone due to the iron opacity bump at log T ~ 5.2 (F-G) and a superficial convective zone due to the ionization of helium (log T ~ 4.6, H-I).

Apart from the large deep opacity bump, the most important characteristic of the N² frequency in massive stars is related to the sharp feature produced by the important radiation pressure around log T ~ 7.0 (C-D). Actually, the increasing contribution of the radiative pressure to the total pressure (with mass, and during the evolution on the MS, due to the increasing central temperature) induces a decreasing adiabatic temperature gradient. Since both gradients (∇_rad and ∇_ad) are already close to one another, the ∇_rad becomes more easily larger than ∇_ad, leading eventually to the appearance of an intermediate convective zone (ICZ) on the post-MS (see Sect. 4.1). However, even before the onset of an ICZ, the discontinuity in the Brunt-Väisälä frequency (C), already present during MS at the upper limit of the ∇_μ region, shows an increasingly small minimum value as the mass increases due to the increasing radiation pressure leading to a smaller coefficient (∂lnρ/∂lnT)_P,μ (Eq. (3)). This lower minimum value in the Brunt-Väisälä frequency will have a similar effect as the ICZ, as we will see in Sect. 4.2.

4. Instability strips

Fig. 2 Instability domain predictions in the sHR diagram for dipole modes (ℓ = 1, see Sect. 4.1). P-mode instabilities are represented by the blue area and the g-mode domain is shown in red. The ATON non-rotating evolutionary tracks are over plotted. The horizontal dashed line indicates the limit above which strange mode oscillations are most likely expected. Black asterisks stand for nine selected models la A to I which will be more deeply investigated in Sect. 4.2. The global properties of these models are given in Table 1.

By investigating modal diagrams, we have determined the instability boundaries of pulsations due to the κ-mechanism in the iron opacity bump. These results, along with the ATON, non-rotating evolutionary tracks are presented in the spectroscopic Hertzsprung-Russell diagram (sHR diagram, Langer & Kudritzki 2014; Castro et al. 2014), separating the case of dipole modes (ℓ = 1, Fig. 2) and higher-degree modes (2 ≤ ℓ ≤ 20, Fig. 3). These figures are complemented with some illustrative, cleaned modal diagrams for 11 selected models, labeled from A to K (Figs. 4 and 5). The location of these models in the sHR diagram are shown by black asterisks in Figs. 2 and 3 and their main parameters are summarized in Table 1. In this case, we have used the dimensional frequency f in cycles per day (instead of the dimensionless frequency ω as in Fig. 1) for a better utility of these diagrams from an observational point of view. The values of τ_dyn (in hours) provided in Table 1 for each model, can be used to transform f into the dimensionless variable ω.

The panels of Figs. 4 and 5 are organized following the location of the selected models in the sHR diagram. In Fig. 4, every row displays models with similar masses (except for model C; decreasing masses from the top to the bottom row) and the models are evolving from the beginning of the MS to the post-MS from left to right (i.e., every column corresponds to a given evolutionary state). Rough evolutionary states are indicated for each model in the right bottom corner of the modal diagram (ZAMS, EMS, MMS, TAMS, and post-MS stand for zero age MS, early MS, middle of MS, terminal age MS, and hydrogen-burning shell models). Finally, Fig. 5 represents models with different masses (8 to 4 M_⊙ from the top to the bottom panel) selected at the same evolutionary stage (MMS, X_c = 0.3).

Modes appear into two groups depending on their frequencies (in particular, see models B, C, and E in Fig. 4):

• modes with high frequencies (small periods) behave mainly likep modes; these modes are excited with low orders;

• modes with low frequencies (large periods) have a general g-mode behavior; these modes are usually excited with high orders.

The limit between these two kinds of modes depends on the frequency but also on the degree ℓ of the mode. For this reason, we do not use the usually accepted criterion in dimensionless frequencies of ω > 1 (resp. ω < 1) for p modes (resp. g modes), and we mainly refer in the following analysis to the high-frequency and low-frequency modes, especially when discussing the results for ℓ > 1.

For a better understanding of the instability domain predictions, we first concentrate on the results for ℓ = 1, and then extend the description of the results to higher-degree modes (2 ≤ ℓ ≤ 20), for which the physics discussed in the ℓ = 1 case applies too.

Fig. 3 Same symbol key as in Fig. 2, except that the instability domain predictions are now done for 2 ≤ ℓ ≤ 20 pulsation modes (see Sect. 4.2).

Fig. 4 Range of frequencies in cycle per day (c/d) for models A to I. Gray dots stand for the frequencies of the stable modes computed, whereas excited mode frequencies are represented by blue bigger dots. From left to right, the rows represent the evolution from the MS to the post-MS: models A to C for high masses (30–50 M⊙), models D to F for 15 M⊙, and models G to I for 6 M⊙ MS models only (see Fig. 2 and Table 1). Rough evolutionary states are indicated for each model in the right bottom corner (ZAMS, EMS, MMS, TAMS, and post-MS stand for zero age MS, early MS, middle of MS, terminal age MS, and hydrogen-burning shell models). We note that the detectability of the pulsation modes in photometry is limited to the low degree (ℓ up to 3–4), whereas higher degree modes should be able to be observed thanks to spectroscopy.

4.1. Instability domains for dipole modes (ℓ = 1)

The results for dipole modes (ℓ = 1) are presented in the sHR diagram in Fig. 2, while the associated modal diagrams for the selected models A–K are highlighted with a vertical gray band in Figs. 4 and 5. We separate the discussion below into the three major instability domains.

We start by discussing the lower mass regime (from 3 to 6 M_⊙) for which instabilities are found for MS stars only. These models are characterized by high-order g modes only (represented by the red area limited by dashed lines in Fig. 2), which basically maps the SPB instability strip for ℓ = 1. The periods related to these modes are of one to three days (see panels I and K in Figs. 4 and 5, respectively).

Moving up in the sHR diagram, we find low-order p and g modes, which are predicted on the MS and also on the early post-MS of massive models. These modes are represented by a blue area limited by dotted lines in Fig. 2. In the intermediate mass range (~10–30 M_⊙), these excited modes with periods of the order of six hours constitute the β Cephei instability strip (e.g., model D, ℓ = 1). For the models with larger masses (M ≥ 30M_⊙), some adiabatic strange modes are also excited with periods of a few hours (models A and B), in addition to the low order p and g modes. These adiabatic strange modes are excited in the iron opacity bump due to the κ-mechanism and their strange mode behavior is related to the amplitude enhancement of the mode trapped into a superficial cavity (e.g., Saio et al. 1998; Glatzel 2009). The nature of these modes has been confirmed by investigating the eigenfunctions. These instabilities extend to the post-MS phase, especially for models with masses up to 40 M_⊙ (with dimensionless frequency ω ~ 3, i.e., periods of the order of hours). In particular, in these post-MS models, the instabilities are essentially due to low-order g modes that present a mixed-mode behavior, that is, they propagate into the central layers as g modes and they have a p-mode behavior in the superficial layers.

Finally, a large number of high-order g modes are expected to be excited mainly during the post-MS phase of stars with masses larger than 9 M_⊙ (red area in Fig. 2 and models C and F in Fig. 4). These modes have periods of about ~5−15 days, and this region of the sHR diagram is often referred to as the extension of the SPB instability strip to higher masses. They are predicted in models which present an ICZ surrounding the radiative core. The ICZ acts as a barrier and prevents the modes from entering the damping core. Indeed, the density contrast between the helium core and the envelope is very high in massive post-MS models. As a result, the Brunt-Väisälä frequency reaches huge values in the core and a large radiative damping occurs (the radiative damping is proportional to N²: see Eq. (3) in Dupret et al. 2009). However, the ICZ, which develops above the radiative core during post-MS evolution, reflects the modes and allows therefore a sufficient excitation by the κ-mechanism in the iron opacity bump (Saio et al. 2006; Gautschy 2009; Godart et al. 2009), in a similar way to the excited modes in SPB stars. In our models, the ICZ appears in the models of 12 M_⊙ (up to higher masses), although g-mode instabilities are already predicted in models of 9 M_⊙. Actually, we found here that the sharp feature in the Brunt-Väisälä frequency discussed in Sect. 3 (at log T ~ 7.1 in Fig. 1) already produces the reflection needed for exciting the g modes if the minimum is low enough, in agreement with Daszyńska-Daszkiewicz et al. (2013), as it is the case for our models with M ≥ 9M_⊙.

The instability strip for these high-order g modes (low frequencies) overlaps with the low-order p- and g-modes instability strip (high frequencies) for the most massive models. The overlapping of both areas appears in purple in Fig. 2. Mainly located on the post-MS (for ~15–30 M_⊙), the instability region reaches the MS for masses of 40–50 M_⊙. When both these frequency ranges are excited at the same time (as is the case at the end of the MS for a 40 M_⊙ star for example, for which the frequency spectrum is similar to model C), a stable region lies between the frequency domains of excited modes (Fig. 4). The amplitude of the modes in that range of intermediate radial orders are larger in the core than in the envelope, and the important radiative damping prevents the κ-mechanism from being sufficiently efficient. Two effects explain this behavior: (1) the evanescent zone is larger in this frequency range which leads to a larger coupling between p and g cavities and thus to a larger radiative damping (Dupret et al. 2008) and; (2) the amplitudes of the eigenfunctions in the superficial layers show a minimum in this frequency range (Dziembowski et al. 1993).

All these theoretical predictions for the dipole instability domains are in good agreement with the work of Saio (2011). For a more extensive description of the properties of the various theoretical instability domains for massive stars for low degree ℓ and for modes excited by other mechanism than the opacity mechanism we refer the reader to Saio (2011, 2015), Saio et al. (2015) and the reviews of Godart et al. (2014), Samadi et al. (2015).

Fig. 5 Same symbol key as in Fig. 4 but for models J to K, selected, from top to bottom, at the same evolutionary state (Xc ~ 0.3) but with decreasing masses from 8 to 4 M⊙ (see Fig. 2 and Table 1).

4.2. Instability domains for 2 ≤ ℓ ≤ 20

As stated in Sect. 1, the motivations that guided us through this work are the investigation of the physical origin of macroturbulent broadening and the interpretation of the observed line-profile variability of OB stars. In this context, we extended our computations to high-degree modes (ℓ> 3), which signatures can be visible in spectroscopy, though they are not detectable in photometry. The instability domains predicted for ℓ = 2 to 20 are shown in Fig. 3 by shaded areas. These instability boundaries are complemented by the information contained in the modal diagrams (as illustrated in Figs. 4 and 5).

We divided our predictions for the frequency ranges and instabilities into two groups. First, a group of modes with high frequencies (small periods), for which ω ~ 2−3, is represented by the blue area limited by dotted lines in Fig. 3. This group of pulsations presents mostly a p-mode behavior and constitutes a rather sparse and discrete spectrum (see modes with f ~ 4.0 c/d for models A and D, f ~ 1.5−2 c/d for models B and E, and modes with f ~ 0.8 c/d for model C). The periods associated to these modes correspond to approximately eight hours for models A and D, approximately 12 h for model B and E, and approximately 24 h for model C, with increasing periods along the evolution. As for dipole modes, these modes are excited with low orders. Secondly, a group of lower frequencies (large periods), for which ω ~1, is represented in the red area with dashed contour lines. These pulsations mainly present a g-mode behavior and are mainly excited with high orders. These modes constitute a much denser spectrum (see the lower frequency group in models B, C, E, F, G, H, and I). The periods range from a couple of hours to ten days. Usually well detached, these two groups of modes reach a roughly common frequency range in the most massive models (see model B), in which the excited frequencies from p and g modes are found very close to each other and the modes behave as mixed modes, in agreement with Balona & Dziembowski (1999; see the explanation in Sect. 4.1).

As expected, the instability domains cover a wider region of the sHR diagram when considering high-degree pulsations. This effect is noticed in low-order p- and g-mode instabilities, for which the instability domain mainly widens towards smaller masses and reaches MS models of 8 M_⊙ (rather than 10 M_⊙ for dipole modes). It is also marked when considering the high-order g-mode domain. Indeed, the instability bands reach hotter T_eff in the whole instability region. In particular, the post-MS gap between low-order p and g modes and high-order g modes is now completely filled in by excited modes in models from 9 to 12−15 M_⊙. In these models, as for example in model E, g modes are only predicted when considering ℓ ≥ 3 which would, in this case, make them difficult to observe in photometry. Furthermore, theoretical computations predict now excited g modes in the high-mass models (up to 70 M_⊙). See for example model B, for which excited quadrupole modes should be detectable in photometry (P ~ 2.5 days), though g modes are not predicted for ℓ = 1. Finally, intermediate mass models (from 5–12 M_⊙) now present instabilities on the MS for 5 to 10 M_⊙. These modes have quite large frequencies (small periods reaching the few hours) – with regard to their g-mode behavior – due to their mixed mode behavior. Here again, while models G and J are characterized by the excitation of a dense spectrum of g modes, these are hardly expected to be detected with photometry, while the situation is more favorable for models H, I, and K. These results show that, with the increasing quality of the photometric datasets offered by the space missions, it should become more and more possible to detect intermediate-degree g modes in some MS massive stars.

Some regions of the sHR diagram remain, however, completely stable when considering ℓ = 1 to 20 modes: the post-MS of models with 3–8 M_⊙, the region close to the ZAMS of 9–70 M_⊙ models, and the much evolved models of 9–70 M_⊙. In particular, for this latter group, it is important to remark that numerous modes are found positive to an instability for models with M> 20 M_⊙ and with log T_eff< 4.3. The eigenfunctions of these modes are, however, sometimes ambiguous and more investigation needs to be performed in order to classify the types of modes encountered (including strange modes⁵) and to identify the excitation mechanisms. Therefore, although we have computed the instabilities until log T_eff = 4.0 and found excited modes in that region, we have decided to restrain the instability domains to the modes excited by the κ-mechanism due to the iron opacity bump only. We note that the spectra of the excited modes encountered in these models were usually less dense than for the earlier evolution (various modes were found in all the frequency range, but less concentrated in a given domain of frequency).

Finally, we remark that the frequency spectrum for lower masses reaches lower-degree modes for decreasing masses at the same evolutionary state (see models J, H to K), and that no modes with ℓ> 1 are predicted for 3 M_⊙. We discuss the global behavior of the high-degree modes when the stars evolve from the MS to the post-MS in the next section.

4.3. Properties of high-degree modes

A global look at modal diagrams presented in Figs. 4 and 5 allows us to characterize the behavior of the high-degree pulsations with the stellar evolution for different initial masses. For M ≳ 9 M_⊙, no pulsations are found close to the ZAMS. Modes with ℓ = 1–10 and high frequencies (i.e., having a p-mode behavior) are the first to appear at some point during the MS: the earlier it is, the more massive the star. Then, as the stars evolve towards lower T_eff in the sHR diagram (i.e., looking from model A to C for higher masses and from model D to F for 15 M_⊙), a large group of modes with lower frequencies (g-mode behavior) appears with a large range of degree ℓ, and the number of excited p modes diminishes (often for all degrees at the same time). As a general trend, the degree ℓ of g modes decreases along the evolution, even reaching ℓ = 1 in the post-MS region. Lastly, the red part of the sHR diagram is characterized by other types of modes which are not taken into account here, as already mentioned in Sect. 4.2.

Concerning lower masses (models G to I), the ZAMS is characterized by SPB type modes with high degrees (modes with ℓ up to 20 are excited, see model G and H). The degree of the mode, along with the frequency range of the excited modes, then decreases over the MS evolution. We note that stars with masses between ~7–9 M_⊙ predicted to pulsate with g modes, only present modes of degrees ℓ larger than three. Finally, models J, H, and K, selected at a given state of evolution (X_c ~ 0.3, Fig. 5), all present SPB-type g modes. At the same evolutionary stage, the largest mass models present high-degree excited modes with periods of the order of hours (models J and H present excited modes with ℓ up to 20). By decreasing the mass of the models (from J to K), we find modes excited with a decreasing ℓ, and the frequency range of instability reaches the usual long periods of the SPB stars observed in photometry, of the order of days (model I). The post-MS phase of evolution is then characterized by pure stable models.

Fig. 6 Perturbed temperature (top panel) and work integral (bottom panel) for three modes selected in model G: one excited mode (see the positive work integral at the surface on the bottom panel) with ω = 1.26 and ℓ = 10 (mode (1); blue solid line) and two stable quadrupole modes with ω = 0.29 (mode (2); dotted orange line) and ω = 1.27 (mode (3), dot-dashed green line).

While high-frequency modes appear and disappear similarly altogether during the evolution, we note two general trends concerning the modes with the lowest frequencies (“g modes”) for increasing ℓ:

• the range of unstable g modes is moving towards higher frequencies;

• the blue boundary of the instability domain of g modes is shifted towards warmer T_eff, overlapping the instability of β Cephei stars.

These two trends can be explained with the following reasoning, based on the two conditions needed for an efficient κ-mechanism. On the one hand, in order to have an efficient κ-mechanism, the amplitudes of the g modes have to be large enough in the opacity bump compared to the g cavity where the radiative damping plays an important role. Therefore, since the amplitudes of the eigenfunctions for g modes are essentially a function of the ratio ω²/ (ℓ(ℓ + 1)) (in the asymptotic limit, a fixed ω²/ (ℓ(ℓ + 1)) corresponds to a given radial order), the condition of the amplitudes implies therefore a given range for the ω²/ (ℓ(ℓ + 1)) ratio. On the other hand, an efficient κ-mechanism also requires that the transition zone (the region in the star where the pulsation period is of the same order as the thermal relaxation time) coincides with the region of the opacity bump. This is the case for a given range of frequencies, independently of the degree ℓ of the mode. These two conditions generate therefore two constraints, one depending on the degree of the mode and the other not. The intersection between these constraints defines the group of unstable g modes in the modal diagrams. To illustrate this effect, we selected three modes in model G:

• (1) a mode with a large degree (ℓ = 10) with f = 5.67 c/d (dimensionless frequency ω = 1.26) which has been selected from amongst the group of unstable modes,

• (2) a stable mode with a low degree (ℓ = 2), and with a lower frequency (f=1.29 c/d and ω = 0.29), which has been selected to fit the same ω²/ (ℓ(ℓ + 1)) ratio as mode (1), and

• (3) a stable quadrupole mode with the same frequency as mode (1) (f=5.70 c/d and ω = 1.27).

Figure 6 shows the eigenfunction of the perturbed temperature (top panel) and the work integral (bottom panel) for these three modes. From Fig. 6, we find that mode (1) presents a large amplitude in the superficial layers and a positive work at the surface: it is excited. However, for a lower-degree mode with the same frequency (as is the case for mode (3)), the eigenfunction takes low values in the superficial layers due to the larger evanescent zone (Dupret et al. 2008), and because of the boundary condition: δP/P = (l(l + 1) /ω²−4−ω²) δr/R (Buta & Smith 1979). As a result, the mode stays stable. Moreover, for a given radial order (i.e., for a fixed ω²/ (ℓ(ℓ + 1)), see modes (1) and (2)), the eigenfunctions of the perturbed temperature present large amplitudes in the superficial layers and, especially, in the region of the opacity bump, due to the iron elements, favoring an efficient κ-mechanism. However, unlike mode (1), mode (2) is found to be stable. Indeed, its frequency being lower, the transition region corresponding to mode (2) is situated deeper into the star. As a result, it is located outside the opacity bump of iron (characterized by a fixed log T), preventing the excitation of the associated mode. In short, the layout of the unstable low-frequency modes in the modal diagrams depends on two conditions controlling the efficiency of the κ-mechanism, that is, the constraint on the amplitudes of the modes depending on ℓ, and the constraint on the location of the transition region, which depends on the frequency itself only: a larger ℓ corresponds therefore to bluer instability domains and to excited modes of higher frequencies.

Before comparing our theoretical predictions with spectroscopic observations in the next section, we remark that the computations performed here are based on a given set of input physics (described in Sect. 2) which affects the exact boundary of the instability domains in the sHR diagram. We will refer to this in more detail in the context of macroturbulence in Sect. 5.

5. Link with observations: macroturbulent broadening

The most commonly considered scenario to explain the existence of a non-rotational line-broadening mechanism – traditionally quoted as macroturbulent broadening – shaping the line profiles of O stars and B Supergiants claims that this broadening is a signature of the collective effect of stellar oscillations (Lucy 1976; Aerts et al. 2009). While there are some empirical hints supporting this hypothesis (Simón-Díaz et al. 2010, 2016; Aerts et al. 2014), the exact physical mechanism driving the type of instabilities required to produce the observed line profiles is still under debate. One obvious possibility refers to heat-driven modes produced by a κ-mechanism. However, other options such as stochastically-excited oscillations driven by turbulent pressure fluctuations initiated either in the sub-surface convective zone (Grassitelli et al. 2015) or by the interface between the convective core and the radiative envelope (Rogers et al. 2013; Aerts & Rogers 2015) have been also recently proposed as potential hypotheses. Actually, as noted by Simón-Díaz et al. (2016), different excitation mechanisms could be contributing, at the same time and to some degree, to the global broadening of the line profiles depending on the mass and the evolutionary stage of the star.

In this section we connect the observations gathered in the framework of the IACOB project (Simón-Díaz et al. 2015) with our predictions for instability domains of high-degree modes driven by the κ-mechanism, in order to investigate the physical origin of the so-called macroturbulent broadening in the O and B star domain. This is a continuation of the work started in Simón-Díaz et al. (2016), where we compared the empirical information extracted from a sample of high-resolution single-snapshot spectra of ≈430 Galactic O and B stars, to the predictions for the dipole instability domains presented in this paper. In the upper part of the sHR diagram, we found an important number of stars which were significantly affected by macroturbulent broadening but were located outside the g-mode instability domains. Furthermore, the O stars and B Sgs located inside the instability domains, populate different regions in the diagram characterized by different type of excited modes. However, a quite homogeneous distribution of profile types (in terms of global shape and relative contribution between the macroturbulent and the rotational broadening) is found for all these stars, despite the expected different effect on the line profiles due to the specific combination of excited modes. All together, these results could be interpreted as a strong empirical challenge to macroturbulent broadening in O stars and B Sgs being a spectroscopic signature of only non-radial modes excited by the κ-mechanism. However, as commented in Simón-Díaz et al. (2016), the situation may improve when including the predictions for high-degree modes or when modifying the input parameters in the stellar model computations, such as the opacities and the metal mixture. These two possibilities are discussed in detail in the next sections.

5.1. High-degree modes vs. macroturbulent broadening

Hereafter, we work under the hypothesis that the bare minimum requirement to end up in a macroturbulent-type profile produced by pulsation modes is a dense frequency spectrum of excited modes (e.g., Aerts et al. 2009). In the context of the κ-mechanism, such spectra can be achieved by the presence of high-order g modes – since they are usually excited in an abundant spectrum – or even by combining the collective effect of a group of modes with a large range of ℓ degrees. Assuming this requirement, we can determine the best models resulting in such a situation by combining the information presented in Fig. 3 (instability domains) and Figs. 4 and 5 (modal diagrams). Basically, stars with pulsational properties similar to those of models E–I, B, and C, are good candidates to produce an important broadening of the line profiles (unlike models A and D). Models E–I, B, and C correspond to regions in Fig. 3 colored in red and purple (low-frequency modes) while models A and D refer to regions colored in blue (high-frequency modes). The region where the existence of an important contribution of the macroturbulent broadening could be explained in terms of high-degree heat-driven modes is, therefore, more extended than when only dipole modes are considered. Given the complexity of the interpretation of the results from the pulsational analysis in the coolest post-MS region of the most massive models, we have completely excluded any representation of the instabilities found in that region (see Sect. 3). However, although some excited high-degree modes are found in this region, the associated frequency spectrum is sparse.

Fig. 7 Location in the sHR diagram of the observational dataset presented in Simón-Díaz et al. (2016): red points stand for the stars presenting an important (or dominant) macroturbulent broadening contribution and green points correspond to stars for which either rotational broadening dominates or both vsini and vmac are below 25 km s-1 (see Simón-Díaz et al. 2016, for further details). The size of the points is linearly scaled with vmac, ranging from 10 to 130 km s-1. Instability domain predictions (as described in Sects. 4.1 and 4.2) for 1 ≤ ℓ ≤ 20 are colored in blue and red.

Having this in mind, we compare in Fig. 7 the empirical distribution of the line-broadening properties of the sample of ≈260 Galactic O and B stars⁶ considered by Simón-Díaz et al. (2016), and our instability domain predictions for heat-driven modes with 1 <ℓ ≤ 20. The sample is divided into two main groups depending on their global line-broadening characteristics, and the size of the symbols is scaled to v_mac, the amount of macroturbulent broadening. The colors used for the instability domains have the same meaning as in Fig. 3.

As commented in Simón-Díaz et al. (2016), most of the stars with a remarkable contribution of the macroturbulent broadening (red points) are concentrated in the uppermost part of the sHR diagram, and the largest values of v_mac are measured for the late-O/early-B Sgs (log L/ L_⊙ ≈ 4.0, and log T_eff ≈ 4.55−4.45). Many of these stars are located in regions where heat-driven modes could potentially explain the observed line profiles. This refers, in particular, to those targets with M ≳ 10 M_⊙ falling inside the instability domain colored in red (corresponding to high-order g modes with a relatively dense global frequency spectrum).

While the situation is improved when including the full range of mode degrees up to ℓ = 20 compared to the case of only considering dipole modes, some problems remain. Indeed, there are still two well-populated regions in Fig. 7 where stars are clearly located outside the instability domains (even taking individual uncertainties into account) but have line profiles dominated by macroturbulent broadening:

• the first region refers to the late-B supergiants. We note, however,that instabilities were found in that region, although they are notdirectly related to the κ-mechanism due to the iron opacity bump (see Sect. 4.2). As a matter of fact, the ϵ-mechanism could play an important role in post-MS models by exciting some high-order g modes which could be trapped inside the H-burning shell (Scuflaire & Noels 1986; Moravveji et al. 2012). Furthermore, in addition to these modes excited by the ϵ-mechanism, adiabatic and non-adiabatic strange modes are also expected for all masses larger than ≈30 M_⊙ (above the dotted line in Fig. 7), but they do not usually lead to a dense frequency spectrum.

• the second region refers to the O stars close to the ZAMS with masses ≈15–40 M_⊙. This latter point could be used, at first glance, as a strong argument to dismiss the possibility that heat-driven modes are the main source of macroturbulent broadening in the O-star and B-supergiant domain. However, before reaching any firm conclusion in this direction, we must explore the effect of considering different assumptions made for the input parameters of the computations on the predicted boundaries of the instability domains.

Finally, we note that there is one region where almost no macroturbulent broadening is found though it is located inside the theoretical instability domain: the low-mass MS stars lying in the SPB instability domain. For this latter group of stars, it is important to notice that, from the observational sample of ≈430 stars presented in Simón-Díaz et al. (2016), we discarded, in the representation of Fig. 7, stars for which it was only possible to derive a maximum limit on the macroturbulence parameter. When taking these stars into account, we find that most of them are located at the beginning of the MS of 6–9 M_⊙ models (where the theoretical computations predict a sparse spectrum of low-order p and g modes), while rather few stars in the theoretical range 4–6 M_⊙ were observed (where high-order g modes are expected).

From these results, it appears that a crucial point for the investigation of the origin of macroturbulence broadening in line profiles resides in the presence (or absence) of 12–40 M_⊙ ZAMS stars with the required pulsational properties. Therefore, we focus in the next section on the efforts needed to excite modes at the beginning of the MS.

5.2. Bringing instability domains towards the ZAMS

As just concluded, any mechanism which aims to produce a macroturbulent-type profile in 12–40 M_⊙ stars must be efficient in their early stage of evolution. As illustrated in the literature, the use of a different set of input parameters (e.g., overshooting parameter, opacities, metal mixture, metallicity) or even evolutionary codes may alter the results of our computations (see, e.g., Miglio et al. 2007a,b; Pamyatnykh 2007; Zdravkov & Pamyatnykh 2008; Martins & Palacios 2013; Castro et al. 2014; Moravveji 2016, for recent illustrations of these effects). The question is thus the following: how can we widen the instability strip towards the ZAMS? While the overshooting parameter will be of no help to stretch the blue boundary, the metallicity, the detailed chemical mixture (e.g., Asplund et al.2009, hereafter AGSS09 vs. Grevesse & Noels1993, hereafter GN93), and opacities are obviously to play a major role in the efficiency of the κ-mechanism (Miglio et al. 2007a; Pamyatnykh 2007; Moravveji 2016).

The computations presented in this paper were performed with the GN93 metal mixture and adopting a metallicity of 0.015, in agreement with AGSS09. The choice of this combination was motivated by the fact that the AGSS09 solar mixture was not yet implemented in the ATON evolutionary code. In order to investigate the effect of the use of a different metal mixture (AGSS09 vs. GN93), metallicity (0.015 vs. 0.020), and opacities (OP vs. OPAL), we started the computation of new models in the 7 to 40 M_⊙ range combining the CLES evolution code from Liège (Scuflaire et al. 2008b, for which both opacity tables and both metal mixtures are implemented) and the non-adiabatic code MAD. The models computed with OP opacities are limited to masses up to 25 M_⊙ due to the narrower temperature and density range in the OP tables. Except for the three sets of parameters indicated above, the rest of the input physics were fixed to those parameters used for our ATON models. While the presentation of these new (still on-going) computations will be the subject of a more extensive study in the future, we focus here, for illustrative purposes, on the predictions for the blue boundary of the p-mode (ω = 1–25), ℓ = 1 instability domain.

Fig. 8 Blue boundary of the instability domains for dipole modes. The blue area limited by a dotted line represents the instability domains already presented (ATON evolutionary code, GN93, Z = 0.015) while green and purple lines stand for the blue boundaries of the instability domains for CLES models computed with the AGSS09 mixture and different metallicities and opacity tables (OPAL and OP, respectively). Green and red points have the same meaning as in Fig. 7 and their sizes are proportional to the macroturbulence velocity.

Figure 8 summarizes the results of these new computations. This figure illustrates how critical is the selection of the three input parameters mentioned above for reaching a definitive conclusion about the heat-driven mode scenario. First, the comparison between the ATON and CLES computations with Z = 0.015 and OPAL allows us to illustrate the effect of the metal mixture. Although no modes are found close enough to the ZAMS yet, the AGSS09 boundary is clearly shifted towards higher T_eff (compared to the GN93 one), and the associated instability strip now gathers some more O stars presenting a dominant macroturbulent broadening contribution. This result is encouraging and brings us one step further. Indeed, the situation improves even more (for a fixed metallicity) when considering the OP opacities (as already pointed out by Miglio et al. 2007a and Pamyatnykh 2007, the combination of AGSS09+OP favors excitation through the κ-mechanism due to the iron opacity bump, much more than the combination of GN93+OPAL). Lastly, the effect of metallicity is clear: the larger the metallicity, the broader the instability domain.

These effects (increasing metallicity and opacity) both broaden the instability region in our calculations and thus help to support the pulsational hypothesis as being the origin of macroturbulent broadening. However, these results must be taken with caution. Firstly, given the most recent value of metallicity (Z = 0.014) derived from the spectroscopic analysis of a sample of ≈30 early B-type stars in the solar vicinity by Nieva & Przybilla (2012, Z = 0.014, higher metallicity (Z = 0.02) favoring a wider instability domain may not be representative of massive stars. Secondly, concerning the opacities, we highlight the recent studies of Salmon et al. (2012) and Turck-Chièze & Gilles (2013), who suggest the need to reassess opacity computations, in particular for nickel, which were shown to have an important effect on the κ-mechanism instability predictions in massive stars (see also Mondet et al. 2015, who revealed a larger mean Rosseland opacity in their OPAS computations of 11% – compared to OPAL and OP – around log T ≃ 6.2 in solar modeling). Moreover, the recent work of Walczak et al. (2015) has shown that the newly developed OPLIB opacity tables from Los Alamos National Laboratory induce an even higher Rosseland mean opacity coefficient at the maximum of the iron opacity bump (see Figs. 1 and 3 of their paper). Lastly, Bailey et al. (2015) found a large disagreement between the experimental measurement of iron opacity and the theoretical computations, for conditions close to that of the base of the solar convective zone. Based on this, Moravveji (2016) showed that a 75% increase in Fe and Ni monochromatic opacities from OP tables extends the β Cephei domain to the ZAMS, in agreement with our results without any metal enhancement presented in Fig. 8. As a result, while the metallicity cannot play a major role in widening the instability strip since the value fixed in the computations already corresponds to the commonly accepted metallicity value for massive stars, the computations of the opacities constitute, therefore, the crucial point to help widen the instability domains towards the ZAMS.

To conclude, our investigation of the instability domains associated with high-degree modes, in combination with the results presented in Fig. 8 and all the studies on the opacity quoted above, show that the heat-driven mode scenario to explain macroturbulent broadening in O and B stars may not be in as bad a situation as when only considering the dipole modes. Of particular importance is the fact that certain assumptions on the opacities can widen the instability domains towards the ZAMS, allowing some of the O stars with macroturbulent-broadening-dominated profiles to lie inside. However, there are still some clues which make us suspect that the heat-driven pulsation modes due to the κ-mechanism in the iron bump cannot explain alone the distribution of line-broadening properties of the line profiles in massive stars. First, in view of panels A and D in Fig. 4, even assuming the opacities resulting from the most favorable conditions, it appears that a dense spectrum of pulsation modes would be very difficult to produce at the beginning of the MS. Second, even if we assume that models similar to those labeled with B and E could be found close enough to the ZAMS, we would need to explain why there is a non-negligible number of narrow line (i.e., with a negligible macroturbulent broadening contribution) MS stars with masses in the range ≈9–20 M_⊙.

We briefly discuss in the next section some additional hypotheses to attempt to explain the origin of the non-rotational extra-broadening in OB stars.

5.3. Which additional driving scenario?

In addition to the spectroscopic signature of the macroturbulent-type profile, direct observations also indicate the presence of pulsations in stars outside or close to the blue instability boundary for masses between 15–40 M_⊙. Indeed, the presence of line-profile variations has been demonstrated to be a common feature among O-type stars, even close to the ZAMS (e.g., Fullerton et al. 1996). For some of them, it has been possible to link this variability with pulsations (e.g., Blomme et al. 2011; Briquet et al. 2011; Mahy et al. 2011). Moreover, there is increasing observational (Degroote et al. 2010) and theoretical evidence (Belkacem et al. 2010; Samadi et al. 2010) for the presence of stochastic oscillations. Furthermore, strange modes (e.g., Noels et al. 2008; Glatzel 2009; Saio 2009), oscillatory convective modes (Saio 2011), and pulsation modes due to the ϵ-mechanism (Scuflaire & Noels 1986), have been suggested to be the explanation of photometric observations (Aerts et al. 2010a; Moravveji et al. 2012). All these mechanisms may help in producing a macroturbulent-type profile, though they usually produce a scarce spectrum of unstable modes. While several other scenarios have been suggested for the stellar oscillation origin, the most viable to this date is stochastically-driven waves induced by turbulent pressure fluctuations in the interface between the convective core and the envelope (Aerts & Rogers 2015) or in the sub-surface convective zone (Grassitelli et al. 2015). In particular, Aerts & Rogers (2015) demonstrated with CoRoT space photometry the observational evidence for the occurrence of convectively driven internal gravity waves in young massive O-type stars giving rise to a macroturbulent-type line profile. Grassitelli et al. (2015) have computed and mapped the turbulent pressure associated to sub-surface convection zones in massive MS and post-MS models. The authors presented a good correlation between ${\mathit{P}}_{\mathrm{turb}}^{\max }\mathit{/}\mathit{P}$ and v_mac. As indicated by Simón-Díaz et al. (2016), the interplay between the turbulent pressure fluctuation hypothesis and the heat-driven modes due to the κ-mechanism scenario seems like a promising explanation for the origin of macroturbulent broadening in the upper part of the sHR diagram.

6. Summary and concluding remarks

In this paper we investigated the pulsation modes driven by the κ-mechanism in the iron opacity bump in massive stars. One of the goals of this paper was to provide new predictions for the instability domains of non-radial modes driven by the κ-mechanism in the opacity bump of iron. The originality of our results comes from the use of a large domain of ℓ degrees for a large range of masses between 3 and 70 M_⊙. Although these modes are not visible with photometric detections, we think our results will be a valuable tool for constraining the new high-quality spectroscopic observations in massive OB stars, as for example, the observational material gathered by the IACOB project. In particular, we provided a global overview of the frequency ranges and characteristics of the modes to expect in the main sequence and the post-main sequence of these massive models. On the main sequence, a few high-frequency modes, with periods around eight hours, are found to be unstable in our theoretical models. These periods increase during the evolution from the main sequence to the post-main sequence, but the frequency spectrum is always characterized by a small number of degree ℓ modes, with ℓ never larger than ten. On the post-main sequence (and even on the main-sequence for the smallest masses), low-frequency modes are excited in a much more abundant spectrum. These modes have periods ranging between a couple of hours to a dozen days. We found that the presence of the low-frequency modes in the evolved models which do not present any intermediate convective zone (appearing for models ≥12 M_⊙), is due to the sharp feature already present in the Brunt-Väisälä frequency in models of 9 M_⊙, in agreement with Daszyńska-Daszkiewicz et al. (2013). We discussed the appearance of these low-frequency unstable modes by considering the constraints on the efficiency of the κ-mechanism. Actually, we found that for an increasing ℓ, the range of unstable low-frequency modes moves towards higher frequencies. Furthemore, by increasing the degree of the modes, the blue boundary of the instability domains for the low-frequency modes is shifted towards the zero-age main sequence, even overlapping the high-frequency instability domain (β Cephei instability strip).

In addition to the usual spectroscopic observations, the emergence in the last few years of the investigation of the so-called macroturbulent broadening has provided a new way to study stellar oscillations which is less expensive in terms of observational time. While the strongest hypothesis for the origin of this macroturbulent-type profile is to relate it to stellar oscillations (e.g., Simón-Díaz et al. 2016), several mechanisms and types of waves and/or modes driving the macroturbulence have been suggested. In this paper we have shown that pulsation modes driven by the κ-mechanism in the opacity bump of iron cannot explain alone the distribution of line-broadening properties of the line profiles in the whole OB star domain, and hence they might not be the dominant agent producing the so-called macroturbulent broadening in O stars and B Supergiants. Indeed, although we considered numerous spectra of modes combining different ℓ degrees, the location of the instability domains for this kind of spectra do not entirely coincide with the location of the high v_mac stars. We investigated how the physics of the models have to be adapted in order to provide such an excitation in the beginning of the main sequence, as we think it is the key point for constraining the origin of the mechanism which is broadening the line profiles. While the major outcome would be achieved considering new computations of the opacities, it is still difficult to obtain such an abundant spectrum close to the zero-age main sequence. Indeed, for a reasonable metallicity in massive stars (Z ~ 0.015), models close to the beginning of the main sequence are found unstable when considering the OP opacities, although the frequency spectrum of these models is scarce. Another additional mechanism should therefore trigger the broadening of the line profiles and it must be efficient in the early evolution of the star, close to the beginning of the main sequence (such as the scenarios of Grassitelli et al. 2015; Aerts & Rogers 2015). Already foreseen by Simón-Díaz et al. (2016), we agree therefore with the idea that the line broadening must be due to a combination of different sources depending on the regions in the HR diagram.

We propose several directions for future investigations: (a) the possibility of extending the instability boundaries for stars with M ≥ 12 M_⊙ towards the zero-age main sequence, with a reasonable enhancement of opacities (especially nickel) although the calculations computed so far in the litterature and in this paper have shown that few modes are excited in these regions, (b) the study of other driving mechanisms for stellar oscillations not accounted for in the present computations, for example, theoretical simulations to confirm that the turbulent pressure fluctuations investigated by Grassitelli et al. (2015) do indeed end up in a macroturbulent profile: in particular, if an excitation mechanism such as pressure fluctuations (which in 1D models correspond to a radial velocity pattern) can give rise to macroturbulent velocity distributions as used here, with a strong tangential component and (c) the extension of the type of comparison between observations and models presented in this paper but including information about line-profile variability as extracted from adequate time-resolved spectroscopic observations (see, e.g., first results in this direction in Simón-Díaz et al. 2016).

Acknowledgments

This work has been funded by the Spanish Ministry of Economy and Competitiveness (MINECO) under grants AYA2010-21697-C05-04, AYA2012-39364-C02-01, and Severo Ochoa SEV-2011-0187, and by the Research Council of KU Leuven under grant GOA/2013/012. We thank C. Aerts, P. Degroote, J. Montalban, E. Moravveji, J. Puls, and J. Telting for interesting and fruitful discussions related to the work.

References

All Tables

All Figures

Fig. 1 Left: modal diagram: range of dimensionless frequencies for a 50 M⊙ MS model (labeled B in Table 1) as a function of the degree ℓ of the mode. Stable modes are represented by gray points while excited modes are shown with solid blue circles. Right: propagation diagram for the same model. The solid red line stands for the Brunt-Väisälä frequency, and the Lamb frequency for three different degrees (ℓ = 1, 3 and 20) is represented by the dotted, dashed, and dot-dashed lines, respectively. The propagation zones for p and g modes with ℓ = 1 are highlighted by the dark gray and light blue shadow regions, respectively. The solid horizontal black line stands for an example of the propagation zones for a mode of ω = 1.0 with ℓ = 3.

In the text

Fig. 2 Instability domain predictions in the sHR diagram for dipole modes (ℓ = 1, see Sect. 4.1). P-mode instabilities are represented by the blue area and the g-mode domain is shown in red. The ATON non-rotating evolutionary tracks are over plotted. The horizontal dashed line indicates the limit above which strange mode oscillations are most likely expected. Black asterisks stand for nine selected models la A to I which will be more deeply investigated in Sect. 4.2. The global properties of these models are given in Table 1.

In the text

	Fig. 3 Same symbol key as in Fig. 2, except that the instability domain predictions are now done for 2 ≤ ℓ ≤ 20 pulsation modes (see Sect. 4.2).
In the text

Fig. 4 Range of frequencies in cycle per day (c/d) for models A to I. Gray dots stand for the frequencies of the stable modes computed, whereas excited mode frequencies are represented by blue bigger dots. From left to right, the rows represent the evolution from the MS to the post-MS: models A to C for high masses (30–50 M⊙), models D to F for 15 M⊙, and models G to I for 6 M⊙ MS models only (see Fig. 2 and Table 1). Rough evolutionary states are indicated for each model in the right bottom corner (ZAMS, EMS, MMS, TAMS, and post-MS stand for zero age MS, early MS, middle of MS, terminal age MS, and hydrogen-burning shell models). We note that the detectability of the pulsation modes in photometry is limited to the low degree (ℓ up to 3–4), whereas higher degree modes should be able to be observed thanks to spectroscopy.

In the text

	Fig. 5 Same symbol key as in Fig. 4 but for models J to K, selected, from top to bottom, at the same evolutionary state (Xc ~ 0.3) but with decreasing masses from 8 to 4 M⊙ (see Fig. 2 and Table 1).
In the text

	Fig. 6 Perturbed temperature (top panel) and work integral (bottom panel) for three modes selected in model G: one excited mode (see the positive work integral at the surface on the bottom panel) with ω = 1.26 and ℓ = 10 (mode (1); blue solid line) and two stable quadrupole modes with ω = 0.29 (mode (2); dotted orange line) and ω = 1.27 (mode (3), dot-dashed green line).
In the text

Fig. 7 Location in the sHR diagram of the observational dataset presented in Simón-Díaz et al. (2016): red points stand for the stars presenting an important (or dominant) macroturbulent broadening contribution and green points correspond to stars for which either rotational broadening dominates or both vsini and vmac are below 25 km s-1 (see Simón-Díaz et al. 2016, for further details). The size of the points is linearly scaled with vmac, ranging from 10 to 130 km s-1. Instability domain predictions (as described in Sects. 4.1 and 4.2) for 1 ≤ ℓ ≤ 20 are colored in blue and red.

In the text

Fig. 8 Blue boundary of the instability domains for dipole modes. The blue area limited by a dotted line represents the instability domains already presented (ATON evolutionary code, GN93, Z = 0.015) while green and purple lines stand for the blue boundaries of the instability domains for CLES models computed with the AGSS09 mixture and different metallicities and opacity tables (OPAL and OP, respectively). Green and red points have the same meaning as in Fig. 7 and their sizes are proportional to the macroturbulence velocity.

In the text