Massive star evolution : rotation, winds, and
overshooting vectors in the Mass-Luminosity plane
Key Words.:stars: massive – evolution – mass loss
email: email@example.com; firstname.lastname@example.org
Context:Massive star evolution is dominated by various physical effects, including mass loss, overshooting, and rotation, but the prescriptions of their effects are poorly constrained, even affecting our understanding of the main sequence.
Aims:We aim to constrain massive star evolution models using the unique testbed eclipsing binary HD 166734 with new grids of MESA stellar evolution models, adopting calibrated prescriptions of overshooting, mass loss, and rotation.
Methods:We introduce a novel tool: the ”mass-luminosity plane” or ” plane”, as an equivalent to the traditional HR diagram, utilising it to reproduce the testbed binary HD 166734 with newly calibrated MESA stellar evolution models for single stars.
Results:We can only reproduce the Galactic binary system with an enhanced amount of core overshooting ( 0.5), mass loss, and rotational mixing. We can utilise the gradient in the plane to constrain the amount of mass loss to 0.5 - 1.5 times the standard Vink et al. (2001) prescriptions, and we can exclude extreme reduction or multiplication factors. The extent of the vectors in the plane leads us to conclude that the amount of core overshooting is larger than is normally adopted in contemporary massive star evolution models. We furthermore conclude that rotational mixing is mandatory to get the nitrogen abundance ratios between the primary and secondary components to be correct (3:1) in our testbed binary system.
Conclusions:Our calibrated grid of models, alongside our new plane approach, present the possibility of a widened main sequence due to an increased demand for core overshooting. The increased amount of core overshooting is not only needed to explain the extended main sequence, but the enhanced overshooting is also needed to explain the location of the upper-luminosity limit of the red supergiants. Finally, the increased amount of core overshooting has – via the compactness parameter – implications for supernova explodibility.
Massive stars with an initial mass above 8 have a diversity of possible evolutionary channels, dictated by the dominant processes acting on their structure. The extent of these dependancies are variant with mass, metallicity, and multiplicity. Stellar winds have a significant impact on the evolution of O-type stars throughout their lives, leading to evolutionary phases involving Luminous Blue Variables (LBV) and Wolf-Rayet (WR) stars. It is also an important factor for dictating their final masses, determining whether a neutron star or black hole is formed in the final stage of evolution, as extensively reviewed by Chiosi & Maeder (1986).
On the main sequence (MS), mass loss via stellar winds has the greatest impact at the highest mass ranges: above 60 mass loss completely dominates the evolution (e.g. Vink & Gräfener 2012; Vink 2015), whilst in the range 30 60 it is one of the important ingredients (e.g. Langer 2012; Groh et al. 2014). At lower masses, (i.e. below 30 ) the evolution is thought to be heavily influenced by rotation (e.g. Maeder & Meynet 2000). Over the explored mass range within this paper, 8 - 60 , we will consider both effects, as well as convective overshooting, which may all play a role.
The extension of the convective core by overshooting is a key structural feature which increases the amount of hydrogen (H) dredged into the core, replenishing its supply, thereby extending the MS lifetime. The parameter which we explore in this study corresponds to the fraction of the pressure scale height H by which particles will continue to travel a distance beyond the convective core boundary. This form of mixing has been explored for decades, with few constraints on its size () in the high-mass range. It has been argued essential for reproducing observations, though lacks evidence for dependancies such as mass (Claret & Torres 2017). Another process which may potentially affect stellar evolution is the presence of a magnetic field, however, Grunhut et al. (2012) show the fraction of magnetic O stars to be just on the order of 7%.
Massive star evolution models are currently not able to fully reproduce observations, even the MS (e.g. Vink et al. 2010; Markova et al. 2018), with many physical processes such as rotation and overshooting yet to be fully understood. The MS width and dependancies remain unresolved, while this stage represents 90% of the overall lifetime – and sets the stage for later evolutionary phases.
Efforts have been made to map out the evolution of massive stars with systematic grids (e.g. Brott et al. 2011a; Ekström et al. 2012) using detailed predictions of chemical abundances, rotation rates, and fundamental parameters such as mass, luminosity, and effective temperature. These models have subsequently been compared to observations for predicting evolutionary stages and characteristics, though due to limitations in both key observations and accurately modelling key physical processes, many assumptions remain, including the amount of core overshooting () that is thought appropriate.
Martins & Palacios (2013) explored a diversity of evolutionary codes (e.g. Ekström et al. 2012; Chieffi & Limongi 2013; Bertelli et al. 2009) in which the implementation of input physics was surveyed, allowing code applicability to be tested, however linear comparisons of physical treatments cannot be drawn due to the variety of prescriptions in different codes. It is clear that all stellar evolution models will have a degree of uncertainty, yet to establish a clear comparison between codes, it would be beneficial if physical implementations are examined with one and the same code. Therefore, we here aim to compute massive star models with both new and existing prescriptions using the same evolutionary code, MESA, ”Modules for Experiments in Stellar Astrophysics” (e.g. Paxton et al. 2011) given its high flexibility and code capabilities, enabling ample comparisons of several key physical processes. Such exploration offers the opportunity for calibrating models with respect to observations.
Due to the variety of possible prescriptions in each code, the evolution of massive O-type stars so far remains model dependent, leaving the MS lifetime ambiguous particularly due to the absence of evidence for objects after the terminal age main sequence (TAMS). Cool B supergiants, the descendants of O-type stars, are less understood, and have yet to be confirmed as core hydrogen or helium burning objects, (Vink et al. 2010). As O-type stars spend the majority of their lifetime on the MS, we would expect a scarcity of B supergiants if they are indeed post-MS objects. However, we observe a (too) large number of these stars (e.g. Garmany & Fitzpatrick 1988), raising the possibility of these objects being MS, core H-burning stars. The existence of a large number of slow rotating B supergiants however, (with 50 km s) is suggestive of an evolved star which has completed the MS phase and been spun down.
Vink et al. (2010) also consider the possibility for bi-stability braking (BSB) as the mechanism by which B supergiants lose their angular momentum, (see also Keszthelyi et al. 2017). If we consider that B supergiants may not represent the end of the H burning phase, this could allow for a wider MS, hypothesised by Vink et al. (2010). This would result in a demand for additional mixing of H in the core, which may be fulfilled by increased overshooting. Vink et al. (2010) addresses that a higher value of would result in a lower critical mass at which BSB would be efficient. Test models show that BSB occurs in present models with 0.335 above a critical mass of 35 in the Large Magellanic Cloud, yet with an increase of overshooting to 0.5 the critical mass drops to 20 for the same metallicity, (Vink et al. 2010).
The determination of for massive stars has been challenging without the aid of astroseismological data for the most massive stars, leading to an array of prescriptions such as the correlation between and log g , Brott et al. (2011a). Many other estimations of have been adapted in stellar evolutionary models leading to a wide variety of potential stellar ages, MS lifetimes and final products, Martins & Palacios (2013). One of the most straightforward ways would be to derive it simply from the MS width, which might potentially be possible from the Galactic Hertzsprung-Russell Diagram (HRD) in Castro et al. (2014) although such Galactic data might be biased, as the unbiased LMC data from the VLT-Flames Survey of massive stars Evans et al. (2005) do not show a gap between O and B supergiants Vink et al. (2010) suggesting a more extended MS.
In this paper, we attempt to constrain the dominant parameters effecting massive star evolution. For this purpose, we present a grid of evolutionary models for two extreme values of 0.1 and 0.5 to illustrate both lower adopted values as well as enhanced overshooting, with varying initial masses, rotation rates, and mass-loss rates, thereby highlighting the sensitivity of stellar models in terms of mixing and mass loss. We introduce the ’Mass-Luminosity Plane’ as an alternative to the HRD to study the key ingredients in massive star evolution on the main sequence (see Fig. 4). Whilst the fundamental stellar parameters of Mass and Luminosity have been plotted logarithmically by Maeder (1983), for example, our version of the plot highlights the independent effects of rotation, overshooting and mass loss on stellar evolution through vectors, with inverted mass on the x-axis providing a useful comparison to the tracks in the HRD.
Weidner & Vink (2010) present an overview of the methods of mass determination for O stars, including the ’mass discrepancy’ often seen between the evolutionary masses and spectroscopic masses. The method of comparing the positions of stars in the HRD with theoretical evolution models (evolutionary masses) has often led to predictions which are systematically higher than the masses derived through stellar spectroscopy (spectroscopic masses), (e.g. Herrero et al. 1992). However, when O stars are found in binary systems, their dynamics can present a model independent mass determination (dynamical masses). Evolutionary masses can present discrepancies amongst themselves when using various theoretical models (e.g. Ekström et al. 2012; Brott et al. 2011a) with differing implementations of rotation, convection and mass loss. Though this is not the widely discussed âmass discrepancyâ problem, it does highlight the necessity of calibrating stellar evolution models to minimise further discrepancies with spectroscopic and dynamical masses (see e.g. Markova et al. 2018). In the case where dynamical masses agree with the spectroscopic masses, we can have faith in the spectroscopic result, thus allowing for calibration of theoretical evolution models. Similar work has been completed by Southworth et al. (2004); Pavlovski et al. (2018); Tkachenko et al. (2014) for detached eclipsing binaries, however these works utilised lower mass stars (up to 15 ) which did not incorporate the interacting effects of mass loss, overshooting and rotation as we do in this study.
We use constraints relative to and to investigate the possible evolutionary paths of a high mass, detached binary, HD 166734, modelled here as a testbed for single star evolution. As previously mentioned, dominant processes take effect at varying mass ranges, yet with dynamical masses of 39.5 and 33.5 , for the primary and secondary respectively, for HD 166734 we may probe the effects of these processes as they interact and overlap. As the spectroscopic masses adeptly agree with the dynamical masses for HD 166734 (Mahy et al. 2017), this system provides a unique opportunity to constrain – and effectively correct – stellar evolution models, whilst for the general case of single massive stars we cannot currently tell if there are issues with spectroscopic masses, leading to mass discrepancies (Markova et al. 2018).
We present a method of producing a calibrated grid of models with an analysis of HD 166734 in section 2.1, with calibration of mixing processes in section 3. We explore a new tool for comparing observations with models in the Mass-Luminosity Plane in section 4 and we provide our final results for HD 166734 in section 5. We present our grid of models alongside a sample of Galactic O-stars in section 6, with further results in Appendix A. Finally, we highlight our conclusions in section 7, with remaining full grid tables in Appendix A.
2.1 MESA : Treatment of convection, mass loss, and rotation
A set of evolutionary models was calculated for massive main-sequence stars with the one-dimensional, stellar evolution code MESA, e.g. Paxton et al. (2011), as a comparison for both the primary and secondary of HD 166734 (see Sect. 8). The extensive capabilities of this code provide a diverse range of available alterations, enabling the user to compare implementations of physical processes with other code treatments. In this paper, we examine the effects of mass loss, convective overshooting, and rotational mixing in terms of fundamental observables such as luminosity, mass, and surface abundances. These models were completed from zero-age main-sequence (ZAMS) to core collapse, unless convergence problems arise in which computations are concluded earlier. We adopt the default metallicity in MESA of Z = 0.02 with the chemical mixture from Grevesse & Sauval (1998) in order to provide direct comparisons with chemical abundances in Galactic observations.
Convection is treated by the mixing length theory where 1.5, with a semi-convection efficiency parameter of 1. The convective core boundary is defined by the Ledoux criterion111 The Ledoux criterion is denoted by , but in chemically homogeneous layers where = the Schwarzschild criterion is effective. , with overshooting succeeding convective mixing at the core boundary, increasing the temperature gradient by implementing a thermal gradient , (e.g. Choi et al. 2016). This method of extending the core is denoted as step-overshooting, which enhances the core by a factor of the pressure scale height . Experiments in the dependancies of this parameter are completed in the following sections.
We later compare with treatments of and rotational mixing from Brott et al. (2011a) and Ekström et al. (2012) grids since these are used extensively in the community. Brott et al. (2011a) present a calibration of the overshooting parameter by comparing the TAMS of 16 models with observations from the FLAMES survey (Evans et al. 2008), suggesting a TAMS at log g 3.2, since this value represents a drop in beyond which a large number of slow rotating B supergiants are located, assumed to be post-MS objects. The model with 0.335 corresponded to the log g 3.2 and has since been used as a static parameter in models to compare against observations of a wide mass range. A lower value of 0.1 is applied for models presented by Ekström et al. (2012), since calibration was completed with lower mass stars of 1.7-2 where convective mixing plays a dominant role compared to that of rotational mixing, hence allowing a linear calibration of convective overshooting without accounting for the more sophisticated treatment of rotational mixing as prescribed in the GENEC code. We adopt step-overshooting for H-burning phases only, as we aim to better understand the MS-width.
Mass-loss rates are adopted from Vink et al. (2001) accounting for metallicity dependencies and the occurrence of the bi-stability jump, an increase of mass loss at 21kK causing effects in the evolution, seen in the HRD. Various factors of this mass loss regime will be tested to determine the possibility of extreme rates. We hence explore a range of multiplication factors of Vink et al. (2001) mass-loss rates from 0.1 to 10 times the standard prescription. Rotation is later applied in our models through a fully diffusive approach with appropriate instabilities such as the Eddington-Sweet circulation, dynamical and secular shear instabilities. We also consider the effects of an internal magnetic field by a Spruit-Taylor dynamo, although we found that this had inconsequential effects on our results. The calibration of our single star models are relevant for evolutionary codes which implement rotational mixing in a similar way, if this process is treated physically different in another code, then results would differ quantitatively, but qualitatively have the same behaviour.
|[ ]||8, 12, 16, 20, 25, 30, 35, 40, 45, 50, 55, 60|
|[ km s]||0, 100, 200, 300, 400, 500|
A systematic grid of models was calculated for comparison with a larger sample, including new prescriptions discussed in Sections 3 and 4. Table 1 shows the range of masses, rotation rates, and overshooting values for which we compose our grid. We choose masses representative for the O-star and early B-star range, with a variety of rotation rates up to break-up speed, and extreme values for to explore the extent of extra mixing. We evolve each model to core collapse, unless convergence problems highlight unlikely solutions. For this purpose, Vink et al. (2001) provides the relevant mass loss prescription, with a factor of unity for all models in the first instance.
2.2 The detached, eclipsing binary HD 166734 : a testbed for massive star evolution
The eclipsing massive binary HD 166734 (see Table 2) provides a unique opportunity to improve physics in stellar evolution models, as Mahy et al. (2017) were able to determine the individual stellar parameters, including their exact positions in the HRD and their dynamical real masses directly. As these dynamical masses were found to be in excellent agreement with their spectroscopic masses, these 2 stars of this massive binary system, enable us to calibrate and correct the evolutionary masses, thereby constraining the relevant physics in the upper HRD for stars above 30-40 . Observations of high-mass eclipsing binaries are sparse, and even more extreme for detached, non-interacting stars which may be treated as evolved single stars. As observations of massive single stars may sometimes highlight discrepancies between spectroscopic and evolutionary masses, we have an ideal opportunity here where the dynamical masses are in agreement with spectroscopic masses, providing a tool for calibrating evolutionary masses and thus evolutionary paths of stars that are massive enough for the physics to be heavily influenced, if not dominated by mass loss via stellar winds.
Though a large fraction of O stars may be present in a binary or multiple system, observations of eclipsing binaries above 30 are extremely rare (see e.g. Bonanos et al. 2004; de Mink et al. 2009; Pfuhl et al. 2014; Gies 2003). Hence the stellar parameters derived by Mahy et al. (2017) have provided a unique opportunity to analyse a non-interacting system which can be treated as single stars. The similar values of for both components may at face value be considered of interest in terms of synchronisation, but as the rotation speeds are lower than the orbital period, Mahy et al. (2017) gave arguments arguing against synchronisation. In addition, we note that the values are close to the inferred macro-turbulent values of 65 10 km s(Mahy et al. 2017) and we therefore urge for caution that the quoted values of are truly the result of rotation (see Simón-Díaz & Herrero 2014). We thus treat the values as upper limits, and we consider the similar values of the 2 components as merely a coincidence.
|[K]||32000 1000||30500 1000|
|log(L/ )||5.840 0.092||5.732 0.104|
|[ ]||39.5 5.4||33.5 4.6|
|[ ]||37.7 29.2||31.8 26.6|
|[ km s]||95 10||98 10|
We utilise this agreement between dynamical and spectroscopic masses, allowing HD 166734 to be treated as an excellent testbed for massive star evolution of the most massive O-type stars. Mahy et al. (2017) analysed the system finding a composition of two supergiant O-type stars in an eccentric 34.5-day orbital period. We recognise that the estimated quantities may be upper limits due to the possibility of macro-turbulence. We can also compare with observed surface N abundances as a secondary assessment of potential rotation rates.
Comparisons to fixed current-day evolutionary sets of models by Brott et al. (2011a) and Ekström et al. (2012) by Mahy et al. (2017) revealed that both sets of models over-predict the evolutionary masses, whilst the secondary star appeared to be more evolved than the primary. We consider this latter finding an artefact of the Mahy et al. (2017) approach, and that in reality it is far more likely that both components formed simultaneously. We can therefore use an ”equal-age assumption” in addition to the exact HRD positions and true current day masses to solve the evolutionary mass discrepancy for both components, and at the same time constrain the relevant physics in this mass range.
Our assumption that this binary has evolved from the same initial stage is important for constraints of the MS width and thus for constraining the overshooting parameter, as well as the determination of rotation rates and possible evolutionary scenarios. As both stars show limited evidence of an evolved nature, we can exclude extreme events in the past such as eruptive mass-loss or binary interactions. Mahy et al. (2017) show surface nitrogen enrichments with a particle fraction [N/H] ratio of 3:1 between the primary and secondary components respectively. We utilise these abundances as evidence for mixing, as well as constraints for the determination of age.
3 Mixing and Mass Loss
3.1 Envelope stripping and nitrogen enrichment.
In developing our initial set of models we aim to minimise interacting physical processes. We start with a set of non-rotating stellar evolution models which exclusively employ mass loss and convective overshooting as mixing processes. In the first instance, initial masses were adopted from Mahy et al. (2017) with 56.1 and 47.4 for the primary and secondary respectively, with varying factors of the mass loss recipe, adopted from Vink et al. (2001), for a range of convective overshooting parameters . We initially attempt to reproduce characteristics of HD 166734 by following analysis from Mahy et al. (2017), with parameters taken from Brott et al. (2011a) and Ekström et al. (2012) grids, though find that these models do not offer solutions where sufficient N enrichment is reached. We hence employ greater mixing through increased factors of mass loss and overshooting.
In reproducing the properties of HD 166734, we can constrain the scenarios which display the 3:1 ratio of [N/H] for the primary to secondary by applying a restriction to the model time. As both stars are assumed to be approximately the same age with this ratio of enrichment, we can exclude the vast majority of possible evolutionary scenarios, i.e. those that do not represent these surface chemical enrichments simultaneously. Accordingly, we do not predict the ages of these stars, but we rather allow for constraints such as surface enrichments, rotation rates, and dynamical masses to provide a solution whereby both stars can reproduce the observables concurrently. Analogous to this, isochrones have not been used here as a method of stellar age determination as we have previously highlighted the sensitivity of model dependancy on these features, thus leading to a wide range of possible ages.
Massive stars produce surface He on the MS by the CNO-cycle, with a rapid increase in 14N by a factor of 10 at the surface when CN-equilibrium is reached. The occurrence of this observational feature has been reviewed widely by Maeder & Meynet (1987), finding that increased convective mixing by overshooting has shown to lower the limit for CN-equilibrium during the MS.
Maeder & Meynet (1987, 1988, 1991) composed grids of evolutionary models based on inputs of mass-loss rates and convective overshooting as the sole mechanisms for chemical mixing. The importance of convective overshooting has been stressed in these early publications as leads to a range of stellar ages, due to the dependence of T at TAMS on , (Maeder & Meynet 1991). Moreover, the MS luminosity increases by 0.9 dex at the reddest point of the MS when overshooting is accounted for leading to increases in age by factors of 1.5 - 2.7.
Meynet et al. (1994) present grids of massive stars with high mass-loss rates since the evolution of the most massive stars is so heavily reliant on the effect of stellar winds. A factor of two enhancement was applied to their mass-loss prescription from Schaller et al. (1992) demonstrating the effects on the evolutionary track presented in a HRD. These results hinted at a metallicity dependancy on mass-loss rates, though also show envelope stripping with increased mass loss leading to evolutionary phases such as WR types and quasi-chemical homogeneous evolution, Meynet et al. (1994). When analysing nitrogen enrichments for these models we find that if surface abundances do increase, it is by a sudden step of a factor of ten, representative of CN-equilibrium. This behaviour applies to factors of 1 - 3 of Vink et al. (2001) mass-loss rates, and overshooting of 0.1 - 0.8. We also note that models with increased overshooting result in earlier enrichment by up to 1Myr, regardless of mass-loss rates. In figure 1 we present the nitrogen enrichments for a sample of models of primary and secondary masses.
We find that chemical mixing of CNO elements by mass loss and overshooting attains CN-equilibrium before any intermediate enrichment occurs. This demonstrates that a combination of increased mass loss and overshooting will result in envelope stripping, whereby fusion products are extensively exposed at the stellar surface. Since this does not provide a solution for reproducing the observed surface enrichments of HD 166734, and moreover any observation with intermediate enrichment, we must explore additional, viable mixing processes, such as rotational mixing.
As we have adopted the method of step-overshooting, physical implications of this may hinder intermediate enrichment since step-overshooting invokes instantaneous homogeneous mixing within the overshooting region, leading to immediate enrichment by a factor of 10 when the envelope is stripped via stellar winds. Therefore, we compared our results with the prescription of exponential overshooting, whereby the length of the scale height is set by a comparable parameter f, but the overshooting region is mixed by a diffusion gradient.
Nevertheless, these results show similar enrichments, as even though intermediate enrichments may be reached through the overshooting region by altering the diffusion coefficient, elements are not mixed intermittently through the envelope from the overshooting layer. Thus another mixing process capable of mixing the chemical elements from the convective layers through the envelope must be implemented in order to match observed enrichments.
Recent studies of massive star observations (e.g. Brott et al. 2011b; Hunter et al. 2008; Maeder 2000) suggest that surface enrichments of CNO products may or may not be a result of rotational mixing. Yet, the necessity of rotational mixing has not been stressed with respect to CN-equilibrium or observed intermediate enrichments. We therefore tested the effects of rotational mixing as a function of surface enrichment, with a set of rotating models of varied initial rotation rates from 100-500 km s. In this set of models we find that a range of intermediate enrichments occurs, also providing solutions for reproducing the 3:1 nitrogen ratios as seen in HD 166734, (Fig.1). The comparison in Fig. 1 illustrates that rotational mixing is essential in reproducing observational surface enrichments, unless another not yet considered mechanism is identified, since previous mixing processes provide either too little or too much mixing leading to insignificant enrichment or CN-equilibrium.
Fig. 1. not only demonstrates the necessity of rotational mixing, but also stresses the importance of enhanced overshooting. In the rotating models of Fig. 1 we see that with an increase in from 0.1 to 0.5, we get much larger surface enrichments which may aid our understanding of the unexplained nitrogen enrichments discussed by Grin et al. (2017). As a significant fraction of the sample cannot be explained by rotational mixing alone, extended overshooting may help towards resolving this problem.
3.2 Rotationally-induced Mass Loss
While analysing a set of rotating models for HD 166734 we approach a problem with respect to interacting processes such as rotation and mass loss, consequently having a non-linear affect on the mass and luminosity. We find that the initial masses sufficient for reproducing the observed luminosities, are excessive when aiming to reach the dynamical masses by the time of observed temperatures or evolutionary phases. We therefore calculated a set of lower initial mass models, yet these diminish the luminosity gradient over time so that current data-points of HD 166734 remain out of reach. Interpreting an initial mass from the observed luminosity allowed for calculation of a possible mass-loss rate that would enable the current dynamical masses to be reached.
Following this method, we find a mass-loss rate of log 5.17, translating to an increase in mass-loss rate by approximately a factor of 3. We therefore completed further models with increased mass-loss rates of a factor of two and three. Though the dynamical masses were now reached, this also lead to a significant drop in luminosity, which correlates to a shallow gradient in the plane (see Fig. 5.), suggesting the observed masses and luminosities could not be reproduced simultaneously, (see Fig. 3). The possibility of rotationally-enhanced mass loss started with one dimensional radiation-driven wind models of Friend & Abbott (1986), who proposed an equatorially enhanced stellar wind as well as an increased mass-loss rate due to a lower effective gravity at the equator. It is this result that is also included in many massive star evolution models (see Heger et al. 2000; Brott et al. 2011a). This same implementation is included in the default MESA settings. The mathematical approach in shown in Eq. (1).
Note that the Geneva group (e.g. Maeder & Meynet 2015; Ekström et al. 2012) employ a slightly different implementation, yet it is based on similar physical principles. Since 1986 there have been many studies of the effects of rotation on radiation-driven wind predictions, with several different levels of sophistication, and different results. Recent 2D modelling by Müller & Vink (2014) encountered cases of equatorial decreases of the mass-loss rate, as well as surface-averaged total mass-loss rates that are lower than for the 1D case. They therefore challenged the implementation of rotationally-enhanced mass-loss in stellar evolution modelling, which is still mostly applied, e.g. it is the default setting in MESA.
Figure 2. highlights the change in initial mass-loss rate due to a change in initial rotation rate from 100 - 300 km s for both a 40 model. As we consider the current enhancement largely as artificial, we explored the difference between disabling rotationally enhanced mass loss (effectively setting 0) and enabling it using the default setting ( 0.43). We thus calculated a series of models with various initial masses, rotation rates, mass-loss rates and overshooting parameters.
4 Mass - Luminosity Plane
When comparing models in Section 3 we find that enhanced mass loss regimes lead to unrealistic luminosities which are too low to reproduce the observed HD 166734 luminosities. We also find that an initial mass representative of the observed luminosities is too high to reproduce the much lower dynamical masses with factor unity of Vink et al. (2001) mass-loss rates. If we aim to simultaneously reproduce the mass and luminosity, we must explore all possible dependancies of these properties:
where varies as a function of mass, and is the mean molecular weight.
The most fundamental characteristics of a star’s evolution are its mass and luminosity. As such, when trying to correlate the theoretical evolution of a star with its observables, these properties are essential. Thanks to analysis of HD 166734 by Mahy et al. (2017), we can reliably utilise the luminosities of both stars determined from bolometric magnitudes to calibrate their evolutionary status. This reasoning is also applicable to the masses of HD 166734, as in this circumstance the dynamical masses agree very well with the derived spectroscopic masses, providing a unique opportunity to constrain the mass-loss rates and physical processes during evolution.
The mass and luminosity of a star are reliant on age and mass-loss rate, so we reach a diversity of possible evolutionary scenarios with respect to mass-loss rates and . Yet we may constrain these solutions through assuming both objects evolved from the same initial starting point, so we can account for primary and secondary masses to be reached at the same time.
Eq. (2) shows that we can increase the luminosity by increased helium abundance. A minor helium enrichment in both the primary and secondary presents the possibility that the initial mass is not required to be insufficiently high to reach the dynamical mass. The observed helium enrichment corresponds to an increase in by 30 or 0.11 dex. Though this offers a potential scenario which would allow for a higher luminosity and lower , it is unlikely that the initial He abundance of HD 166734 is enriched rather than having been exposed as fusion products at the surface during hydrogen burning. We consider this solution unlikely.
Alternatively, the observed luminosity could be higher than would be required for the relevant initial mass due to the evolutionary phase at which these stars are currently undergoing. Beyond the TAMS, we observe an increased luminosity as models evolve to cooler temperatures. If HD 166734 was in fact composed of helium burning objects, the observed luminosity could be explained by this increased post-TAMS. Yet when comparing our models with the observed T’s, we note that both objects remain too hot to be post-MS objects, regardless of , thus excluding later evolutionary phases as a viable solution.
We observe some models which reach the dynamical mass of the primary due to higher mass-loss rates relative to the higher initial mass, though even these models must be excluded due to the observed T since the dynamical mass is only reached during the bi-stability regime, at a much cooler temperature than observed. Scaling factors and dependancies between , and present a complex situation to break into linear effects.
We constrain our models with HD 166734 observations by utilising a variety of plots for consistency between mass, luminosity, temperature and age, see figure 3. We explore the HRD position and compare this with the spectroscopic HRD (sHRD), which removes uncertainties with distance and luminosity. Simultaneously, we correlate ages of the primary and secondary with dynamical masses and mass-loss rates. Figure 3 illustrates the relevant plots in which HD 166734 was compared to in order to concurrently reproduce the observed masses and luminosities.
Maeder (1986) discusses the complexity of mixing processes that apply to stellar evolution and the disentanglement required to fully understand the linear effects of each process. Mass loss is thought to dredge up fusion products to the surface while diminishing the core mass, extending the MS- lifetime at the expense of the He-burning lifetime. In this respect, stellar winds behave similarly to convective overshooting or rotational mixing, even in extreme cases where a star may evolve quasi-chemically homogeneously due to extensive mixing.
Earlier models which solely employ mass loss as the mixing process may present a simpler solution to understanding the full effects of this process. Though it has been stressed that in order to reproduce observations such as in the 34 open clusters from Maeder & Mermilliod (1981), an extended MS was required leading to conclusions that overshooting is required as an additional mechanism of mixing.
Similarly, in section 3 we have emphasised the necessity of rotational mixing in reproducing observed surface abundances. Therefore, since overshooting, mass loss and rotation have similar effects on the MS-lifetime and appearance of CNO-products at the stellar surface, a method of separating these processes must be developed.
Challenges in reproducing masses and luminosities simultaneously remained while comparing HRD’s and mass-age plots. It was consequently thought to be more insightful to compare our models by mass and luminosity directly. Interpreting behavioural characteristics of physical processes in this way has opened a diversity of information on luminosity and mass, as illustrated in Fig. 4.
Figure 4 highlights the key features in the Mass-Luminosity Plane. As the star evolves with time, the vector of mass and luminosity increases in length, since main-sequence stars increase in luminosity due to hydrogen burning. In this sense the plot is similar to the HRD where time can be interchanged with temperature, since we also follow the vector length with respect to temperature, reaching characteristics such as the bi-stability jump. Figure 4 demonstrates the evolution of a theoretical model to a particular age or temperature, by which we can compare this point with observations (e.g. an observed effective temperature).
We note that the gradient of this vector is reliant on the mass-loss rate, or in this case factors of the mass loss prescription from Vink et al. (2001). Unsurprisingly, this feature becomes more prominent with higher masses, e.g. 60 compared to that of a 20 model, (see Fig. 5). We find that the position of the vector at a given evolutionary phase or temperature can only be further extended in length by increased rotation or overshooting , since greater mixing leads to higher luminosities (see Fig. 6 and 7), which will have a higher mass-loss rate and subsequently a lower mass.
When analysing our grid of models for the mass range 8-60 we found a set of features in the plane which provide fundamental boundaries to stellar evolutionary models. Figure 4 illustrates one of these boundaries by a red solid ’forbidden’ region, by which the mass-luminosity relation (see Eq. 2) sets the initial mass and luminosity. As a result of this relationship, stellar evolution models cannot lie within the red ’forbidden’ region.
Similarly, if the length of the vector in the plane increases not only with time, but also temperature (as in the HRD), then we can adjust the length of our model based on an observed temperature. Thus we set an initial position and a final position in the plane for our models based on observed stellar parameters such as mass, luminosity and temperature. We can then utilise these positions to better understand processes such as rotation, mass loss and overshooting, since these all have an affect on our now ’measured’ vector length.
Figures 5, 6 and 7 each illustrate a process which influences the length or gradient of our vector in the plane. Figure 5 demonstrates that the mass loss rate will dictate a steep or shallow gradient, which again must be reached with the initial and final positions determined by the boundaries shown (i.e. the black line representing the forbidden region, and observations illustrating the final position). Figure 6 shows the possibility of extending the length of the vector by increasing the initial rotation rate, hence enhancing the luminosity. Finally, we can further extend the vector length by overshooting as represented in figure 7 if rotation can be constrained through other methods such as and surface enrichments.
The range of explored factors of Vink et al. (2001) mass-loss prescription can be seen in Fig. 5 for models with initial masses 40 and 60 . As we would expect, the factor of mass-loss rate has a much larger effect at 60 than the 40 . We find that due to the ’forbidden’ region highlighted in Fig. 4, the gradients of models with 2-3 times the Vink et al. (2001) prescription are much too shallow to reach observed initial luminosities of a 60 star for example.
Fig. 6 illustrates an increase in luminosity by 0.1 dex for an increase in rotation of 200-400 km s. We find models with initial rotation rates of 100 km sand 200 km sare indistinguishable in the plane, though a notable increase in luminosity occurs above 200 km s. We use the TAMS here as a reference point (blue dots) for each model demonstrating the effects of increased mixing by rotation or overshooting.
5 Observational Constraints
5.1 HD 166734 Parameter Space
To determine the initial parameters of the system HD 166734, we have computed a collection of models which adapt our methods from sections 3 and 4, for a variety of initial masses, mass-loss rates, and rotation rates. Due to constraining observations we have reproduced dynamical masses, luminosities, and surface nitrogen abundances based on a selection of parameters.
Since there are multiple solutions to the current evolutionary stage, we present a parameter space in which the system can be reproduced within observational errors. This is necessary as following models with increased rotation or overshooting will lead to higher luminosities. For example, models with higher mass-loss rates will require lower initial masses and thus lower rotation rates.
We can reject extreme factors of Vink et al. (2001) mass-loss rates due to the initial mass boundary in figure 4, such that we can reproduce the system with factors 0.5 - 1.5 of the Vink et al. (2001) recipe. For initial masses of 55 - 60 for the primary and 42 - 47 for the secondary, we find a range of relevant overshooting parameters of 0.3 - 0.5 and 0.1 - 0.4 for the primary and secondary respectively. We also stress that when calibrating our theoretical models, we ensure that the factor of mass loss recipe (Vink et al. 2001) remains constant between the two objects to reach the most reliable solution.
When fixing the mass-loss prescription to a factor unity of Vink et al. (2001), we predict initial masses of 55 and 45 for the primary and secondary respectively. Initial rotation rates have been selected such that observed surface N abundances are reproduced, with 250 km sand 120 km sfor the primary and secondary respectively. Having fixed the mass-loss rate and rotation rates of our models, we utilise the plane to measure the necessary overshooting required to reach the observed mass, luminosity and effective temperature of the primary and secondary respectively. We discover greater values of required to reproduce these stellar parameters with the primary adopting = 0.3 0.1 and the secondary requiring extra mixing of = 0.5 0.1 in order to reach the observed luminosity.
5.2 Applications from analysis
We can further constrain the evolution of HD 166734 due to constraints provided by and [N/H] abundances, after we constrained in the plane. We seem to require a larger amount of core overshooting required for the secondary star than for the more luminous primary. As the primary initial mass is of the order of 60 , effects such as envelope inflation (Gräfener et al. 2012; Sanyal et al. 2015) and mass loss may potentially effect the primary’s stellar radius to a larger extent than it would for the secondary. Therefore, instead of arguing for an inverse mass dependence of the , we remain conservative, and consider the determination of the secondary star as more secure than that of the primary.
5.3 Galactic observational sample
We aimed to consolidate our results from sections 3-5 by overlaying our calibration models for HD 166734 with a sample of 30 Galactic O stars from Markova et al. (2018) to ensure our calibration is representative of a larger sample, and not unique to our selected testbed HD 166734 only. Markova’s analysis provided photospheric and wind parameters, including rotation rates, and surface N abundances by applying the model atmosphere code FASTWIND (Puls et al. 2005) to optical spectroscopy. Table 3 provides the key parameters we explored. We compared these Galactic data to our grid of models with the aim to constrain treatments of rotation and convection and contrast this with treatments from other evolutionary codes.
|HD/CPD||[kK]||log(L/ )||[ ]|
|HD 64568a||48.0 1.5||5.80||48.5 17.9||8.18|
|HD 46223||43.51.5||5.58||38.9 14.4||8.58|
|HD 93204a||40.51.0||5.70||60.9 22.5||7.78|
|CPD-59 2600a||40.0 1.0||5.40||40.3 14.9||7.78|
|HD 93843a||39.0 1.5||5.91||64.123.8||7.98|
|HD 91572a||38.5 1.0||5.35||32.712.1||8.37|
|HD 93222||38.0 1.0||5.36||35.213.0||7.98|
|HD 97848||36.5 1.0||5.03||19.67.2||8.38|
Over-plotting our models with the Galactic sample from Markova et al. (2018) not only allowed evolutionary masses, as derived from the spectroscopic HRD, to be contrasted with the masses derived from the standard HRD, but also allowed a comparison to be made between our grids of models and the Brott et al. (2011a) -like tracks. For a sample of representative masses 20 , 40 and 60 of our grid, and the applied parameters from Brott et al. (2011a) models (see Table 4) for , and in order to test our new prescriptions discussed in sections 3-4, denoted here as Brott-like models. Note that we make comparisons based on two rotation rates for both Brott-like models and our new grid.
In figures 10 and 11, we present our tracks in blue and grey for 100 km sand 250 km srespectively, and Brott et al. (2011a) -like tracks in red and pink with 100 km sand 250 km srespectively. We find a diminished discrepancy between with our new models when compared to that of Brott et al. (2011a) parameters by approximately 0.1dex as a result of increased luminosities with increased , as well as the absence of the . This discrepancy is noted in Markova et al. (2018, pg. 12) as a systematic difference in Ekström et al. (2012) models whereas Brott et al. (2011a) models appear 10-20% less massive in HRD’s compared to sHRD’s for masses above 30 .
While exploring the possibility of a reduced discrepancy between evolutionary masses derived from the HRD and sHRD, we compared luminosities deduced by Markova et al. (2018) with the recent Gaia DR2 distance estimates, Gaia Collaboration et al. (2018). We found discrepancies in the distances of our sample and therefore luminosities when using the newly calculated distances through Bailer-Jones et al. (2018). However, when considering potential errors due to reddening, we found the reddening error to be larger than the already substantial error in the Gaia distances. Therefore, final answers would require a new spectroscopic analysis with proper consideration for reddening parameters and Gaia DR2 distances, which lies beyond the scope of this study.
6 Grid analysis
Our systematic grid has been completed for two extreme values of 0.1 and 0.5, since analysis from HD 166734 invokes an argument for increased overshooting of 0.5, and 0.1 for allowing comparisons with previously published grids such as Ekström et al. (2012) models. For this reason, we show our tracks in figures 10 and 11, computed with our larger prescription of 0.5. Table 1. highlights the parameter space in which we compose our grid, compared with models from Ekström et al. (2012) and Brott et al. (2011a) in table 4. We identify the key variances as extra mixing by overshooting of up to 0.5, as well as decreased mass-loss by excluding rotationally-induced mass-loss. In order for our results to have relevance beyond the MS, we must ensure that observations of later evolutionary phases can be matched with our parameters.
|Code||Brott et al. (2011a) Parameters||Ekström et al. (2012) Parameters||This work : MESA|
|[ ]||5 - 60||8 - 120||8, 12, 16, 20, 25, 30, 35, 40, 45, 50, 55, 60|
|[ km s]||0 - 550||0.4||0, 100, 200, 300, 400, 500|
6.1 RSG upper luminosity limit
Red supergiants (RSGs) have been observed at luminosities up to log L/ 5.5 - 5.8 (Levesque et al. 2005; Humphreys & Davidson 1994), with an observed cut-off after which RSGs are not created. In analysing our set of models, we compare final stages of evolution to RSG, BSG, or possible WR-type evolution (bluewards). We find that having fixed our mass-loss rates for clarity, overshooting has the dominant effect on the maximum mass/luminosities at which RSGs are formed. Note that the treatment of convection in the outer layer as well as the adopted mass-loss regime for this evolutionary phase will also affect the position of the RSG.
If both and are fixed at lower values such as 0.1, with standard Vink et al. (2001) mass-loss rates, then RSGs are formed at luminosities of up to log L / = 6.0, even for masses up to 60 . Since observations of RSGs suggest a lower cut-off in the range of 5.5-5.8 dex, a higher value of is required to match observations. Models which have adopted 0.5 remain blue above log L/ 5.8 without evolving to RSGs, in agreement with the Humphreys-Davidson limit.
We also examine final evolutionary phases for models from section 4 with factors of between 0.5 and 1.5, for 0.1 and 0.5, finding that regardless of mass-loss rate (within our accepted parameter range), models which adopt of 0.1 result in RSG evolution even at 60 . Figure 12. represents the observed luminosity cut-off for RSG evolution when implementing an enhanced overshooting of 0.5.
6.2 Compactness parameter
Enhancing the overshooting parameter has repercussions for the final fate of our models. The consequences of our grid results impact final mass estimates as well as compactness parameters which may dictate black hole and neutron star formation. We include our estimates of the compactness parameter for all final models, in which post-MS evolution is set as standard in MESA, as a function of and initial rotation rate. O’Connor & Ott (2011) quantifies the compactness of a presupernova stellar core as seen in Eq. (3) where has been selected as the relevant mass within which the iron core density gradient may be defined. The parameter thus denotes how easily a presupernova stellar core will explode, with a low value leading to a more likely solution where the star will explore rather than collapse to form a black hole. Sukhbold & Woosley (2014) find dependancies in the treatment of convection, including overshooting, with the ”explodability” of presupernova models computed with MESA.
We note that with an extended convective core ( 0.5), and thus MS- lifetime, it is more difficult to form black holes than for 0.1 at the same mass range, as shown in Fig. 13. A clear correlation with rotation rate is not reached, but can be compared via the representative colour-bar. However, we note that rapidly rotating models with 500 km s have a very low and may explode more easily, regardless of .
Renzo et al. (2017) present values of 0.25 for the mass range 15-30 with varying wind parameters, which are analogous to our results for a similar mass range with 0.5. However we reach values of 0.9 with 0.1 for the mass range 20-30 , suggesting this level of convective mixing may be less desirable for a high chance of explodability.
7 Discussion and Conclusion
We presented a calibrated grid of evolutionary models for the mass range 8 - 60 , with initial rotation rates 0 - 500 km s, and two values of overshooting 0.1 and 0.5, (Higgins & Vink 2018). These models have been constrained based on results of our testbed eclipsing binary HD 166734. We found that rotational mixing is necessary for reproducing the observed intermediate nitrogen enrichments, after first having explored the possibility that this could be achieved by overshooting and mass loss alone. In particular, we developed a method of reproducing the eclipsing binary HD 166734 based on the fundamental properties of mass and luminosity, utilising a new tool: the mass-luminosity plane, the plane. It presents extensive information about the dominant physical processes for various mass ranges.
- First of all, the plane allowed us to exclude very large increases or reductions in the standard mass-loss rates, via the gradient in the plane. More specifically, we can exclude mass-loss factors that lie beyond 0.5 - 1.5 times the Vink et al. (2001) prescriptions.
- Secondly, the extension of the data in the plane forced us to conclude that an additional process is required. Therefore, we favour large overshooting values of order 0.5. The reproduced evolution of our testbed high-mass binary HD 166734 required this enhanced mixing by rotation and overshooting in order to increase the luminosity to that of the observed primary and secondary luminosities.
- Rotational mixing proves intrinsically necessary as the process whereby nitrogen is dredged to the surface in any intermediate quantity. Though the process has been widely researched in the last few decades, the importance of reproducing observed surface abundances such as in HD 166734 has not been sufficiently emphasised. We confirm that alternative mechanisms such as convection and mass loss cannot alone reproduce observed surface nitrogen abundances.
- Finally, we disfavour the application of rotationally-induced mass loss, in agreement with results from 2D computations of Müller & Vink (2014), as interacting processes artificially altering the initial mass-loss rate, leads to an entangled set of processes which cannot be separately constrained. The evolution of HD 166734 cannot be reproduced with the inclusion of this theory.
If we compare observations to our new prescriptions of overshooting, mass loss, and rotation, we now open the possibility of an extended MS width, reinforcing the argument of B supergiants still being core H-burning objects (Vink et al. 2010).
- Bailer-Jones et al. (2018) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58
- Bertelli et al. (2009) Bertelli, G., Nasi, E., Girardi, L., & Marigo, P. 2009, A&A, 508, 355
- Bonanos et al. (2004) Bonanos, A. Z., Stanek, K. Z., Udalski, A., et al. 2004, ApJ, 611, L33
- Brott et al. (2011a) Brott, I., de Mink, S. E., Cantiello, M., et al. 2011a, A&A, 530, A115
- Brott et al. (2011b) Brott, I., Evans, C. J., Hunter, I., et al. 2011b, A&A, 530, A116
- Castro et al. (2014) Castro, N., Fossati, L., Langer, N., et al. 2014, A&A, 570, L13
- Chieffi & Limongi (2013) Chieffi, A. & Limongi, M. 2013, ApJ, 764, 21
- Chiosi & Maeder (1986) Chiosi, C. & Maeder, A. 1986, ARA&A, 24, 329
- Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102
- Claret & Torres (2017) Claret, A. & Torres, G. 2017, ApJ, 849, 18
- de Mink et al. (2009) de Mink, S. E., Cantiello, M., Langer, N., & Pols, O. R. 2009, Communications in Asteroseismology, 158, 94
- Ekström et al. (2012) Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146
- Evans et al. (2008) Evans, C., Hunter, I., Smartt, S., et al. 2008, The Messenger, 131, 25
- Evans et al. (2005) Evans, C. J., Smartt, S. J., Lee, J.-K., et al. 2005, A&A, 437, 467
- Friend & Abbott (1986) Friend, D. B. & Abbott, D. C. 1986, ApJ, 311, 701
- Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, ArXiv e-prints [\eprint[arXiv]1804.09365]
- Garmany & Fitzpatrick (1988) Garmany, C. D. & Fitzpatrick, E. L. 1988, ApJ, 332, 711
- Gies (2003) Gies, D. R. 2003, in IAU Symposium, Vol. 212, A Massive Star Odyssey: From Main Sequence to Supernova, ed. K. van der Hucht, A. Herrero, & C. Esteban, 91
- Gräfener et al. (2012) Gräfener, G., Owocki, S. P., & Vink, J. S. 2012, A&A, 538, A40
- Grevesse & Sauval (1998) Grevesse, N. & Sauval, A. J. 1998, Space Sci. Rev., 85, 161
- Grin et al. (2017) Grin, N. J., Ramírez-Agudelo, O. H., de Koter, A., et al. 2017, A&A, 600, A82
- Groh et al. (2014) Groh, J. H., Meynet, G., Ekström, S., & Georgy, C. 2014, A&A, 564, A30
- Grunhut et al. (2012) Grunhut, J. H., Wade, G. A., & MiMeS Collaboration. 2012, in Astronomical Society of the Pacific Conference Series, Vol. 465, Proceedings of a Scientific Meeting in Honor of Anthony F. J. Moffat, ed. L. Drissen, C. Robert, N. St-Louis, & A. F. J. Moffat, 42
- Heger et al. (2000) Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368
- Herrero et al. (1992) Herrero, A., Kudritzki, R. P., Vilchez, J. M., et al. 1992, A&A, 261, 209
- Higgins & Vink (2018) Higgins, E. R. & Vink, J. S. 2018, ArXiv e-prints [\eprint[arXiv]1810.12924]
- Humphreys & Davidson (1994) Humphreys, R. M. & Davidson, K. 1994, PASP, 106, 1025
- Hunter et al. (2008) Hunter, I., Brott, I., Lennon, D. J., et al. 2008, ApJ, 676, L29
- Keszthelyi et al. (2017) Keszthelyi, Z., Puls, J., & Wade, G. A. 2017, A&A, 598, A4
- Langer (2012) Langer, N. 2012, ARA&A, 50, 107
- Levesque et al. (2005) Levesque, E. M., Massey, P., Olsen, K. A. G., et al. 2005, ApJ, 628, 973
- Maeder (1983) Maeder, A. 1983, A&A, 120, 113
- Maeder (1986) Maeder, A. 1986, in IAU Symposium, Vol. 116, Luminous Stars and Associations in Galaxies, ed. C. W. H. De Loore, A. J. Willis, & P. Laskarides, 287–299
- Maeder (2000) Maeder, A. 2000, New A Rev., 44, 291
- Maeder & Mermilliod (1981) Maeder, A. & Mermilliod, J. C. 1981, A&A, 93, 136
- Maeder & Meynet (1987) Maeder, A. & Meynet, G. 1987, A&A, 182, 243
- Maeder & Meynet (1988) Maeder, A. & Meynet, G. 1988, A&AS, 76, 411
- Maeder & Meynet (1991) Maeder, A. & Meynet, G. 1991, A&AS, 89, 451
- Maeder & Meynet (2000) Maeder, A. & Meynet, G. 2000, A&A, 361, 159
- Maeder & Meynet (2015) Maeder, A. & Meynet, G. 2015, in IAU Symposium, Vol. 307, New Windows on Massive Stars, ed. G. Meynet, C. Georgy, J. Groh, & P. Stee, 9–19
- Mahy et al. (2017) Mahy, L., Damerdji, Y., Gosset, E., et al. 2017, A&A, 607, A96
- Markova et al. (2018) Markova, N., Puls, J., & Langer, N. 2018, ArXiv e-prints [\eprint[arXiv]1803.03410]
- Martins & Palacios (2013) Martins, F. & Palacios, A. 2013, A&A, 560, A16
- Meynet et al. (1994) Meynet, G., Maeder, A., Schaller, G., Schaerer, D., & Charbonnel, C. 1994, A&AS, 103
- Müller & Vink (2014) Müller, P. E. & Vink, J. S. 2014, A&A, 564, A57
- O’Connor & Ott (2011) O’Connor, E. & Ott, C. D. 2011, ApJ, 730, 70
- Pavlovski et al. (2018) Pavlovski, K., Southworth, J., & Tamajo, E. 2018, MNRAS, 481, 3129
- Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3
- Pfuhl et al. (2014) Pfuhl, O., Alexander, T., Gillessen, S., et al. 2014, ApJ, 782, 101
- Puls et al. (2005) Puls, J., Urbaneja, M. A., Venero, R., et al. 2005, A&A, 435, 669
- Renzo et al. (2017) Renzo, M., Ott, C. D., Shore, S. N., & de Mink, S. E. 2017, A&A, 603, A118
- Sanyal et al. (2015) Sanyal, D., Grassitelli, L., Langer, N., & Bestenlehner, J. M. 2015, A&A, 580, A20
- Schaller et al. (1992) Schaller, G., Schaerer, D., Meynet, G., & Maeder, A. 1992, A&AS, 96, 269
- Simón-Díaz & Herrero (2014) Simón-Díaz, S. & Herrero, A. 2014, A&A, 562, A135
- Southworth et al. (2004) Southworth, J., Maxted, P. F. L., & Smalley, B. 2004, MNRAS, 351, 1277
- Sukhbold & Woosley (2014) Sukhbold, T. & Woosley, S. E. 2014, ApJ, 783, 10
- Tkachenko et al. (2014) Tkachenko, A., Degroote, P., Aerts, C., et al. 2014, MNRAS, 438, 3093
- Vink (2015) Vink, J. S. 2015, in Astrophysics and Space Science Library, Vol. 412, Very Massive Stars in the Local Universe, ed. J. S. Vink, 77
- Vink et al. (2010) Vink, J. S., Brott, I., Gräfener, G., et al. 2010, A&A, 512, L7
- Vink et al. (2001) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574
- Vink & Gräfener (2012) Vink, J. S. & Gräfener, G. 2012, ApJ, 751, L34
- Weidner & Vink (2010) Weidner, C. & Vink, J. S. 2010, A&A, 524, A98