Explaining the differences in massive star models from various simulations

The chemical composition of the Sun is often used as a yardstick in computing the metal content of other stars. However, the exact value remains inconclusive and has undergone several revisions since 2004 (see Basu, 2009; Asplund et al., 2021, for an overview). Therefore, different stellar models often make use of different abundance scales.

The BPASS, PARSEC, Geneva and MIST models base their chemical compositions on the Sun, while the BoOST models use a mixture tailored to the sample of massive stars from the FLAMES survey (Evans et al., 2005) with $Z_{Gal}=0.0088$. For stellar winds and opacity calculations, BoOST models use $Z=0.017$ from Grevesse et al. (1996) as the reference solar metallicity. The BPASS models use solar abundances from Grevesse & Noels (1993) with  Z_⊙$=0.02$. The PARSEC models follow Grevesse & Sauval (1998) with revisions from Caffau et al. (2011) and  Z_⊙$=0.01524$. Geneva models use Asplund et al. (2005) abundances with Ne abundance from Cunha et al. (2006) and  Z_⊙$=0.014$ while, finally, the MIST models base their abundance on Asplund et al. (2009) with  Z_⊙$=0.01428$.

For the purpose of comparison here, we use $Z=0.014$ models for each set except for BoOST where we use the Galactic composition, $Z=0.0088$ as the closest match. Iron is an important contributor to metallicity, as numerous iron transition lines dominate both opacity and mass-loss rates, therefore directly affecting the structure of massive stars (e.g., Puls et al., 2000). The stellar models compared here have similar iron content, with the normalised number density ($A_{\rm Fe}={\rm log(N_{Fe}/N_{H})+12.0})$ ranging from $7.40$ (for BoOST models) to $7.54$ (for MIST models).

Stellar mass is a key determinant of a star’s life and evolutionary outcome. It can, however, change as stars lose their outer layers in the form of stellar winds, and through interactions with a binary companion. Consequently, mass loss can affect not only the structure and chemical composition of the star, but is also important in determining its final state (Renzo et al., 2017).

For massive stars, the effects of mass loss are even more pronounced. The mass loss experienced by hot massive stars (O type stars and Wolf–Rayet stars) is known to be line-driven (Lamers & Cassinelli, 1999) while that of cool massive stars (red supergiants, Levesque, 2017) is suggested to be dust-driven. Both types of mass loss are an intensively studied subject. However, the complexity of the problem of atomic and molecular transitions in the wind together with the rarity of stars at these high masses means that the model assumptions are usually based either on a few observations (a small sample of stars) or on what we know about the wind properties of low-mass stars.

All models in this study follow Vink et al. (2000, 2001) for hot wind-driven mass-loss. The PARSEC, Geneva, BPASS and MIST models follow de Jager et al. (1988) for cool dust-driven mass-loss and Nugis & Lamers (2000) for mass loss in the naked helium star phase. The Geneva models further switch to Crowther (2000) for hydrogen-rich stars with ${\rm logT_{eff}/K\leq 3.7}$. They also use the maximum of Vink et al. (2001) and Gräfener & Hamann (2008) for stars with surface hydrogen mass fraction between 0.3 and 0.05, before switching to Nugis & Lamers (2000) when the surface hydrogen mass fraction falls below 0.05. The BoOST models follow Nieuwenhuijzen & de Jager (1990) for cool winds and Hamann et al. (1995, reduced by a factor of 10) for stars with surface hydrogen mass fraction $<\,0.3$. For computing mass-loss rates of stars with surface hydrogen mass fraction between 0.3 and 0.6, BoOST models linearly interpolate between the mass-loss rates of Vink et al. (2001) and Hamann et al. (1995, reduced by a factor of 10).

To account for the dependence of the mass-loss rates on the chemical composition, BoOST, MIST and BPASS models scale the mass-loss rates by a factor of $Z^{0.85}$ (Vink et al., 2001) ¹¹1Note that for some models the factor $Z^{0.69}$ is quoted, depending on whether the dependence of the terminal velocity on Z is explicitly considered or not.. Geneva and PARSEC models also use additional mass-loss as described in Section 3.

Internal mixing processes such as convection and overshooting play an important role in determining both the structure and evolution of massive stars (e.g., see Sukhbold & Woosley, 2014). Similar to mass loss, these processes represent another major source of uncertainty in massive stellar evolution (Schootemeijer et al., 2019; Kaiser et al., 2020). In 1D stellar evolution codes convection is modelled using the Mixing-Length Theory (Böhm-Vitense, 1958, MLT) in terms of the mixing length parameter ${\rm\alpha_{MLT}}$. However, 3D simulations suggest that convection in massive stars might be more sophisticated and turbulent than described by MLT (Jiang et al., 2015).

The BoOST, Geneva, PARSEC and BPASS models used here follow standard MLT (Cox & Giuli, 1968) for convective mixing with mixing length parameter $\alpha_{\rm MLT}$ = $(1.5,1.6,1.74,2.0)$ respectively. MIST follows a modified version of MLT given by Henyey et al. (1965) with $\alpha_{\rm MLT}=1.82$. Convective boundaries in PARSEC, Geneva and BPASS models are determined using the Schwarzschild criterion (Schwarzschild, 1958). BoOST and MIST use the Ledoux criterion (Ledoux, 1947) for determining convective boundaries with semiconvective mixing parameters of 1.0 and 0.1, respectively. For determining convective core overshoot, Geneva and BoOST use step overshooting with overshoot parameter $\alpha_{ov}$ = $(0.1,0.335)$. MIST uses exponential overshooting following Herwig (2000) with $\alpha_{ov}$ = 0.016. PARSEC uses overshoot from Bressan et al. (1981) with $\alpha_{ov}$ = $0.5$. BPASS uses the overshoot prescription from Pols et al. (1998) with $\alpha_{ov}$ = 0.12. For comparison, the rough equivalent in the step overshooting formalism would be 0.2, 0.25 and 0.4 for the MIST, PARSEC and BPASS models respectively (See Choi et al., 2016; Pols et al., 1998; Bressan et al., 2012, for details of each method.)

MIST and PARSEC also include small amounts of overshoot associated with convective regions in the envelope. However, apart from modifying surface abundances, envelope overshoot has a negligible effect on the evolution of the star (Bressan et al., 2012). Rotational mixing also plays an important role in the evolution of massive stars. In fact, the calibration of the free parameters in the stellar codes is often based on their rotating models. For simplicity we only compare non-rotating models for PARSEC, MIST, Geneva and BPASS in this study. Although, for BoOST, in the absence of non-rotating models for stars more massive than 60 M_⊙, we do use slowly rotating (100 km s^-1) models. As shown by Brott et al. (2011a), this small difference in the initial rotation rate is not relevant from the point of view of the overall evolutionary behaviour.

Major input parameters used in each set of models are summarized in Table 1.

The evolution of stars can be easily represented through tracks on the Hertzsprung-Russell (HR) diagram, depicting the evolutionary paths followed by a series of stars. Figure 1 presents the HR diagram of stars of various initial masses from the five simulation approaches. The observational analogue to the Eddington-limit is the Humphreys–Davidson limit or HD-limit (Humphreys & Davidson, 1979). Since the luminosity of a star depends on its mass, more massive stars also are more luminous. This means they can easily exceed the Eddington-limit, develop density inversions and require the use of numerical solutions as discussed in Section 3.

From Figure 1 we see that the tracks of the 25 M_⊙ (or 24 M_⊙ in some cases) stars agree well during most of the evolution. This is because stars of this mass do not exceed the Eddington limit and are thus not affected by the related numerical treatments. The minor differences in their tracks are due to the difference in physical inputs (Section 2) between the simulations. For example, the differences in the position of the main-sequence (MS) hook feature in the HR diagram arise due to the varied extent of convective overshoot used in each set of models.

A 40 M_⊙ star is clearly affected by the numerical treatment employed during the post-main sequence phase of its evolution, as evidenced by the difference in the tracks in the HR diagram shown in Figure 1. More massive stars, i.e, those with initial masses 80/85 M_⊙ and 120/125  M_⊙, can exceed the Eddington-limit in their envelopes while on the MS and therefore their simulations differ significantly from each other. At these masses, the mass-loss rates can be as high as 10^-3 $-$ 10^-4 M_⊙ yr^-1, completely dominating over every other physical ingredient in determining the evolutionary path. While all tracks have been computed with similar prescriptions for wind mass loss (cf. Section 2), the actual rates can be strongly modified by the numerical methods adopted by each code in response to numerical instabilities (Section 3), resulting in vast differences in the tracks.

The ionization released by a stellar population in e.g. a cluster or galaxy is influenced by the contribution of the most massive stars (Topping & Shull, 2015). As shown above however, these are the stars for which the simulations give the most diverse predictions.

To demonstrate this effect, we calculate the ionizing radiation emitted by a simple stellar population, supposing a Salpeter (1955) initial mass function with an upper mass of 120 M_⊙ and a star-forming region of 10⁷ M_⊙ total mass – which is aimed to represent either a typical starburst galaxy, or a young massive cluster in the Milky Way. In Table 2 we list the Lyman photon flux predicted by the individual stellar models analysed here. To simplify the population synthesis calculations, the table provides time averaged values, that is, the photon number flux emitted over the whole evolution is divided by the lifetime. This way the emission coming from the population is estimated by simply weighting the time averaged values by the initial mass function (i.e. without needing to follow the time evolution of the modelled cluster or galaxy). The results of these simple population syntheses are also reported in Tab. 2. In the absence of spectral synthesis models computed for all five sets (cf. Wofford et al., 2016), we have opted to simply use black body estimation. To correct for optically thick winds, we follow the method explained in Chapter 4.5.1 of Szécsi (2016, which relies on ).

We find that, in terms of how much Lyman flux is emitted by a given synthetic population, the model predictions can differ as much as $\sim$18 percent between simulations. This supports earlier findings (e.g. Topping & Shull, 2015) that relying on the ionizing properties of massive stars from evolutionary models should be done with caution. Indeed one should keep in mind that the behaviour of the most massive models, those that dominate the radiation profile of any star-forming region, is weighted with large uncertainties – the source of which is the treatment of the Eddington limit, explained in Section 3.

The radial expansion of a star plays a significant role in determining the nature of binary interaction as it can lead to episodes of mass transfer in close interacting binaries. Recent studies indicate that most of the massive stars occur in binaries (Sana et al., 2012; Moe & Di Stefano, 2017). Therefore, predictions of stellar radii become even more important for determining the binary properties of massive stars.

Figure 2 shows the maximum radial expansion achieved by massive stars from each simulation. For stars with initial mass up to 30 M_⊙, all simulations predict the formation of a red supergiant. The maximum difference in the radius predictions here is $\lesssim$1000 R_⊙. For higher initial masses, the predictions for maximum stellar radii become more divergent as proximity to the Eddington-limit increases and numerical treatments adopted by each code modify the mass-loss rates.

The greatest difference in the maximum radius predictions ($\gtrsim$1000 R_⊙) occurs for stars with initial masses between 40 and 100 M_⊙. Above $\sim$ 100 M_⊙, stars have even higher mass-loss rates which can completely strip a star of its envelope before it can become a red-supergiant. Such stars evolve directly towards the naked helium star phase and have much smaller radii. Therefore, for stars with initial masses more than 100 M_⊙, the difference between the maximum radius predictions by each simulation reduces to $\lesssim$100 R_⊙. The predictions in this mass range seem to further converge into two main groups: PARSEC and MIST represent one group predicting smaller radii compared to the second group which consists of BoOST and BPASS models. The maximum radii in this mass range are predicted by the Geneva models. This is due to the difference in the mass loss rates adopted during the naked helium phase of the star, as explained in the next section.

Stellar evolutionary models provide an easy way of estimating the properties of stellar remnants such as black holes and neutron stars, which are needed in many fields including supernova studies (e.g. Aguilera-Dena et al., 2018; Raithel et al., 2018), gamma-ray bursts progenitors (e.g. Yoon et al., 2006; Szecsi, 2017), and gravitational-wave event rate predictions (e.g. Stevenson et al., 2019; Mapelli et al., 2020).

Following Belczynski et al. (2010) we show in Figure 3 how the uncertainties in the models we compare here also pose a challenge for the predictions of remnant properties. Remnant masses have been calculated from the carbon-oxygen (CO) core mass and the total mass of the star at the end of the core helium-burning phase using the prescription of Belczynski et al. (2008) (same as the StarTrack prescription in Fryer et al., 2012).

The Geneva models do not provide information on core masses in their publicly available dataset. Therefore, the final CO-core mass for these models has been taken from Georgy et al. (2012, for stars with initial mass up to 120 M_⊙) and Yusof et al. (2013, for stars with initial mass more than 120 M_⊙). Note that in the Geneva models from Yusof et al. (2013), CO-core mass is defined as the mass of the core where the sum of mass fraction of carbon and oxygen exceeds 75 percent, and is different to the definition used by Georgy et al. (2012).

The remnant masses are heavily influenced by the modelling assumptions (cf. Section 2) and the numerical methods (cf. Section 3) especially above M_ini $=$ 40 M_⊙ where the most massive black holes are predicted. For stars with initial masses between 9 and 120 M_⊙, we find that the mass of the black holes predicted by the different sets of models can differ by $\sim$ 20 M_⊙. The maximum black hole mass varies from about 20 M_⊙ for BPASS models to about 35 M_⊙ for MIST and PARSEC models, and 32 M_⊙ for the BoOST models. BPASS models consistently predict the lowest values of remnant mass for most of the massive stars while predictions from the BoOST, MIST and PARSEC models peak at 60, 90 and 120 M_⊙ before flattening out at the higher initial masses.

For more massive stars, the variation in the remnant mass between BPASS, BoOST, MIST and PARSEC models reduces to about 10 M_⊙. An interesting behaviour is shown by the Geneva models, which predict one of the lowest remnant masses for stars up to 120 M_⊙, reaching a maximum of only 28 M_⊙ at 120 M_⊙. These models however, predict the highest values of remnant masses beyond 120 M_⊙. At their farthest, for a model with initial mass of 200 M_⊙, the predictions between the Geneva models and other models can be as high as 30 M_⊙. Similar variability is found in the core mass and the final total mass of the star (from which the remnant masses have been calculated).

The significantly higher remnant masses predicted by the Geneva models for stars with initial mass beyond 120 M_⊙can be explained as follows. Due to their high luminosity, stars more massive than 100–120 M_⊙ rapidly lose mass during the main-sequence phase and directly evolve towards the naked helium phase (cf. Figure 1). The mass-loss prescriptions from both Nugis & Lamers (2000) and Hamann et al. (1995) predict the highest mass-loss rates for stars occur during this naked helium phase. BPASS, MIST and PARSEC models switch to mass-loss rates from Nugis & Lamers (2000) when the mass fraction of hydrogen at the surface of the star falls below 0.4, 0.4 and 0.5 respectively. BoOST models linearly interpolate between the mass-loss rates of Vink et al. (2001) and Hamann et al. (1995, reduced by a factor of 10) for stars with surface hydrogen mass fraction between 0.3 and 0.6, before completely switching to the mass-loss rate from Hamann et al. (1995, reduced by a factor of 10) for stars with surface hydrogen mass fraction less than 0.3.

Geneva models, on the other hand, switch to using mass-loss rates from Nugis & Lamers (2000) only when the mass fraction of hydrogen at the surface of the star falls below 0.05. For stars with surface hydrogen mass fraction between 0.3 and 0.05, they use the maximum of Vink et al. (2001) and Gräfener & Hamann (2008), which predict lower mass-loss rates compared to both Nugis & Lamers (2000) and Hamann et al. (1995) (See section 2.1 of Yusof et al., 2013). Thus, they do not lose as much mass as other models during the naked helium star phase and end with higher total mass and thus with the higher remnant mass.

Note that the Belczynski et al. (2008) prescription is one of several methods for predicting the remnant properties of the stars. Other methods for calculating the remnant masses may predict higher or lower values. For example, the remnant mass calculations based on the binding energy of the star (e.g. see Eldridge et al., 2017) are generally lower than those predicted here. However, all the models we study have near-solar metallicity and therefore rather high mass-loss rates, none of them stays massive enough at the end of their lives to undergo pair-instability.