© ESO 2019

1. Introduction

Massive stars with an initial mass above 8 M_⊙ have a diversity of possible evolutionary channels that are dictated by the dominant processes acting on their structure. The extent of these dependencies 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 and 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 M_⊙ mass loss completely dominates the evolution of O-type stars (e.g. Vink & Gräfener 2012; Vink 2015), whilst in the range 30 M_⊙ <  M <  60 M_⊙ mass loss is one of the important ingredients (e.g. Langer 2012; Groh et al. 2014). At lower masses, (i.e. below ∼30 M_⊙) the evolution is thought to be heavily influenced by rotation (e.g. Maeder & Meynet 2000). Over the explored mass range within this paper, i.e. 8–60 M_⊙, we will consider the effects of mass loss and rotation, as well as convective overshooting, which may all play a role in the evolution of these stars.

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 α_ov which we explore in this study corresponds to the fraction of the pressure scale height H_p by which particles continue to travel a distance l_ov beyond the convective core boundary. This form of mixing has been explored for decades, with few constraints on its size (α_ov) in the high-mass range. It has been argued as essential for reproducing observations, although evidence is lacking for dependencies such as mass (Claret & Torres 2017). Another process that may potentially affect stellar evolution is the presence of a magnetic field, however, Grunhut et al. (2012) showed 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 dependencies 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 (α_ov) 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 have a degree of uncertainty, yet to establish a clear comparison between codes, it would be beneficial to examine physical implementations with one and the same code. Therefore, in this work we aim to compute massive star models with both new and existing prescriptions using the same evolutionary code, Modules for Experiments in Stellar Astrophysics (MESA; 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 that these objects are MS, core H-burning stars. The existence of a large number of slow-rotating B supergiants however, (with υ sini ≲ 50 km s⁻¹) is suggestive of an evolved star that has completed the MS phase and been spun down.

Vink et al. (2010) also considered 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) addressed that a higher value of α_ov would result in a lower critical mass at which BSB would be efficient. Test models show that BSB occurs in present models with α_ov = 0.335 above a critical mass of 35 M_⊙ in the Large Magellanic Cloud (LMC), yet the critical mass drops to 20 M_⊙ for the same metallicity with an increase of overshooting to α_ov = 0.5 (Vink et al. 2010).

The determination of α_ov 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 υ sini and log g (Brott et al. 2011a). Many other estimations of α_ov 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 approaches 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). However, the Galactic data from this study might be biased compared with the unbiased LMC data from the VLT-Flames Surveys of massive stars Evans et al. (2005, 2011) as they do not show a gap between O and B supergiants (Vink et al. 2010), thereby 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 α_ov = 0.1 and 0.5 to illustrate both lower adopted values and 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 MS (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) presented 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. Although 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 Tka for detached eclipsing binaries, however these works utilised lower mass stars (up to ∼15 M_⊙) which did not incorporate the interacting effects of mass loss, overshooting and rotation as we do in this study.

We use constraints relative to α_ov and Ṁ to investigate the possible evolutionary paths of a high mass, detached binary, HD 166734, modelled in this work as a test-bed for single star evolution. As previously mentioned, dominant processes take effect at varying mass ranges, yet with dynamical masses of 39.5 M_⊙ and 33.5 M_⊙, 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 Sect. 2.1 and a calibration of mixing processes in Sect. 3. We explore a new tool for comparing observations with models in the mass-luminosity plane in Sect. 4 and we provide our final results for HD 166734 in Sect. 5. We present our grid of models alongside a sample of Galactic O-stars in Sect. 6, and provide further results in Appendix A. Finally, we highlight our conclusions in Sect. 7; the remaining full grid plots are available in Appendix A.

2. Methodology

2.1. MESA: Treatment of convection, mass loss, and rotation

A set of evolutionary models was calculated for massive MS stars with the 1D, stellar evolution code MESA, for example Paxton et al. (2011), as a comparison for both the primary and secondary of HD 166734 (see Sect. 2.2). 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 arose in which computations were concluded earlier. We adopt the default metallicity in MESA of Z = 0.02 with the chemical mixture from Grevesse & Sauval (1998) to provide direct comparisons with chemical abundances in Galactic observations.

Convection is treated by the mixing length theory where α_MLT = 1.5, with a semi-convection efficiency parameter of α_semi = 1. The convective core boundary is defined by the Ledoux criterion¹, in which overshooting succeeds convective mixing at the core boundary, increasing the temperature gradient ∇_T by implementing a thermal gradient ∇_rad (e.g. Choi et al. 2016). This method of extending the core is denoted as step-overshooting, which enhances the core by a factor α_ov of the pressure scale height H_p. Experiments in the dependencies of this parameter are completed in the following sections.

We then compared our grid with treatments of α_ov 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) presented a calibration of the overshooting parameter by comparing the TAMS of 16 M_⊙ 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 υ sini beyond which a large number of slow-rotating B supergiants are located, assumed to be post-MS objects. The model with α_ov = 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 α_ov = 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 M_⊙ where convective mixing plays a dominant role compared to that of rotational mixing. Hence this allowed for 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. We tested various factors of this mass-loss regime to determine the possibility of extreme rates. We hence explored a range of multiplication factors of Vink et al. (2001) mass-loss rates from 0.1 to 10 times the standard prescription. We later applied rotation in our models through a fully diffusive approach with appropriate instabilities such as the Eddington-Sweet circulation, dynamical and secular shear instabilities. We also considered 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 as physically different in another code, then the results would differ quantitatively, but qualitatively have the same behaviour.

A systematic grid of models was calculated for comparison with a larger sample, including new prescriptions discussed in Sects. 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 α_ov to explore the extent of extra mixing. We evolved each model to core collapse, unless convergence problems highlighted unlikely solutions. For this purpose, Vink et al. (2001) provided 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 test-bed 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 two 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 M_⊙. Observations of high-mass eclipsing binaries are sparse, and even more extreme for detached, non-interacting stars that may be treated as evolved single stars. Since observations of massive single stars may sometimes highlight discrepancies between spectroscopic and evolutionary masses. In this case, we have an ideal opportunity because 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 M_⊙ are extremely rare (see e.g. Bonanos et al. 2004; de Mink et al. 2009; Pfuhl et al. 2014; Gies et al. 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 υ sini for both components may at face value be considered of interest in terms of synchronisation, but Mahy et al. (2017) have argued against synchronisation because the rotation speeds are lower than the orbital period. In addition, we note that the υ sini 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 υ sini are truly the result of rotation (see Simón-Díaz & Herrero 2014). We thus treat the υ sini values as upper limits, and we consider the similar values of the two components as merely a coincidence.

We utilise this agreement between dynamical and spectroscopic masses, allowing HD 166734 to be treated as an excellent test-bed 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 υ sini 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 and determining the 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) showed 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 that 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 M_⊙ and 47.4 M_⊙ 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 α_ov. We initially attempted 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. We found however that these models do not offer solutions in which 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 that 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 dependency 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 ¹⁴N 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 α_ov as the sole mechanisms for chemical mixing. The importance of convective overshooting has been stressed in these early publications as α_ov leads to a range of stellar ages, due to the dependence of T_eff at TAMS on α_ov, (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) presented 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 dependency 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 α_ov of 0.1–0.8. We also note that models with increased overshooting result in earlier enrichment by up to 1 Myr, regardless of mass-loss rates. In Fig. 1 we present the nitrogen enrichments for a sample of models of primary and secondary masses.

Fig. 1. Surface nitrogen abundances as a function of stellar age for extreme values of αov and M⊙. The blue lines represent rotating 40 M⊙ models with an initial rotation rate of 200 km s−1, αov = 0.1 (dash-dotted), and αov = 0.5 (solid). The red lines show the corresponding non-rotating models for the same mass and values of αov respectively.

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 results 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.

Figure 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 α_ov from 0.1 to 0.5, we get much larger surface enrichments that 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 discovered a problem with respect to interacting processes such as rotation and mass loss, which consequently have a non-linear effect 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 the 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. We now reached the dynamical masses; this also led to a significant drop in luminosity, which correlates to a shallow gradient in the M−L 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 1D radiation-driven wind models of Friend & Abbott (1986), who proposed an equatorially enhanced stellar wind and an increased mass-loss rate due to a lower effective gravity at the equator. This result 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), i.e.

(1)

We note that the Geneva group (e.g. Maeder et al. 2015; Ekström et al. 2012) employed 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 and 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; for example 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 M_⊙, 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.

Fig. 2. Mass-loss rate as a function of stellar age for comparison of rotationally induced mass loss. The red solid lines represent models with default MESA settings of rotationally induced Ṁ (see Eq. (1)), with an initial mass of 40 M⊙ and rotation rates of 100 km s−1 and 300 km s−1. All other processes have been set to default values to avoid conflict in our analysis. The dashed blue lines show the corresponding models with ξ = 0, for the same mass and rotation rates respectively.

4. Mass–luminosity plane

When comparing models in Sect. 3 we find that enhanced mass-loss regimes lead to unrealistic luminosities that 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 dependencies of these properties,

(2)

where α varies as a function of mass and μ is the mean molecular weight.

The most fundamental characteristics of the evolution of a star 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 α_ov. Yet we may constrain these solutions by 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.

Equation (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 L_init by ≈30% or 0.11 dex. This offers a potential scenario that would allow for a higher luminosity and lower M_init, nonetheless 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_eff’s, we note that both objects remain too hot to be post-MS objects, regardless of α_ov, thus excluding later evolutionary phases as a viable solution.

We observe some models that 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_eff since the dynamical mass is only reached during the bi-stability regime at a much cooler temperature than observed. Scaling factors and dependencies between Ṁ, α_ov and υ sini 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 Fig. 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 comparing HD 166734 to reproduce the observed masses and luminosities concurrently.

Fig. 3. Evolution of a 40 M⊙ model with an initial rotation rate of 100 km s−1 and αov = 0.1. This is shown for a variety of conventional plots such as the luminosity and mass as a function of stellar age (left upper and lower), as well as in a standard HRD (right upper), and finally mass as a function of effective temperature.

Maeder et al. (1986) discussed the complexity of mixing processes that apply to stellar evolution and the disentanglement required to understand the linear effects of each process fully. 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. Although it has been stressed that 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 Sect. 3 we 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 the HRD 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.

Fig. 4. Illustration of the mass-luminosity plane with a typical evolutionary track entering the ZAMS at the red dot, evolving along the black arrow. The dotted vector suggests how increased rotation and/or convective overshooting may extend the M−L vector. The curved dashed line represents the gradient at which mass-loss rates affect this M−L vector. The red solid region represents the boundary set by the mass-luminosity relationship, and as such is forbidden.

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 MS stars increase in luminosity due to hydrogen burning. In this sense the M−L plot is similar to the HRD in that 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, for example 60 M_⊙ compared to that of a 20 M_⊙ 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 α_ov, since greater mixing leads to higher luminosities (see Figs. 6 and 7), which have a higher mass-loss rate and subsequently a lower mass.

Fig. 5. Evolution of both 40 M⊙ and 60 M⊙ models with initial rotation rates of 100 km s−1 and αov = 0.1 are shown for a variety of factors of Vink et al. (2001) mass-loss rate (0.1–3 times), demonstrating the gradients for each model represented by the green and red solid lines in the M−L plane.

Fig. 6. Evolution of a 40 M⊙ model with a factor of unity of the mass-loss prescription and αov = 0.1 shown for a variety of initial rotation rates from 100–500 km s−1. The length of the M−L vector at a given evolutionary stage can be extended via increased rotation as shown by the blue dots corresponding to TAMS for each model.

Fig. 7. Evolution of a 40 M⊙ model with an initial rotation rate of 100 km s−1 and a factor of unity for the mass-loss prescription shown for a variety of overshooting αov = 0.1, 0.3, 0.5. The blue dots correspond to TAMS for each model, demonstrating the increase in luminosity or decrease in mass for an increase in overshooting of αov = 0.1 − 0.5. This illustrates the possible further extension of a vector in this plane by extending αov.

When analysing our grid of models for the mass range 8–60 M_⊙ we found a set of features in the M−L plane that 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 M−L 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 M−L 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–7 each illustrate a process which influences the length or gradient of our vector in the M−L plane. Figure 5 demonstrates that the mass-loss rate dictates 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 Fig. 7 if rotation can be constrained through other methods such as υ sini 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 M_⊙ and 60 M_⊙. As we would expect, the factor of mass-loss rate has a much larger effect at 60 M_⊙ than the 40 M_⊙. We find that due to the “forbidden” region highlighted in Fig. 4, the gradients of models with two-three times the Vink et al. (2001) prescription are much too shallow to reach observed initial luminosities of a 60 M_⊙ star for example.

Figure 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 s⁻¹ and 200 km s⁻¹ are indistinguishable in the M−L plane, although a notable increase in luminosity occurs above 200 km s⁻¹. We use the TAMS 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 computed a collection of models that adapt our methods from Sects. 3 and 4, for a variety of initial masses, mass-loss rates, α_ov 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 leads to higher luminosities. For example, models with higher mass-loss rates requires 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 Fig. 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 M_⊙ for the primary and 42–47 M_⊙ 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 M_⊙ and 45 M_⊙ for the primary and secondary respectively. Initial rotation rates have been selected such that observed surface N abundances are reproduced with 250 km s⁻¹ and 120 km s⁻¹ for the primary and secondary, respectively. Having fixed the mass-loss rate and rotation rates of our models, we utilise the M−L plane to measure the necessary overshooting required to reach the observed mass, luminosity and effective temperature of the primary and secondary. We discover greater values of α_ov required to reproduce these stellar parameters with the primary adopting α_ov = 0.3 ± 0.1 and the secondary requiring extra mixing of α_ov = 0.5 ± 0.1 to reach the observed luminosity.

Figure 8 illustrates the evolution of the selected models which simultaneously reproduce observed luminosities and dynamical masses at an age of ∼3 Myrs. Observed nitrogen abundances are reproduced in Fig. 9, showing the observed 3:1 ratio of the primary to secondary.

Fig. 8. HD 166734 constrained models for primary and secondary in the M−L plane (left) and HRD (right).

Fig. 9. Nitrogen enrichments from models in Fig 8. Observations of HD 166734 highlight the desired (3:1) ratio.

5.2. Applications from analysis

We can further constrain the evolution of HD 166734 due to constraints provided by υ sini and [N/H] abundances, after we constrain α_ov in the M−L plane. We seem to require a larger amount of core overshooting for the secondary star than for the more luminous primary. As the primary initial mass is of the order of 60 M_⊙, effects such as envelope inflation (Gräfener et al. 2012; Sanyal et al. 2015) and mass loss may potentially effect the stellar radius of the primary to a larger extent than it would for the secondary. Therefore, instead of arguing for an inverse mass dependence of the α_ov, we remain conservative, and consider the α_ov 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 Sects. 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 test-bed HD 166734 only. The analysis by Markova et al. (2018) 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 we contrast this with treatments from other evolutionary codes.

Over-plotting our models with the Galactic sample from Markova et al. (2018) not only allowed us to contrast evolutionary masses as derived from the spectroscopic HRD with the masses derived from the standard HRD, but also allowed us to compare our model grids with the prescriptions from Brott et al. (2011a). A sample of representative masses 20 M_⊙, 40 M_⊙ and 60 M_⊙ of our grid, and the applied parameters from Brott et al. (2011a) models (see Table 4) for α_ov and Ṁ_{rot, boost} that test our new prescriptions discussed in Sects. 3 and 4, are denoted in this work as Brott-like models. Note that we make comparisons based on two rotation rates for both Brott-like models and our new grid.

In Figs. 10 and 11, we present our tracks in blue and grey for 100 km s⁻¹ and 250 km s⁻¹, respectively, and Brott et al. (2011a) – like tracks in red and pink with 100 km s⁻¹ and 250 km s⁻¹, respectively. We find a diminished discrepancy between M_{evol, sHRD}M_{evol, HRD} with our new models when compared to that of Brott et al. (2011a) parameters by approximately 0.1 dex as a result of increased luminosities with increased α_ov, as well as the absence of the Ṁ_{rot, boost}. 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 the HRD compared to the sHRD for masses above 30 M_⊙.

Fig. 10. Comparison of parameters taken from Brott et al. (2011a) in red/pink with our new prescription in blue/grey, contrasted alongside our adopted Galactic sample of O stars from Markova et al. (2018).

Fig. 11. Comparison of models as described in Fig. 10. Data from Markova et al. (2018) are shown as black triangles in the form of a sHRD.

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 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 α_ov = 0.1 and 0.5, since analysis from HD 166734 invokes an argument for increased overshooting of α_ov = 0.5, and α_ov = 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 Figs. 10 and 11, which are computed with our larger α_ov prescription of α_ov = 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 α_ov = 0.5, and 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.

6.1. Red supergiant upper luminosity limit

Red supergiants (RSGs) have been observed at luminosities up to log L/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. We note that the treatment of convection in the outer layer and the adopted mass-loss regime for this evolutionary phase also affects the position of the RSG.

If both α_ov and Ṁ are fixed at lower values such as α_ov = 0.1, with standard Vink et al. (2001) mass-loss rates, then RSGs are formed at luminosities of up to log L/L_⊙ = 6.0, even for masses up to 60 M_⊙. Since observations of RSGs suggest a lower cut-off in the range ≈5.5–5.8 dex, a higher value of α_ov is required to match observations. Models which have adopted α_ov = 0.5 remain blue above log L/L_⊙ = 5.8 without evolving to RSGs, in agreement with the Humphreys-Davidson limit.

We also examine final evolutionary phases for models from Sect. 4 with factors of Ṁ between 0.5 and 1.5, for α_ov = 0.1 and 0.5, finding that regardless of mass-loss rate (within our accepted parameter range), models which adopt α_ov of 0.1 result in RSG evolution even at 60 M_⊙. Figure 12 represents the observed luminosity cut-off for RSG evolution when implementing an enhanced overshooting of α_ov = 0.5.

Fig. 12. Rotating models with αov = 0.5 and a factor of 1.5 Ṁ, highlighting the post-MS evolution to the red for log L up to 5.8, but above this limit remaining in the blue at TAMS.

6.2. Compactness parameter

Enhancing the overshooting parameter α_ov 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 ζ_2.5 for all final models, in which post-MS evolution is set as standard in MESA, as a function of α_ov and initial rotation rate. O’Connor & Ott (2011) quantified the compactness of a presupernova stellar core as seen in Eq. (3), where M = 2.5 M_⊙ is selected as the relevant mass within which the iron core density gradient may be defined. The parameter ζ_2.5 thus denotes how easily a presupernova stellar core explodes; with a low value leading to a more likely solution in which the star explodes rather than collapsing to form a black hole. Sukhbold & Woosley (2014) found dependencies in the treatment of convection, including overshooting, with the explodability of presupernova models computed with MESA,

(3)

We note that with an extended convective core (α_ov = 0.5), and thus MS lifetime, it is more difficult to form black holes than for α_ov = 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 υ_init = 500 km s⁻¹ have a very low ζ_2.5 and may explode more easily, regardless of α_ov.

Fig. 13. Final value of compactness parameter for the initial mass range 20–60 M⊙ with αov = 0.1 and 0.5, rotation rates of 0–500 km s−1. The circles represent models with αov = 0.1 and the triangles represent αov = 0.5. Rotation rates are shown as a fraction of the maximum 500 km s−1 in the relevant colour bar.

Renzo et al. (2017) presented values of ζ_2.5 ≥ 0.25 for the mass range 15–30 M_⊙ with varying wind parameters, which are analogous to our results for a similar mass range with α_ov = 0.5. However we reach values of ζ_2.5 ≥ 0.9 with α_ov = 0.1 for the mass range 20–30 M_⊙, suggesting this level of convective mixing may be less desirable for a high chance of explodability.

7. Discussion and conclusions

We presented a calibrated grid of evolutionary models for the mass range 8–60 M_⊙, with initial rotation rates 0–500 km s⁻¹, and two values of overshooting α_ov = 0.1 and α_ov = 0.5 (Higgins & Vink 2018). These models have been constrained based on results of our test-bed eclipsing binary HD 166734. We found that rotational mixing is necessary to reproduce 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 known as the mass-luminosity plane, the M−L plane. This tool presents extensive information about the dominant physical processes for various mass ranges.

• First of all, the M−L plane allows us to exclude very large increases or reductions in the standard mass-loss rates, via the gradient in the M−L 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 M−L plane forces us to conclude that an additional process is required. Therefore, we favour large overshooting values of order α_ov = 0.5. The reproduced evolution of our test-bed high-mass binary HD 166734 required this enhanced mixing by rotation and overshooting 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. Even 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 that 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).

Acknowledgments

We thank Andreas Sander and the referee for constructive comments, and we acknowledge the MESA developers for facilitating the use of their publicly available code.

References

Appendix A: Calibrated grid of models

Fig. A.1. Grid results for the non-rotating models in the mass range 8–60 M⊙ employing αov = 0.1 (solid lines) and αov = 0.5 (dashed lines).

Fig. A.2. Grid results for the rotating models with initial rotation rates of 100 km s−1 for the mass range as mentioned in Fig. A.2, αov = 0.1 (solid lines) and αov = 0.5 (dashed lines).

Fig. A.3. Grid results for the rotating models with initial rotation rates of 200 km s−1 for the mass range 8–60 M⊙ employing αov = 0.1 (solid lines) and αov = 0.5 (dashed lines).

Fig. A.4. Grid results for the rotating models with initial rotation rates of 300 km s−1 for the mass range as mentioned in Fig. A.2, αov = 0.1 (solid lines) and αov = 0.5 (dashed lines).

Fig. A.5. Grid results for the rotating models with initial rotation rates of 400 km s−1 for the mass range as mentioned in Fig. A.2, αov = 0.1 (solid lines) and αov = 0.5 (dashed lines).

Fig. A.6. Grid results for the rotating models with initial rotation rates of 500 km s−1 for the mass range as mentioned in Fig. A.2, αov = 0.1 (solid lines) and αov = 0.5 (dashed lines).

All Tables

All Figures

	Fig. 1. Surface nitrogen abundances as a function of stellar age for extreme values of αov and M⊙. The blue lines represent rotating 40 M⊙ models with an initial rotation rate of 200 km s−1, αov = 0.1 (dash-dotted), and αov = 0.5 (solid). The red lines show the corresponding non-rotating models for the same mass and values of αov respectively.
In the text

Fig. 2. Mass-loss rate as a function of stellar age for comparison of rotationally induced mass loss. The red solid lines represent models with default MESA settings of rotationally induced Ṁ (see Eq. (1)), with an initial mass of 40 M⊙ and rotation rates of 100 km s−1 and 300 km s−1. All other processes have been set to default values to avoid conflict in our analysis. The dashed blue lines show the corresponding models with ξ = 0, for the same mass and rotation rates respectively.

In the text

	Fig. 3. Evolution of a 40 M⊙ model with an initial rotation rate of 100 km s−1 and αov = 0.1. This is shown for a variety of conventional plots such as the luminosity and mass as a function of stellar age (left upper and lower), as well as in a standard HRD (right upper), and finally mass as a function of effective temperature.
In the text

Fig. 4. Illustration of the mass-luminosity plane with a typical evolutionary track entering the ZAMS at the red dot, evolving along the black arrow. The dotted vector suggests how increased rotation and/or convective overshooting may extend the M−L vector. The curved dashed line represents the gradient at which mass-loss rates affect this M−L vector. The red solid region represents the boundary set by the mass-luminosity relationship, and as such is forbidden.

In the text

	Fig. 5. Evolution of both 40 M⊙ and 60 M⊙ models with initial rotation rates of 100 km s−1 and αov = 0.1 are shown for a variety of factors of Vink et al. (2001) mass-loss rate (0.1–3 times), demonstrating the gradients for each model represented by the green and red solid lines in the M−L plane.
In the text

	Fig. 6. Evolution of a 40 M⊙ model with a factor of unity of the mass-loss prescription and αov = 0.1 shown for a variety of initial rotation rates from 100–500 km s−1. The length of the M−L vector at a given evolutionary stage can be extended via increased rotation as shown by the blue dots corresponding to TAMS for each model.
In the text

Fig. 7. Evolution of a 40 M⊙ model with an initial rotation rate of 100 km s−1 and a factor of unity for the mass-loss prescription shown for a variety of overshooting αov = 0.1, 0.3, 0.5. The blue dots correspond to TAMS for each model, demonstrating the increase in luminosity or decrease in mass for an increase in overshooting of αov = 0.1 − 0.5. This illustrates the possible further extension of a vector in this plane by extending αov.

In the text

	Fig. 8. HD 166734 constrained models for primary and secondary in the M−L plane (left) and HRD (right).
In the text

	Fig. 9. Nitrogen enrichments from models in Fig 8. Observations of HD 166734 highlight the desired (3:1) ratio.
In the text

	Fig. 10. Comparison of parameters taken from Brott et al. (2011a) in red/pink with our new prescription in blue/grey, contrasted alongside our adopted Galactic sample of O stars from Markova et al. (2018).
In the text

	Fig. 11. Comparison of models as described in Fig. 10. Data from Markova et al. (2018) are shown as black triangles in the form of a sHRD.
In the text

	Fig. 12. Rotating models with αov = 0.5 and a factor of 1.5 Ṁ, highlighting the post-MS evolution to the red for log L up to 5.8, but above this limit remaining in the blue at TAMS.
In the text

	Fig. 13. Final value of compactness parameter for the initial mass range 20–60 M⊙ with αov = 0.1 and 0.5, rotation rates of 0–500 km s−1. The circles represent models with αov = 0.1 and the triangles represent αov = 0.5. Rotation rates are shown as a fraction of the maximum 500 km s−1 in the relevant colour bar.
In the text

	Fig. A.1. Grid results for the non-rotating models in the mass range 8–60 M⊙ employing αov = 0.1 (solid lines) and αov = 0.5 (dashed lines).
In the text

	Fig. A.2. Grid results for the rotating models with initial rotation rates of 100 km s−1 for the mass range as mentioned in Fig. A.2, αov = 0.1 (solid lines) and αov = 0.5 (dashed lines).
In the text

	Fig. A.3. Grid results for the rotating models with initial rotation rates of 200 km s−1 for the mass range 8–60 M⊙ employing αov = 0.1 (solid lines) and αov = 0.5 (dashed lines).
In the text

	Fig. A.4. Grid results for the rotating models with initial rotation rates of 300 km s−1 for the mass range as mentioned in Fig. A.2, αov = 0.1 (solid lines) and αov = 0.5 (dashed lines).
In the text

	Fig. A.5. Grid results for the rotating models with initial rotation rates of 400 km s−1 for the mass range as mentioned in Fig. A.2, αov = 0.1 (solid lines) and αov = 0.5 (dashed lines).
In the text

	Fig. A.6. Grid results for the rotating models with initial rotation rates of 500 km s−1 for the mass range as mentioned in Fig. A.2, αov = 0.1 (solid lines) and αov = 0.5 (dashed lines).
In the text