The excess of cool supergiants from contemporary stellar evolution models defies the metallicity-independent Humphreys-Davidson limit
Avishai Gilkis, Tomer Shenar, Varsha Ramachandran, Adam S. Jermyn, Laurent Mahy, Lidia M. Oskinova, Iair Arcavi, Hugues Sana
MMNRAS , 1–13 (2021) Preprint 8 February 2021 Compiled using MNRAS L A TEX style file v3.0
The excess of cool supergiants from contemporary stellar evolution modelsdefies the metallicity-independent Humphreys-Davidson limit
Avishai Gilkis (cid:63) , Tomer Shenar † , Varsha Ramachandran , Adam S. Jermyn ,Laurent Mahy , Lidia M. Oskinova , Iair Arcavi , and Hugues Sana The School of Physics and Astronomy, Tel Aviv University, Tel Aviv 6997801, Israel Institute of Astrophysics, KU Leuven, Celestijnenlaan 200 D, 3001 Leuven, Belgium Institut f¨ur Physik und Astronomie, Universit¨at Potsdam, Karl-Liebknecht-Str. 24 /
25, D-14476 Potsdam, Germany Center for Computational Astrophysics, Flatiron Institute, 162 5th Ave, New York, NY 10010, USA CIFAR Azrieli Global Scholars program, CIFAR, Toronto, Canada
ABSTRACT
The Humphreys-Davidson (HD) limit empirically defines a region of high luminosities (log ( L / L (cid:12) ) (cid:38) .
5) and low e ff ectivetemperatures ( T e ff (cid:46)
20 kK) on the Hertzsprung-Russell Diagram in which hardly any supergiant stars are observed. Attemptsto explain this limit through instabilities arising in near- or super-Eddington winds have been largely unsuccessful. Usingmodern stellar evolution we aim to re-examine the HD limit, investigating the impact of enhanced mixing on massive stars. Weconstruct grids of stellar evolution models appropriate for the Small and Large Magellanic Clouds (SMC, LMC), as well asfor the Galaxy, spanning various initial rotation rates and convective overshooting parameters. Significantly enhanced mixingapparently steers stellar evolution tracks away from the region of the HD limit. To quantify the excess of over-luminous stars instellar evolution simulations we generate synthetic populations of massive stars, and make detailed comparisons with cataloguesof cool ( T e ff ≤ . ( L / L (cid:12) ) ≥ .
7) stars in the SMC and LMC. We find that adjustments to the mixingparameters can lead to agreement between the observed and simulated red supergiant populations, but for hotter supergiants thesimulations always over-predict the number of very luminous (log ( L / L (cid:12) ) ≥ .
4) stars compared to observations. The excessof luminous supergiants decreases for enhanced mixing, possibly hinting at an important role mixing has in explaining the HDlimit. Still, the HD limit remains unexplained for hotter supergiants.
Key words: stars: evolution – stars: massive
The upper-right part of the Hertzsprung-Russell diagram (HRD)features a stark absence of observed stars (Fig. 1), a phenomenontermed the Humphreys-Davidson (HD) limit (Humphreys & David-son 1979). While the luminous blue variables (LBVs) venture intothis region during outbursts, with a few exceptions, cool supergiants(CSGs), which comprise red, yellow, and blue supergiants (RSGs,YSGs, BSGs) with e ff ective temperatures T e ff (cid:46) . ( L max / L (cid:12) ) (cid:38) . M i (cid:38)
30 M (cid:12) ) arenot observed. In contrast, hundreds of main sequence progenitorswith M i (cid:38)
30 M (cid:12) , ranging all the way to ≈
150 M (cid:12) and perhapsmore, were directly observed (Crowther et al. 2010; Bestenlehneret al. 2011; Almeida et al. 2017; Shenar et al. 2017; Tehrani et al.2019; Mahy et al. 2020). This implies two possibilities: either starswith M i (cid:38)
30 M (cid:12) skip the CSG phase altogether, or they experiencethis phase very briefly, making it observationally rare.Despite various attempts, no evolutionary models are currently ca- (cid:63) [email protected] † [email protected] pable of reproducing the observed absence of stars beyond the HDlimit (Davies, Crowther & Beasor 2018; Schootemeijer & Langer2018). The HD limit has consequences not only for our understand-ing of the evolution of the progenitors of Wolf-Rayet (WR) stars andblack holes (BHs), but also for estimates of the likelihood of binaryinteraction in the upper-mass end of stars. The reliability of our pre-dictions of gravitational-wave (GW) events are thus severely limitedas long as the HD limit has not been su ffi ciently understood.It is commonly assumed that the HD limit is a result of radia-tive instability, tracing a region in which stars reach their Edding-ton luminosity and become LBVs (e.g., Lamers & Fitzpatrick 1988;Glatzel & Kiriakidis 1993). However, attempts to prove this bymeans of direct radiative transfer calculations in stellar models havehad only limited success. Ulmer & Fitzpatrick (1998) showed thatthe so-called modified Eddington limit mimics the shape of the HDlimit, but that it lies roughly one magnitude above the observed HDlimit. More recently, Sanyal et al. (2017) showed that stars with M i (cid:38)
30 M (cid:12) and a solar metallicity content ( Z = Z (cid:12) ) reach the Ed-dington limit in their interiors and undergo envelope inflation. How-ever, both these studies and others clearly predict that L max shouldgrow with decreasing Z . This is because the opacities in the inte- © a r X i v : . [ a s t r o - ph . S R ] F e b Gilkis, Shenar, Ramachandran, et al.
T*/kK
ZAMS V rot,ini ~180 km/s Ramachandran et al. 2019Hainich et al. 2015 (WR single)Shenar et al. 2016 (WR binary)Castro et al. 2018Bouret et al. 2013Hunter et al. 2008Dufton et al. 2019 Davies et al. 2018Neugent et al. 2010Kalari et al. 2018Trundle et al. 2004,05
H-D limit
He-ZAMS
SMC
WR OB BSG YSGRSG
R40(LBV)
AB8a AB8b l og ( L / L ) T*/kK
ZAMS V rot,ini ~100 km/s 6 M7 M9 M12 M15 M20 M25 M40 M60 M80 M100 M150 M Ramachandran et al. 2018a,bSchneider et al. 2018Hunter et al. 2008Hainich et al. 2014 (WR single)Shenar et al. 2019 (WR binary)Crowther et al. 2002Tramper et al. 2015 Davies et al. 2018,Neugent et al. 2012Neugent et al. 2012Urbaneja et al. 2017Humphreys et al. 2016
LMC
WR OB
BSG/LBV
YSGRSG
H-D limit
He-ZAMS l og ( L / L ) Figure 1.
HRD positions of populations of massive stars in the SMC (left panel) and LMC (right panel), based on analyses of apparently single and binary WRstars (Hainich et al. 2014, 2015; Shenar et al. 2016, 2019; Crowther et al. 2002; Tramper et al. 2015), YSGs (Neugent et al. 2010, 2012), RSGs (Davies et al.2018; Neugent et al. 2012), LBVs (Humphreys et al. 2016; Kalari et al. 2018), BSGs (Trundle et al. 2004; Trundle & Lennon 2005; Urbaneja et al. 2017), andpopulations of OB-type stars (Ramachandran et al. 2018a,b, 2019; Schneider et al. 2018; Hunter et al. 2008; Castro et al. 2018; Dufton et al. 2019; Bouret et al.2013). The typical error bar is shown in the top left corner. Evolutionary tracks (black solid lines), accounting for rotation with (cid:51) rot , init ≈
100 km s − for theLMC and ≈
180 km s − for the SMC (Brott et al. 2011; K¨ohler et al. 2015), are labeled with their initial mass. The red dashed lines mark the observational limitreported by Humphreys & Davidson (1979). The diagrams are largely complete for supergiants with T e ff (cid:46) . rior of stars are strongly correlated with the content of metals withinthem. The lower Z is, the weaker the outward radiative pressure be-comes. Hence, from this perspective, environments with a lower am-bient Z are expected to allow the stability of more massive stars, andhigher L max .This prediction does not seem to be confirmed by observations.Humphreys & Davidson (1979) originally reported a maximum RSGluminosity at around log ( L / L (cid:12) ) ≈ .
7, later revised slightly down-wards to 5 . Z -environment of the Triangu-lum galaxy (M33).This trend continues with the Small and Large Magellanic Clouds(SMC, LMC), which have metallicities of Z ≈ . , . Z (cid:12) , respec-tively (Korn et al. 2000; Hunter et al. 2007; Trundle et al. 2007).Figure 1 shows the HRD positions of analysed massive stars in theSMC and LMC, adapted from Ramachandran et al. (2018b, 2019).These populations are thought to be complete for the WRs as wellas the luminous (log ( L / L (cid:12) ) (cid:38) .
7) YSGs and RSGs, but far fromcomplete for the OB-type stars. Figure 1 illustrates strongly the ab-sence of RSGs with log ( L / L (cid:12) ) (cid:38) . ( L max / L (cid:12) ) (cid:46) . L max in the SMC comparedto the LMC. Davies et al. (2018) illustrated why this is very likely a genuine physical fact rather than a statistical one, and how thisstands in tension with recent rotating and non-rotating Geneva evo-lution models (Ekstr¨om et al. 2012; Georgy et al. 2013).The challenge of explaining the HD limit solely by radiative in-stabilities might imply that other processes are involved. It was al-ready recognized in the past that rotationally induced mixing canstrongly hinder the redward evolution of massive stars towards theRSG phase (e.g., Maeder & Meynet 2000). However, initial rotationrates in the excess of ≈
200 km s − (at which strong rotational mix-ing starts) seem to be inconsistent with observed rotation rates ofmassive stars (Ram´ırez-Agudelo et al. 2015, 2017; Ramachandranet al. 2019), and models using realistic rotation rates do not repro-duce the apparently Z -independent HD limit (Davies et al. 2018).Higgins & Vink (2020) recently explored this problem by inves-tigating the impact of mixing parameters on RSG stellar models.They found that enhanced semi-convection can reproduce the ob-served HD limit. However, a solution to the HD limit problem needsto consider YSGs and cool BSGs as well, as we argue in our study.In this paper we tackle the puzzle of the HD limit from a di ff erentangle, quantifying the duration spent beyond the HD limit and theexpected number of stars for various stellar evolution models. Forthis purpose we study models with enhanced mixing and constructsynthetic populations based on these models. Moreover, unlike pre-vious studies that considered RSGs alone, we consider simultane-ously RSGs, YSGs, and cool BSGs ( T e ff (cid:46) . ff ect of enhanced MNRAS , 1–13 (2021) he cool supergiant problem and the HD limit Table 1.
Initial compositions in terms of mass fractions for stellar evolutioncalculations. SMC LMC MWH 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Z . . . mixing on stellar evolution tracks. In Section 4 we make compar-isons between observations and synthetic populations that we gen-erate from our stellar evolution tracks. In Section 5 we examine theimpact of modelling assumptions. We discuss our results in relationto previous studies, and summarise, in Section 6. We use the Modules for Experiments in Stellar Astrophysics code(
MESA , version 10398, Paxton et al. 2011, 2013, 2015, 2018) toevolve stellar models with 39 di ff erent zero-age main sequence(ZAMS) masses between M ZAMS = (cid:12) and M ZAMS =
107 M (cid:12) .Three di ff erent initial compositions are used, appropriate for SMC,LMC and Milky Way (MW) stars, as listed in Table 1 (followingHainich et al. 2019). The equation of state employed by
MESA is a blend of the follow-ing equations of state: OPAL (Rogers & Nayfonov 2002), SCVH(Saumon, Chabrier & van Horn 1995), HELM (Timmes & Swesty2000), and PC (Potekhin & Chabrier 2010). Radiative opacitiesare primarily from OPAL (Iglesias & Rogers 1993, 1996), withlow-temperature data from Ferguson et al. (2005) and the high-temperature, Compton-scattering dominated regime according toBuchler & Yueh (1976). Electron conduction opacities follow Cas-sisi et al. (2007).We use the built-in
MESA nuclear reaction network approx21 . TheJoint Institute for Nuclear Astrophysics REACLIB reaction rates(Cyburt et al. 2010) are used, with additional tabulated weak reac-tion rates (Fuller, Fowler & Newman 1985; Oda et al. 1994; Lan-ganke & Mart´ınez-Pinedo 2000) and screening via the prescriptionsof Salpeter (1954), Dewitt, Graboske & Cooper (1973), Alastuey &Jancovici (1978) and Itoh et al. (1979). The formulae of Itoh et al.(1996) are used for thermal neutrino loss rates.
For cool phases ( T e ff ≤
10 000 K), the mass-loss prescription of deJager, Nieuwenhuijzen & van der Hucht (1988) is employed, regard-less of the surface hydrogen mass fraction X s . Wind mass loss is ac-cording to Vink, de Koter & Lamers (2001) for hot ( T e ff ≥
11 000 K) hydrogen-rich phases ( X s ≥ .
4) of the evolution, most notablythe main sequence. If the surface hydrogen mass fraction is lowbut non-negligible (0 . ≤ X s < .
4) then the empirical mass-lossrate relation of Nugis & Lamers (2000) is used. For hot hydrogen-deficient models ( X s < .
1) we follow Yoon (2017) and Woosley(2019), whose prescriptions are based on the mass-loss rates ofHainich et al. (2014) and Tramper, Sana & de Koter (2016). For10 000 K < T e ff <
11 000 K we interpolate between the two regimesdescribed above. The mass-loss rate is multiplied by a factor η w . Wemostly use η w =
1, but in one set of models we explore the impactof boosted mass loss on our results by taking η w = M ∝ Z . . The fit by de Jager et al. (1988),which prevails during most of the CSG phases, has no implicit de-pendence on metallicity. The models have initial rotation velocities at the equator of V i = − . The shellular approximation where the angu-lar velocity is constant over isobars (Meynet & Maeder 1997) isused for rotating models in MESA (Paxton et al. 2013). Mixing pro-cesses induced by rotation are implemented in a di ff usion approx-imation (Paxton et al. 2013) and principally follow Heger, Langer& Woosley (2000). Transport of angular momentum and chemicalmixing caused by internal magnetic fields follows the new prescrip-tion of Fuller, Piro & Jermyn (2019).The e ffi ciency of rotational mixing is set by two parameters, theratio of the turbulent viscosity to the di ff usion coe ffi cient f c (inputparameter am D mix factor in MESA ), and the ratio between the ac-tual molecular weight gradient and the value used for computingthe mixing coe ffi cients f µ (input parameter am gradmu factor in MESA ). These parameters are calibrated to give the observed nitro-gen enhancement for evolved stars. We use the calibration of Hegeret al. (2000), f c = /
30 and f µ = .
05 for most models. However,as discussed by Chie ffi & Limongi (2013), this calibration is notunique (also Potter, Tout & Eldridge 2012). In one set of models weuse alternative values of f c = . f µ = ffi & Limongi2013).Convective regions are defined by the Ledoux criterion (except forone set of models which employs the Schwarzschild criterion) andtreated according to Henyey, Vardya & Bodenheimer (1965) witha mixing-length parameter of α MLT = .
5. Semiconvective mixingin regions which are Ledoux stable but Schwarzschild unstable fol-lows Langer, Fricke & Sugimoto (1983) with an e ffi ciency param-eter α sc = α sc = ffi ciency pa-rameter of α th =
1. We use the so-called MLT ++ implementation fore ffi cient energy transport in convective regions (Paxton et al. 2013).Mixing above the convective core boundary is extended in two ap-proaches. First, for most models, a step overshoot approach is taken(e.g. Shaviv & Salpeter 1973; Maeder & Meynet 1987), where theconvective region is extended by a fraction α ov of the pressure scaleheight H P . The lowest value for α ov that we use is α ov = . α ov = . α ov = . MNRAS000
1. We use the so-called MLT ++ implementation fore ffi cient energy transport in convective regions (Paxton et al. 2013).Mixing above the convective core boundary is extended in two ap-proaches. First, for most models, a step overshoot approach is taken(e.g. Shaviv & Salpeter 1973; Maeder & Meynet 1987), where theconvective region is extended by a fraction α ov of the pressure scaleheight H P . The lowest value for α ov that we use is α ov = . α ov = . α ov = . MNRAS000 , 1–13 (2021)
Gilkis, Shenar, Ramachandran, et al.
Table 2.
Modelling assumptions for all sets of models computed. α sc f c f µ η w criterion1 Ledoux α ov = . /
30 0 .
05 12 Ledoux α ov = .
335 1 1 /
30 0 .
05 13 Ledoux α ov = . /
30 0 .
05 14 Ledoux α ov = . /
30 0 .
05 15 Ledoux α ov = . /
30 0 .
05 16 Ledoux α ov = . /
30 0 .
05 17 Ledoux f ov , JTC /
30 0 .
05 18 Ledoux α ov = . /
30 0 .
05 19 Ledoux α ov = .
335 100 1 /
30 0 .
05 110 Ledoux α ov = . /
30 0 .
05 111 Ledoux α ov = . /
30 0 .
05 112 Ledoux α ov = . /
30 0 .
05 113 Ledoux α ov = . /
30 0 .
05 114 Ledoux f ov , JTC
100 1 /
30 0 .
05 115 Schwarzschild α ov = . − /
30 0 .
05 116 Ledoux α ov = .
335 100 0 . . α ov = .
335 100 1 /
30 0 .
05 2 ditional higher values of α ov = . α ov = α ov = . . Ina second approach we use an exponential core overshooting (e.g.Herwig 2000), where the mixing e ffi ciency decays smoothly outsidethe core, rather than dropping abruptly as in the step overshoot ap-proach. We follow the prescription of Jermyn, Tout & Chitre (2018)to compute the fraction f ov used for the decay scale, f ov = − v c / c s ) + (cid:16) ( H P / r c ) − (cid:17) , (1)where v c / c s is the mass-averaged convective core Mach number, H P is taken at the top of the core, and r c is the convective core radius.The value of f ov is updated according to equation (1) at the end ofevery evolution step. This prescription corresponds to a meridionalcirculation driven by anisotropy in the heat flux emerging from theconvective core. The anisotropy is rotationally-induced, but becausethe convective turnover time is so long in the core this e ff ect saturatesat slower angular velocities than any we consider here. A total of 5967 models were evolved, for 39 initial masses ( M ZAMS =
4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 19, 20, 22, 23, 25,27, 29, 31, 33, 36, 39, 42, 45, 48, 52, 56, 60, 64, 69, 74, 80, 86,92, 99, and 107 M (cid:12) ), 3 initial rotation velocities ( V i = − ), 3 compositions (SMC, LMC, and MW), and 17sets of modelling assumptions as detailed in Table 2. The evolu-tion reached the end of core helium burning for the lower masses( M ZAMS <
10 M (cid:12) ) or core carbon burning for the higher masses( M ZAMS ≥
10 M (cid:12) ). We present several sub-sets of the results to high-light our main findings.Figure 2 shows evolutionary tracks for LMC and MW compo-sitions and an initial rotation velocity of V i =
200 km s − , for α sc =
100 and step overshooting with α ov = .
335 and α ov = . These parameter values might not be physical, but rather used as proxy toinvestigate the issue at hand. dots, to illustrate the duration of di ff erent phases, and the relative ex-pected number of stars at each part of the evolution for each initialmass. It can be seen that for α ov = .
335 numerous models spend along time beyond the HD limit, while for α ov = . α ov = . α ov = . ff ect of overshooting. To test the modelsa statistical analysis is discussed in Section 4.Figure 3 shows evolutionary tracks for LMC and MW composi-tions and an initial rotation velocity of V i =
100 km s − (two timesslower compared to the models in Figure 2), for the same mixingparameters as in Figure 2 ( α sc = α ov = .
335 or α ov = . Γ Edd (computed by
MESA ,taking into account the gas opacity). Compared to the tracks in Fig-ure 2, the models with α ov = . ) and final positions on the HRD . For α ov = . α ov = .
2, the TAMS moves redward for the lower masses, butfor the higher masses it can move bluewards. There are less pre-SN red supergiants for α ov = .
2, though for such a large extent ofovershooting lower initial masses might give rise to CCSNe whoseprogenitors are red supergiants.The duration of the evolutionary stage when a star is located onthe HRD above the HD limit is shown in Figure 4 for various param-eters. For each composition, initial mass, and overshooting scheme,there are 3 di ff erent initial rotation velocities, and we take the longesttime from the three tracks. For α ov = . > ∼
100 000 yr beyond the HD limit. The highest masses with MWcomposition do not cross the HD limit thanks to mass removal bywinds (besides the models with α ov = . ff ect is not present in the SMC and LMCmodels, except for those with significantly enhanced overshooting.With increasing α ov , the LMC and MW models spend less time be-yond the HD limit.For the SMC the models with 30 < ∼ M ZAMS / M (cid:12) < ∼
80 always crossthe HD limit and spend a long time beyond it. This is because cross-ing the HD limit is a result of the combination of mixing and massloss. With increased mixing, envelope material is used as fuel in thecore, both increasing the core mass and luminosity and thereforethe mass loss, while at the same time decreasing the envelope masswhich needs to be removed by winds.The models which use the overshooting prescription of Jermynet al. (2018) give quantitatively similar results to those with α ov = f ov varying slightly with time and initial mass but generallyclose to f ov ≈ . f ov relative to what Claret & The TAMS is defined as the point where the central hydrogen mass fractiondrops below 0 . The final position is marked only for models which ignited carbon in theircentre, and are therefore considered to be CCSN progenitors.MNRAS , 1–13 (2021) he cool supergiant problem and the HD limit Figure 2.
Hertzsprung-Russell diagrams for models with a semiconvective mixing e ffi ciency of α sc =
100 and an overshoot parameter of α ov = .
335 (left)and α ov = . − . Initial masses are in the range4 M (cid:12) ≤ M ZAMS ≤
107 M (cid:12) . Models every 50 000 yr are marked. The thick magenta line maps the HD limit.
Figure 3.
Hertzsprung-Russell diagrams for models with a semiconvective mixing e ffi ciency of α sc =
100 and an overshoot parameter of α ov = .
335 (left)and α ov = . − . Initial masses are in the range4 M (cid:12) ≤ M ZAMS ≤
107 M (cid:12) . The line colour follows the surface Eddington factor. Black crosses show the points in the evolution where core hydrogen burningends, and red asterisks mark the end of the evolution, when core carbon burning ends. The thick magenta line maps the HD limit. MNRAS000
107 M (cid:12) . The line colour follows the surface Eddington factor. Black crosses show the points in the evolution where core hydrogen burningends, and red asterisks mark the end of the evolution, when core carbon burning ends. The thick magenta line maps the HD limit. MNRAS000 , 1–13 (2021)
Gilkis, Shenar, Ramachandran, et al.
Figure 4.
Time spent beyond the HD limit as function of ZAMS mass forSMC (top), LMC (middle) and MW (bottom) composition, for di ff erent coreovershooting prescriptions. For each initial mass, the maximal value out ofall initial rotation rates is taken. Points labeled as f ov , JTC use the exponen-tial overshooting coe ffi cient f ov as function of core properties described byJermyn et al. 2018. Torres (2017) infer from observations. Furthermore, while the usageof f ov is convenient, their mixing mechanism does not exhibit an ex-ponentially decaying geometry like the f ov prescription assumes andso this implementation does not accurately describe the physics. Ourresults imply that there is a strong motivation to further investigatethe mixing mechanism proposed by Jermyn et al. (2018). We compile the luminosities and temperatures of all known CSGsin the Magellanic Clouds with a luminosity of log ( L / L (cid:12) ) = . ff ective Figure 5.
Exponential decay scale as fraction of the pressure scale heightas function of time for several models with MW composition and an initialrotation velocity of V i =
100 km s − as calculated by using equation (1). temperatures of ≈ . T e ff ≤ < T e ff < ≤ T e ff ≤
12 500 K, respectively.For the LMC, we cross-match the RSG list of Davies et al. (2018)with the RSG-YSG list of Neugent et al. (2012). For targets thatappear in both compilations, we adopt temperatures and luminosi-ties from Davies et al. (2018), who derived these parameters froma complete spectral energy distribution (SED) fitting. When temper-atures are not specified, we use calibrations between spectral typesand temperatures by Tabernero et al. (2018) to derive the tempera-ture. Since both studies claim to be complete for log ( L / L (cid:12) ) ≥ . . − . < Bp − Rp < . G < . E B − V = .
09 (e.g. Fitzpatrick & Garmany 1990). We then cross-matched our list with the SIMBAD catalogue to retrieve spectraltypes using the V izier X- match service. Main references are Sand-uleak (1970), Ardeberg et al. (1972), Stock et al. (1976), Evans et al.(2006), and Urbaneja et al. (2017). All identified targets were clas-sified before and have spectral types consistent with BSGs, and themajority of those were included in previous spectroscopic analysesof CSGs in the LMC.All stars with spectral types earlier than B7 in our final list wereremoved, including a few WR stars. Finally We included four LBVsfrom the compilation given by Smith (2019).For all remaining objects, we extracted radial velocities (RVs) andproper motions (PMs) from the SIMBAD database (Wenger et al.2000). The sources of the RVs were predominantly the Gaia DR2catalogue (Gaia Collaboration 2018), Massey & Olsen (2003), Neu-gent et al. (2012), Fehrenbach (1972), and Fehrenbach & Duflot(1982). PMs originate in the Gaia DR2 catalogue for all sources butten, for which they are retrieved from Gaia DR1 (Lindegren et al.2016). The mean PM is 1 .
78 mas yr − with a standard deviation of0 . − , which reflects the measurement limit of Gaia. Thereare 15 outliers with PMs larger than 4 mas yr − within their respec-tive errors, and they are omitted from our sample to ensure that wedo not include foreground Galactic objects.When available, T e ff and log L values for the cool BSGs were MNRAS , 1–13 (2021) he cool supergiant problem and the HD limit adopted from Urbaneja et al. (2017) and Smith (2019). Otherwise,we used spectral-type calibrations by Fitzpatrick & Garmany (1990)to derive the e ff ective temperatures and estimated the extinction pa-rameters based on the expected intrinsic colours. We then used bolo-metric corrections following Flower (1996) and Torres (2010), as-suming a distance of 49 .
97 kpc (Pietrzy´nski et al. 2013). The finallist for the LMC comprises 375 stars: 265 RSGs, 39 YSGs, and 71cool BSGs (four of which are LBVs).For the SMC, we repeat this procedure using the RSGs listed byDavies et al. (2018) and the YSGs listed by Neugent et al. (2010).The RVs, PMs, and spectral types are again extracted using SIM-BAD. The RVs and PMs originate predominantly from the GaiaDR2 catalogue, but also from Massey & Olsen (2003), Neugent et al.(2010), and Gonz´alez-Fern´andez et al. (2015). The spectral types areretrieved from Feast et al. (1960), Dubois et al. (1977), Humphreys(1983), Lennon (1997), and Dufton et al. (2000). Again, all objectsearlier than B7 are removed. We identify about ten outliers withPM > − , which are removed from our sample. When avail-able, T e ff and log L values for the cool BSGs were adopted fromDufton et al. (2000). Otherwise, we use calibrations by Evans &Howarth (2003) to derive the e ff ective temperatures, make the sameassumptions as Neugent et al. (2010) regarding the reddening andthe distance towards the SMC, and use the same relations as aboveto derive the bolometric corrections and luminosities. The final listcomprises 179 stars: 140 RSGs, 7 YSGs, and 32 cool BSGs.We note that accurate derivations of T e ff and log L should rely onthe fitting of SEDs. However, given the statistical nature of our study,the calibrations used above should be su ffi cient for our purpose (e.g.Neugent et al. 2010, 2012). We construct synthetic populations by generating random initialmasses according to a Salpeter IMF. The initial rotation velocity ischosen according to the observed distributions for the LMC (Ra-machandran et al. 2018b) and the SMC (Ramachandran et al. 2019).For both the initial mass and velocity, the nearest values availablein our models are used to chose a stellar evolution track to follow,rather than interpolating between tracks. This results in 58% of SMCmodels being assigned an initial rotation velocity of 100 km s − , 31%getting 200 km s − and 11% with 300 km s − . For the LMC the cor-responding percentages are 80%, 19% and 1%.The stellar age is chosen according to a uniform distribution,corresponding to a constant star-formation rate (SFR) . The stellarproperties are interpolated from the evolutionary track according tothe generated stellar age. If the generated age is longer than the life-time of the computed stellar evolution track, the star is discarded.Random stars are generated until the combined number of RSGsand YSGs ( T e ff < . ≤ log ( L / L (cid:12) ) < . For each of the 17 sets of modelling assumptions (Table 2) we gen-erate 25 random realisations of such populations to get an error esti- This is a reasonable assumption as we are interested in young stars, thoughin general the SFR is not constant in the SMC and the LMC (Section 5.4).
Figure 6.
Example synthetic populations of stars with SMC (top) and LMC(bottom) initial compositions. mate for the computed numbers. We count the number of stars gen-erated which are over-luminous, i.e. those with log ( L / L (cid:12) ) ≥ . α ov , although for e ffi cient semiconvective mixing ( α sc = α ov = α ov = .
335 and α sc =
1. The excess in the simulations increaseswhen including also the YSGs, especially as there is a non-negligiblenumber of such stars in the simulated populations, while the num-ber of observed luminous YSGs is small. When including also coolBSGs, i.e. all stars with T e ff ≤
12 500 K, the excess of over-luminousstars in the simulated populations becomes even larger. An excess inthe lower luminosities for BSGs is acceptable as we do not expectthe sample to be complete for these temperatures, but higher lumi-nosity stars (log ( L / L (cid:12) ) ≥ .
4) would definitely be observed, andtherefore the simulations cannot be taken to properly account for thestellar population.In Figure 9 we show the ratio between the number of RSGs with4 . ≤ log ( L / L (cid:12) ) < . MNRAS000
4) would definitely be observed, andtherefore the simulations cannot be taken to properly account for thestellar population.In Figure 9 we show the ratio between the number of RSGs with4 . ≤ log ( L / L (cid:12) ) < . MNRAS000 , 1–13 (2021)
Gilkis, Shenar, Ramachandran, et al.
Figure 7.
Number of supergiant stars with log ( L / L (cid:12) ) ≥ . ffi ciency, as well as models employing the Schwarzschild stability criterion instead of Ledoux, for the SMC (left) and the LMC (right). For α sc = T e ff ≤ T e ff < T e ff ≤
12 500 K)supergiant stars. and employing MLT ++ results in very luminous stars being hotter(Klencki et al. 2020). In the present study we focus on the most lumi-nous stars, whose numbers are small compared to those with lowerluminosities which mostly a ff ect the YSG / RSG (or BSG / RSG) ratio.Supergiants of all colours need to be taken into account, to separatethe issue of the number ratio of di ff erent colours from the issue ofthe HD limit. Alongside mixing, continuous or eruptive mass loss is a key param-eter that is responsible for stripping a star of its H-rich envelope.The larger the mass loss, the less likely the star is to cross the HDlimit. Motivated by the apparent Z -independence of the HD limit, wechose to focus on mixing here. However, to shed more light on theinterplay between mass loss and mixing in this context, we also pro-vide a set of models with mass-loss rates that are boosted by a factorof two. Recent theoretical and empirical determinations of mass-loss rates during the OB and the RSG phases (e.g., Ramachandran et al.2019; Bj¨orklund et al. 2020; Beasor et al. 2020) rather suggest thatstandard prescriptions such as those used here already lead to anoverestimation of the mass loss. Nevertheless, considering the poorunderstanding of eruptive mass loss, we explore the impact of anincreased mass loss. We do this for the case of α ov = .
335 and α sc = Z -independence of theHD limit. To conclude, while continuous mass loss will play a role inshaping the HD limit, there is little support that it alone can explainthe observations. It is well possible that a Z -independent rapid phaseof mass loss (e.g., in a CSG or LBV phase) needs to be invoked toavoid the discrepancy between observations and theory, but a con-sistent physical framework for implementing it is still missing. MNRAS , 1–13 (2021) he cool supergiant problem and the HD limit Figure 8.
Number of supergiant stars as function of luminosity, for a sample of five di ff erent sets of mixing parameters compared to the observed distributionsfor the SMC (left) and the LMC (right). The top panels show the number of RSGs ( T e ff ≤ T e ff < T e ff ≤
12 500 K) supergiant stars.
The main uncertainty in stellar evolution modelling on which wefocus in this work is mixing in stellar interiors. Therefore, in ad-dition to the e ffi ciency of semiconvection and the extent of over-shooting, the e ff ects of a couple of other mixing assumptions weretested. While we mostly used the Ledoux criterion for defining theconvective boundary, in one set of models (with α ov = . ff ects not only the boundary of the convective core, but alsointermediate regions between the core boundary and the stellar pho-tosphere and the position of models in the HRD (Georgy, Saio &Meynet 2014). The analysis of our synthetic populations generatedfrom stellar evolution tracks which employed the Schwarzschild cri-terion yields a somewhat reduced excess of over-luminous CSGscompared to the analysis with the Ledoux criterion and the sameovershooting extent. Similarly to our results with boosted mass loss(Section 5.1), the excess remains non-negligible, especially whenconsidering the entire relevant temperature range.Up until now we have discussed mixing only for convective re-gions and near their boundaries. In models of rotating stars, consid-erable mixing occurs also in radiative regions, owing to various in-stabilities (e.g., Heger et al. 2000). In one set of models we changed the rotational mixing parameters in radiative regions to f c = . f µ =
1, as described in Section 2.3. Our analysis of the synthetic pop-ulations generated from tracks with these alternative mixing param-eters yields almost no change in the CSG excess (Fig. 7). This mightbe a result of the possible degeneracy of these parameters (Chie ffi &Limongi 2013).As an additional test of the sensitivity to the implementation ofrotational mixing, we also compare our results to the Geneva trackswith MW and SMC compositions (Ekstr¨om et al. 2012; Georgy et al.2013). In Figure 10 we show the time spent beyond the HD limit forseveral of our modelling assumptions and for the Geneva models,for which rotational mixing is treated as an advective-di ff usive pro-cess, compared to the di ff usion approximation adopted by MESA . TheGeneva models spend time periods beyond the HD limit rather sim-ilar to our tracks computed with
MESA . While a quantitative compar-ison using synthetic populations will give a more definitive answer,we surmise that the computed excess of CSGs will be similar withthe Geneva tracks.
MNRAS000
MNRAS000 , 1–13 (2021) Gilkis, Shenar, Ramachandran, et al.
Table 3.
Excess of over-luminous stars in synthetic populations.overshooting α sc N excess , RSG + YSG N excess , RSG
SMC α ov = . ±
15 79 ± α ov = .
335 1 20 ± ± α ov = . ± ± α ov = . ± ± α ov = . ± ± α ov = . ± ± f ov , JTC ± ± α ov = . ± ± α ov = .
335 100 10 ± ± α ov = . ± ± α ov = . ± ± α ov = . ± ± α ov = . ± ± f ov , JTC
100 7 ± ± α ov = . ±
15 103 ± α ov = .
335 1 41 ± ± α ov = . ± ± α ov = . ± ± α ov = . ± ± α ov = . ± ± f ov , JTC ± ± α ov = . ± ± α ov = .
335 100 44 ± ± α ov = . ± ± α ov = . ± ± α ov = . ± ± α ov = . ± ± f ov , JTC
100 14 ± ± Table 4.
Excess of over-luminous stars in additional synthetic populations.model assumptions N excess , RSG + YSG N excess , RSG
SMC Schwarzschild 4 ± ± f c = . f µ = ± ± η w = ± ± ± ± f c = . f µ = ± ± η w = ± ± In principle, the presence of a close companion (orbital period P (cid:46) ff ected by binary interactions, the e ff ect of stellar multiplic-ity is not trivial. If the binary separation shrinks as the binary massincreases, then one may expect binary interactions to help resolve Figure 9.
Number of YSG stars relative to the number of RSG stars in syn-thetic populations as function of the overshooting extent, for two values ofthe semiconvective mixing e ffi ciency, for the SMC (top) and the LMC (bot-tom). the discrepancy. However, solid evidence for this at the upper-massend is currently lacking. Addressing this question in future studiesshould be helpful to quantitatively constrain the impact of multiplic-ity on the HD limit. There are strong indications that the recent star-formation history inthe Magellanic Clouds, especially in the SMC, was not continuousbut is rather characterised by various peaks over the past 100 Myr(Indu & Subramaniam 2011; Ramachandran et al. 2019; Fulmeret al. 2020). As such, this can have an important impact on thesynthetic populations obtained, and hence on the final conclusions.However, it is unlikely that star-formation history alone can explainthe discrepancy in the SMC and the LMC simultaneously. Moreover,the presence of WR stars in both galaxies that span a substantial lu-minosity range (Hainich et al. 2014, 2015; Shenar et al. 2016, 2019)suggests that a similar distribution of CSGs across the luminosityrange could be anticipated. Given the many uncertainties, repeatingthis investigation with detailed star-formation histories is beyond thescope of our work, but should be explored in future studies.
We evolve numerous grids of stellar evolution tracks, for MW, LMCand SMC compositions, with a variety of mixing prescriptions. Wefind that enhanced mixing diverts stellar evolution tracks from the“forbidden region” defined by the HD limit (Figs. 2-3). Based uponthese grids of stellar models we construct synthetic populationswith initial compositions appropriate for the LMC and SMC, andmake quantitative comparisons with the observed populations ofcool ( T e ff ≤ . MNRAS , 1–13 (2021) he cool supergiant problem and the HD limit Figure 10.
Time spent beyond the HD limit as function of ZAMS mass forSMC (top) and MW (bottom) composition, for stellar evolution tracks asdescribed in the inset.
We find that enhanced mixing reduces the excess of over-luminousCSGs in our simulated stellar populations. While for RSGs theredoes not seem to be a severe problem, the tension between observa-tions and simulations increases with increasing the upper tempera-ture cuto ff . We can therefore consider the existence of a “Cool Su-pergiant Problem” as the apparent mismatch between observed su-pergiants and the results of stellar evolution models.When considering only RSGs, we find that the dependence ofthe excess on the overshooting extent is non-monotonic for e ffi cientsemiconvective mixing ( α sc = ffi -cient mixing in regions of semiconvection to account for the proper-ties of supergiants in the LMC and SMC. Schootemeijer et al. (2019)also claim that convective overshooting can be constrained by theproperties of the populations, such as the ratio between RSGs andBSGs. We di ff er from Schootemeijer et al. (2019) in two regards:( i ) We consider much higher values of α ov , beyond an apparent ex-tremum around α ov ≈ .
5; ( ii ) Schootemeijer et al. (2019) focus onthe bulk of the supergiant population, while we investigate the starswith the highest luminosities, which are a small fraction of the over-all population.Higgins & Vink (2020) suggest that the HD limit can be explainedwith decreased overshoot mixing, α ov = .
1, though they consideronly the red part of the evolution. We show that for red supergiantstaking α ov = . α sc =
100 indeed gives a reasonableaccount of the supergiant populations in the LMC and SMC (Figs.7-8). This is explained by the significant e ff ect of e ffi cient semicon-vection in increasing the e ff ective surface temperatures of the stellarmodels. Thus, for e ffi cient semiconvection, core He-burning starstend to appear as YSGs or BSGs instead of RSGs. However, thisdoes not solve the supergiant excess beyond the HD limit, but merely“sweeps it under the carpet”. When we consider the entire tempera-ture range to which the horizontal HD limit applies ( T e ff (cid:46) . ff ects the formation of WR stars. Enhanced mixing reduces themass threshold for a star to remove its envelope through winds, andtherefore a ff ects the necessity of binary interactions in forming WRstars (Shenar et al. 2020), as well as decreasing the single-star evo-lution WR luminosity threshold. For example, Figs. 2-3 show thatincreasing α ov from 0 .
335 to 1 . ( L / L (cid:12) ) ≈ ≈ . ≈ . ≈ .
25 for the MW.Our results also have implications for energetic transient astro-physical phenomena. Firstly, the initial mass threshold for core-collapse supernovae is lower for enhanced mixing. Secondly, the po-sition in the HRD where stellar models end their lives is substantiallyhotter for enhanced mixing. This has implications for the so-called“Red Supergiant Problem”, which is the apparent discrepancy be-tween the most luminous progenitor of a Type IIP supernova andthe most luminous known RSGs (Smartt et al. 2009; Horiuchi et al.2014; Davies & Beasor 2018, 2020; Kochanek 2020). Our modelswith lower overshooting values exhibit such a problem, though forenhanced mixing the more massive stars do not die as RSGs, in ac-cordance with observations.The formation of black hole binaries (which are progenitors ofgravitational wave events) depends on the expansion of massive starsand their interaction with a companion. Klencki et al. (2020) discussthe role of metallicity in stellar expansion and interaction, as withlower metallicity stars are generally more compact because of thelower gas opacity. Smaller stellar radii imply that a smaller frac-tion of stars will experience significant binary interactions, such ascommon envelope evolution, that lead to short-period binary blackholes which will merge quickly enough to produce observable grav-itational wave events. Our models with enhanced mixing expand tosmaller radii in general, therefore a ff ecting the evolution towardsmerging black holes. Since a large fraction of massive stars interactwith a companion during their evolutionKlencki et al. (2020) also discuss the role of MLT ++ in relationto the HD limit. According to Klencki et al. (2020), the more limitedexpansion of massive stellar models employing MLT ++ is favorablein terms of explaining the HD limit, and it might be more accuratefor rather massive ( M > ∼
50 M (cid:12) ) stars. We note that even though weuse the favorable MLT ++ prescription we still find an excess of over-luminous stars in the temperature range relevant for the HD limit.It is important to stress that our study focuses on the upper-massend. For example, applying high α ov values as invoked here cansuppress blue loops (Stothers & Chin 1991; Walmswell, Tout &Eldridge 2015) on the HRD (Fig. 2) and prevent lower-mass starsfrom becoming Cepheids (Anderson et al. 2014, 2016), in contrast toobservations. However, studies of intermediate-mass stars suggest amass-dependent convective overshoot extent (Claret & Torres 2016,2017, 2018, 2019). These studies suggest an increase up to α ov ≈ . M ≈ (cid:12) , plateauing afterwards. A similar trend is explained byJermyn et al. (2018). Overshooting is less well-constrained for thehigher masses which are relevant for our study and for the existenceof the HD limit. Models of massive stars make various calibrationsfor the overshooting extent, with results ranging from α ov = . . (cid:12) (Ekstr¨om et al. 2012) to α ov = . M > ∼
30 M (cid:12) (Higgins & Vink2019). So while very high overshoot values cannot be applied to thelowest masses in our grids ( M = (cid:12) ), higher mass stars might MNRAS000
30 M (cid:12) (Higgins & Vink2019). So while very high overshoot values cannot be applied to thelowest masses in our grids ( M = (cid:12) ), higher mass stars might MNRAS000 , 1–13 (2021) Gilkis, Shenar, Ramachandran, et al. exhibit behaviour appropriate for enhanced mixing, with few realconstraints on the extent.In conclusion, we propose that internal mixing in massive starsmight play an important part in explaining the empiric HD limit.Enhanced mixing prevents the redward evolution of stellar mod-els towards or beyond the HD limit. We do not suggest to adopt ahigher overshooting parameter, but rather that our results hint at adeficiency in the modelling of mixing in stellar interiors (e.g., Aertset al. 2019; Schootemeijer et al. 2019). The extent of core overshoot-ing might be highly mass-dependent, or rotational mixing is moree ffi cient for slowly rotating stars. Moreover, it is well possible thatthe final explanation relies on a multitude of mechanisms (e.g., mix-ing, mass loss, multiplicity, star formation history; see Sect. 5). Wetherefore encourage future investigations of this problem that ad-dress these various mechanisms. For now, internal mixing in massivestars remains an unresolved issue in stellar modelling with broad im-plications, and the origin of the Humphreys-Davidson limit remainsuncertain. ACKNOWLEDGMENTS
We acknowledge the constructive report by our referee Dr. CyrilGeorgy. This research has made use of the SIMBAD database andcross-match service operated at CDS, Strasbourg, France. The re-search has received funding from the European Research Council(ERC) under the European Union’s Horizon 2020 research and in-novation programme with grant number 772225: MULTIPLES. AGand TS would like to thank Danny Lennon and Ben Davies for help-ful input and exchanges. The Flatiron Institute is supported by theSimons Foundation. IA is a CIFAR Azrieli Global Scholar in theGravity and the Extreme Universe Program and acknowledges sup-port from that program, from the European Research Council (ERC)under the European Union’s Horizon 2020 research and innovationprogram (grant agreement number 852097), from the Israel ScienceFoundation (grant numbers 2108 /
18 and 2752 / DATA AVAILABILITY STATEMENT
The data underlying this article will be shared on reasonable requestto the corresponding author.
REFERENCES
Aerts C., Mathis S., Rogers T. M., 2019, ARA&A, 57, 35Alastuey A., Jancovici B., 1978, ApJ, 226, 1034Almeida L. A., et al., 2017, A&A, 598, A84Anderson R. I., Ekstr¨om S., Georgy C., Meynet G., Mowlavi N., Eyer L.,2014, A&A, 564, A100Anderson R. I., Saio H., Ekstr¨om S., Georgy C., Meynet G., 2016, A&A,591, A8Ardeberg A., Brunet J. P., Maurice E., Prevot L., 1972, A&AS, 6, 249Beasor E. R., Davies B., Smith N., van Loon J. T., Gehrz R. D., Figer D. F.,2020, MNRAS, 492, 5994Bestenlehner J. M., et al., 2011, A&A, 530, L14Bj¨orklund R., Sundqvist J. O., Puls J., Najarro F., 2020, arXiv e-prints, p.arXiv:2008.06066Bouret J. C., Lanz T., Martins F., Marcolino W. L. F., Hillier D. J., DepagneE., Hubeny I., 2013, A&A, 555, A1Brott I., et al., 2011, A&A, 530, A115 Buchler J. R., Yueh W. R., 1976, ApJ, 210, 440Cassisi S., Potekhin A. Y., Pietrinferni A., Catelan M., Salaris M., 2007, ApJ,661, 1094Castro N., Oey M. S., Fossati L., Langer N., 2018, ApJ, 868, 57Chie ffi A., Limongi M., 2013, ApJ, 764, 21Claret A., Torres G., 2016, A&A, 592, A15Claret A., Torres G., 2017, ApJ, 849, 18Claret A., Torres G., 2018, ApJ, 859, 100Claret A., Torres G., 2019, ApJ, 876, 134Crowther P. A., Dessart L., Hillier D. J., Abbott J. B., Fullerton A. W., 2002,A&A, 392, 653Crowther P. A., Schnurr O., Hirschi R., Yusof N., Parker R. J., Goodwin S. P.,Kassim H. A., 2010, MNRAS, 408, 731Cyburt R. H., et al., 2010, ApJS, 189, 240Davies B., Beasor E. R., 2018, MNRAS, 474, 2116Davies B., Beasor E. R., 2020, MNRAS, 493, 468Davies B., Crowther P. A., Beasor E. R., 2018, MNRAS, 478, 3138Dewitt H. E., Graboske H. C., Cooper M. S., 1973, ApJ, 181, 439Drout M. R., Massey P., Meynet G., 2012, ApJ, 750, 97Dubois P., Jaschek M., Jaschek C., 1977, A&A, 60, 205Dufton P. L., McErlean N. D., Lennon D. J., Ryans R. S. I., 2000, A&A, 353,311Dufton P. L., Evans C. J., Hunter I., Lennon D. J., Schneider F. R. N., 2019,A&A, 626, A50Ekstr¨om S., et al., 2012, A&A, 537, A146Evans C. J., Howarth I. D., 2003, MNRAS, 345, 1223Evans C. J., Lennon D. J., Smartt S. J., Trundle C., 2006, A&A, 456, 623Feast M. W., Thackeray A. D., Wesselink A. J., 1960, MNRAS, 121, 337Fehrenbach C., 1972, A&A, 19, 427Fehrenbach C., Duflot M., 1982, A&AS, 48, 409Ferguson J. W., Alexander D. R., Allard F., Barman T., Bodnarik J. G.,Hauschildt P. H., He ff ner-Wong A., Tamanai A., 2005, ApJ, 623, 585Fitzpatrick E. L., Garmany C. D., 1990, ApJ, 363, 119Flower P. J., 1996, ApJ, 469, 355Fuller G. M., Fowler W. A., Newman M. J., 1985, ApJ, 293, 1Fuller J., Piro A. L., Jermyn A. S., 2019, MNRAS, 485, 3661Fulmer L. M., Gallagher J. S., Hamann W.-R., Oskinova L. M., Ramachan-dran V., 2020, A&A, 633, A164Gaia Collaboration 2018, VizieR Online Data Catalog, p. I / , 1–13 (2021) he cool supergiant problem and the HD limit Kalari V. M., Vink J. S., Dufton P. L., Fraser M., 2018, A&A, 618, A17Kippenhahn R., Ruschenplatt G., Thomas H. C., 1980, A&A, 91, 175Klencki J., Nelemans G., Istrate A. G., Pols O., 2020, A&A, 638, A55Kochanek C. S., 2020, MNRAS, 493, 4945K¨ohler K., et al., 2015, A&A, 573, A71Korn A. J., Becker S. R., Gummersbach C. A., Wolf B., 2000, A&A, 353,655Lamers H. J. G. L. M., Fitzpatrick E. L., 1988, ApJ, 324, 279Langanke K., Mart´ınez-Pinedo G., 2000, Nuclear Physics A, 673, 481Langer N., Fricke K. J., Sugimoto D., 1983, A&A, 126, 207Lennon D. J., 1997, A&A, 317, 871Levesque E. M., Massey P., Olsen K. A. G., Plez B., Josselin E., Maeder A.,Meynet G., 2005, ApJ, 628, 973Lindegren L., et al., 2016, A&A, 595, A4Maeder A., Meynet G., 1987, A&A, 182, 243Maeder A., Meynet G., 2000, A&A, 361, 159Mahy L., et al., 2020, A&A, 634, A118Massey P., Evans K. A., 2016, ApJ, 826, 224Massey P., Olsen K. A. G., 2003, AJ, 126, 2867Meynet G., Maeder A., 1997, A&A, 321, 465Moe M., Di Stefano R., 2017, ApJS, 230, 15Neugent K. F., Massey P., Ski ff B., Drout M. R., Meynet G., Olsen K. A. G.,2010, ApJ, 719, 1784Neugent K. F., Massey P., Ski ff B., Meynet G., 2012, ApJ, 749, 177Nugis T., Lamers H. J. G. L. M., 2000, A&A, 360, 227Oda T., Hino M., Muto K., Takahara M., Sato K., 1994, Atomic Data andNuclear Data Tables, 56, 231Ohnaka K., Driebe T., Hofmann K. H., Weigelt G., Wittkowski M., 2008,A&A, 484, 371Paxton B., Bildsten L., Dotter A., Herwig F., Lesa ff re P., Timmes F., 2011,ApJS, 192, 3Paxton B., et al., 2013, ApJS, 208, 4Paxton B., et al., 2015, ApJS, 220, 15Paxton B., et al., 2018, ApJS, 234, 34Pietrzy´nski G., et al., 2013, Nature, 495, 76Potekhin A. Y., Chabrier G., 2010, Contributions to Plasma Physics, 50, 82Potter A. T., Tout C. A., Eldridge J. J., 2012, MNRAS, 419, 748Ramachandran V., Hainich R., Hamann W. R., Oskinova L. M., Shenar T.,Sander A. A. C., Todt H., Gallagher J. S., 2018a, A&A, 609, A7Ramachandran V., Hamann W. R., Hainich R., Oskinova L. M., Shenar T.,Sander A. A. C., Todt H., Gallagher J. S., 2018b, A&A, 615, A40Ramachandran V., et al., 2019, A&A, 625, A104Ram´ırez-Agudelo O. H., et al., 2015, A&A, 580, A92Ram´ırez-Agudelo O. H., et al., 2017, A&A, 600, A81Rogers F. J., Nayfonov A., 2002, ApJ, 576, 1064Salpeter E. E., 1954, Australian Journal of Physics, 7, 373Sana H., et al., 2012, Science, 337, 444Sanduleak N., 1970, Contributions from the Cerro Tololo Inter-AmericanObservatory, 89Sanyal D., Langer N., Sz´ecsi D., -C Yoon S., Grassitelli L., 2017, A&A, 597,A71Saumon D., Chabrier G., van Horn H. M., 1995, ApJS, 99, 713Schneider F. R. N., et al., 2018, A&A, 618, A73Schootemeijer A., Langer N., 2018, A&A, 611, A75Schootemeijer A., Langer N., Grin N. J., Wang C., 2019, A&A, 625, A132Shaviv G., Salpeter E. E., 1973, ApJ, 184, 191Shenar T., et al., 2016, A&A, 591, A22Shenar T., et al., 2017, A&A, 598, A85Shenar T., et al., 2019, A&A, 627, A151Shenar T., Gilkis A., Vink J. S., Sana H., Sander A. A. C., 2020, A&A, 634,A79Smartt S. J., Eldridge J. J., Crockett R. M., Maund J. R., 2009, MNRAS, 395,1409Smith N., 2019, MNRAS, 489, 4378Stock J., Osborn W., Ibanez M., 1976, A&AS, 24, 35Stothers R. B., Chin C.-W., 1991, ApJ, 374, 288Tabernero H. M., Dorda R., Negueruela I., Gonz´alez-Fern´andez C., 2018,MNRAS, 476, 3106 Tehrani K. A., Crowther P. A., Bestenlehner J. M., Littlefair S. P., PollockA. M. T., Parker R. J., Schnurr O., 2019, MNRAS, 484, 2692Timmes F. X., Swesty F. D., 2000, ApJS, 126, 501Torres G., 2010, AJ, 140, 1158Tramper F., et al., 2015, A&A, 581, A110Tramper F., Sana H., de Koter A., 2016, ApJ, 833, 133Trundle C., Lennon D. J., 2005, A&A, 434, 677Trundle C., Lennon D. J., Puls J., Dufton P. L., 2004, A&A, 417, 217Trundle C., Dufton P. L., Hunter I., Evans C. J., Lennon D. J., Smartt S. J.,Ryans R. S. I., 2007, A&A, 471, 625Ulmer A., Fitzpatrick E. L., 1998, ApJ, 504, 200Urbaneja M. A., Kudritzki R. P., Gieren W., Pietrzy´nski G., Bresolin F., Przy-billa N., 2017, AJ, 154, 102Vink J. S., de Koter A., Lamers H. J. G. L. M., 2001, A&A, 369, 574Vink J. S., Brott I., Gr¨afener G., Langer N., de Koter A., Lennon D. J., 2010,A&A, 512, L7Walmswell J. J., Tout C. A., Eldridge J. J., 2015, MNRAS, 447, 2951Wenger M., et al., 2000, A&AS, 143, 9Wolf B., 1972, A&A, 20, 275Woosley S. E., 2019, ApJ, 878, 49Yoon S.-C., 2017, MNRAS, 470, 3970de Jager C., Nieuwenhuijzen H., van der Hucht K. A., 1988, A&AS, 72, 259MNRAS000