Wind modelling of very massive stars up to 300 solar masses

# Wind modelling of very massive stars up to 300 solar masses

###### Key Words.:
Stars: early-type – Stars: mass-loss – Stars: winds, outflows – Stars: evolution
1

The stellar upper-mass limit is highly uncertain. Some studies have claimed there is a universal upper limit of 150. A factor that is often overlooked is that there might be a significant difference between the present-day and the initial masses of the most massive stars – as a result of mass loss. The upper-mass limit may easily supersede 200. For these reasons, we present new mass-loss predictions from Monte Carlo radiative transfer models for very massive stars (VMS) in the mass range 40-300 , and with very high luminosities 6.0 7.03, corresponding to large Eddington factors . Using our new dynamical approach, we find an upturn or “kink” in the mass-loss versus dependence, at the point where the model winds become optically thick. This coincides with the location where our wind efficiency numbers surpass the single-scattering limit of , reaching values up to . In all, our modelling suggests a transition from common O-type winds to Wolf-Rayet characteristics at the point where the winds become optically thick. This transitional behaviour is also revealed with respect to the wind acceleration parameter, , which starts at values below 1 for the optically thin O-stars, and naturally reaches values as high as 1.5-2 for the optically thick Wolf-Rayet models. An additional finding concerns the transition in spectral morphology of the Of and WN characteristic He ii line at 4686Å. When we express our mass-loss predictions as a function of the electron scattering Eddington factor alone, we obtain an vs. dependence that is consistent with a previously reported power law (Vink 2006) that was based on our previous semi-empirical modelling approach. When we express  in terms of both  and stellar mass, we find optically thin winds and   for the range 0.4 0.7, and mass-loss rates that agree with the standard Vink et al. recipe for normal O stars. For higher values, the winds are optically thick and, as pointed out, the dependence is much steeper,   . Finally, we confirm that the effect of on the predicted mass-loss rates is much stronger than for the increased helium abundance (cf. Vink & de Koter 2002 for Luminous Blue Variables), calling for a fundamental revision in the way stellar mass loss is incorporated in evolutionary models for the most massive stars.

## 1 Introduction

The prime aim of this paper is to investigate the mass-loss behaviour of very massive stars (VMS) with masses up to 300 that approach the Eddington limit. Mass loss from hot massive stars is driven by radiative forces on spectral lines (Lucy & Solomon 1970; Castor, Abbott & Klein 1975; CAK). CAK developed the so-called force multiplier formalism in order to treat all relevant ionic transitions. This enabled them to simultaneously predict the wind mass-loss rate, , and terminal velocity, , of O-type stars. Although these predictions provided reasonable agreement with observations, they could account for neither the high wind efficiencies of the denser Of stars with their strong He ii 4686Å lines, nor that of the even more extreme Wolf-Rayet (WR) stars.

This discrepancy had been proposed as due to the neglect of multi-line scattering (Lamers & Leitherer 1993, Puls et al. 1996). Using a global energy Monte Carlo approach (Abbott & Lucy 1985, de Koter et al. 1997) in which the velocity law was adopted – aided by empirical constraints – Abbott & Lucy (1985) and Vink et al. (2000) provided mass-loss predictions for galactic O stars including multi-line scattering. This appeared to solve the wind momentum problem for the denser O-star winds. Mass-loss rates were obtained that were a factor 3 higher than for cases in which single scattering was strictly enforced.

Historically, the situation for the WR stars was even more extreme. Here, values of 10 had been reported (e.g. Barlow et al. 1981). With the identification of major wind-clumping effects on the empirical mass-loss rates (Hillier 1991; Moffat et al. 1994; Hamann & Koesterke 1998), these numbers should probably be down-revised to values of 3. Although it has been argued that WR winds are also driven by radiation pressure (Lucy & Abbott 1993, Springmann 1994, Gayley & Owocki 1995, Nugis & Lamers 2002, Gräfener & Hamann 2005), the prevailing notion is still that these optically thick outflows of WR stars, where the sonic point of the accelerating flow lies within the pseudo or false photosphere, are fundamentally different from the transparent line-driven O-star winds (e.g. Gräfener & Hamann 2008).

For O-type stars Müller & Vink (2008) recently suggested a new parametrization of the line acceleration, expressing it as a function of radius rather than of the velocity gradient (as in CAK theory). The implementation of this new formalism improves the local dynamical consistency of Monte Carlo models that originally imposed a velocity law. Not only do we find fairly good agreement with observed terminal velocities (see also Muijres et al. 2011b), but as our method naturally accounts for multi-line scattering, it is also applicable to denser winds, such as those of WR stars.

Still adopting a velocity stratification, Vink & de Koter (2002) and Smith et al. (2004) predicted mass-loss rates for Luminous Blue Variables (LBVs), and showed that is a strong function of the Eddington factor for these objects. They also found that, despite their extremely large radii, even LBV winds may develop pseudo-photospheres under special circumstances: when they find themselves in close proximity to the bi-stability and Eddington limits.

In this paper, our aim is to study the mass-loss behaviour of stars as they approach the Eddington limit (see also Gräfener & Hamann 2008). We do this in a systematic way by targeting VMS in the range 40-300 . A pilot study was performed by Vink (2006) who found a steep dependence of on , finding , but this was obtained using the earlier global energy approach, in which the velocity stratification was adopted, rather than our new dynamically-consistent approach that we explore in the following.

Both approaches have their pros and cons. In the semi-empirical approach, the terminal velocity constraint (where / is constant) might aid the modelling at imposing the correct wind structure, and as long as the adopted velocity law is close to the correct one, it provides the most accurate mass-loss rates. The prime advantage of this approach is that any missing physics that could affect the mass-loss rate might be balanced by employing the empirical (and probably close to realistic) terminal-wind velocity. The second approach, however, is theoretically more appealing, as it enforces local dynamical consistency, and one thus no longer needs to rely on free parameters. Ultimately, one would aspire to predict accurate mass-loss rates and terminal velocities simultaneously from first principles via the second approach explored here.

The stellar upper-mass limit is highly controversial. On purely statistical grounds, some investigators have claimed there is a universal upper limit of 150 (e.g. Weidner & Kroupa 2004, Oey & Clarke 2005, Figer 2005). However, a physical factor that is often overlooked concerns the possibility of a significant difference between the present-day and the initial masses of the most massive stars, as a result of strong mass loss. In other words, the initial masses of the most massive stars may be significantly above 150, possibly superseding 200  (e.g. Figer et al. 1998, Crowther et al. 2010, Bestenlehner et al. 2011). The issue of the upper mass-limit will remain uncertain as long as there is only limited quantitative knowledge of mass loss in close proximity to the Eddington limit.

Our aim is thus to explore wind models of stars with masses up to 300, using a well-established methodology that has been extensively tested against observations for lower mass common O-type stars. VMS have been proposed as leading to the production of intermediate mass (of the order of 100) black holes that have been suggested to be at the heart of ultra-luminous X-ray sources (Belkus et al. 2007 and Yungelson et al. 2008). Clearly, the success of such theories depends critically on the applied mass-loss rates. The present study may help advance these theories.

Our paper is organized as follows. In Sect. 2 we briefly describe the Monte Carlo mass-loss models, before presenting the parameter space considered in this study (Sect. 3). The mass-loss predictions (Sect. 4) are followed by a description of the spectral morphology of the Of-WN transition in terms of the characteristic He ii 4686Å line in Sect. 5. Subsequently, we compare our wind model parameters against empirical values for the most massive stars in the Arches cluster (Sect. 6.1), as well as theoretical models (Sect. 6.2) of CAK and Gräfener & Hamann (2008), before ending with a summary in Sect. 7.

## 2 Monte Carlo models

Mass-loss rates are calculated with a Monte Carlo method that follows the fate of a large number of photon packets from below the stellar photosphere throughout the wind. The core of our approach is related to the total loss of radiative energy that is coupled to the momentum gain of the outflowing material. Since the absorptions and scatterings of photons in the wind depend on the density in the wind, hence on the mass-loss rate, one can obtain a consistent model where the momentum of the wind material equals the transferred radiative momentum. We have recently improved our dynamical approach (Müller & Vink 2008, Muijres et al. 2011b) and are now able to predict simultaneously with and the wind structure parameter . The essential ingredients and assumptions of our approach have been discussed more extensively in Abbott & Lucy (1985), de Koter et al. (1997), and Vink et al. (1999). Here we provide a brief summary.

The Monte Carlo code mc-wind uses the density and temperature stratification from a prior model atmosphere calculation performed with isa-wind (de Koter et al. 1993, 1997). These model atmospheres account for a continuity between the photosphere and the stellar wind, and describe the radiative transfer in spectral lines by adopting an improved Sobolev treatment. The chemical species that are explicitly calculated (in non-LTE) are H, He, C, N, O, S, and Si. The iron-group elements, which are crucial for the radiative driving and the calculations, are treated in a generalized version of the “modified nebular approximation” (e.g. Schmutz 1991). However, we performed a number of test calculations in which we treated Fe explicitly in non-LTE. These tests showed that differences with respect to the assumption of the modified nebular approximation for Fe were small. Therefore, we decided to treat Fe in the approximate way, as was done in our previous studies.

The line list used for the MC calculations consists of over of the strongest transitions of the elements H - Zn extracted from the line list constructed by Kurucz & Bell (1995). The wind was divided into 90 concentric shells, with many narrow shells in the subsonic region, and wider shells in supersonic layers. For each set of model parameters, a certain number of photon packets was followed, typically .

Other assumptions in our modelling involve wind stationarity and spherical geometry. The latter seems to be a good approximation, as the vast majority of O-type stars show little evidence of significant amounts of linear polarization (Harries et al. 2002, Vink et al. 2009). Nevertheless, asphericity has been found in roughly half the population of LBVs (Davies et al. 2005, 2007), although those polarimetry results have been interpreted as the result of small-scale structure or “clumping” of the wind, rather than of significant wind asymmetry.

With respect to wind clumping, it has been well-established that small-scale clumping of the outflowing gas can have a pronounced effect on the ionization structure of both O-star and Wolf-Rayet atmospheres (e.g. Hillier 1991). This has lead to a downward adjustment of empirical mass-loss rates, by factors of up to three (e.g. Moffat et al. 1994, Hamann & Koesterke 1998, Mokiem et al. 2007, Puls et al. 2008), and possibly even more (Bouret et al. 2003; Fullerton et al. 2006). In addition, clumping may have a direct effect on the radiative driving, therefore on predicted mass-loss rates (e.g. Gräfener & Hamann 2008). The subtle issues of both clumping and porosity on the predicted mass-loss rates have recently been investigated by Muijres et al. (2011a). Whilst it was found that the impact on  can be high for certain clumping prescriptions, the overall conclusion was that clumping does not affect the wind properties of O-type dwarfs and supergiants in a dramatic way for moderate clumping factors and porosity. Stars close to the Eddington limit, however, may be much more susceptible to – even modest – clumping (Shaviv, 1998, 2000; van Marle et al. 2008). Another relevant factor that might affect the results concerns the interplay of radiation with rotation, as the Eddington factor for a rotating star depends explicitly on the rotation rate (Langer 1997, Maeder & Meynet 2000). As a consequence, the effective Eddington limit could be reduced by rotation, and might even become the dominant factor.

In the present set of computations, we do not account for the effects of wind clumping or rotation, but it should be kept in mind that these effects might play a quantitative role.

## 3 Parameter space and model applicability

Stars approach the Eddington limit when gravity is counterbalanced by the radiative forces; i.e. . Photons can exert radiative pressure through bound-free, free-free, electron scattering, and bound-bound interactions. In early-type stars hydrogen, the dominant supplier of free electrons, is fully ionized. Therefore is essentially independent of distance and constitutes a fixed number for each model. Because of this useful property, which provides a well-defined and simple quantitative handle, we opt to discuss our results in terms of . We discuss this choice in more detail in Sect. 3.1

The dependence of the mass-loss rate on  represents a non-trivial matter because depends on both the mass and the stellar luminosity . To properly investigate the effect of high on mass-loss predictions, we first need to establish the relevant part of parameter space in terms of , , and . For a fully ionized2 plasma,  equals

 Γe≡gegNewton=L⋆σe4πcGM⋆=10−4.813(1+X)(%$L⋆$L⊙)(M⋆%$M⊙$)−1 , (1)

The luminosities are chosen in such a way that in combination with the stellar mass the desired  value is obtained. It is principally the effective temperature that sets the ionization stratification in the atmosphere, and thus determines which lines are most active in driving the wind. As a result,  affects the predicted mass-loss rate. For most parts of this paper, we investigate the influence of  for a fixed stellar temperature of 50 000 K. The  dependence is studied separately in Sect. 4.4. All models are for the solar metallicity from Anders & Grevesse (1989) and with the element-to-element distribution from Allen (1973). The prime reason we use these older abundances rather than the newer (and lower) solar or B-type star abundances is to be able to directly compare the new results to the older Vink et al. (2000) rates. We note however that it is the element of Fe that sets the mass-loss rate, and this element has not changed. The low (CNO) and intermediate-mass elements however dominate the outer wind, where the terminal wind velocity is set. Nonetheless, even a substantial decrease in these abundances (by several tens of percent) is not expected to lower the terminal wind velocities significantly, since the terminal wind velocity has been found to depend only weakly on metallicity (Leitherer et al. 1992).

We divided our model stars into three different groups according to their characteristics. The first group comprises objects that have relatively common O-star masses in the range 40-70 . The second group of objects are rather high-mass stars within the “observable range” of 70-120 . They might be close to the Eddington limit already early-on on the main sequence because of their intrinsically high luminosity. The third group involves very massive stars in the mass range 120-300 . They are near the Eddington limit for the same reason as the second group. So far, there is a lack of compelling observational evidence of any such stars in the present-day universe; however, we note that Crowther et al. (2010) have suggested a revision of the upper-mass limit to 300 .

The bulk of the models in our grid have been chosen such that the behaviour of mass loss as a function of and can be studied separately. The grid is presented in Table 1. We note that the combinations are intentionally rather extreme to assure high  values. The reason is to specifically map that part of parameter space where physically the most extreme winds are expected to appear.

### 3.1 Model applicability regime

With respect to the potential limitations of our modelling approach, we make one rather stringent assumption in the manner the (sub-) photospheric density structure is set-up. In the deepest layers of the model atmosphere (with 1km/s), we assume that the run of density is provided by the equation of motion using , so we apply . In reality , as well as being depth-dependent as a result of bound-bound, bound-free, and free-free processes. Notably, the opacities from millions of weak iron lines may contribute significantly, but they are largely neglected in the deep layers of our models.

Nugis & Lamers (2002) highlight the importance of the iron peak opacities in deep photospheric layers for the initiation of Wolf-Rayet winds (see also Heger & Langer 1996). This approach was subsequently included in models by Gräfener & Hamann (2005, 2008) for WC and WNL stars. They find that the presence of these opacity bumps may locally cause to approach unity, leading to the formation of optically thick winds. In our Monte Carlo approach, we trace the radiative driving of the entire wind, and as most of the energy is transferred in the supersonic part of the outflow, we are less susceptible to the details of the (sub)photospheric region. However, this also means that we do not treat these deep regions self-consistently. This implies that we can (and we will) compute model atmospheres with values of very close to one. This strategy has the advantage of allowing us to explore the transition from transparent to dense stellar winds. As our models do capture the full physics in the layers around and above the sonic point, we argue that they correctly predict the qualitative behaviour of dense winds, but that for one of our optically thick wind models would correspond to a model with smaller if the ionic contributions were included in the deepest parts of the atmosphere. This “shift” in is not fixed but would depend on the sonic point temperature and density. From the behaviour of the Rosseland mean opacity, we would expect the size of the shift to increase at higher and higher temperatures.

If exceeds unity at some depth in the subphotospheric part of the atmosphere, a density inversion is expected to occur for the static case, i.e. for increasing radial distance from the centre, the density very near the domain where 1 is anticipated to increase. This is encountered in studies of stellar structure and evolution, but it is unclear what really happens in nature. The potential effects may involve strange-mode pulsations (e.g. Glatzel & Kiriakidis 1993), subsurface convection (Cantiello et al. 2009), or an inflation of the outer stellar envelope (e.g. Ishii et al. 1999). These processes tend to occur only when is very close to unity or above it (see e.g. Petrovic et al. 2006). In assessing the outcome of our computations, we find that at the results behave rather oddly. Though we show the wind results over the entire  range, we only quantify the mass-loss rates up to this value of . This boundary is indicated by a vertical dashed line in all relevant figures.

## 4 Results

### 4.1 Mass-loss predictions at high Γe

Table 1 lists our mass-loss predictions for all three considered mass ranges. Most columns are self-explanatory, but we note that the effective escape velocity (7th column) is defined as . The predicted wind terminal velocities, mass-loss rates, wind efficiency numbers, and wind acceleration parameter are given in columns (8), (9), (10), and (11), respectively. For comparison column (12) lists the mass-loss values from the standard mass-loss recipe of Vink et al. (2000) where was held fixed at unity.

The predicted mass-loss rates (column 9) are shown in Fig. 1. Different symbols are used to identify the different mass ranges. The figure shows that increases with . This is in qualitative agreement with the luminosity dependence of the standard mass-loss recipe of Vink et al. (2000), derived from a set of models with   0.4. Analogous to the results from the standard Vink et al. (2000) recipe, Fig. 1 suggests that there is an additional mass-loss dependence on mass, as for fixed   the higher mass stars have higher mass-loss rates. This finding confirms that mass-loss rates cannot solely be described by a dependence on luminosity or Eddington factor. This will be discussed further in Sect. 4.2.

When comparing columns (9) and (12) from Table 1, it can be noted that our new high mass-loss predictions tend to be larger than those determined using the standard Vink et al. (2000) recipe. In order to quantify these differences, we divide the new mass-loss rates over those determined using the Vink et al. (2000) recipe (using the derived terminal wind velocities as input), and show the results in Fig. 2. For the range 0.7, the differences are small. However, for values of exceeding 0.7, the new and the old results diverge sharply. The maximum difference reaches a factor of five, which is similar in magnitude to what was reported previously for LBVs (Vink & de Koter 2002) and WR stars (Vink & de Koter 2005). We note that, although these prior results were based on global energy consistency with fixed /, where the velocity stratification was adopted, the reason for the differences revealed in Fig. 2 is that we probe a different part of parameter space.

We now turn our attention to the wind velocity structure. We first inspect the associated terminal wind velocity predictions. Figure 3 shows the behaviour of terminal wind velocity versus . The highest values are reached for the highest mass stars and exceed 5000 km/s. As expected drops with . In the range 0.4-0.95, the terminal wind velocity divided over the escape velocity is of the order 3-4, which is similar to the values for common O-type stars (Muijres et al. 2011b), where is in the range 30-40 kK, and the wind velocities are closer to 3000 km/s (see Sect. 4.4).

We next turn our attention to the other wind velocity structure parameter, , which describes how rapidly the wind accelerates. The predicted values of are depicted in Fig. 4, but does not show a significant dependence on stellar mass. For up to 0.7, values are near unity, in accordance with the dynamical consistent models of Pauldrach et al. (1986), Müller & Vink (2008), and Muijres et al. (2011b). However, when exceeds 0.7 and approaches unity, steadily rises to values of about 1.7. These higher values are supposedly more commensurate in Wolf-Rayet stars (see e.g. Ignace et al. 2003), and it is reassuring to find that our models naturally predict this transition, without the use of any free parameter.

In all, our results suggest a natural extension from O-type mass loss to more extreme WR behaviour for increasing . An upturn in the behaviour is found at 0.7. Inspection of our models reveals a change from optically thin to optically thick wind models at the position where we obtained the kink in the mass-loss versus relationship. Closer scrutiny of our model output also revealed that here the character of the Fe line driving changes. Whilst various ionization stages of Fe contribute at all models, we find that for the optically thin models at the low end, just one ionization state of Fe dominates the relevant part of the wind driving domain (from just below the sonic point to about half the terminal velocity). By contrast, for the optically thick models at the high end, two or more ionization stages of Fe contribute to the primary driving regime.

### 4.2 Γe dependence of mass loss

In order to determine the dependence of the mass-loss rate on , we could simply fit the datapoints to a power law:

 ˙M∝Γep (2)

Using the semi-empirical approach Vink (2006) found to be equal to 5, and our dynamically consistent results provide the same slope here. However, in order to also take the mass dependence into account we divide the mass-loss rates by , and show the results in Fig. 5. We fit the data with the following power-law

 ˙M∝M⋆qΓep (3)

Below   0.95, Fig. 5 shows two mass-loss regimes, divided by a boundary at   0.7. This is not only the point where the slope of the mass-loss versus relation changes, but also where the wind efficiency parameter surpasses the single scattering limit (see below). Upon further inspection of our models, we find that as long as   0.7 the winds are optically thin, implying that the sonic point of the outflowing material lies outside the photosphere, whilst the winds become optically thick – with the photosphere moving outside of the sonic point – for  values above 0.7.

We derive two independent mass-loss relationships for the two separate regimes.

For 0.4 0.7 we find

 ˙M∝M0.68Γe2.2

For 0.7 0.95 we determine that

 ˙M∝M0.78Γe4.77

We first note that we intentionally do not provide equations here, as we expect the absolute values of these mass-loss predictions for the high models to be underpredicted. In Sect. 6.1, we see that our predicted wind terminal velocities are a factor 2-4 higher than found empirically. If we had employed our older semi-empirical approach, and assumed empirical – i.e. a factor 2-4 lower – terminal wind velocities, we would have predicted higher mass-loss rates. Indeed, the pilot study of Vink (2006) that was based on the semi-empirical approach provided such higher mass-loss rates.

Secondly, we note that the exact value of the transition value is model-dependent. We have already discussed our potential model deficiencies in the deepest layers in Sect. 3.1, but we should note that several other factors might also play a role; in particular, higher mass-loss rates – as would be obtained from our semi-empirical approach – would shift the kink to lower values. We note that the latter could occur in the case the wind is clumped. As long as porosity effects are small, wind clumping is expected to increase the momentum transfer in our MC models (Muijres et al. 2011a). Although future dynamical-consistent modelling of clumped winds is required to test this, wind clumping could potentially increase our predicted mass-loss rates, and subsequently decrease the value of the kink.

The above mass-loss relationships can easily be transformed using Eq. (1), leading to the entirely analogous mass-loss relationships and Eq. (4.2) to . Interestingly, if one subsequently applies a mass-luminosity relationship for classical (He-rich) WR stars of Maeder & Meynet (1987) or for very massive H-rich stars such as that of Yungelson et al. (2008), with for both cases, it follows that . This appears to be in good accord with the radio mass-loss rate relation for classical WR stars with measured masses from binaries by Abbott et al. (1986). It also agrees with the versus stellar mass relationship of that has been applied in WR evolution models by Langer (1989).

### 4.3 Increased wind efficiency close to the Eddington limit?

In order to learn whether radiation-driven mass-loss rates continue to increase with increasing  or reach a maximum in  instead, it is insightful to consider the wind efficiency parameter . We show the predicted values of in Fig. 6. As the symbols denote different mass ranges, the small scatter on the datapoints shows that is not very sensitive to stellar mass. At values of 0.5 we find wind efficiency numbers of order 1, in accordance with standard Vink et al. (2000) models. However, when approaches unity, rises in a curved manner to values as high as 2.5. Such high values are more commensurate with Wolf-Rayet winds than with common O star winds, and these results thus confirm a natural extension from common O-type mass loss to more extreme WR behaviour. In Sect. 6.1, we find that our predicted wind-terminal velocities are higher than the empirical values. As we note that an overprediction of the wind velocity is likely offset by a mass-loss rate underprediction by a similar amount, we argue that the combination of these two quantities, i.e. their product constituting the parameter, might be less affected by model deficiencies than either of these quantities would be individually.

The maximum mass loss in our models up to  = 0.95 is log = -3.8. This is the mass-loss rate that is retrieved for the most extreme models in our grid. Owocki et al. (2004) have investigated the mass loss of stars that formally exceed their Eddington limit and show that the expected mass loss falls well below the values required to account for the mass that is lost during LBV giant eruptions, such as that of Carinae in the 1840s. Interestingly, they introduce a porosity-moderated continuum driven mass loss that might account for the huge mass-loss rates associated with LBV eruptions (which may be of the order of 1/yr).

### 4.4 Effect of Teff on high Γe models

To establish whether there is an additional temperature dependence on , we varied  over the range 50-30 kK for selected  models, with mass-loss predictions presented in Fig. 7 and terminal wind velocities shown in Fig. 8. As we wish to stay above the temperature of the bi-stability jump (which starts at  values below 27.5 kK; see Vink et al. 2000), we restrict our  range to a minimum value of 30 kK. We find that for all  values  is not a strong function of temperature. In terms of the terminal velocity dependence, Fig. 8 shows a rather steep dependence on temperature, with dropping by a factor two. This is merely a reflection of the escape velocity dropping by a similar factor of two over the temperature range under consideration.

### 4.5 Effect of the helium abundance on high Γe models

To establish the existence of a potential helium dependence on , we computed additional models across the entire  region, setting the hydrogen abundance to zero and increasing the helium abundance accordingly. The results are listed in Table 2 and shown in Fig. 9. The mass-loss rates are similar to those of H-rich models for objects with the same  (see Table 1). This is not too surprising given that the indirect effects of different continuum energy distributions for H-rich versus H-poor are rather subtle (Vink & de Koter 2002). However, when the He-rich results are plotted in Fig. 9 they lie above the H-rich models. For equal luminosity and  the masses of the He-rich models are lower since  is a function of the chemical composition through (see Eq. 1); is lower for He-rich models, therefore the mass must be lowered to keep  constant. Similar to the H-rich models, there appears to be an upturn in the mass loss vs.  dependence for models at about   0.7.

With respect to the terminal velocity and -dependence, we do not find any significant differences between H-rich and He-rich models (see Table 1 versus Table 2).

## 5 Spectral morphology: the characteristic He 4686 Angstrom line

In the previous section, we provided evidence for a natural transition in the mass-loss- exponent, as well as in the velocity parameter and wind-efficiency from moderate “optically thin wind” cases to “optically thick wind” cases for objects that find themselves above 0.7. We have inspected our models and confirmed that for 0.7 the sonic velocity is reached outside the photosphere, whilst the stars form a pseudo-photosphere for 0.7.

We expect that the occurrence of a pseudo-photosphere has a consequence for the spectral morphology of the stars in question. We might suspect that the transition 0.7 is the point where the spectral morphology of normal O stars changes from the common O and Of-types into a WN-type spectrum. The spectral sequence involving the Of/WN stars has a long history (e.g. Conti 1976, Walborn et al. 1992, de Koter et al. 1997, Crowther & Dessart 1998) but it still has to be placed into a theoretical context. Figure 10 shows a sequence for the predicted He ii 4686Å lines for three gradually increasing values of : 0.70 (model 24), 0.84 (model 15), and 0.93 (model 20), respectively. These models have been selected to be objects with a constant luminosity of , and we simply lowered the mass from 120  to 100, to 90. It is insightful to note that, although the first spectrum below the transition already shows some emission – characteristic of Of stars – the line-flux is rather modest in comparison to what is found for the next two cases with values exceeding the critical value of 0.7. These objects show very strong and broad He ii 4686Å emission lines that are more characteristic of Of/WN or “slash” stars, progressing towards the Wolf-Rayet stars of the nitrogen sequence (WN).

These models thus indicate that the observed spectral transition from Of to WN corresponds to a transition from relatively low  to high  values (and larger ) for WN stars. This assertion is based not only on the higher predicted mass-loss rates themselves, but also on the finding that the mass-loss behaviour (as a function of  ) changes at 0.7. We note that the increasing Heii 4686Å equivalent width (EW) amounts to EW values of 2, 7, 20, respectively (for these unclumped models).

## 6 Discussion

### 6.1 Comparison with empirical mass-loss rates

Comparing our new mass-loss predictions against observed mass-loss rates is a non-trivial undertaking, as high objects are scarce. The largest sample of such potentially high objects that involves state-of-the-art modelling analysis is probably that of the Arches cluster by Martins et al. (2008). They provided stellar and wind properties (accounting for wind clumping) of 28 of its brightest members from -band spectroscopy. Roughly half of their sample comprises O4-O6 supergiants whilst the other half includes H-rich WN7-9 stars.

It is not possible to quote direct mass-loss predictions, because the Martins et al. analysis did not yield object masses. However, on the basis of the high stellar luminosities, with ) up to 6.3, these objects were suggested to be consistent with initial masses of up to 120. For the O4-6 supergiant population, luminosity values are in the range ) 5.75–6.05, consistent with initial masses . For this mass and luminosity range, 0.2 – comfortably within our low- regime. Assuming the current mass of these objects is about the same as their initial mass, our mass-loss formula yields values of 6.1, which is in reasonable agreement with the lower end of the Martins et al. mass-loss rates for their O4-O6 I objects.

The second group of Martins et al. objects comprise the WN7-9 objects. If we again assume their current masses can directly be inferred from the observed luminosities, we find 0.4 and mass-loss rates . Even if the helium abundances of these objects are increased, these properties would not result in a pronounced emission profile of the He ii line, i.e. a profile shape that is typical of late-WN stars. For this to happen the mass-loss rates need to be higher by at least a factor of a few. This seems to require high . In the framework of our models, this could be achieved by lowering the mass. However as the non-electron contribution in our models is not self-consistently treated in the high regime, we refrain from providing quantitative assessments.

As the O4-O6 supergiants from the Martins et al. (2008) analysis have effective temperatures in the range 32-40 kK, we expect to fall in the range 2500-3500 km s (see Fig. 8), which reasonably agrees with the upper end of the Arches O4-O6 supergiant stars. However, for the late WN stars, the terminal velocities presented by Martins et al. drop to values as low as 800-1600 km s, which is a factor 2-4 lower than we predict. We identify a number of possible reasons for this discrepancy. One option could be that the K-band spectral fits of Martins et al. (2008) yield terminal velocities that are too low (note that no ultra violet P Cygni blue edges are available for these obscured objects), but a more plausible reason is that as a result of our modelling assumptions (no rotation, smooth winds, etc.), we overpredict the wind terminal velocity. As the overprediction in terminal velocity is a factor 2-4, this would naturally translate in a mass-loss underestimate of a factor 3-5.

Although the effects of rotation on our mass-loss predictions are beyond the scope of this investigation, rotation may become relevant once one wishes to compare the non-rotating predictions to observations of objects that might rotate at a relevant rate. In particular, given that Gräfener et al. (2011) discuss the possibility that the WNh stars in the Arches cluster may evolve close to being chemically homogeneous. To first order, one might expect the effect of rotation to lower the effective gravity; i.e. one could replace the mass by the effective mass , in which case one would anticipate to scale as . In a similar vein, one might expect the terminal velocity to scale with the escape velocity, i.e. 3 (see e.g. Gayley 2000; Puls et al. 2008 and references therein). In Sect. 4.5, we computed some test models in which we lowered the stellar mass, hence the effective gravity, so as to keep  constant (in order to study a potential helium dependence on ). These models indeed showed higher mass-loss rates for lower effective gravity, but can we quantify the effective-gravity effect?

If we were to attribute the offset between the observed and predicted wind-terminal velocity of a factor 2 to stellar rotation, we should have an of 0.85 to “match” theory to observations. In that case, the mass-loss rate would be enhanced by a factor of for values in the range 0.5-0.6. Nevertheless, this could significantly change for lower temperatures. These high values of the order of 0.85 seem to be rather high when we consider recent stellar evolution models (Brott et al. 2011; Friedrich et al. in prep.), but they are not out of line with respect to spectral observations of LBVs (Groh et al. 2006), objects that are presumably in close proximity to the Eddington limit.

In all, we interpret our high wind terminal velocities as a sign that some physics is missing in our high models. Therefore, we refrain from using our dynamically determined mass-loss rates in a quantitative way at the optically-thick high end. For normal O-type stars, for low and moderate values, we achieve much better agreement between observed and predicted wind terminal velocities (Muijres et al. 2011b), and for this regime we have a much higher confidence in the accuracy of the absolute mass-loss rates, as long as O-type winds are not extremely porous, in which case mass-loss rates could drop significantly (Muijres et al. 2011a).

### 6.2 Comparison to other models

#### Comparison to CAK and other O-type star mass-loss models

We now wish to compare our results with previous model predictions. In this paper, we have investigated the mass-loss behaviour at high  for an extensive grid of models, and we revealed the existence of two mass-loss regimes. Moreover, we have found that mass-loss rates are dependent on both  and stellar mass (or stellar luminosity) and that the shape of these dependencies is well-described by a power law.

We remind the reader of classical CAK theory, where the mass-loss rate is found to be proportional to

 ˙M ∝ L (Γe1 − Γe)1 − αα (4)

where is a force multiplier parameter expressing the importance of optically thin lines to the total ensemble of lines. It is generally found that is 2/3 for galactic O-type stars (Puls et al. 2008) and assumed to be constant throughout the atmosphere. In reality, however, is depth-dependent (Vink 2000, Kudritzki 2002, Gräfener & Hamann 2005, Muijres et al. 2011b), which is captured better by an alternative representation of the line acceleration (Müller & Vink 2008). Nevertheless, the classical CAK formalism – as described by Eq. 4 – already shows a dependence on both and , and one could rewrite this mass-loss dependence as a function of and , using a mass-luminosity relation.

In the standard Vink et al. (2000) mass-loss parametrization , which can be reorganized to . This is the type of mass-loss parametrization that is currently employed in modern evolutionary computations (see e.g Meynet & Maeder 2003, Palacios et al 2005, Limongi & Chieffi 2006, Eldridge & Vink 2006, Vink et al. 2010, and Brott et al. 2011).

In this paper, we obtain much steeper  vs. dependencies, in agreement with our previous models for constant-luminosity LBVs (Vink & de Koter 2002, Smith et al. 2004). For the “low” range considered here we obtained “modest” dependencies of , but for the “high” regime, we found which involves a much steeper dependence on than any CAK-type mass-loss relationship provides. We note that such a steep dependence on the Eddington limit agrees with radiation-driven wind models of Vink (2006) and Gräfener & Hamann (2008), whilst Gräfener et al. (in prep.) also provide empirical evidence for such a strong mass-loss dependence on the Eddington parameter.

We emphasize that, through the use of the Vink et al. (2000) theoretical mass-loss recipe, most current stellar models already include the effect of positive mass-loss feedback (contrary to recent claims by Smith & Conti 2008). This effect describes how the mass-loss rate increases with the Eddington parameter. However, as we here obtain much steeper  vs. dependencies, it is likely that the mass-loss feedback effect that is currently employed in the stellar evolution models may not be sufficient for certain areas of the Hertzsprung-Russell diagram. We thus concur with the notion of Smith & Conti (2008) that new stellar evolution computations that take this effect into account properly are desirable

#### Comparison to alternative Wolf-Rayet mass-loss models

We also compare our models to the optically thick wind models for Wolf-Rayet stars, such as the critical-point analysis by Nugis & Lamers (2002) and the hydrodynamical model atmosphere analysis of Gräfener & Hamann (2008). As there is a significant qualitative difference between our Monte Carlo approach and the optically thick wind approaches, a meaningful quantitative comparison is a non-trivial undertaking (see section 3.1).

First, we quantitatively compare our versus dependence to the WNL star mass-loss dependence suggested by Gräfener & Hamann (2008). For the models in our grid at = 50 kK, we find very good agreement with the Gräfener & Hamann (2008) mass-loss rates and also find that the power-law slope of our dependence is very similar. However, the onset of WR-type behaviour occurs earlier, i.e. for lower , in the models by Gräfener & Hamann.

In Sect. 3, we discussed the possibility of such a shift in , because the actual Eddington parameter is expected to be affected by free-free and bound-free contributions and peaks in the iron opacity. By comparison with OPAL opacity tables (Iglesias & Rogers 1996), we estimate an increase of by in the region of the sonic point, assuming the location of the sonic point remains unaffected. This value corresponds roughly to a shift in between our relation and that by Gräfener & Hamann for typical parameters of Galactic WNL stars ( kK, ). This could be considered a maximum shift as we may slightly overestimate the line force near the sonic point by applying the Sobolev approximation (Pauldrach et al. 1986). However, a change in affects the atmospheric structure and therefore the location of the sonic point, consequently the effect on is hard to establish.

The Gräfener & Hamann mass-loss rates also display a strong temperature dependence, with . The actual size of the shift in is thus strongly dependent on the specific stellar parameters. Our Monte Carlo models suggest a much smoother dependence on (see Fig. 7) as long as we stay above the location of the predicted bi-stability jump, where the mass-loss properties jump drastically (Vink et al. 1999, Pauldrach & Puls 1990).

What we wish to emphasize is that both modelling approaches show an versus dependence that is much stronger than any additional mass or luminosity dependence. Where the two distinct mass-loss prescriptions differ is in the treatment of wind clumping, as well as the value of for the onset of WR-type mass loss behaviour. The exact location of this transition is of paramount importance for the evolution of the most massive stars. Ultimately, this should be testable with comparisons to observational data when sufficient objects are available in the appropriate range. This will be a crucial aim of the VLT Flames Tarantula Survey ( Evans et al. 2011, Bestenlehner et al. in prep.).

## 7 Summary

We presented mass-loss predictions from Monte Carlo radiative transfer models for very massive stars in the mass range 40-300 and with Eddington factors  in the range 0.4–1.0. An important outcome is that when winds become optically thick the spectral and mass-loss properties change. This transitional behaviour can be summarized as follows

• We find a transition from common O-type stars to more extreme Wolf-Rayet behaviour when  exceeds a critical value.

• The way in which the mass-loss rate depends on  in the range 0.4 0.7 is   , where rates are found to be consistent with the standard Vink et al. (2000) mass-loss rates.

• At the  dependence shows a “kink”; i.e. the slope is steeper for objects closer to the Eddington limit. Here the slope becomes   . This slope agrees with WNL models by Gräfener & Hamann (2008).

• When approaches unity, the wind efficiency number rises in a curved manner to values as high as 2.5. Such high values are more commensurate with Wolf-Rayet winds than with common O stars winds, and these results thus confirm a natural extension from common O-type mass loss to more extreme WR behaviour.

• This transitional behaviour is also found in terms of the wind acceleration parameter , which naturally reaches values as high as 1.5

• The spectral morphology of the He ii line at 4686Å changes gradually as a function of . This links the spectral sequence O-Of-Of/WN-WN to a transition of optically thin to optically thick winds.

• The mass-loss rate is found to be only modestly dependent on the effective temperature for the range of 30 to 50 kK.

• Last but not least, we highlight the fact that for fixed  the He abundance only has a minor effect on the predicted rate of mass loss (cf. Vink & de Koter 2002 for LBVs). This contradicts how O-type, to LBV-type, to WR-type mass-loss transitions are employed in massive star evolution models. We thus call for fundamental changes in the way mass loss is included in stellar evolution models for objects in close proximity to the Eddington limit.

###### Acknowledgements.
We would like to thank the anonymous referee for providing constructive comments.

### Footnotes

1. offprints: Jorick S. Vink, jsv@arm.ac.uk
2. In reality,  changes once hydrogen recombines (which starts below 30 000 K), or when the hydrogen-to-helium surface abundance changes, relevant for classical WR stars.
3. The issue of rotation in radiation-driven wind is highly complex. All existing studies (e.g. Friend & Abbott 1986, Bjorkman & Cassinelli 1993, Petrenz & Puls 2000, Pelupessy et al. 2000, Cure & Rial 2004, Madura et al. 2007) have only treated part of the problem, but not tackled the combined multi-dimensional, high (), and high aspects.

### References

1. Abbott D.C., Lucy L.B. 1985, ApJ 288, 679
2. Abbott D.C., Bieging, J.H., Churchwell, E., Torres, A.V., 1986, ApJ 303, 239
3. Allen C.W., 1973, Astrophysical quantities, University of London, Athlone Press
4. Anders, E. & Grevesse N., 1989, GeCoA 53, 197
5. Baraffe I., Heger A., Woosley S.E. 2001, ApJ 550, 890
6. Barlow M.J., Smith L.J., Willis A.J., 1981, MNRAS 196, 101
7. Belkus H., Van Bever J., Vanbeveren D., 2007, ApJ 659, 1576
8. Bestenlehner J.M., Vink J.S., Gräfener G., et al., A&AL, in press
9. Bjorkman J.E., Cassinelli J.P., 1993, ApJ 409, 429
10. Bouret S.-C., Lanz T., Hillier D.J., 2003, ApJ 595, 1182
11. Bromm V., Kudritzki R.P., Loeb A., 2001, ApJ 552, 464
12. Brott I. et al., 2011, A&A, arXiv1102.0530
13. Cantiello M., Langer N., Brott I., de Koter A., Shore S. N., Vink J. S., Voegler A., Lennon D. J., Yoon S.-C., 2009, A&A 499, 279
14. Castor J., Abbott D.C., Klein R.I., 1975, ApJ 195, 157
15. Conti P.S., 1976, MSRSL 9, 193
16. Crowther P.A., Dessart L., 1998, MNRAS 296, 622
17. Crowther P.A., Schnurr O., Hirschi R., Yusof N., Parker R.J., Goodwin S.P., Kassim H.A., 2010, MNRAS in press
18. Cure M., Rial D.F., 2004, A&A 428, 545
19. Davies B., Oudmaijer R.D., Vink J.S., 2005, A&A 439, 1107
20. Davies B., Vink J.S., Oudmaijer R.D., 2007, A&A 469, 1045
21. de Koter, A., Schmutz, W., Lamers, H.J.G.L.M., 1993, A&A 277, 561
22. de Koter A., Heap S.R., Hubeny I., 1997, ApJ 477, 792
23. Eldridge, J.J., & Vink, J.S. 2006, A&A, 452, 295
24. Evans C.J., Taylor W.D., Henault-Brunet V., 2011, A&A, arXiv1103.5386
25. Figer D.F., Najarro F., Morris M., et al., 1998, ApJ 506, 384
26. Figer D.F., 2005, Nature 434, 192
27. Friend D.B., Abbott D.C., 1986, ApJ 311, 701
28. Fullerton, A. W., Massa, D. L., Prinja, R. K. 2006, ApJ 637, 1025
29. Gayley K.G., 2000, ApJ 529, 1019
30. Gayley K.G., Owocki S.P., Cranmer S.R., 1995, ApJ 442, 296
31. Glatzel, W., Kiriakidis, M., 1993, MNRAS 263, 375
32. Gräfener G.,& Hamann W.-R., 2005, A&A 432 633
33. Gräfener G., Hamann W.-R. 2008, A&A 482, 945
34. Gräfener G., Vink J.S., de Koter A., Langer N., 2011, A&A submitted
35. Groh J.H., Hillier D.J., Damineli A., 2006, ApJ 638, 33
36. Hamann W.-R., & Koesterke L., 1998 A&A 335, 1003
37. Harries T.J., Howarth I.D., Evans C.J., 2002, MNRAS 337, 341
38. Heger A., & Langer N., 1996, A&A 315, 421
39. Hillier D.J., 1991, A&A 247, 455
40. Humphreys R.M., Davidson K. 1994, PASP 106, 1025
41. Iglesias, C.A., Rogers, F.J., 1996, ApJ 464, 943
42. Ignace R., Quigley M.F., Cassinelli J.P., 2003, ApJ 596, 538
43. Ishii, M., Ueno, M., Kato, M., 1999, PASJ 51, 417
44. Kudritzki R.-P., 2002, ApJ 577, 389
45. Kurucz R.L. & Bell B. Atomic Line Data Kurucz CD-ROM No. 23. Cambridge, Mass.: Smithsonian Astrophysical Observatory, 1995
46. Lamers H.J.G.L.M., Leitherer C., 1993, ApJ 412, 771
47. Langer, N., 1997, ASPC 120, 83
48. Langer, N., 1989, A&A 220, 135
49. Leitherer C., Robert C., Drissen L., 1992, ApJ 401, 596
50. Limongi, M., Chieffi, A. 2006, ApJ 647, 483
51. Lucy L.B., & Solomon P.M., 1970, ApJ 159, 879
52. Lucy L.B., Abbott D.C., 1993, ApJ 405, 738
53. Madura T.I., Owocki S.P., Feldmeier A., 2007, ApJ 660, 687
54. Maeder A., & Meynet G., 1987, A&A 182, 243
55. Maeder A., & Meynet G., 2000, A&A 361, 159
56. Martins F., Hillier D.J., Paumard T., Eisenhauer F., Ott T., Genzel R., 2008, A&A 478, 219
57. Meynet G., & Maeder A., 2003, A&A 404, 975
58. Moffat A.F.J., & Robert C., 1994, ApJ 421, 310
59. Mokiem M.R., de Koter A., Vink J.S. 2007, et al., A&A 473, 603
60. Muijres L., de Koter A., Vink J.S., et al., 2011a, A&A 526, 32
61. Muijres L., Vink J.S., de Koter A., et al., 2011b, A&A submitted
62. Müller P.E., Vink J.S., 2008, A&A 492, 493
63. Nugis T., & Lamers H.J.G.L.M., 2002, A&A 389, 162
64. Oey M.S., Clarke C.J., 2005, ApJ 620, 43
65. Owocki S.P., Gayley K.G., Shaviv N.J., 2004, ApJ 616, 525
66. Palacios A., Arnould M., Meynet G., 2005, A&A 443, 243
67. Pauldrach A.W.A., Puls J., Kudritzki R.P., 1986, A&A 164, 86
68. Pauldrach, A.W.A. & Puls, J. 1990, A&A 237, 409
69. Pelupessy I., Lamers H.J.G.L.M., Vink J.S., 2000, A&A 359, 695
70. Petrenz P., Puls J., 2000, A&A 358, 956
71. Petrovic J., Pols O., Langer N., 2006, A&A 450, 219
72. Puls J., Kudritzki R.P., Herrero A., et al., 1996, A&A 305, 171
73. Puls J., Vink J.S., Najarro F., 2008, A&ARv 16, 209
74. Schmutz W., 1991, in: “Stellar Atmospheres: Beyond Classical Models”, eds. Crivellari L., Hubeny I., Hummer D.G., NATO ASI Series C, Vol. 341, 191
75. Shaviv N. J.,1998, ApJ 494, L193
76. Shaviv N. J., 2000, ApJ 532, L137
77. Smith N., Conti P.S., 2008, ApJ 679, 1467
78. Smith N., Vink J.S., de Koter A., 2004, ApJ 615, 475
79. Springmann, U., 1994, A&A 289, 505
80. van Marle A. J., Owocki S. P., Shaviv N. J., 2008, MNRAS 389, 1353
81. Vink J.S., de Koter A., Lamers H.J.G.L.M. 1999, A&A 345, 109
82. Vink J.S., de Koter A., Lamers H.J.G.L.M. 2000, A&A 362, 295
83. Vink J.S., 2000, PhD thesis at Utrecht University
84. Vink J.S., de Koter A., Lamers H.J.G.L.M. 2001, A&A 369, 574
85. Vink J.S., de Koter A. 2002, A&A 393, 543
86. Vink J.S., de Koter A. 2005, A&A 442 587
87. Vink, J.S. 2006, ASPC 353, 113 (astro-ph/0511048)
88. Vink J.S., 2009, in: “Eta Carinae”, eds R.M. Humphreys & K. Davidson, Springer, astro-ph/0905.3338
89. Vink J.S., Davies B., Harries T.J., Oudmaijer R.D., Walborn N.R., 2009, A&A 505, 743
90. Vink J. S., Brott I., Gräfener G., Langer N., de Koter A., Lennon D. J., 2010, A&A 512, L7
91. Walborn N.R., Ebbets D.C., Parker J.W., Nichols-Bohlin J., White R.L., 1992, ApJ 393, 13
92. Weidner C., Kroupa P., 2004, MNRAS 348, 187
93. Yungelson L.R., van den Heuvel E.P.J., Vink J.S., Portegies Zwart S.F., de Koter A., 2008, A&A 477, 223
248364