Explaining the differences in massive star models from various simulations

The evolution of massive stars is the basis of several astrophysical investigations, from predicting gravitational-wave event rates to studying star formation and stellar populations in clusters. However, uncertainties in massive star evolution present a significant challenge when accounting for these models' behaviour in stellar population studies. In this work, we present a comparison between five published sets of stellar models from the BPASS (Binary Population and Spectral Synthesis), BoOST (Bonn Optimized Stellar Tracks), Geneva, MIST (MESA Isochrones and Stellar Tracks), and PARSEC (PAdova and TRieste Stellar Evolution Code) simulations at near-solar metallicity. The different sets of stellar models have been computed using slightly different physical inputs in terms of mass-loss rates and internal mixing properties. Moreover, these models also employ various pragmatic methods to overcome the numerical difficulties that arise due to the presence of density inversions in the outer layers of stars more massive than 40 M_⊙. These density inversions result from the combination of inefficient convection in the low-density envelopes of massive stars and the excess of radiative luminosity to the Eddington luminosity. We find that the ionizing radiation released by the stellar populations can change by up to 18 per cent, the maximum radial expansion of a star can differ between 100 and 1600 R_⊙, and the mass of the stellar remnant can vary up to 20 M_⊙ between the five sets of simulations. We conclude that any attempts to explain observations that rely on the use of models of stars more massive than 40 M_⊙ should be made with caution.