© ESO, 2012

1. Introduction

In the standard picture for the evolution of the most massive stars (with M ≳ 30   M_⊙), O stars have been proposed to evolve through successive luminous blue variable (LBV) and Wolf-Rayet (WR) phases before exploding as hydrogen (H) free supernovae (SNe) Ibc (e.g., Langer et al. 1994; Meynet & Maeder 2003; Heger et al. 2003). This evolutionary path is however still open as LBVs have more recently been suggested to be the direct progenitors of some H-rich type II SNe (Kotak & Vink 2006; Smith et al. 2007; Gal-Yam & Leonard 2009). Additional renewed interest in WR stars arises from their connection to long-duration gamma ray bursts (long GRBs, cf. Yoon & Langer 2005; Woosley & Heger 2006).

What LBVs and WR stars have in common is their proximity to the Eddington limit. Not only is this thought to be instrumental for their mass-loss behavior (Vink & de Koter 2002; Gräfener & Hamann 2008; Vink et al. 2011; Gräfener et al. 2011), but it may also have key consequences for their stellar structures. Using modern OPAL opacities, Ishii et al. (1999), Petrovic et al. (2006), and Yungelson et al. (2008) studied the interior structure of massive stars, revealing a “core-halo” configuration that involves a relatively small convective core and an extended radiative envelope. This is referred to as the “inflation” of the outer envelope.

It has been known for many decades that the observed radii of WR stars are almost an order of magnitude larger than canonical WR stellar structure models indicate. This is oftentimes attributed to their so-called pseudo-photospheres, which concern an “effective” photosphere several times larger than the hydrostatic radius, as a result of a dense stellar outflow (e.g., Crowther 2007). Pseudo-photospheres have also been discussed in the context of LBVs, although the issue is still under debate (e.g., Smith et al. 2004).

An alternative explanation for the large radii of WR stars may involve envelope inflation. The issue is particularly relevant in the context of WR stars as GRB progenitors, as there have been several suggestions that the required radii of GRB progenitors are up to an order of magnitude (~1   R_⊙ vs.  ~ 10   R_⊙) smaller than the radii of observed WR stars (Modjaz et al. 2009; Cui et al. 2010). The WR radius issue is thus highly relevant for GRB modeling.

For LBVs the issue of their radii is equally interesting. The defining property of LBVs concerns their “S Doradus” cycles. On timescales of years to decades LBVs are seen to vary between effective temperatures of  ~30 kK (early B spectral type) to  ~8 kK (early F). Due to the lack of convincing counter-evidence, these S Dor excursions in the upper Hertzsprung-Russell diagram have generally been assumed to occur at constant bolometric luminosity (Humphreys & Davidson 1994), but recent studies have challenged this behavior for a couple of objects (Groh et al. 2009b; Clark et al. 2009). Either way, the issue of the increased LBV radii during S Dor cycles remains un-questioned, but a satisfying explanation for it has yet to be provided (Vink 2009). Outer envelope inflation may turn out to be an interesting explanation for S Dor variations, as will be detailed in the following.

The paper is organized as follows. In Sect. 2 we present numerical stellar structure models that show the envelope inflation effect. In Sect. 3 we describe the effect analytically, including a new instability limit. We also provide a recipe that can be used to estimate the inflation effect for arbitrary stars with given (observed) parameters. The implications of our results are discussed in Sect. 4, with focus on WR stars and LBVs. The conclusions are summarized in Sect. 5.

2. Model computations

In this section we present stellar structure models for chemically homogeneous stars close to the Eddington limit. To model the conditions in classical WN-stars (i.e., WR stars in the phase of core He-burning that have lost their H-rich envelope), we use pure He-models, while models with hydrogen resemble the conditions in massive WNh stars, and LBVs. At this point we note that we only use the assumption of chemical homogeneity, to keep the models as simple as possible. The occurrence of an envelope inflation is not connected to this assumption, as it only affects the outermost layers of the star.

In Sect. 2.1 we give an overview of our numerical approach for computing the stellar structure, and in Sect. 2.2 we describe the envelope inflation, as it occurs in our models. In Sect. 2.3, we present a model grid to investigate the occurrence of this effect in the HR diagram. Finally, in Sect. 2.4 we discuss the important influence of density inhomogeneities on our models.

2.1. Numerical method

Our models are computed with a code that integrates the stellar structure equations using a shooting method. The numerical methods are described in the textbook by Hansen & Kawaler (1994). The code is based on an example program that is distributed with the book, but is completely re-written, and updated with OPAL opacities (Iglesias & Rogers 1996). As the chemical structure of the stars is fixed, a set of four differential equations needs to be solved. These are the equations of hydrostatic equilibrium (1)mass conservation (2)energy conservation (3)and energy transport (4)where ∇ = ∇_rad for the case of purely radiative energy transport, and ∇_rad > ∇ > ∇_ad for the convective case. As our models are pure stellar structure models, i.e. they include no time-dependent terms, the energy generation rate ϵ in Eq. (3) does not take changes in the gravitational energy, due to contraction/expansion of the star into account. This omission has no consequences for the physics of the low-density envelopes discussed in this work, but can potentially affect the L/M ratio of WR stars in late evolutionary stages (see Sect. 3.5.2, for a detailed discussion).

Furthermore, we neglect dynamic terms due to mass loss, i.e., we only concentrate on the hydrostatic case. The dynamic case has been discussed by Petrovic et al. (2006). For the case of WR stars with strong winds, Petrovic et al. found that dynamic effects may inhibit the formation of inflated envelopes. We will address this point in more detail in Sect. 3.5.3.

The choice of the outer boundary condition, particularly the temperature at the outer boundary, turns out to be important for the formation of inflated envelopes. Here, we prescribe ρ, and T at the outer boundary of the star in a way, that resembles reasonable values for the sonic points of stars with strong, radiatively driven winds.

An envelope inflation only occurs in our models if we choose outer boundary temperatures below  ~70 kK, well below the temperature of the hot Fe-opacity peak. Such temperatures are expected at the sonic points of optically thick, radiatively driven winds, such as the late-type WN stars described by Gräfener & Hamann (2008). The density ρ is chosen such that the resulting mass loss rate, (5)lies in a realistic range.

We note that our envelope solutions are very robust with respect to the detailed choice of ρ, and T, as long as the temperature T lies clearly below the regime of the hot Fe-opacity peak, which has its maximum at  ~160 kK. The effect of changes in the outer boundary condition is discussed in more detail in Sect. 3.5.1.

The stellar temperatures T _⋆  given in this work, are effective temperatures related to the outer boundary radius R _⋆  of our models, i.e., . For static atmospheres these values will be almost identical to classical effective temperatures T_eff. In the presence of dense, extended stellar winds, our T _⋆  values will resemble the “effective core temperatures”  T _⋆ , that are typically inferred from models for expanding atmospheres (cf. Gräfener et al. 2002).

Fig. 1 Density vs. radius for our 23   M⊙ He model (cf. Table 1).

Fig. 2 Pressure vs. radius. Gas pressure Pgas, radiation pressure Prad, and total pressure Ptot = Pgas + Prad.

2.2. Stellar structure models with inflated envelopes

As a typical example for a model with an inflated envelope, we present a homogeneous stellar structure model for a 23 M_⊙ helium star with solar metallicity (X = 0, Y = 0.98, Z = 0.02). In agreement with previous works by Ishii et al. (1999), and Petrovic et al. (2006), this model develops an extended low-density envelope, with a density inversion close to the surface (cf. Fig. 1). The star thus forms a shell around the actual stellar core, with a radius that is 3–4 times larger than the core radius.

Figure 2 shows that the inflated region is dominated by radiation pressure. In this model, the gas pressure P_gas only accounts for a fraction of  ~4 × 10^-4 of the total pressure P. In contrast to this, P_gas and P_rad are roughly equal in the stellar core. The core is thus only fairly close to the Eddington limit, as reflected by an Eddington parameter of Γ_e ≈ 0.4 (cf. Table 1). Here Γ_e denotes the “classical” Eddington parameter, which is related to the Thomson opacity κ_e due to free electrons, i.e., Γ_e = κ_eL_rad/(4πcGM). The Γ_e given in Tables 1 and 2, are computed for a fully ionized plasma with a given hydrogen mass fraction X, and L_rad = L. Γ_e then simplifies to (6)In contrast, the total Eddington factor Γ = κL_rad/(4πcGM), which is related to the total opacity κ, approaches unity in the inflated zone. This is shown in Fig. 3, where we compare the total luminosity within our model with the Eddington luminosity L_Edd = 4πcGM/κ (the condition L_rad = L_Edd is thus equivalent to Γ = 1). In Fig. 3 we can identify three zones within the star. In the inner, convective core, the total luminosity L is larger than the Eddington luminosity L_Edd. The main reason for this is that the energy transport is almost entirely convective, and thus L_rad is very low. Above this region follows the radiative envelope with L < L_Edd. On top of this follows the inflated zone, where L_Edd almost equals L. Γ is thus extremely close to one, almost throughout the complete outer envelope of the star. This is surprising, because Γ is a function of the opacity κ, which is a strongly varying function of ρ, and T. The conditions in the inflated envelope thus demand for a mechanism that adjusts Γ very precisely to one.

To achieve a large radial extension it is necessary that the density scale height H is comparable to the stellar core radius R_c. For our example model we have R_c/H = GM(1 − Γ)/(a²R_c) ≈ 3 × 10³(1 − Γ). From Fig. 2 we can infer that (1 − Γ) ≈ P_gas/P ≈ 4 × 10^-4 (cf. Eq. (15)). Due to the low ratio of gas pressure to radiation pressure, Γ is thus just close enough to one, to achieve a large radial extension.

A qualitative idea about the reason for the envelope inflation can be obtained from Fig. 4, where we plot the gas pressure P_gas, and the total pressure P, vs. the radiation pressure P_rad. This diagram is particularly useful, as for stars close to the Eddington limit, P is of the same order of magnitude as P_rad, i.e., in a logarithmic diagram the stellar structure just follows a straight line (in accordance with the Eddington standard model (Eddington 1918) for which P_rad/P = const is adopted). For P_gas the situation is more complicated. In the stellar core, P_gas follows a similar relation as P. In the extended envelope, P_gas however drops significantly, and rises again after a certain minimum pressure P_min is reached. As demonstrated in Fig. 5, the region where P_gas drops corresponds precisely to the location of the Fe-opacity peak in the P_rad–P_gas plane.

More explicitly, starting from the stellar core, the solution follows a straight line with P_gas/P_rad = const. This naturally leads into a region in the P_rad–P_gas plane, where the opacity κ increases significantly, due to the presence of the Fe-opacity peak. When κ reaches κ_Edd, this leads to Γ → 1. To avoid a super-Eddington situation, with Γ > 1, either L_rad, or κ need to be reduced. L_rad could in principle be reduced by convection, however, in our case the involved densities turn out to be too low, i.e., convection is inefficient (cf. Sect. 3.5.4). The only possibility to reduce κ, is thus to lower the density. For this reason the density has to drop in our model, until the tip of the Fe-peak is reached (cf. Fig. 5). At the same time, however, this leads to a situation where P_gas ≪ P_rad, and thus to Γ → 1 (because, as already mentioned above (1 − Γ) ≈ P_gas/P). This mechanism thus forces the solution to follow a path with Γ = 1, i.e., with κ(ρ,T) = κ_Edd.

As P_rad, and P_gas are functions of ρ, and T, we can express κ in the form κ(P_gas,P_rad). This is actually done in Fig. 5, where the colors indicate κ(P_gas,P_rad)/κ_Edd on a logarithmic scale, as obtained from the OPAL opacity tables (Iglesias & Rogers 1996). In this plot we can see that the solution within the low-density zone is completely dominated by the topology of the opacity κ, i.e., by material properties. As soon as the solution reaches low densities close to the Fe-opacity peak, it just follows a path with Γ = 1, until a minimum density ρ_min is reached. This minimum is located at the tip of the Fe-opacity peak, typically at P_rad ≈ 10⁶   dyn   cm^-2, corresponding to a temperature of  ~ 150 kK. After the opacity peak is passed, the solution still follows Γ = 1, but now with increasing P_gas, leading to a density inversion.

Fig. 3 Total stellar luminosity L(r), compared to the Eddington luminosity LEdd = 4πcGM/κ.

Fig. 4 Stellar structure in the Prad–Pgas plane. Gas pressure Pgas, radiation pressure Prad, and total pressure Ptot = Pgas + Prad are plotted vs. Prad throughout the whole star.

To maintain the described envelope inflation, it is thus necessary that 1) the star is close enough to the Eddington limit, so that Γ → 1; and 2) convection is inefficient. If the envelope inflation occurs, its properties are determined by the topology of the Fe-opacity peak in the P_rad–P_gas plane, i.e., by material properties. In Sect. 3, this will help us to find an analytical description of the inflation effect. It also implies a dependence on material properties like the metallicity Z, as discussed by Ishii et al. (1999), and Petrovic et al. (2006). In Sect. 2.4, we will discuss an additional important effect, namely the influence of density inhomogeneities within the inflated zone, which is based on a similar principle.

2.3. A grid of homogeneous stellar structure models with envelope inflation

To assess the effect of the envelope inflation in the HR diagram, we have computed a grid of homogeneous stellar structure models with hydrogen mass fractions X = 0.0, 0.2, 0.4, and 0.7, that covers a luminosity range between log (L/L_⊙) = 5.0, and 6.5. The results are compiled in Tables 1 and 2. In Fig. 6, we show the position of our models in the HR diagram, together with theoretical zero-age main sequences (ZAMS) that do not take the inflation effect into account.

At low luminosities, our He-burning models (with X = 0.0) closely resemble the He-ZAMS from Langer (1989b), and our H-burning models (X = 0.2–0.7) the ZAMS from Schaller et al. (1992). Towards higher luminosities, our models develop inflated envelopes, i.e., they show lower T _⋆ . With increasing luminosity, T _⋆  decreases rapidly until, above a certain luminosity, we do not find static solutions anymore. The reason for the occurrence of this important limit will be explained in Sect. 3.2.

In our present models, the envelope inflation sets in for Eddington factors Γ_e between 0.3 (X = 0.0), and 0.4 (X = 0.7), i.e., within a relatively narrow range. The fact that it sets in earlier for H-deficient stars is a consequence of the fact that the L/M ratios of these objects are higher (although κ_e is lower for lower H abundances). In Sects. 3.2, and 3.3 we will show that this behavior, as well as the described upper luminosity limit, are a consequence of the strength, and shape of the Fe-opacity peak in the P_rad–P_gas plane.

2.4. The influence of density inhomogeneities

For our pure He models in Table 1, the described envelope extension starts at luminosities of log (L/L_⊙) ≈ 5.5, and fully develops for log (L/L_⊙) ≈ 5.8 (cf. Fig. 6). As we will discuss in Sect. 4.1, the largest part of the Galactic H-free Wolf-Rayet stars, however, seems to show inflated envelopes for luminosities just below this value. In the following we will show that the inflation in this parameter range can be explained by density inhomogeneities. Such inhomogeneities are expected to arise in the outer envelopes of WR stars, e.g. due to strange-mode instabilities (Glatzel & Kaltschmidt 2002; Glatzel 2008).

If material is inhomogeneous, or clumped, the density within clumps is larger than the mean density. As a result, the opacity κ has to be evaluated for an increased density, instead of the mean density (note that κ denotes the mass absorption coefficient in cm²/g). Here we assume a constant clumping factor D within the inflated envelope, and evaluate κ(ρ,T) for the enhanced density D × ρ, instead of the mean density ρ. This is equivalent to assuming a medium with clumps of a constant density D × ρ, and an inter-clump medium which is void. In this picture the clumps would thus only fill a fraction f = 1/D of the total volume. The same definition of clumping factors D, and filling factors f is commonly used for the modeling of inhomogeneous stellar winds (e.g. Hamann & Koesterke 1998).

In the clumping approximation used here, it is assumed that no radiative flux gets lost in between clumps, e.g. by shadowing/porosity effects. Such a situation could be achieved by optically thin clumps, or by a clump geometry (such as large shells, or pancakes) that is not sensitive to shadowing effects.

Fig. 5 Envelope solution in the Prad–Pgas plane. The colors indicate the logarithm of the Eddington factor Γ, for given Prad and Pgas, according to the OPAL opacity tables (Iglesias & Rogers 1996), for our pure He model with 23   M⊙. The envelope solution almost precisely follows a path with Γ = 1. The radii of 2, 4, and 6   R⊙ within the inflated envelope are indicated. Arrows indicate the slope of the solution, as expected from Eq. (12).

Following our argumentation from the previous section, large envelope extensions, i.e., large scale heights H, are reached for Eddington factors Γ that are very close to one. We have shown that this is achieved for the low values of P_gas at the tip of the Fe-opacity peak, as P_gas/P ≈ (1 − Γ). By the assumption of clumping, the opacity peak is shifted towards even lower mean densities, i.e., towards lower (mean) values of P_gas. In this way, the envelope inflation is further enhanced (a more concise explanation of this effect is given in Sect. 3, Eqs. (19), (20) (29)).

In Fig. 6 we show He-ZAMS models that are computed with clumping factors D = 4, and D = 16. As expected, the envelope extension occurs much earlier in these models. In Sect. 4.1 we will show that these models cover the observed HR diagram positions of the Galactic H-free WR stars, i.e., we can explain the observed temperatures of these objects with moderate clumping factors. The adopted values for D are in notable agreement with spectroscopically determined clumping factors for the winds of WR stars, e.g. by Hamann & Koesterke (1998).

3. Analytical description of the envelope inflation

In the previous section we have shown that the properties of inflated envelopes are chiefly determined by the fact that the solution has to follow a path with Γ = 1 in the P_rad–P_gas plane (cf. Fig. 5). In the following we elaborate on this, and derive an analytical description of this process. As we concentrate on the physics of the inflated envelopes alone, the resulting relations do not depend on the internal structure of the star, i.e., they are generally applicable to stars with given (observed) parameters. In Sect. 3.1, we start with the equations describing the envelope structure. In Sect. 3.2 we derive analytical expressions for the radial extension of the envelope ΔR, and in Sect. 3.3 we present a recipe to estimate ΔR based on observed/adopted stellar parameters. Furthermore, we derive an estimate for the mass of the surrounding shell, ΔM, in Sect. 3.4. Finally, we discuss in Sect. 3.5, to which extent the inflation effect may be affected by model assumptions.

3.1. Inflated low-density envelopes near the Eddington limit

Fig. 6 HR diagram for pure helium models (red), and models containing hydrogen (blue).

In this section we concentrate on the physics of the outermost stellar envelope, where m = M, and L = const. With these assumptions, we eliminate two equations from the stellar structure equations (Eqs. (1)–(4)), so that only the equation of hydrostatic equilibrium (Eq. (1)), and the energy transport equation (Eq. (4)) are left.

In Sect. 2.2 we concluded, that two conditions are mandatory to maintain an inflated envelope. 1) the star needs to be close enough to the Eddington limit, so that Γ = 1 can be reached; and 2) convective energy transport needs to be inefficient. Under condition 2) (the validity of this condition is discussed in Sect. 3.5.4), the energy transport equation (Eq. (4)) becomes purely radiative, i.e., L = L_rad, and ∇ = ∇_rad, with (7)After substituting Eqs. (1), and (2) in Eq. (4), we thus obtain the energy transport equation in the form (8)If this relation is written in terms of P_rad = (a/3)T⁴ we obtain (9)for the (outward directed) radiative acceleration g_rad. Together with the equation of hydrostatic equilibrium (10)we obtain by division of Eqs. (10) and (9) (11)or (12)Notably, ∂P_gas/∂P_rad depends only on the Eddington factor Γ, i.e., via Γ/Γ_e = κ/κ_e on Γ_e, and κ. This is in line with our previous finding that the envelope inflation is fully determined by the classical Eddington factor Γ_e (as given by Eq. (6)), which is a function of stellar parameters, and κ, which depends on material properties.

In Fig. 5 we compare our numerical solution from Sect. 2.2 (blue line) with the slopes ∂P_rad/∂P_gas (black arrows) following from Eq. (12), in the P_rad–P_gas plane. The background colors indicate the value of Γ = Γ_e(κ/κ_e) for our model, where κ has been extracted from the OPAL opacity tables (Iglesias & Rogers 1996).

According to Eq. (12), the point at which the envelope solution crosses the Eddington limit, i.e. where Γ = 1, needs to be an extremum in P_gas. Figure 5 nicely shows, that our numerical solution indeed approaches the Eddington limit, and crosses it at the point of the lowest possible gas pressure. If the Eddington limit would be crossed earlier, at a higher density, the gas pressure had to increase so fast, that extremely high densities would be required at the outer boundary. The only way to reach reasonably low temperatures, and densities at the outer boundary, thus leads around the Fe-opacity peak, via a low gas pressure, i.e., via low densities. The formation of a low-density envelope with a density inversion is thus a natural consequence of the topology of the Fe opacity peak in the P_rad–P_gas plane, in combination with Eq. (12).

For the example in Fig. 5, we find that P_gas/P_rad ≈ 4    ×    10^-4 at the tip of the Fe-opacity peak. In Sect. 2.2 we already discussed that this implies that Γ → 1, as P_gas/P ≈ (1 − Γ). This relation follows from Eq. (11), if we assume a slowly changing ratio P_rad/P. In this case we have ∂(lnP_rad)/∂(lnP) ≈ 1 (cf. Fig. 4). It thus follows that (13)which leads to (14)or (15)As discussed in Sect. 2.2, this small value of (1 − Γ) leads to a density scale height of the order of the stellar radius, i.e., to an envelope inflation.

The occurrence of an envelope inflation is thus intimately connected to the fact that P_gas ≪ P_rad at the tip of the Fe-opacity peak. As noted earlier, the location of this point in the P_rad–P_gas plane is only determined by the topology of the opacity peak, i.e., by material properties. In the following we will take advantage of this fact, to obtain analytical estimates of the radial envelope extension ΔR.

3.2. Computation of the radial extension ΔR

To compute the radial extension ΔR, we integrate Eqs. (9) and (10) from the outer boundary of the inflated zone, at the “envelope radius” R_e, to the inner boundary, at the “core radius” R_c. As this integration is performed within the inflated zone, where P = P_rad, and κ = κ_Edd, Eqs. (9) and (10) degenerate to one equation, that describes the change of P_rad, and thus also the change of the temperature, in hydrostatic equilibrium. (16)or (17)The integration is performed from R_e, where P ≡ P_e, to R_c, where P ≡ P_c, leading to (18)The left hand side of this equation can be written in the form (19)where ΔP = P_c − P_e, and is the mean of 1/ρ with respect to P_rad. In Sect. 3.3 we will show that ΔP, and can be estimated alone on the basis of opacity tables. For given W_c, the radial envelope extension can then be computed from Eq. (18), which becomes (20)where the latter equality uses a definition of W in terms of the observable, outer-envelope radius R_e(21)The expression for the radial extension ΔR = R_e − R_c for given W_c, or W_e follows directly from Eq. (20) For the W_c, and W_e as determined from our models (Tables 1 and 2) these relations are fulfilled very precisely.

The parameters W_c/e are of essential importance for the inflation effect. As u = 3P/ρ in a radiation dominated gas (where u denotes the specific energy per volume), is related to the specific energy per gram of material within the inflated zone. According to Eq. (18), it denotes the increase of the specific energy, from the outer to the inner edge of the inflated envelope, divided by 3. The parameters W_c/e indicate the ratio between this energy increase, and the gravitational energy per gram of material, at R_c/e.

Notably, Eq. (20) imposes an instability limit for W_c → 1, for which R_e → ∞. This limit coincides with the limit that we have described in Sect. 2.3, where we found that no hydrostatic solutions exist, if certain luminosities are exceeded. Stable solutions only exist for W_c < 1, which is equivalent to (24)As is largely fixed by material properties, Eq. (24) imposes a limit on the gravitational energy at the core radius R_c. If the gravitational energy becomes too small (i.e., R_c becomes too large), the stellar envelope becomes energetically unbound, and there exists no hydrostatic envelope solution. This limit might be of profound importance for the occurrence of instabilities close to the Eddington limit, as observed, e.g. for LBVs (cf. Sect. 4.2).

The occurrence of the envelope inflation is thus mainly due to the fact that the Fe-opacity peak forces the envelope to low densities ρ, while the temperature of the opacity peak is fixed, i.e., u = const. This leads to a very high specific energy u/ρ, i.e., the envelope becomes less gravitationally bound.

3.3. A recipe to compute ΔR

Fig. 7 Plot of 1/ρ vs. Prad. The dashed horizontal lines indicate 1/ρmin, and 1/(10   ρmin). The intersections of the latter with the envelope solution represent boundary pressure Pc, and Pe, as indicated by vertical dashed lines. The grey shaded area illustrates that the integral ∫   dP/ρ ≈ ΔP/(2ρmin).

The remaining step for the computation of ΔR is to obtain estimates of the parameter W (Eqs. (19), (21)), which depends on the ratio M/R (where R denotes either R_c, or R_e), and the ratio . While M, and R are given stellar parameters, needs to be determined from opacity tables.

Based on Eqs. (18), (19) we have (25)In Fig. 7 we illustrate that, due to the topology of the Fe-opacity peak, the integral on the right hand side is roughly equal to ΔP/(2ρ_min), i.e., (26)with f_ρ ≈ 2. Numerically obtained values for f_ρ from our model computations confirm this value (cf. Tables 1 and 2). W can thus be computed on the basis of ΔP = P_c − P_e, and ρ_min.

To extract these values from the opacity tables, we make use the fact that the envelope solution follows almost precisely a path with Γ = 1, i.e., a path with κ = κ_Edd in the P_rad–P_gas plane. From this path we extract the minimum gas pressure P_min, at the tip of the Fe-opacity peak (cf. Fig. 5). We use the density at this point as our reference density ρ_min. For given P_min, we define the boundaries of the inflated envelope as the two adjacent points on the path with κ = κ_Edd, where P_gas = f_P·P_min (in this work we use f_P = 10). The radiation pressure at these two points gives us P_c, and P_e. Our choice of f_P = 10 has a moderate influence on the resulting values, will however not change the overall results.

The procedure described above thus provides ρ_min, and ΔP = P_c − P_e for given κ_Edd, or Γ_e (note that Γ_e(κ_Edd/κ_e) = 1), i.e., we can now compute W, and ΔR. To determine ΔP/ρ_min over a wide parameter range, we employ a root finding algorithm to extract ΔP, and ρ_min from the OPAL tables (Iglesias & Rogers 1996), as described above. The value of ΔP/ρ_min depends on the chemical composition, i.e., X and Z, and on Γ_e.

At this point we define a function Q(X,Z,Γ_e), that only depends on ΔP/ρ_min(27)Adopting a Galactic metallicity of Z = 0.02, we determine a fitting relation for Q(X,Γ_e), based on the numerically determined values of ΔP/ρ_min. The resulting relation has the form (28)The fits are displayed in Fig. 8, where the crosses denote the numerically determined values of Q, for the range X = 0–0.7, and the blue lines denote the fitting relation Eq. (28), for X =  0, and X = 0.7.

Fig. 8 Fitting relation for Q(X,Z = Z⊙,Γe). Black crosses indicate numerically determined values of Q, for X = 0 – 0.7. Blue lines indicate our fitting relation Eq. (28), for X = 0 (top), and X = 0.7 (bottom).

The values of Q determined this way, lie between Q ~ 1, and Q ~ 150 (cf. Fig. 8), with Γ_e = 0.3–0.6. For lower Q (i.e., lower Γ_e), the extension of the Fe-peak becomes so small that it does not cover a factor f_P = 10 in P_gas. Under these conditions an envelope inflation may still occur, however only of moderate strength. For Q larger than  ~150, ρ_min becomes so small that the range of the opacity tables is exceeded. In this case we still expect a very strong envelope inflation, but are just not able to infer Q from the tables.

For given Q, W can be determined from Eq. (26). If we additionally allow for clumping with a clumping factor D (cf. Sect. 2.4), we have , i.e., we can express W in the form (29)with f_ρ ≈ 2.

To compute the envelope extension ΔR for a star with given parameters M, L, X, and R (where R may denote R_c or R_e), the following steps thus need to be performed.

• 1)

  For given M, L, and X, Γ_e can be computed from Eq. (6).

• 2)

  For given X, and Γ_e, Q can be estimated from Eq. (28).

• 3)

  For given Q, W can be determined from Eq. (29), adopting reasonable values for f_ρ, and D. Note that the stability limit Eq. (24) can be expressed in the form Q/(f_ρ/D) < M/R (in solar units).

• 4)

  Finally, for given W, the radial extension ΔR can be computed from Eq. (22), or Eq. (23).

E.g., for our He models in Table 1 (with X = 0), masses are typically of the order of 15   M_⊙, and core radii are  ~ 1   R_sun, i.e., M/R ~ 15 in solar units. This means that an unstable situation, with W = 1, is reached for Q/(f_ρ/D) ~ 15 (cf. Eq. (29)). For f_ρ = 2, and D = 1, this corresponds to Q ~ 30. According to Eq. (28) this value is reached for Γ_e ~ 0.4, in agreement with our models with D = 1 in Table 2. For a large clumping factor D = 16, a value of W = 1 is already reached for Q ~ 2, suggesting an instability at Γ_e ~ 0.3, in agreement with our models with D = 16 in Table 2.

Fig. 9 HR diagram for pure He models (red), compared with our analytical relations. Long-dashed lines: analytical relation with fρ = 2, short-dashed line: analytical relation with fρ = 1.5.

In Fig. 9 we compare our analytical predictions, based on core radii R_c from Langer (1989b), with our numerical models. Long-dashed lines indicate predicted HR diagram positions within our analytical formalism for clumping factors D = 1, 4, and 16, and an adopted f_ρ = 2. The short-dashed line indicates the results for f_ρ = 1.5, and D = 1. Note that the use of the He-ZAMS models by Langer introduces inconsistencies. Nevertheless, our formalism describes the envelope inflation very well, in particular its dependence on the clumping factor D. Also the uncertainty due to the adopted f_ρ turns out to have only moderate influence on the results.

3.4. The mass of the surrounding shell ΔM

The mass ΔM of the dense shell surrounding the inflated zone (cf. Fig. 1) may be of interest, e.g. in relation to the energy budget, and the timescales of variations in ΔR.

To compute ΔM we use the following assumptions. 1) The radial extension of the shell is small compared to the total stellar radius, i.e., r = R_e; 2) the mass of the shell is small compared to the stellar mass, i.e., m = M; 3) the pressure P at the stellar surface is small compared to the pressure at the outer boundary of the inflated region P_e.

Under these assumptions, the equation of hydrostatic equilibrium (Eq. (1)) can be integrated directly, and we obtain (30)The mass of the shell thus simply follows from the fact that the pressure force at the rim of the low density region needs to be in equilibrium with the gravitational force on the shell. Again, P_e depends on material properties. Using P_e ≈ 0.4...0.9 × 10⁶   dyn/cm² from our models (cf. Tables 1 and 2), we obtain (31)ΔM thus mainly depends on the radius of the star. This is in agreement with our model computations (Tables 1 and 2) where we find ΔM as low as 10^-9   M_⊙ for compact WR stars, whereas our larger, LBV-type models reach ΔM up to 10^-2   M_⊙.

3.5. Dependence on model assumptions

The occurrence of an envelope inflation, involving a density inversion, has been reported in many previous works. Examples are the Wolf-Rayet models by Schaerer (1996), and models for very massive main sequence stars, as well as hydrogen-deficient carbon stars, by Saio et al. (1998). What all these stars have in common, is their proximity to the Eddington limit. Especially the latter work, which discusses the connection of inflated envelopes to the strange-mode instability, demonstrates the universality of the principles discussed here, for stars covering the range of 0.9–110 M_⊙.

On the other hand, the inflation effect has not been detected in many standard works on stellar evolution. E.g., in the rotating evolutionary models by Meynet & Maeder (2003), the 120   M_⊙ track enters the H-free WR stage with a mass of 23.2   M_⊙, and an effective core temperature of T_⋆ = 123 kK. According to our results in Table 1, we would expect a considerable envelope inflation, with T_⋆ ~ 60   kK, for the same object. In this section we discuss possible sources of such discrepancies, such as the influence of the imposed outer boundary condition (Sect. 3.5.1), dynamical terms on a secular timescale (Sect. 3.5.2), dynamical terms introduced by mass-loss (Sect. 3.5.3), and the treatment of convection and turbulence (Sect. 3.5.4).

3.5.1. The effect of the outer boundary condition

Fig. 10 The effect of ρ and T at the outer boundary, on the envelope solution. The solid blue line indicates our standard 23   M⊙ He-model in the Prad–Pgas plane. The dashed lines indicate solutions for which the boundary values for the density ρ, and the temperature T have been changed. The inflated zone between 2 and 6 R⊙ (as indicated by the radius marks) is mainly affected by the choice of the temperature T at the outer boundary. Temperatures T ≳ 80 kK, can lead to a cut-off of the inflated zone.

In Fig. 10 we illustrate the influence of the adopted boundary values for ρ, and T, on the envelope solution. The blue line indicates our standard 23   M_⊙ He-model, with boundary values T = 50 kK, and ρ = 5 × 10^-10   g   cm^-2. In addition we show six test models (dashed lines), for which we have changed ρ to values of 5 × 10^-8   g   cm^-2 and 5 × 10^-12   g   cm^-2, for temperatures of T = 50, 80, and 160 kK.

Firstly, the envelope solution turns out to be very robust with respect to changes of the outer boundary. For increasing P_rad, i.e., towards the interior of the star, all test models merge into the standard solution. However, boundary temperatures in the region of the Fe-opacity peak, i.e., temperatures that would lie within the inflated zone, lead to a cut-off of the inflated envelope. In Fig. 10, the inflated zone is indicated by radius marks (2–6   R_⊙). Models with outer boundary temperatures T ≳ 80 kK display envelopes where the inflated zone ends close to the respective outer boundary temperature, i.e. the radii are significantly reduced.

Boundary conditions that take the back-warming effect of an optically thick wind into account, as the ones discussed by Schaerer (1996), may thus inhibit an envelope inflation (Schaerer indeed detected a very similar effect in his model computations). In the same way, strong stellar winds with high temperatures at the wind base, as the ones described by Gräfener & Hamann (2005) for early WC subtypes, will prevent the formation of inflated envelopes. This is not surprising, as such winds are initiated by the hot Fe-opacity peak, i.e., the Fe-opacity peak lies within the wind. We thus expect a strong influence of the detailed wind physics on the inflation effect.

3.5.2. Dynamical terms on a secular timescale

The contraction/expansion of a star during its evolution can lead to energy sources/sinks, in addition to the nuclear energy production rate ϵ, which is included in our models (Eq. (3)). To estimate the importance of this effect, we compare the gravitational energy E_g, which is gained, or lost on a timescale τ, with the stellar luminosity L, and estimate in this way a limiting timescale (32)Here ΔM is the mass of the layers involved in the contraction/expansion, and R is a typical radius.

Because Wolf-Rayet stars are very compact, and evolve on relatively short timescales, such gravitational terms may play a role in their evolution. However, as the masses of the inflated envelopes discussed here are very small, these terms play no role for the inflation effect. For our 23.7 M_⊙ He-model, which represents the most extreme case of an inflated envelope, we estimate τ = 0.5 d, based on the parameters from Table 1, and R = R_c. We thus expect no influence of secular changes on the inflation effect. In Sect. 4.2 we will discuss the case of LBVs, for which gravitational terms indeed seem to play a role on timescales of the observed S Dor-type variability.

For the complete radiative envelope of a typical WR star (here we use the 15   M_⊙ model from Table 1, with ΔM = 1.5   M_⊙, and R = 0.4   R_⊙) we estimate τ = 5000 yr. Particularly in the end phases of the WR evolution, gravitational terms may thus play an important role in determining the L/M ratio of WR stars, which is an important input for our formalism.

3.5.3. Dynamical terms due to mass-loss

Because of the low densities, mass loss may introduce considerable velocity fields within the inflated envelopes, via the equation of continuity Ṁ = 4πρνr². The effect of the resulting dynamical terms on the envelope structure has been investigated by Petrovic et al. (2006). They found that strong mass loss can inhibit the formation of inflated envelopes, if the velocities within the inflated zone exceed the local escape speed, i.e., when ν(dν/dr) ≈ ν²/r has the same order of magnitude as GM/r². This condition imposes an upper limit for the mass-loss rate (33)where we have taken into account that the maximum velocity is reached close to R_m = (R_e + R_c)/2, where ρ ≈ ρ_min. An inspection of the dynamical term that would result from the density structure within our own models, suggests however that this limit should be lowered by  ~0.4 dex.

In Table 3 we have compiled conservative estimates of Ṁ₀, based on ρ_min = 5 × 10^-11   g/cm³, for a set of WR, and LBV models. For the WR models with the lowest temperatures, which may be subject to an envelope inflation, we find log (Ṁ₀/M_⊙yr^-1) =  −4.2... −3.7. For comparison, spectroscopically determined mass-loss rates by Hamann et al. (2006) reach values of up to log (Ṁ/M_⊙yr^-1) =  −4.3... −4.1, dependent on L. In view of the involved uncertainties, it is thus not clear whether the envelope inflation of WR stars is affected by dynamical effects.

Fig. 11 Herzsprung-Russel diagram of the Galactic WN stars, and the LBV AG Car. Red/blue symbols indicate observed HR diagram positions of WN stars from Hamann et al. (2006) (blue: with hydrogen (X > 0.05); red: hydrogen-free (X < 0.05)). Black symbols indicate the HR diagram positions of AG Car throughout its S Dor Cycle from 1985–2003, according to Groh et al. (2009b). Large symbols refer to stars with known distances from cluster/association membership. The symbol shapes indicate the spectral subtype (see inlet). Arrows indicate lower limits of T⋆ for stars with strong mass loss. The observations are compared to stellar structure models with/without hydrogen from this work (blue/red lines, cf. Fig. 6), and for AG Car (X = 0.36, black line).

3.5.4. Convection and turbulence

As already noted in Sects. 2.2, and 3.1, the inefficiency of convection is a mandatory condition for the occurrence of an envelope inflation. As the convective efficiency may be affected by details of the numerical implementation, or effects that are not taken into account in the standard mixing-length approach, we give basic physical arguments, why the convective efficiency should be low for the cases discussed in this work.

Generally, the strong increase in opacity near the iron bump will lead to the onset of convection through the usual Schwarzschild stability criterion. But for the models here such convection should be quite inefficient in the sense that it can carry only a small fraction of the local stellar flux F. To see this, note that a robust upper limit to the flux carried by convection is set by the free-streaming of internal energy U at the local sound speed ν_s, F_fs = ν_sU. In general, U = U_gas + U_rad, but in the luminous stars here, and particularly in the region around the iron bump, the radiative energy dominates¹. The associated maximum fraction of stellar flux F that could be carried by convection is thus (34)where T_c is the effective temperature for the (pre-inflated) core radius R_c. The last equality applies the characteristic iron bump temperature T ≈ 1.5 × 10⁵ K, with associated sound speed ν_s ≈ 50 km s^-1. This shows directly that for core effective temperatures appropriate for Wolf-Rayet stars, convection can carry only a small fraction of the total stellar flux, viz. 0.3% for T_c = 10⁵ K, and about 5% for T_c = 50   000 K. This implies that, even when convection sets in, the bulk of the stellar flux is still transported by radiative diffusion, thus justifying this basic simplifying assumption in our analysis of envelope inflation.

We note that, due to the dominance of radiation pressure, and radiation energy in the inflated zone, the contributions from internal gas pressure/energy are generally negligible with respect to F_fs. The same holds for turbulent pressure/energy, as P_turb = (2/3)u_turb = (1/3)ρa²ℳ² is of the same order of magnitude as gas pressure/energy (cf., Maeder 2009, p. 105). Here u_turb denotes the turbulent energy density, a the sonic speed, and ℳ the Mach number, which is expected to be lower than one.

For low convective efficiency, the structure of inflated envelopes is solely determined by Eq. (16), which is (in first order approximation) equivalent to the condition that Γ = 1. The latter only follows from material properties, so that the addition of turbulent pressure has almost no effect on the resulting envelope structure.

4. Implications

Let us next discuss the implications of the predicted envelope inflation for stars close to the Eddington limit, i.e., for WR stars and LBVs. Both types of objects are extremely luminous, and are thought to reside close to the Eddington limit. In Sect. 4.1, we focus on the still unexplained WR locations in the HR diagram. In Sect. 4.2, we discuss the possible connection of the radius inflation to the S Doradus radius changes characteristic of LBVs.

4.1. The radius problem for Wolf-Rayet stars

In Fig. 11, we show an HR diagram for galactic WR stars, as obtained from spectral analyses by Hamann et al. (2006). As already highlighted by these authors there exists a radius problem for H-free WR stars, namely that these objects are not located on the theoretical He-main sequence. According to a population synthesis performed by Hamann et al., almost all WR stars with a hydrogen surface abundance X < 0.05 are expected to display stellar temperatures T_⋆ in excess of 100 kK. However, in the observed HR diagram they cover a broad temperature range from 40–140 kK, with the majority below 100 kK (red symbols in Fig. 11). The observed WR radii are thus larger than predicted by stellar structure models (e.g., Langer 1989b).

At this point we note that the stellar temperatures T_⋆, as determined by Hamann et al. (2006), are not classical effective temperatures (related to τ = 2/3), but effective temperatures, related to the (hydrostatic) surface radius R_⋆, i.e., (35)with the Stefan-Boltzmann constant σ_SB. R _⋆  is the inner boundary radius of the employed atmosphere models, and is typically located at a much higher optical depth of τ_max ≈ 20. As long as R _⋆  is located in the hydrostatic part of the wind, its value does not depend on the detailed choice of τ_max². The T _⋆  determined by Hamann et al. thus reflect actual surface radii, which are not affected by the optical depth of the strong stellar winds of WR stars. The low T_⋆ in Fig. 11 thus display a true inconsistency with respect to the stellar structure models.

This is particularly surprising, as H-free WR stars are very simple objects. Due to their large convective cores, and high mass-loss rates they are expected to be close to chemical homogeneity (Langer 1989a). Even if the material in their He-burning cores is partly processed to carbon and oxygen, the mean molecular weight stays almost constant throughout the star. Moreover, the strong temperature sensitivity of He-burning guarantees an almost fixed temperature in their centers. The main unknown for the determination of the stellar structure is thus the opacity. In fact, with the advent of the OPAL opacity data (Iglesias & Rogers 1996) and their strong Fe-opacity peaks, Ishii et al. (1999) were the first to predict the large envelope inflation that we also describe in this work.

Ishii et al. (1999) could however not explain the HR diagram positions of these stars. In the galaxy, most H-free WR stars are located at luminosities below log (L/L_⊙) ≈ 5.8, reaching down to log (L/L_⊙) ≈ 5.3 (cf. Fig. 11). Ishii et al. demonstrated that unrealistically high metallicities, up to Z = 0.1, would be needed to achieve a significant radius inflation in this regime. An important result of the present work is that the same effect can be achieved by the assumption of an inhomogeneous (clumped) structure within the inflated zone (cf. Sect.2.4). For moderate clumping factors, up to D = 16, our models cover the observed parameter range (red (dashed) lines in Figs. 6 and 11). Notably, similar clumping factors are determined spectroscopically, for the winds of WR stars (e.g., Hamann & Koesterke 1998).

The metallicity dependence, as predicted by Ishii et al. (1999) and Petrovic et al. (2006), may however still play an important role in the determination of the temperatures of WR stars, as their subtype distribution seems to depend on Z (e.g., Crowther 2007). Later spectral subtypes are typically found at higher Z, which is in line with a stronger radius inflation due to the stronger Fe-opacity peak at higher Z. This effect may however also depend on the increased mass loss at higher Z (see Crowther & Hadfield 2006). We further note that, as discussed in Sect. 3.5.3, the mass-loss limit by Petrovic et al. (2006) may inhibit an envelope inflation for WR stars with strong mass-loss.

4.2. Radius changes of S Doradus-type LBVs

A defining property of LBVs are S Doradus-type radius variations on timescales of months to years. The inflation effect provides a reasonable explanation for radius variations on such timescales, as only small portions of the stellar envelope are affected. The underlying reason for the time dependence may be manifold. E.g., Georgiev et al. (2011) recently proposed that the variability of HD 5980, an LBV-like star that temporarily shows a WR-type spectrum, may be triggered by a transition from a hot, to a cool wind base. This picture fits very nicely into the framework of the inflation effect. As discussed in Sect. 3.5.1, we expect that a hot, optically thick wind would inhibit a radial inflation, i.e., a wind transition as proposed by Georgiev et al. could indeed lead to substantial radius changes.

In the following we discuss the case of AG Car, one of the best-studied galactic S Doradus LBVs (e.g., Leitherer et al. 1994; Stahl et al. 2001; Vink & de Koter 2002). More recently, AG Car has been studied in detail throughout its full S Dor cycle (1985–2003), by Groh et al. (2009b, 2011). During this period the star increased its radius from 58.5   R_⊙ to 120.4   R_⊙, corresponding to a change in the stellar temperature from T_⋆ = 26   450 K to 16 650 K. Such radius and temperature changes are characteristic aspects of the S Dor cycle, but they remain as yet unexplained by stellar structure calculations (see however Stothers & Chin 1993). For AG Car the luminosity was found to decrease from log (L/L_⊙) = 6.17 to 6.0 between the S Dor visual minimum and maximum phase. Furthermore, Groh et al. (2009b) determined a considerable hydrogen deficiency for AG Car, with a hydrogen surface mass fraction of X = 0.36.

We have computed a series of homogeneous stellar structure models for the relevant parameter range, with a hydrogen mass fraction of X = 0.36, and masses ranging from 60–73   M_⊙ (cf. Table 2). These models display a similar luminosity range as AG Car, and partly show very low T_⋆ , down to  ~17 000 K, due to a substantial radius inflation. A comparison with the observed HR diagram positions of AG Car is shown in Fig. 11. Notably, during the minimum phases of the S Dor cycle, AG Car (indicated by black symbols) shows a very good agreement with our homogeneous models (black dashed line), i.e., its position in the HR diagram can be explained by a chemically homogeneous star with M ~ 70   M_⊙, and a substantial radius inflation.

The core of such a star would be relatively small, with R_c ~ 14   R_⊙, corresponding to an effective core temperature T_c ~ 53 kK (cf. Table 2). Due to the envelope inflation, the star would display a core-halo structure, as described in Sect. 2.2. According to Eq. (31), the mass ΔM, of the outer shell around the inflated envelope, is expected to change throughout the S Dor cycle. For the smallest radius (58.5   R_⊙), we obtain ΔM ~ 10^-4   M_⊙, while ΔM ~ 2 × 10^-3   M_⊙ for the largest extension (120.4   R_⊙). The energy to lift this material from the core radius to the stellar surface is ΔE_g ~ ΔMMG/R_c = 3.8 × 10⁴⁶   erg. This value compares well with the “missing energy” due to the decrease in luminosity ΔE_lum = 1.83 × 10³⁹   erg/s × 575d = 9.1 × 10⁴⁶   erg (cf. Groh et al. 2009b). We thus find that within the inflated envelope, the S Dor-type variability introduces dynamical effects as the ones discussed in Sect. 3.5.2.

The good agreement with our models suggests that the extension of the envelope of AG Car may be explained by the inflation effect. The mass-loss rate of AG Car (log (Ṁ/M_⊙yr^-1) <  −4.2, Groh et al. 2009b) lies well below the upper mass-loss limit discussed in Sect. 3.5.3 (cf. Table 3), implying a stable, static configuration that is not affected by dynamical effects due to mass loss. Based on our chemically homogeneous models, AG Car lies close to the Eddington limit, with an Eddington factor of Γ_e = 0.39. According to Eq. (28), this results in a Q parameter of  ~9, and with Eq. (29) we obtain W_c ~ (9/f_ρ)(R_c/R_⊙)/(M/M_⊙) ~ 0.8 (where we have adopted f_ρ = 2.2 from Table 2). The star is thus very close to the instability limit (W_c = 1), as discussed in Sect. 3.2.

One possible mechanism that could lead to an unstable situation involving W_c > 1, is the reduction of the effective stellar mass through rotation. Groh et al. (2011) found that the observed rotational velocity of AG Car declines from ν_rotsini = 220 km s^-1 in the minimum phase, to 90 km s^-1 for the largest extension. For these values the ratio of centrifugal to gravitational acceleration at the stellar surface decreases from , to 0.07/sin²i. The corresponding reduction of the effective stellar mass by 20% in the minimum phase, could thus just lead to a situation where W_c ≈ 1. In this case the outer stellar envelope would become unstable, and could start expanding. The reduction of the rotational velocity during the expansion would lead to W_c < 1, and a re-stabilization of the envelope.

Groh et al. (2009a, 2011) proposed a somewhat similar scenario for AG Car, where the Eddington limit is exceeded due to rotation. In this work, we have demonstrated that already for moderate Eddington factors of Γ_e ~ 0.4³ AG Car may form an inflated stellar envelope that is dominated by radiation. Within this envelope the star adjusts very precisely to the Eddington limit (Γ = 1). In our scenario, the instability that potentially induces the S Dor cycle is due to the fact that the specific internal energy due to radiation exceeds the gravitational potential, with the stellar envelope thus becoming unbound.

We want to point out that the formation of inflated envelopes does not depend on the assumption of chemical homogeneity in our numerical models. The instability follows from our formalism in Sect. 3, and is thus independent of the internal structure of the star. For AG Car, the observed luminosity (log (L/L_⊙) = 6.17), and surface hydrogen mass fraction (X = 0.36) imply a mass of 74   M_⊙ for the case of chemical homogeneity (cf. Gräfener et al. 2011). In reality, AG Car is likely more chemically enriched in the core, i.e., its mass is lower, and its core radius larger. This would lead to an even earlier onset of the inflation effect, and its instability. As the evolutionary status of LBVs is not well known, it is however not possible to obtain reliable mass estimates for this type of object.

5. Conclusions

In the present work, we have discussed the possible formation of radially inflated stellar envelopes for stars approaching the Eddington limit, as originally proposed by Ishii et al. (1999), and Petrovic et al. (2006). In addition to numerical models, we could provide an analytic description of this process, leading to the discovery of a new instability limit, and a clumping dependence. Both effects turn out to be of profound importance. While the former may be connected to S Doradus-type instabilities in LBVs, the latter can account for the large observed radii of Galactic WR stars, and thus resolve the long standing WR radius problem.

Within our analytical approach we could describe the inflation effect in terms of a dimensionless parameter W (Eqs. (28) and (29)) that characterizes the ratio between internal and gravitational energy at the base of the envelope. W can be estimated on the basis of opacity tables and a combination of stellar parameters. The envelope inflation occurs when W approaches one. For W ≥ 1, we find that the envelope becomes gravitationally unbound, i.e., there exists no static solution. A prerequisite for envelope inflation is the proximity of the star to the Eddington limit, with Eddington factors of Γ_e > 0.3–0.35, at solar metallicity. For a given chemical composition and Γ_e, the resulting stellar radius depends on the ratio M/R, and on the degree of (in)homogeneity of the material within the envelope (characterized by a clumping factor D).

Envelope inflation has a strong impact on the effective temperatures of stars in the upper HR diagram. For Wolf-Rayet stars it can account for the “radius problem”, i.e., that the observed stellar temperatures of H-free WR stars are much lower than predicted by canonical stellar structure models. To reproduce the observed HR diagram positions of these objects, it is necessary to assume that the material within the envelope is clumped. The resulting clumping factors lie in the range of D = 1–16, in notable agreement with typical clumping factors detected in the winds of WR stars (e.g., Hamann & Koesterke 1998). As the envelope inflation is expected to depend on metallicity (Ishii et al. 1999; Petrovic et al. 2006), it may also account for the observed Z-dependence of the spectral subtype distribution of WR stars, i.e., that WR stars in high-Z environments show lower effective temperatures (e.g., Crowther 2007).

For the LBV AG Car we suggest that the observed HR diagram position of a  ~70   M_⊙ star, with a relatively compact core, might be subject to substantial envelope inflation. During its S Dor cycle, AG Car is found to increase its radius by a factor of two, while its luminosity decreases by a factor of 1.5. The luminosity decrease might be explained as a result of the formation of a dense outer shell with a mass of  ~2 × 10^-3   M_⊙ as predicted by our models. The energy needed to lift this material from the core to the outer radius is of the same order of magnitude as the “missing energy” due to the decrease in luminosity. Similar to the mechanism proposed by Groh et al. (2009a, 2011), the S Dor variability of AG Car may be due to the effect of stellar rotation, which reduces the effective stellar mass, and thus leads to an instability with W > 1, that forces the star to expand.

We conclude that the stellar effective temperatures in the upper HR diagram are potentially strongly affected by the inflation effect. The peculiar structure of the inflated envelopes,

with an almost void, radiatively dominated region, beneath a thin and dense shell could mean that many, rather very compact stars, are hidden below inflated envelopes, and thus display much lower effective temperatures to the observer. Also compact, fast rotating stars may be hidden this way. This may particularly affect massive stars just before the final collapse, and the question of whether they form a SN Ibc, or a GRB (for the latter smaller core radii are inferred than the observed WR radii, cf. Cui et al. 2010). The complex envelope structure may also affect the early X-ray afterglow of GRBs (e.g. Li 2007).

Acknowledgments

We thank A. Maeder for a very constructive referee report, and S.-C. Yoon for helpful discussions about evolutionary effects in WR stars.

References

All Tables

All Figures

	Fig. 1 Density vs. radius for our 23   M⊙ He model (cf. Table 1).
In the text

	Fig. 2 Pressure vs. radius. Gas pressure Pgas, radiation pressure Prad, and total pressure Ptot = Pgas + Prad.
In the text

	Fig. 3 Total stellar luminosity L(r), compared to the Eddington luminosity LEdd = 4πcGM/κ.
In the text

	Fig. 4 Stellar structure in the Prad–Pgas plane. Gas pressure Pgas, radiation pressure Prad, and total pressure Ptot = Pgas + Prad are plotted vs. Prad throughout the whole star.
In the text

Fig. 5 Envelope solution in the Prad–Pgas plane. The colors indicate the logarithm of the Eddington factor Γ, for given Prad and Pgas, according to the OPAL opacity tables (Iglesias & Rogers 1996), for our pure He model with 23   M⊙. The envelope solution almost precisely follows a path with Γ = 1. The radii of 2, 4, and 6   R⊙ within the inflated envelope are indicated. Arrows indicate the slope of the solution, as expected from Eq. (12).

In the text

	Fig. 6 HR diagram for pure helium models (red), and models containing hydrogen (blue).
In the text

	Fig. 7 Plot of 1/ρ vs. Prad. The dashed horizontal lines indicate 1/ρmin, and 1/(10   ρmin). The intersections of the latter with the envelope solution represent boundary pressure Pc, and Pe, as indicated by vertical dashed lines. The grey shaded area illustrates that the integral ∫   dP/ρ ≈ ΔP/(2ρmin).
In the text

	Fig. 8 Fitting relation for Q(X,Z = Z⊙,Γe). Black crosses indicate numerically determined values of Q, for X = 0 – 0.7. Blue lines indicate our fitting relation Eq. (28), for X = 0 (top), and X = 0.7 (bottom).
In the text

	Fig. 9 HR diagram for pure He models (red), compared with our analytical relations. Long-dashed lines: analytical relation with fρ = 2, short-dashed line: analytical relation with fρ = 1.5.
In the text

Fig. 10 The effect of ρ and T at the outer boundary, on the envelope solution. The solid blue line indicates our standard 23   M⊙ He-model in the Prad–Pgas plane. The dashed lines indicate solutions for which the boundary values for the density ρ, and the temperature T have been changed. The inflated zone between 2 and 6 R⊙ (as indicated by the radius marks) is mainly affected by the choice of the temperature T at the outer boundary. Temperatures T ≳ 80 kK, can lead to a cut-off of the inflated zone.

In the text

Fig. 11 Herzsprung-Russel diagram of the Galactic WN stars, and the LBV AG Car. Red/blue symbols indicate observed HR diagram positions of WN stars from Hamann et al. (2006) (blue: with hydrogen (X > 0.05); red: hydrogen-free (X < 0.05)). Black symbols indicate the HR diagram positions of AG Car throughout its S Dor Cycle from 1985–2003, according to Groh et al. (2009b). Large symbols refer to stars with known distances from cluster/association membership. The symbol shapes indicate the spectral subtype (see inlet). Arrows indicate lower limits of T⋆ for stars with strong mass loss. The observations are compared to stellar structure models with/without hydrogen from this work (blue/red lines, cf. Fig. 6), and for AG Car (X = 0.36, black line).

In the text