© The Authors 2022

1. Introduction

Stars more massive than about 9 M_⊙ are key to several astrophysical processes. During their lives, they enrich their surroundings with ionizing flux and nuclear-processed material while altering the dynamics of their host systems. At the end of their lives, these massive stars again expel a copious amount of radiation and metal-rich matter in the form of supernovae, leaving behind compact remnants: neutron stars and black holes. Furthermore, mergers of these compact remnants result in gravitational wave (Abbott et al. 2016, 2017) emission and can also lead to the formation of rare elements (Kasen et al. 2017). Therefore, a better understanding of how these stars evolve is crucial in comprehending their contribution to the evolution of star clusters and galaxies.

The evolution of massive stars is typically modeled using one-dimensional (1D) stellar evolution codes. In the last few decades, these stellar evolution codes have progressed a lot and so has our understanding of massive stars. Together with the advances in our observing capabilities (Evans et al. 2011; Simón-Díaz et al. 2015; Wade et al. 2014; Abbott et al. 2016), the development of sophisticated numerical methods for simulating physical processes and newer input data in the form of opacity tables and nuclear reaction rates has led to the development of modern and improved stellar structure and evolution codes (Langer 2012; Ekström et al. 2020).

Even with these new capabilities, 1D modeling of massive stars is limited by a number of approximations. Several evolutionary properties of massive stars such as mass-loss rates (Smith 2014; Renzo et al. 2017), nuclear reaction rates (Heger et al. 2002; Fields et al. 2018), and rotation (Heger et al. 2000; Maeder 2009) remain uncertain. The high mass of these stars makes it feasible for sophisticated physical processes to operate in these stars, but their short lives makes it difficult to obtain observational constraints.

One such process is the treatment of convective transport of energy through the mixing length theory (MLT; Böhm-Vitense 1958). In this theory, energy is transported through fluid elements supported by buoyancy forces. These elements travel over a radial distance known as the mixing length, after which they dissolve in their surroundings. MLT assumes hydrostatic equilibrium in stars and only depends on local conditions (i.e., local values of pressure, density, etc.), without taking other parts of the star into account. Time-dependency and nonlocal treatments are included through ad hoc methods such as convective overshoot, semiconvection, and diffusion (e.g., Renzini 1987; Kippenhahn et al. 2012).

The simplicity of MLT makes it a popular choice for many stellar evolution codes. While MLT gives fairly good results in the deep interiors of stars where density is high and convection is nearly adiabatic (with negligible radiative losses), its limitations start becoming apparent in low-density environments where the convection is highly superadiabatic and prone to radiative losses (Maeder 2009). For example, in the low-density envelopes of massive stars, convection, as given by MLT, is inefficient. Furthermore, changes in the elemental opacity as the star evolves can cause the radiative luminosity to exceed the local Eddington luminosity and can lead to the formation of density and gas pressure inversions in the subsurface layers.

The proximity to the Eddington luminosity and the associated density inversions have been attributed as the source of several instabilities in massive stars, for example, the dynamical instability (Stothers & Chin 1993), the convective instability (Langer 1997), and the strange-mode instability (Saio et al. 1998, 2013). Observationally, these have been linked to stellar variability phenomena such as stochastic low-frequency photometric variability (Pedersen et al. 2019; Bowman 2020), spectroscopic macroturbulence (Simón-Díaz & Herrero 2014; Simón-Díaz et al. 2017), and episodic mass ejection behavior in luminous blue variables (LBVs: Bestenlehner et al. 2014; Gräfener 2021).

From a numerical perspective, the presence of density inversions in the inflated envelopes of supergiant stars requires 1D stellar evolution codes to take prohibitively short time steps (on the order of hours and minutes), leading to convergence issues (Maeder 1987; Alongi et al. 1993; Paxton et al. 2013). Evolving stars past these numerical difficulties in the supergiant phase has been a long-standing challenge for the 1D stellar evolution approach. As shown in Agrawal et al. (2022; hereafter Paper I), stellar evolution codes often resort to numerous pragmatic solutions to evolve stars, such as enhancing the convective efficiency (e.g., Ekström et al. 2012) or limiting temperature gradients such that the density gradient is always positive (e.g., Chen et al. 2015). While these solutions help evolve stars through numerically difficult phases of evolution, they can also modify their surface behavior, such as radius evolution and mass-loss rates.

The different solutions used by 1D codes and their interplay with other physical parameters of massive stars can alter the dynamics of stellar evolution, thereby adding a potential bias to any study aiming to determine the properties of these input parameters from the evolution of a star. While other physical uncertainties in the evolution of massive stars have received considerable attention in a number of studies, the impact of numerical issues has not been explored to the same extent.

In Paper I, the comparison of stellar models from five different codes revealed large differences in the evolutionary behavior of stars more massive than 40 M_⊙ around solar metallicity (at Z = 0.014). The maximum radial expansion predicted by different stellar models varied by more than 1000 R_⊙ and the predictions of remnant mass varied by 20 M_⊙. However, the stellar models also had different physical inputs besides the pragmatic treatment of density inversions arising due to the Eddington luminosity, making it difficult to untangle the impact of this process from other inputs. There is therefore a need for a systematic study of the impact of these solutions within a single code and single set of assumptions.

In this work, we perform a study of the impact of density inversions on the evolution of massive stars up to 110 M_⊙ using consistent input parameters. Since Paper I focused on solar metallicity stars only, we chose a metallicity ten times lower than solar here to demonstrate the impact of density inversions at a metallicity relevant to the progenitors of current gravitational wave observations (e.g., Stevenson et al. 2017). As we show here, the different pragmatic solutions used by 1D codes can have a non-negligible impact on the evolutionary properties of massive stars. These differences are important, as they can help us explain the formation of gravitational wave progenitors and other observations of stellar populations.

The paper is organized as follows. We provide an overview of the physics related to density inversions in Sect. 2. In Sect. 3, we describe our standard or default set of models with MESA and discuss the effect of density inversions on their completeness. We recompute the models that fail to reach the end of evolution, using three different pragmatic solutions in Sect. 4 and compare the impact of these solutions in predicting the different evolutionary properties of massive stars in Sect. 5. In Sect. 6, we compare the stellar models from Sect. 4 with the observations of massive stars and conclude our study in Sect. 7.

2. Physics of density inversions in massive stars

In this section, we describe the conditions for the formation of density and gas pressure inversions in stellar envelopes and their impact on modeling the evolution of massive stars. For a spherically symmetric star containing mass m(r) inside radius r and with radiative opacity κ(r) and density ρ(r), the luminosity that can be carried by radiative transport of energy is given by

$$\begin{array}{c}\hfill L_{\mathrm{rad}}(r)=-\frac{4\pi r^{2}c}{\rho (r)\kappa (r)}\frac{\mathrm{d}{P}_{\mathrm{rad}}}{\mathrm{d}r},\end{array}$$(1)

where P_rad denotes the radiation pressure and c is the speed of light.

The Eddington luminosity gives the maximum value of luminosity that can be transported by radiation while maintaining hydrostatic equilibrium (Eddington 1926). The expression for the Eddington luminosity is given by

$$\begin{array}{c}\hfill L_{\mathrm{Edd}}(r)=\frac{4\pi cGm(r)}{\kappa (r)},\end{array}$$(2)

where G represents the gravitational constant. The total pressure P in the star is the sum of the radiation pressure, P_rad and the gas pressure P_gas. Using Eqs. (1) and (2) and the equation of hydrostatic equilibrium dP/dr = −Gm(r)ρ(r)/r², the ratio of radiative luminosity to the Eddington luminosity (Eddington factor) can be defined as

$$\begin{array}{c}\hfill \frac{L_{\mathrm{rad}}}{L_{\mathrm{Edd}}}=\frac{\mathrm{d}{P}_{\mathrm{rad}}}{\mathrm{d}P}=[1+\frac{\mathrm{d}{P}_{\mathrm{gas}}}{\mathrm{d}{P}_{\mathrm{rad}}}].\end{array}$$(3)

Normally, the luminosity transported by radiation (L_rad) is less than the Eddington luminosity (L_Edd) inside a star, and the density and gas pressure of stellar material decrease with the stellar radius (dρ/dr <  0, dP_gas/dr <  0). However, during the evolution of the star, changes in the ionization states of various elements can lead to an increase in the opacity κ(r). In the low density (ρ ≪ 1 g cm⁻³) radiation pressure dominated (P_rad/P ≈ 1) envelopes of massive stars, an increase in opacity can reduce the local Eddington luminosity below the radiative luminosity $(\frac{L_{\mathrm{rad}}}{L_{\mathrm{Edd}}}>1)$. Since dP_rad/dr is always less than 0, when $\frac{L_{\mathrm{rad}}}{L_{\mathrm{Edd}}}>1$, Eq. (3) implies dP_gas/dr >  0, that is, a gas pressure inversion.

From the ideal gas equation of state, P_gas can be expressed as a function of ρ and P_rad. Using Eq. (3), the condition for density inversion (dρ/dr >  0) can therefore be written as

$$\begin{array}{c}\hfill \frac{L_{\mathrm{rad}}}{L_{\mathrm{Edd}}}>{[1+{\left(\frac{\partial P_{\mathrm{gas}}}{\partial P_{\mathrm{rad}}}\right)}_{\rho }]}^{-1}\end{array}$$(4)

(Joss et al. 1973; Paxton et al. 2013).

Since ${\left(\frac{\partial P_{\mathrm{gas}}}{\partial P_{\mathrm{rad}}}\right)}_{\rho }$ is a always positive, $\frac{L_{\mathrm{rad}}}{L_{\mathrm{Edd}}}>1$ implies a density inversion in stars.

The effect of density inversions on the evolution of massive stars is complex and remains an active field of research (Mihalas 1969; Érgma 1971; Langer 1997; Owocki et al. 2004; Cantiello et al. 2009; Sanyal et al. 2015). Gräfener et al. (2012) found that an increase in the radiative luminosity near the opacity peak due to the iron-group elements (at ∼2 × 10⁵ K) leads to the formation of an “inflated envelope” containing density inversions (also see Ishii et al. 1999; Petrovic et al. 2006; Köhler et al. 2015). In this state, the star has an extended radiative envelope with a relatively small convective core. As pointed out by Sanyal et al. (2015), these inflated stars are different from the classical red supergiants, as envelope inflation does not require hydrogen shell burning and can even occur while the star is on the main sequence.

For hydrogen-rich stars, the formation of the inflated envelope reduces the opacity and therefore the Eddington factor, L_rad/L_Edd, stays less than or close to the unity. However, as the envelope expands, its outer layers become sufficiently cool to encounter even steeper opacity bumps due to helium ionization (at ∼2 × 10⁴ K) and hydrogen recombination (at ∼10⁴ K). Further expansion no longer reduces L_rad/L_Edd, and a density inversion forms in the outer layers of the star.

The presence of these inversions in the convectively-inefficient and radiation-dominated envelope of massive stars (P_gas/P ≪ 1) leads to numerical instabilities (Maeder 1987; Maeder et al. 1992) and an additional mode of energy transport (such as turbulent convection) may be needed to carry the energy (Jiang et al. 2015; Schultz et al. 2020). Without proper treatment of these instabilities, 1D stellar evolution codes find it difficult to obtain the solution of the stellar structure equations. They are, therefore, forced to adopt exceedingly small time steps and struggle to complete the evolution of these stars (Paxton et al. 2013).

3. Stellar models: The standard set

Module for Experiments in Stellar Astrophysics (MESA: Paxton et al. 2011, 2013, 2015; Paxton et al. 2018, 2019) is an advanced one dimensional stellar structure and evolution code. MESA solves the coupled differential equations of the stellar structure simultaneously with the energy transport equation for all mass cells from the surface to the center. In this work, we make use of version 11701 of MESA.

3.1. Physical ingredients

Convection is treated as a diffusive process with convective boundaries determined by the Ledoux criterion. To account for convection due to varying opacity in layers, we use the Henyey et al. (1965) treatment of MLT with mixing length parameter, α_MLT = 1.82. Semiconvection is implemented using Langer et al. (1985) with α_{semiconvection} = 1.0 while a small amount of Thermohaline mixing is used (α_th = 2.0) following Kippenhahn et al. (1980). For convective overshooting we follow the calibration adopted by Brott et al. (2011): step overshoot with f = 0.335 and f0 = 0.05. We only apply overshooting in the core and not in the convective envelope.

Mass loss is modeled using the commonly used MESA “Dutch” scheme (Glebbeek et al. 2009) with the mass loss scaling factor, η_Dutch = 1.0. The scheme follows Vink et al. (2000, 2001) for hot winds, de Jager et al. (1988) for cool winds and Nugis & Lamers (2000) for Wolf-Rayet winds.

Opacities are calculated using OPAL (Iglesias & Rogers 1993, 1996) Type I and Type II (at the end of hydrogen burning) using Asplund et al. (2009) photospheric abundances. Electron screening is included for both the weak and strong regimes using the “extended” option in MESA, which computes the screening factors by extending the classic Graboske et al. (1973) method with that of Alastuey & Jancovici (1978), and adopting plasma parameters from Itoh et al. (1979) for strong screening. Boundary conditions at the surface of the star are calculated using the “Eddington_grey” option of MESA, which uses the Eddington T-tau integration (Eddington 1926) to obtain temperature and pressure in the outer layer of the stars.

Nuclear reaction networks can have a nontrivial effect on the evolution of massive stars (Farmer et al. 2016), especially the inclusion of elements like iron and nickel (Nabi et al. 2019). Hence, we make use of an extensive grid of isotopes for the nuclear reaction network in our models to closely follow the evolution of massive stars. The network has been chosen to match globular cluster observations and has 72 elements, including Mg, Li and Fe. The reaction rates are determined using the Jina Reaclib database (Cyburt et al. 2010).

We use high spatial resolution with mesh_delta_coeff = 0.5 and the maximum relative cell size max_dq = 5 × 10⁻⁴. Model-to-model structure variation is kept modest with varcontrol between 7 × 10⁻⁴ and 10⁻³, and additional constraints are used to limit time steps where necessary.

We compute the evolution of stars in the mass range of 10–110 M_⊙ at metallcity, Z  =  0.00142 ([Fe/H] = −1). All stars start from the pre-main-sequence (PMS) with uniform composition and the goal is to evolve each model until the point of carbon depletion in the core (X_c  ≤  10⁻²). Following Choi et al. (2016), we adopt solar scaled abundances from Asplund et al. (2009), with Solar metallicity Z_⊙ = Z_{⊙, protosolar} = 0.0142. The primordial helium abundance Y_p is taken to be 0.249 while the protosolar helium abundance Y_{⊙, protosolar} = 0.2703. Initial hydrogen (X) and helium abundances (Y) for PMS models are calculated using the following formula

$$\begin{array}{cc}& Y=Y_{\mathrm{p}}+\left(\frac{Y_{\odot ,\mathrm{protosolar}}-Y_{\mathrm{p}}}{Z_{\odot ,\mathrm{protosolar}}}\right)Z,\hfill \end{array}$$(5)

$$\begin{array}{c}\hfill \\ \hfill & X=1-Y-Z.\hfill \end{array}$$(6)

Models are evolved without mass loss from the PMS until the zero-age main sequence (ZAMS), defined as when the central hydrogen abundance reduces by 1 percent of the initial value. Evolution is then restarted (including mass loss) using the stellar model saved at ZAMS, and continues until either the termination condition is reached (X_c ≤ 10⁻²) or the end of 96 CPU-hrs (running parallel on four cores for 24 h).

The set of models evolved using the physical inputs described above fail to reach the end of carbon burning in the core within the allocated time (24 h) for stars more massive than 20 M_⊙. Hereafter, we refer to this set as the “standard set” and label it as “NoModifier” in the figures, since the models were computed without any modifications to facilitate their evolution. The results for the standard set are unaffected by the increase in temporal and spatial resolution.

3.2. Results

Figure 1 presents the stellar tracks in the Hertzsprung-Russell (HR) diagram for our standard set of models. In the left panel, the tracks are colored by the maximum of L_rad/L_Edd, while the right panel shows the minimum of P_gas/P_total for each stellar model. We see that in these models L_rad/L_Edd approaches unity and can even exceed it by factors of a few during the evolution as the stars encounter opacity bumps in their envelopes. These opacity bumps are due to the partial ionization states of iron and helium, as well as the recombination of hydrogen. However, for 10 and 20 M_⊙ stars, P_gas/P_total remains high enough (>0.5) for their evolution to proceed uninterrupted (see Sect. 2).

Fig. 1. Hertzsprung-Russell (HR) diagram showing stellar models from the standard set colored by the maximum of Lrad/LEdd (left panel) and the minimum of Pgas/Ptotal (right panel). The blue cross marks the end of core hydrogen burning and the red cross marks the end of core helium burning (where applicable). The brown dashed line in the left panel denotes the position of the observational Humphreys & Davidson (1979) limit beyond which few stars are observed, while the gray dashed line signifies the luminosities of the brightest red supergiants as inferred by Davies et al. (2018). Large values of Lrad/LEdd and small values of Pgas/Ptotal in the envelopes of stars with initial masses 30 M⊙ and above causes the evolution of these stars to become halted at log Teff/K ≈ 3.7, and their models fail to even finish core helium burning.

Models with initial masses between 30 and 110 M_⊙ fail during core helium burning at similar effective temperatures, log T_eff/K ≈ 3.7. In these models, the stellar envelope inflates in response to the iron-opacity peak, thereby preventing density inversion. However, as the star expands and cools, the minima in gas pressure fraction at the opacity peak decreases as highlighted in the right panel of Fig. 1. The evolution of these models progresses smoothly until stars encounter the opacity peak due to hydrogen and helium in their subsurface layers at log T_eff/K ≈ 4.0. The envelope inflation there is not sufficient to prevent the radiative luminosity from exceeding the Eddington luminosity, which leads to density inversions.

To elaborate on the conditions in the stellar interior, the temperature profiles for a 110 M_⊙ star at the ZAMS (at log T_eff/K = 4.77), and at the end of the track (at log T_eff/K = 3.72, corresponding to the final model reached after 24 h of computation), are shown in Figs. 2 and 3. In both these Figures the x-axis is limited to temperatures less than 5 × 10⁶ to show the impact of the various opacity bumps on the evolution of the stellar model.

Fig. 2. Temperature profile of a 110 M⊙ star at the ZAMS evolved using the standard set. Opacity peaks due to partial ionization states of iron are present at 1.5 × 106 and 1.8 × 105 K. However, Lrad/LEdd is less than 1 throughout the star and the evolution of the star proceeds smoothly. For clarity, only the outermost layers of the star (with T ≤ 5 × 106 K) containing the various opacity peaks have been plotted. See Sect. 3.2 for details of each panel.

Fig. 3. Temperature profile of the outermost layers of a 110 M⊙ star at the end of the simulation, evolved using the standard set. The high luminosity of the star combined with the peak in opacity due to hydrogen and helium ionization around 104 K causes Lrad/LEdd to exceed unity. This causes density and gas pressure inversions at 3 × 104 K. However, the low density of the environment renders convection inefficient and the star struggles to evolve despite reaching supersonic convective velocities.

In panel (a) of Fig. 2, L_rad/L_Edd shows a small increase corresponding to the opacity peaks due to partial ionization of iron at 1.8 × 10⁵ and 1.8 × 10⁵ K (panel b). In response to the increase in L_rad/L_Edd, gas pressure dips a little (visible as a minimum in P_gas/P_total in panel e), although, density consistently decreases outward (shown in panel c). The bottom panel (g) shows the difference between the actual temperature gradient and the adiabatic temperature gradient inside the star. This difference, ∇_T − ∇_ad, is known as superadiabaticity (see Sect. 4.3 for details). At the location of the opacity peak, superadiabaticity is positive but small and the convective velocity (v_conv) is less than the isothermal sound velocity (v_sound) (panel f), signifying efficient convection. The specific entropy (S/N_Ak_B) in panel (d) at the base of the convective region is small and the evolution of the star proceeds smoothly.

In panel (b) of Fig. 3, in addition to the iron opacity peak, opacity peaks corresponding to partial ionization of helium and hydrogen recombination can be seen at 3.5 × 10⁴ and 10⁴ K. In contrast to Fig. 2, the increase in opacity due to hydrogen and helium ionization increases L_rad/L_Edd above one, causing the density inversion near the surface of the star. However, the high value of superadiabaticity (∇_T − ∇_ad ≈ 10) renders convection inefficient and prone to radiative losses. Consequently, the specific entropy shows a steep increase (S/N_Ak_B ≥ 300) while the gas pressure fraction becomes very low (P_gas/P_total ≲ 0.01) at the base of the convective envelope near the iron opacity bump. The convective velocity increases to increase the amount of flux convection can carry, becoming more than the local sound speed. Nonetheless, convection remains inefficient and the high value of entropy causes issues with the convergence of the models (Schwab 2019). Struggling to find a solution, MESA retries with smaller time steps until they become on the order of days.

The evolution of the star is essentially halted until it is able to get rid of the density inversion as mass loss slowly chips away the outer convective layers. However, evolving the star this way requires a lot of computational time (e.g., it takes ≈200 CPU-hrs for a 40 M_⊙ star) which might not be feasible.

4. Stellar models: The pragmatic solutions

In the absence of efficient convection, the radiation-dominated envelopes of massive stars in 1D modeling are prone to numerical difficulties (Stothers & Chin 1979; Maeder 1987). Therefore, stellar evolution codes adopt various pragmatic solutions to compute the evolution of massive stars beyond these numerically difficult points (Alongi et al. 1993; Ekström et al. 2012; Paxton et al. 2013). The exact solution differs from code to code; however, they can be summed up into three main categories: using higher mixing length to increase the efficiency of the mixing process, using higher mass-loss rates to remove layers with numerical instabilities, or suppressing the numerical instability by limiting the temperature gradient and thereby suppressing density inversions. We explore each of them in detail in the following subsections.

4.1. Enhancing mixing length

In the convective transport of energy by MLT, the mixing length, l travelled by a fluid element before dissolving in the surroundings is given in terms of the local pressure scale height H_p and mixing length parameter α_MLT such that l = H_v × α_MLT ¹. Therefore a higher value of l implies better convective transport of energy and can reduce the value of the Eddington factor, thereby preventing numerical instabilities.

In the low-density envelopes of massive stars, density scale height, H_ρ, is much greater than the pressure scale height, H_p. Thus, using H_ρ instead of H_p in calculating mixing length results in larger l and can help overcome density inversion in the stellar models (e.g., Érgma 1971; Ekström et al. 2011). Similarly, for a given H_p, a higher value of α_MLT also implies a higher l and again better convective efficiency. We use the latter approach in this work. Beginning with α_MLT = 1.82 (used in the standard set), we compute a series of stellar models with α_MLT = 3.0, 3.64, 4.0, 5.0, 5.46, 7.28, 8.0 for stars with initial mass 30 M_⊙ and above.

We find that for α_MLT≥ 5.46, which is three times the value used in the standard set, the models are able to evolve without any numerical instabilities until carbon depletion in the core.

The top panel of Fig. 4 shows the evolutionary tracks evolved with α_MLT= 5.46. From the figure we see that L_rad/L_Edd still exceeds one in the stellar envelope, however, this does not limit the time steps of the computation of the stellar models. The reason for this can be understood from Fig. 5, where we show the temperature profile of 110 M_⊙ star at a similar location in the HR diagram where the computation for the 110 M_⊙ star from the standard set became stuck. In this Figure (as well as in Figs. 6–7) the x-axis is limited to temperatures less than 5 × 10⁶ to show the impact of the various opacity bumps on the evolution of the stellar model. Similar to Fig. 3, the stellar profile contains opacity peaks due to ionization of iron, helium and hydrogen (panel b), resulting in excess L_rad/L_Edd (panel a) and the density and gas pressure inversions near the surface (panel c and e). However, higher convective velocities resulting from higher α_MLT imply that the convective fluid element travels faster and transports more energy before it leaks out due to radiative losses. The superadiabaticity in the outer layers is also smaller compared to the standard case (∇_T − ∇_ad ≈ 2), meaning radiative losses are also smaller. The overall convective flux is, therefore, higher than the standard case. The gas pressure fraction is non-negligible (P_gas/P_total >  0.01) and the specific entropy is small (S/N_Ak_B <  100). Thus the time steps remain large enough to efficiently compute the evolution of the star until the end of carbon burning.

Fig. 4. HR diagrams showing the evolutionary tracks graded with the maximum of Lrad/LEdd for stars in the mass range 30–110 M⊙, computed with the model variations described in Sect. 4. The top panel shows the tracks computed with enhanced mixing (αMLT = 5.46), the middle panel shows tracks with enhanced mass loss (ηDutch = 8.0 whenever L exceeds LEdd) and the bottom panel shows tracks computed with MLT++. Similar to Fig. 1, blue and red crosses mark the end of core hydrogen burning and core helium burning while the brown and gray dashed lines signify the Humphreys & Davidson (1979) limit and the Davies et al. (2018) limit respectively.

Fig. 5. Temperature profile of the outermost layers of a 110 M⊙ star computed with the mixing length parameter, αMLT = 5.46. Similar to Fig. 3, a density inversion can be seen in panel (c). However, higher convective velocity (panel f) and smaller superadiabaticity (panel g) resulting from the higher value of αMLT helps to keep the specific entropy small (panel d) and time steps large enough to evolve the model without numerical instabilities.

Fig. 6. Temperature profile of the outermost layers of a 70 M⊙ star computed with enhanced mass loss, as described in Sect. 4.2. Opacity peaks due to the ionization states of hydrogen, helium and iron and associated density inversions can be seen in panel (b) and panel (c) respectively. However, high mass-loss rates remove the outer layers containing the hydrogen and helium opacity peaks, moderating the specific entropy at the base of the convective envelope and the model is able to evolve until completion with reasonably large time steps.

Fig. 7. Temperature profile of the outermost layers of a 110 M⊙ star computed with the MLT++ method of MESA. The method artificially reduces the difference between ∇T and ∇ad in the stellar envelope, acting as an source of additional envelope mixing (see Sect. 4.3 for details). The small envelope mixing used here is not enough to suppress density inversions completely (panel c) but it does help to keep the specific entropy small (panel d), gas pressure fraction non-negligible (panel e) and therefore time steps reasonable.

Increasing convective efficiency in this way helps compute the evolution of stars until core carbon burning. However, it also changes the effective temperature of these stars and makes them appear bluer in the HR diagram (Maeder 1987; Kippenhahn et al. 2012). Comparing the tracks with enhanced mixing length (top panel of Fig. 4), with the standard set (left panel of Fig. 1), we find that stellar models with enhanced mixing length are indeed limited to log T_eff/K ≈ 3.73 in the HR diagram which is greater than the minimum log T_eff/K ≈ 3.63 reached by models in the standard set. We discuss this further in Sect. 5.

4.2. Enhancing mass loss

Winds of hot massive main-sequence stars are radiation-driven winds. Due to their high luminosity, massive stars can generate a high number of photons which can be scattered via ions, transferring momentum with them which then accelerates material outward that can escape the gravitational potential of the star. As the star evolves to the cool red-supergiant phase, these winds transition to being dust-driven where they are generated by the interaction of photons with dust grains instead of ions. Massive stars can lose a substantial amount of mass through stellar winds, even their entire hydrogen envelope to become naked helium stars. The mass loss in naked helium stars is again driven by radiation pressure and can be higher by a factor of ten or more than their hydrogen-rich counterparts (Smith & Tombleson 2015)

Although the mass encompassed in the outer layers of the star which contains the density inversions is small (< -pagination1 M_⊙), the classical mass-loss rates described above are unable to get rid of it as the convective envelope quickly expands to lower temperatures at which hydrogen and helium opacity peaks are present, and density inversions are thereby sustained. However, if the mass-loss rates are high enough, convection will not be able to replenish the envelope and the star will contract and move to higher effective temperatures. High mass-loss rates can therefore help remove outer layers in the stellar model containing the density inversion, helping the stellar model avoid numerical instabilities (Petrovic et al. 2006; Cantiello et al. 2009).

The mass-loss rates for massive stars are highly uncertain (Renzo et al. 2017; Björklund et al. 2021) and can be affected by instabilities and processes other than those described at the beginning of this section. For example, Gräfener & Hamann (2008) and Vink (2011) have shown the contribution of optically thick (clumped) winds in the presence of subsurface opacity bumps in massive stars. Moreover, massive stars exceeding the classical Eddington limit can also experience episodes of much stronger mass-loss rates (up to 10⁻³ M_⊙ yr⁻¹), known as LBV eruptions or super-Eddington winds (Lamers & Fitzpatrick 1988; Humphreys & Davidson 1994; Smith et al. 2004). Despite being the stronger contributor to mass-loss rates for massive stars, the exact rates and the mechanism behind the LBV eruptions remain disputed (Puls et al. 2008; Smith 2014; Owocki 2015) and therefore most stellar evolution models often exclude their contributions in the mass-loss rates.

To determine the impact of enhanced wind mass loss on the convergence properties of models that fail to evolve in the standard set, we recomputed these models with increased wind mass-loss rates for different evolutionary phases until the models are able to evolve without numerical difficulties. We find that setting the mass loss scaling factor to η_Dutch ≥ 8.0, that is, at least eight times the mass loss used in the standard set, whenever the stellar luminosity exceeds the local Eddington luminosity in the envelope, leads to smooth evolution of the stellar models with initial mass greater than or equal to 30 M_⊙.

We also find that this mass-loss enhancement is only required for stars with a hydrogen-rich envelope (X_surf ≥ 0.4). Although naked helium stars can also exceed the Eddington limit and can develop density inversions in their envelopes, their evolution proceeds uninterrupted for our models without any numerical instabilities.

The models with enhanced mass loss as described above, with η_Dutch = 8.0, are shown in the middle panel of Fig. 4. A 110 M_⊙ star with enhanced mass loss is able to evolve to the end of carbon burning, while the maximum of L_rad/L_Edd remains close to unity throughout the evolution. Here, the high mass-loss rates and subsequent changes in surface composition due to the removal of H-rich layers reduce the inflation effect (Gräfener 2021) and cause the star to lose its envelope before the steep opacity peaks due to hydrogen and helium can develop. However, for models less massive than 70 M_⊙, the maximum of L_rad/L_Edd can be up to 8, similar to models in the standard set. The high values of L_rad/L_Edd do lead to the formation of density inversions, as shown in the temperature profile of the 70 M_⊙ star in Fig. 6. The enhanced mass-loss rates (up to 10⁻⁴ M_⊙ yr⁻¹) are unable to prevent the envelope inflation and the models do encounter hydrogen and helium opacity peaks in the cooler outer layers. However, the mass-loss rates are high enough to remove these outer layers containing density inversions before the specific entropy at the base of the envelope becomes too large, or the time steps become too small. Thus, all models are able to evolve smoothly to completion without any difficulty.

Enhancing mass loss this way helps compute the evolution of massive stars all the way through to carbon depletion in their core. However, it also influences the structure and evolutionary properties of the stars as explained in Sect. 5. As we discuss further in Sect. 6, such enhanced mass-loss rates are in tension with recent observational and theoretical estimates for massive stars (e.g., Björklund et al. 2021; Hawcroft et al. 2021), but may be interpreted as being due to eruptive LBV mass loss (e.g., Smith & Owocki 2006; Groh et al. 2020).

4.3. Suppressing density inversions

According to the MLT, convection sets in when the temperature gradient of the surrounding material is greater than the gradient interior to the moving element: ∇_T >  ∇_element. For convection to transport the maximum possible energy, the element should move adiabatically, that is without dissipating energy in the surroundings or ∇_T ∼ ∇_ad.

The difference between ∇_T and ∇_ad is termed as the superadiabaticity and is a measure of the efficiency of convection (Maeder 2009; Kippenhahn et al. 2012). For efficient convection, superadiabaticity is positive but close to 0. It means that the element loses hardly any energy as it traverses the mixing length. However, a higher value of superadiabaticity (≳10⁻²) implies that the element suffers energy losses as it travels, and by the time it reaches the end of the mixing length it is left with hardly any energy, thereby, rendering the convective transport of energy quite inefficient.

In the envelopes of massive stars superadiabaticity is the order of unity and the convective transport of energy given by MLT is inefficient, providing fertile ground for density inversions to form and be sustained. Some authors even consider density inversions to be nonphysical, possibly an artifact of the MLT and 1D stellar evolution (Ekström et al. 2012) while others consider the possible presence of some unknown mixing mechanism that helps stars get rid of density inversions in nature (Paxton et al. 2013). Either way, since these density inversions are associated with numerical instabilities in the models of massive stars, many stellar evolution codes suppress them either by limiting ∇_T or by reducing ∇_T to make it closer to ∇_ad. This reduces the superadiabaticity, thereby making convection efficient and helping stellar models overcome numerical instabilities.

MESA uses the latter approach of reducing ∇_T in the stellar envelope through a method known as MLT++². In this method, whenever the superadiabaticity exceeds a predefined threshold gradT_excess_f1, MESA decreases ∇_T to make it closer to ∇_ad. The amount of decrease is given by the combination of the parameter gradT_excess_f2 which can be defined by the user, and the parameter gradT_excess_alpha which is calculated based on the maximum of L_rad/L_Edd and the minimum of P_gas/P_total (see Appendix A for details). In general, a smaller gradT_excess_f2 implies a larger reduction in the superadiabaticity and more efficient convective transport of energy.

We find that using the default values of the MLT++ parameters completely suppresses density inversions but it also gives unrealistic values of luminosity for the most massive stars in our set. Therefore, we test the models in the standard set with different combinations of parameters in MLT++, as described in Appendix A. Compared to the default values of the MLT++ parameters, we find that using a smaller value of gradT_excess_f2=10⁻¹, and therefore a smaller reduction in the superadiabaticity but occurring more frequently inside the star (with λ₁ = 0.6 and β₁ = 0.05), is sufficient for the smooth evolution of the models without any numerical instabilities or inaccuracies.

The stellar models evolved using MLT++ are shown in the bottom panel of Fig. 4. The evolutionary paths of models computed with MLT++ are quite similar to models evolved with enhanced mixing length. Although, unlike the models with enhanced mixing length, models with MLT++ are not limited to log T_eff/K ≈ 3.73 and evolve to lower effective temperatures. The temperature profile of a 110 M_⊙ star evolved with MLT++ (Fig. 7) again shows a similar behavior compared to the temperature profile of the 110 M_⊙ star evolved with enhanced mixing length (Fig. 5), at similar log L and log T_eff. However, MLT++ artificially reduces ∇_T, such that the superadiabaticity and the specific entropy at the base of the convective envelope remain small despite having lower convective velocity compared to the 110 M_⊙ model with enhanced mixing length.

5. Impact on the stellar properties

Each of the three solutions described in Sect. 4 help compute the evolution of massive stellar models in the standard set until the end of carbon burning. However, in the process they also modify the evolutionary pathway of the stars and impact their evolutionary outcome. In this section, we compare the set of stellar models obtained with the minimum numerical enhancement from each solution and determine the impact of these solutions on the structure and evolution of the massive stellar models.

5.1. Structure of the star

Figure 8 shows the Kippenhahn plot of a 110 M_⊙ star–depicting regions within the star by mass as a function of time in the period leading up to the end of the run–for each of the models computed with the pragmatic solution described in Sect. 4. The evolution of the 110 M_⊙ model evolved using MLT++ is similar to the model with enhanced mixing length but quite different to the evolution of the model computed with enhanced mass loss.

Fig. 8. Kippenhahn plots showing the structure of a 110 M⊙ star evolved with enhanced mixing (left), with MLT++ (middle) and with enhanced mass loss (right). The Y-axis represents the mass coordinate inside the star while the X-axis represents the time remaining in the life of the star before the end of the run is reached. Green, purple and red hatching mark the regions with convection, overshooting and semiconvective mixing, respectively. In all three panels, the helium core boundary is the outermost location where the hydrogen mass fraction is < 0.01, while the helium mass fraction is ≥0.01, and is represented by the blue dashed line. Similarly, the carbon core boundary is defined as the outermost location where the hydrogen and helium mass fraction are < 0.01 while the carbon mass fraction is ≥0.01. It is represented as the red dashed line.

All three models start with a 90 M_⊙ convective core, (shown by the green hatching) accompanied by a 3 M_⊙ overshoot region outside the core, (shown by the purple hatching). The convective core decreases in size as the star evolves through the main sequence.

In models computed with enhanced mixing length and MLT++, thin strips of convection, two close to the surface and the third at about 90 M_⊙, are formed as the star encounters the hydrogen, helium and iron opacity bumps in the final 10⁶ yr of its evolution. The envelope inflation due to the iron opacity bump causes the star to become a red supergiant before core hydrogen burning can finish. The star suffers higher mass-loss rates as a red supergiant that can be – seen as a rapid decline in the total mass of the star (black solid line) – until a 60 M_⊙ helium core is formed – depicted by a blue dotted line in the figure. The star continues to lose mass, as the core helium burning and the hydrogen shell burning ensues, although at a comparatively lower rate. It ultimately loses its envelope to form a naked helium star in the final 10⁵ yr of its evolution, before finally forming a roughly 50 M_⊙ carbon core in the end.

In the model computed with enhanced mass loss, mass-loss rates are quite high (> 10⁻⁵ M_⊙ yr⁻¹) during the main-sequence evolution. The convective core shrinks rapidly in response to high mass-loss rates. The star loses its outer layers before density inversion due to hydrogen and helium opacity bumps can happen, becoming a naked helium star shortly after a 50 M_⊙ helium core is formed. The final product is a 40 M_⊙ naked helium star with a 35 M_⊙ carbon core.

The similarities in the evolution of the models using MLT++ and enhanced mixing can be understood as follows. In models with enhanced mixing, convection is already efficient in the core, owing to its high density and increasing α_MLT (and therefore the mixing length) hardly makes any difference. However, the density in the subsurface of layers of the star can be ≈10⁻¹⁰ g cm⁻³ and convection is highly inefficient, and therefore increasing α_MLT leads to more efficient convective transport of energy in the stellar envelope.

For models using MLT++ to suppress density inversions, the story is similar. Reducing the temperature gradient (superadiabaticity) prevents radiative losses from the convective cells, making them more efficient at transporting energy. However, near the center convection is nearly adiabatic and the value of superadiabaticity is small ( <  10⁻⁴), therefore MLT++ is not applicable there. Thus models with high α_MLT and MLT++ produce similar core structures.

For models with extra mass loss the evolution is quite different from the first two cases as the star loses quite a lot of mass even on the main sequence. This same trend continues for the post-main-sequence evolution. Therefore, it has the least of both helium and carbon core masses. Similar to the other two solutions, the 110 M_⊙ model again loses all its envelope and ends up as a naked helium star during the core helium burning phase.

Interestingly, the maximum mass-loss rate encountered by the 110 M_⊙ star with enhanced mass loss is lower than the models with enhanced mixing and with MLT++, as shown in Fig. 9. However, models with enhanced mass loss by construction have higher mass-loss rates during most of the main sequence, therefore they become a naked helium star without ever undergoing the red-supergiant phase where the peak in mass-loss rates usually occurs (due to the large radius of the star). Models with enhanced mixing and MLT++ lose their envelope later in the evolution as they encounter high mass-loss rates as a red supergiant and therefore end up with higher total mass compared to the model with enhanced mass loss.

Fig. 9. Variation in mass-loss rates with time for a 110 M⊙ star evolved with enhanced mixing length (blue line), enhanced mass loss (pink line) and with MLT++ (orange line). While the maximum mass-loss rate encountered by the 110 M⊙ model with enhanced mass loss is less than the maximum mass-loss rate encountered by the models evolved with the other two solutions, it still ends up being the least massive of all in the end. See Sect. 5.1 for an explanation.

5.2. Final mass and remnant properties

Massive stars are expected to end their lives in supernovae, leaving behind compact remnants (neutron stars and black holes) in a core-collapse or a pulsational-pair instability supernova or undergoing complete disruption in a pair-instability supernova. Their demise as a supernova explosion is also important for modifying the chemical and energy makeup of their surroundings, paving way for the formation of new generations of stars. Furthermore, stellar remnants are important for studies across a wide spectrum. They are the progenitors of X-ray binaries, gamma rays bursts, and compact binary mergers that lead to gravitational waves. Hence, it is important to quantify the differences between the properties of the remnants formed by the massive star models.

There are many prescriptions available in the literature that relate the final properties of the star with the mass of the remnant it would form (e.g., Eldridge & Tout 2004; O’Connor & Ott 2011; Fryer et al. 2012; Ertl et al. 2016). Here we use the Belczynski et al. (2008) prescription, which is the same as the StarTrack prescription in Fryer et al. (2012), to calculate the mass of the stellar remnants for each set of models. The method uses the total mass and the core mass of the star at the end of carbon burning to calculate the mass of the remnant (we refer the interested reader to Sect. 6.1 of Agrawal et al. 2020, for further details of the method).

Following Belczynski et al. (2010), we plot the remnant mass of the stars as a function of their initial mass as given by the three sets of stellar models computed using the pragmatic solutions in Fig. 10. The top panel of the figure shows the remnant mass of the stars while the bottom panel shows the final total mass of the star and the mass of the carbon-oxygen core.

Fig. 10. Initial versus final stellar masses for the three sets of stellar models that were computed using the pragmatic solutions. The top panel shows the mass of stellar remnants as a function of their initial mass, as predicted by the Belczynski et al. (2008) prescription for the sets of stellar models computed using the pragmatic solutions. The bottom panel shows the final total mass (solid line) and the carbon-oxygen core mass (dashed line) that were used in calculating the remnant mass. For stars more massive than 60 M⊙ the differences in the remnant mass can be up to 14 M⊙.

For all the sets of models, the Belczynski et al. (2008) prescription predicts the formation of black holes with masses in the range of 13–50 M_⊙. The difference in the mass of black holes as predicted by each set is less than 2 M_⊙ for stars with initial masses up to 60 M_⊙, with the exception of the 30 M_⊙ star where the model with enhanced mixing length predicts a significantly higher remnant mass (22 M_⊙ black hole) compared to other two sets (13 M_⊙ black hole). For stars more massive than 60 M_⊙ the curve diverges rapidly, and the difference between the remnant masses can be up to 14 M_⊙. A similar trend can be seen in the total mass and core masses for each set. For the 30 M_⊙ star, the model with enhanced mixing length predicts a higher final total mass but lower carbon-oxygen core mass compared to the other two sets. This is because a larger mixing length in the model leads to a more compact star with lesser mass loss. Thus, the 30 M_⊙ model is able to retain most of its envelope and ends up with the final total mass of 26 M_⊙. Similarly, the origin of differences in the black hole mass predictions can be traced back to the mass-loss rates experienced by each model which themselves are dependent on the surface properties of the star.

In models with enhanced mass loss, stars with initial masses greater than 50 M_⊙ lose their hydrogen envelopes due to mass-loss and become Wolf–Rayet stars (cf. Fig. 4). As mentioned in Sect. 3, we use the mass-loss prescription from Nugis & Lamers (2000) for Wolf–Rayet stars. Whilst more modern Wolf–Rayet mass-loss prescriptions exist (e.g., Gräfener & Hamann 2005; Vink & de Koter 2005), and may lead to different absolute mass-loss rates, we expect that the relative differences between the models shown in Fig. 10 should be robust to this difference (see for example Belczynski et al. 2010).

5.3. The maximum radial expansion

Figure 11 shows the maximum radial expansion achieved by the stars during their evolution, computed with the three different pragmatic solutions. The maximum difference in the maximum radial expansion is ∼2000 R_⊙ which occurs between models with enhanced mass loss and models with MLT++, for the most massive 110 M_⊙ star in the set. For models with enhanced mixing length and MLT++ which appear to undergo similar evolutionary paths and final fates, the difference in maximum radial expansion can still be up to 1000 R_⊙, especially for stars in the mass range 30–80 M_⊙. Even for the 110 M_⊙ star, where the difference between models with enhanced mixing length and MLT++ appears to be the least, the maximum radius can differ by 500 R_⊙. This has important implications for the binary evolution of the star, as radial proximity determines the episodes of mass transfer in close binary systems.

Fig. 11. Maximum radial expansion achieved by a star during its lifetime as a function of its initial mass, as predicted by the three sets of stellar models computed using the pragmatic solutions. The difference in predictions of the maximum radial expansion by each model varies between 500–2000 R⊙, and can have important implications if the star is in a binary system.

These differences in the stellar radii are again the result of the pragmatic solutions used in each set. For models with initial masses greater than 60 M_⊙ and computed with enhanced mass loss, high mass-loss rates strip the envelope of the stars before they can become a red supergiant. Thus these stars evolve directly toward the naked helium star phase and show the least radial expansion. For stars with enhanced mixing length, a higher mixing length means the fluid element is more efficient at transporting energy through convection. This decreases the size of the convective cells in the envelope and the model remains more compact and at higher effective temperatures than the models from the other two sets.

An interesting case is presented by models with MLT++ where reducing the temperature gradient can have the same effect as excess envelope mixing (Sabhahit et al. 2021). For stars up to 50 M_⊙, models using MLT++ closely mimic the behavior of models with enhanced mass loss, although beyond 50 M_⊙, the excess envelope mixing in models with MLT++ makes them resemble more closely the models with enhanced mixing.

6. Clues from observations: The Humphreys-Davidson limit

It is well-established that the evolution of massive stars is highly uncertain. It is riddled with many physical and numerical problems. Therefore the need for pragmatic solutions arises. To make matters worse, it is difficult to say which solution is closer to reality as there is an apparent scarcity of massive stars in the region of the HR diagram where numerical issues with density inversions occur.

First characterized by Humphreys & Davidson (1979), the Humphreys-Davidson (HD) limit defines the region in the HR diagram above log L/L_⊙ = 5.8 where very few massive stars in the Galaxy have been observed to date. Using recent observations of red supergiants in the Large Magellanic Cloud (LMC) and the Small Magellanic Cloud (SMC), Davies et al. (2018) found the limiting luminosity of the brightest red supergiants to be log L/L_⊙ = 5.5 in both galaxies, slightly lower than previous studies. The LMC and SMC are lower metallicity environments than the Milky Way, closer to the metallicity used in this work. The findings of Davies et al. (2018) also suggest that the maximum red-supergiant luminosity does not vary strongly with metallicity.

While the absence of massive stars as red supergiants beyond the HD limit may not help trace the exact evolutionary path of massive stars, it does provide an important clue that stars more massive than about 40 M_⊙ do not spend much time in the HD region. In recent years, several studies (e.g., Castro et al. 2018; Kaiser et al. 2020; Vink et al. 2021) have tried to constrain the different mixing mechanisms (such as semiconvection and overshooting), mass-loss rates, and binary properties of massive stars by using stellar models evolved with different physical inputs to reproduce the HD limit. However, they also employ pragmatic solutions for instabilities due to density inversions to compute the complete evolution of the stellar track. Using these pragmatic solutions interferes with other physical inputs and adds an implicit bias in the computation of models. As these pragmatic solutions can be different across different stellar evolution codes, results obtained with them might not reflect the true value of the physical parameter being constrained for the massive stars.

Studies that do not employ any pragmatic solutions (e.g., Klencki et al. 2018) are usually limited to stars less massive than 30 M_⊙ or to the evolutionary phases before numerical instabilities arise in more massive stars. While this helps avoid bias due to pragmatic solutions, it also inhibits exploring the late-stage evolution of massive stars, such as the end of core helium burning and beyond. Thereby, affecting the potential studies involving stellar remnants and transients.

Some studies also show that the HD limit can also be reproduced by stellar models by just using these pragmatic solutions. For example, using just MLT++ as the source of excess envelope mixing in massive stars up to 50 M_⊙, Sabhahit et al. (2021) were able to reproduce the lack of massive stellar models beyond the HD limit and the Davies limit at Galactic (Z = 0.017), LMC (Z = 0.008) and SMC (Z = 0.004) metallicity.

Recently, Gilkis et al. (2021) showed that using significantly enhanced mixing parameters in stellar models can reduce the time spent by stars beyond the HD limit. Following Gilkis et al. (2021), in Fig. 12 we show the amount of time stars spend beyond the HD limit as a function of initial mass for each of our sets of models. We see that the models computed with enhanced mass loss spend the least time (except for a 40 M_⊙ stellar model) while most of the massive stars computed with enhanced mixing and MLT++ can spend between 2 × 10⁵ and 4 × 10⁵ yr in the HD region. Therefore, models with enhanced mass loss may appear to be closest to observations (or lack of observations) of massive stars. This, however, has serious implications for gravitational wave observations, as the maximum black hole mass predicted by this set of the model is just 35 M_⊙ (see Fig. 10).

Fig. 12. Time spent beyond the HD limit as a function of the initial mass of the star as given by each set of stellar models computed using the pragmatic solutions.

To unravel this problem, we plot the evolutionary tracks, colored according to the mass-loss rates, from the set computed with the enhanced mass loss in Fig. 13. Despite the enhancement near the Eddington limit, the maximum mass-loss rates for models with enhanced mass loss seem to be consistent with the typical mass-loss rates for massive stars (Smith 2014). The maximum mass-loss rate of 1.8 × 10⁻³ M_⊙ yr⁻¹ is experienced by stars in the 40–60 M_⊙ mass range during the red-supergiant phase of their evolution. However, these rates only last for a maximum of 7 × 10³ yr. For stars more massive than 60 M_⊙, the maximum mass-loss rates are an order of magnitude lower, about 3 × 10⁻⁴ M_⊙ yr⁻¹, and last between 3 × 10²–2 × 10⁴ yr. The lower mass-loss rates encountered by more massive stars are a consequence of the enhanced mass loss during their main-sequence evolution. These stars experience mass-loss rates of about 10⁻⁵ M_⊙ yr⁻¹ for more than 3 × 10⁵ yr. Therefore, they lose a significant portion of their envelope while on the main sequence and do not undergo the red-supergiant phase of evolution.

Fig. 13. HR diagram showing the mass-loss rates for stars in the mass range 30–110 M⊙, computed with the enhanced mass loss described in Sect. 4.2.

Whether massive stars can retain these high mass-loss rates for prolonged periods of time is currently questionable. Recent (theoretical and observed) estimates of steady-state mass-loss from massive O stars (M  ≳  50 M_⊙) suggest that the mass-loss rates for these stars should be lower, rather than higher, than those typically used in computing single star models (Beasor et al. 2020; Vink & Sander 2021; Björklund et al. 2021; Hawcroft et al. 2021). However, these massive stars also experience eruptive mass-loss during LBV phases, which can be important for the mass-loss history of these stars (Smith & Owocki 2006), with effective mass-loss rates up to 10⁻⁴ M_⊙ yr⁻¹. Such mass-loss is often neglected in one-dimensional stellar models, but can be approximately included as an average, enhanced mass-loss rate, as we do here (e.g., Groh et al. 2020; Sabhahit et al. 2022). Interaction with a companion in a binary system can also lead to higher mass-loss rates, but again the smaller radii predicted by these models (Fig. 11) make binary interactions less likely. Therefore, the validity of mass-loss rates in models computed with extra mass-loss, or in general any of the pragmatic solutions used by the codes cannot be ascertained at present.

Mixing processes such as rotation, semiconvection, and overshooting have also been shown to have a large impact on the red-supergiant phase of evolution for massive stars and on the reproducibility of the HD limit (Schootemeijer et al. 2019; Higgins & Vink 2020; Gilkis et al. 2021; Sabhahit et al. 2021). Further, many studies suggest that magnetic fields can have a significant effect on the subsurface convection regions of the stars and hence on the density inversions. For example, Jiang et al. (2018) have shown that the turbulent velocity fields around the iron opacity peak can be escalated in the presence of magnetic fields. We have not explored the role of semiconvection, rotation, magnetic fields and binarity in this work, each of which we acknowledge can have a significant impact on the evolution of massive stars. Evolving massive stellar models with these properties is important and will be a part of future work.

7. Conclusion

In this work, we performed a systematic study of the density inversions arising due to the super-Eddington layers in the envelope of the massive stars and their impact on computing the evolution of massive stars with the 1D stellar evolution code MESA. Using commonly used input parameters for massive stellar models, we computed the nonrotating evolutionary tracks of 10–110 M_⊙ stars for Z = 0.00142. Models with initial masses between 30 and 110 M_⊙ in this standard set fail to reach the end of core carbon depletion (X_c ≤ 10⁻²) as the time steps become too small (on the order of days) as models encounter density inversions in the envelope. We find that this is due to large specific entropy and a low gas pressure fraction at the iron opacity bump and a high Eddington factor at the hydrogen and helium opacity bumps in the stellar envelope.

We recomputed the models for 30–110 M_⊙ stars using three pragmatic solutions: enhancing the mixing length parameter, enhancing mass-loss rates, and suppressing density inversions by adopting MLT++. These three solutions are a general form of those used by different stellar evolution codes to resolve numerical instabilities associated with density inversions in models of massive stars. With each solution, we were able to compute the evolution of massive stars up to carbon depletion in the core.

To determine the least possible impact each solution can have on the evolution of the star, we picked the sets with the minimum enhancement for each solution. With these, we find that density inversion regions can still form in these models, but the evolution of models proceeds smoothly, that is to say without time steps becoming too short.

We compare the stellar models evolved using each solution. Even with sets of models with minimum numerical enhancement from each solution, the remnant mass of the stars can vary by up to 14 M_⊙, while the maximum radial expansion achieved by stars can vary by up to 2000 R_⊙ between the sets. These differences are important for comparing stellar models with observations and the possible feedback on the evolution of massive stars. For example, the abovementioned differences can have huge implications for studies involving binary interactions of stars and stellar remnants.

We also find that the differences in the various evolutionary properties predicted by the set of models computed with enhanced mixing length and those using MLT++ are quite small. However, these models show a difference of 1000 K in the minimum effective temperatures achieved by the star as a red supergiant (cf. Sect. 4.3). The effective temperature is critical for determining the spectral properties and surface abundances of the stars (Davies et al. 2013). Thus, this discrepancy in the effective temperature can be significant for studies involving the abundance properties of massive stars, such as galactic chemical composition studies.

In commonly used 1D stellar evolution models, a combination of these solutions is used and not necessarily limited to their minimum value. Thus, the differences in the evolutionary properties of massive stars can he higher than what we get here, as shown in Paper I.

Here we have considered an intermediate value of metallicity (Z= 1.42 × 10⁻³), as this is the metallicity where progenitors of current gravitational wave observations are more likely to form. Due to the dependence of opacity on chemical composition and the diminishing effect of opacity peaks at lower metallicities (Sanyal et al. 2017), numerical issues related to the proximity to the Eddington limit become less prominent at lower metallicities. For example, using the same set of input parameters as the standard set (Sect. 3), we find that stellar models fail to evolve for initial masses beyond 20 M_⊙ at higher (solar) metallicity (Z = 1.42 × 10⁻²), but they evolve smoothly for initial masses up to 100 M_⊙ at lower metallicity (Z = 1.42 × 10⁻⁴).

Multidimensional stellar models suggest that convection can be more turbulent and nonlocalized than assumed by MLT (Kupka 2009). Three-dimensional modeling of the density inversion regions in massive stars by Jiang et al. (2015, 2018) shows a complex interplay of convective and radiative transport dependent on the ratio of the photon diffusion time to the dynamical time and a smaller convective flux compared to 1D codes. There are ongoing efforts to improve the treatment of convective mixing in 1D codes using results from 3D simulations (Mosumgaard et al. 2018; Arnett et al. 2019; Schultz et al. 2020).

As observations of massive stars become more accessible, and as our ability to compute 3D models of massive stars improves, the problem of density inversions might be resolved in the future. Meanwhile, it is important to be aware of and acknowledge the impact different pragmatic solutions can have on the evolutionary model sequences of massive stars.

Acknowledgments

We thank Ross Church, Jan Eldridge, Debashis Sanyal, Matteo Cantiello, and Mathieu Renzo for useful comments and discussions. We also thank Pablo Marchant for the script (https://github.com/orlox/mkipp) to make Kippenhahn diagrams. P.A., J.H. and S.S. acknowledge the support by the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), through project number CE170100004. SS is a recipient of an ARC Discovery Early Career Research Award (DE220100241). D.Sz. has been supported by the Alexander von Humboldt Foundation. This work made use of the OzSTAR high performance computer at Swinburne University of Technology. OzSTAR is funded by Swinburne University of Technology and the National Collaborative Research Infrastructure Strategy (NCRIS). This research was funded in part by the National Science Center (NCN), Poland under grant number OPUS 2021/41/B/ST9/00757. For the purpose of Open Access, the author has applied a CC-BY public copyright license to any Author Accepted Manuscript (AAM) version arising from this submission.

References

Appendix A: Details of MLT++

In MLT++, MESA makes convection efficient by reducing superadiabaticity in the stellar envelope. In the regions where superadiabaticity is greater than threshold gradT_excess_f1, it decreases ∇_T to make it closer to ∇_ad. The fraction of decrease is determined by the parameter gradT_excess_alpha or α and is calculated based on the value of (cf. equation 38 of Paxton et al. 2013):

$$\begin{array}{c}\hfill {\lambda }_{max}\equiv max\left(\frac{L_{\mathrm{rad}}}{L_{\mathrm{Edd}}}\right)\text{and}{\beta }_{min}\equiv min\left(\frac{P_{\mathrm{gas}}}{P}\right).\end{array}$$(A.1)

For each stellar model, MESA computes α by comparing λ_max with the thresholds λ₁ and λ₂ and β_min with β₁ and β₂ using the following conditions.

If λ_max ≥ λ₁, then

$$\begin{array}{c}\hfill \alpha =\{\begin{array}{cc}1\hfill & {\beta }_{\min }\le {\beta }_1\hfill \\ \frac{{\beta }_1+d\beta -{\beta }_{\min }}{d\beta }\hfill & {\beta }_1<{\beta }_{\min }<{\beta }_1+d\beta \hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}.\end{array}$$(A.2)

If λ_max ≥ λ₂, then

$$\begin{array}{c}\hfill \alpha =\{\begin{array}{cc}1\hfill & {\beta }_{\min }\le {\beta }_{\mathrm{limit}}\hfill \\ \frac{{\beta }_{\mathrm{limit}}+d\beta -{\beta }_{\min }}{d\beta }\hfill & {\beta }_{\min }<{\beta }_{\mathrm{limit}}+d\beta \hfill \\ 0\hfill & {\beta }_{\min }\ge {\beta }_{\mathrm{limit}}+d\beta \hfill \end{array}.\end{array}$$(A.3)

If λ_max >  λ₂ − dλ, then

$$\begin{array}{c}\hfill \alpha =\{\begin{array}{cc}1\hfill & {\beta }_{\min }\le {\beta }_2\hfill \\ \frac{{\lambda }_{\max }+d\lambda -{\lambda }_2}{d\lambda }\hfill & {\beta }_2<{\beta }_{\min }<{\beta }_2+d\beta \hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}.\end{array}$$(A.4)

The net fraction of decrease is determined by a combination of alpha and user defined gradT_excess_f2 or f₂ using the following equation,

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{net}}=f_2+(1-f_2)(1-\alpha ).\end{array}$$(A.5)

The excess fraction is then subtracted from ∇_T, to give reduced ∇_T, new as,

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{\mathrm{T},\mathrm{new}}={\alpha }_{\mathrm{net}}\times {\mathrm{\nabla }}_{\mathrm{T}}+(1-{\alpha }_{\mathrm{net}})\times {\mathrm{\nabla }}_{\mathrm{ad}}.\end{array}$$(A.6)

The default values for the different thresholds in the Equations A.2–A.4 are: λ₁ = 1.0 and λ₂ = 0.5, β₁ = 0.35 and β₂ = 0.25, dλ = 0.1 and dβ = 0.1. gradT_excess_f1 defaults to 10⁻⁴ while the default value of gradT_excess_f2 or f₂ is 10⁻³. These values can be redefined by the user.

Equations A.2–A.4 yield a value of α between 0 and 1. The maximum fraction of ∇_T used in calculating ∇_T, new is limited to f₂ (when α = 1), with smaller f₂ implying larger contribution of ∇_ad in the equation A.6 and therefore larger reduction in superadiabaticity.

We tested our models with for different values of f₂, λ₁ and β₁. The HR diagram for 110 M_⊙ stellar model, calculated with the four different combinations of parameters in MLT++ is shown in Figure A.1. We find that using a small reduction in superadiabaticity with f₂=10⁻¹ and setting λ₁ = 0.6 and β₁ = 0.05 helps the stellar models reach completion timely and without any numerical inaccuracies or difficulty.

Fig. A.1. Evolutionary tracks of a 110 M⊙ star evolved with four different MLT++ parameter combinations. Apart from the changes indicated in the figure legend, default values for all other MLT++ parameters have been used. A negative value of λ1 implies that the maximum reduction in superadiabaticity is applied throughout the track. The tracks corresponding to the default values of MLT++ and to λ1 = −1.0 are indistinguishable here and predict unrealistic luminosities.

All Figures

Fig. 1. Hertzsprung-Russell (HR) diagram showing stellar models from the standard set colored by the maximum of Lrad/LEdd (left panel) and the minimum of Pgas/Ptotal (right panel). The blue cross marks the end of core hydrogen burning and the red cross marks the end of core helium burning (where applicable). The brown dashed line in the left panel denotes the position of the observational Humphreys & Davidson (1979) limit beyond which few stars are observed, while the gray dashed line signifies the luminosities of the brightest red supergiants as inferred by Davies et al. (2018). Large values of Lrad/LEdd and small values of Pgas/Ptotal in the envelopes of stars with initial masses 30 M⊙ and above causes the evolution of these stars to become halted at log Teff/K ≈ 3.7, and their models fail to even finish core helium burning.

In the text

Fig. 2. Temperature profile of a 110 M⊙ star at the ZAMS evolved using the standard set. Opacity peaks due to partial ionization states of iron are present at 1.5 × 106 and 1.8 × 105 K. However, Lrad/LEdd is less than 1 throughout the star and the evolution of the star proceeds smoothly. For clarity, only the outermost layers of the star (with T ≤ 5 × 106 K) containing the various opacity peaks have been plotted. See Sect. 3.2 for details of each panel.

In the text

Fig. 3. Temperature profile of the outermost layers of a 110 M⊙ star at the end of the simulation, evolved using the standard set. The high luminosity of the star combined with the peak in opacity due to hydrogen and helium ionization around 104 K causes Lrad/LEdd to exceed unity. This causes density and gas pressure inversions at 3 × 104 K. However, the low density of the environment renders convection inefficient and the star struggles to evolve despite reaching supersonic convective velocities.

In the text

Fig. 4. HR diagrams showing the evolutionary tracks graded with the maximum of Lrad/LEdd for stars in the mass range 30–110 M⊙, computed with the model variations described in Sect. 4. The top panel shows the tracks computed with enhanced mixing (αMLT = 5.46), the middle panel shows tracks with enhanced mass loss (ηDutch = 8.0 whenever L exceeds LEdd) and the bottom panel shows tracks computed with MLT++. Similar to Fig. 1, blue and red crosses mark the end of core hydrogen burning and core helium burning while the brown and gray dashed lines signify the Humphreys & Davidson (1979) limit and the Davies et al. (2018) limit respectively.

In the text

Fig. 5. Temperature profile of the outermost layers of a 110 M⊙ star computed with the mixing length parameter, αMLT = 5.46. Similar to Fig. 3, a density inversion can be seen in panel (c). However, higher convective velocity (panel f) and smaller superadiabaticity (panel g) resulting from the higher value of αMLT helps to keep the specific entropy small (panel d) and time steps large enough to evolve the model without numerical instabilities.

In the text

Fig. 6. Temperature profile of the outermost layers of a 70 M⊙ star computed with enhanced mass loss, as described in Sect. 4.2. Opacity peaks due to the ionization states of hydrogen, helium and iron and associated density inversions can be seen in panel (b) and panel (c) respectively. However, high mass-loss rates remove the outer layers containing the hydrogen and helium opacity peaks, moderating the specific entropy at the base of the convective envelope and the model is able to evolve until completion with reasonably large time steps.

In the text

Fig. 7. Temperature profile of the outermost layers of a 110 M⊙ star computed with the MLT++ method of MESA. The method artificially reduces the difference between ∇T and ∇ad in the stellar envelope, acting as an source of additional envelope mixing (see Sect. 4.3 for details). The small envelope mixing used here is not enough to suppress density inversions completely (panel c) but it does help to keep the specific entropy small (panel d), gas pressure fraction non-negligible (panel e) and therefore time steps reasonable.

In the text

Fig. 8. Kippenhahn plots showing the structure of a 110 M⊙ star evolved with enhanced mixing (left), with MLT++ (middle) and with enhanced mass loss (right). The Y-axis represents the mass coordinate inside the star while the X-axis represents the time remaining in the life of the star before the end of the run is reached. Green, purple and red hatching mark the regions with convection, overshooting and semiconvective mixing, respectively. In all three panels, the helium core boundary is the outermost location where the hydrogen mass fraction is < 0.01, while the helium mass fraction is ≥0.01, and is represented by the blue dashed line. Similarly, the carbon core boundary is defined as the outermost location where the hydrogen and helium mass fraction are < 0.01 while the carbon mass fraction is ≥0.01. It is represented as the red dashed line.

In the text

Fig. 9. Variation in mass-loss rates with time for a 110 M⊙ star evolved with enhanced mixing length (blue line), enhanced mass loss (pink line) and with MLT++ (orange line). While the maximum mass-loss rate encountered by the 110 M⊙ model with enhanced mass loss is less than the maximum mass-loss rate encountered by the models evolved with the other two solutions, it still ends up being the least massive of all in the end. See Sect. 5.1 for an explanation.

In the text

Fig. 10. Initial versus final stellar masses for the three sets of stellar models that were computed using the pragmatic solutions. The top panel shows the mass of stellar remnants as a function of their initial mass, as predicted by the Belczynski et al. (2008) prescription for the sets of stellar models computed using the pragmatic solutions. The bottom panel shows the final total mass (solid line) and the carbon-oxygen core mass (dashed line) that were used in calculating the remnant mass. For stars more massive than 60 M⊙ the differences in the remnant mass can be up to 14 M⊙.

In the text

	Fig. 11. Maximum radial expansion achieved by a star during its lifetime as a function of its initial mass, as predicted by the three sets of stellar models computed using the pragmatic solutions. The difference in predictions of the maximum radial expansion by each model varies between 500–2000 R⊙, and can have important implications if the star is in a binary system.
In the text

	Fig. 12. Time spent beyond the HD limit as a function of the initial mass of the star as given by each set of stellar models computed using the pragmatic solutions.
In the text

	Fig. 13. HR diagram showing the mass-loss rates for stars in the mass range 30–110 M⊙, computed with the enhanced mass loss described in Sect. 4.2.
In the text

Fig. A.1. Evolutionary tracks of a 110 M⊙ star evolved with four different MLT++ parameter combinations. Apart from the changes indicated in the figure legend, default values for all other MLT++ parameters have been used. A negative value of λ1 implies that the maximum reduction in superadiabaticity is applied throughout the track. The tracks corresponding to the default values of MLT++ and to λ1 = −1.0 are indistinguishable here and predict unrealistic luminosities.

In the text