© The Authors 2023

1. Introduction

Numerous difficulties arise when stars are modeled. The treatment of convection has long been one of these difficulties. The classical way of dealing with convection in one-dimensional (1D) stellar evolution codes is the following: (1) find the boundaries of the convective zone by applying an instability criterion, (2) empirically increase the size of (some of) the convective zones (a process often called “overshooting”), and (3) compute a thermal gradient to be applied inside the convective zone. For the first point, two criteria are widely used in the literature to determine the stability of a thermally stratified medium with the gravity as the restoring force: the so-called Schwarzschild criterion, which writes (for stability against convection)

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{\text{rad}}<{\mathrm{\nabla }}_{\text{ad}},\end{array}$$(1)

where ${\mathrm{\nabla }}_{\text{rad}}={\left(\frac{\text{d}lnT}{\text{d}lnP}\right)}_{\text{rad}}$ is the radiative thermal gradient, and ${\mathrm{\nabla }}_{\text{ad}}={\left(\frac{\text{d}lnT}{\text{d}lnP}\right)}_{\text{ad}}$ is the adiabatic thermal gradient; and the Ledoux criterion, which writes (also for stability)

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{\text{rad}}<{\mathrm{\nabla }}_{\text{ad}}+\frac{\phi }{\delta }{\mathrm{\nabla }}_{\mu },\end{array}$$(2)

where $\phi ={\left(\frac{\partial ln\rho }{\partial ln\mu }\right)}_{T,P}$, $\delta =-{\left(\frac{\partial ln\rho }{\partial lnT}\right)}_{T,\mu }$, and ${\mathrm{\nabla }}_{\mu }=\left(\frac{\text{d}ln\mu }{\text{d}lnP}\right)$. It is possible that some layers inside the star are Ledoux-stable but Schwarzschild-unstable. In case of a thermally dissipative medium, however, Kato (1966) showed that Eq. (2) is reduced to Eq. (1) due to the onset of oscillatory convection. For this reason, the Schwarzschild criterion is most of the time preferred to the Ledoux criterion in stellar evolution computations. Another option is to apply a partial mixing of the layers encompassed by the two criteria. This is called “semiconvection” (see Langer et al. 1983). In this framework, models computed with the Ledoux criterion correspond to a totally inefficient semiconvection, while models computed with the Schwarzschild criterion correspond to an infinitely efficient semiconvection.

Point 2 was introduced in stellar evolution codes as a necessity to reproduce some observational features, such as the width of the main sequence (MS) of open clusters (see e.g., Maeder & Mermilliod 1981; Maeder & Meynet 1989; Martinet et al. 2021). Different implementations for this additional mixing (co)exist in stellar evolution codes: penetrative overshoot (e.g., Maeder & Mermilliod 1981; Zahn 1991), diffusive overshoot (e.g., Freytag et al. 1996; Herwig et al. 2006), entrainment (Meakin & Arnett 2007; Cristini et al. 2017; Scott et al. 2021), or more recently deduced from 2D fully compressible time-implicit simulations (Baraffe et al. 2023). A thorough discussion of these various implementations and their link to hydrodynamics simulation can be found in Viallet et al. (2015).

Finally, point 3 requires the computation of a thermal gradient to be applied inside the convective regions. For deep convection, the timescale for radiative exchanges between the convective cells and the environment are longer than the convective instability timescale, and the process is very close to being adiabatic, thus the adiabatic gradient ∇_ad can be applied. In convective regions closer to the surface, more complex theories need to be used. A usual choice in this case is the mixing-length theory (Böhm-Vitense 1958).

To obtain a complete picture of convection in stellar interiors, 3D hydrodynamics simulations are required. Convection in stars is very turbulent: the estimated Reynolds number is as high as 10⁹ (Arnett & Meakin 2016). The mass flux is therefore highly asymmetric, requiring a multi-dimensional treatment. Efforts to do this have been made during the past decade for various physical environments (e.g., Meakin & Arnett 2007; Cristini et al. 2017; Woodward et al. 2015; Jones et al. 2017; Horst et al. 2021). However, all of these simulations were limited to a small region and/or a short time compared to the size and lifetime of a star. These results can nevertheless be used to build new algorithms that can obtain a better agreement between the 1D stellar evolution codes and the 3D convection simulations. Some attempts have been made (e.g., Scott et al. 2021), but these procedures are not yet fully mature and still need to be improved for routine use in stellar evolution calculations.

It is not clearly established which of the criteria for convection should be used in stellar evolution computations. Moreover, when the Ledoux criterion is used, the efficiency of semiconvection is not known either. Since the pioneering works of Stothers & Chin (1973, 1975, 1976), very little progress has been made. Lawlor et al. (2015) simulated and compared the evolution and supernova (SN) light curves of Population III stars computed with either criterion and found Schwarzschild stars to be much larger and cooler than their Ledoux counterparts. Chun et al. (2018) calibrated the mixing-length parameter with red supergiants (RSGs) of different metallicities and predicted a decrease in the hydrogen content of stellar envelopes with increasing metallicity when using the Schwarzschild criterion, but not the Ledoux criterion. Schootemeijer et al. (2019) computed grids of Ledoux models at the metallicity of the Small Magellanic Cloud (SMC; Z = 0.002) and found that more efficient semiconvection increases the ratio of blue to red supergiants. Kaiser et al. (2020) computed non-rotating models of 15, 20, and 25 M_⊙ at solar metallicity to study the relative impact of various convective parameters. They found the strength of convective boundary mixing to be more effective than the choice of criterion for convection, except for determining the initial location of intermediate convective zones. Anders et al. (2022) performed a 3D hydrodynamical simulation of a convective zone adjacent to a semiconvective region. They reported that both criteria gave the same results when the evolutionary timescale exceeds the convective-overturn timescale (e.g., during the main sequence), but that differences can emerge when the evolution occurs on a rapid timescale (e.g., the Hertzsprung gap is crossed). We pursue our efforts and to try to determine which criterion reproduces the observed features of massive stars better (Saio et al. 2013a; Georgy et al. 2014b, 2021). To do this, we provide the first detailed comparison of Schwarzschild and Ledoux criterion models over a wide mass range at solar metallicity, including the effects of rotation.

In Sect. 2, we present the main physical ingredients of the stellar models. We discuss the impact of changing the convective criterion on the cores and convective structures of stars in Sect. 3. Section 4 presents the impact on various aspects of post-MS stellar evolution. We compare the properties of stars at the end of their evolution in Sect. 5. We present models for population synthesis in Sect. 6. We compare our results to previous works on the subject and to an observed population of evolved massive stars in Sect. 7. The discussion and conclusions are given in Sect. 8.

2. Ingredients of the models

The models computed with the Schwarzschild criterion were presented in Ekström et al. (2012). The convective regions are determined starting from the centre of the star and going outwards. In each shell of the model, the adiabatic and the radiative gradient are computed and compared, which determines whether the shell is to be radiative or convective. During central H- and He-burning, the size of the core is artificially increased by a length corresponding to 10% of the local pressure scale height at the edge of the formal core (overshooting). This is not done for subsequent burning phases or for convective shells. In the internal convective zones, the thermal gradient is imposed to be adiabatic, and the chemical species are instantaneously mixed (so that the chemical composition of the convective zones is always homogeneous). In the envelope, the classical mixing-length theory (Böhm-Vitense 1958) is used to obtain the thermal gradient, with a mixing-length parameter α = 1.6. For models more massive than 40 M_⊙, the mixing-length parameter is set to 1, and the mixing-length is computed according to the density scale height instead of the pressure scale height. Moreover, the turbulent pressure is accounted for by adding an acoustic flux term in the mixing-length formulation (see Maeder 1987).

The Ledoux models are new models that are computed with the same physical ingredients as the Schwarzschild models, except that the μ-gradient is computed and accounted for when a shell is to be determined as radiative or convective. All the other ingredients are kept the same. In particular, the treatment of overshooting and of mixing is the same, and the mass-loss prescriptions are the same as well.

Rotation is also treated in the same way in both sets of models. The full advection-diffusion equation for the transport of angular momentum is solved during the main sequence (see e.g., Eggenberger et al. 2008, for more details). The horizontal diffusion coefficient that intervenes in this formalism was proposed by Zahn (1992), and the diffusion coefficient associated with the shear mixing was proposed by Maeder (1997). All the rotating models shown in this paper have an initial rotation rate υ_ini/υ_crit = 0.4, where the critical velocity υ_crit is computed as in Georgy et al. (2011, see their definition of υ_crit,1).

Mass-loss is applied to all of our models according to different recipes, which are summarised here (see also Ekström et al. 2012). For non-Wolf-Rayet (WR) stars hotter than log(T_eff) = 4.0, the mass-loss rates by Vink et al. (2000, 2001) are applied. In this case, a correction factor accounting for the effects of rotation on the mass loss is also applied (see e.g., Georgy et al. 2011). For models up to 9 M_⊙ and when log(T_eff) < 3.8, the Reimers (1977) rates are used, with the parameter η = 0.6. For more massive models, the mass-loss rate is a linear fit of observational data from Sylvester et al. (1998) and van Loon et al. (1999; see also Crowther 2001). When the surface abundance of H drops below X_s < 0.3 and the effective temperature is above T_eff > 4.0, we switch to WR-kind mass-loss rates (Nugis & Lamers 2000; Gräfener & Hamann 2008), unless they are lower than the Vink et al. (2000) rates. Finally, we apply the de Jager et al. (1988) rates when the above recipes are not applicable.

For massive RSGss, some layers in the external layers, very close to the surface, reach luminosities far above the Eddington luminosity. To facilitate the computation in this case, we therefore increase the mass-loss rate as described above by a factor of 3 when the over-Eddington luminosity is overcome by a factor of 5 for models more massive than 20 M_⊙.

3. Convective criterion and internal structure of the models

Changing the convective criterion impacts the sizes of the convective zones and thus the chemical structure of a star. Before discussing the impacts on the observable outputs of the stellar models using different criteria, we discuss here the direct impact on various internal properties of our stellar models. In Table 1, we present for each model the masses of the He- and CO cores at the end of the central H- and He-burning phases, respectively, the maximum extent (in Lagrangian mass coordinates) of the intermediate convective zone (ICZ) before the beginning of core He-burning, and the deepest reach of the outer convective zone (OCZ) at different evolutionary stages.

We define the cores in the following way: the He core is the region inside the mass coordinate where at the first time the stellar layers are scanned from the surface towards the interior the mass fraction of helium Y > 0.9. The boundary of the CO core is where for a similar scanning, it crosses Y < 10⁻².

The masses of the He cores at the end of the MS phase obtained with the two different criteria show very small differences. For non-rotating models, the differences are always smaller than 3%¹. The same difference is obtained by comparing rotating models. These differences are blurred by other uncertainties pertaining to the models and are consistent with no effect due to the change in the convective criterion. This is consistent with the fact that in phases during which the convective core size decreases in mass, no differences are expected when the Schwarzschild or the Ledoux criterion are applied because there is no mean molecular weight gradient just above the core.

The only exception is the case of the rotating 120 M_⊙ model. The core mass obtained with the Schwarzschild criterion (33 M_⊙) is much lower than the core mass obtained with the Ledoux model (55.6 M_⊙). The structure evolution of the two rotating 120 M_⊙ models shows intermediate convective shells that appear just above the core in the Schwarzschild models and quickly merge with the core. This grants more fuel to the core, prolongs the main sequence by ∼100 kyr, and causes greater mass-loss of the star.

Rotation increases the mass of the He core at the end of the MS phase. This increase is very similar for the Schwarzschild and Ledoux criteria. It is 5–6% for the 9 M_⊙, reaches a maximum for a mass about 40 M_⊙ where it amounts to 26–27%, and then decreases because mass loss by stellar winds becomes important and thus blurs the effects of rotation. The rise in the He-core mass with initial mass is linked to the higher efficiency of rotational mixing when the initial mass increases (see the discussion in Sect. 3.1 of Maeder & Meynet 2000).

Figure 1 shows how the mass of the convective core evolves during the core He-burning phase in the non-rotating stellar models for initial masses between 9 and 15 M_⊙. The rotating cases (not shown here) exhibit very similar qualitative behaviours as those shown in Fig. 1: Their core masses increase during He-burning, and there are frequent breathing pulses² that inject helium from the H-burning shell into the core. Interestingly, the convective cores are more massive in the 12 and 15 M_⊙ models when the Ledoux criterion is applied. At first sight, this appears to be counter-intuitive. A stricter criterion for instability (Ledoux) should produce smaller convective cores than a less strict criterion (Schwarzschild). However, changing the criterion also affects the sizes of other convective zones in the models, and the changes in these other convective zones in turn affect the size of the convective core (see below).

Fig. 1. Relative mass of the convective core as a function of central He (from 1 to 0) during the He-burning phase for the non-rotating Ledoux (red) and Schwarzschild (blue) models between 7 and 15 M⊙.

The CO-core masses (see Table 1) present modest differences for masses below and including 15 M_⊙. The largest differences in the higher-mass star range are due to differences in the mass lost by stellar winds.

As expected, when the Schwarzschild criterion is used, larger intermediate convective zones associated with the H-burning shell appear. This has important consequences for the evolution of the He core, the occurrence of a blue loop (Maeder & Meynet 2001), and for the surface abundances of stars, which evolve back to the blue after having lost a large amount of mass during the red supergiant stage (Saio et al. 2013b; Georgy et al. 2014b, 2021). The ICZ tends to increase the contribution of the H-burning shell to the total luminosity (as it transports energy more efficiently), and as a result, the star reacts to this by decreasing the energy generated in the He-burning core. The Schwarzschild core has a lower central temperature than the Ledoux core and is therefore less dense. The Schwarzschild core has the same extent in radius as the Ledoux core, but its mass is lower, as shown in the bottom panels of Fig. 1.

A convective envelope does not appear in the high-mass range (above 60 M_⊙) as a result of mass loss, which prevents the stars from evolving in the red part of the Hertzsprung–Russell diagram (HRD). The outer convective zones for the 40 M_⊙ models are very small. The 32 M_⊙ initial mass is a transition case, without a convective envelope for the rotating models and with a significant envelope, reaching down to 25.5 M_⊙ for Schwarzschild and 16 M_⊙ for Ledoux, in the non-rotating case. For masses below and including 25 M_⊙, all the models have significant convective envelopes during the core He-burning phase. Models between 12 and 25 M_⊙ computed with the Ledoux criterion have deeper outer convective zones, but this is not the case for the 7 and 9 M_⊙ (see Sect. 4.2).

It might appear surprising that adopting a criterion that hinders convection can increase the extent of the outer convective zone. This extent also strongly depends on the position of the star in the HRD during the core He-burning phase, however: A hotter position during that phase decreases the extent of the outer convective zone. The position in the HRD at which core He-burning occurs is very sensitive to the variation in the abundances above the core as well as to the presence or absence of an intermediate convective zone. The intermediate convective zone typically tends to make the star more compact (i.e., blue), which then reduces the size of the convective envelope. In this manner, a criterion favouring convection (Schwarzschild) can produce a larger ICZ and consequently a less extended OCZ.

4. Evolutionary tracks and lifetimes

We show the evolutionary tracks for all models in Fig. 2. Main-sequence tracks are as expected to be similar in the Schwarzschild and Ledoux models. The main difference occurs in the rotating 120 M_⊙ models, as mentioned in Sect. 3.

Fig. 2. Evolutionary tracks of our stellar models in the HRD, colour-coded by the central helium mass fraction Yc. The grey shaded area represents the Cepheid instability strip (Tammann et al. 2003).

After the MS phase, the largest differences between the tracks (whether rotating or not) computed with the Schwarzschild and Ledoux criteria occur for the mass range between 7 and 40 M_⊙. In the upper mass range (between 60 and 120 M_⊙), stellar winds are the main factor governing the evolution of the stars. In the next paragraphs, we discuss in more detail how changing the convection criterion affects the first crossing of the Hertzsprung gap (Sect. 4.1), the blue loops (Sect. 4.2), the properties of the Cepheids (Sect. 4.3), the properties of stars ending their evolution as red supergiants (Sect. 4.4), and the properties of stars ending their evolution as yellow (YSG), blue supergiants (BSG), or as Wolf-Rayet stars (Sect. 4.5).

4.1. First crossing of the Hertzsprung gap

Figure 3 shows the evolution of the effective temperature as a function of the mass fraction of helium at the centre Y_c for all the models with initial masses between 7 and 40 M_⊙.

Fig. 3. Effective temperature towards the end of the main sequence (Yc > 0.8) and core He-burning for non-rotating (solid curves) and rotating (dashed curves) Ledoux (red) and Schwarzschild (blue) models. Upper row: masses between 7 and 15 M⊙. Lower row: masses between 20 and 40 M⊙. In each row, the vertical axes in all panels share the same limits, which are indicated in the leftmost panels.

Models with masses below and including 12 M_⊙ all begin their core He-burning phase when the star is a red supergiant, with log(T_eff [K]) ∼ 3.6. For the 15–25 M_⊙ mass range, the Ledoux criterion favours a quick crossing of the Hertzsprung gap, whether the models are computed with or without rotation. This is likely to be related to the fact that the use of the Ledoux criterion restrains the apparition and/or the extension of the intermediate convective zone associated with the H-burning shell, as shown in Fig. 4, which shows Kippenhahn diagrams for the post-MS phases of the four 15 M_⊙ models. This tends to produce steeper gradients of H and He in that zone and favours helium ignition in the core in the red supergiant stage (Farrell et al. 2022). The models of 15 and 20 M_⊙ with the Schwarzschild criterion (as well as the non-rotating 25 M_⊙ model) show a different behaviour. These models ignite helium in their cores when the star still has a high effective temperature (log T_eff [K] > 4.0), whether it rotates or not. This is linked to the larger intermediate convective zone in these models when the Schwarzschild criterion is used. A tentative explanation is that convection is more effective at transporting energy than radiation, and therefore, the presence of this ICZ allows the increased luminosity generated by the contraction of the core to reach the surface more easily. In the Ledoux models, however, energy is transported by radiation, and a larger fraction of it is deposited into the envelope itself, causing its rapid expansion and a drop in its luminosity. In both cases, the core contraction timescales are similar; the ICZ only affects the timescale of the envelope expansion. This means that at the onset of core He-burning, the Schwarzschild star remains in a blue part of the HRD.

Fig. 4. Kippenhahn diagrams showing the post-MS phases for the four models computed with an initial mass of 15 M⊙. The vertical axis shows the mass coordinate, and the numbered lines show isoradius levels. Left column: Schwarzschild models. Right column: Ledoux models. Upper row: non-rotating models. Lower row: rotating models. The intermediate convective zone above the H-burning shell is clear in the Schwarzschild models and is associated with the slower expansion of the envelope.

The cases of the more massive models (equal to or higher than 25 M_⊙) result from intricate interactions between convection, rotation, and mass loss. In general, helium ignition in the core of these models always occurs at an effective temperature below log T_eff [K] = 4.0. The only exception is the non-rotating 25 M_⊙ Schwarzschild model.

The evolution of the luminosity during the first crossing varies depending on the convective criterion. In general, the Ledoux models for masses between 15 and 20 M_⊙ experience a stronger decrease in luminosity before going up the Hayashi line than Schwarzschild models because these models are in a stronger radiative disequilibrium than the Schwarzschild models. Because they expand more rapidly, a larger amount of energy per unit time is tapped from the gravitational contraction of the core to expand the envelope. This tends to decrease the amount of energy that is radiated away. This decrease is to some extent a consequence of the time taken to cross the Hertzsprung gap. If it occurs on a sufficiently long timescale, the decrease in luminosity is modest. When it occurs on a rapid timescale, the luminosity drop is marked. To summarise, the Schwarzschild criterion favours the apparition of an ICZ above the H-burning shell after the main sequence, and for the reasons already described just above, this favours He ignition in the core when the star is still in the blue part of the HRD. The crossing of the Hertzsprung gap occurs on a longer timescale (nuclear, about 500 kyr) than with the Ledoux criterion (thermal, about 30 kyr). With the Schwarzschild criterion, the envelope has more time to radiate the excess energy produced by the contraction of the core away, and the Hertzsprung gap crossing occurs at roughly constant luminosity. The Ledoux criterion suppresses the ICZ, implying a shorter timescale for the Hertzsprung gap crossing. A large part of the energy radiated by the core contraction is absorbed by the envelope. As a result, the luminosity at the surface decreases significantly during the crossing.

4.2. Blue loops

Figure 3 shows that for initial masses above 20 M_⊙ the effective temperature can increase during core He-burning. These phases are not considered as blue loops for two reasons, however. First, these stars end their lifetimes at high effective temperatures as BSG or WR stars. Second, this evolution is mainly driven by mass loss and not by some internal evolution due to changes in the hydrostatic structure. In the rest of this section, we thus focus on the mass domain between 7 and 15 M_⊙.

For the lower-mass models, blue loops, characterised by an increase in effective temperature (blueward evolution) around Y_c ∼ 0.4 − 0.3, appear and are followed by a decrease (redward evolution) at the end of core He-burning. The physics of the blue loops was discussed in detail in previous papers (see e.g. Lauterborn et al. 1971; Maeder & Meynet 2001; Anderson et al. 2014). Many aspects of the stellar models influence the blue loops: for instance the rate of the ¹⁴N(p, γ)¹⁵O reaction (Halabi et al. 2012), of the ¹²C(α, γ)¹⁶O reaction (Brunish & Becker 1990), the undershooting below the outer convective zone (Lai & Li 2011), the helium and metal content (Bono et al. 2000), and rotation (Zhao et al. 2023). Recently, the occurrence of blue loops was used to constrain the magnetic moment of massive neutrinos (Mori et al. 2020) or the impact of axions (Friedland et al. 2013). Stellar models accounting for these effects undergo an additional loss of energy, and this tends to suppress the blue loops. This demonstrates that if the evacuation of energy from the core is facilitated, then this favours the suppression of the blue loops. Here, we focus on only one aspect, namely the impact of the convective criterion.

It has been shown that the blue loops may be linked to subtle differences in the abundance distributions near the H-burning shell, and that excess helium above the shell will suppress the blue loop (see e.g., Walmswell et al. 2015, who propose that the increase in mean molecular weight caused by having helium instead of hydrogen above the shell is the main culprit in suppressing the blue loop). Lauterborn et al. (1971) have proposed a criterion based on the gravitational potential of the core in order to decide whether a model produces a blue loop. Stars whose core has a gravitational potential higher than a given critical value remain along the Hayashi line, while those whose core has a gravitational potential lower than this critical value develop a blue loop. Interestingly, the quantity to be compared to the critical value is the core potential only if there is a steep chemical gradient just above the He core or near the H-burning-shell. Otherwise, the gravitational potential of the core has to be multiplied by a factor that is larger for a milder gradient, thus inducing the model to stay along the Hayashi line.

The non-rotating 7 and 9 M_⊙ models computed with the Ledoux criterion show no blue loops, while the 7 and 9 M_⊙ models with the Schwarzschild criterion have well-developed loops. Because these two masses behave very similarly in the features of their overall evolution, we focus on the 7 M_⊙ models, and the results qualitatively apply to the 9 M_⊙ models. Figure 5 shows the ⁴He profiles at different times before and during core He-burning for the two non-rotating 7 M_⊙ models. In the Schwarzschild model, the outer convective zone has equalised the abundances from the surface down to a mass near 1.6 M_⊙, reducing the extent of the zone with a chemical composition gradient (see the line corresponding to ⁴He_c where the gradient is between the mass coordinates 1.3 and 1.6 M_⊙). This occurs at the very beginning of the core He-burning phase. In the Ledoux model, the base of the outer convective zone always remains above the mass coordinate of 2.1 M_⊙ and the zone of the chemical composition gradient is larger (from 1.3 to 2.1 M_⊙), thus milder, than in the Schwarzschild model. This acts to prevent the formation of a blue loop. The Ledoux model also has an excess of helium above the H-burning shell (at 1.7–2.1 M_⊙), which suppresses the blue loop (see again Walmswell et al. 2015). The Ledoux models have this excess helium but the Schwarzschild models lack it because the mean molecular weight gradient above the He-burning shell stabilises this region against convection when the Ledoux criterion is used. As a result, the convective envelope does not reach down far enough to equalise the helium abundance to the surface value, leaving a helium excess above the shell.

Fig. 5. 4He profiles at different times before and during core He-burning for the two non-rotating 7 M⊙ models.

Models of 7 and 9 M_⊙ with rotation present less marked differences when the convective criterion is changed from the Schwarzschild to the Ledoux criterion. Rotational mixing blurs the differences that arise from different choices of the criterion for convection, and chemical composition gradients as well as helium abundances above the H-burning shell are similar for those models.

For the non-rotating 12 M_⊙, we do not observe any blue loop or any difference between the Schwarzschild and the Ledoux models. The only 12 M_⊙ model showing a blue loop is the rotating Ledoux model. This model would thus predict the most luminous Cepheids of those discussed in this paper (excluding the 15 M_⊙ models that cross the Cepheid instability strip very briefly).

4.3. Cepheid properties

Here, we discuss results concerning stars that cross the Cepheid instability strip (shown as the grey shaded region in Fig. 2 and in the small panels of Figs. 6 and 8). The initial masses of these stars are between 7 and 15 M_⊙.

Fig. 6. Durations of the different Cepheid phases for masses at which at least one of the four (Ledoux, Schwarzschild, rotating, non-rotating) models has a Cepheid phase. We distinguish between the first crossing (redward crossing of the Hertzsprung gap after the MS), the second crossing (blueward during the blue loop), and the third crossing (redward during the blue loop). These three moments are indicated in the small top right inset by coloured lines and dots with corresponding colours. Left panel: non-rotating model. Right panel: rotating model.

Figure 6 shows the durations of the different Cepheid phases. As is already well known, the first crossing is very short in all cases: It is much shorter than the blue loop when one occurs. As a result, most observed Cepheids should be undergoing a blue loop rather than crossing the Hertzsprung gap for the first time. For non-rotating models, the Schwarzschild criterion in general favours a longer duration of the Cepheid phase. The 7 and 9 M_⊙ Schwarzschild models undergo blue loops and the Ledoux models do not. Moreover, the 15 M_⊙ Schwarzschild models start core He-burning at an effective temperature of log(T_eff [K]) ∼ 4 and thus subsequently evolve on a helium nuclear timescale, whereas the Ledoux stars begin core He-burning as red supergiants and thus cross on a much shorter Kelvin-Helmholtz timescale.

For rotating models, the situation is less different between the Schwarzschild and Ledoux models, although the case of the 9 M_⊙ shows very significant differences: the Schwarzschild model spends 110 kyr as a Cepheid, while the Ledoux model spends only 20 kyr. We find no obvious correlation between the overall duration of the Cepheid phase and the stellar mass.

Figure 7 shows the surface abundances of ¹⁴N before and during the blue loop. In non-rotating models, there is no nitrogen surface enrichment during the MS phase, and as a result, no model shows nitrogen enrichment during the first crossing. Non-rotating models show an increase in surface nitrogen when they go through the RSG stage, where some nitrogen is dredged up to the surface through convection. This dredge-up occurs before the blue loop and is thus visible during the blue loops for the 7 and 9 M_⊙ models. For rotating models, the impact of rotational mixing during the MS is apparent in the surface ¹⁴N during the first crossing as it is larger than the initial value for all models. Interestingly, in the nitrogen surface abundance of the 9 M_⊙ Schwarzschild models, the blue-loop value obtained for the non-rotating model is equivalent to that given by the rotating model before the blue loop. This shows that the surface nitrogen abundance by itself is a poor indicator of the stellar phase (before or after the blue loop), unless its surface velocity is measured (see more on this point below). The nitrogen surface abundance increases with the initial mass, and it is slightly higher for Schwarzschild than for Ledoux models. It is higher in rotating models than in non-rotating models during the blue loop.

Fig. 7. Surface abundances of 14N before (cyan bar) and during (yellow bar) the blue loop for stars between 7 and 15 M⊙. The dotted line indicates the initial surface abundance. Left panel: non-rotating model. Right panel: rotating model.

Figure 8 shows the surface velocities of rotating stars at the three moments described above, as well as at the hottest point of the blue loop. In the 9 M_⊙ Schwarzschild model, the surface rotation speed along the blue loop strongly increases during blueward evolution (from 10 km s⁻¹ up to nearly 100 km s⁻¹ or nearly 70% of the critical velocity). Interestingly, it decreases very fast during the redward evolution after reaching the hottest point of the blue loop. The strong surface acceleration arises because during the RSG phase, an extended convective envelope is present. This convective zone rotates as a solid body, which means that the angular momentum is highest at the outer border of the convective zone. When the star contracts, the angular momentum accumulated by convection in the outer layers produces a rapid acceleration. The star loses angular momentum by winds during the blue loop, and it has a larger radius after than before (the effective temperature is similar, but the luminosity has increased). As a result, when it evolves back to the red, its situation is not the same as before the blue loop: Its surface velocity is lower. The Ledoux model shows a qualitatively similar behaviour, but the different evolution of the outer convective zone makes the evolution of the surface rotation slightly different.

Fig. 8. Surface velocities of rotating stars during the first crossing of the Hertzsprung gap and at three different times of the blue loop for stars between 7 and 15 M⊙. These four moments are indicated in the small top right inset by coloured lines and dots with corresponding colours. Left panel: equatorial velocity υsurf. Right panel: ratio of the equatorial to the critical velocity υsurf/υcrit.

4.4. Red supergiants

In this and the next section, we examine the impact of the convective criterion on different subclasses of post-MS stars. These subclasses are defined in Table 2.

Here, we discuss whether the convective criterion affects the maximum luminosity of red supergiants. Figure 2 shows that the maximum luminosity reached for red supergiants is not much affected by the convective criterion, but it is affected by rotation. Rotating models produce an upper limit for the RSG luminosity of about log(L [L_⊙]) = 5.4, while the non-rotating models extend this upper limit to values of about log(L [L_⊙]) = 5.7.

When we now compare the upper luminosity of red supergiants that remain red until the end of their evolution and thus will be progenitors of a type II SN, then the upper limits are lower than those in the previous paragraph because mass loss by stellar winds induces a blueward evolution and thus produces yellow or blue progenitors in the higher-luminosity range (see Sect. 4.5). This upper limit in luminosity for red supergiant core-collapse progenitors is about log(L [L_⊙]) = 5.2 for the non-rotating Schwarzschild models, log(L [L_⊙]) = 4.95 for the non-rotating Ledoux models, log(L [L_⊙]) = 5.1 for rotating Schwarzschild models, and log(L [L_⊙]) = 5.05 for the rotating Ledoux models. We note that the models that spend a larger fraction of their core He-burning phase in the red supergiant phase have slightly lower upper limits. This is expected because their mass is reduced by strong RSG winds for a longer duration. This tends to lower the minimum mass above which stars later evolve blueward. For instance, in the non-rotating case (with the largest difference in maximum RSG luminosity), the most luminous Schwarzschild RSG has an initial mass of 20 M_⊙ while for the Ledoux RSG, the initial mass is 15 M_⊙ because the 20 M_⊙ model becomes blue after its mass-loss episode.

4.5. Yellow and blue supergiants, luminous blue variables, and Wolf-Rayet stars

The mass domain between 20 and 40 M_⊙ is a transition domain between the stars that evolve into the red supergiant stage and remain a red supergiant until the end of their lifetimes and the stars, for instance the 60 M_⊙ models, that never evolve into a RSG, but become a luminous blue variable before evolving into a WR phase. In this transition mass domain, the star can end its evolution as a red, yellow, or blue supergiant, or become a Wolf-Rayet star (in the latter case after having been a red supergiant for a while). The tracks show a complex behaviour resulting from intricate effects involving convection, rotational mixing, and mass loss.

Below an initial mass M_RSG ∼ 15 M_⊙, stars end their lives as red supergiants (see Sect. 4.4). Stars initially between M_RSG and M_RSG-WR ∼ 20 − 25 M_⊙ end their lives as a yellow or blue supergiants. For initial masses between M_RSG-WR and M_WR ∼ 40 M_⊙, stars end their lives as Wolf-Rayet stars after a red supergiant stage. We note that this mass range of stars that experience both an RSG and a WR phase is expected to be narrow because at the moment, apart from Westerlund-1 (see e.g., Clark et al. 2005; Crowther et al. 2006; Negueruela et al. 2010; Beasor et al. 2021, we present isochrones of our models compared with observed data of massive evolved stars in Westerlund-1 in Sect. 7.2), there is no single-age stellar population in which both red supergiants and Wolf-Rayet stars are observed at the same time. Finally, stars with initial masses above M_WR enter the WR phase without previously becoming a RSG.

These limits appear to be more sensitive to rotation than to the convective criterion (at least for the initial rotation speed considered here). In general, rotation tends to decrease the limits M_RSG-WR and M_WR (see more below). Because M_WR is shifted to lower values than M_RSG-WR in rotating models, this implies that rotation suppresses or at least disfavours the production of WR stars with RSG progenitors (at least for a single-star evolution).

Mass loss by stellar winds (which increases when the luminosity increases and the effective temperature decreases) begins to become a dominant feature in this mass domain (it becomes an even more important feature for higher initial masses). As a numerical example, the 20 M_⊙ models lose significant amounts of mass during the post-MS phase (see Table 3), ending core He-burning with a total mass of about 7–9 M_⊙ (for comparison, the 15 M_⊙ models reach the end of He-burning at around 11–13 M_⊙). The masses of the 25 M_⊙ models at the end of He-burning are about 9–10 M_⊙ (this is higher than for the 20 M_⊙ models, but more mass has been lost by the 25 M_⊙ models). The 32 and 40 M_⊙ models experience even stronger mass loss (losing up to 22 and 28 M_⊙ for the 32 and 40 M_⊙ models, respectively). They all end core He-burning at high (log T_eff [K] > 4.4) effective temperatures, and some of them lose so much mass that they become Wolf-Rayet stars.

While all the stars including and above 20 M_⊙ lose large amounts of mass, the mass of the core remains an increasing function of the initial mass. In other words, the 20 M_⊙ models reach the end of He-burning with a lower total mass than the 15 M_⊙ models, but their CO cores are more massive (3.65–4.44 M_⊙ compared with 2.11–2.66 M_⊙). Furthermore, the models between 20 and 40 M_⊙ experience much stronger surface nitrogen enrichment than those below 20 M_⊙.

The mass domain above this transition mass range is dominated by mass loss. Not many differences are therefore visible between Schwarzschild and Ledoux models. The most noticeable difference occurs for the 120 M_⊙ models, as discussed in Sect. 3.

Figure 9 shows the durations of the different Wolf-Rayet phases (stacked in chronological order) for the masses where at least one of the four models reaches a WR phase. A first striking difference appears between the rotating and non-rotating models. Rotation produces longer WR phases. This is in line with previous studies (Meynet & Maeder 2003; Georgy et al. 2012). The time spent as WR stars is longer for rotating than for non-rotating models because rotation drives additional mass loss and mixing. More mixing implies that less mass loss is required to uncover hydrogen-poor layers. Rotation also favours the WNL phase, decreases WNE and WC durations, and prevents the WO phase. Interestingly, rotation increases the WNC-phase duration.

Fig. 9. Durations of the different Wolf-Rayet phases for masses at which at least one of the four (Ledoux, Schwarzschild, rotating, and non-rotating) models has a Wolf-Rayet phase. Left panel: non-rotating models. Right panel: rotating models.

Above and including 60 M_⊙ the Schwarzschild and Ledoux models spend similar times (within 6%) as Wolf-Rayet stars because mass loss is the dominating effect here. This agrees with the fact that the differences between Ledoux and Schwarzschild models are mostly due to the ICZ. They do not appear in the most massive stars due to the strong mass loss. As a result, the choice of criterion for convection intervenes as a second-order effect.

For the mass range between 25 and 40 M_⊙, the differences induced by changing from the Ledoux to Schwarzschild criterion are more marked. Without rotation, the Ledoux criterion produces Wolf-Rayet stars for lower initial masses than the Schwarzschild criterion. These Ledoux models spend a longer time of their core He-burning phase as an RSG than the Schwarzschild models. As a result, mass loss is stronger for the former, and they reach the WR phase earlier (as a result, the WR phase lasts longer).

5. Final properties of the models

5.1. Core masses

Figure 10 shows the CO-core masses and the mass ratios of the He and CO core as functions of the initial mass for the four sets of models at the end of their evolution. The masses of the cores are defined as in Table 1 (here, we discuss their values at the last computed stage, however). The ratio of the mass of the He core and that of the CO core decreases from 2 to roughly 1.3 between 7 and 32 M_⊙ and then remains constant around 1.2–1.3, indicating that the He-burning shell sits relatively farther out from the C-burning core in lower-mass models than in higher-mass models.

Fig. 10. Mass of the 4He core (left panel) and mass ratio of the He and CO cores (right panel) at the last computed stage as a function of the initial mass for the different models. Stars (blue: Schwarzschild, and red: Ledoux) represent non-rotating models, and crosses (green: Schwarzschild, and orange: Ledoux) show rotating models.

The differences among models of the same mass in the left panel of Fig. 10 show that rotating models have a 10% larger core on average (except for the rotating 32, 40, and 120 M_⊙ Schwarzschild models) than non-rotating models. There is no clear trend in the differences between Schwarzschild and Ledoux models, however. The impact of rotation is stronger than that of changing the convective criterion.

The strongest impact of changing from the Ledoux to Schwarzschild criterion appears for the 7 and 9 M_⊙ stellar model. The ratio of He- to CO-core mass can vary from 1 to 2, depending on the criterion. This mass domain covers the transitions between stars that would produce white dwarfs and those that would produce neutron stars at the end of their lifetimes. This transition has been studied for instance by Siess (2006, 2007). This final fate depends on whether ignition of carbon occurs in degenerate, mildly degenerate, or non-degenerate conditions. A star whose central region is for a large part in mildly degenerate conditions can still succeed in igniting (off-centre) carbon, but this requires a more massive CO core than in non-degenerate conditions. It is therefore expected that a larger part of the central regions for models with the lowest ratio of the He- to the CO-core mass lies in the degenerate domain. Only models that ignite carbon in non-degenerate or very mildly degenerate conditions (i.e., the rotating 9 M_⊙ models) show high ratios of the He to CO core.

5.2. Properties of the compact object progenitors

We list integrated chemical abundances (in mass fractions) of ¹H, ⁴He, ¹²C, and ¹⁶O in the models at the end of core C-burning (core He-burning for stars that do not reach the end of C-burning) in Table 4. This moment is close enough to the end of evolution that the abundances of the isotopes shown in the ejecta will not change dramatically (see e.g., the appendix table in Limongi & Chieffi 2018). Furthermore, explosive nucleosynthesis during the supernova is not expected to significantly affect the yields of these four elements. We computed the remnant types and masses using the models from Patton & Sukhbold (2020), whose predictions concerning the remnant type (neutron star or black hole) were based on the ‘Ertl criterion’ (Ertl et al. 2016, 2020). For each star, we took the mass of the CO core M_CO (defined, like in Table 1, as the region within which the mass fraction of helium Y < 10⁻²) and the mass fraction of ¹²C in the core X_¹²C at the end of He-burning. When M_CO < 1.4 M_⊙, then we considered the remnant to be a white dwarf with a mass of M_WD = M_CO. When 1.4 M_⊙ < M_CO < 2.5 M_⊙, the remnant is a neutron star with a mass of M_NS = 1.4 M_⊙. When 2.5 M_⊙ < M_CO < 10 M_⊙, we consulted the top panel of Fig. 3 of Patton & Sukhbold (2020) to determine whether the star will explode or implode (without ejecting anything). If it explodes, then we estimated the baryonic³ mass of the resulting neutron star to be M_NS = M₄, where M₄ is defined in Patton & Sukhbold (2020, we use their Table 1 which provides the values of M₄ for a range of (M_CO, X_¹²C) at the end of core He-burning). If it implodes, then the mass of the resulting black hole is M_BH = M_fin. Finally, if 10 M_⊙ < M_CO < 30 M_⊙, the outcome is a black hole, with M_BH = M_fin. This equality between the black hole mass and the final mass of the star may be an overestimation, but the underlying assumption is that the entire star collapses into the black hole and that no matter is ejected. Any value between M_CO and M_fin would be a reasonable estimate for M_BH. The case M_CO > 30 M_⊙ could yield a pulsational or regular pair-instability supernova, but this applies to none of our stars. We then integrated the abundances of ¹H, ⁴He, ¹²C, and ¹⁶O above the remnant (between its mass coordinate and the surface). If the remnant is a black hole, we list the integrated chemical abundances above the mass coordinate of the CO core. We grouped all other elements into the quantity 𝒵 (usually, the metallicity Z includes ¹²C and ¹⁶O, and therefore, we used the calligraphic 𝒵 = 1 − X − Y − X(¹²C)−X(¹⁶O)). We also provide the mass, effective temperature, and luminosity of the final models.

Using the criteria defined above on M_CO, as well as the grid of CO cores with varying ¹²C core mass fractions computed by Patton & Sukhbold (2020), we can predict the type of compact object that would remain after the end of stellar evolution. All of the 7 and 9 M_⊙ models would yield white dwarfs (WD in Table 4). Of the 12, 15, and 20 M_⊙ models, the rotating 15 M_⊙ Schwarzschild model would become a black hole, while all the others would become neutron stars. In the 25 M_⊙ stars, the Ledoux models become black holes and the Schwarzschild models become neutron stars. All stars 32 M_⊙ and above become black holes.

The composition of the envelopes shows that, as expected for stars that reach the WR phase, the models above 20 M_⊙ do not have any ¹H at the end of their evolution. More generally, the integrated abundances of both ¹H and ⁴He decrease with increasing initial mass.

Conversely, the quantity of metals (including ¹²C and ¹⁶O) in the envelope is an increasing function of the initial mass. We find no striking consistent difference between the Schwarzschild and Ledoux models. The differences mostly concern the integrated ¹²C and ¹⁶O abundances: in some cases, one model produces more ¹²C and less ¹⁶O than its counterpart, but there is no clear-cut effect.

The case of the 25 M_⊙ models is interesting because the two Ledoux models are predicted to become black holes and the Schwarzschild models to become neutron stars. There is a striking difference in the compositions of their envelopes: the former are much richer in ⁴He (Y_env of 0.58 and 0.74 against 0.20 and 0.33 for the Schwarzschild models). The latter, conversely, are much richer in ¹²C and ¹⁶O (with the most marked differences being in the ¹⁶O abundance: X(¹⁶O)_env of 0.525 and 0.417 against 0.188 and 0.105 for the Ledoux models).

In the 32 and 40 M_⊙ models, the rotating 32 M_⊙ and non-rotating 40 M_⊙ Schwarzschild models present very different ¹²C and ¹⁶O abundances from the other models of the same mass. The former is much richer in ¹²C and ¹⁶O (and thus poorer in ⁴He) than the other 32 M_⊙ models; the latter is much poorer in ¹²C and ¹⁶O (and thus richer in ⁴He) than the other 40 M_⊙ models.

From the ¹H and ⁴He content of the models at the last computed stage, we can infer the type of core-collapse supernova that is expected to occur for all the exploding models. All the models of 12 and 15 M_⊙ contain significant amounts of ¹H (about 50% of the envelope mass), and therefore, we expect type II SNe from them (except for the rotating Schwarzschild 15 M_⊙ model, which we expect to directly collapse into a black hole). Of the 20 M_⊙ stars, only the non-rotating Schwarzschild model is expected to lead to a type II SN because it still contains a non-negligible amount of ¹H. We expect the other 20 M_⊙ models to die in type Ib SNe because they have no ¹H left, but still contain ⁴He. The same holds for the 25 M_⊙ Schwarzschild models, for which we expect an explosion.

Finally, we theoretically predict the remnant masses. The most massive neutron star comes from the non-rotating 20 M_⊙ Schwarzschild model, with a baryonic mass M_NS = 1.77 M_⊙. The least massive black hole results from the non-rotating 25 M_⊙ Ledoux model, with M_BH = 9.3M_⊙. Acknowledging that taking M_BH = M_fin may be an overestimation of the black hole mass, we can estimate the absolute minimum predicted black hole mass from that of the CO cores (taking M_BH = M_CO). In this case, the 15 M_⊙ rotating Schwarzschild model yields the lowest-mass black hole, with M_BH = 2.68M_⊙.

6. Stellar population synthesis

6.1. The population synthesis models

The population synthesis results we present in this work were computed using the SYCLIST code (Georgy et al. 2014a), using the stellar models from Ekström et al. (2012, computed with the Schwarzschild criterion), and from Saio et al. (2013a) and Georgy et al. (2014b, 20 and 25 M_⊙ models computed with the Ledoux criterion) as input. We added all the other Ledoux models that we computed for the current paper in order to cover masses from 7 to 120 M_⊙. SYCLIST generates a stellar population with initial masses sampled according to the Salpeter IMF, and it saves the state of the population at each requested time step. To mimic continuous star formation rates, we summed the populations at each time step, and we then counted stars and classified them into subtypes according to some of their properties. The subtypes are defined in Table 2.

We generated two types of clusters. One type consisted of a burst of star formation at t = 0, and the other type had continuous star formation for 60 Myr. The second type allowed us to simulate a stationary state and derive equilibrium relative abundances of stellar subtypes. We also generated isochrones of our populations to compare them to the observed populations of evolved massive stars in Westerlund-1.

6.2. Case of a burst of star formation

Figure 11 shows the number of blue and red supergiants (top panels) and Wolf-Rayet stars (bottom panels) alive at each time in the case of a burst of star formation, normalised by the initial number of stars. For non-rotating stars, the Ledoux clusters produce far fewer blue and many more red supergiants than the Schwarzschild clusters. This is attributable to the different crossing durations of the Hertzsprung gap above 15 M_⊙ (the Schwarzschild models start core He-burning as BSGs, the Ledoux models start as RSGs), and to a larger extent (because of the IMF slope) to the blue loop behaviour, which is present in Schwarzschild stars between 7 and 12 M_⊙ but is absent in Ledoux models.

Fig. 11. Number of stars of certain types alive at each time, normalised by the total initial number of stars generated in a synthetic cluster for the case of a burst of star formation at t = 0. The solid lines represent stars computed with the Schwarzschild criterion, and the dashed lines show the same for the Ledoux criterion. Top panels: blue and red supergiants. The vertical axis is logarithmic, and the time axis extends to 60 Myr. Bottom panels: Wolf-Rayet stars. The vertical axis is linear and in units of 10−3, and the time axis extends to 12 Myr (there are no WR stars later on). Left panels: non-rotating models. Right panels: rotating models.

This is not so much the case for rotating stars because, as mentioned previously, rotation tends to mitigate the differences between the two sets of models, and the blue loop behaviour is much more similar in rotating models. The difference in these models occurs just before 20 Myr, which is when the Ledoux 12 M_⊙ model undergoes its blue loop while the Schwarzschild model does not.

Finally, we note that the non-rotating Wolf-Rayet populations (especially the WNL stars) are quite different, with Ledoux clusters clearly hosting more WR stars than Schwarzschild clusters. This is expcected because the WR stage is reached for non-rotating Ledoux stars of lower masses (down to 25 M_⊙), and it lasts longer for these lower masses (see the left panel of Fig. 9) than in the Schwarzschild models. WNE, WNC, and WO type stars are very rare in both types of clusters, and the differences in WC number densities are not significant.

6.3. Stationary regime populations

We now discuss populations in the stationary regime. Table 5 shows the number ratios of different classes of stars in this equilibrium state, where we simulated populations with a constant continuous star formation episode lasting 60 Myr. The ratios of post-MS to MS stars are very similar in Schwarzschild and Ledoux models: we find differences of 5% (non-rotating) and 9% (rotating), and the Schwarzschild post-MS to MS ratio is higher than the Ledoux ratio because Ledoux models have a slightly longer main sequence and shorter He-burning phases. The ratios of BSGs to RSGs are higher for the Schwarzschild populations regardless of the rotation scheme, although the difference is more marked for populations of non-rotating stars (the BSG to RSG ratio is almost five times higher for Schwarzschild than for Ledoux populations); this echoes what we mentioned in Sect. 6.2, and we attribute this difference to the different blue loop behaviour and to the durations of the Hertzsprung gap crossing.

The WR to RSG ratio is quite interesting in the non-rotating case. Ledoux populations produce more RSGs, but also more WR stars than Schwarzschild populations. The overproduction is higher for the WR stars, as WR/RSG_Ledoux > WR/RSG_{Schwarzschild}. Finally, while the populations produce similar numbers of WC stars (see the lower panels of Fig. 11), the ratio of WC to WN stars is three times higher for the non-rotating Schwarzschild population because the Ledoux population contains more WN(L) stars.

Figures 12 (non-rotating) and 13 (rotating) show density maps in the HRD of populations generated with continuous star formation. The densities are higher in the lower parts of the diagram because more lower-mass stars are generated due to the Salpeter IMF we used. The comparison of populations both with and without rotation shows that the MS regions are the same for the Schwarzschild and Ledoux models.

Fig. 12. Density maps in the HRD of the Schwarzschild and Ledoux non-rotating populations in the case of continuous star formation (stationary regime). We add the tracks of the 7, 15, 25, and 60 M⊙ models. Left: Schwarzschild model. Right: Ledoux model.

We note for the non-rotating populations that they contain a region of higher density that corresponds to the blue loops of the Schwarzschild models, but this region is absent from the Ledoux models. We also see differences in the supergiant stars above log(L [L_⊙]) ∼ 4.5 because the 20–40 M_⊙ models experience a blue and/or redward evolution in diverse parts of the HRD.

For the rotating populations, the blue loop stars are present in both sets of models, but their high-density regions extend to lower surface temperatures for the Ledoux models (log T_eff [K] ∼ 3.85) than for the Schwarzschild models (log T_eff [K] ∼ 3.95). This would not be noticeable from the HRD tracks alone (the loops even extend farther bluewards for the rotating Ledoux stars), which means that the Ledoux models spend less time at the blue edge of the loop than the Schwarzschild models do. Figure 13 shows that the higher-density region of the blue loops corresponds to the hottest point along the loop of 7 M_⊙ Schwarzschild stars, but the equivalent point of the Ledoux-star loop extends to a higher temperature. The high-luminosity part of the blue-loop region reaches higher luminosities for the Ledoux models (log(L [L_⊙]) ∼ 4.6) than for the Schwarzschild models (log(L [L_⊙]) ∼ 4.3) because the threshold for blue loops is between 12 and 15 M_⊙ for the former and between 9 and 12 M_⊙ for the latter. We also see differences in the higher luminosity supergiants, which we explain in a similar way to those of the non-rotating models.

Fig. 13. Density maps in the HRD of the Schwarzschild and Ledoux rotating populations in the case of continuous star formation (stationary regime). We add the tracks of the 7, 15, 25, and 60 M⊙ models. Left: Schwarzschild model. Right: Ledoux model.

7. Comparison with previous works and observations

7.1. Comparison with previous theoretical works

While this paper presents the first grid that spans such a wide mass range in the comparison of Schwarzschild and Ledoux models at solar metallicity, with and without rotation, it is not the first to study the impact of various convective parameters on stellar evolution. In this section, we compare our conclusions with a few such previous theoretical works.

Schootemeijer et al. (2019) used MESA (Paxton et al. 2011, 2013, 2015, 2018) to compute models between 9 and 100 M_⊙, with and without rotation, at SMC metallicity (Z = 0.002). They used the Ledoux criterion and varied parameters related to the strength of overshooting and efficiency of semiconvective mixing. The models relevant to our comparison are those with overshooting parameter α_ov = 0.11, and the most extreme values of the semiconvection parameter α_sc = 0.01 and 300 (or 100 because they discussed more results related to that value). They found that very efficient semiconvection (large α_sc, Schwarzschild in the current paper) leads to more time being spent by stars as blue supergiants. For instance, these models start to burn helium in the core at a higher effective temperature than those with inefficient semiconvection. Figure 2 (top panel for α_ov = 0.11) of Schootemeijer et al. (2019) shows that blue loops occur in the lower-mass range for efficient semiconvection, but not when it is inefficient. All of these observations lead to a higher ratio of blue to red supergiants when semiconvective mixing is efficient. Qualitatively, this agrees very well with what we find. The specific values for the mass range of the blue loops, the effective temperatures at the beginning of core He-burning, and the ratio of blue to red supergiants differ from what we predict, but that is to be expected because their models were computed at SMC metallicity (which is almost an order of magnitude lower than ours).

Kaiser et al. (2020) also used MESA to compute non-rotating models of 15, 20, and 25 M_⊙ at solar metallicity. They computed models with both the Schwarzschild and Ledoux criteria and used a free parameter f_CBM, where CBM stands for convective boundary mixing, to vary the amount of additional mixing at the convective boundary (larger f_CBM corresponds to more mixing). This parameter intervenes in a diffusive model rather than the penetrative overshoot that we used in this study. They found that the main-sequence evolution is not affected by the choice of criterion for convective stability as long as CBM is included (in our case, as long as overshooting is included). They also discussed the intermediate convective zones and found that ICZs that are larger and remain for longer mean that stars continue to be blue until the ICZ recedes. Their stars computed with the Schwarzschild criterion have stronger (meaning larger and longer-living) ICZs than those computed with the Ledoux criterion. As a result, their Schwarzschild stars start core He-burning as blue supergiants and the Ledoux stars start as red supergiants. Their 15 M_⊙ Ledoux models also cross the HRD much faster after the main sequence, and this is correlated with a drop in surface luminosity. Overall however, their results are more affected by the value of f_CBM than by the choice of Schwarzschild and Ledoux criterion.

Anders et al. (2022) performed a 3D hydrodynamical simulation of a convective zone adjacent to a semiconvective (Ledoux-stable but Schwarzschild-unstable) region. They found that overshooting mixes the semiconvective region, increasing the size of the convective zone by a process they called entrainment. After a few thousand convective-overturn times, the Ledoux and Schwarzschild criteria predict the same convective boundary. The consequence is that as the title of Anders et al. (2022) states, ‘Schwarzschild and Ledoux are equivalent on evolutionary timescales’. This is the case during the main sequence (when the evolutionary timescale is much longer than the convective-overturn timescale). However, when the two timescales are of similar order, Ledoux predicts the instantaneous location of the boundary, whereas Schwarzschild provides its location in a stationary state. This agrees with what we find: that the largest differences are introduced during the post-MS expansion, which occurs on a rapid timescale.

7.2. Comparison with observations

We briefly present here a comparison of our four sets of models with the evolved massive stars in the Westerlund-1 (Wd1) cluster. Wd1, at a distance of ∼4 kpc (Beasor et al. 2021), is the best-studied young cluster in the galactic disk, and it contains both RSG, YSG, and WR stars. Its metallicity is unknown because the gas associated with its formation has been dispersed. Based on its age, however, we may expect that its metallicity could be higher than that in our study (which is why Yusof et al. 2022 used it as a comparison for their super-solar metallicity (Z = 0.020) models). We nevertheless chose to also compare our models to Wd1 because it is the best-studied cluster containing both cool supergiants and Wolf-Rayet stars. We again used SYCLIST (Georgy et al. 2014a), this time, to generate isochrones for our four grids of models. We used the data from Yusof et al. (2022) for the observed evolved stars in Wd1.

Figure 14 shows isochrones for our four sets of grids and the observed data of evolved massive stars in Wd1. While no isochrone perfectly fits the observed distribution of stars, we find that the non-rotating Ledoux isochrone at an age of 6.3 Myr best incorporates the presence of RSG and WR stars, as well as the spread in luminosity of the RSGs. The 10 Myr Ledoux isochrones (both rotating and non-rotating) fit the observed YSGs. We recall the possibility that a cluster formation event lasted a few million years, typically between 10 and 6 Myr ago, and that real stellar populations have varied rotation rates. Moreover, multiple star evolution and cluster dynamics will impact the observable properties of stars and induce a scatter in their luminosities, and this can explain some of the observed data.

Fig. 14. HRD showing isochrones at log(Age [yr]) = 6.7,  6.8,  6.9,  and 7.0 (corresponding to ages of 5, 6.3, 7.9, and 10 Myr) for our four sets of grids, compared with observations of evolved massive stars in the Westerlund-1 cluster. We use the data from Yusof et al. (2022) for the observed stars. Left column: Schwarzschild models. Right column: Ledoux models. Upper row: non-rotating models. Lower row: rotating models.

8. Conclusions and discussion

In the above sections, we discussed the possible effect of changing the convective criterion on stellar model outputs with and without rotation. In many aspects, the effects appear to be relatively minor compared to those of rotation or mass-loss. We summarise the main results from this work below.

• The main-sequence phase is not much affected by the choice of the convection criterion.

• The Schwarzschild criterion does not necessarily lead to more extended convective zones. It does facilitate the formation of an extended intermediate convective zone at the end of the main sequence (especially for stars with initial masses between 15 and 32 M_⊙), however, which gives rise to significant differences in the post-MS evolution of these stars.

• One such difference occurs during the first crossing of the HRD, where stars computed with the Schwarzschild criterion can start core He-burning as BSG, whereas the Ledoux stars start as RSG. This has important implications in the subsequent mass loss of the models. This in turn affects the ratio of blue to red supergiants as well as the way in which case B mass transfer occurs in close binary systems. Because the Schwarzschild criterion causes longer crossing times, it predicts more nitrogen-poor Cepheids (those undergoing their first crossing of the HRD) than the Ledoux criterion. Observations of nitrogen-poor Cepheids would definitely tip the scale in favour of the Schwarzschild criterion (or efficient semiconvection).

• The occurrence of blue loops is affected by the choice of a Schwarzschild and Ledoux criterion: non-rotating 7–9 M_⊙ Ledoux stars do not have them, but Schwarzschild stars do. This is probably due to different chemical profiles and helium abundances near and above the H-burning shell, and it also affects the ratio of blue to red supergiants.

• The duration of and the surface velocities (for rotating stars) during the Cepheid phase are influenced by the choice of the convection criterion. The Schwarzschild criterion predicts longer-lasting and faster-rotating Cepheids than the Ledoux criterion. The Ledoux criterion predicts the most massive Cepheids.

• Red supergiants computed with the Ledoux criterion reach a lower maximum luminosity than those computed with the Schwarzschild criterion.

• The impact of the criterion for convection on stars with initial masses between 20 and 40 M_⊙ is complex because of the interplay with rotation and mass loss. For instance, the Ledoux criterion produces non-rotating WR stars for lower initial masses than the Schwarzschild criterion, but the inverse occurs for rotating stars.

• The impact of changing the convective criterion on stars whose initial masses are above 60 M_⊙ is modest because their evolution is dominated by mass loss.

• The main differences relative to the endpoint of stellar evolution occur for stars with initial masses of 25–40 M_⊙. For the 25 M_⊙ the Ledoux criterion predicts the final compact objects to be black holes, while the Schwarzschild criterion predicts neutron stars. If these stars still undergo a supernova, the Ledoux criterion predicts much stronger helium signatures than the Schwarzschild criterion because the amount of helium in their envelope is larger.

• Our population synthesis models imply different number ratios of evolved stars depending on the convection criterion, especially for non-rotating stars: The Schwarzschild criterion predicts much higher ratios of BSGs to RSGs and WCs to WNs than the Ledoux criterion. Conversely, the Ledoux criterion predicts a higher ratio of WRs to RSGs.

We recall, however, that the comparison presented in this work was made for a given choice of mass-loss rates, step overshoot, physics of rotation, and metallicity, and only for single-star evolution. Without performing detailed evolutionary calculations, it is not possible to know how the picture we presented in this work would change when one, a few, or all of the physical ingredients would change.

Acknowledgments

The authors would like to thank the anonymous referee for their very insightful comments which greatly improved the quality of this manuscript. The authors have received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 833925, project STAREX).

References

All Tables

All Figures

	Fig. 1. Relative mass of the convective core as a function of central He (from 1 to 0) during the He-burning phase for the non-rotating Ledoux (red) and Schwarzschild (blue) models between 7 and 15 M⊙.
In the text

	Fig. 2. Evolutionary tracks of our stellar models in the HRD, colour-coded by the central helium mass fraction Yc. The grey shaded area represents the Cepheid instability strip (Tammann et al. 2003).
In the text

	Fig. 3. Effective temperature towards the end of the main sequence (Yc > 0.8) and core He-burning for non-rotating (solid curves) and rotating (dashed curves) Ledoux (red) and Schwarzschild (blue) models. Upper row: masses between 7 and 15 M⊙. Lower row: masses between 20 and 40 M⊙. In each row, the vertical axes in all panels share the same limits, which are indicated in the leftmost panels.
In the text

Fig. 4. Kippenhahn diagrams showing the post-MS phases for the four models computed with an initial mass of 15 M⊙. The vertical axis shows the mass coordinate, and the numbered lines show isoradius levels. Left column: Schwarzschild models. Right column: Ledoux models. Upper row: non-rotating models. Lower row: rotating models. The intermediate convective zone above the H-burning shell is clear in the Schwarzschild models and is associated with the slower expansion of the envelope.

In the text

	Fig. 5. 4He profiles at different times before and during core He-burning for the two non-rotating 7 M⊙ models.
In the text

Fig. 6. Durations of the different Cepheid phases for masses at which at least one of the four (Ledoux, Schwarzschild, rotating, non-rotating) models has a Cepheid phase. We distinguish between the first crossing (redward crossing of the Hertzsprung gap after the MS), the second crossing (blueward during the blue loop), and the third crossing (redward during the blue loop). These three moments are indicated in the small top right inset by coloured lines and dots with corresponding colours. Left panel: non-rotating model. Right panel: rotating model.

In the text

	Fig. 7. Surface abundances of 14N before (cyan bar) and during (yellow bar) the blue loop for stars between 7 and 15 M⊙. The dotted line indicates the initial surface abundance. Left panel: non-rotating model. Right panel: rotating model.
In the text

Fig. 8. Surface velocities of rotating stars during the first crossing of the Hertzsprung gap and at three different times of the blue loop for stars between 7 and 15 M⊙. These four moments are indicated in the small top right inset by coloured lines and dots with corresponding colours. Left panel: equatorial velocity υsurf. Right panel: ratio of the equatorial to the critical velocity υsurf/υcrit.

In the text

	Fig. 9. Durations of the different Wolf-Rayet phases for masses at which at least one of the four (Ledoux, Schwarzschild, rotating, and non-rotating) models has a Wolf-Rayet phase. Left panel: non-rotating models. Right panel: rotating models.
In the text

	Fig. 10. Mass of the 4He core (left panel) and mass ratio of the He and CO cores (right panel) at the last computed stage as a function of the initial mass for the different models. Stars (blue: Schwarzschild, and red: Ledoux) represent non-rotating models, and crosses (green: Schwarzschild, and orange: Ledoux) show rotating models.
In the text

Fig. 11. Number of stars of certain types alive at each time, normalised by the total initial number of stars generated in a synthetic cluster for the case of a burst of star formation at t = 0. The solid lines represent stars computed with the Schwarzschild criterion, and the dashed lines show the same for the Ledoux criterion. Top panels: blue and red supergiants. The vertical axis is logarithmic, and the time axis extends to 60 Myr. Bottom panels: Wolf-Rayet stars. The vertical axis is linear and in units of 10−3, and the time axis extends to 12 Myr (there are no WR stars later on). Left panels: non-rotating models. Right panels: rotating models.

In the text

	Fig. 12. Density maps in the HRD of the Schwarzschild and Ledoux non-rotating populations in the case of continuous star formation (stationary regime). We add the tracks of the 7, 15, 25, and 60 M⊙ models. Left: Schwarzschild model. Right: Ledoux model.
In the text

	Fig. 13. Density maps in the HRD of the Schwarzschild and Ledoux rotating populations in the case of continuous star formation (stationary regime). We add the tracks of the 7, 15, 25, and 60 M⊙ models. Left: Schwarzschild model. Right: Ledoux model.
In the text

Fig. 14. HRD showing isochrones at log(Age [yr]) = 6.7,  6.8,  6.9,  and 7.0 (corresponding to ages of 5, 6.3, 7.9, and 10 Myr) for our four sets of grids, compared with observations of evolved massive stars in the Westerlund-1 cluster. We use the data from Yusof et al. (2022) for the observed stars. Left column: Schwarzschild models. Right column: Ledoux models. Upper row: non-rotating models. Lower row: rotating models.

In the text