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, BoOST, Geneva, MIST, and PARSEC 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$_\odot$. 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 percent, the maximum radial expansion of a star can differ between 100-1600 R$_\odot$, and the mass of the stellar remnant can vary up to 20 M$_\odot$ 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$_\odot$ should be made with caution.