OBSERVATIONAL CONSEQUENCES OF TURBULENT PRESSURE IN THE ENVELOPES OF MASSIVE STARS

Massive stars are of key importance for the enrichment of the interstellar medium (Kobayashi et al. 2006), for regulating star formation (Mac Low & Klessen 2004), and as progenitors of supernovae and gamma-ray bursts (Langer 2012). However, various physical processes acting in massive stars are as yet not well understood, therefore preventing the current stellar models to self-consistently explain several general observational properties, such as the mass-discrepancy (Herrero et al. 1992; Markova & Puls 2015), the presence of the Humphreys–Davidson limit (Humphreys & Davidson 1979), and the position of the terminal-age main sequence in the HR diagram (e.g., Vink et al. 2010; Castro et al. 2014). Radiation-pressure-dominated layers, envelope inflation, clumpy stellar winds, and dynamical instabilities are features that commonly occur in these objects, but that have not yet been thoroughly investigated. Some of these phenomena have been connected to the presence of turbulent motions in the outer layers (e.g., de Jager 1984; Cantiello et al. 2009). Despite its simplicity, the mixing-length theory (MLT; Böhm-Vitense 1958) for convection has been very successful in describing the main features of non-adiabatic turbulent stellar envelopes (e.g., Trampedach et al. 2014). However, the MLT is expected to have shortcomings when the convective velocities approach the speed of sound, as it is expected in the envelopes of very luminous stars (Sanyal et al. 2015). In particular, it has been argued by de Jager (1984) and Maeder (2009) that turbulent pressure may become important in this situation. It is the aim of the present paper to discuss the role of turbulent pressure in the envelopes of massive stars, investigate its effects on the stellar structure, and focus on its potential observational signatures.

We use the Lagrangian one-dimensional hydrodynamic stellar evolution code BEC (Heger et al. 2000; Petrovic et al. 2005; Brott et al. 2011) for computing massive star models. It treats convection following the MLT with a mixing-length parameter of $\alpha =1.5$ (Brott et al. 2011). The opacity of stellar matter is interpolated from the OPAL opacity tables (Iglesias & Rogers 1996), and mass loss by stellar wind is adopted following the prescription of Vink et al. (2001).

We modified the BEC code in order to include the turbulent pressure term in the stellar envelopes, in particular, in the convective zones associated with the iron opacity peak (FeCZ), and with the partial ionization of helium (HeCZ) and hydrogen (HCZ). We neglect the turbulent pressure in the convective core, where it is typically more than seven orders of magnitude smaller than the ideal gas pressure. We included the turbulent pressure ${P}_{\mathrm{turb}}$ and the turbulent energy density ${e}_{\mathrm{turb}}$ in the momentum, energy transport, and energy conservation equations, respectively, following (Canuto & Mazzitelli 1991; Jiang & Huang 1997; Stothers 2003; Maeder 2009; Trampedach et al. 2014)

$$\begin{eqnarray}&&{P}_{\mathrm{turb}}=\zeta \rho {v}_{{\rm{c}}}^{2},\quad {e}_{\mathrm{turb}}=\displaystyle \frac{3}{2}\displaystyle \frac{{P}_{\mathrm{turb}}}{\rho },\end{eqnarray} \tag{ 1 }$$

where ρ is the local mass density, ζ is a parameter chosen to be $\zeta =1/3$ for isotropic turbulence (Stothers 2003; Maeder 2009), and v_c is the local convective velocity. As the convective velocities approach the local sound speed, we limit the convective velocities to this value, i.e.,

$$\begin{eqnarray}&&{v}_{{\rm{c}}}\leqslant {c}_{{\rm{s}}},\quad {c}_{{\rm{s}}}^{2}=\displaystyle \frac{{k}_{{\rm{B}}}T}{\mu {m}_{{\rm{H}}}}=\displaystyle \frac{{P}_{\mathrm{gas}}}{\rho },\end{eqnarray} \tag{ 2 }$$

where T is the local temperature, k_b is the Boltzmann constant, μ is the mean molecular weight, and m_H is the proton mass. Here, we use the isothermal sound speed, as the layers under consideration are characterized by a small ratio of the local thermal-to-dynamical timescale (${\tau }_{\mathrm{th}}/{\tau }_{\mathrm{dyn}}\lt 1$). As a consequence of Equation (2), the turbulent pressure in our models can not exceed a value of one-third of the local gas pressure.

We computed a set of stellar models and stellar evolutionary tracks for non-rotating stars with Galactic metallicity in the initial mass range 7–80 ${M}_{\odot }$, including the physics described in Section 2. All input and physics parameters were chosen as in Brott et al. (2011).

3.1. Evolutionary Tracks

The evolutionary tracks of our newly computed models are plotted in Figure 1 and superposed to those obtained by Brott et al. (2011). Throughout most of the evolution, the effect of turbulent pressure and energy on the luminosity and the surface temperature of the models is small, such that our new evolutionary tracks nearly coincide with those from Brott et al. (2011). According to the color scheme in Figure 1, which indicates the maximum fraction of turbulent pressure inside the stellar models, the HR diagram can be roughly divided into three areas: the hot (${T}_{\mathrm{eff}}\gt {10}^{4}$ K) and luminous ($\mathrm{log}(L/L\odot )\gt 4.5$) stars, the cool stars (${T}_{\mathrm{eff}}\lt {10}^{4}$ K) of any luminosity, and the hot stars with $\mathrm{log}(L/L\odot )\lt 4.5$.

Zoom In Zoom Out Reset image size

In the latter regime, the turbulent pressure does not exceed a fraction of the total pressure of a few tenths of a percent anywhere inside the stellar models. On one hand, the iron opacity peak in these models is located deep inside the envelope in nearly adiabatic layers, which leads to relatively low convection velocities. On the other hand, the models are too hot to contain a hydrogen ionization zone. Finally, they do show a He II partial ionization zone close to the surface, which is, however, not vigorous enough to play a significant role (Cantiello et al. 2009).

In the cool supergiant region, the turbulent pressure fraction inside the star can be very significant. The maximum value is typically 25%–30%, while the highest possible value, 33%, is reached for the most luminous models. The computed transonic convective velocities, of the order of 20–30 km s⁻¹, arise in the HCZ very close to the surface, where high opacities induce high local Eddington factors, density inversions, and high degrees of superadiabaticity (Sanyal et al. 2015). Such a high turbulent pressure fraction is possible here since ideal gas pressure is by far the dominant contribution to the total pressure in the envelopes of the cool supergiants. Even though a significant fraction of the total pressure arises from the turbulent motion, the evolutionary tracks are only shifted by few tens of Kelvin toward lower temperatures, compared to the tracks without turbulent pressure, which corresponds to a radius increase by a few percent. This is due to the fact that the region that contains high turbulent pressure contains very little mass (see below).

In the hot and very luminous stars, the turbulent pressure can account for up to ∼5% of the total pressure. The evolutionary tracks are not significantly different from those of Brott et al. (2011), showing displacements of the order of hundreds of Kelvin toward lower effective temperatures. Turbulent pressure becomes more important as the stars expand during their main-sequence evolution, as the convective velocities increase. It achieves a maximum in the O supergiant regime, and then decreases in the B supergiant regime as the iron opacity peak moves deeper inside, but turbulent pressure remains significant for surface temperature above ∼10,000 K. The maximum turbulent pressure in these models occurs within the FeCZ, where high local Eddington factors are achieved, giving rise to envelope inflation and density inversions (Sanyal et al. 2015). It does not reach as high fractions as in the cool supergiant because of the predominance of radiation pressure in the hot star envelopes.

3.2. The Structure of the Outer Layers

Following the results shown in Section 3.1, we investigate the structure of the 20 ${M}_{\odot }$ and 60 ${M}_{\odot }$ models in more detail. Figures 2 and 3 show the relative fraction of the turbulent pressure as function of the optical depth in the stellar envelopes of both models throughout the evolution.

Zoom In Zoom Out Reset image size

In the 20 ${M}_{\odot }$ model (Figure 2), we find that for ${T}_{\mathrm{eff}}$ $\geqslant \;{10}^{4}\;{\rm{K}}$ the turbulent pressure fraction has its maximum in the FeCZ, located at an optical depth of $\;\mathrm{log}(\tau )\;$ $\approx \;2.5-4$. The turbulent pressure fraction is about 1.8% at ${T}_{\mathrm{eff}}$ $\simeq \;25,000\;$ K, and the FeCZ moves deeper inside the star as the stellar model expands during its evolution. Once the stellar model reaches effective temperatures well below 10,000 K, the HCZ arises. This convective region, within which the maximum turbulent pressure fraction arrives at about 25%–30%, reaches the stellar surface in stellar models with ${T}_{\mathrm{eff}}\lt 8000$ K and covers an extended range of optical depths. The bottom panel of Figure 2 compares the density profiles and the radial extent of the convection zones for a model with ${T}_{\mathrm{eff}}$ $\simeq \;7000\;{\rm{K}}$ and turbulent pressure included in the same model where turbulent pressure was disregarded.

In the envelope of the 60 ${M}_{\odot }$ model, we find convective zones associated with the partial ionization of iron and helium, while the model does not become cool enough to show hydrogen recombination.

As shown in Figure 3, the FeCZ is located at an optical depth of $\;\mathrm{log}(\tau )\;$ $\gt \;1.5$, with ${P}_{\mathrm{turb}}/P$ rising to ∼5% for ${T}_{\mathrm{eff}}$ $\simeq \;40,000\;$ K. As the star evolves, the convective region moves deeper inside the star. Figure 3 shows that due to the turbulent pressure, the star increases its radius by few percent, which leads to a slightly reduced density and an increase of the radial extent of the iron convection zone.

Zoom In Zoom Out Reset image size

In particular, the iron opacity peak at $T\;\approx \;$ 200,000 K is known to induce a variety of dynamical phenomena at the stellar surface, e.g., pulsations in β-Cephei and slowly pulsating B stars (SPB; Pamyatnykh 1999; Miglio et al. 2007), and stochastically excited traveling waves generated by turbulent motions in the FeCZ (Goldreich & Kumar 1990; Cantiello et al. 2009; Belkacem et al. 2010; Samadi et al. 2010; Mathis et al. 2014) leading to a subsonic small-scale velocity field at the surface that has been proposed to be the physical origin of the so-called "microturbulence."

The spectra of luminous OB stars are known to be also affected by the so-called macroturbulent broadening, an extra line-broadening, usually ad hoc, associated with large-scale (compared to the line-forming region) motions at the surface (see Markova et al. 2014, Simón-Díaz & Herrero 2014 and references therein). Similarly to the case of microturbulence, convection might play a significant role in the origin of macroturbulence as well. This view is supported by the work of Sundqvist et al. (2013), who showed that macroturbulence is suppressed in strongly magnetic massive stars, where the magnetic field is expected to at least partially inhibit convection. A similar effect was found for intermediate-mass magnetic chemically peculiar stars (Ryabchikova et al. 1997) and for spots in late-type stars (Strassmeier 2009). In this context, vigorous envelope convection in the temperature range of the hydrogen recombination (Figure 2) may be responsible for the non-thermal (macroturbulent) broadening observed in red supergiants (Collet et al. 2007; Carney et al. 2008).

With this motivation, we pursue the hypothesis that the relative strength of turbulent pressure in the sub-surface convective zones is related to the appearance and strength of macroturbulence at the stellar surface. We investigate the case of the luminous OB stars, where the turbulent pressure constitutes up to ≈5% of the total pressure in the FeCZ.

We make use of the results from the quantitative spectroscopic analysis of the rich sample of spectra compiled by the IACOB project (Simón-Díaz et al. 2011, 2015). In particular, we benefit from the derived values of surface temperature (${T}_{\mathrm{eff}}$), gravity (log(g)), projected rotational velocity (v sin(i)), and macroturbulent velocity (v_macro) for a sample of ∼300 Galactic OB stars used in Simón-Díaz (2015).⁴

In Figure 4, we compare the observationally derived macroturbulent velocities to the maximum fraction of turbulent pressure in our models in the spectroscopic HRD (sHR; Langer & Kudritzki 2014). Stars presenting a clear signature of macroturbulence in their line profiles (i.e., a v_macro larger than 50 km s⁻¹) are marked by a bigger black circle in Figure 4. Interestingly, those stars are located mainly in regions of the sHR where the turbulent pressure yields the highest contribution to the total pressure in the FeCZ. At places where the models predict a small contribution from the turbulent pressure, only very few stars show an unambiguously high macroturbulent velocity. We discuss these exceptions at low ${\mathcal{L}}$ in the next section.

Zoom In Zoom Out Reset image size

The agreement becomes even more striking looking at Figure 5, where the observed macroturbulent velocities are directly plotted against the maximum fraction of turbulent pressure in our models. This plot reveals a clear correlation between macroturbulent velocity and turbulent pressure starting from ${P}_{\mathrm{turb}}^{\mathrm{max}}/P\approx 0.005$, with a Spearmanʼs rank correlation coefficient of 0.812.

Zoom In Zoom Out Reset image size

Figure 4 revealed about 10 stars with log(${\mathcal{L}}/{{\mathcal{L}}}_{\odot })\lt 3.3$ for which our models predict very small turbulent pressure contributions, but which unambiguously show macroturbulence at a level above 50 km s⁻¹ (see also Figure 5). Interestingly, we find that all these stars are located inside or very near to the region where stars are expected to be pulsationally unstable to high-order g-modes (Pamyatnykh 1999; Miglio et al. 2007; see also Simón-Díaz 2015), which is drawn as a gray band in Figure 4. Indeed, Aerts et al. (2009) showed that the collective effect of high-order non-radial pulsations may produce a velocity field in the spectra of hot stars that resembles closely what is observed as macroturbulence in the stars discussed in Section 4.

Consequently, we interpret these stars as affected by a macroturbulent broadening that can be explained via a heat-driven pulsational origin. However, the homogeneity of the spectroscopic signatures of macroturbulence over the whole effective temperature range in the luminous stars calls for a single dominant mechanism to produce it (S. Simón-Díaz 2015, in preparation). This is not the case when considering only classical (κ-mechanism) instability domains, which do not cover the full region where most of the stars showing a macroturbulent velocity field are observed (see M. Godart et al. 2015, in preparation; Simón-Díaz 2015). Shiode et al. (2013) and Aerts & Rogers (2015) consider gravity waves originating in the convective core as the cause of macroturbulence. Whereas Shiode et al. (2013) find that the surface velocity fluctuations do not exceed 10 m s⁻¹ even in their most optimistic case, based on their massive star models, Aerts & Rogers (2015) concluded that it might explain the macroturbulence observed in O-stars, based on 2D nonlinear simulations of convection-driven waves in a modified 3 ${M}_{\odot }$ model.

Within our scenario, the Reynolds stresses associated with turbulent pressure induce uncorrelated turbulent pressure fluctuations in the form

$$\begin{eqnarray}&&\delta {P}_{\mathrm{turb}}\sim \rho {v}_{{\rm{c}}}^{2},\end{eqnarray} \tag{ 3 }$$

(Goldreich & Kumar 1990; Grigahcène et al. 2005; Lecoanet & Quataert 2013; Shiode et al. 2013), where $\delta {P}_{\mathrm{turb}}$ is the Lagrangian pressure perturbation associated with the convective motions. Such stochastic fluctuations at the percent level can produce a strong local deviation from hydrostatic equilibrium, and will thus naturally excite high-order pulsations in the range of eigenmodes, which are closest to the spectrum of the fluctuations.

This also would imply that in the luminous stars with log(${\mathcal{L}}/{{\mathcal{L}}}_{\odot })\gt 3.3$ in Figure 4, the macroturbulence may be signifying high-order pulsations, which are, however, excited by turbulent pressure fluctuations rather than through the κ-mechanism or by strange mode instability. If so, on one hand, we can expect that linear pulsation analyses that include the Reynolds stress tensor, as in, e.g., Dupret et al. (2004, 2005) or Antoci et al. (2014), may uncover that stars in a large fraction of the red and orange colored region in Figure 4 are unstable to high-order g-mode pulsations. On the other hand, as the pressure fluctuations in these stars, as predicted by our simple analysis, can amount up to 5% of the total pressure, it is conceivable that in linear stability analyses, which require the growth of the instability from infinitesimal perturbations, an instability is not detected in all stellar models where high-order g-mode can be excited through finite amplitude pressure perturbations.

We implemented the effect of the turbulent motion of convective eddies in a simple formalism in the momentum and energy equations of our stellar evolution code. By comparing with previous computations (Brott et al. 2011), we find that the turbulent pressure does not alter the stellar structure significantly. However, we find maximum turbulent pressure contribution of up to 5% and 30% in our models for OB supergiants and cool red supergiants, respectively. By comparing the maximum turbulent pressure contribution in our models with spectroscopically derived macroturbulent velocities for a large sample of Galactic OB stars (Simón-Díaz 2015), we find both quantities to be strongly correlated.

Several less luminous stars, in which the turbulent pressure is thought to be small, nevertheless show high macroturbulent velocities. These are located close to or inside the region where linear pulsation analysis predicts high-order g-mode pulsations, arguing therefore for κ-mechanism pulsations, and not turbulent pressure fluctuations, as the origin of the macroturbulence phenomenon, which is in line with previous suggestions (Aerts et al. 2009; Simón-Díaz 2015).

We argue that the turbulent pressure fluctuations in hot luminous stars can excite such high-order pulsations, most likely non-radial g-modes, which may explain the occurrence of macroturbulence in stars that are found outside of the currently predicted pulsational instability domains. This view is in agreement with the indication that macroturbulence can be suppressed in strongly magnetic stars, given that such a field might effectively inhibit convective motions in the FeCZ (Sundqvist et al. 2013).

At the moment, turbulent pressure fluctuations appear to be the only mechanism that may excite high-order oscillations in luminous stars (log(L/${L}_{\odot }$ $)\gt 4.5$) in the wide effective temperature regime for which strong macroturbulence is observed.

L.G. is member of the International-Max-Planck-Research-School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne. L.F. acknowledges financial support from the Alexander-von-Humboldt foundation. S.S.-D. acknowledges funding by the Spanish Ministry of Economy and Competitiveness (MINECO) under the grants AYA2010-21697-C05-04, AYA2012-39364-C02-01, and Severo Ochoa SEV-2011-0187. Based on observations made with the Nordic Optical Telescope, operated by NOTSA, and the Mercator Telescope, operated by the Flemish Community, both at the Observatorio del Roque de los Muchachos (La Palma, Spain) of the Instituto de Astrofísica de Canarias. Moreover, L.G. and L.F. thank A. Miglio for useful discussions.