© ESO, 2015

1. Introduction

Massive stars are powerful engines and strongly affect the evolution of star-forming galaxies throughout cosmic time (Bresolin et al. 2008). In particular the most massive stars produce copious amounts of ionising photons (Doran et al. 2013), emit powerful stellar winds (e.g. Kudritzki & Puls 2000; Smith 2014) and, in their final explosions, are suspected to produce the most energetic and spectacular stellar explosions, as pair-instability supernovae (Kozyreva et al. 2014), superluminous supernovae (Gal-Yam et al. 2009), and long-duration gamma-ray bursts (Larsson et al. 2007; Raskin et al. 2008).

Massive main-sequence stars, which we understand here as those that undergo core hydrogen burning, have a much higher luminosity than the Sun, as they are known to obey a simple mass-luminosity relation, L ~ M^α, with α> 1. However, whereas this relation is very steep near the solar mass (α ≃ 5), it is shown in Kippenhahn & Weigert (1990) that α → 1 for M → ∞. Indeed, Köhler et al. (2015) find α ≃ 1.1 for M = 500 M_⊙.

Since the Eddington factor is proportional to L/M, it is debated in the literature whether main-sequence stars of higher and higher initial mass eventually reach the Eddington limit (Langer 1997; Crowther et al. 2010; Maeder et al. 2012). The answer is clearly: yes, they do. Even when only electron scattering is considered as a source of radiative opacity, the Eddington limit corresponds to a luminosity-to-mass ratio of $\mathrm{ℛ:}\mathrm{=}\log \left(\frac{\mathit{L}}{\hspace{0.17em}{\mathit{L}}_{\mathrm{\odot }}}\mathit{/}\frac{\mathit{M}}{\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot }}}\right)\mathrm{\simeq }\mathrm{4.6}$ (Langer & Kudritzki 2014) for hot stars with a solar helium abundance. This is extremely close to the ℛ-values obtained for models of supermassive stars, where this ratio is nearly mass independent (Fuller et al. 1986; Kato 1986). In fact, Kato (1986) showed that zero-age main-sequence (ZAMS) models computed only with electron scattering opacity reach the Eddington limit at a mass of about ~ 10⁵M_⊙.

Whether supermassive stars exist is an open question. Also the mass limit of ordinary stars is presently uncertain (Schneider et al. 2014), however, there is ample evidence for stars with initial masses well above 100 M_⊙ in the local Universe. A number of close binary stars have been found with component initial masses above 100 M_⊙ (Schnurr et al. 2008, 2009; Taylor et al. 2011; Sana et al. 2013). Crowther et al. (2010) proposed initial masses of up to 300 M_⊙ for several stars in the LMC, based on their luminosities. Bestenlehner et al. (2014) identified more than a dozen stars more massive than 100 M_⊙ from the sample of ~ 1000 OB stars near 30 Doradus, which are analysed in the frame of VLT-Flames Tarantula Survey (Evans et al. 2011). The hydrogen-rich stars among them have measured ℛ-values of up to 4.3. In hot stars with finite metallicity, the ion opacities can easily exceed the electron scattering opacity (Iglesias & Rogers 1996). It is thus to be expected that the true Eddington limit, which accounts for all opacity sources, is located at ℛ-values of 4.3 or below. This implies that these stars should in fact also have reached, or exceeded their Eddington limit.

In this paper, we explore the question of massive main-sequence stars reaching, or exceeding the Eddington limit from the theoretical side. We show by means of detailed stellar models, as described in Sect. 2, that even all stars more massive than ~ 40 M_⊙ are found to reach the Eddington limit. In Sect. 3, we demonstrate the need to properly define a local Eddington factor in the stellar interior, which we then use in Sect. 4 to show that when it exceeds the critical value of one, the stellar envelope becomes inflated. We show further in Sects. 5 and 6 that super-Eddington conditions can lead to density inversions and induce convection. We compare our results to previous studies in Sect. 7, and relate them to observations in Sect. 8, before summarising our conclusions in Sect. 9.

2. Stellar models

The grids of the stellar models we used have been published in Brott et al. (2011) and Köhler et al. (2015). We consider only the core hydrogen burning models computed with LMC metallicity. Each stellar evolution sequence computed by Brott et al. (2011) and Köhler et al. (2015) consists of typically 2000 individual stellar models. However, the full amount of data defining a stellar model is only stored for a few dozen time points per sequence, in non-regular intervals. We analyse those stored models. This scheme has the disadvantage that the density of models in the investigated parameter space is not always as high as it should be ideally. Still, as shown below, it allows for a thorough sampling of the considered parameter space, and it is fully consistent with the results already published.

The stellar models were computed with a state-of-the-art, one-dimensional hydrodynamic implicit Lagrangian code (BEC), which incorporates latest input physics (for details, see Braun 1997; Yoon et al. 2006; Brott et al. 2011; Köhler et al. 2015, and references therein). Convection was treated as in the standard non-adiabatic mixing length approach (Böhm-Vitense 1958; Kippenhahn & Weigert 1990) and a mixing length parameter of α = l/H_p = 1.5 (Langer 1991) was adopted, with l and H_p being the mixing length and the pressure scale height, respectively. This value of the mixing length parameter leads to a good representation of the Sun (Suijs et al. 2008), whereas its calibration to multi-dimensional hydrodynamic models shows that it tends to decrease towards lower gravities (Trampedach et al. 2014; Magic et al. 2015). The convective velocities were limited to the local value of the adiabatic sound speed. The contribution of turbulent pressure (de Jager 1984) was neglected, since it is not expected to be important in determining the stellar hydrostatic structure (Stothers 2003). Indeed, our recent study which includes turbulent pressure (Grassitelli et al. in preparation) shows that e.g. for an 80 M_⊙ evolutionary sequence, the stellar radius is increased over that of models without turbulent pressure by at most a few per cent at any time during its main-sequence evolution. Rotational mixing of chemical elements, following Heger et al. (2000), and transport of angular momentum by magnetic fields due to the Spruit-Taylor dynamo were also included (Spruit 2002). The efficiency parameters f_c and f_μ for rotational mixing were set to 0.0228 and 0.1, respectively (Brott et al. 2011). Radiative opacities were interpolated from the OPAL tables (Iglesias & Rogers 1996). The opacity enhancement due to Fe-group elements at T ~ 200 kK plays a vital role in determining the envelope structure in our stellar models. Even though flux-mean opacities are appropriate to study the momentum balance near the stellar photosphere, in the following we only consider the Rosseland mean opacities, which are thought to behave very similarly to the flux-mean opacities especially at an optical depth larger than one.

The outer boundary condition of the stellar models corresponds to a plane-parallel grey atmosphere model on top of the photosphere. In other words, the effective temperature was used as a boundary condition at a Rosseland optical depth of 2/3. The adopted stellar wind mass-loss recipe leads to small, but finite outflow velocities in the outermost layers, which induces a slight deviation from hydrostatic equilibrium.

The mass-loss prescription from Vink et al. (2000, 2001) was employed to account for the winds of O- and B-type stars. Moreover, parameterized mass-loss rates from Nieuwenhuijzen & de Jager (1990) were used on the cooler side of the bi-stability jump, i.e. at effective temperatures less than 22 000 K, if the Nieuwenhuijzen & de Jager (1990) mass-loss rate exceeded that of Vink et al. (2000, 2001). Wolf-Rayet (WR) type mass-loss was accounted for using the empirical prescription from Hamann et al. (1995) divided by a factor of 10 (Yoon et al. 2006), when the surface helium mass fraction became greater than 70%.

Evolutionary sequences of massive stars, with and without rotation, were computed up to an initial mass of 500 M_⊙, starting with LMC composition. The initial mass fractions of hydrogen, helium and metals were taken to be 0.7391, 0.2562, and 0.0047 respectively, in accordance with the observations of young massive stars in the LMC (Brott et al. 2011).

3. The Eddington limit

The Eddington limit refers to the condition in which the outwards radiative acceleration in a star balances the inwards gravity, in hydrostatic equilibrium. It is a concept that is thought to apply at the stellar surface, in the sense that if the Eddington limit is exceeded, a mass outflow should arise (Eddington 1926; Owocki et al. 2004). If we denote the gravity as g = GM/r² and the radiative acceleration mediated through the electron scattering opacity as g_rad = κ_eL/ 4πr², then the classical Eddington factor Γ_e is defined as ${\mathrm{\Gamma }}_{\mathrm{e}}\mathrm{:}\mathrm{=}\frac{\mathit{L}}{{\mathit{L}}_{\mathrm{Edd}}}\mathrm{=}\frac{{\mathit{g}}_{\mathrm{rad}}}{\mathit{g}}\mathrm{=}\frac{{\mathit{\kappa }}_{\mathrm{e}}\mathit{L}}{\mathrm{4}\mathit{\pi cGM}}\mathit{,}$(1)where L,M, and κ_e are the luminosity, mass, and electron-scattering opacity, respectively, with the physical constants having their usual meaning. The classical Eddington parameter Γ_e therefore does not depend on the radius r as the inverse r² scaling in both g_rad and g cancel out. Whereas Γ_e is often convenient to consider, it provides a sufficient instability criterion to stars, but not a necessary one because usually the true opacity exceeds the electron scattering opacity significantly and also contributes to the radiative force.

As it turns out below, even when the Rosseland mean opacities are used, the models we analyse practically never reach the Eddington limit at their surface. Therefore, we instead consider the Eddington factor in the stellar interior as ${\mathrm{\Gamma }}^{\mathrm{\prime }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{:}\mathrm{=}\frac{\mathit{L}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{L}}_{\mathrm{Edd}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathrm{=}\frac{\mathit{\kappa }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{L}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{4}\mathit{\pi cGM}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathit{,}$(2)where M(r) is the Lagrangian mass coordinate, κ(r) is the Rosseland mean opacity, and L(r) is the local luminosity (Langer 1997). However, Γ′(r) > 1 also does not provide a stability limit in the stellar interior because the stellar layers turn convectively unstable following Schwarzschild’s criterion when Γ′(r) → 1 (Joss et al. 1973; Langer 1997). As the luminosity transported by convection does not contribute to the radiative force, we subtract the convective luminosity in the above expression and redefine the Eddington factor as $\mathrm{\Gamma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{:}\mathrm{=}\frac{{\mathit{L}}_{\mathrm{rad}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{L}}_{\mathrm{Edd}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathrm{=}\frac{\mathit{\kappa }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{\left(}\mathit{L}\mathrm{\right(}\mathit{r}\mathrm{)}\mathrm{-}{\mathit{L}}_{\mathrm{conv}}\mathrm{(}\mathit{r}\mathrm{\left)}\mathrm{\right)}}{\mathrm{4}\mathit{\pi cGM}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathrm{·}$(3)For example, near the stellar core where convective energy transport is highly efficient, Γ(r) stays well below unity in spite of Γ′(r) ≫ 1 and no instability, i.e. departure from hydrostatic equilibrium, occurs (see Fig.  1). In the rest of the paper we refer to Γ(r) as the Eddington factor unless explicitly specified otherwise.

Even with this definition, a super-Eddington layer inside a star does not necessarily lead to a departure from hydrostatic equilibrium or a sustained mass outflow. In the outer envelopes of massive stars, non-adiabatic conditions prevail and convective energy transport is highly inefficient, which pushes Γ(r) close to (or above) one. We find that the stellar models counteract this kind of a super-Eddington luminosity by developing a positive gas pressure gradient, thus restoring hydrostatic equilibrium (Langer 1997; Asplund 1998). In these situations, the canonical definition of L_edd being the maximum sustainable radiative luminosity locally in the stellar interior (in hydrostatic equilibrium) breaks down and loses its significance. As we see below, the radiative luminosity beneath the photosphere can be up to a few times the Eddington luminosity.

In Fig.  1, the behaviour of Γ and Γ′ is shown along with the electron-scattering Eddington factor Γ_e in a 285 M_⊙ non-rotating stellar model, which provides an educative example (see Appendix D for further examples). As explained above, Γ′ and Γ_e are significantly greater than one in the convective core of the star. The indicated sub-surface convection zones are caused by the opacity peaks at T ~ 1.5 × 10⁶ K (deep iron bump) and at T ~ 2 × 10⁵ K (iron bump). Near the bottom of the inflated envelope (r/r_core ≳ 1; see Sect. 4 for the definition of r_core), Γ approaches one and the Fe opacity bump drives convection. An extended region with Γ ≈ 1 follows. A thin shell very close to the photosphere contains the layers with a positive density gradient and with Γ > 1.

The stellar models have been computed with a hydrodynamic stellar evolution code, however, because of the large time steps required for stellar evolution calculations, non-hydrostatic solutions are suppressed by our numerical scheme. The resulting hydrostatic structures are still valid solutions of the hydrodynamic equations (see Heger et al. 2000; Kozyreva et al. 2014, for the equations to be solved). Models computed with time steps small enough to resolve the hydrodynamic timescale reveal that some, and potentially many of our models are pulsationally unstable, as will be shown in a forthcoming paper. However, in the cases analysed so far, the pulsations saturate and do not lead to a destruction or ejection of the inflated envelopes. In this respect, we consider our analysis of the hydrostatic equilibrium structures as useful.

3.1. Effect of rotation on the Eddington limit

The effect of the centrifugal force on the structure of rotating stellar models has been studied by a number of groups in the past, including Heger et al. (2000) and Maeder & Meynet (2000). This is done by describing the models in a 1D approximation where all variables are taken as averages over isobaric surfaces (Kippenhahn & Thomas 1970). The stellar structure equations are modified to include the effect of the centrifugal force (Endal & Sofia 1976). The equation of hydrostatic equilibrium becomes $\frac{\mathrm{d}\mathit{P}}{\mathrm{d}\mathit{m}}\mathrm{4}\mathit{\pi }{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{f}}_{\mathit{P}}\frac{\mathit{GM}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{=}\mathrm{0}\mathit{,}$(4)and the radiative temperature gradient in the energy transport equation (in the absence of convection) takes the form ${\mathrm{\nabla }}_{\mathrm{rad}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{16}\mathit{\pi acG}}\frac{\mathit{\kappa PL}}{\mathit{M}{\mathit{T}}^{\mathrm{4}}}\frac{{\mathit{f}}_{\mathit{T}}}{{\mathit{f}}_{\mathit{P}}}\mathit{,}$(5)where the quantities f_P and f_T have the same definition as in Heger et al. (2000). Consequently, the Eddington luminosity gets modified as: ${\mathit{L}}_{\mathrm{Edd}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi cGM}}{\mathit{\kappa }}\frac{{\mathit{f}}_{\mathit{P}}}{{\mathit{f}}_{\mathit{T}}}\mathrm{·}$(6)However, the Eddington factor, $\mathrm{\Gamma }\mathrm{=}\frac{{\mathit{L}}_{\mathrm{rad}}}{{\mathit{L}}_{\mathrm{Edd}}}\mathrm{=}\frac{\mathrm{\nabla }}{{\mathrm{\nabla }}_{\mathrm{rad}}}\frac{\mathit{L}}{{\mathit{L}}_{\mathrm{Edd}}}\mathrm{=}\frac{\mathrm{4}\mathit{a}}{\mathrm{3}}\frac{{\mathit{T}}^{\mathrm{4}}\mathrm{\nabla }}{\mathit{P}}$(7)does not have any explicit dependence on f_P and f_T because the factor f_P/f_T cancels out. Therefore formally, the Eddington factor remains unaffected by rotation. Of course, if the internal evolution of a rotating model is changed, for example, by rotational mixing, its Eddington factor is still different from that of the corresponding non-rotating model.

Of course, real stars are three-dimensional and the centrifugal force must affect the hydrostatic stability limit. However, this is expected to be a function of the latitude at the stellar surface, and in a 2D view, the effect is largest at the equator (Langer 1997). To first order, the critical luminosity L_c to unbind matter at the stellar equatorial surface becomes ${\mathit{L}}_{\mathrm{c}}\mathrm{=}{\mathit{L}}_{\mathrm{Edd}}\left(\mathrm{1}\mathrm{-}{\left(\frac{{\mathit{v}}_{\mathrm{rot}}}{{\mathit{v}}_{\mathrm{Kep}}}\right)}^{\mathrm{2}}\right)\mathit{,}$(8)where v_rot and v_Kep are the stellar equatorial rotation velocity and the corresponding Keplerian value, respectively. However, to compute the effect reliably, the stellar deformation due to rotation as well as the effect of gravity darkening need to be accounted for (Maeder & Meynet 2000; Maeder 2009). To do this realistically for stars near the Eddington limit requires 2D calculations at least.

The implication is that the effect of rotation on the critical stellar luminosity cannot be properly described through the models analysed here. Those models see the same critical luminosity as if rotation was absent. Since mixing of helium in these models is very weak for rotation rates below those required for chemically homogeneous evolution, most of the rotating models evolve very similar to the non-rotating models (Brott et al. 2011; Köhler et al. 2015), and thus merely serve to augment our database.

3.2. The maximum Eddington factor

In our stellar models, we have determined the maximum Eddington factor Γ_max over the whole star, i.e Γ_max: = max [Γ(r)]. The maximum Eddington factor Γ_max generally occurs in the outer envelopes of our models, where convective energy transport is much less efficient than in the deep interior. The variation of Γ_max across the upper HR diagram is shown in Fig. 2 for all analysed core hydrogen burning models that have Γ_max> 0.9.

Three distinct regions with Γ_max> 1 can be identified in Fig. 2, which can be connected to the opacity peaks of iron, helium, and hydrogen. When one of these opacity peaks is situated sufficiently close to the stellar photosphere, the densities in these layers are so small that convective energy transport becomes inefficient. As a consequence, super-Eddington layers develop that are stabilised by a positive (i.e. inwards directed) gradient in density and gas pressure (see Sect. 5 below). The envelope inflation, which occurs when Γ_max approaches one, is discussed in Sect. 4.

Figure 3 shows Γ_max as a function of the effective temperature for all our models with Γ_max> 0.9. The models that have the hydrogen opacity bump close to their photosphere can obtain values of Γ_max as high as ~ 7. This manifests itself as a prominent peak around T_eff ≈ 5.5 kK. The inset panel shows the much weaker peaks in Γ_max due to partial ionisation zones of Fe and HeII, at T/ kK ~ 200 and 50, respectively. The peak caused by the Fe opacity bump may extend to hotter effective temperatures and apply to hot, hydrogen-free Wolf-Rayet stars, which are not part of our model grid.

For stars above about 125 M_⊙, Γ_max reaches one, even on the ZAMS. This is demonstrated in Fig. 4 which shows both Γ_max and effective temperature as a function of mass for the non-rotating stellar models. As these models evolve away from the ZAMS to cooler temperatures, super-Eddington layers develop in their interior. The blue curve shows a maximum effective temperature of 57 000 K at about 200 M_⊙, beyond which it starts decreasing with further increase in mass. This behaviour is related to the phenomenon of inflation, which is discussed in detail in Sect. 4. However, the “effective temperature” at the base of the extended, inflated envelope, T_eff,core (see Sect. 4), still increases with mass in the whole considered mass range.

3.3. The spectroscopic HR diagram

Figure 5 shows the location of our analysed models in the spectroscopic HR diagram (sHRD, Langer & Kudritzki 2014; Köhler et al. 2015) where instead of the luminosity, $L\mathrm{:}\mathrm{=}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathit{/}\mathit{g}$ is plotted as a function of the effective temperature. The quantity L can be measured for stars without knowing their distance. Moreover, we have log (L/ L_⊙) = ℛ (cf. Sect. 1), such that L is directly proportional to the Eddington factor Γ_e as ${\mathrm{\Gamma }}_{\mathrm{e}}\mathrm{=}\frac{{\mathit{\kappa }}_{\mathrm{e}}\mathit{L}}{\mathrm{4}\mathit{\pi cGM}}\mathrm{=}\frac{{\mathit{\kappa }}_{\mathrm{e}}\mathit{\sigma }{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}}{\mathit{cg}}\mathrm{=}\frac{{\mathit{\kappa }}_{\mathrm{e}}\mathit{\sigma }}{\mathit{c}}L\mathit{,}$(9)where g is the surface gravity and the constants have their usual meaning. Therefore one can directly read off Γ_e (right Y-axis in Fig. 5) from the sHRD. Massive stellar models often evolve with a slowly increasing luminosity over their main-sequence lifetimes. Therefore, where in a conventional HR diagram models with very different Γ_max might cluster together (see Fig. 2), they separate out nicely in the sHRD since L∝ L/M. The effect of the opacity peaks on the maximum Eddington factor (Γ_max) at temperatures corresponding to the three partial ionisation zones (Fe, HeII and H) is seen more clearly in the sHRD in Fig. 5 compared to the ordinary HR diagram (Fig. 2).

We find that for our ZAMS models the electron scattering opacity is κ_e ≈ 0.34, while the true photospheric opacity κ_ph is around 0.5. Therefore it is expected that the true Eddington limit (Γ = 1) is achieved at about Γ_e = 0.7 for stellar models that retain the initial hydrogen abundance at the photosphere. Therefore in Fig. 5 we have drawn two horizontal lines corresponding to Γ_e = 0.7, one assuming the initial hydrogen mass fraction X = 0.74 (green line) and the other assuming X = 0 (red line). While models with helium-enriched photospheres exceed the green line comfortably, even the most helium-enriched models (rotating or otherwise) stay below the red line.

From Köhler et al. (2015), we know that models with log L/ L_⊙> 4.4 are all hydrogen-deficient, either due to mass-loss or rotationally induced mixing, as both processes lead to an increasing L/M-ratio (cf. their Fig. 18). Figure 5 thus demonstrates that the models that contain super-Eddington layers due to the partial ionisation of helium all have hydrogen deficient envelopes, i.e. they are correspondingly helium-enriched.

Figure 5 reveals that the electron scattering Eddington factor Γ_e is not a good proxy for the maximum true Eddington factor (Γ_max) obtained inside the star. For example, along the horizontal line log L/ L_⊙ = 4.3, corresponding to Γ_e ≃ 0.5, Γ_max varies from well below one to values near seven at the cool end. However, we note that below 30 000 K helium and hydrogen recombine, the gas is not fully ionised any more. The line opacities of helium and hydrogen become important, which causes the increase in Γ_max (see Figs. 3 and 14).

3.4. Surface Eddington factors and the location of Γ_max

The optical depth where Γ_max is reached gives an idea of how deep in the stellar interior the layer with the highest Eddington factor is located. We investigate this in Fig. 6, which shows that Γ_max is located at largely different optical depths in different types of models. While the maximum Eddington factors occur generally at optical depths below ~ 10 000, we see that in the three effective temperature regimes identified by the super-Eddington peaks in Figs. 2 and 3, Γ_max can even be located at an optical depth of ~ 10 or below.

For example, when the tracks above log L/L_⊙ = 6.2 approach effective temperatures of ~ 30 kK, Γ_max is located at the Fe-peak, which is deep inside the envelope (τ ≃ 1400). The models at this stage become helium-rich (Y_s ≳ 70%) and the Wolf-Rayet mass-loss prescription is applied. Once these tracks turn bluewards in the HR diagram, the position of Γ_max jumps to the helium opacity peak, which is located much closer to the stellar surface. Consequently, we find three orders of magnitude of difference between these two types of models with similar effective temperature and luminosity.

When considering the surface Eddington factors in the sHRD (Fig. 7), we see that only the models with Γ(R_⋆) > 0.98 have log L/ L_⊙ values of more than 4.6. As discussed above, these models, which started on the main-sequence with initial masses above 300 M_⊙ are extremely helium-rich and may correspond to the most extreme late-type WN stars (Sander et al. 2014). As shown in Fig. 8, they exceed the Eddington limit by just a few per mille, which is possible because of the high assumed mass-loss rates that imply a slight deviation from hydrostatic equilibrium near the stellar surface (cf. Sect. 3 above). However, these models are those in which our assumption of an optically thin wind might break down (see Fig. 7 in Köhler et al. 2015). Since the inclusion of an optically thick outflow may lead to changes of the temperature and density structure near the surface, the surface Eddington-factors for these particular models are not reliable.

In summary, we find, on the one hand, that many of our models contain layers at optical depths between a few and a few thousand in which the Eddington factor exceeds the critical value of one. On the other hand, for none of our models we can conclude that the Eddington limit is reached very near to, or at the surface, where for the vast majority we can even exclude that this happens. This finding leads to a shift in the expectation of the response of stars that reach the Eddington limit during their evolution. We might not expect direct outflows driven by super-Eddington luminosities, but instead internal structural changes, in particular envelope inflation.

4. Envelope inflation

Inflation of massive, luminous stars refers to the formation of extended, extremely dilute stellar envelopes. An example of an inflated model is shown in Fig. 9. The red shaded region is the non-inflated core and the blue shaded region is what we refer to as the inflated envelope. In the example, the model is inflated by 60% of its core radius (defined below). In the presented model, the inflated envelope only contain a small fraction of a solar mass, i.e. ≈ 10^-5M_⊙.

Envelope inflation is inherently different from classical red supergiant formation. The latter occurs after core hydrogen exhaustion, as a consequence of vigorous hydrogen shell burning. This process expands all layers above the shell source, which usually comprise several solar masses in massive stars, and it also operates in low-mass stars, such that no proximity to the Eddington limit is required. The mechanism of envelope inflation that we discuss already works during core hydrogen burning, i.e. even on the ZAMS for sufficiently luminous stars (cf. Fig. 4). Previous investigations have suggested that inflation is related to the proximity of the stellar luminosity to the Eddington luminosity (Ishii et al. 1999; Petrovic et al. 2006; Gräfener et al. 2012) in the envelopes of massive stars with a high luminosity-to-mass ratio (≳ 10⁴L_⊙/M_⊙). The amount of mass contained in an inflated envelope is usually very small. As we shall see below, inflation, in extreme cases, can also produce core hydrogen burning red supergiants.

We define inflation in our models through Δr/r_core: = (R_⋆−r_core) /r_core, with r_core being the radius at which inflation starts and R_⋆, the photospheric radius. Since the densities in inflated envelopes are small, the dominance of radiation pressure in these envelopes is much larger than it is in the main stellar body. We define a model to be inflated if β(r), which is the ratio of gas pressure to total pressure, reaches a value below 0.15 in the interior of a model. The radius at which β goes below 0.15 for the first time from the centre outwards is denoted as r_core, i.e. the start of the inflated region. The remaining extent of the star until the photosphere (R_⋆−r_core) is considered the inflated envelope.

We emphasize that our choice of the threshold value for β is arbitrary and not derivable from first principles. However, we have verified that this prescription identifies inflated stars in different parts of the HR diagram very well (cf. Appendix D). As β → 0 for M → ∞, our criterion may fail for extreme masses, however, the mass averaged value of β for the most massive ZAMS model that we analysed (500 M_⊙) is 0.3. A threshold value of 0.15 thus appears adequate for the present study. As an example, let us consider a typical inflated model, shown in Appendix A. The value of β in Fig. A.4 decreases sharply at the base of the inflated envelope, to around 0.01. Even if the β threshold is varied by 30%, i.e. 0.15 ± 0.045, the non-inflated core radius r_core changes by only 4%. This goes to show that for clearly inflated models, the value of r_core is insensitive to the threshold value of β.

We furthermore performed a numerical experiment which is suited to show that the core radii identified as described above are indeed robust. We chose an inflated 300 M_⊙ model, and then increased the mixing length parameter α such that convection becomes more and more efficient. As shown in Fig. B.1, as a result the extent of the inflated envelope decreased without affecting the model structure inside the core radius, which thus remained independent of α. For α = 40 convection became nearly adiabatic, inflation almost disappeared, and the fact that the photospheric radius in this case became very close to the core radius validates our method of identifying r_core.

Figure 10 shows the amount of inflation as defined above, for all our models that fulfil the inflation criterion, in the HR diagram. It reveals that overall, inflation is larger for cooler temperatures. This is not surprising since inflation does not appear to change the stellar luminosity and must therefore induce smaller surface temperatures. We also see that inflation is larger for more luminous stars, which is expected because the Eddington limit is supposed to play a role (see below). We also find inflation increases along the evolutionary tracks of the most massive stars that turn back from the blue supergiant stage, which in this case is due to the shrinking of their core radii. A distinction between the inflated and the non-inflated models is made by drawing the black lines in Fig. 10. They are drawn such that no model is inflated below the solid line and all the models are inflated above the dotted line. In between these two lines we find a mixture of both inflated and non-inflated models. We find that essentially all models above log (L/L_⊙) ≃ 5.6 are inflated. Consequently, stars above ~ 40 M_⊙ inflate during their main-sequence evolution.

Figure 11 shows the inflation factor as function of the stellar effective temperature for our inflated models. Whereas inflation increases the radius of our hot stars by up to a factor of 5, the cool supergiant models can be inflated by a factor of up to 40. We refer to Appendix A for the detailed structures of several inflated models.

In Fig. 12, we take a look at inflation as a function of the core effective temperature T_eff,core defined as ${\mathit{T}}_{\mathrm{eff}\mathit{,}\mathrm{core}}\mathrm{=}\frac{\mathit{L}}{\mathrm{4}\mathit{\pi \sigma }{\mathit{r}}_{\mathrm{core}}^{\mathrm{2}}}\mathit{,}$(10)where L refers to the surface luminosity and the constants have their usual meaning. We can see that even our coolest models have high core effective temperatures in the sense that if their inflated envelopes were absent their stellar effective temperatures would have been higher than 20 000 K. Those stars with stellar effective temperatures below ~50 000 K contain the He II ionisation zone within their envelopes, and stars with stellar effective temperatures below ~10 000 K also contain the H/He I ionisation zone. However, as revealed by the density and temperature structure of these models (cf. Appendix D), the temperature at the bottom of the inflated envelope is always about 170 000 K, and thus corresponds to the temperature of the iron opacity peak. We conclude that the iron opacity is at least in part driving the inflation of all the stars. For those with cool enough envelopes, helium and hydrogen are likely relevant in addition.

4.1. Why do stellar envelopes inflate?

As suggested earlier, the physical cause of inflation in a given star may be its proximity to the Eddington limit. Figure 13 shows the correlation between inflation and Γ_max for our models. As expected, we find that our stellar models are not inflated when Γ_max is significantly below 1, and they are all inflated for Γ_max> 1. Indeed, the top panel of Fig. 13 gives the clear message that the Eddington limit, in the way it is defined in Sect. 3, is likely connected with envelope inflation.

Comparing Fig. 13 (top panel) to Fig. 11 shows that inflation increases up to T_eff ≈ 5500 K. Thereafter, T_eff and Γ_max decrease and the stars keep getting bigger without significant changes in r_core, and, hence, inflation still increases. However, the drop in inflation for the coolest models shows an opposite trend. This is because the non-inflated core radius r_core now moves outwards (increases) such that inflation (Δr/r_core) decreases even though R_⋆ keeps increasing (cf. definition of r_core in Sect. 4).

In the zoom-in view shown in the lower panel of Fig. 13, we see some models being inflated for Γ_max in the range ~ 0.9...1. Partly, this may be because of the arbitrariness in our definition of inflation. The exact value of Δr/r_core depends somewhat on the choice of the threshold value of β to characterize inflation (cf. Appendix D), i.e. the models with Δr/r_core ≲ 2 and Γ_max< 1 may be at the borderline of inflation. The models with Δr/r_core ≲ 2 but Γ_max> 1 are all very hot (T_eff ≳ 40 000 K) and in those models, the inflation is intrinsically small, but generally unambiguous.

Still, we see a significant number of models below the Eddington limit (Γ_max< 1), which show a quite prominent inflation, i.e. which have a radius increase due to inflation of more than a factor of five. We investigated this kind of a model by artificially increasing its mass-loss rate above the critical value Ṁ_crit (Petrovic et al. 2006), such that the inflated envelope was removed (cf. Sect. 4.2). We then found that, on turning down the mass-loss rate to its original value, the model regained its initial inflated structure with Γ_max< 1. However, Γ_max = 1 was reached and exceeded in the course of our experiment. We conclude that a stellar envelope may remain inflated even if the condition Γ_max = 1 is not met any more in the course of evolution, but that Γ_max ≳ 1 may be required to obtain inflation in the first place.

In contrast to earlier ideas of a hydrodynamic outflow being triggered when the stellar surface reaches the Eddington limit (Eddington 1926; Owocki et al. 2004), in our models this never happens. When the properly defined Eddington limit is reached inside the envelope, its outermost layers expand hydrostatically and produce inflation. Two possibilities arise in this process. When the star approaches the Eddington limit, the ensuing envelope expansion leads to changes in the temperature and the density structure. Consequently, the envelope opacity can either increase or decrease. Figure 14 shows that the effect of expansion generally leads to a reduced opacity such that the expansion is indeed alleviating the problem. The star then expands until the Eddington limit is just not exceeded any more, which is the reason why we find so many inflated models with Γ_max ≃ 1.

Figure 14 shows the OPAL opacities for hydrogen-rich composition for various constant values of R as function of temperature, where R = ρ/ (T/ 10⁶)³. Kippenhahn & Weigert (1990) showed that for constant β = P_gas/P and constant chemical composition, R as a function of spatial co-ordinate inside the star is a constant. Thus, for un-inflated models, the opacity curves in Fig. 14 may closely represent the true run of opacity with temperature inside the star. In the inflated models, β is dropping abruptly at the base of the inflated envelope, which means that the opacity is jumping from a curve with a higher R-value to one with a lower R-value at this location. That is, the opacity is smaller everywhere in the inflated envelope compared to the situation where inflation would not have happened.

For the chemical composition given in Fig. 14 and assuming β ≡ const., we find $\mathit{R}\mathrm{\simeq }\mathrm{1.8}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{-5}\frac{\mathit{\beta }}{\mathrm{1}\mathrm{-}\mathit{\beta }}\mathit{,}$(11)such that if β drops from 0.5 in the bulk of the star to 0.1 in the inflated envelope, R drops by one order of magnitude. The corresponding reduction in opacity can be significant, i.e. up to about a factor of two.

When upon expansion the envelope becomes cool enough for another opacity bump to come into play, the problem of not exceeding the Eddington limit might not be solvable this way. Instead, when a new opacity peak is encountered in the outer part of the envelope, super-Eddington conditions occur, i.e. layers with Γ_max> 1 (cf. Figs. 2 and 5), along with a strong positive gas pressure (and density) gradient (cf. Sects.  3 and 5). This is most extreme when the envelopes become cool enough (T_eff ≲ 8000 K) such that the hydrogen ionisation zone is present in the outer part of the envelope, where Eddington factors of up to seven are achieved.

4.2. Influence of mass-loss on inflation

One might wonder about the sustainability of the inflated layers against mass-loss, which is an important factor in the evolution of metal-rich massive stars. Petrovic et al. (2006) estimated that the inflated envelope cannot be replenished when the mass-loss rate exceeds a critical value of $\mathit{Ṁ}\mathrm{crit}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{r}}_{\mathrm{core}}^{\mathrm{2}}\hspace{0.17em}{\mathit{\rho }}_{\min }\sqrt{\frac{\mathit{GM}}{{\mathit{r}}_{\mathrm{core}}}}\hspace{0.17em}\mathit{,}$(12)where M and r_core stand for the stellar mass and the un-inflated radius, respectively, and ρ_min is the minimum density in the inflated region. Petrovic et al. (2006) found Ṁ_crit ~ 10^-5M_⊙ yr^-1 for a massive hydrogen-free Wolf-Rayet star of 24 M_⊙. However, for a typical inflated massive star on the main-sequence (see Fig. A.1), this critical mass-loss rate is of the order 10^-3...10^-1M_⊙ yr^-1. These high mass-loss rates are only expected in luminous blue variable (LBV)-type giant eruptions. The mass-loss rates applied to our models are several orders of magnitude smaller (cf. Köhler et al. 2015).

The mass-loss history of four evolutionary sequences without rotation are shown in Fig. 15. Even the 500 M_⊙ model never exceeds a mass-loss rate ~ 5 × 10^-4M_⊙ yr^-1. The critical mass-loss rate for all models shown in Fig. 15 is much higher than the actual mass-loss rates applied. Whereas Ṁ_crit typically exceeds Ṁ by a factor of 1000 for the inflated models in the 50 M_⊙ sequence, it exceeds that of the 500 M_⊙ sequence by a factor of 3...100. It is thus not expected that mass-loss prevents the formation of the inflated envelopes in massive stars near the Eddington limit. In fact, it may be difficult to identify a source of momentum that might drive such strong mass-loss (Shaviv 2001; Owocki et al. 2004). Gräfener et al. (2011) in the Milky Way and Bestenlehner et al. (2014) in the LMC found a steep dependence of the mass-loss rates on the electron-scattering Eddington factor Γ_e for very massive stars, but they do not find mass-loss rates that substantially exceed 10^-4M_⊙ yr^-1.

As many of the models analysed here may be pulsationally unstable, the mass-loss rates may be enhanced in this case. Grott et al. (2005) show that hot stars near the Eddington limit may undergo mass-loss due to pulsations, although extreme mass-loss rates are not predicted. For very massive cool stars, on the other hand, Moriya & Langer (2015) find that pulsations may enhance the mass-loss rate to values of the order of 10^-2M_⊙ yr^-1. These extreme values could prevent the corresponding stars to spend a long time on the cool side of the Humphreys-Davidson limit. A detailed consideration of this issue is beyond the scope of the present paper.

5. Density inversions

An inflated envelope can be associated with a “density inversion” near the stellar surface, i.e. a region where the density increases outwards. An example is shown in Fig. 9. In hydrostatic equilibrium, Γ(r) > 1 implies $\frac{\mathrm{d}{\mathit{P}}_{\mathrm{gas}}}{\mathrm{d}\mathit{r}}\mathit{>}\mathrm{0}$, and thus $\frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}\mathit{r}}\mathit{>}\mathrm{0}$. As a consequence, all the models that have layers in their envelopes exceeding the Eddington limit show density inversions. The criterion for density inversion can be expressed as (Joss et al. 1973; Paxton et al. 2013) $\frac{{\mathit{L}}_{\mathrm{rad}}}{{\mathit{L}}_{\mathrm{Edd}}}\mathit{>}{\left[\mathrm{1}\mathrm{+}{\left(\frac{\mathit{\partial }{\mathit{P}}_{\mathrm{gas}}}{\mathit{\partial }{\mathit{P}}_{\mathrm{rad}}}\right)}_{\mathit{\rho }}\right]}^{-1}\mathit{,}$(13)where P_gas, P_rad and ρ stand for the gas pressure, radiation pressure, and density, respectively. A density inversion gives an inwards force and acts as a stabilizing agent for the inflated envelopes. In the above inequality, P_gas is assumed to be a function of ρ and T only, i.e. the mean molecular weight μ is assumed to be constant. Density inversions might also be present in low-mass stars like the Sun where they are caused by the steep increase of μ around the hydrogen recombination zone (cf. Érgma 1971).

Figure 16 identifies our core hydrogen burning models, which contain a density inversion. The quantity Δρ/ρ represents the strength of the density inversion normalized to the minimum density attained in the inflated zone. We can identify three peaks in Δρ/ρ at T_eff/ kK ~ 55, 25, and 5.5 (see also Fig. 17), which coincides exactly with the three T_eff-regimes in which models exceed the Eddington limit (cf. Fig. 3). The maximum of the density inversions in the three zones is related to the relative prominence of the three opacity bumps of Fe, HeII, and H respectively, as shown in Fig. 17.

However, an inflated model is not necessarily accompanied by a density inversion. This is depicted clearly in Fig. 18 where we investigate the correlation between inflation and density inversion (this can also be seen by comparing Figs. 10 to 16). Figure 18 shows many models that are even substantially inflated but do not develop a density inversion. The three peaks in the distribution of density inversions of Fig. 17 also show up distinctly in this plot at the three characteristic effective temperatures (shown in colour). Models that show a density inversion do always show some inflation. This is less obvious from Fig. 18 because the hottest models show the smallest amount of inflation (Fig. 11).

The stability of density inversions in stellar envelopes has been a matter of debate for the last few decades but there has been no consensus on this issue yet (see Maeder 1992). There have been early speculations by Mihalas (1969) while studying red supergiants that a density inversion might lead to Rayleigh-Taylor instabilities (RTI), resulting in “elephant trunk” structures washing out the positive density gradient. However, as rightly pointed out by Schaerer (1996), RTI does not develop since the effective gravity g_eff = g(1−Γ) acting on the fluid elements is directed outwards in the super-Eddington layers, which contain the density inversion. Kutter (1970), on the other hand, claimed that a hydrodynamic treatment of the stellar structure equations prevents any density inversion and instead leads to a steady mass outflow. However, this claim is refuted by the present work, since our code solves the 1D hydrodynamic stellar structure equations, in agreement with previous hydrodynamical models by Glatzel & Kiriakidis (1993) and Meynet (1992). Stothers & Chin (1973) suggested that density inversions will lead to strong turbulent motions instead of drastic mass-loss episodes. However, these layers are unstable to convection, so turbulence is present in any case.

Additionally, Glatzel & Kiriakidis (1993) argued in favour of a sustainability of density inversions in the sense that they can be viewed as a natural consequence of strongly non-adiabatic convection, and pointed out that the only plausible way to suppress density inversions is to use a different theory of convection. The only instability expected from simple arguments therefore is convection, which is in line with Wentzel (1970) and Langer (1997).

Still, Ekström et al. (2012) and Yusof et al. (2013) recently considered density inversions as “unphysical”. Density inversions have been suppressed in their models by replacing the pressure scale height in the Mixing Length Theory with the density scale height (cf. Sect. 7). This approach has often been adopted by many investigators in the past to prevent numerical difficulties. As the density scale height tends to infinity when a density inversion starts to develop, this measure tends to enormously increase the convective flux in the relevant layers. It is doubtful whether in reality the convective flux can be increased so much, as the ratio of the local thermal to the local dynamical timescale in the relevant layers is much smaller than one, such that convective eddies lose their thermal energy much faster than they rise, and thus hardly transport any energy at all. Multi-dimensional hydrodynamical simulations would help to settle this issue. We briefly return to this point in Sect. 7.

6. Sub-surface convection

We also studied the convective velocities in the sub-surface convection zones associated with the opacity peaks in our stellar models (Cantiello et al. 2009). We measure these velocities in units of either the isothermal or the adiabatic sound velocity, i.e. c_s,ad and c_s,iso, respectively, which we compute as ${\mathit{c}}_{\mathrm{s}\mathit{,}\mathrm{ad}}\mathrm{=}\sqrt{\frac{\mathit{\gamma P}}{\mathit{\rho }}}$(14)and${\mathit{c}}_{\mathrm{s}\mathit{,}\mathrm{iso}}\mathrm{=}\sqrt{\frac{{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\mathit{\mu }}}\mathrm{=}\sqrt{\frac{{\mathit{P}}_{\mathrm{gas}}}{\mathit{\rho }}}\mathit{,}$(15)where γ is the adiabatic index, P is the total pressure, ρ is the density, μ is the mean molecular weight, T is the temperature, and k_B is the Boltzmann constant. We define M_iso as the maximum ratio of the convective velocity over the isothermal sound speed in the stellar envelope, and M_ad correspondingly using the adiabatic sound speed.

The true sound speed is in between the adiabatic and isothermal speed, closer to the first one in the inner parts of the star, and closer to the second in the inflated stellar envelope (cf. Sect. 5). In Figs. 19 and 20, we show the values of M_iso and M_ad for our models in the HR diagram. Whereas the convective velocities are always smaller than the adiabatic sound speed, Fig. 19 shows that the isothermal sound speed can be exceeded locally in our models by a factor of a few. The convective velocity and sound speed profiles for an extreme model are presented in Appendix C.

Supersonic convective velocities (adiabatic or isothermal, depending on the physical conditions in the envelope) may not be realistic and are outside the frame of the standard Mixing Length Theory. Therefore, in some of our models, the convective velocities, and thus the convective energy transport, may have been overestimated. A limitation of the velocities to the adequate sound speed is expected to reduce the convective flux, which might lead to further inflation of the stellar envelope.

The cool models, with the strongest inflation, have relatively smaller values of M_iso (compared to the hot WR-type models) but large values of M_ad (Fig. 20) in the sub-surface convection zones. This is primarily because of the fact that while c_s,ad depends on the total pressure, c_s,iso only depends on the gas pressure. In the very outer layers of the cool, luminous models, β → 1 and hence P_gas ≈ P_tot. In these situations, c_s,ad and c_s,iso are only a factor $\sqrt{\mathit{\gamma }}$ apart, where γ is the adiabatic index.

We find that the convective energy transport is not always negligible in the inflated models (cf. Sect. 4.1). We therefore evaluate the amount of flux that is actually carried by convection in the inflated envelopes of our models. We define the quantity η(M_iso) as the fraction of the total flux carried by convection in the stellar envelope, at the location where the isothermal Mach number is the largest. This quantity is plotted as a function of the effective temperature in Fig. 21. It is evident from this figure that η(M_iso) need not be small for stellar envelopes to be inflated. However, the hotter a model is the lower its η(M_iso)-value at a given luminosity (see Fig. 22). For models hotter than T_eff ≈ 63 kK (e.g. the hydrogen-free He stars), η(M_iso) indeed goes towards zero (Grassitelli et al., in prep.). The behaviour of the quantity η(M_iso) in the HR diagram is shown in Fig. 22.

7. Comparison with previous studies

7.1. Stellar atmosphere and wind models

Since the Eddington limit was thought to be reached in massive stars near their surface (cf. Sect. 3), several papers have investigated this using stellar atmosphere calculations. Lamers & Fitzpatrick (1988) investigated the Eddington factors in the atmospheres of luminous stars in the temperature-gravity diagram, while Ulmer & Fitzpatrick (1998) did so in the HR diagram. Both studies took the full radiative opacity into account. While their technique did not allow them to reach Eddington factors of one or more, Ulmer and Fitzpatrick found that model atmospheres with a maximum Eddington factor of 0.9 are located near the observed upper luminosity limit of stars in the Large Magellanic Cloud (Humphreys & Davidson 1979; Fitzpatrick & Garmany 1990).

The main feature in the lines of constant Eddington factors in the HR diagram found by Ulmer and Fitzpatrick is a drop from T_eff ≃ 60 000 K to 15 000 K. This may correspond to the drop in the maximum Eddington factor seen in our models in the same temperature interval (cf. Figs. 2 and 5). Note that the peak around T_eff ≃ 30 000 K in Fig. 3 only corresponds to helium-rich models, which are not considered by Ulmer and Fitzpatrick.

On the other hand, neither inflation nor super-Eddington layers or density inversions are reported by Ulmer and Fitzpatrick, or, to the best of our knowledge, from any hot, main-sequence star model atmosphere calculation so far. One reason might be that many model atmospheres only include a rather limited optical depth range (e.g. up to τ = 100 in Ulmer and Fitzpatrick), such that the iron opacity peak is often not included in the model. Additionally, the computational methods employed might not allow for a non-monotonic density profile.

Given the ubiquity of inflation for models above log  L/L_⊙> 5.5 or M> 50 M_⊙ in the LMC, and a correspondingly lower limit in the Milky Way due to its higher iron content, it is helpful to construct model atmospheres that include this effect and identify its observational signatures. As the density profiles of these kinds of atmospheres near the photosphere are significantly different from those in non-inflated atmospheres, these signatures may indeed be expected.

Asplund (1998) gives a thorough analysis of the Eddington limit in cool star atmosphere models. Indeed, he finds super-Eddington layers and density inversions in his models, and gives arguments for the physically appropriate nature of these phenomena. He also discusses the effects of stellar winds on these features, and finds they may be suppressed by extremely strong winds, but not by winds with mass-loss rates in the observed range. Asplund does not find inflation in his models, again, arguably because the iron opacity peak is not included in his model atmospheres, which appears essential even for our models with cool effective temperatures.

Owocki et al. (2004) and van Marle et al. (2008) studied the winds of stars that reach or exceed the Eddington limit at their surface. As we have shown above, this condition is generally not found in our models (cf. Figs. 7 and 8). However, it may occur in helium-rich stars (see again Fig. 7) and hydrogen-free Wolf-Rayet stars (cf. Heger & Langer 1996), as well as in stars that deviate from thermal or hydrostatic equilibrium. Noticeably, Owocki et al. (2004) find that the mass-loss rates in this case are still quite limited, because of the energy loss attributed to lifting the wind material out of the gravitational potential (see, Heger & Langer 1996).

We want to emphasize in this context that the Eddington limit investigated in the quoted models as well as in our own may be different from the true Eddington limit, because of a number of effects that are all related to the opacity of the stellar matter in the stellar envelopes. One is that convection, which is necessarily present in the layers near or above the Eddington limit, may induce density inhomogeneities or clumping that can alter the effective radiative opacity (Shaviv 1998). In fact, depending on the nature of the clumping, the opacity may be enhanced (Gräfener et al. 2012) or reduced (Owocki et al. 2004; Ruszkowski & Begelman 2003; Muijres et al. 2011). Furthermore, such opacity calculations are tedious, and even in the currently used opacities, important contributions might still be missing.

Finally, the effect of stellar rotation on the stability limit in the atmospheres especially of hot stars is clearly important (Langer 1997, 1998; Maeder & Meynet 2000). However, it adds another dimension to this difficult problem and is therefore generally not included (cf. Sect. 7.2).

7.2. Stellar interior models

The peculiar core-halo density structure of inflated stars was first pointed out by Stothers & Chin (1993), after Iglesias et al. (1992) found the large iron bump in the opacities near 170 000 K. Further studies pointing out this phenomenon comprise Ishii et al. (1999), Petrovic et al. (2006), Gräfener et al. (2012), and Köhler et al. (2015). Conceivably, inflation may be present in further models of very massive stars, but often no statements on the presence or absence of this phenomenon are made in these papers.

For example, the models for very massive stars by Yusof et al. (2013) only discuss the electron-scattering Eddington factor in their models. On the ZAMS, their models are hotter and more compact than those of Köhler et al. (2015) analysed here, which implies that inflation is either weaker or absent. This difference might be due to the different treatment of convection in the sub-surface convective zones, where Yusof et al. (2013) assume the mixing length to be proportional to the density scale height instead of the standard pressure scale height. This prohibits the formation of density inversions, and since the density scale height tends to infinity when a model attempts to establish a density inversion, the convection may transport an arbitrarily large energy flux in this scheme. While the physics of convection introduces one of the biggest uncertainties in the atmospheres of stars close to the Eddington limit, efficient convective energy transport in inflated envelopes appears unlikely (cf. Sect. 5).

A suppression of inflation may have significant consequences for the evolution of massive stars, as the stellar models stay bluer and as a result have lower mass-loss rates and lower spin-down rates. The final fates of these non-inflated stars is significantly different compared to inflated stars (see Köhler et al. 2015, for a detailed discussion).

Gräfener et al. (2012) find inflation which, for their models without clumping, correspond well to those of Petrovic et al. (2006) for the Wolf-Rayet case, and to our unpublished solar metallicity main-sequence models, which show a bending of the ZAMS to cool temperatures for M ≳ 100 M_⊙. The models of Ishii et al. (1999) also agree very well. Including the work of Stothers & Chin (1993), we conclude that the effect of inflation in models of massive main-sequence stars is found in at least four independent stellar structure codes, with three of them quantitatively producing very similar results.

As pointed out above, massive star evolutionary models, which include effects of rotation, are being produced routinely these days (cf. Maeder & Meynet 2010; Langer 2012; Chieffi & Limongi 2013), but an investigation of the effect of stellar rotation on the stability limit in the atmospheres of hot stars requires the construction of two-dimensional stellar models.

8. Comparison with observations

8.1. The VFTS sample

A prime motivation of Köhler et al. (2015) for computing the evolutionary models for the very massive stars analysed here was to provide a theoretical framework for the VLT Flames-Tarantula Survey (VFTS, Evans et al. 2011). Within VFTS, multi-epoch spectral data of about 700 early B and 300 O stars are being analysed through detailed model atmosphere calculations. Within this effort, Bestenlehner et al. (2014) and McEvoy et al. (2015) derived the physical properties of more than 50 very massive stars, with luminosities log (L/L_⊙) > 5.5. We confront the models of Köhler et al. (2015) with this sample in Fig. 23.

Two sets of model data are included in Fig. 23, one that uses the effective temperatures of the Köhler et al. models directly, and a second one where the effective core temperature is used as defined in Sect. 4 (cf. Eq. (10)). The latter approximates the surface temperature of our models if inflation was completely absent. An example calculation presented in Appendix B, where inflation in a 300 M_⊙ is suppressed by increasing the mixing length parameter, shows that this approximation is indeed quite good. The ZAMS is also drawn for both sets of models. Note that while wind effects are clearly seen in the spectra of all stars in the sample, the optical depth of their winds is expected not to exceed τ = 2 (cf. Fig. 7 in Köhler et al. 2015) until the stars become very helium rich at their surface. Therefore, the effective temperatures derived from the observations need not be corrected for optically thick winds.

As shown in Bestenlehner et al. (2014), the hottest stars in their sample follow the ZAMS of the Köhler et al. models very closely, well into the regime of inflation. One might expect un-evolved stars to the left of the Köhler et al. ZAMS if inflation was not present. In that case, the stars above log L/L_⊙ ≃ 6.2 might spend a significant fraction of their lifetime on the hot side of the Köhler et al. ZAMS. The absence of these hot stars, however, does not conclusively argue that inflation does exist in nature, i.e. even without inflation, the star formation history in 30 Doradus might preclude the existence of stars like this, or they may be hidden in their natal cloud because of their youth (Yorke 1986; Castro et al. 2014).

On the cool side, Fig. 23 shows an absence of observed stars for log L/L_⊙ ≳ 6.15 and T_eff ≲ 35 000 K. As the evolutionary models predict about 30% of the core hydrogen burning to take place at T_eff ≲ 35 000 K in this luminosity regime, which may indicate that the inflation in the models of Köhler et al. is too strong. On the other hand, again, the absence of correspondingly cool stars 30 Doradus may be a result of the local star formation history.

In the luminosity range below, at 5.5 ≲ log L/L_⊙ ≲ 6.15, stars as cool as T_eff ≲ 15 000 K are observed for which McEvoy et al. (2015) concluded that they are still core hydrogen burning objects. The observed stars are somewhat cooler than the coolest core effective temperature of our models, which may argue in favour of inflation in real stars. Note that the life time of stars in the regime T_eff< 20 000 K is only 10% for the Köhler et al. models in the considered luminosity range.

In summary, as the stellar evolution models for these high masses are still quite uncertain, we do not find possible to argue for or against inflation being present in the observed stars considered here based in Fig. 23. In fact, it is intriguing that most of the observed stars are found in the regime where the inflated and non-inflated models overlap. Nevertheless, the observed sample above log L/L_⊙ ≃ 5.5 might constitute the best test case, since according to our models, the envelopes of all of them are expected to be strongly affected by the Eddington limit. Model atmosphere calculations for these stars, which include inflation, might shed new light on this question.

8.2. Further possible consequences of inflation

S Doradus type variability

Gräfener et al. (2012) argue that the S Doradus type variability of LBVs may be related to the effect of inflation, and, in particular, focussed on the case of AG Car (Groh et al. 2009). They propose that an instability sets in when their 70 M_⊙ chemically homogeneous hydrostatic stellar model is highly inflated (≈ 120 R_⊙) by virtue of which the inflated layer becomes gravitationally unbound and a mass-loss episode follows.

In contrast to this idea, our inflated hydrodynamic stellar models do not show any signs of this kind of instability. This could be because of various simplifying assumptions made in the models of Gräfener et al. (2012), in particular, their neglect of the convective flux in the inflated envelopes (cf. Sect. 6). Nevertheless, the physics of convection is very complex in these envelopes, and our results do not necessarily imply that instabilities do not occur.

In fact, considering the hot edge of the S Doradus variability strip, according to Smith et al. (2004), Fig. 22 shows that it roughly separates the models with a low maximum convective efficiency (η_max< 0.2) from those with a higher convective efficiency. If these high fluxes could not be achieved in these envelopes (cf. Sect. 6), a dynamical instability might well be possible.

The hot edge of the S Doradus variability in the HR diagram also coincides quite well with the borderline separating mildly (Δr/r_core< 1) from strongly (Δr/r_core> 1) inflated models. Comparing this with the observed distribution of very massive stars in Fig. 23, which indicates that essentially no stars are found far to the cool side of this line, could indicate again that strongly inflated envelopes are indeed unstable and might lead to S Doradus type variability and an increased time-averaged mass-loss rate.

LBV eruptions

Glatzel & Kiriakidis (1993) speculated that strange mode pulsations might be responsible for the LBV phenomenon. These pulsations are characterized by very short growth times (~ τ_dyn) and small brightness fluctuations roughly of the order ~ 10...100 mmag (Glatzel et al. 1999; Grott et al. 2005). However, the mass contained in the pulsating envelopes of their models is negligible compared to the stellar mass, and the associated brightness variations cannot explain the humongous luminosity variations observed in LBV eruptions.

We have seen in Sect. 5 that an inflated envelope often produces a density inversion. These density inversions have been repeatedly proposed as a source of instabilities giving rise to eruptive mass-loss in LBVs (Maeder 1989; Maeder & Conti 1994; Stothers & Chin 1993). Given our results, we consider it unlikely that a density inversion can be the sole cause of LBV eruptions. Density inversions are a generic feature present in a multitude of our models (see Fig. 16), while the LBV phenomenon is quite rare. Furthermore, the density inversions in our models are found very close to the surface of the star with very small amounts of mass above the inversion.

Given our results, inflation per se appears unlikely to cause LBV eruptions, again, because it occurs too abundantly in our models, and also because the mass of the inflated envelope is generally very small. However, Fig. 24 reveals that this is not so for our coolest models. Whereas for most models the mass of the inflated envelope is smaller than ~ 10^-3M_⊙, intriguingly it rises to several solar masses in the models that have effective temperatures below ~ 10 000 K. These cool models, of which detailed examples are presented in Appendix D, also show the highest Eddington factors (Fig. 3) and the strongest inflation (Fig. 11). This behaviour is seen in the mass range of 40...100 M_⊙, which corresponds well to the masses of observed LBVs.

A key feature in our cool models with massive inflated envelopes is visible in Fig. 21. As the opacity in the hydrogen recombination zone becomes very large, effectively blocking any radiation transport, convective efficiencies of the order of η ≃ 1 are needed to transport the stellar luminosity through this zone. Also, in the iron convection zone at the bottom of the envelope, a high convective efficiency (η ≃ 0.5) is found in these models. As shown in Figs. 19 and 20, this requires sonic or even supersonic convective velocities in the framework of the standard MLT as implemented in our code. It thus appears conceivable that in reality convection is less efficient in this kind of a situation, e.g. because of viscous dissipation. When a star enters this region with a massive inflated envelope, throughout which the stellar luminosity can neither be transported by radiation nor by convection, hydrostatic equilibrium is possible any more, and the loosely bound inflated envelope may by dynamically ejected. We believe that this scenario may relate to major LBV eruptions.

As proposed by Langer (2012), a rapid evolutionary timescale of a star may be required to obtain LBV outbursts in addition to the star reaching the Eddington limit. If evolution on the thermal timescale is rapid enough, this might produce LBVs after core hydrogen exhaustion. This rapid evolution may relate to most of the observed LBVs in our Galaxy and the Magellanic Clouds, as well as LBVs after core helium exhaustion, which may address the recently accumulated evidence of LBV outbursts in immediate supernova progenitors (Smith & Arnett 2014).

Supernova shock break-out

A recent study by Moriya et al. (2015) concluded that inflated stellar models can help to explain the extended rise time of the shock break-out signal from the Type Ib supernova SN2008D (Soderberg et al. 2008). In this scenario, the shock break-out occurs deep inside the inflated envelope and consequently the rise time is determined by the radiative diffusion time of the envelope and not the light crossing time. They also noted that more of these kinds of events, if observed in the future, might serve as indicators of inflated supernova progenitors.

Whereas the above result supports the idea that hydrogen-free Wolf-Rayet stars may possess inflated envelopes (Petrovic et al. 2006; Gräfener et al. 2012), LBVs have also been considered to be immediate progenitors of supernovae (Kotak & Vink 2006; Groh et al. 2013). In the realm of high cadence supernova surveys, this opens up the possibility to also test the existence of envelope inflation in hydrogen-rich stars through supernova shock break-out observations.

9. Discussion and conclusions

We investigated the internal structures of the massive star models computed by Brott et al. (2011) and Köhler et al. (2015) using a 1D hydrodynamical stellar evolution code, with particular emphasis on the Eddington limit. We find that the conventional idea of sufficiently massive stars reaching the Eddington limit at the stellar surface is not reproduced by our core hydrogen burning models, not even at 500M_⊙ (cf. Figs. 7 and 8). Instead, we find a suitably defined Eddington limit inside the star (Eq. (3)) is reached by models with log (L/L_⊙) ≳ 5.6 (Fig. 2), which leads to sub-surface convection, envelope inflation (Fig. 13), and possibly to pulsations. Many of our models even exceed this Eddington limit, and in the extreme case of red supergiants even by factors of up to seven (Fig. 3), with the consequence that strong density inversions develop such that hydrostatic equilibrium is maintained (Fig. 17).

In the analysed models, whose initial composition is chosen to match that of the LMC, all stars above ~ 40 M_⊙ reach the Eddington limit in their envelopes. As iron opacities are mainly responsible for this phenomenon, we expect that this mass limit is higher at a lower metallicity, and similarly lower for massive stars in our galaxy. Furthermore, there may be two groups of stars for which this limit comes down even further. First, the centrifugal force in rapidly rotating stars may lead to similar conditions in the envelope layers near the stellar equator, i.e. to a strong latitude dependence of inflation. Perhaps, this could give rise to the so-called B[e] supergiants, which show a slow and dense equatorial wind and a fast polar wind at the same time (Zickgraf et al. 1985). Second, the mass losing stars in interacting close binary systems evolve to much higher L/M-values than corresponding single stars (Langer & Kudritzki 2014), and are therefore expected to reach the Eddington limit for much lower initial masses.

The stability of the inflated envelopes is not investigated here, but many of them are likely to be pulsationally unstable (Glatzel & Kiriakidis 1993; Saio et al. 1998, Sanyal et al., in prep.). If so, it is expected that the pulsations will lead to mass-loss enhancements (e.g. Moriya & Langer 2015), or to the loss of the inflated envelope. In the latter case, the envelope is expected to re-grow unless the achieved time average mass-loss rate exceeds the high critical mass-loss rate (Sect. 4.2). We find that in our coolest models, the mass contained in the inflated envelopes can reach several solar masses (Fig. 24), and speculate that their dynamical loss may resemble LBV major eruptions (cf.  Sect. 8.2). Consequently, even though reaching or exceeding the Eddington limit may not immediately lead to strong outflows in stars, clearly the mass-loss rate of the stars is strongly affected, in the sense that the mass-loss is significantly enhanced one way or another.

It will be crucial to test observationally whether luminous, main-sequence stars indeed possess inflated stellar envelopes. This possibility has not yet been investigated with stellar atmosphere models for hot stars. Perhaps the best candidates are the S Doradus variables (Gräfener et al. 2012), which appear in the part of the HR diagram where our models predict a radius inflation by more than a factor of two (cf. Fig. 10).

Finally, we note that besides massive stars, the Eddington limit is relevant to various other types of stars, as luminous post-AGB star, X-ray bursts, Novae, R Corona Borealis stars, and accreting compact objects. It may be interesting to assess to what extent similar phenomena as found in this work might play a role in these objects.

Online material

Appendix A: Interior structure of a 85 M_⊙ stellar model

Appendix B: Effect of efficient convection on inflation

Knowing that convective flux is proportional to the mixing length, we show here (Fig. B.1) that by increasing the mixing length parameter α in an inflated 300 M_⊙ model near the ZAMS, the inflation gradually goes away and what we are left with is an almost non-inflated star, whose radius is well-approximated by core radius r_core of the inflated model.

Appendix C: Convective velocity profile in a WR model

In Fig. C.1, the convective velocity is shown as a function of radius in a massive (147 M_⊙) WR-type (Y_S = 0.89) stellar model.

The variation of isothermal and adiabatic sound speeds are also plotted for comparison. The convective velocities exceed the local isothermal sound speed in the envelope where conditions are non-adiabatic, i.e. the thermal adjustment time is short. In these models, turbulent pressure becomes important (which is not taken into account in our models) as well as standard MLT fails to be a good approximation for modelling convection.

Appendix D: Representative models

The profiles of different relevant physical quantities are shown for a few selected stellar models at five distinct effective temperatures corresponding to the three peaks in Γ_max and the two troughs in between the peaks (cf. Fig. 2).

Acknowledgments

We thank the referee, Raphael Hirschi, for useful comments that improved the manuscript. We are grateful to S. C. Yoon for many fruitful discussions. We also thank C. J. Evans, D. Szécsi, and J. S. Vink for helpful comments on the manuscript. L.G. thanks the International Max Planck Research School for Astronomy and Astrophysics at Bonn.

References

All Figures