A systematic study of super-Eddington layers in the envelopes of massive stars

TL;DR

This study investigates the effects of numerical solutions on modeling super-Eddington layers in massive stars, revealing significant variations in stellar evolution outcomes like expansion and remnant mass.

Contribution

It systematically quantifies how different pragmatic numerical treatments influence the evolutionary properties of massive stars near the Eddington limit.

Findings

Maximum stellar radius varies by up to 2000 R$_\odot$ due to numerical solutions.

Remnant mass differences can reach up to 14 M$_\odot$.

Numerical difficulties cause models with initial masses ≥30 M$_\odot$ to fail before core helium burning ends.

Abstract

The proximity to the Eddington luminosity has been attributed as the cause of several observed effects in massive stars. Computationally, if the luminosity carried through radiation exceeds the local Eddington luminosity in the low-density envelopes of massive stars, it can result in numerical difficulties, inhibiting further computation of stellar models. This problem is exacerbated by the fact that very few massive stars are observed beyond the Humphreys-Davidson limit, the same region in the Hertzsprung-Russell diagram where the aforementioned numerical issues relating to the Eddington luminosity occur in stellar models. One-dimensional stellar evolution codes have to use pragmatic solutions to evolve massive stars through this computationally difficult phase. In this work, we quantify the impact of these solutions on the evolutionary properties of massive stars. We used the stellar evolution code MESA with commonly used input parameters for massive stellar models to compute the evolution of stars in the initial mass range of 10-110 M$_\odot$ at one-tenth of solar metallicity. We find that numerical difficulties in stellar models with initial masses greater than or equal to 30 M$_\odot$ cause these models to fail before the end of core helium burning. Recomputing these models using the same physical inputs but three different pragmatic solutions to treat the numerical instability, we find that the maximum radial expansion achieved by stars can vary by up to 2000 R$_\odot$, while the remnant mass of the stars can vary by up to 14 M$_\odot$ between the sets. These differences can have implications on studies such as binary population synthesis.