EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 707 (March 2026) A&A, 707 (2026) A11 Full HTML
Open Access
Issue
A&A
Volume 707, March 2026
Article Number A11
Number of page(s) 18
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/202557766
Published online 26 February 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1. Introduction

Massive stars (Mini ≳ 8 M⊙), although relatively rare, are the powerful engines of the Universe. Through their intense ionizing radiation, strong stellar winds, Luminous Blue Variable (LBV) eruptions, and supernova explosions, they shape their surroundings and are the regulators of star formation and sources of chemical enrichment within their host galaxies (e.g., Maeder 1983; Dray et al. 2003; Hopkins et al. 2014; Ramachandran et al. 2018; Crowther 2019). Despite their significance, the evolution of massive stars, particularly at the high-mass end, where their impact is the strongest, remains poorly understood. The predicted evolution of these stars is highly dependent on the assumptions about internal mixing processes and mass-loss rates applied during different evolutionary stages (e.g., Shenar et al. 2020a; Gilkis et al. 2021; Sabhahit et al. 2022; Zapartas et al. 2025).

One of the major challenges facing current stellar evolution models is their tendency to overproduce extremely luminous red supergiant (RSG) stars (e.g., Brott et al. 2011; Georgy et al. 2013; Choi et al. 2016; Limongi & Chieffi 2018; Costa et al. 2025) beyond the empirical Humphreys-Davidson (HD) limit (log L/L⊙ ≳ 5.5), where no stars are observed (Humphreys & Davidson 1979; Davies et al. 2018). This issue has been the subject of numerous studies exploring the effects of mixing processes and different RSG mass-loss rates on the synthetic RSG populations across various metallicities (Schootemeijer & Langer 2018; Klencki et al. 2020; Higgins & Vink 2020; Zapartas et al. 2025; Decin et al. 2024). However, a solution to this problem remains elusive.

Luminous Blue Variables are massive stars that have evolved into a region in the Hertzsprung-Russell diagram (HRD) where their outer atmosphere is barely gravitationally bound, leading to enhanced mass-loss rates with possibly episodic outbursts. During the LBV phase, a star can lose several solar masses, making it a critical component of massive stellar evolution that could potentially alleviate some of the tension of the overluminous RSG problem. However, our understanding of the LBV phase is limited. Its impact on the evolution of massive stars is unclear (Smith 2026), and due to the absence of a consistent framework for modeling LBV eruptions, it is often ignored in detailed stellar evolution calculations. Early pioneering modeling attempts tried to incorporate the effect of LBV eruptions in stellar evolution calculations, by increasing the mass-loss rates of stars massive enough to undergo an LBV phase as soon as they enter a specific evolutionary stage (such as the post-main sequence, Vanbeveren 1991) or to enhance the mass-loss in regions unstable against strange-mode pulsations (Langer et al. 1994). In recent years, progress has been made in understanding the physics, impact, and modeling of LBV mass ejections (Grassitelli et al. 2021; Mukhija & Kashi 2024; Cheng et al. 2024). Nevertheless, a physical model that prevents stellar evolution simulations from surpassing the HD limit for a significant duration, especially at low metallicity (Z), is still lacking.

A third major issue is the difficulty to explain the existence of the least luminous, apparently single WR stars, particularly in metal-poor galaxies (Shenar et al. 2020a). The strength of UV radiation-driven stellar winds weakens (e.g., Vink et al. 2001; Mokiem et al. 2007; Hainich et al. 2017; Pauli et al. 2025), which leads to the expectation that the minimum luminosity of WR stars formed by single stars should increase with decreasing Z. Consequently, stellar evolution models predict that the faintest WR stars in low-metallicity environments are all post-interaction binary products (e.g., Smith 2014; Massey et al. 2021; Pauli et al. 2022). Contrarily, dedicated multi-epoch surveys of the least-luminous WR stars in low-Z environments strongly suggest that these stars are currently single (Foellmi et al. 2003; Schootemeijer et al. 2024). To explain the presence of these faint single WR stars, single-star stellar evolution calculations require enhanced mass-loss rates (e.g., Yoon 2017; Sabin et al. 2022). Yet this approach is in contrast to empirical mass-loss measurements, which over the past few decades have exposed lower mass-loss rates for hot stars than those commonly assumed in standard mass-loss prescriptions (Hainich et al. 2017; Shenar et al. 2019; Pauli et al. 2025).

The Large (LMC) and Small Magellanic Clouds (SMC) are two nearby galaxies with low metallicities (i.e., Fe-group elements) of ZLMC = 0.5 Z⊙ and ZSMC = 0.14 Z⊙, respectively (Vink et al. 2023). These galaxies have been subject to numerous observing campaigns, and their massive star content is mostly known. In particular, their RSG populations (Davies et al. 2018; Yang et al. 2023) as well as the WR population, including their binary properties (Neugent et al. 2018; Bartzakos et al. 2001a,b; Foellmi et al. 2003; Hainich et al. 2014, 2015; Shenar et al. 2019), are mostly known. This makes these galaxies ideal test beds to benchmark stellar evolution models.

In this work, we present a new physically motivated and empirically calibrated method for modeling Eddington-limit induced mass ejections in 1D stellar evolution codes. We investigate how these mass ejections influence massive star populations. We compare the resulting synthetic populations with the observed massive star populations of the LMC and SMC.

Section 2 outlines the physical assumptions incorporated into our stellar evolution calculations, the new model for the Eddington-limit induced mass ejections, and the assumptions used to generate the synthetic stellar populations. The resulting synthetic populations are presented in Sect. 3. Section 4 discusses their sensitivity to the assumptions made during the modeling process as well as their implications for understanding stellar populations. Our results are summarized in Sect. 5.

2. Methods

2.1. Stellar evolution modeling

We employed the Modules for Experiments in Stellar Astrophysics (MESA) code, version 24.08.1 (Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023), to calculate stellar evolution models at LMC (ZLMC = 6.17 × 10−3) and SMC (ZSMC = 2.44 × 10−3) metallicity. Our models cover initial masses spanning from Mini ≈ 10 M⊙  −  400 M⊙ (log(Mini/M⊙) = 1–2.65) in steps of Δlog(Mini/M⊙) = 0.05. For simplicity, we assume that all stars are born with a moderate rotation velocity of vrot,ini = 100 km s−1 (i.e., where the observed vrot-distribution peaks, e.g., see Ramírez-Agudelo et al. 2015; Ramachandran et al. 2019). In total, 33 single-star models per initial metallicity and per input physics variation (see Sect. 2.3) are calculated.

A detailed list of the input physics is provided in Appendix A, and our source files can be downloaded from Zenodo1. Here, only the most relevant physics are mentioned. Our models incorporate semiconvective mixing with αsc = 1 (Langer et al. 1983; Schootemeijer et al. 2019), and step core overshooting with fov = 0.345 and fov,  0 = 0.01 (Brott et al. 2011) during core hydrogen and helium burning. To account for stabilizing chemical gradients that may limit the growth of the overshooting regions, we include the Brunt-Vaisala frequency and set the stabilizing Brunt composition gradient to B = 0.1. As a result, in most of our calculations, overshooting during core helium burning occurs primarily in the very early stages, before a steep chemical gradient between the convective core and the envelope forms.

During the late evolutionary stages of the most massive stars (Mini ≳ 40 M⊙), with helium core masses exceeding 13 M⊙, we encountered convergence issues. To address this, we employed for these models after core hydrogen depletion the superadiabatic reduction method of MESA (Jermyn et al. 2023, Section 7.2). This method enhances energy transport locally in layers close to the Eddington limit, resulting in less impact on the final stellar structure compared to MLT++. For the superadiabatic reduction, we set the energy transport enhancement factors to α1 = α2 = 3.5, which allows our models to converge within a reasonable time frame while minimizing deviations from standard mixing length theory. Since this option is only turned on after core hydrogen depletion, all of our core hydrogen burning models approaching the Eddington limit locally are affected by envelope inflation (Sanyal et al. 2015).

In our standard setup, we account for stellar winds as follows. For hot stars with radiative envelopes, we apply the recent empirical mass-loss recipe from Pauli et al. (2025), which depends on the classical Eddington factor Γ​e and the metallicity Z. However, this approximation is only valid for optically thin winds, and it breaks down for stars approaching the Eddington limit. Therefore, when logΓ​e > −0.4, we switch to the empirical Wolf-Rayet (WR) mass-loss recipe from Shenar et al. (2020b) if the mass-loss rate is higher than the one from Pauli et al. (2025). To avoid abrupt changes in the mass-loss rates that could lead to numerical issues, we linearly interpolate between the hot star and WR mass-loss recipes in the range from logΓ​e = −0.5  to   − 0.4. For RSGs, we adopt the empirical mass-loss recipe from Yang et al. (2023), which has been shown to best reproduce the RSG population in the SMC (Zapartas et al. 2025). While the launching mechanism of RSG winds remains a topic of debate, recent theoretical models suggest that these winds are powered by turbulence near the stellar surface and are therefore Z-independent (Kee et al. 2021). This is further supported by observational evidence, which shows that even at extremely low Z, no RSG with log(L/L⊙)≳5.5 is observed (Davies et al. 2018; McDonald et al. 2022; Schootemeijer et al. 2026). To switch to an RSG wind in our models, we require that the layers where hydrogen can recombine (T < 10 000 K) are convective and account for more than 0.1% of the total hydrogen envelope mass. Further, we account for Eddington-limit induced mass ejections as associated with LBVs or the radial pulsations of hydrodynamically unstable RSG envelopes when approaching the Eddington limit (e.g., Heger et al. 1997; Yoon & Cantiello 2010; Bronner et al. 2025; Laplace et al. 2025).

2.2. Eddington-limit induced mass ejections

A connection of the LBV phenomenon with the Eddington limit is evident since the pioneering work of Lamers & Fitzpatrick (1988) and Ulmer & Fitzpatrick (1998), who showed that the empirical HD-limit is essentially recovered by the Eddington limit in stellar model atmospheres. This connection is strengthened by 1D stellar structure and evolution models (Ishii et al. 1999; Gräfener et al. 2012; Köhler et al. 2015; Sanyal et al. 2015; Grassitelli et al. 2021), which found that stars approaching their Eddington limit locally in their envelopes inflate their radius to allow for a sufficient radiative energy transport, a phenomenon typical for hydrostatic codes that cannot be recovered by 3D radiation-hydrodynamic calculations (Jiang et al. 2015, 2018; Poniatowski et al. 2021).

In contrast to Cheng et al. (2024), who model eruptive mass loss in 1D massive star models based on the local super-Eddington luminosity, we choose an empirically calibrated estimate of the global instability of the inflated stellar envelopes. In the absence of systematic stability analyses of such envelopes, we assume here that a stronger inflation leads to a higher chance for instability. A possible instability mechanism is provided by strange mode pulsations (Gautschy et al. 1990) which have been shown to occur in supergiant stellar models near the Eddington limit (Kiriakidis et al. 1993; Sanyal & Langer 2013; Lovekin & Guzik 2014), and which have growth time scales that are close to the hydrodynamic time scale of the envelope. To approximate Eddington-limit induced mass ejections, we assume that a star can inflate to a certain degree without becoming unstable, but that there must be a threshold value after which a star will expel its envelope until it returns to a “stable” inflation (i.e., when inflation is again below the given threshold). In our models, we fully account for the effect of envelope inflation during core hydrogen burning. During core helium burning, we apply the local superadiabatic reduction method from MESA for numerical stability, which still allows for significant inflation.

To determine the extent of inflation R/Rcore in a stellar evolution model with radius R, we need to compute the radius that the star would have if inflation did not occur, which we designate Rcore. This radius can technically be recovered by artificially increasing the efficiency of convective energy transport, for instance, by increasing the mixing-length parameter (Sanyal et al. 2015; Nathaniel et al. 2025). However, this is computationally expensive. Instead, we follow Sanyal et al. (2015), who showed that Rcore can be well approximated by the location in the stellar envelope at which the gas pressure (Pgas) contribution to the total pressure (P) drops below Pgas/P ≃ 0.15. As shown by Sanyal et al. (2015), the result is not sensitive on the threshold value. Increasing the threshold to 0.30 impacts Rcore by less than 5%. We confirm that this also holds for our model calculations.

To calibrate the extent of inflation beyond which Eddington-limit induced mass ejections should be triggered, we calculated stellar evolution models up to core hydrogen depletion using the physical setup described in Sect. 2.1. In Fig. 1, we compare the model tracks to the observed stellar populations in the HRD, for the LMC and the SMC. Contours indicate the amount of inflation in our models across the HRD. Interestingly, the slope of lines with fixed inflation values is parallel to the HD limit, supporting our hypothesis that a certain amount of inflation in the stellar envelope might be stable. One can see that core hydrogen-burning stars exceed the HD limit when their inflation reaches R/Rcore ≈ 2 (see also Sanyal et al. 2017), which we adopt as the triggering criterion for Eddington-limit induced mass ejections for such stars.

Thumbnail: Fig. 1. Refer to the following caption and surrounding text. Fig. 1.

Hertzsprung-Russell diagrams of the massive star population (symbols) in the LMC (left) and SMC (right). The different symbols mark the positions of observed RSG (small red plusses, Davies et al. 2018; Yang et al. 2023), YSG (crosses, Neugent et al. 2010, 2012), LBVs (filled triangles with quiescent and outburst positions, Humphreys et al. 2016; Kalari et al. 2018), OB stars (open circles and squares, Evans et al. 2004; Trundle et al. 2004; Trundle & Lennon 2005; Hunter et al. 2008; Bestenlehner et al. 2011; Bouret et al. 2013; Urbaneja et al. 2017; Castro et al. 2018; Schneider et al. 2018; Ramachandran et al. 2018, 2019; Dufton et al. 2019; Bestenlehner 2020; Bouret et al. 2021; Rickard et al. 2022; Pauli et al. 2023; Bernini-Peron et al. 2024; Gómez-González et al. 2025; Alkousa et al. 2025), partially stripped stars (white large plusses Pauli et al. 2022; Rickard & Pauli 2023; Ramachandran et al. 2023, 2024), and WR stars (white stars for apparent single, gray stars for binaries Crowther et al. 2002; Hainich et al. 2014, 2015; Tramper et al. 2015; Shenar et al. 2016, 2019). The empirical HD limit (Humphreys & Davidson 1979) is indicated by the hatched regions. Solid black and gray lines show stellar evolution tracks of the Model no-Eject until core hydrogen depletion. For clarity only the tracks with initial masses Mini = 40 M⊙, 63 M⊙, 100 M⊙, 158 M⊙, and 251 M⊙ are labeled and colored in black. Dotted black lines mark the H-ZAMS and He-ZAMS. Contours in the background highlight the degree of inflation covering R/Rcore = 1 – 3. Inflations of R/Rcore = 1.6, 2.2, and 2.8 for stars during hydrogen burning, and R/Rcore = 1.05, 1.1, and 1.15 for the He-ZAMS models are marked by colored lines and numbers.

Since it is uncertain whether this criterion would also hold for helium stars, such as the classical WR stars, and because for those, we need to employ the superadiabatic reduction method from MESA, we also calculated the extent of inflation for stars on the He-ZAMS, which we indicate in Fig. 1. From inspection, it is clear that, regardless of metallicity, the maximum inflation for helium stars appears to be R/Rcore = 1.1, as no classical WR stars are observed above this threshold (see also Figure 18 of Köhler et al. 2015).

Based on these empirical limits, we defined the maximum allowed inflation in our models as

R max / R core = 2 α ( X H ) , Mathematical equation: $$ \begin{aligned} R_{\rm max}/R_{\rm core} = 2\,\alpha (X_{\rm H}), \end{aligned} $$(1)

with a scaling factor, α, depending on the surface hydrogen abundance, XH, as

α ( X H ) = { 1 X H > = 0.2 2.125 X H + 0.575 X H < 0.2 . Mathematical equation: $$ \begin{aligned} \alpha (X_{\rm H}) = {\left\{ \begin{array}{ll} 1&X_{\rm H}> = 0.2\\ 2.125\,X_{\rm H}+0.575&X_{\rm H} < 0.2 . \end{array}\right.} \end{aligned} $$(2)

As shown by Petrovic et al. (2006), mass loss with rates above a critical value can reduce the size of the inflated envelopes. Therefore, to approximate the mass lost during an ejection episode, we adopted the critical mass-loss rate (Ṁeject) of Petrovic et al. (2006) needed to prevent a star from further inflating as

M ˙ eject = 4 π R core 2 ρ min G M R core , Mathematical equation: $$ \begin{aligned} \dot{M}_{\rm eject} = 4\pi \,R_{\rm core}^2\,\rho _{\rm min}\,\sqrt{\dfrac{\mathrm{G} M}{R_{\rm core}}}, \end{aligned} $$(3)

where M is the total mass of the star and ρmin is the minimum density in the inflated layers. To avoid extreme jumps in the mass-loss rates, we increased the mass-loss rate linearly between R/Rcore = α(XH)⋅1.9 and Rmax/Rcore. For illustrative purposes, we display in Fig. E.1 the stellar evolution tracks calculated with the Eddington-limit induced mass ejections and indicate their ejection episodes. One can see that the starting points of the ejection episode in the models overlap with the observed LBVs in the LMC and SMC.

We performed test calculations where we varied the inflation threshold for the maximum inflation by 50%. While this slightly alters the number of models that exceed the HD limit and the maximum WR luminosity, it does not significantly affect the total time spent during the OB or WR phase of a specific model. This gave us confidence that our chosen approach allows for the impact of mass ejections, such as during an LBV phase, on stellar populations to be studied without being overly sensitive to the choice of the threshold parameters.

We also note that the critical mass-loss rate, and thus a star’s subsequent evolution, is only weakly sensitive to the choice of stellar-wind prescription during hot-star evolutionary phases, since typical variations among commonly used OB- or WR-star mass-loss recipes primarily affect the integrated mass loss but are not high enough to significantly alter the local density structure and thus prevent the effect of envelope inflation, especially at low metallicity. However, when adopting different mass-loss prescriptions, the empirical limits used to define the maximum allowed degree of inflation may need to be recalibrated.

2.3. Model variations

Throughout this paper, we investigate how different choices of the input physics affect the synthetic stellar populations of the SMC and LMC. To do so, we created the following sets of stellar evolution models:

  • Model Eject: This is our standard model. It uses the input physics as described in the previous sections, including the modeling of Eddington-limit induced mass ejections. These models assume a constant star-formation history (SFH). We created an additional synthetic population for the SMC called SMC-SFH, where we assume that star formation ended 3.5 Myr ago to match the observed population better.

  • Model no-Eject: This model has the same setup as detailed in Sect. 2.1, but without the Eddington-limit induced mass ejections.

  • Model Binary: This model is a combination of our Model Eject with the LMC binary grids from Marchant (2016) and Pauli et al. (2022). We provide more details in Sect. 2.4.2.

  • Model RSG: Within this model, we switch to the RSG mass-loss recipe of Nieuwenhuijzen & de Jager (1990) that is frequently employed in stellar evolution models. Note that these mass-loss rates are weaker than those from (Yang et al. 2023). We still use our Eddington-limit induced mass ejection formalism here.

  • Model sc10: This model has the same setup as in Model Eject, but we changed the efficiency of semiconvective mixing to αsc = 10.

A discussion of including or excluding our Eddington-limit induced mass ejections as well as the role of binaries are provided in Sects. 4.1 and 4.2, while the impact of choosing a different physical setup is discussed in Appendix D.

2.4. Population synthesis

2.4.1. Single stars

To test our physical setup, we create synthetic populations of single stars per initial metallicity and per model variation, using the different single-star grid of detailed MESA models (see Sects. 2.1 and 2.3). The resulting synthetic populations are compared to the observed populations of the LMC and SMC. For our each of our single-star populations, each of the 33 detailed single-star models of a grid is assigned a statistical weight based on its formation probability and the duration of each time step. The probability of forming a star with a specific mass is determined by the initial mass function (IMF; Salpeter 1955). In this work, we adopt a Kroupa (2001) IMF and, for simplicity, a constant star formation rate (SFR) of SFRLMC = 0.1 M⊙ yr−1 for the LMC (Harris & Zaritsky 2009, corrected for a Kroupa IMF) and SFRSMC = 0.05 M⊙ yr−1 for the SMC (Rubele et al. 2018). The time spent in a specific evolutionary phase is determined by the detailed evolutionary models and does not rely on any external assumptions. More details on how the weights for individual models are calculated can be found in Appendix B.1.

2.4.2. Binary stars

Most massive stars are born in binary or higher-order multiple systems (Sana et al. 2012, 2013, 2014, 2025; Dunstall et al. 2015; Frost et al. 2025; Villaseñor et al. 2025). If they have sufficiently close separations, the components interact during their lifetimes. Calculating a large grid of binary models spanning the full initial period and initial mass ratio parameter space is beyond the scope of this paper. However, to test the impact of binarity, we combine our single-star grids with already calculated MESA binary evolutionary grids at LMC metallicity. For the initial primary mass range from M1,  ini ≈ 10 M⊙  −  25 M⊙, we employ the binary evolutionary grid presented in Marchant (2016), and for the range M1,  ini ≈ 28 M⊙  −  89 M⊙, we use the grid from Pauli et al. (2022), which reproduces the observed luminosity distribution of the LMC WR population. These grids have the same mass spacing (Δlog(Mini/M⊙) = 0.05) as our single-star grids. The models from Marchant (2016) span an initial period range from log(Porb,  ini/d) = 0.15  −  3.5, while the more massive star models from Pauli et al. (2022) that can evolve to larger stellar radii span initial periods from log(Porb,  ini/d) = 0.15  −  4.0. Both grids cover the full initial period space where binary interactions occur. We use a period spacing of Δlog(Porb,  ini/d) = 0.05. Furthermore, both grids cover initial mass ratios in the range from qini = 0.25  −  0.95, with a spacing of Δqini = 0.05. In total, these binary grids encompass about 18 000 detailed MESA models. To model mass-transfer, both binary model grids employ the “contact” scheme of MESA, allowing, as the name suggests, for contact phases. Furthermore, it is assumed that the accretion of mass is rotationally limited, leading to quasi-conservative mass-transfer in short-period systems and non-conservative mass-transfer in wide-period systems. Note that these grids use some different physical assumptions compared to our single-star models, such as mass-loss rates, and thus, detailed comparisons should be made with caution.

To create the binary population, we assume a binary fraction of 100% and initial orbital periods ranging from Porb,  ini = 1.4 d–100 yr (log(Porb,  ini/d) = 0.15–4.55. From the binary grids of Marchant (2016) and Pauli et al. (2022), we know that binary interactions, such as mass exchange, occur only in systems with Porb,  ini ≲ 3000 d, while stars in systems with longer periods evolve effectively as single stars. Following our assumed initial orbital period distribution, this translates to an initial interacting binary fraction of ≈70%, which aligns with observations of massive stars (Sana et al. 2012, 2025; Villaseñor et al. 2025). For all models that do not undergo binary interactions, we use our single-star models, meaning that for the synthetic binary population 40 000 different model combinations were used. Since our Eddington-limit induced mass ejections set an upper limit on a star’s radius, we check if a primary star’s model at a given initial mass and initial orbital period exceeds the radius of the corresponding single star model. If the primary’s radius grows larger, we assume that the star will be stripped by the Eddington-limit induced mass ejections rather than by mass transfer, and treat the evolution of the system as two independent single stars. For binary systems with initial masses beyond the range covered by the binary models, we assume that if the stars fit within the orbit upon their formation, their evolution can be approximated by single-star models. At these high masses, binary interactions have a less important role, as the maximum radius of a star with Mini = 100 M⊙ is limited to Rmax ≈ 58 R⊙ due to the Eddington-limit induced mass ejections, which prevents most mass-transfer phases.

As with the single-star population, each binary model is assigned a statistical weight. The initial primary masses are assumed to follow the Kroupa (2001) IMF, and we adopt the same SFR as for the single star population of the LMC. For simplicity, we assume a flat mass-ratio distribution (i.e., consistent with the distribution reported by Sana et al. 2012). The initial orbital period distribution is proportional to (log Porb,  ini)−0.55, following Sana et al. (2012). More details on the calculation can be found in Appendix B.2.

The binary grids used here assume several criteria for mass-transfer stability, such as L2 overflow during the contact phase, inverse mass-transfer, a limit on the maximum accretion onto the companion based on their ability to expel non-accreted material from the system, and when the mass-transferrate exceeds a high limit, of 0.1 M⊙ yr−1 (for more details, see Marchant 2016; Pauli et al. 2022). If mass transfer is flagged as unstable, we assume that the two stars will merge. If both stars are still on the main sequence (MS+MS merger), we calculate the mass of the merger as Mmerger,  ini = (1 − Φ)(M1,  current + M2,  current), where the fraction of mass lost during the merger is Φ = 0.3q/(1 + q)2 (Glebbeek et al. 2013; Wang et al. 2022). The merger product is then approximated by the closest model with the respective initial mass from our single-star grid. If one of the stars has already exhausted H in its core, we assume that the merger product should also have a helium core and an H-rich envelope (He-MS merger). To estimate the mass of the merger in such cases, we assume that during the potential common envelope phase, half of the envelope of the core-He burning star is lost. The mass of the merger is then Mini,  merger = (3M1 − M1,  he − core)/2 + M2, and we select the nearest single-star model from our grid after core hydrogen depletion. This is a simplification, and in reality, these mergers should have overmassive H-rich envelopes that would end up as blue supergiants.

To assess the number of stars in our binary population and compare them with the observed population, we only count the primary stars, the secondaries where the primary star has already died, and the mergers. Secondaries with a primary companion are assumed to be outshone by the companion and should also not be included in the number of observed stars.

2.5. Definitions of distinct stellar groups

For a consistent comparison of observations to the synthetic populations, we apply selection criteria to both populations to differentiate between the various evolutionary phases. The evolutionary phases considered in this work include WR stars, along with their multiple subtypes, as well as O-type stars, very massive stars (VMSs), and supergiants that are blue (BSG), yellow (YSG), and red (RSG). Figure 2 provides an overview of the regions occupied by the different stellar groups. A detailed explanation on how stars are assigned to each group is provided in Appendix C.

Thumbnail: Fig. 2. Refer to the following caption and surrounding text. Fig. 2.

Sketch of the regions in the HRD where the different stellar groups are located. For reference, in the background, the stellar content of the LMC is shown. The symbols have the same meaning as in Fig. 1.

3. Results

In this section, the properties and characteristics of the synthetic population using our new Eddington-limit induced mass ejection formalism, namely Model Eject, are described. A summary of the number of stars during the distinct evolutionary phases of the observed and all synthetic populations is provided in Table 1. Figure 3 shows HRDs comparing the observed populations of the LMC and SMC with the synthetic populations generated using our Model Eject.

Thumbnail: Fig. 3. Refer to the following caption and surrounding text. Fig. 3.

Hertzsprung-Russell diagrams of the massive star population (symbols) in the LMC (left) and SMC (right) compared to the synthetic population of the Model Eject (contours). A selection of the same stellar evolution models is shown as solid lines, color coded by their surface abundances. Small black dots on the evolution tracks mark equidistant timesteps of Δt = 30 000 yr. The hatched region indicates the empirical HD limit (Humphreys & Davidson 1979). The different symbols mark different stellar groups and have the same sources as mentioned in Fig. 1. In the SMC’s population, we added the WO star DR1 from IC 1613 (Tramper et al. 2013), a galaxy that has an SMC-like metallicity and only this WR star.

Table 1.

Summary of the number of stars in a specific stellar group for the different stellar populations.

3.1. The LMC population

For the LMC, the synthetic single-star population contains about 2650 O-type stars, in agreement with the estimated 2780 O-type stars in the LMC (see footnote in Table 1). The predicted 235 RSGs are also in good agreement with the 262 observed RSGs. Figure 4 compares the predicted and empirical luminosity distributions of BSGs, YSGs, and RSGs in both the LMC and SMC. Not only is the total number of RSGs well matched, but their luminosity distribution is also closely reproduced.

Thumbnail: Fig. 4. Refer to the following caption and surrounding text. Fig. 4.

Empirical luminosity distributions for the different cool supergiants (diamonds) in the LMC (left) and SMC (right) compared to the predictions from the Model Eject (hatched bins). The error margins represent the upper and lower bounds of the 84.13% confidence interval of a Poisson distribution, accounting for small numbers (Israel 1968). The sources of the observed stars are the same as in Fig. 1.

From the HRD shown in Fig. 3, one can see that stars with Mini ∼ 25 M⊙ – 50 M⊙ lose their H-rich envelopes during the RSG phase and evolve into WR stars. In these models, after the RSG phase, when most of the H-rich envelope has been removed and the envelope switches from convective to radiative, the inflation criterion for the mass ejections is triggered (see the 28 M⊙ and 40 M⊙ tracks shown in Figs. E.1). This outburst typically removes less than one solar mass, and occurs in a region close to the S Dor instability strip where the faintest LBVs are observed. Martin et al. (2025) examined the spatial distribution of LBVs in the LMC and report that the faint LBVs in the LMC are associated with evolved stars, and that they are most likely core-He burning, which is in agreement with our model predictions.

Models with Mini ≳ 50 M⊙ undergo an ejection episode already during core-H burning, preventing them from crossing the HD limit. Upon entering this phase, the models fall out of thermal equilibrium and reach mass-loss rates of log(Ṁeject/(M⊙ yr−1)) ≈ −2 to −3. These high rates quickly stop the stellar models from inflating, such that they can regain thermal equilibrium. Following the initial rapid mass ejection, the models remain close to their inflation limit (R/Rcore = 1.9) with lower but still substantial rates of log(Ṁ/(M⊙ yr−1)) ≈ −4.5, comparable to luminous WNh stars located to the right of the ZAMS. During the mass ejection episode, large fractions of the H-rich envelope are shed. Depending on initial mass, ≈20 M⊙  −  300 M⊙ can be lost.

After core-H depletion, the models expand again (regardless of whether or not MESA’s superadiabatic reduction method has been used), triggering a second mass ejection episode. This leads to an intense but short-lived mass-loss episode in which the remaining H-rich envelope (less than a few solar masses) is removed, producing hot H-free WN stars that soon evolve into WC stars. While this outburst phase may be too brief to be observed directly, it should leave behind a nebula around the WR star. Interestingly, three WR stars in the LMC are known to have nebulae (Stock et al. 2011): one is a single hot early-type WN, one a single WC, and one a WC in a binary. Our model predictions agree with the WR temperatures and luminosities. However, there is no proof that these stars have gone through an Eddington-limit induced mass ejections, as the material could also be leftovers from a RSG or common envelope phase.

Deficiencies of our single-star population are the lack of OB supergiants near the terminal age main-sequence as well as the BSG and YSG populations when compared to the observed ones. Sana et al. (2013) reported a lack of binaries in the O9.7 I, suggesting that they could be post-interaction binaries or merger products. Alternatively, the BSG might be explainable when using different assumptions on mixing (Gilkis et al. 2021). We discuss this in more detail in Sect. 4 and Appendix D. At the high-mass end, the predicted number of 28 VMSs differs by a factor of two from the observed value of 13. At such high masses, the population is very sensitive to the assumed form of the IMF and SFH. Small changes in these assumptions can easily account for the difference of the factor of two.

The synthetic single-star population of the LMC contains about 132 WR stars, which is close to the observed number of 122 WRs. The different WR subtypes serve as an excellent probe to test the stellar physics used in models, as they are sensitive to the assumptions made about mass-loss and mixing. Figure 5 shows the luminosity distributions of the WR subtypes in the LMC and SMC. Our population covers roughly the luminosity distribution of the observed WR stars and can explain the so far unexplained faintest single WR stars (Shenar et al. 2020a, and references therein). There are only small discrepancies: a few too many H-poor WN stars around log(L/L⊙)≈6.1, which originate from the core hydrogen-burning WN population after the Eddington-limit induced mass ejections, and too few bright H-free WNs around log(L/L⊙)≈5.8, which originate from stars experiencing mass ejections. Both of these features, as well as the overabundance of VMSs in the synthetic population compared to observations, might suggest that our assumptions on the IMF or SFH are too simplistic and thus lead to such artifacts. The low number of predicted fainter H-free WN stars in the LMC population may be related to our choice of when to transition from the hot star mass-loss recipe to the WR mass-loss recipe. However, shifting this transition to lower Γe values would lead to an overproduction of WN stars in the SMC, which is inconsistent with observations. Another crucial factor is the impact of binary interactions, which can lead to a more efficient removal of the H-rich envelope, making it easier for the WN stars to become H-free. Pauli et al. (2022) report that binary interactions could be the primary mechanism behind the WN and WC populations in the LMC. We discuss this further in Sect. 4.2. For the WC/WO population in the LMC, our models predict that only stars that undergo a mass ejection episode evolve into WC/WO-type stars. Note that the brightest WC/WO stars in the LMC are most likely binaries, leading to too high luminosities when applying the bolometric correction (Bartzakos et al. 2001a,b; Pauli et al. 2022).

Thumbnail: Fig. 5. Refer to the following caption and surrounding text. Fig. 5.

Empirical luminosity distributions for the different WR subtypes (diamonds) in the LMC (left) and SMC (right) compared to the predictions from the Model Eject population (hatched bins). The error margins mark the 84% confidence intervals of a Poissonian distribution accounting for small numbers (Israel 1968). In the synthetic population, we differentiate between core H-burning (light dotted bins) and core He-burning stars (bold dotted bins). For definitions of the individual WR subtypes, refer to Appendix D. The Of/WN stars are not included in the H-rich WR subtype, as they do not meet the criteria used in our models to identify the WR phase. For the LMC’s WC/WO population, not all stars have been studied in detail using stellar atmosphere models, and for a large set of the WCs, we had to estimate their luminosity from their V-band photometry (see Section 4 of Pauli et al. 2022). Bins containing only luminosities derived from a detailed spectral analysis are marked by dark blue diamonds, those only containing luminosities derived from the V-band photometry are marked by light blue diamonds, and those containing both are marked by half dark and half light blue diamonds. Since some of these stars might be binaries, the luminosities estimated from the V-band should be treated with caution and considered as upper limits. In the SMC’s WC/WO population, we added the WO star DR1 from IC 1613 (Tramper et al. 2013).

3.2. The SMC population

The HRD for SMC metallicity, shown in Fig. 3, reveals a somewhat different picture compared to the LMC. The synthetic population predicts about 1800 O-type stars, roughly three times as many as are observed. This discrepancy primarily arises from the absence of stars on the first half of the main sequence for initial masses above Mini ≳ 25 M⊙. This phenomenon was studied in more detail by Schootemeijer et al. (2021), who conclude that changes in the SFH or the IMF alone cannot explain the observed discrepancy. They conclude that another process is likely at work, and propose that the absent stars might be hidden in their birth clouds.

The RSG population, mainly originating from stars with Mini ≲ 25 M⊙, contains 131 stars in the synthetic population, closely matching the observed 124 RSGs and their luminosity distribution (Fig. 4). At SMC metallicity, OB-star winds are weaker, leading to less mass removal from the H-rich envelope during core-H burning. Consequently, more envelope mass must be lost during the RSG phase, raising the minimum initial mass required to form WR stars to Mini ≳ 28 M⊙. Interestingly, the luminosity of our faintest WR star aligns with the faintest observed single WR in the SMC (Hainich et al. 2015; Schootemeijer et al. 2024), which previous single-star models could not explain.

As in the LMC, models with Mini ≳ 50 M⊙ experience a mass ejection episode, preventing most stars from exceeding the HD limit. Envelope inflation is Z-dependent, causing stars to expand less rapidly at lower metallicity. This can be seen by a few models that surpass the HD limit for less than a few thousand years. In total, one star is predicted beyond the HD limit, consistent with observations. Given the weak Z dependence of the Eddington-limit induced mass ejection prescription, it might still be possible that below SMC metallicity stars can be found beyond the HD limit. However, this number should still be much smaller compared to models that ignore Eddington-limit induced mass ejections (see Sect. 4.1). Furthermore, the synthetic SMC population predicts 17 VMSs, whereas none have been observed. Although we are working within a parameter regime where low-number statistics are relevant, the absence of any observed stars that should be the brightest in the galaxy is concerning.

Examining the WR subtype luminosity distributions, displayed in Fig. 5, the synthetic population generally reproduces the observed distributions of the H-rich and H-poor WNs but fails at the H-free WNs and WC/WO stars. In general, it also overpredicts the total number of WR stars by a factor of three, which is mostly due to the prediction of too many bright (massive) stars. From the H-rich WN population, two key points emerge. First, the synthetic population tends to predict too many bright WR stars. Second, our models struggle to reproduce the faint H-rich WN stars. This discrepancy arises from our criterion for defining the WR phase. If one looks at the Mini = 28 M⊙ track in Fig. 3, we see that it passes through the region where H-rich WN stars are observed and has the fitting surface H-abundance. For the H-poor WN population, the synthetic model matches the shape of the observed population well, although it is slightly overabundant, especially at the high luminosity end, which is populated by stars going through a mass ejection episode. The additional mass loss efficiently strips the H-rich envelope in the most massive stars, and predicts 4 H-free WNs and 4 WC/WO stars. While the SMC does not contain any H-free WNs, it does host one luminous WO-type star in a binary system. The position of this star matches the upper luminosity range covered by our evolutionary tracks. Interestingly, the evolutionary tracks shown in Fig. 3 indicate that the WC/WO phase stars extend down to luminosities of about log(L/L⊙)≈5.5. IC 1613, a small galaxy with SMC-like metallicity (Tautvaišienė et al. 2007), hosts one WR star (DR1), a WO-type star with a luminosity of log(L/L⊙) = 5.68 ± 0.10 (Tramper et al. 2013). So far, this star has not been explained by single-star evolution models using mass-loss recipes that match typical observational data. We believe that our Eddington-limit induced mass ejection implementation offers a new possibility to explain this star’s evolutionary history.

3.3. The SMC population with changed SFH

The overabundance of O-type stars, VMSs, and WR stars, each arising from stars with Mini ≳ 28 M⊙, in the synthetic SMC population, may suggest that our assumption of a constant SFH could be erroneous. To explore this further, we created a new synthetic population, assuming that star formation in the SMC stopped 3.5 Myr ago. We chose this stopping condition as it matches well the shape of the observed stars in the SMC and is roughly consistent with the age estimates of the bright O-stars in the young cluster NGC 346 (Rickard & Pauli 2023). The HRD and WR luminosity distribution for this new population are shown in Fig. 6.

Thumbnail: Fig. 6. Refer to the following caption and surrounding text. Fig. 6.

Hertzsprung-Russell diagram of the SMC population (left) and luminosity distributions of the WR subtypes in the SMC (right) compared to the model predictions of our Model Eject, where star-formation stopped 3.5 Myr ago. In the SMC’s WC/WO population, we added the WO star DR1 from IC 1613 (Tramper et al. 2013). The symbols and colors are the same as in Figs. 3 and 5.

From the HRD of the synthetic population, illustrated in the left panel of Fig. 6, it is immediately apparent that the age cut removes all VMSs (their evolutionary lifetime is much shorter than our age cut), bringing the synthetic population into better agreement with observations. Additionally, the region near the ZAMS, for stars with Mini ≳ 25 M⊙, is no longer populated, and the contours of the synthetic population now match the observed distribution much more closely. Our synthetic population now contains about 890 O-type stars, which is closer to the estimated number of 500 O-type stars, but still not in full agreement. Upon closer inspection of the HRD shown in Fig. 6, an absence of O-stars in the spectroscopic sample can be seen near the predicted turn-off point located around log(T/kK)≈4.5 and log(L/L⊙)≈5.7. However, this absence is not present in the photometric Gaia sample of Schootemeijer et al. (2021, their Figure 5), suggesting it might be a coincidence that these stars are missing in our spectroscopic sample. Alternatively, this lack of stars could reflect a binarity effect, where the most luminous stars are binary products, namely post-interaction binaries or mergers. If it were a binary effect, it could also imply that our SFH is still too active and thus might lead to an overestimation of the O-type stars.

The WR luminosity distribution, shown in the right panel of Fig. 6, changed substantially compared to the constant-SFR case. The luminous H-rich and H-poor WNs are no longer present in our synthetic population. With the adjusted SFH, the distribution of H-poor stars is well matched, though we are now missing a few bright WR stars. Interestingly, due to the age cut, only one H-free WN star, which is created as a result of the Eddington-limit induced mass ejections during the evolution of the most massive stars in our model grid, is expected to exist in the SMC. Furthermore, our models still predict three WC/WO stars within the observed luminosity range of the WO star in the SMC. Note that this observed WO star is in a close binary, which might indicate that it is the result of binary interaction rather than Eddington-limit induced mass ejections. However, DR1, the apparently single WO star in IC 1613, still fits this scenario.

4. Discussion

4.1. Neglecting Eddington-limit induced mass ejections

An HRD for a single-star population computed with the same physical setup as before, but without including Eddington-limit induced mass ejections (Model no-Eject), is shown in Fig. 7. In this case, all of our massive star models evolve beyond the HD limit, which conflicts with observations.

Thumbnail: Fig. 7. Refer to the following caption and surrounding text. Fig. 7.

Hertzsprung-Russell diagrams of the massive star population (symbols) in the LMC (left) and SMC (right) compared to the synthetic population of the Model no-Eject (contours). The symbols and colors are the same as in Fig. 3.

In Fig. 8, we compare the luminosity distribution of the observed cool supergiants to the predictions of different models. Model no-Eject still reproduces the RSG population fairly well. Stars with Mini ≲ 25 M⊙ (i.e., those not affected by envelope inflation) evolve similarly in both the Eject and no-Eject models, as they never experience a mass ejection episode. Stars with Mini ≈ 25 M⊙−40 M⊙ still lose most of their H-rich envelope during the RSG phase, but without the mass ejections present in Model Eject (see Fig. E.1). Instead, they must shed these layers through their stellar winds, which is reflected in fewer WR stars.

Thumbnail: Fig. 8. Refer to the following caption and surrounding text. Fig. 8.

Empirical luminosity distributions for the different cool supergiants (diamonds) in the LMC (left) and SMC (right) compared to the predictions from the models with different physical setups (bins). The sources of the observed stars are the same as in Fig. 1.

For Mini ≳ 50 M⊙, the neglect of Eddington-limit induced mass ejections allows the stars to expand beyond the HD limit while still in core hydrogen burning, spending most of their time as YSG. Since we switch to the RSG mass-loss rate as soon as a certain fraction of the envelope becomes convective and contains recombined hydrogen, it is applied for these massive models already during the YSG phase. The mass-loss rates from Yang et al. (2023) increase rapidly with luminosity, preventing these stars from expanding and reaching temperatures as low as those of observed RSGs. The applied mass-loss rates during the cool supergiant phase effectively remove the H-rich layers from these massive stars, leaving them with surface H-abundances typically below XH ≲ 0.3 by the end of their lives. Additionally, since we use the locally working superadiabatic reduction method in MESA instead of MLT++, these stars can remain inflated even with such low surface H-abundances instead of contracting and becoming WR stars. The resulting overproduction of YSGs is evident in Fig. 8: the model predicts 23 YSGs in the LMC and 16 in the SMC beyond the HD limit (log(L/L⊙)≳5.5 (Davies et al. 2018). However, no star is seen in this region of the HRD in any of these galaxies, nor in any other galaxies with a similar metallicity (Schneider et al. 2025), underscoring this tension.

4.2. The role of binarity in stellar populations

Large multi-epoch observing campaigns in the Galaxy, LMC, and SMC have established that most massive stars belong to a close binary system (e.g., Sana et al. 2025, and references therein). In binary systems, interactions such as mass transfer or common-envelope evolution can drastically alter the stars’ subsequent evolution and ultimate fate. Pauli et al. (2022) computed an extensive grid of binary evolution models and concluded that the majority of WR stars in the LMC are likely formed via mass transfer. However, their models neglected Eddington-limit induced mass ejections and adopted the RSG mass-loss recipe from Nieuwenhuijzen & de Jager (1990) with a Z dependence, which effectively prevented most single stars from evolving into WR stars (cf. Appendix D.1). Conversely, Shenar et al. (2020a) questioned whether binary interaction is the dominant WR formation channel at low Z, noting that the observed WR binary fraction appears independent of metallicity at roughly 30%  −  40% (e.g., van der Hucht 2001; Foellmi et al. 2003; Deshmukh et al. 2024).

In this work, we assess the role of binary interactions across different stellar groups using our new Eddington-limit induced mass ejection prescription. Instead of calculating new binary model grids, which is computationally expensive, we combine our single-star models with the binary evolution models of Pauli et al. (2022) and Marchant (2016) (for more details, see Sect. 2.4.2). The resulting predictions for the number of stars in each stellar group are listed in Table 1.

4.2.1. The OB and VMS population

An HRD comparing the synthetic binary population at LMC metallicity with observations is shown in the left panel of Fig. 9. Overall, it is quite similar to the single-star population (Fig.3), apart from the appearance of a He-star population located below the WR stars. The predicted number of O stars is about 3200, a factor of 1.2 times higher than in the single-star models. This number is sensitive to both the assumed post-interaction evolution of primaries, secondaries, and merger products, as well as the way one counts the stars that are seen by an observer. Given that the difference remains within a few tens of percent, we consider our models to still provide a good representation of the LMC population. Interestingly, because the model weights are assigned according to the IMF of the primary mass and the secondaries are assumed to follow a uniform mass-ratio distribution, the inclusion of binaries reduces the predicted number of VMSs to 19, bringing the models at the high mass end into closer agreement with observations.

Thumbnail: Fig. 9. Refer to the following caption and surrounding text. Fig. 9.

Hertzsprung-Russell diagram of the LMC population (left), luminosity distributions of the WR subtypes (center), and luminosity distribution of the cool supergiants (right) compared to the model predictions of our Model Binary.

4.2.2. The WR population

The LMC’s synthetic binary population contains 120 WR stars, which is remarkably similar to the single-star case. In the central panel of Fig. 9, we compare the luminosity distribution of the observed WR subtypes to our model predictions. Different hatch patterns indicate WR stars formed via binary interaction (stable mass-transfer) or self-stripping (including non-interacting binaries and mergers). Caution is needed when interpreting WR stars brighter than log(L/L⊙)≳6.1, as in our models, these can only form through self-stripping. This limitation arises because the adopted binary grids cover initial primary masses only up to M1,  ini = 89 M⊙, meaning that interactions in more massive systems with orbits short enough to initiate mass transfer (Porb,  ini ≲ 50 d) are not properly modeled.

The luminosity distribution of H-poor WN stars is well reproduced. Interestingly, most of these stars (≈80%) originate from non-interacting binaries or mergers. In contrast, in the H-free WN population, the self-stripped fraction drops to 65%, particularly at low luminosities where binary interactions dominate. This is expected since single stars must first remove their extended H-rich envelopes before entering a WR phase, whereas binary mass transfer can efficiently strip most of the H-rich layers regardless of initial mass, allowing stars to more easily evolve into H-free WNs, especially at lower initial masses.

In the H-free WN luminosity distribution, a distinct peak at log(L/L⊙)≈5.8 appears, which is dominated by self-stripped stars and not visible in the observed LMC’s WR population. This feature is a by-product of our Eddington-limit induced mass ejection treatment. By limiting the maximum inflation, all stars with Mini ≳ 89 M⊙ end their lives in the same HRD region (cf. the 89 M⊙ and 150 M⊙ tracks in Fig. 9, left panel), leading to an overpopulation. Since binary evolution is not modeled for the most massive stars, this peak may be an artifact of our simplified approach to generating a binary population and might vanish in properly calculated binary evolution models using our setup.

For WC/WO stars, our model predicts 16 objects compared to the 26 observed. Since we are only calculating our models until core He depletion, it might be that we accidentally cut off a part of this population. However, given that core C burning takes only about 10% of the time a star spends core He burning, this would marginally increase the numbers and might only help explain the faint WO stars seen in the LMC.

Above, we have only discussed the fraction of WR stars formed through binary interactions. This is not equivalent to the binary fraction of WR stars. Schnurr et al. (2008) observed that only ≈30% of WN stars in the LMC are in binaries with P ≲ 200 d. Since we do not know the current period for all systems, we try to estimate the binary fraction using the initial orbital period. Using this approach, we obtain that 43% of the WN stars should be in a binary system with P ≲ 200 d. The discrepancy in the binary fractions, or at least a large part of it, can be explained by observational biases that were not taken into account when reporting the observed WR binary fraction. Overall, this agreement provides additional support for our model, which lifts parts of the long-standing question of how low-Z WR stars form.

4.2.3. The cool supergiant population

Binary interaction and binary evolution are often associated with mass transfer, which is considered a shortcut to forming WR stars at the expense of skipping the RSG phase and thus causes there to be fewer RSGs in a binary population. In contrast to these expectations, our binary population contains 270 RSGs, which is slightly more stars than in the single-star population. The RSGs in our population originate both from mergers and ejected secondaries. This can be seen in the right panel of Fig. 9, which illustrates the luminosity distribution of the cool supergiant population, separated into contributions from non-interacting binaries, secondaries, main-sequence mergers, and post-main-sequence mergers. The fraction of RSGs originating from mergers decreases with increasing luminosity: while ≈25% of the RSGs at log(L/L⊙)≲5.4 are merger products, this fraction drops sharply below < 5% at higher luminosities. This behavior reflects the stability criteria for mass transfer adopted in our binary grids: lower-mass donors often undergo unstable mass transfer, whereas mass transfer from higher-mass donors is typically stable, thus leading to fewer mergers at higher initial masses.

For the BSG and YSG populations, our models predict the same deficit as seen in the single-star population. This is expected, since we approximate mergers with our single-star models. However, the majority of mergers in our models (≈85%) occur after the primary has finished core-H burning. A merger of a post-main-sequence star and a main-sequence companion has a different core-to-envelope ratio, drastically impacting a star’s further evolution. He-MS mergers are expected to evolve differently from “normal” single stars and spend more or even all time as BSGs or YSGs instead of as RSGs (e.g., Podsiadlowski et al. 1990; Vanbeveren et al. 2013; Justham et al. 2014; Schneider et al. 2025). This effect could also help reduce the modest RSG excess we see in the synthetic population. This hypothesis is further supported by recent BLOeM campaign results in the SMC (Shenar et al. 2024), which indicate that nearly all BAF supergiants are apparently single (Patrick et al. 2025), a fraction significantly lower than the ones measured for their main-sequence progenitors (Sana et al. 2025; Villaseñor et al. 2025). This suggests that the BAF population may contain a large fraction of merger products.

We caution that these results and the discussion are based on several simplifying assumptions, including (i) all systems are disrupted after the primary’s death, (ii) our merger mass estimates, and (iii) the approximation of mergers by single-star tracks. Changing one of these could substantially affect the predicted cool supergiant population and its luminosity distribution.

4.3. Comparison to other LBV models

The literature on LBV mass-loss prescriptions is sparse, limiting the scope for detailed comparisons with other implementations. Grassitelli et al. (2021) performed time-dependent hydrodynamic stellar evolution calculations of a 73 M⊙ model to demonstrate, on the example of AG Car, that LBV eruptions are linked to opacity and inflated envelopes near the Eddington limit. Their approach reproduces both the quiescent and outburst phases of LBVs, but is computationally expensive, taking over 100 000 timesteps for the LBV phase alone. Nevertheless, we can roughly compare the physical properties of their Galactic model with our 71 M⊙ LMC model. In the hydrodynamic models, the transition from quiescence to outburst occurs at T = 20–25 kK, with mass-loss rates of log(ṀLBV,Gras/M⊙ yr−1) = −2.9 to −3.3 during the outburst. Following our Eddington-limit induced mass ejection prescription, our 71 M⊙ model enters the rapid mass ejection episode at T = 22 kK and reaches an average log(Ṁeject/M⊙ yr−1) ≈ −3.1, which is in agreement with the model from Grassitelli et al. (2021).

In the time-dependent hydrodynamic model, the LBV eruptions recur every 12 years. Our prescription treats the Eddington-limit induced mass ejections, and thus a LBV phase, as a continuous outflow, which may underestimate the actual time spent in the LBV state. However, because the LBV phase in our models proceeds on a timescale comparable to the thermal timescale, even a tenfold increase in duration with intervening quiescent phases would have little effect on subsequent stellar evolution and would only impact the predicted number of LBVs.

Cheng et al. (2024) used the MESA code to estimate the energy available in layers with super-Eddington luminosities that could be used to unbind the stellar envelope. They present stellar models at solar and SMC metallicity. In their prescription, the fraction of excess energy available for envelope ejection is treated as a free parameter between 0 and 1. Applying this formalism can yield mass-loss rates up to log(ṀLBV,Cheng/M⊙ yr−1) ≈ −2, comparable to our models.

Their Galactic-metallicity models using the maximum possible ejection efficiency exhibit a bimodal behavior: stars with Mini = 25 M⊙−60 M⊙ undergo eruptive mass loss during the YSG phase, preventing them from evolving into RSGs (even at luminosities where RSGs are observed) while stars with Mini ≳ 70 M⊙ lose mass too early, turning back before reaching the S Dor instability strip. In their Galactic models with lower ejection efficiencies as well as in their SMC-metallicity models, the eruptions are insufficient to remove most of the H-rich envelope, allowing the stars to remain beyond the HD limit for extended periods. This is in contrast to both our models and those of Grassitelli et al. (2021), which do not spend a significant fraction of their evolution beyond the HD limit except during Eddington-limit induced mass ejections.

5. Summary

In this work, we have presented a new mass-loss prescription to account for mass ejections during the evolution of luminous massive stars approaching the Eddington limit, such as those expected during an LBV phase. Within our models, the mass ejections are triggered when stellar models inflate their envelope as a consequence of reaching the Eddington limit beyond a threshold value. The HD limit and the WR populations of the LMC and SMC were used as benchmarks to calibrate the maximum inflation a star can have before going through a mass ejection episode, thus effectively removing inflation.

From our models, Eddington-limit induced mass ejections can be triggered during two main phases. Stars more massive than Mini ≳ 50 M⊙ enter an LBV-like phase during the main sequence. Before the mass ejection episode, they have mass-loss rates comparable to those of the early O supergiants, while after the mass ejections, they have mass-loss rates comparable to the WNh stars; both were observed in the same region of the HRD. Stars with Mini ≈ 25 M⊙−40 M⊙ evolve first into RSGs. As soon as they expel their H-rich envelope, radiation pressure dominates in the envelope, which is associated with envelope inflation. After the RSG phase, these stars experience short mass ejections in a region similar to those of the faintest observed LBVs.

Our new Eddington-limit induced mass ejection model exhibits only a slight Z dependence, as the opacity governs how much a star inflates. This Z dependence is evident when comparing the LMC and SMC populations, although it remains consistent with observed data. Due to the weak Z dependence, one would still expect that at very low metallicities, only a few stars would reside beyond the HD limit. These stars would be observable as cool supergiants rather than RSGs.

Despite their simplicity, the models using our new physical setup successfully replicate many of the features observed in stellar populations, even at low Z, and they provide a definitive improvement compared to models without Eddington-limit induced mass ejections. Within the models, the Eddington-limit induced mass ejections occur in regions of the HRD where LBVs are observed. These mass ejections explain the lack of stars beyond the HD limit, including the upper limit of RSG luminosities. Furthermore, the number of observed O-type stars, WRs, and RSGs in the LMC and SMC can be matched. When using strong Z-independent RSG winds, such as those from Yang et al. (2023), we can explain the faintest observed single WR stars; however, it is only when we use the Eddington-limit induced mass ejections that we can reproduce the observed low-Z WO-type stars in the SMC and IC 1613.

We combined our single-star models with large detailed grids of binary evolution models at LMC metallicity. According to our binary population, binary interactions still play a fundamental role in forming WR stars, especially for the H-free WN stars. However, a large fraction of the WR stars can be formed by stars in non-interacting binaries and mergers. One of the limitations of our Eddington-limit induced mass ejection model is that all stars ≳150 M⊙ converge to the same track, leading to an excess of luminous H-free WR stars. Our binary population assumes that all stars are born in binary systems, with 70% of the O-stars found in close binaries where they can interact. Interestingly, our new model predicts that about less than 40% of the WR stars are in close binary systems, which is in agreement with observations. This further underlines the importance of Eddington-limit induced mass ejections in stellar populations.

Acknowledgments

DP acknowledges financial support from the FWO in the form of a junior postdoctoral fellowship No. 1256225N. AP acknowledges support from the FWO under grant agreement No. 11M8325N (PhD Fellowship), and K209924N, K223124N, K1A4925N (Travel Grants). PM acknowledges support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 101165213/Star-Grasp), and from the Fonds Wetenschappelijk Onderzoek (FWO) senior postdoctoral fellowship number 12ZY523N.

References

  1. Aadland, E., Massey, P., Hillier, D. J., et al. 2022, ApJ, 931, 157 [NASA ADS] [CrossRef] [Google Scholar]
  2. Aguilera-Dena, D. R., Langer, N., Antoniadis, J., & Müller, B. 2020, ApJ, 901, 114 [NASA ADS] [CrossRef] [Google Scholar]
  3. Alkousa, T., Crowther, P. A., Bestenlehner, J. M., et al. 2025, A&A, 699, A314 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Antoniadis, K., Bonanos, A. Z., de Wit, S., et al. 2024, A&A, 686, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Barlow, M. J., & Hummer, D. G. 1982, IAU Symp., 99, 387 [Google Scholar]
  6. Bartzakos, P., Moffat, A. F. J., & Niemela, V. S. 2001a, MNRAS, 324, 18 [NASA ADS] [CrossRef] [Google Scholar]
  7. Bartzakos, P., Moffat, A. F. J., & Niemela, V. S. 2001b, MNRAS, 324, 33 [NASA ADS] [CrossRef] [Google Scholar]
  8. Beasor, E. R., Davies, B., Smith, N., et al. 2023, MNRAS, 524, 2460 [Google Scholar]
  9. Bernini-Peron, M., Sander, A. A. C., Ramachandran, V., et al. 2024, A&A, 692, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Bestenlehner, J. M. 2020, MNRAS, 493, 3938 [Google Scholar]
  11. Bestenlehner, J. M., Vink, J. S., Gräfener, G., et al. 2011, A&A, 530, L14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Bouret, J. C., Lanz, T., Martins, F., et al. 2013, A&A, 555, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Bouret, J. C., Martins, F., Hillier, D. J., et al. 2021, A&A, 647, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Breysacher, J. 1981, A&A, 43, 203 [Google Scholar]
  15. Bronner, V. A., Laplace, E., Schneider, F. R. N., & Podsiadlowski, P. 2025, A&A, 703, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Castro, N., Oey, M. S., Fossati, L., & Langer, N. 2018, ApJ, 868, 57 [NASA ADS] [CrossRef] [Google Scholar]
  18. Cheng, S. J., Goldberg, J. A., Cantiello, M., et al. 2024, ApJ, 974, 270 [NASA ADS] [CrossRef] [Google Scholar]
  19. Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102 [Google Scholar]
  20. Claret, A., & Torres, G. 2016, A&A, 592, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Costa, G., Shepherd, K. G., Bressan, A., et al. 2025, A&A, 694, A193 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Cox, J. P., & Giuli, R. T. 1968, Principles of stellar structure (Gordon and Breach) [Google Scholar]
  23. Crowther, P. A. 2019, Galaxies, 7, 88 [Google Scholar]
  24. Crowther, P. A., Dessart, L., Hillier, D. J., Abbott, J. B., & Fullerton, A. W. 2002, A&A, 392, 653 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Davies, B., Crowther, P. A., & Beasor, E. R. 2018, MNRAS, 478, 3138 [NASA ADS] [CrossRef] [Google Scholar]
  26. Decin, L. 2021, ARA&A, 59, 337 [NASA ADS] [CrossRef] [Google Scholar]
  27. Decin, L., Richards, A. M. S., Marchant, P., & Sana, H. 2024, A&A, 681, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Deshmukh, K., Sana, H., Mérand, A., et al. 2024, A&A, 692, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Doran, E. I., Crowther, P. A., de Koter, A., et al. 2013, A&A, 558, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Dray, L. M., Tout, C. A., Karakas, A. I., & Lattanzio, J. C. 2003, MNRAS, 338, 973 [NASA ADS] [CrossRef] [Google Scholar]
  31. Dufton, P. L., Evans, C. J., Hunter, I., Lennon, D. J., & Schneider, F. R. N. 2019, A&A, 626, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Dunstall, P. R., Dufton, P. L., Sana, H., et al. 2015, A&A, 580, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Evans, C. J., Crowther, P. A., Fullerton, A. W., & Hillier, D. J. 2004, ApJ, 610, 1021 [NASA ADS] [CrossRef] [Google Scholar]
  34. Foellmi, C., Moffat, A. F. J., & Guerrero, M. A. 2003, MNRAS, 338, 360 [NASA ADS] [CrossRef] [Google Scholar]
  35. Frost, A. J., Sana, H., Le Bouquin, J. B., et al. 2025, A&A, 701, A171 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Gautschy, A., Glatzel, W., Gautschy, A., & Glatzel, W. 1990, MNRAS, 245, 597 [Google Scholar]
  37. Georgy, C., Ekström, S., Eggenberger, P., et al. 2013, A&A, 558, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Gilkis, A., Vink, J. S., Eldridge, J. J., & Tout, C. A. 2019, MNRAS, 486, 4451 [CrossRef] [Google Scholar]
  39. Gilkis, A., Shenar, T., Ramachandran, V., et al. 2021, MNRAS, 503, 1884 [NASA ADS] [CrossRef] [Google Scholar]
  40. Glebbeek, E., Gaburov, E., Portegies Zwart, S., & Pols, O. R. 2013, MNRAS, 434, 3497 [NASA ADS] [CrossRef] [Google Scholar]
  41. Gómez-González, V. M. A., Oskinova, L. M., Hamann, W. R., et al. 2025, A&A, 695, A197 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Gräfener, G., Vink, J. S., de Koter, A., & Langer, N. 2011, A&A, 535, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Gräfener, G., Owocki, S. P., & Vink, J. S. 2012, A&A, 538, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Grassitelli, L., Langer, N., Mackey, J., et al. 2021, A&A, 647, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Hainich, R., Rühling, U., Todt, H., et al. 2014, A&A, 565, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Hainich, R., Pasemann, D., Todt, H., et al. 2015, A&A, 581, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Hainich, R., Shenar, T., Sander, A., Hamann, W. R., & Todt, H. 2017, IAU Symp., 329, 171 [Google Scholar]
  48. Harris, J., & Zaritsky, D. 2009, AJ, 138, 1243 [NASA ADS] [CrossRef] [Google Scholar]
  49. Heger, A., Jeannin, L., Langer, N., & Baraffe, I. 1997, A&A, 327, 224 [NASA ADS] [Google Scholar]
  50. Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368 [NASA ADS] [CrossRef] [Google Scholar]
  51. Higgins, E. R., & Vink, J. S. 2020, A&A, 635, A175 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Hopkins, P. F., Kereš, D., Oñorbe, J., et al. 2014, MNRAS, 445, 581 [NASA ADS] [CrossRef] [Google Scholar]
  53. Humphreys, R. M., & Davidson, K. 1979, ApJ, 232, 409 [Google Scholar]
  54. Humphreys, R. M., Weis, K., Davidson, K., & Gordon, M. S. 2016, ApJ, 825, 64 [NASA ADS] [CrossRef] [Google Scholar]
  55. Hunter, I., Lennon, D. J., Dufton, P. L., et al. 2008, A&A, 479, 541 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Ishii, M., Ueno, M., & Kato, M. 1999, PASJ, 51, 417 [NASA ADS] [Google Scholar]
  57. Israel, M. H. 1968, Caltech Space Radiation Lab. Internal Rept. No. 4 [Google Scholar]
  58. Jermyn, A. S., Tout, C. A., & Chitre, S. M. 2018, MNRAS, 480, 5427 [NASA ADS] [CrossRef] [Google Scholar]
  59. Jermyn, A. S., Bauer, E. B., Schwab, J., et al. 2023, ApJS, 265, 15 [NASA ADS] [CrossRef] [Google Scholar]
  60. Jiang, Y.-F., Cantiello, M., Bildsten, L., Quataert, E., & Blaes, O. 2015, ApJ, 813, 74 [Google Scholar]
  61. Jiang, Y.-F., Cantiello, M., Bildsten, L., et al. 2018, Nature, 561, 498 [NASA ADS] [CrossRef] [Google Scholar]
  62. Jin, H., Langer, N., Lennon, D. J., & Proffitt, C. R. 2024, A&A, 690, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Justham, S., Podsiadlowski, P., & Vink, J. S. 2014, ApJ, 796, 121 [NASA ADS] [CrossRef] [Google Scholar]
  64. Kalari, V. M., Vink, J. S., Dufton, P. L., & Fraser, M. 2018, A&A, 618, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Kee, N. D., Sundqvist, J. O., Decin, L., de Koter, A., & Sana, H. 2021, A&A, 646, A180 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Kippenhahn, R., Ruschenplatt, G., & Thomas, H. C. 1980, A&A, 91, 175 [Google Scholar]
  67. Kiriakidis, M., Fricke, K. J., & Glatzel, W. 1993, MNRAS, 264, 50 [NASA ADS] [Google Scholar]
  68. Klencki, J., Nelemans, G., Istrate, A. G., & Pols, O. 2020, A&A, 638, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Klencki, J., Istrate, A., Nelemans, G., & Pols, O. 2022, A&A, 662, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  70. Köhler, K., Langer, N., de Koter, A., et al. 2015, A&A, 573, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  72. Lamers, H. J. G. L. M., & Fitzpatrick, E. L. 1988, ApJ, 324, 279 [CrossRef] [Google Scholar]
  73. Langer, N. 1989, A&A, 210, 93 [NASA ADS] [Google Scholar]
  74. Langer, N., Fricke, K. J., & Sugimoto, D. 1983, A&A, 126, 207 [NASA ADS] [Google Scholar]
  75. Langer, N., Hamann, W.-R., Lennon, M., et al. 1994, A&A, 290, 819 [NASA ADS] [Google Scholar]
  76. Laplace, E., Bronner, V. A., Schneider, F. R. N., & Podsiadlowski, P. 2025, ArXiv e-prints [arXiv:2508.11088] [Google Scholar]
  77. Limongi, M., & Chieffi, A. 2018, ApJS, 237, 13 [NASA ADS] [CrossRef] [Google Scholar]
  78. Lovekin, C. C., & Guzik, J. A. 2014, MNRAS, 445, 1766 [NASA ADS] [CrossRef] [Google Scholar]
  79. Maeder, A. 1983, A&A, 120, 113 [NASA ADS] [Google Scholar]
  80. Magg, E., Bergemann, M., Serenelli, A., et al. 2022, A&A, 661, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Marchant, P. 2016, Ph.D. Thesis, Rheinischen Friedrich-Wilhelms-Universität Bonn [Google Scholar]
  82. Martin, J. C., Humphreys, R. M., & Davidson, K. 2025, ApJ, 994, 159 [Google Scholar]
  83. Massey, P., Neugent, K. F., Dorn-Wallenstein, T. Z., et al. 2021, ApJ, 922, 177 [NASA ADS] [CrossRef] [Google Scholar]
  84. McDonald, S. L. E., Davies, B., & Beasor, E. R. 2022, MNRAS, 510, 3132 [NASA ADS] [CrossRef] [Google Scholar]
  85. Menon, A., Ercolino, A., Urbaneja, M. A., et al. 2024, ApJ, 963, L42 [NASA ADS] [CrossRef] [Google Scholar]
  86. Mokiem, M. R., de Koter, A., Vink, J. S., et al. 2007, A&A, 473, 603 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Mukhija, B., & Kashi, A. 2024, ApJ, 974, 124 [NASA ADS] [CrossRef] [Google Scholar]
  88. Nathaniel, K., Langer, N., Simón-Díaz, S., et al. 2025, A&A, 702, A197 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Neugent, K. F., Massey, P., Skiff, B., et al. 2010, ApJ, 719, 1784 [NASA ADS] [CrossRef] [Google Scholar]
  90. Neugent, K. F., Massey, P., Skiff, B., & Meynet, G. 2012, ApJ, 749, 177 [NASA ADS] [CrossRef] [Google Scholar]
  91. Neugent, K. F., Massey, P., & Morrell, N. 2018, ApJ, 863, 181 [NASA ADS] [CrossRef] [Google Scholar]
  92. Nieuwenhuijzen, H., & de Jager, C. 1990, A&A, 231, 134 [NASA ADS] [Google Scholar]
  93. Patrick, L. R., Lennon, D. J., Najarro, F., et al. 2025, A&A, 698, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. Pauli, D., Langer, N., Aguilera-Dena, D. R., Wang, C., & Marchant, P. 2022, A&A, 667, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Pauli, D., Oskinova, L. M., Hamann, W. R., et al. 2023, A&A, 673, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Pauli, D., Oskinova, L. M., Hamann, W. R., et al. 2025, A&A, 697, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
  98. Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [Google Scholar]
  99. Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
  100. Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [Google Scholar]
  101. Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [Google Scholar]
  102. Petrovic, J., Pols, O., & Langer, N. 2006, A&A, 450, 219 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  103. Podsiadlowski, P., Joss, P. C., & Rappaport, S. 1990, A&A, 227, L9 [NASA ADS] [Google Scholar]
  104. Poniatowski, L. G., Sundqvist, J. O., Kee, N. D., et al. 2021, A&A, 647, A151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  105. Ramachandran, V., Hamann, W. R., Hainich, R., et al. 2018, A&A, 615, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. Ramachandran, V., Hamann, W. R., Oskinova, L. M., et al. 2019, A&A, 625, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  107. Ramachandran, V., Klencki, J., Sander, A. A. C., et al. 2023, A&A, 674, L12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  108. Ramachandran, V., Sander, A. A. C., Pauli, D., et al. 2024, A&A, 692, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  109. Ramírez-Agudelo, O. H., Sana, H., de Mink, S. E., et al. 2015, A&A, 580, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  110. Renzo, M., Zapartas, E., Justham, S., et al. 2023, ApJ, 942, L32 [NASA ADS] [CrossRef] [Google Scholar]
  111. Rickard, M. J., & Pauli, D. 2023, A&A, 674, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  112. Rickard, M. J., Hainich, R., Hamann, W. R., et al. 2022, A&A, 666, A189 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  113. Rubele, S., Pastorelli, G., Girardi, L., et al. 2018, MNRAS, 478, 5017 [NASA ADS] [CrossRef] [Google Scholar]
  114. Sabhahit, G. N., Vink, J. S., Higgins, E. R., & Sander, A. A. C. 2022, MNRAS, 514, 3736 [NASA ADS] [CrossRef] [Google Scholar]
  115. Sabin, L., Gómez-Llanos, V., Morisset, C., et al. 2022, MNRAS, 511, 1 [CrossRef] [Google Scholar]
  116. Salpeter, E. E. 1955, ApJ, 121, 161 [Google Scholar]
  117. Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444 [Google Scholar]
  118. Sana, H., de Koter, A., de Mink, S. E., et al. 2013, A&A, 550, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  119. Sana, H., Le Bouquin, J. B., Lacour, S., et al. 2014, ApJS, 215, 15 [Google Scholar]
  120. Sana, H., Shenar, T., Bodensteiner, J., et al. 2025, Nat. Astron., 9, 1337 [Google Scholar]
  121. Sander, A. A. C., Lefever, R. R., Josiek, J., et al. 2025, Nat. Astron., submitted [arXiv:2508.18410] [Google Scholar]
  122. Sanyal, D., & Langer, N. 2013, Massive Stars: From alpha to Omega, 97 [Google Scholar]
  123. Sanyal, D., Grassitelli, L., Langer, N., & Bestenlehner, J. M. 2015, A&A, 580, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  124. Sanyal, D., Langer, N., Szécsi, D., Yoon, S.-C., & Grassitelli, L. 2017, A&A, 597, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  125. Schneider, F. R. N., Sana, H., Evans, C. J., et al. 2018, Science, 359, 69 [NASA ADS] [CrossRef] [Google Scholar]
  126. Schneider, F. R. N., Laplace, E., & Podsiadlowski, P. 2025, A&A, 700, A253 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  127. Schnurr, O., Casoli, J., Chené, A., Moffat, A. F. J., & St-Louis, N. 2008, MNRAS, 389, L38 [Google Scholar]
  128. Schootemeijer, A., & Langer, N. 2018, A&A, 611, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  129. Schootemeijer, A., Langer, N., Grin, N. J., & Wang, C. 2019, A&A, 625, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  130. Schootemeijer, A., Langer, N., Lennon, D., et al. 2021, A&A, 646, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  131. Schootemeijer, A., Shenar, T., Langer, N., et al. 2024, A&A, 689, A157 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  132. Schootemeijer, A., Götberg, Y., Langer, N., et al. 2026, A&A, in press, https://doi.org/10.1051/0004-6361/202557675 [Google Scholar]
  133. Shenar, T., Hainich, R., Todt, H., et al. 2016, A&A, 591, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  134. Shenar, T., Sablowski, D. P., Hainich, R., et al. 2019, A&A, 627, A151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  135. Shenar, T., Gilkis, A., Vink, J. S., Sana, H., & Sander, A. A. C. 2020a, A&A, 634, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  136. Shenar, T., Sablowski, D. P., Hainich, R., et al. 2020b, A&A, 641, C2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  137. Shenar, T., Bodensteiner, J., Sana, H., et al. 2024, A&A, 690, A289 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  138. Smith, N. 2014, ARA&A, 52, 487 [NASA ADS] [CrossRef] [Google Scholar]
  139. Smith, N. 2026, Encyclopedia Astrophys., 2, 508 [Google Scholar]
  140. Stock, D. J., Barlow, M. J., & Wesson, R. 2011, MNRAS, 418, 2532 [NASA ADS] [CrossRef] [Google Scholar]
  141. Tautvaišienė, G., Geisler, D., Wallerstein, G., et al. 2007, AJ, 134, 2318 [CrossRef] [Google Scholar]
  142. Testor, G., Schild, H., & Lortet, M. C. 1993, A&A, 280, 426 [NASA ADS] [Google Scholar]
  143. Tramper, F., Gräfener, G., Hartoog, O. E., et al. 2013, A&A, 559, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  144. Tramper, F., Straal, S. M., Sanyal, D., et al. 2015, A&A, 581, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  145. Trundle, C., & Lennon, D. J. 2005, A&A, 434, 677 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  146. Trundle, C., Lennon, D. J., Puls, J., & Dufton, P. L. 2004, A&A, 417, 217 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  147. Ulmer, A., & Fitzpatrick, E. L. 1998, ApJ, 504, 200 [NASA ADS] [CrossRef] [Google Scholar]
  148. Urbaneja, M. A., Kudritzki, R. P., Gieren, W., et al. 2017, AJ, 154, 102 [NASA ADS] [CrossRef] [Google Scholar]
  149. van der Hucht, K. A. 2001, New Astro. Rev., 45, 135 [Google Scholar]
  150. Vanbeveren, D. 1991, A&A, 252, 159 [NASA ADS] [Google Scholar]
  151. Vanbeveren, D., Mennekens, N., Van Rensbergen, W., & De Loore, C. 2013, A&A, 552, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  152. Villaseñor, J. I., Sana, H., Mahy, L., et al. 2025, A&A, 698, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  153. Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  154. Vink, J. S., Mehner, A., Crowther, P. A., et al. 2023, A&A, 675, A154 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  155. Wang, C., Langer, N., Schootemeijer, A., et al. 2022, Nat. Astron., 6, 480 [NASA ADS] [CrossRef] [Google Scholar]
  156. Yang, M., Bonanos, A. Z., Jiang, B., et al. 2023, A&A, 676, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  157. Yoon, S.-C. 2017, MNRAS, 470, 3970 [NASA ADS] [CrossRef] [Google Scholar]
  158. Yoon, S.-C., & Cantiello, M. 2010, ApJ, 717, L62 [NASA ADS] [CrossRef] [Google Scholar]
  159. Zapartas, E., de Wit, S., Antoniadis, K., et al. 2025, A&A, 697, A167 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

Appendix A: Additional information on the assumed input physics

Convection is modeled using the Ledoux criterion and standard mixing length theory, following Cox & Giuli (1968), with a mixing length parameter αmlt = 1.5. Thermohaline mixing is included with an efficiency factor of αth = 1 (Kippenhahn et al. 1980). Rotational mixing in MESA is included as a diffusive process. Following Heger et al. (2000), we include the effects of dynamical and secular shear instabilities, the Goldreich-Schubert-Ricke instability, and Eddington-Sweet circulations. The efficiency factors for rotational mixing are calibrated according to Brott et al. (2011), with values set to fc = 1/30 and fμ = 0.1.

The initial chemical compositions were as follows: for H and He, we used XH = 0.7343 and XHe = 0.2594 in the LMC models, and XH = 0.7452 and XHe = 0.2522 in the SMC models. For C, N, O, Ne, Na, Mg, Al, Si, S, Cl, Ar, Ca, Cr, Fe, and Ni, we adopted the recent baseline abundances for the LMC and SMC as reported by Vink et al. (2023). For elements not specified in these abundances, we assumed they can be scaled down from the solar abundances provided by Magg et al. (2022). The resulting metallicities for our LMC and SMC grids were ZLMC = 0.006177 and ZSMC = 0.0024386, respectively. For reference, the solar metallicity value of Magg et al. (2022) is Z⊙ = 0.012823.

To properly model the nuclear reactions of massive stars, we use a nuclear network inspired by the stellar evolution code STERN (Heger et al. 2000; Jin et al. 2024). We include 53 reaction rates for the p-p chains, CNO, Ne-Na, and Mg-Al burning, which are important already during core hydrogen and helium burning in VMSs.

Appendix B: Calculating statistical weights for the population synthesis

B.1. Single star population

The total number of stars in a star population can be calculated as

N tot = ∫ 0.01 400 A ξ ( m ) d m ≈ A ∑ i = 1 n ξ ( m i ) Δ m i , Mathematical equation: $$ \begin{aligned} N_{\rm tot} = \int _{0.01}^{400}A\,\xi (m)dm\approx A \sum _{i = 1}^n\xi (m_i)\Delta m_i\,, \end{aligned} $$(B.1)

where mi is the initial mass of a stellar model, Δmi is the difference in the initial mass between the successive stellar models (which is not constant, given the logarithmic spacing of our grid), n is the maximum number of models in our grid, A the normalization factor, and ξ(m) the Kroupa (2001) IMF given by

ξ ( m ) = { ( m 0.08 M ⊙ ) − 0.3 0.01 < m < 0.08 ( m 0.08 M ⊙ ) − 1.3 0.08 ≤ m < 0.5 ( 0.5 M ⊙ 0.08 M ⊙ ) − 1.3 ( m 0.5 M ⊙ ) − 2.3 0.5 ≤ m . Mathematical equation: $$ \begin{aligned} \xi (m) = {\left\{ \begin{array}{ll} \left(\dfrac{m}{0.08\,{M_{\odot }}}\right)^{-0.3}&0.01<m<0.08\\ \left(\dfrac{m}{0.08\,{M_{\odot }}}\right)^{-1.3}&0.08\le m < 0.5\\ \left(\dfrac{0.5\,{M_{\odot }}}{0.08\,{M_{\odot }}}\right)^{-1.3}\left(\dfrac{m}{0.5\,{M_{\odot }}}\right)^{-2.3}&0.5\le m\\ \end{array}\right.}. \end{aligned} $$(B.2)

Each stellar track is split into time steps of Δt = 1000 yr and it is assumed that both galaxies had constant star-formation rates (SFRs) with SFRLMC = 0.1 M⊙ yr−1 in the LMC (Harris & Zaritsky 2009, corrected for a Kroupa IMF) and SFRSMC = 0.05 M⊙ yr−1 in the SMC (Rubele et al. 2018) within the last 27.5Myr (i.e., the age of our oldest model). The predicted number of stars within an initial mass interval [m, m + Δm] and age interval [t, t + Δt] after a simulation time of Tsim is

Δ N = A ξ ( m ) Δ m Δ t / T sim . Mathematical equation: $$ \begin{aligned} \Delta N = A\,\xi (m)\,\Delta m\,\Delta t/T_{\rm sim}. \end{aligned} $$(B.3)

To calculate the normalization factor A, we use the condition that the total mass used to form stars within the simulation time is given by

M tot = SFR · T sim = ∫ 0.01 400 m ′ A ξ ( m ′ ) d m ′ . Mathematical equation: $$ \begin{aligned} M_{\rm tot}=\mathrm{SFR} \,\cdot T_{\rm sim} = \int _{0.01}^{400}m^{\prime }A\,\xi (m^{\prime })\mathrm{d} m^{\prime }. \end{aligned} $$(B.4)

Hence, the number of stars per initial mass bin Δm per age step Δt (as long as t + Δt is less than the age of a star with the corresponding initial mass) in our simulation is

Δ N = SFR · ξ ( m ) Δ m Δ t ∫ 0.01 400 m ′ ξ ( m ′ ) d m ′ . Mathematical equation: $$ \begin{aligned} \Delta N = \dfrac{\mathrm{SFR} \cdot \xi (m)\,\Delta m\,\Delta t}{\int _{0.01}^{400}m^{\prime }\,\xi (m^{\prime })\mathrm{d} m^{\prime }}. \end{aligned} $$(B.5)

B.2. Binary star population

The total number of stars in a population consisting only of binary stars is given by

N tot = ∫ 0.01 400 A ξ ( m 1 ) d m 1 ∫ 0 1 f ( q ) d q ≈ A ∑ i = 1 n ξ ( m 1 , i ) Δ m 1 , i Δ q , Mathematical equation: $$ \begin{aligned} N_{\rm tot} = \int _{0.01}^{400}A\,\xi (m_1)dm_1\int _{0}^{1}f(q)\,dq\approx A\,\sum _{i = 1}^n\xi (m_{1,i})\Delta m_{1,i}\Delta q\,, \end{aligned} $$(B.6)

where ξ(m) is the Kroupa (2001) IMF (see Eq. B.2), q = m2/m1 is the initial mass ratio, and f(q) the assumed mass-ratio distribution. Here we assume a uniform mass-ratio distribution, yielding f(q) = q0 = 1.

Similar to the calculation of the single-star population, we can use the total mass enclosed in the stars to calculate the normalization factor by

M tot = B · SFR · T sim = A ∫ 0.01 400 ( m 1 + q m 1 ) ξ ( m 1 ) f ( q ) d q d m 1 . = 1.5 A · ∫ 0.01 400 m 1 ξ ( m 1 ) d m 1 . Mathematical equation: $$ \begin{aligned} M_{\rm tot}=B\cdot \mathrm{SFR} \,\cdot T_{\rm sim}&= A\,\int _{0.01}^{400}(m_1+q\,m_1)\,\xi (m_1)f(q)dq\,\mathrm{d} m_1. \nonumber \\&= 1.5\,A\cdot \int _{0.01}^{400}m_1\,\xi (m_1)\,\mathrm{d} m_1. \end{aligned} $$(B.7)

The SFR used in this work is calculated using theoretical isochrones, photometry, and some assumptions to correct for the intrinsic distributions, such as the IMF. To correct the SFR of the LMC (Harris & Zaritsky 2009) for an intrinsic binary fraction of 100%, we need to employ a correction factor of B = 2.

For the initial orbital periods, we assume all stars are born in binaries with Porb,  ini = 1 d–100 yr. Assuming that this range encompasses all possible binary configurations, the probability of finding a star in a binary whose initial period was in the interval [log P, log P + Δlog P] is

p ( log P , log P + Δ log P ) = ∫ log P log P + Δ log P π ( log P ′ ) d log P ′ ∫ 0 4.55 π ( log P ″ ) d log P ″ . Mathematical equation: $$ \begin{aligned} p(\log P,\log P+\Delta \log P) = \dfrac{\int _{\log P}^{\log P+\Delta \log P}\pi (\log P^{\prime })d\log P^{\prime }}{\int _0^{4.55}\,\pi (\log P^{\prime \prime })\,d\log P^{\prime \prime }}. \end{aligned} $$(B.8)

Following Sana et al. (2012), we adopt the intrinsic period distribution π(log P) = (log P)−0.55.

Analogous to the single-star case, the number of binaries formed per initial initial primary mass interval [m1, m1 + Δm1], per time interval [t, t + Δt], per initial mass ratio interval [q + Δq], and per logarithmic initial orbital period interval [log P, log P + Δlog P] is given by

Δ N = 2 · SFR ξ ( m 1 ) Δ m 1 Δ q ∫ log P log P + Δ log P π ( log P ′ ) d log P ′ Δ t 1.5 ∫ 0.01 400 m 1 ′ ξ ( m 1 ′ ) d m 1 ′ ∫ 0 4.55 π ( log P ″ ) d log P ″ . Mathematical equation: $$ \begin{aligned} \Delta N = \dfrac{2\cdot \,\mathrm{SFR} \,\xi (m_{\rm 1})\Delta m_{\rm 1}\,\Delta q\,\int _{\log P}^{\log P+\Delta \log P}\pi (\log P^{\prime })d\log P^{\prime } \,\Delta t}{1.5\int _{0.01}^{400}m_1^{\prime }\xi (m_1^{\prime })dm_1^{\prime }\int _0^{4.55}\pi (\log P^{\prime \prime })d\log P^{\prime \prime }}\,. \end{aligned} $$(B.9)

Appendix C: Detailed definitions of the distinct stellar groups

Wolf-Rayet stars exhibit optically thick winds, which lead to strong emission lines in their spectra. These emission features in the observations are used for classification. However, stellar evolution models do not output the spectrum of a star. To identify when a stellar model enters the WR phase, we follow the method described Aguilera-Dena et al. (2020) and calculate the optical depth of a stellar model using

τ = − κ M ˙ 4 π R ( v ∞ − v 0 ) ln ( v ∞ v 0 ) Mathematical equation: $$ \begin{aligned} \tau =\dfrac{-\kappa \dot{M}}{4\pi R(v_\infty -v_0)}\,\ln \left(\dfrac{v_\infty }{v_0}\right) \end{aligned} $$(C.1)

from Langer (1989), where κ is approximated by the electron scattering opacity κe = 0.2 (1 + XH) cm2 g−1, Ṁ is the mass-loss rate, R is a stellar model’s radius, v0 is the expansion velocity at the photosphere, which is on the order of 20 km s−1, and v∞ is the terminal wind velocity, which is assumed to scale with the escape velocity vesc as

v ∞ = n · v esc = n 2 G M R ( 1 − Γ ) . Mathematical equation: $$ \begin{aligned} v_\infty =n\cdot v_{\rm esc}=n\sqrt{\dfrac{2GM}{R}(1-\Gamma )}. \end{aligned} $$(C.2)

Here, Γ is the Eddington factor and n is a scaling relation. While n might vary during the evolution (e.g., see Vink et al. 2001), we assume n = 1.3, which provides a good approximation for hot O-type stars and WN-type stars (Vink et al. 2001; Gräfener et al. 2011). Following Pauli et al. (2022), who calibrated the τ-criterion to the observed population of the WR stars in the LMC, we assume that a stellar model is in a WR phase as soon as τ ≥ 0.75 and log(T/K)≥4.6 (T ≳ 40 000 K).

Wolf-Rayet stars have multiple subtypes, each characterized by emission lines from different elements, namely nitrogen, carbon, and oxygen. These subtypes are associated with distinct evolutionary phases and serve as excellent benchmarks for comparison of model predictions with observed populations. In this work, we define the following categories for the nitrogen-rich WN-type stars: H-rich WNs that have XH > 0.4, H-poor WNs with 0.05 < XH ≤ 0.4, and H-free WNs are those WR stars with XH ≤ 0.05 and XC < 0.005. While observational WC and WO-type stars are classified by their different strengths in carbon and oxygen lines, recent studies indicate that the different emission strengths are due to a temperature effect rather than an abundance effect (Aadland et al. 2022; Sander et al. 2025). Therefore, we do not differentiate between these two classes, and within our models, we define the WC/WO phase as soon as XC ≥ 0.005.

The O-star phase in our models is assigned to all stars that are on the right side of the hydrogen zero age main sequence (H-ZAMS), have log(T/K)≥4.5 (T ≳ 31 500 K), and an optically thin wind (i.e., τ < 0.5). The temperature cut corresponds to a spectral type of about O9.5. Note that we apply this temperature cut also to the observed dataset rather than using the spectral classification.

As VMSs, we define all stars with a luminosity above log(L/L⊙)≥6.5 (i.e., Mini ≳ 150 M⊙). These stars, even though extremely rare in stellar populations, are dominating the integrated spectra of actively star-forming galaxies and are therefore of particular interest when trying to mimic real stellar populations. Note that the VMSs are also included in the WR and O-type stars sample.

To investigate how our synthetic populations perform, we additionally consider the evolutionary phases of BSG, YSG, and RSG. Comparing these stars provides valuable insights into the impact of mixing efficiencies as well as the performance of the assumed mass-loss rates. The RSG phase is assigned to every synthetic star that is cooler than log(T/K)≤3.68 (T ≲ 4500 K), the YSG phase for all stars with 3.68 < log(T/K)≤3.87 (4500 ≲ T ≲ 7500 K), and the BSG phase for stars with 3.87 < log(T/K)≤4.09 (7500 ≲ T ≲ 12 500 K). Since our grids only cover stars with initial masses above 10 M⊙, they do not cover the complete cool supergiant population. As a consequence of varying physical assumptions, the maximum luminosity of the RSG of a star with an initial mass of 10 M⊙ also changes. To ensure consistency, we apply a conservative luminosity cut at log(L/L⊙)≥4.7 for cool supergiants in both the synthetic and observed populations.

Appendix D: Additional variations of the input physics

D.1. Changing the RSG mass-loss rates

One of the biggest uncertainties in modeling massive star evolution, besides the Eddington-limit induced mass ejections, is the mass lost via the winds of RSG. Comparing RSG mass-loss rates derived using different model assumptions (such as the launching mechanism of the wind) shows that the measured RSG mass-loss rate can differ by up to two orders of magnitude (e.g., Decin 2021; Yang et al. 2023; Beasor et al. 2023; Antoniadis et al. 2024; Decin et al. 2024). These differences arise partly from the method of modeling itself (e.g., see figure 13 in Antoniadis et al. 2024). When comparing different recent mass-loss prescriptions (e.g., Beasor et al. 2023; Yang et al. 2023; Decin et al. 2024), it becomes evident that most of them align within a factor of a few for luminosities above log(L/L⊙)≳5 (i.e., the potential progenitors of WR stars) and that stronger mass-loss rates compared to the “classic” Nieuwenhuijzen & de Jager (1990) are observed. Still in many publicly available stellar evolution grids, the mass-loss recipe of Nieuwenhuijzen & de Jager (1990) is commonly employed, sometimes even with a Z dependence. Here, we present our Model RSG set, where we have switched from the Yang et al. (2023) mass-loss recipe to that of Nieuwenhuijzen & de Jager (1990), but without any Z dependence.

The HRDs of the LMC and SMC populations compared to the synthetic Model RSG population are shown in Fig. D.1. Since Eddington-limit induced mass ejections are still included, stellar models with Mini ≳ 50 M⊙ continue to avoid the HD limit and evolve into WR stars. One can see that these stars only produce the most luminous WR stars in the population, as well as the majority of WC/WO-type stars. For stars with Mini ≈ 25 M⊙ – 40 M⊙, which lost their H-rich envelopes during the RSG phase in the Model Eject, most of them can no longer strip off their envelopes when employing the RSG mass-loss rates from Nieuwenhuijzen & de Jager (1990) and thus do not evolve into WR stars. This challenges stellar evolution in explaining the faintest observed single WR stars in both the LMC and SMC.

Thumbnail: Fig. D.1. Refer to the following caption and surrounding text. Fig. D.1.

Hertzsprung-Russell diagrams of the massive star population (symbols) in the LMC (left), and SMC (right), compared to the synthetic population of the Model RSG (contours). The symbols and colors have the same meaning as in Fig. 1.

In Fig. 8, we compare the luminosity distributions of BSG, YSG, and RSG in the LMC and SMC to those of the Model RSG population. It is evident that the high-luminosity bins (log(L/L⊙)≳5.5) of the RSG contain too many stars. This is because these stars never lose their H-rich envelopes, and they spend all of their core He-burning lives as RSGs, contrary to the stars in Model Eject, which spend a significant fraction of core He burning as WR stars. The number of BSG and YSG is underpredicted at low luminosities (log(L/L⊙)≲5.5) and is overpredicted in the model at high luminosities (log(L/L⊙)≳5.5). We conclude that in order to achieve a better agreement with the observed stellar populations, stronger RSG mass-loss rates are favored, as this seems to be the only way to explain the observed faintest single low-Z WR stars.

D.2. More efficient semiconvective mixing

All the aforementioned single-star models fail to reproduce the observed BSG and YSG populations in the LMC and SMC. Previous studies have shown that semiconvective mixing has a strong influence on the late evolutionary stages of massive stars, as it shortens the Hertzsprung gap, affects the He-core mass, and can lead to blue loops during which the star appears as a BSG or YSG (e.g., Schootemeijer et al. 2019; Gilkis et al. 2019; Klencki et al. 2020). In Fig. D.2, the HRD for Model sc10, which adopts the same physical setup as Model Eject but with a tenfold increase in the efficiency of semiconvective mixing, is shown. The impact on stars with initial masses Mini ≈ 14 M⊙−40 M⊙ is striking: instead of remaining in the RSG phase, these stars evolve through the expected blue loops.

Thumbnail: Fig. D.2. Refer to the following caption and surrounding text. Fig. D.2.

Hertzsprung-Russell diagrams of the massive star population (symbols) in the LMC (left) and SMC (right) compared to the synthetic population of the Model sc10 (contours). The symbols and colors have the same meaning as in Fig. 1.

The luminosity distribution of cool supergiants for this model is shown in Fig. 8. While the blue-loop evolution successfully populates both the BSG and YSG regions, most of the predicted stars are found at log(L/L⊙)≳5.5, a regime where observations show few to no objects. Moreover, the Model sc10 predicts more YSGs than BSGs, opposite to what is observed. Both discrepancies may stem from our adopted semiconvection and overshooting parameters. As shown by Schootemeijer et al. (2019), stars at SMC metallicity only undergo blue loops for specific combinations of overshooting and semiconvection. It is still unclear whether these mixing parameters, such as overshooting, may vary with mass and what this implies for the cool supergiant population (e.g., Claret & Torres 2016; Jermyn et al. 2018; Gilkis et al. 2019).

A further issue of the Model sc10 is that the stars spend substantial time as BSGs or YSGs and thus do not experience strong RSG winds. This prevents them from shedding their H-rich envelopes and evolving into WR stars, which is again conflicting with observations. Interestingly, recent multi-epoch observations of BAF supergiants in the SMC suggest that most, maybe all, of them are apparently single (Patrick et al. 2025). This raises the possibility that many observed BSGs and YSGs are in fact merger products or ejected accretors (Menon et al. 2024). It has been shown that accretors that gain significant mass are not fully rejuvenated and exhibit different core-to-envelope ratios and internal structures (e.g., Klencki et al. 2022; Renzo et al. 2023; Pauli et al. 2023), which can naturally lead to blue-loop evolution without invoking highly efficient semiconvective mixing. A similar evolutionary pathway is expected for post-main-sequence merger products (e.g., Vanbeveren et al. 2013; Justham et al. 2014). Consequently, the inability of single-star models to match the observed BSG and YSG populations may not be a major concern.

Appendix E: Additional figures

Thumbnail: Fig. E.1. Refer to the following caption and surrounding text. Fig. E.1.

Hertzsprung-Russell diagrams showing the tracks of our Model Eject grid for LMC (left) and SMC metallicity (right). The Eddington-limit induced mass ejection episodes of our models are indicated by bold gray outlines around the tracks. Black dots on the tracks indicate equidistant timesteps of Δt = 30 000 yr. For comparison, the known LBVs of the LMC and SMC are shown as triangles during quiescence (left triangle) and outburst (right triangle) (Humphreys et al. 2016; Kalari et al. 2018).

All Tables

Table 1.

Summary of the number of stars in a specific stellar group for the different stellar populations.

In the text

All Figures

Thumbnail: Fig. 1. Refer to the following caption and surrounding text. Fig. 1.

Hertzsprung-Russell diagrams of the massive star population (symbols) in the LMC (left) and SMC (right). The different symbols mark the positions of observed RSG (small red plusses, Davies et al. 2018; Yang et al. 2023), YSG (crosses, Neugent et al. 2010, 2012), LBVs (filled triangles with quiescent and outburst positions, Humphreys et al. 2016; Kalari et al. 2018), OB stars (open circles and squares, Evans et al. 2004; Trundle et al. 2004; Trundle & Lennon 2005; Hunter et al. 2008; Bestenlehner et al. 2011; Bouret et al. 2013; Urbaneja et al. 2017; Castro et al. 2018; Schneider et al. 2018; Ramachandran et al. 2018, 2019; Dufton et al. 2019; Bestenlehner 2020; Bouret et al. 2021; Rickard et al. 2022; Pauli et al. 2023; Bernini-Peron et al. 2024; Gómez-González et al. 2025; Alkousa et al. 2025), partially stripped stars (white large plusses Pauli et al. 2022; Rickard & Pauli 2023; Ramachandran et al. 2023, 2024), and WR stars (white stars for apparent single, gray stars for binaries Crowther et al. 2002; Hainich et al. 2014, 2015; Tramper et al. 2015; Shenar et al. 2016, 2019). The empirical HD limit (Humphreys & Davidson 1979) is indicated by the hatched regions. Solid black and gray lines show stellar evolution tracks of the Model no-Eject until core hydrogen depletion. For clarity only the tracks with initial masses Mini = 40 M⊙, 63 M⊙, 100 M⊙, 158 M⊙, and 251 M⊙ are labeled and colored in black. Dotted black lines mark the H-ZAMS and He-ZAMS. Contours in the background highlight the degree of inflation covering R/Rcore = 1 – 3. Inflations of R/Rcore = 1.6, 2.2, and 2.8 for stars during hydrogen burning, and R/Rcore = 1.05, 1.1, and 1.15 for the He-ZAMS models are marked by colored lines and numbers.
In the text
Thumbnail: Fig. 2. Refer to the following caption and surrounding text. Fig. 2.

Sketch of the regions in the HRD where the different stellar groups are located. For reference, in the background, the stellar content of the LMC is shown. The symbols have the same meaning as in Fig. 1.
In the text
Thumbnail: Fig. 3. Refer to the following caption and surrounding text. Fig. 3.

Hertzsprung-Russell diagrams of the massive star population (symbols) in the LMC (left) and SMC (right) compared to the synthetic population of the Model Eject (contours). A selection of the same stellar evolution models is shown as solid lines, color coded by their surface abundances. Small black dots on the evolution tracks mark equidistant timesteps of Δt = 30 000 yr. The hatched region indicates the empirical HD limit (Humphreys & Davidson 1979). The different symbols mark different stellar groups and have the same sources as mentioned in Fig. 1. In the SMC’s population, we added the WO star DR1 from IC 1613 (Tramper et al. 2013), a galaxy that has an SMC-like metallicity and only this WR star.
In the text
Thumbnail: Fig. 4. Refer to the following caption and surrounding text. Fig. 4.

Empirical luminosity distributions for the different cool supergiants (diamonds) in the LMC (left) and SMC (right) compared to the predictions from the Model Eject (hatched bins). The error margins represent the upper and lower bounds of the 84.13% confidence interval of a Poisson distribution, accounting for small numbers (Israel 1968). The sources of the observed stars are the same as in Fig. 1.
In the text
Thumbnail: Fig. 5. Refer to the following caption and surrounding text. Fig. 5.

Empirical luminosity distributions for the different WR subtypes (diamonds) in the LMC (left) and SMC (right) compared to the predictions from the Model Eject population (hatched bins). The error margins mark the 84% confidence intervals of a Poissonian distribution accounting for small numbers (Israel 1968). In the synthetic population, we differentiate between core H-burning (light dotted bins) and core He-burning stars (bold dotted bins). For definitions of the individual WR subtypes, refer to Appendix D. The Of/WN stars are not included in the H-rich WR subtype, as they do not meet the criteria used in our models to identify the WR phase. For the LMC’s WC/WO population, not all stars have been studied in detail using stellar atmosphere models, and for a large set of the WCs, we had to estimate their luminosity from their V-band photometry (see Section 4 of Pauli et al. 2022). Bins containing only luminosities derived from a detailed spectral analysis are marked by dark blue diamonds, those only containing luminosities derived from the V-band photometry are marked by light blue diamonds, and those containing both are marked by half dark and half light blue diamonds. Since some of these stars might be binaries, the luminosities estimated from the V-band should be treated with caution and considered as upper limits. In the SMC’s WC/WO population, we added the WO star DR1 from IC 1613 (Tramper et al. 2013).
In the text
Thumbnail: Fig. 6. Refer to the following caption and surrounding text. Fig. 6.

Hertzsprung-Russell diagram of the SMC population (left) and luminosity distributions of the WR subtypes in the SMC (right) compared to the model predictions of our Model Eject, where star-formation stopped 3.5 Myr ago. In the SMC’s WC/WO population, we added the WO star DR1 from IC 1613 (Tramper et al. 2013). The symbols and colors are the same as in Figs. 3 and 5.
In the text
Thumbnail: Fig. 7. Refer to the following caption and surrounding text. Fig. 7.

Hertzsprung-Russell diagrams of the massive star population (symbols) in the LMC (left) and SMC (right) compared to the synthetic population of the Model no-Eject (contours). The symbols and colors are the same as in Fig. 3.
In the text
Thumbnail: Fig. 8. Refer to the following caption and surrounding text. Fig. 8.

Empirical luminosity distributions for the different cool supergiants (diamonds) in the LMC (left) and SMC (right) compared to the predictions from the models with different physical setups (bins). The sources of the observed stars are the same as in Fig. 1.
In the text
Thumbnail: Fig. 9. Refer to the following caption and surrounding text. Fig. 9.

Hertzsprung-Russell diagram of the LMC population (left), luminosity distributions of the WR subtypes (center), and luminosity distribution of the cool supergiants (right) compared to the model predictions of our Model Binary.
In the text
Thumbnail: Fig. D.1. Refer to the following caption and surrounding text. Fig. D.1.

Hertzsprung-Russell diagrams of the massive star population (symbols) in the LMC (left), and SMC (right), compared to the synthetic population of the Model RSG (contours). The symbols and colors have the same meaning as in Fig. 1.
In the text
Thumbnail: Fig. D.2. Refer to the following caption and surrounding text. Fig. D.2.

Hertzsprung-Russell diagrams of the massive star population (symbols) in the LMC (left) and SMC (right) compared to the synthetic population of the Model sc10 (contours). The symbols and colors have the same meaning as in Fig. 1.
In the text
Thumbnail: Fig. E.1. Refer to the following caption and surrounding text. Fig. E.1.

Hertzsprung-Russell diagrams showing the tracks of our Model Eject grid for LMC (left) and SMC metallicity (right). The Eddington-limit induced mass ejection episodes of our models are indicated by bold gray outlines around the tracks. Black dots on the tracks indicate equidistant timesteps of Δt = 30 000 yr. For comparison, the known LBVs of the LMC and SMC are shown as triangles during quiescence (left triangle) and outburst (right triangle) (Humphreys et al. 2016; Kalari et al. 2018).
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.