Eruptive mass loss in massive stars is known to occur but the mechanisms are not yet well understood. One proposed physical explanation appeals to opacitydriven superEddington luminosities in stellar envelopes. Here we present a 1D model for eruptive mass loss and implement this model in the MESA stellar evolution code. The model identifies regions in the star where the energy associated with a local superEddington luminosity exceeds the binding energy of the overlaying envelope. The material above such regions is ejected from the star. Stars with initial masses of 10100 M SUBSUB at solar and SMC metallicities are evolved through core helium burning with and without this new eruptive massloss scheme. We find that eruptive mass loss of up to 10SUP2SUP M SUBSUByrSUP1SUP can be driven by this mechanism and occurs in a vertical band on the HR diagram between inlineformula mmlmath overflowscrollmmlmn3.5mmlmnmmlmommlmommlmilogmmlmimmlmo stretchyfalsemmlmommlmsubmmlmrowmmlmiTmmlmimmlmrowmmlmrowmmlmieffmmlmimmlmrowmmlmsubmmlmrowmmlmo stretchytruemmlmommlmrowmmlmi mathvariantnormalKmmlmimmlmo stretchyfalsemmlmommlmommlmommlmn4.0mmlmnmmlmath inlineformula. This predicted eruptive mass loss prevents stars of initial masses 20 M SUBSUB from evolving to become red supergiants RSGs with the stars instead ending their lives as blue supergiants and offers a possible explanation for the observed lack of RSGs in that mass regime.
['Stellar physics', 'Stellar winds', 'Stellar evolution', '1621', '1636', '1599', 'Astrophysics - Solar and Stellar Astrophysics', 'Astrophysics - Astrophysics of Galaxies']
A Model for Eruptive Mass Loss in Massive Stars
{'A Model for Eruptive Mass Loss in Massive Stars': '<!-- image -->  \nShelley J. Cheng , 1, 2 Jared A. Goldberg , 2 Matteo Cantiello , 2, 3 Evan B. Bauer , 1 Mathieu Renzo , 4, 2 and Charlie Conroy 1  \n1 Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA 2 Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA 3 Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA 4 Steward Observatory, University of Arizona, 933 N Cherry Avenue, Tucson, AZ 85721, USA', 'ABSTRACT': 'Eruptive mass loss in massive stars is known to occur, but the mechanism(s) are not yet wellunderstood. One proposed physical explanation appeals to opacity-driven super-Eddington luminosities in stellar envelopes. Here, we present a 1D model for eruptive mass loss and implement this model in the MESA stellar evolution code. The model identifies regions in the star where the energy associated with a locally super-Eddington luminosity exceeds the binding energy of the overlaying envelope. The material above such regions is ejected from the star. Stars with initial masses 10 -100 M ⊙ at solar and SMC metallicities are evolved through core helium burning, with and without this new eruptive mass-loss scheme. We find that eruptive mass loss of up to ∼ 10 -2 M ⊙ yr -1 can be driven by this mechanism, and occurs in a vertical band on the HR diagram between 3 . 5 ≲ log( T eff / K) ≲ 4 . 0. This predicted eruptive mass loss prevents stars of initial masses ≳ 20 M ⊙ from evolving to become red supergiants, with the stars instead ending their lives as blue supergiants, and offers a possible explanation for the observed lack of red supergiants in that mass regime.', '1. INTRODUCTION': "Massive stars ( ≳ 10 M ⊙ ) are important for understanding the evolution and structure of many astrophysical systems. When studying galaxies, massive stars are important since they impact the overall density and movement of matter which alters galaxy structure, dynamics and evolution (e.g., see review by Tacconi et al. 2020). Additionally, mass loss from massive stars and their explosions play an important role in galaxy feedback (e.g., Heckman et al. 1990; Ceverino & Klypin 2009; Geen et al. 2015; Steinwandel & Goldberg 2023). At smaller scales, quantifying the intrinsic mass loss from massive stars can help constrain the masses of compact objects (Belczynski et al. 2010; Dominik et al. 2012; Renzo et al. 2017; Banerjee et al. 2020; Moriya & Yoon 2022; Fragos et al. 2023), study gravitational wave sources (Ziosi et al. 2014; Mapelli 2016; Uchida et al. 2019; Vink et al. 2021), and better understand the circumstellar environment responsible for observed signatures of interaction in early-time SN emission (e.g.,  \nshelley.cheng@cfa.harvard.edu  \nSchlegel 1990; Pastorello et al. 2007; Morozova et al. 2017; Bruch et al. 2023; Jacobson-Gal'an et al. 2024).  \nObservationally, massive stars are known to exhibit high levels of intrinsic mass loss (e.g., Wood et al. 1983, 1992; Smith 2014; Meynet & Maeder 2003). In single massive stars, mass loss can generally occur in three ways; a radiation-driven wind, a pulsation-induced dustdriven wind, and through eruptions. One of the first analytical radiation-driven wind models was the CAK model (Castor et al. 1975), which predicted a mass loss rate on the order of 10 -6 M ⊙ yr -1 for a O5-type main sequence star. The CAK model describes how the forces due to metal line absorption of photons can produce a mass outflow, and predicts mass loss rates and wind terminal velocities (see, e.g., Lamers & Cassinelli 1999; Owocki et al. 1988; Smith 2014; Cur'e & Araya 2023). While the mass loss rate predictions of CAK were in general agreement with observations (de Jager et al. 1988; Nieuwenhuijzen & de Jager 1990; Nugis & Lamers 2000), the emissivity of the recombination processes that form emission lines depends on the square of density, so small density inhomogeneities or 'clumps' in the wind can cause an over-prediction of mass loss rates if they are not properly modeled. Subsequent work by Crowther et al. (2002); Figer et al. (2002); Hillier et al. (2003); Massa  \net al. (2003); Evans et al. (2004); Puls et al. (2006); Beasor & Davies (2018); Beasor et al. (2020) have determined that CAK-based mass loss rates should be reduced by a factor of 2 -10 (Crowther et al. 2002; Hillier et al. 2003; Fullerton et al. 2006; Puls et al. 2008; Smith 2014).  \nPulsation-induced dust-driven winds are likely associated with red supergiants (RSGs) and AGB stars, where atmospheric shock waves caused by stellar pulsations and possibly convection - are thought to lift gas to low enough temperatures to form dust which then collisionally couples with photons and accelerates outwards as a wind (e.g., Liljegren et al. 2016; Willson 2000; van Loon et al. 2005; Goldberg et al. 2022; Chiavassa et al. 2024). Observations have suggested wind mass loss rates of up to 10 -3 M ⊙ yr -1 in RSGs (van Loon et al. 2005), but mass loss rates for RSGs are highly uncertain (Smith 2014) and potentially overestimated (e.g., Beasor et al. 2020, 2023; Decin et al. 2024). Many common prescriptions assume dust-enshrouded RSGs with higher mass loss rates (e.g., Georgy et al. 2012; Ekstrom et al. 2012; Goldman et al. 2017), and the choice of mass loss model can significantly alter the evolutionary track of a RSG, with lower mass loss rates allowing the star to remain in the red but higher mass loss rates driving the star towards the blue (Smith 2014; Renzo et al. 2017; Beasor & Davies 2018; Beasor et al. 2020; Josiek et al. 2024) after spending time in the RSG phase. Observationally, RSGs exhibit an (apparent) maximum luminosity of log( L/L ⊙ ) ≈ 5 . 5, which corresponds to a maximum initial mass of around 20 -30 M ⊙ (Massey et al. 2000; Heger et al. 2003; Levesque et al. 2009; Ekstrom et al. 2012; Dorn-Wallenstein et al. 2023). Post main-sequence mass loss has been suggested as a possible theory explaining this, particularly involving stars that exceed the Eddington luminosity limit (e.g., Puls et al. 2008). However, the specific nature of this mass loss is not well-understood, and the physical mechanism causing the dearth of RSGs above 20 -30 M ⊙ remains an open problem. Alternative explanations other than eruptive mass loss include envelope inflation (e.g., Sanyal et al. 2015, 2017) and binary interactions (e.g., Smith 2014; Morozova et al. 2015; Margutti et al. 2017).  \nAt even higher masses, luminous blue variables (LBVs) have shown very high mass loss rates of up to 0 . 1 -5 M ⊙ yr -1 (Humphreys & Davidson 1984; Smith 2017). This high rate of mass loss seen in LBVs is observed to only occur for a short period of time, and as such is referred to as 'eruptive massive loss.' Despite a few pioneering theoretical works (e.g., Owocki et al. 2004; Smith & Owocki 2006), efforts to form a comprehensive model for mass loss of massive stars have tra-  \ntionally focused on developing scaling formulae based on empirical rates (e.g., de Jager et al. 1988; Nieuwenhuijzen & de Jager 1990; Vink et al. 2001; van Loon et al. 2005). Additionally, Mauron & Josselin (2011) and Renzo et al. (2017) compare these scaling formulae for solar-metallicity non-rotating stars between 15 -35 M ⊙ . However, these steady wind mass loss schemes do not attempt to and cannot reproduce the high rates of observed eruptive massive loss. Therefore, there has been some work that scales and tunes the theoretical wind mass-loss schemes using observations to reproduce observed rates of mass loss (e.g., Ekstrom et al. 2012; Dessart et al. 2013; Meynet et al. 2015; Vink & Sabhahit 2023). However, these schemes do not provide a physically motivated model that self-consistently captures the effect of eruptive mass loss.  \nIn this work, we begin with an overview of existing literature on radiation-pressure-dominated envelopes of massive stars from one-dimensional (1D) and threedimensional (3D) models that motivates our work (Section 2). We then develop an energy-based model for eruptive mass loss and implement it within the 1D stellar evolution software MESA (Section 3), and apply our model to massive stars of differing masses (between 10 -100 M ⊙ ) and metallicities (at 0 . 2 Z ⊙ and 1 Z ⊙ ). We explore the predicted rate of eruptive mass loss (Section 4), and finally discuss our interpretation of our results and conclude in Section 5.", '2. RADIATION-PRESSURE-DOMINATED ENVELOPES OF MASSIVE STARS': "Energy is transported in stars primarily through radiative diffusion and convection, including in the radiation-pressure-dominated extended outer regions of massive star envelopes, where both mechanisms are present. Due to opacity peaks in the envelope (from recombination of ions of iron, helium, and hydrogen), the local radiative luminosity may exceed the Eddington limit in layers of the envelope below sub-Eddington regions (Jiang et al. 2015, 2018). Therefore, in these envelopes, density and gas pressure inversions may occur which are thought to contribute to outburst-like 'eruptive mass loss' (Joss et al. 1973; Humphreys & Davidson 1994; Paxton et al. 2013). Specifically, the gas pressure increases outwards when radiative luminosity is superEddington (when L rad > L Edd ), and under specific conditions, a density inversion may occur even if gas pressure is not inverted; see Joss et al. (1973) and Paxton et al. (2013) for more details.  \nThe envelopes of massive stars therefore present both physical and numerical challenges due to the presence of inversions and super-Eddington radiative luminos-  \nity. 1D numerical treatments can fail to find a hydrostatic equilibrium solution in these challenging regions (e.g., Paxton et al. 2013). This is the case when using mixing-length-theory (MLT, e.g., Bohm-Vitense 1958; Cox & Giuli 1968) or even time-dependent 1D formulations (e.g., Kuhfuss 1986; Jermyn et al. 2023). Realistic 3D simulations show the limits of the 1D treatment, stressing the importance of the interplay between radiative transport and the complex landscape of turbulent fluctuations in the stellar envelope (Jiang et al. 2015). However these calculations are computationally very expensive, which has prevented a detailed study of the full parameter space relevant to massive stars.  \nIn the stellar evolution code MESA used in this work, several modifications to MLT aimed at these regions in the envelopes of massive stars are available, namely MLT++ (Paxton et al. 2013) and an implicit method ( use superad reduction ; Jermyn et al. 2023). These methods involve artificially enhancing the efficiency of convection by reducing the superadiabaticity ( ∇-∇ ad ) in regions nearing the Eddington limit, which usually helps the code find a hydrostatic equilibrium solution. The first of these methods, MLT++, was introduced in Paxton et al. (2013). While MLT++ allowed evolution of massive stars to core-collapse, it lead to different efficiencies of energy transport in radiation-dominated stellar envelopes and the resulting evolution differs to that of convective energy transport by MLT. Subsequently, an implicit method was introduced in Jermyn et al. (2023) that produces more local and robust results in closer agreement to that from MLT compared to MLT++ (when comparing evolutionary tracks on the HR diagram from main sequence to RSG phase, see Figure 17 in Section 7.2 of Jermyn et al. 2023). Nevertheless, while these 1D recipes can help improve numerical convergence and evolve massive stars further, they are not directly physically motivated.  \nMass loss in massive stars has also been studied using 3D radiation-hydrodynamics (RHD) models. Indeed, 3D RHD models find enhanced convection in these radiation- and turbulent-pressure-dominated envelopes (e.g., Jiang et al. 2015, 2018; Schultz et al. 2023). Jiang et al. (2015) and Jiang et al. (2018) simulated the envelope of massive stars at specific points in their evolution, and determined that non-steady mass loss outbursts potentially consistent with LBVs can be driven by helium and iron opacity peaks. Their predicted mass loss rates are significantly lower than those inferred from observations (Humphreys & Davidson 1984; Smith 2017), possibly due to the fact that line-driving and dust formation are not included in the models. The importance of helium and iron opacity peaks in mediating continuum-  \nopacity-driven mass-loss is not surprising, as the dependence of opacity on temperature and density near opacity peaks can lead to radiation hydrodynamic instabilities (Prendergast & Spiegel 1973; Kiriakidis et al. 1993; Blaes & Socrates 2003; Su'arez-Madrigal et al. 2013).  \nAs photons carrying energy diffuse outwards through the envelopes, they may encounter regions where photon diffusion occurs on a faster timescale than dynamical, convective, and acoustic processes (radiatively inefficient convection). This occurs at optical depths lower than the critical value τ crit ≡ c/c s (e.g, Jiang et al. 2015), where c is the speed of light and c s is the sound speed. This can also be derived from balancing the diffusive radiation flux and the maximum convective flux in the radiation pressure-dominated regime (for more details see Jiang et al. 2015). This definition of critical optical depth is equivalent to the more general definition presented by Jermyn et al. (2022) (see also discussions in Sections 2 and 4.4 of Goldberg et al. 2022).  \nThe role of τ crit in determining convective regions where radiation nonetheless carries significant flux (often termed radiatively inefficient, or sometimes just 'inefficient' convection) has been demonstrated in 3D simulations. In particular, Jiang et al. (2015) simulated the radiation-dominated envelopes of massive stars in a plane parallel geometry and found that in certain regions close to the surface of the star more shallow than a critical optical depth τ < τ crit , convection is unable to carry the super-Eddington flux at the iron opacity peak (which is convectively unstable and locally superEddington). When τ ≫ τ crit , convection is efficient (in the radiative sense), with the averaged radiation advection flux consistent with the calculation of convection flux following MLT and sub-Eddington radiative acceleration. In contrast, radiation carries non-negligible flux for τ < τ crit , and a very significant fraction of the flux for τ ≪ τ crit , despite the convective motion. Even while τ ≪ τ crit , the envelope continues to behave in a turbulent and non-stationary manner, with shocks forming and driving large changes in density, and density inversions dissipating (due to convection) and reforming cyclically. The density inversion at the iron opacity peak is ultimately retained in the time-averaged density profile and the radiative acceleration is larger than the gravitational acceleration, creating an unbinding effect, though a mass loss or mass loss rate estimate were unable to be determined since the simulation neglected line driving and the global stellar geometry.  \nJiang et al. (2018) simulated large wedges of radiationdominated envelopes in spherical geometry, and also explored different regions of the H-R diagram compared to Jiang et al. (2015). They found that a convective insta-  \nbility at the iron opacity peak can lead to the formation of a strong helium opacity peak. Upon formation, the helium opacity peak causes local radiative acceleration to exceed the gravitational acceleration by an order of magnitude, leading to most of the overlaying mass being unbound and thereby triggering an outburst up to an instantaneous mass loss rate of up to 0 . 05 M ⊙ yr -1 (which is still ∼ 10 -100 times smaller than observed by e.g., Smith 2014). Additionally, after the envelope has settled into a steady state, convection causes large envelope oscillations during which mass is lost at a rate of ≈ 5 × 10 -6 M ⊙ yr -1 .  \nHere, instead of enhancing energy transport artificially as in MLT++ (Paxton et al. 2013) and the implicit method (Jermyn et al. 2023), we introduce a 1D model that converts local super-Eddington radiative luminosity into an average global outflow and explore the effect of this mass loss on stellar evolution.", '3.1. MESA': "We used MESA version r23.05.1 (Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023) to evolve stars between the masses of 10 -100 M ⊙ (at every 1 M ⊙ mass increment) at 1 Z ⊙ = 0 . 0142 (following Houdek & Gough 2011; Asplund et al. 2021) and 0 . 2 Z ⊙ (SMC metallicity). We used the Grevesse & Sauval (1998) metallicity abundance ratios (re-scaled to Z ⊙ = 0 . 0142), and the approx21.net nuclear network (Timmes 1999). Non-rotating models were allowed to evolve from the pre-main-sequence until the end of core helium burning (defined as central helium less than a mass fraction of 10 -4 ). All models evolve through core helium ignition, but due to aforementioned numerical issues (see Section 2), not all models reach the end of core helium burning.  \nFor convection, we used the MLT treatment following Cox & Giuli (1968) with an α MLT of 1 . 73 following Jones et al. (2013); Li et al. (2019) based on fits to the Sun (Trampedach & Stein 2011). We account for convective overshooting (Herwig 2000; Paxton et al. 2011) with a combination of step and exponential overshooting for, respectively, the top of the hydrogen burning core ( overshoot f = 0 . 345 and overshoot f0 = 0 . 01) and the top of all other convective cores ( overshoot f = 0 . 01 and overshoot f0 = 0 . 005). This closely follows the setup used to produce Figure 18 in § 7 . 2 of Jermyn et al. (2023). Our inlists are provided at this repository: DOI: 10.5281/zenodo.13306749.  \nThe MESA equation of state (EOS) combines OPAL (Rogers & Nayfonov 2002), SCVH (Saumon et al. 1995), FreeEOS (Irwin 2004), HELM (Timmes & Swesty 2000),  \nPC(Potekhin & Chabrier 2010), and Skye (Jermyn et al. 2021) EOSes. Radiative opacities combines OPAL (Iglesias & Rogers 1993, 1996) and data from Ferguson et al. (2005) and Poutanen (2017). Electron conduction opacities are from Cassisi et al. (2007). Nuclear reaction rates are from JINA REACLIB (Cyburt et al. 2010), NACRE (Angulo et al. 1999) and Fuller et al. (1985); Oda et al. (1994); Langanke & Mart'ınez-Pinedo (2000). Screening is included via the prescription of Chugunov et al. (2007). Thermal neutrino loss rates are from Itoh et al. (1996).", '3.2. Mass loss model': "In MESA , mass loss is implemented through wind massloss schemes. We employ MESA 's ' Dutch ' prescription for T eff > 4000 K, which combines the models of Vink et al. (2001) and Nugis & Lamers (2000) as described in Glebbeek et al. (2009) and captures the mass loss rate due to line-driven stellar winds ˙ M wind . For T eff < 4000 K, we adopt the more recent empirical mass loss rates of RSGs from Decin et al. (2024). This replaces the de Jager et al. (1988) or Nieuwenhuijzen & de Jager (1990) mass loss rates typically used in MESA 's ' Dutch ' prescription which are considered to overestimate mass loss rates (see Section 1). The strength of the Dutch wind is modulated by the Dutch scaling factor , which is uncertain. We choose a Dutch scaling factor of 0 . 8 (following e.g. Glebbeek et al. 2009; Fields & Couch 2021; Johnston et al. 2024; Pereira et al. 2024; Shiber et al. 2024) throughout our evolution.  \nTo account for eruptive mass loss, we create a model based on the following energy argument (see also Section 2). Motivated by Jiang et al. (2015), we define the transition between convectively-inefficient ( τ < τ crit ) and convection-dominated ( τ > τ crit ) energy transport at a critical optical depth τ crit :  \nτ crit ≡ c c s , (1)  \nwhere c s is the total sound speed. This is similar to Equation 1 in Jiang et al. (2015), except we adopt the total sound speed rather than the gas isothermal sound speed as a conservative estimate 1 for the effective energy transport by dynamical processes. The critical optical depth τ crit can be found by equating the photon diffusion timescale and the acoustic timescale for a radial distance  \n∆ r in the star, for which  \nacoustic timescale = ∆ r c s (2)  \nphoton diffusion timescale = ∆ r c/τ (3)  \nwhere c s , ∆ r , and τ are local values of, respectively, total sound speed, radial distance, and optical depth. Thus, τ crit is a local value determined at every radius. In the stellar envelope, τ = τ crit is only true at one point since both c s and τ are monotonic functions, and therefore τ = τ crit defines a unique location above which τ < τ crit everywhere in the outer envelope.  \nWe also follow Jiang et al. (2015) and Jiang et al. (2018)'s interpretation of the balance between radiative acceleration and gravitational acceleration in driving outbursts of mass loss. Specifically, in the τ < τ crit region, convection transports the convective luminosity L conv , and radiative diffusion transports the radiative luminosity L rad , but only up to the Eddington limit L Edd , defined as:  \nL Edd = 4 πGcM enc κ , (4)  \nwith M enc and κ ( r ) being, respectively, the enclosed mass and local opacity.  \nFor any locally super-Eddington radiative luminosity (i.e., for any L rad exceeding L Edd ), we regard the excess radiative luminosity above L Edd as available to do work to unbind the envelope and contribute to eruptive mass loss. More formally, since photons are able to couple with matter when τ ≳ 2 / 3, the 2 / 3 ≲ τ < τ crit region is of interest if the Eddington limit is locally exceed by the radiative luminosity. In this case, the photons in the radiatively super-Eddington regions in 2 / 3 ≲ τ < τ crit can energetically couple with matter and drive mass loss. Therefore we consider any excess radiative luminosity above Eddington ( L excess = L rad -L Edd > 0) in regions of the stellar model with τ < τ crit as available to drive mass loss. The ratio of the radiative, convective, and total luminosity to the Eddington luminosity can be expressed as Γ Edd , rad = L rad /L Edd , Γ Edd , conv = L conv /L Edd , and Γ Edd , tot = L tot /L Edd respectively.  \nWe choose to adopt the local dynamical time t dyn , since eruptive mass loss must be dynamical, to define a local excess energy E excess as  \nE excess = E excess ( r ) = L excess t dyn , (5)  \nwhere t dyn is  \nt dyn = t dyn ( r ) = √ r 3 Gm ( r ) , (6)  \nFigure 1. Profiles of luminosity, energy, and opacity as a function of mass above each radius for a M ZAMS = 30 M ⊙ , Z = 1 Z ⊙ star. The model chosen corresponds to time when log T eff / K = 3 . 76, log L/L ⊙ = 5 . 45 on the H-R diagram just after the start of core helium burning (center helium mass fraction of 0 . 98). For all panels, the horizontal axis is defined as M above = M -m ( r ). The dashed black and red lines respectively show τ crit and r mass loss as defined in Equations 1 and 11. The x-intercept of the red dashed line shows the total mass lost (see Equation 10). Top panel: local radiative, Eddington (Equation 4), and excess luminosities. Second panel: total excess energy (Equation 7) and binding energy (Equation 8). Third panel: local excess energy (Equation 5). Last panel: local opacity, showing the opacity peak at around M above ≈ 4 × 10 -4 M ⊙ due to both hydrogen and helium ionization transitions.  \n<!-- image -->  \nand r and m ( r ) are, respectively, the radius and mass coordinates. E excess is a local value at each radius r of the star.  \n<!-- image -->  \nFigure 2. Stellar temperature profiles showing locally super-Eddington radiative luminosity due to opacity peaks. Left panels: a M ZAMS = 30 M ⊙ , Z = 1 Z ⊙ model with a opacity peak due to both hydrogen and helium ionization transitions. Left panels correspond to time when log T eff / K = 3 . 76, log L/L ⊙ = 5 . 45 on the H-R diagram. Right panels: a M ZAMS = 60 M ⊙ , Z = 1 Z ⊙ model with an iron opacity peak. Right panels correspond to time when log T eff / K = 4 . 60, log L/L ⊙ = 5 . 80 on the H-R diagram. Top panels: Eddington ratios; ratios of radiative (purple) and convective (green) luminosity against Eddington luminosity. Here, radiative luminosity or convective luminosity refer to, respectively, the luminosity carried by radiative diffusion or convection. Unity is shown in a dashed black line, and the critical optical depth is indicated by a dashed red line (see text in Section 3.) Bottom panels: local opacity in each zone throughout the star (purple) and Eddington opacity (dashed black) 4 πGcM/L . Only locally super-Eddington radiative luminosity (in green) is considered in our eruptive mass loss model.  \n<!-- image -->  \nFor τ < τ crit , we find the total excess energy at every radius r from  \nE tot , excess = E tot , excess ( r ) = ∫ m ( r ) m ( r ( τ crit )) E excess dm ∫ m ( r ) m ( r ( τ crit )) dm . (7)  \nThis total excess energy is a cumulative integral from τ crit outwards, and represents a mass-average of the excess energy in the super-Eddington region. E tot , excess has a value at every radius r up to the surface and is only defined for τ < τ crit , since this is the only region we consider as possibly able to drive mass loss (see Section 2). We note that the meshing in MESA is very close to uniform in the region of interest ( τ < τ crit ), which renders the weighting effectively irrelevant.  \nWe compare this total excess energy with the gravitational binding energy of the overlaying material, which for the mass external to radius r is expressed as  \nE bind = E bind ( r ) = -∫ M m ( r ) Gm r dm , (8)  \nfor which G is the gravitational constant, m and r are the local mass and radius, dm is the local mass increment in question (at radius r ), and M is the mass at  \nthe surface (total mass). We conservatively choose not to include other sources of energy such as thermal energy, recombination energy, and kinetic energy, which are all positive energies that decrease E bind . Thus, we expect our estimated eruptive mass loss rates to be conservative.  \nFor every radius r where τ < τ crit , we define a local energy difference between the binding energy and the total excess energy as  \nE diff = E diff ( r ) = E tot , excess -| E bind | . (9)  \nWe find the inner-most radius r mass loss greater than the radius of τ crit above which E diff > 0, and then find the difference between the surface mass M and the mass at r mass loss to give ∆ M . This can be expressed as  \n∆ M = M -m ( r mass loss ) , (10)  \nwhere  \nr mass loss ≡ min( r ) > r ( τ crit ) for which E diff > 0 . (11)  \nWe emphasize that, at each timestep, E bind , E tot , excess , and E diff are functions of radius r . We show  \nFigure 3. Kippenhahn diagram for M ZAMS = 30 M ⊙ , Z = 1 Z ⊙ stars. The age x-axis focuses on the time of core He burning. This is the same model shown in Figure 1. Top panel: modeled with stellar wind and eruptive mass loss. Bottom panel: modeled only with stellar wind. Blue contours show the rate of nuclear energy generation, with colorbar on the right.  \n<!-- image -->  \nhow local luminosities (radiative, Eddington, and excess), local excess energy, and total excess energy (a cumulative integral) are accounted for when determining mass loss for our model in Figure 1, which shows a representative model of a M ZAMS = 30 M ⊙ and Z = 1 Z ⊙ star with mass loss due to a hydrogen/helium opacity peak (also shown in the left panel of Figure 2). Figure 2 shows examples of how the opacity peaks (lower panels) can cause the star's radiative luminosity to locally exceed the Eddington limit (green line in top panels).  \nFinally, we define the mass loss rate as  \n˙ M erupt = ξ ∆ M t dyn , local , (12)  \nwhere ξ is a free parameter we introduce to tune our model. t dyn , local is the local dynamical time at r mass loss ; this is the local dynamical time at the smallest radius in the star for which both E tot , excess > | E bind | (where the purple line exceeds the green line in the second panel of Figure 1) and τ < τ crit (to the right of the dashed  \nblack line in Figure 1). We explore two values of ξ in our modelling; ξ = 0 . 1 and ξ = 1 . 0.  \nFor our total mass loss model, we combine the wind mass-loss model with our novel eruptive mass loss model:  \n˙ M = ˙ M erupt + ˙ M wind . (13)  \nThis mass loss model is implemented as a custom hook in run star extras.f90 and is calculated at each timestep of a star's evolution. We note that although no explicit metallicity dependence is included in our eruptive mass loss model, the mass loss rate is expected to be indirectly impacted by changes in stellar structure due to metallicity (e.g., Xin et al. 2022) such as the effective temperature as well as the depth and strength of the helium and iron opacity peaks.  \nUsing this model, we find that our eruptive mass loss model behaves as expected; mass loss occurs due to locally super-Eddington radiative luminosity at the hydrogen/helium and iron opacity peaks. As shown in Figure 2, these opacity peaks occur in the τ < τ crit region (left of the red dashed line). In particular, the radiative luminosity (marked in green lines and denoted by Γ Edd , rad ) can be locally super-Eddington (above the black dashed line) and correspond to the locations of opacity peaks. This locally super-Eddington radiative luminosity can drive eruptive mass loss and, when ˙ M erupt is sufficiently large, can alter the structure of the star by preventing the formation of an outer convective envelope.", '4. RESULTS': 'In this section, we show our results for eruptive mass loss. Figure 3 presents Kippenhahn diagrams (mass profile over time) of stars modelled with and without eruptive mass loss in addition to stellar winds. As shown in the green hatches in the bottom panel for the star modelled without eruptive mass loss, a convective Hrich envelope is present above the helium-burning region. However, for the model with eruptive mass loss, this convective layer never forms due to eruptive mass loss at 6 . 3 Myr (in agreement with Schootemeijer et al. 2024).  \nFigure 4 more clearly shows the timing of this eruptive mass loss by presenting the mass loss rates due to stellar winds and eruptive mass loss over time, as well as the relative importance of ˙ M erupt compared to ˙ M wind . We see that the mass loss at around 6 . 3 Myr is indeed primarily due to eruptive mass loss, with ˙ M erupt exceeding ˙ M wind by two orders of magnitude.  \nAs a result of this eruptive mass loss, the convective layer is not formed and the star no longer evolves cool enough to become a RSG. This can be seen in Figure 5,  \nFigure 4. Mass loss rate, mass loss ratio, and effective temperature over time of a M ZAMS = 30 M ⊙ star with ξ = 1 . 0 and Z = 1 Z ⊙ . The middle panel shows the ratio between the mass loss rate due to eruptive mass loss and the stellar wind mass loss rate. The dashed black line marks when the ratio is unity.  \n<!-- image -->  \nWe present our full mass loss results for our 10 -100 M ⊙ stars at Z = 1 Z ⊙ in the top panels of Figure 6, which shows the evolution tracks of stars that span the full mass range on an H-R diagram. The mass loss rate due to eruptive mass loss is colored in yellow, green, and blue contours. We take particular note that significant mass loss (green and blue contours) occurs in a vertical band between 3 . 5 < log( T eff / K) < 4 . 0 as the stars evolve past the main sequence towards the red supergiant phase. Most of these models experiencing eruptive mass loss retain a surface hydrogen mass fraction between 0 . 25 -0 . 5, with surface helium and nitrogen mass fractions between 0 . 5 -0 . 75 and 1 -1 . 8 × 10 -3 respectively. The mass loss in this band is caused by the opacity peak in the stellar envelope due to both hydrogen and helium ionization transitions, as shown in Figure 2 where the peaks in Eddington ratios coincide with the hydrogen/helium opacity peak at ∼ 10 4 K (see discussion in Jiang et al. 2018). This result agrees with the 3D models in Jiang et al. (2018) where helium opacity peaks in the envelope of massive stars were found to drive mass loss.  \nwhich shows the evolutionary tracks of the two stars featured in Figure 3 on a H-R diagram. For the model with eruptive mass loss (in purple), the star remains hotter throughout its evolution and does not evolve to become a RSG. This is contrasted with the model with only stellar winds (in green), which evolves cooler to become a RSG. The left panel of Figure 5 compares models run with (purple) and without (green) eruptive mass loss in addition to stellar winds. The right panel of Figure 5 samples the evolutionary track every 1000 years, which allows better correspondence to where we expect stars to be observed. Based on this, we can see that the star modelled with eruptive mass loss is likely to be observed as a blue supergiant (BSG) whereas the star modelled with only stellar winds is likely to be observed as a RSG. Additionally, these sampled tracks capture how the star evolves across the Hertzsprung gap very rapidly in ∼ 10 , 000 years. We note that our results are not sensitive to the choice of RSG wind model since eruptive mass loss prevents evolution to the low temperatures applicable to RSG wind models.  \nComparing the left and right panels in Figure 6, we see that while all the stars evolved with ξ = 0 . 1 (left column) evolve to become RSGs, stars with mass ≳ 20 M ⊙ evolved with ξ = 1 . 0 (right column) do not become RSGs. This indicates that the presence of eruptive mass loss significantly alters the evolution of stars with mass > 20 M ⊙ , suggesting that eruptive mass loss may be a physical mechanism behind the lack of RSGs above 20 -30 M ⊙ as proposed by Massey et al. (2000); Heger et al. (2003); Levesque et al. (2009); Ekstrom et al. (2012) and the excess of BSGs (Bellinger et al. 2023). We show the change in the evolutionary tracks of our stars modeled with eruptive mass loss compared to those modeled only with stellar wind in Figure 7. From the top right panel of Figure 7, we confirm that for M ZAMS ≳ 20 M ⊙ , models evolved with eruptive mass loss at ξ = 1 . 0 do not evolve to temperatures cooler than log( T eff / K) ≈ 3 . 65 and therefore do not become RSGs.  \nAdditionally, for stars with initial masses greater than 50 M ⊙ , some mass loss of order 10 -5 M ⊙ yr -1 at higher log( T eff / K) > 4 . 5 K is predicted by our eruptive mass loss model (seen in top left of Figure 6). This is due to the iron opacity peak in the stellar envelope at ∼ 1 . 8 × 10 5 K (Paxton et al. 2013; Jiang et al. 2015, see Figure 2), and follows the results from 3D models in Jiang et al. (2015). However, mass loss in the log( T eff / K) > 4 . 5 region is dominated by stellar winds rather than eruptive mass loss. This can be seen in Figure 8, which shows the fractional increase in overall mass loss due to eruptive mass loss (total mass loss/wind  \n<!-- image -->  \nFigure 5. Evolutionary tracks of M ZAMS = 30 M ⊙ , Z = 1 Z ⊙ stars. Left panel: H-R tracks for wind only (green) and eruptive mass loss + wind (purple) models. From the differing green and purple tracks, eruptive mass loss prevents evolution to RSG phase. Right panel: evolutionary tracks sampled to 1 model every 1000 years to show where stars are expected to be observed.  \n<!-- image -->  \nmass loss). Regions in the H-R diagram in which eruptive mass loss dominates substantially over stellar winds are marked in blue, which only occurs at lower temperatures. Therefore, from Figure 8, we determine that our eruptive mass loss is dominant over stellar winds in the 3 . 5 < log( T eff / K) < 4 . 0 vertical band region.  \nSimilarly, as shown in the bottom panels in Figures 6, 7, and 8, a lower metallicity of Z = 0 . 2 Z ⊙ ≈ Z SMC also predicts a lack of red supergiants above ∼ 20 M ⊙ at ξ = 1 . 0. Since lower metallicities lead to lower opacities, increasing L Edd and thereby reducing Γ Edd , rad = L rad /L Edd , some of the numerical challenges involving locally super-Eddington regions in the star (see Section 2) are avoided and the high-mass Z = 0 . 2 Z ⊙ ≈ Z SMC models are able to evolve further than those with Z = 1 Z ⊙ . The lower panels of Figure 6 shows that at both metallicities, stars evolved with ξ = 1 . 0 exhibit higher rates of eruptive mass loss of up to ∼ 10 -2 M ⊙ year -1 whereas those evolved with ξ = 0 . 1 feature correspondingly lower mass loss rates of up to ∼ 10 -3 M ⊙ year -1 . Unsurprisingly, the region above the Humphreys-Davidson (HD) limit exhibits high levels of eruptive mass loss, as this observational limit corresponds to regions on the H-R diagram where surface stellar luminosities transition to become superEddington. Additionally, this region above the HD limit also exceeds the de Jager limit (de Jager 1984), wherein supersonic turbulent motions in the photosphere may have a destabilizing effect. Models at Z = 1 Z ⊙ exhib-  \nmass loss rates by up to an order of magnitude than those at Z = 0 . 2 Z ⊙ for M < 40 M ⊙ . For M > 40 M ⊙ , the Z = 0 . 2 Z ⊙ models were able to evolve further and therefore exhibited higher maximum eruptive mass loss rates by at least an order of magnitude compared to Z = 1 Z ⊙ models. The maximum eruptive mass loss rate across all models is ≈ 3 × 10 -2 M ⊙ year -1 (occurring for ξ = 1 . 0), which is consistent with observations (see Section 1).', '5. DISCUSSION AND CONCLUSIONS': "We studied eruptive mass loss in massive stars of varying masses (between 10 -100 M ⊙ ) at two different metallicities ( Z SMC ≈ 0 . 2 Z ⊙ and 1 Z ⊙ ) with an implementation of an energy-based eruptive mass loss model in the stellar evolution code MESA . Our model's mass loss rate is determined by considering photon diffusion through regions in the outer stellar envelope of these massive stars. As parts of the outer envelope may be locally super-Eddington, the excess (above Eddington) energy output from radiation may act on the fluid mass faster than dynamical, convective, and acoustic processes and exceed the binding energy of the overlaying layers, leading to mass loss.  \n<!-- image -->  \n<!-- image -->  \nFigure 6. Eruptive mass loss rates across the H-R diagram. Left column: ξ = 0 . 1. Right column: ξ = 1 . 0. Top row: Z = 1 Z ⊙ . Bottom row: Z = Z SMC = 0 . 2 Z ⊙ . Grey lines show the evolutionary tracks of stars from the beginning of main sequence until the end of the model run. Models were run at every 1 M ⊙ increments but only selected models are shown here for clarity. Dashed grey lines show two empirical tracks of Quiescent LBVs and LBVs in Outburst for reference. Solid red lines show the Humphreys-Davidson limit for reference. Colored contours show the rate of eruptive mass loss experienced by the star, with legend on the right.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 7. Evolutionary tracks of models with and without eruptive mass loss. Left column: ξ = 0 . 1. Right column: ξ = 1 . 0. First row: Z = 1 Z ⊙ . Second row: Z = Z SMC = 0 . 2 Z ⊙ . All axes and labels follow Figure 6.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 8. Ratio of total mass loss against wind mass loss rates across the H-R diagram. Left column: ξ = 0 . 1. Right column: ξ = 1 . 0. Top row: Z = 1 Z ⊙ . Bottom row: Z = Z SMC = 0 . 2 Z ⊙ . Colored contours here show the ratio between total mass loss and wind mass loss, with legend on the right. Cooler colors (blues) indicate that eruptive mass loss is more dominant over stellar wind. All axes, labels, and other lines follow Figure 6.  \n<!-- image -->  \nWe find that, at ξ = 1 . 0 (a model efficiency of one for the eruptive mass loss) and at both metallicity choices, stars with initial masses greater than 20 M ⊙ no longer evolve to become red supergiants. This suggests that eruptive mass loss could be a physical mechanism that contributes to the observed lack of RSGs above 20 -30 M ⊙ , despite recent careful work determining relatively low mass-loss rates in well-characterized isolated RSGs (Beasor et al. 2020, 2023; Decin et al. 2024).  \nAlternative explanations other than eruptive mass loss for the lack of RSGs above 20 -30 M ⊙ include envelope inflation (e.g., Sanyal et al. 2015, 2017) and binary interactions (e.g., Smith 2014; Morozova et al. 2015; Margutti et al. 2017). In the case of envelope inflation, Sanyal et al. (2015) ran 1D hydrodynamic stellar evolution models on core hydrogen burning massive stars and found that stars with initial mass M > 40 M ⊙ locally exceeded the Eddington limit in the partial ionisation zones of iron, helium, or hydrogen. The stars featured density inversions in their envelopes and experienced an increase of radius of up to a factor of ∼ 40, which alters the temperature and density structure such that envelope opacity decreases until the Eddington limit is no longer exceeded. While the authors found that the inflated envelopes remained bound, they noted that their models may have overestimated convective energy transport due to their implementation of MLT (therefore reducing L rad ) and that multi-dimensional hydrodynamical simulations are needed to fully determine the stability of density inversions in stellar envelopes. Additionally, the authors speculatively mentioned that eruptions are a possible consequence of their very coolest models. Alternatively, binary interactions such as envelope stripping can lead to the formation of Wolf-Rayet stars that may never enter the RSG phase (e.g., Dray & Tout 2003; Meynet et al. 2011; Grafener et al. 2012). However, mass-transferring binaries leave behind accretors and mergers, many of which are expected to evolve to the RSG phase (e.g., Podsiadlowski et al. 1993; Aldering et al. 1994; Maund et al. 2004; Claeys et al. 2011; Yoon et al. 2017; Schneider et al. 2020; Renzo & Gotberg 2021; Renzo et al. 2023).  \nAdditionally, our work may offer a possible solution to the 'Blue Supergiant Problem' wherein many blue supergiants (BSGs) having been observed despite the expectation of rapid evolution through the BSG phase. BSGs are hot, luminous stars with luminosities up to log( L/L ⊙ ) ≈ 6 . 0 and temperatures up to log( T eff ) / K ≈ 4 . 4 (e.g., HD 5980 by Pickering & Fleming 1901; Cannon & Pickering 1918; Koenigsberger et al. 2010). It is generally thought that BSGs are either due to changes in stellar structure via mixing during evolution across  \nthe Hertzsprung gap (e.g., Iben 1993; Sugimoto & Fujimoto 2000; de Mink et al. 2009, 2013; Schootemeijer et al. 2019), or post-main-sequence interactions such as mass transfer or mergers (e.g., Podsiadlowski et al. 1992; Braun & Langer 1995; Justham et al. 2014; de Mink et al. 2014). The eruptive mass loss presented in this work can be regarded as a possible intrinsic post-mainsequence behavior that may contribute to the formation of BSGs.  \nThe energy-based model for eruptive mass loss presented in this work draws parallels with the 'geyser model' presented in Maeder (1992). In this model, the mass loss rate during a shell ejection event in a massive star is estimated by assuming the mass above super-Eddington peaks is ejected and accounting for the inward-movement of these peaks at the local thermal timescale. Our model does not explicitly account for this effect since it would operate at timescales shorter than our code timestep. As such, the 'geyser model' may further enhance eruptive mass loss, though our model may approximately account for an average of such effects.  \nWe emphasize that due to the numerical challenges of modelling massive stars in 1D (see Section 2), our models were only evolved far enough into core helium burning to probe the rightmost section of the HR diagram where red supergiants are found. Nevertheless, our work demonstrates that a physically-motivated energybased mass loss model can adequately describe eruptive mass loss in massive stars. Additionally, we note that our work here treated the mass loss model's efficiency factor ξ as a free parameter and explored the implications of ξ = 1 . 0 and ξ = 0 . 1. An observational study aimed at tuning this efficiency factor using LBV mass loss observations will be conducted in future work.", 'ACKNOWLEDGMENTS': 'SJC and CC acknowledge support from NSF grant AST-131576. The Flatiron Institute is supported by the Simons Foundation.', 'APPENDIX': "For RSG winds, we choose to adopt Decin et al. (2024) rather than one of the schemes included in MESA 's Dutch prescription (e.g., de Jager et al. 1988; Nieuwenhuijzen & de Jager 1990) for which mass loss rates are overestimated (see Section 1). While doing so alters stellar evolution in the RSG phase, it does not affect our eruptive mass loss result: the left panel of Figure 9 compares models run with Decin et al. (2024) and with de Jager et al. (1988) mass loss for T eff < 4000 K. The Decin et al. (2024) winds only affects the very late RSG phase, and does not impact the model with eruptive wind in the case shown (in purple) since such models do not evolve cool enough to enter the Decin et al. (2024) regime. We note that although Decin et al. (2024) produces mass loss rates calibrated to careful observations in well-characterized individual stellar populations (in broad agreement with Beasor et al. (2023)'s correction to the Beasor et al. (2020) rates), it contains in it a dependence on the initial mass of stars as estimated from cluster characteristics, which need not correspond to the initial mass of a stellar model in MESA . A more proper treatment correcting for this when comparing models to stellar populations might include a variable scaling fac-  \ntor akin to the dutch\\_scaling\\_factor for the ' Dutch ' stellar wind prescription.  \nFigure 9. Evolutionary tracks of M ZAMS = 30 M ⊙ , Z = 1 Z ⊙ stars. H-R tracks for wind only (green) and eruptive mass loss + wind (purple) models. The models marked with 'Decin' contain the modification described in Section 3, while those using the old de Jager et al. (1988) rates are marked with 'de Jager'. From the differing green and purple tracks, eruptive mass loss prevents evolution to RSG phase.  \n<!-- image -->", 'REFERENCES': 'Aldering, G., Humphreys, R. M., & Richmond, M. 1994, AJ, 107, 662, doi: 10.1086/116886 Angulo, C., Arnould, M., Rayet, M., et al. 1999, NuPhA, 656, 3, doi: 10.1016/S0375-9474(99)00030-5 Asplund, M., Amarsi, A. M., & Grevesse, N. 2021, A&A, 653, A141, doi: 10.1051/0004-6361/202140445 Banerjee, S., Belczynski, K., Fryer, C. L., et al. 2020, A&A, 639, A41, doi: 10.1051/0004-6361/201935332 Beasor, E. R., & Davies, B. 2018, MNRAS, 475, 55, doi: 10.1093/mnras/stx3174 Beasor, E. R., Davies, B., Smith, N., et al. 2020, MNRAS, 492, 5994, doi: 10.1093/mnras/staa255 -. 2023, MNRAS, 524, 2460, doi: 10.1093/mnras/stad1818 Belczynski, K., Bulik, T., Fryer, C. L., et al. 2010, ApJ,  \n714, 1217, doi: 10.1088/0004-637X/714/2/1217  \nBellinger, E. P., de Mink, S. E., van Rossem, W. E., & Justham, S. 2023, arXiv e-prints, arXiv:2311.00038, doi: 10.48550/arXiv.2311.00038 Blaes, O., & Socrates, A. 2003, ApJ, 596, 509, doi: 10.1086/377637 Bohm-Vitense, E. 1958, ZA, 46, 108 Braun, H., & Langer, N. 1995, A&A, 297, 483 Bruch, R. J., Gal-Yam, A., Yaron, O., et al. 2023, ApJ, 952, 119, doi: 10.3847/1538-4357/acd8be Cannon, A. J., & Pickering, E. C. 1918, Annals of Harvard College Observatory, 91, 1 Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M. 2007, ApJ, 661, 1094, doi: 10.1086/516819 Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157, doi: 10.1086/153315'}