Eruption of the Envelope of Massive Stars by Energy Injection with Finite Duration

A significant fraction of supernovae show signatures of dense circumstellar material (CSM). While multiple scenarios for creating a dense CSM exist, mass eruption due to injection of energy at the base of the outer envelope is a likely possibility. We carry out radiation hydrodynamical simulations of eruptive mass loss from a typical red supergiant progenitor with initial mass of $15\ M_\odot$, for the first time focusing on the timescale of the injection as well as energy. We find that not only sufficient injection energy but also sufficient rate of energy injection per unit time, $L_{\rm{min}} \sim 8\times 10^{40}$ erg s$^{-1}$ in this particular model, is required for eruption of unbound CSM. This result suggests that the energy injection rate needs to be greater than the binding energy of the envelope divided by the dynamical timescale for the eruption. The density profile of the resulting CSM, whose shape was analytically and numerically predicted in the limit of instantaneous energy injection, similarly holds for a finite injection timescale. We discuss our findings in the framework of proposed mass outburst scenarios, specifically wave-driven outbursts and common envelope ejection.


INTRODUCTION
Many massive stars explode at the end of stellar evolution, which are observed as supernovae (SNe). SNe display a large variety of spectra and light curves, mainly reflecting the composition of the progenitors and their environments. Although emission lines of SNe are usually broad, some SNe show very narrow hydrogen lines in the spectra. These kinds of SNe are classified as Type IIn, first named by Schlegel (1990). The existence of the narrow lines strongly suggests that there is a dense circumstellar medium (CSM) slowly moving outside the star.
Assuming that the CSM originates from the SN progenitor, there should be violent mass loss at the late stage of the stellar evolution that forms dense CSM before the SN (e.g., Chugai et al. 2004). The mass loss rates of SNe IIn progenitors are inferred to be in the range of 0.026 − 0.12 M yr −1 (Kiewe et al. 2012), 10 −4 − 10 −2 M yr −1  or larger than 10 −3 M yr −1 (Moriya et al. 2014). These values are much more extreme than standard radiation-driven stellar winds, which generally have mass loss rates oḟ M 10 −5 M yr −1 (e.g., Vink et al. 2001;Smith 2014). It is also observationally known that a large fraction of Type IIn SNe progenitors temporarily exhibited outbursts just before (∼ 1-10 years) the SN explosion (e.g., Ofek et al. 2014), implying that these events may be responsible for creating the dense CSM.
Despite the increasing number of observations of SNe IIn, the theoretical mechanism of such violent mass loss is not fully understood. There are various scenarios suggested to explain such mass loss. For example, Smith & Arnett (2014) proposed unsteady nuclear burning due to the turbulent convection at the late stage of stellar evolution. Chevalier (2012) proposed the release of gravitational energy due to the common envelope evolution with a binary companion, likely a neutron star or a black hole. Woosley et al. (2007) proposed pulsational pair instability arising from the production of electron-positron pairs in the core of very massive stars.
For many of the models raised here, eruption of the progenitor's envelope is triggered by injection of energy into the bottom of the envelope. Such process was recently studied by Kuriyama & Shigeyama (2020) with radiation hydrodynamical simulations, and Linial et al. (2021) with analytical modeling. They showed that there is a minimal injection energy required for mass loss, that depends on the binding energy of the progenitor. However, these calculations have assumed the energy injection to occur instantaneously, i.e. in a much shorter time compared to the envelope's dynamical timescale. Since such an ideal approximation would not always hold, it would be important to consider the effect of having a finite injection duration. A longer duration would make the gravitational pull from the center more important, which would reduce the mass and energy of the eruption. Hence we expect there to be a constraint on the duration on the injection, as well as the energy, for mass eruption to occur.
Modelings of continuous energy injection onto the stellar envelope were recently considered by several works (Leung & Fuller 2020;Leung et al. 2021a,b), in the context of wave-driven mass loss. In Leung & Fuller (2020), the amount and asymmetry of the mass loss was investigated by two-dimensional hydrodynamic simulations of hydrogen-rich progenitors. However, the simulations used in their work did not take into account radiation transfer. This may overestimate the mass loss, because loss of energy due to the escape of radiation around shock breakout is not included. In Leung et al. (2021a), the dependence of mass loss on the injection energy and duration was investigated taking into account the radiation transfer effect. This work focused on parameterizing the energy injection onto a hydrogen-poor progenitor in order to reproduce a specific Type Ic SN 2018gep. Leung et al. (2021b) improved the modelling of wave heating and investigated the mass loss, but this work also focused on hydrogen-poor progenitors which are irrelevant to SNe IIn. Moreover, common to these works is that the specific conditions necessary for mass loss to occur were not calculated. A model-agnostic criterion of mass loss, as well as its dependence on the injection energy and duration, would be useful to constrain these parameters from observations of SNe with massive CSM and test the proposed models of violent mass loss.
In this work, we calculate the mass eruption process with radiation hydrodynamical simulations, varying the duration of the energy injection. We simulate the mass eruption of a red supergiant (RSG) progenitor, and compare the amount of mass loss as well as the final density profile of the resulting CSM. We find that not only a sufficient injection energy but also a sufficiently short timescale, within the order of weeks to months for RSGs, is required. This paper is organized as follows. In Sect. 2, we introduce our progenitor model generated with MESA and the simulations done in this work. In Sect. 3, we present the results of our calculations, with focus on the amount of mass lost. In Sect. 4, we conclude and discuss our results in the context of several proposed mechanisms for mass eruption.

METHODS
We use the mass eruption part (Kuriyama & Shigeyama 2020) of the open-source code CHIPS  to simulate the eruption of the envelope of a massive star years before core-collapse. In this section we describe the basic setup and key features of our calculations.

Progenitor Model
We adopt a RSG progenitor with an initial mass of 15 M , generated in Tsuna et al. (2021) (their R15) using the revision 12778 (example_make_pre_ccsn test suite) of MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019 and released as one of the sample progenitor models in CHIPS. The basic parameters of our RSG model are listed in Table 1. As in Kuriyama & Shigeyama (2020), since we are interested in the (partial) ejection of the envelope, we extract the hydrogen-rich envelope of this progenitor as the computational region. We inject energy into the inner boundary at r inner ≈ 1.7 × 10 12 cm, which correspond to the bottom of the envelope. We use the progenitor that was evolved until core-collapse, but the change in the envelope over the last decade is found to be negligible.
The HELMHOLTZ equation of state 1 (Timmes & Swesty 2000) is used to calculate thermodynamic quan-tities from the internal energy density and ρ, which are obtained from the time integration of the hydrodynamical equations. For the range not covered by this equation of state, the pressure is assumed to be the sum of that of ideal gas and that of black-body radiation, i.e., where R denotes the gas constant, and µ the mean molecular weight. When solving the hydrodynamical equations, the boundary conditions are given as follows: v inner = 0, r inner = 1.657 × 10 12 cm, p outer = 0. (14) Here the subscript inner (outer) indicates quantities at the base (outer edge) of the envelope. In contrast to Kuriyama & Shigeyama (2020) that considered a fixed energy injection timescale much shorter than the envelope's dynamical time, we consider the injection timescale as a free parameter. The energy is injected at a constant rate as where E inj is the total injected energy. The energy is equally injected in the innermost 10 cells of the computational region. Following Tsuna et al. (2021) and Takei et al. (2021), we characterize E inj by the parameter f inj , which is the ratio of E inj to the initial total binding energy of the outer envelope E outer in Table 1. The earlier works by Kuriyama & Shigeyama (2020) found that for RSGs f inj of a few 10% or more results in (partial) ejection of the envelope. We investigate 22 different parameter sets of (f inj , ∆t inj ) in this range, as shown in Table 2, and follow the hydrodynamical evolution of the envelope for 5 years.

RESULTS
The energy injection generates an outward shock wave that propagates in the envelope of the progenitor (see e.g. Fig. 3 of Kuriyama & Shigeyama 2020). When the shock approaches the star's surface it breaks out, and a part of the envelope is ejected to form a CSM. In this section we present results of our simulations, focusing on the shock propagation, mass ejection and the final density profile of the CSM.

Shock Propagation in the Star
We extract the time evolution of the shock from the output files of the CHIPS code 2 that record the hydrodynamical quantities at each cell. The evolution of the 2 resultXX.txt, where XX is an integer from 00 to 99.
shocks for various parameter sets is shown in Figure 1. We have fit the radius of the shock as a function of time r shock (t) with a third-order polynomial using the numpy.polyfit module, which results in an error within 0.3% in radius.
Due to the energy injection a shock is formed and accelerated near the base of the envelope, but it gradually decelerates as more material is swept up as shown in Figure 1. Figure 1 shows that the shock propagation in the star initially follows a power law of r shock ∝ t 0.6 , especially for ∆t inj 's much shorter than the dynamical timescale of the envelope This dependence is found to be consistent with the Sedov-Taylor self-similar solution (Taylor 1950;Sedov 1959), where the expansion of the blast wave for our initial RSG model (with a density profile of approximately ρ ∝ r −1.7 near the inner boundary) should follow r shock ∝ t 0.61 .
As the shock approaches the surface, the upstream density quickly decreases and the shock is slightly reaccelerated. Eventually the shock reaches the breakout radius, where the optical depth τ satisfies v shock = dr shock /dt ∼ c/3τ and the internal energy stored in the shock downstream starts to be lost by radiative diffusion. This happens 1-3 months after energy injection, depending on the model parameters.

Mass Eruption
The material released from the star upon shock breakout expands and reduces the pressure by converting its internal energy to kinetic energy. When the pressure in the material becomes negligible, the motion is exclusively controlled by the central star's gravity. At this point, the material is composed of an unbound part at the outer region and a bound part at the inner region that eventually falls back towards the central star. Thus we focus on the unbound CSM. Figure 2 shows the dependence of the mass of the unbound material on ∆t inj , obtained from our simulations in the range f inj = 0.3-0.8. It can be seen that while the mass is rather insensitive to ∆t inj at low values, the energy injection with a ∆t inj longer than a few ×10 6 s significantly reduces the amount of ejected envelope. This is because the gravitational force from the center becomes important when ∆t inj is comparable to the dynamical timescale at the outer edge of the hydrogen envelope (see Eq. (16)). A longer ∆t inj would result in deceleration of the shock, hence a lower amount of mass  to be accelerated to velocities sufficient for being gravitationally unbound.
Since the total energy of the envelope after energy injection (1 − f inj )(−E outer ) is still negative, a very long ∆t inj should result in no mass ejection. As shown in Figure 2, we find that the maximum value ∆t inj,max of ∆t inj for mass eruption, is of order months and shorter for smaller f inj . Figure 3 shows the relation between f inj and ∆t inj,max , and illustrates that there is a minimum luminosity required for mass loss to occur. The relation is linear with a slope of ∆t inj,max /f inj = 6 × 10 6 s, which is close to the dynamical time scale t dyn ∼ 8 × 10 6 s. This suggests that the following equation roughly holds: Figure 3. Relation between finj and ∆tinj,max obtained from our simulations (orange dots and black crosses). The orange points (black crosses) are where unbound CSM is found to form (not form) in our calculations. The relation is approximately linear (blue line) at large finj with slope ∆tinj,max/finj = 6 × 10 6 [s], while it transitions at finj ≈ 0.2 (Einj ≈ 10 47 erg), the minimum finj required for mass loss. We predict mass eruption to occur only in the shaded region of the parameter space.
In the case of our RSG progenitor, the constraint of ∆t inj,max /f inj for mass loss corresponds to an injection luminosity of 7.9 × 10 40 erg s −1 .
In order for the shock to reach the surface of the star, a certain amount of energy is required, because some of the energy involved in the shock is lost due to ascending the gravitational potential during its propagation and escape of radiation around breakout. This condition should appear as a different component from the blue line in Figure 3. Linial et al. (2021) and Kuriyama & Shigeyama (2020) found that there exists a minimal energy required for RSG progenitors to achieve eruptive Figure 4. The shock and photon diffusion velocities (c/3τ ) and opacity near the photosphere of the star as functions of the normalized depth from the pohotosphere. The two shocks are generated with the same finj = 0.5 but with different ∆tinj = 10 4 s and 2 × 10 6 s. The opacity is calculated from equation (8) using the density and temperature profile of the star just before shock breakout. mass loss, of the order of 10 46 -10 47 erg. We repeat the simulation for lower values of f inj , with an instantaneous energy injection of ∆t inj = 10 3 s and check the existence of minimum injection energy (E inj 8 × 10 46 erg) as shown in Figure 3 which is roughly consistent with the both works. By simulating f inj just above this lower limit, we additionally found that this appears as a hard cutoff near the lower limit, shown by the green line in Figure 3.
We conclude that there exists both a minimum energy and minimum luminosity for mass loss to occur. In Sect. 4 we consider the implications of our results in the context of several models proposed for mass eruption.

Analytical Investigations of Erupted Mass
The value of the mass of the erupted material can not only be obtained by numerically following the mass eruption for many years like done above, but can also be roughly inferred by simpler analytical approaches. We introduce such inferences and check the consistency between the numerical results.
First, there exists a lower limit of the ejecta mass, depending only on the progenitor, because mass ejection occurs only when the material at and beyond the breakout radius swept up by the shock can obtain a velocity exceeding its escape velocity (Linial et al. 2021). Assuming that the opacity κ is constant near the surface, this mass is expressed as where C ≈ 2. This coefficient C is actually a function of the depth from the photosphere (see also Eq. (21)). For the outer layer of RSGs, the opacity is however sensitive to the temperature due to the partial ionization of hydrogen. Figure 4 shows the radial dependence of the shock velocity, photon diffusion velocity and the opacity κ (Eq. (8)), around the breakout radius where τ ∼ c/3v shock .
The shock velocity (nearly horizontal solid lines) is of the order of 100 km s −1 for all cases, and thus the breakout radius is close to the photosphere (r ≈ 0.97R ). Around this radius the opacity rapidly changes from ≈ 10 cm 2 g −1 to < 1 cm 2 g −1 . If we set κ = 1 cm 2 g −1 as a fiducial value, we then obtain the lower limit of the mass as This result is consistent with the mass loss from the hydrodynamics simulations, where we obtain the lowest finite mass of the unbound CSM to be 3.2 × 10 −2 M from a model with a parameter set of f inj = 0.3 and ∆t inj = 1 × 10 6 s.
Second, one can infer the unbound mass from the evolution of the shock velocity depicted in Figure 1. The material that crosses the shock is heated, accelerated, and acquires an outward velocity upon mass eruption. The final speed v eje,fin of this material is determined by the structure of the envelope as well as the shock velocity, and can be written as (e.g., Matzner & McKee 1999;Ro & Matzner 2013 where both of C and v shock are functions of the normalized depth from the photosphere of the star x 0 = 1 − r shock /R * . For an envelope with a polytropic index 3/2 appropriate for RSGs, an approximate fit to C(x 0 ) is given as which we note is independent of our model parameters. From our polynomial fit of r shock (t), we can obtain v shock = dr shock /dt and thus v eje,fin (x 0 ). An important assumption is that gravity from the central star is neglected when determining the parameter C. We can nevertheless compare this v eje,fin (x 0 ) with the escape velocity v esc = 2Gm/r at the same depth, and determine whether (and roughly how much) mass can be ejected.
To illustrate this, Figure 5 shows a comparison between v eje,fin (solid lines) and v esc (dotted line) for two parameter sets. The figure indicates that a longer ∆t inj would result a narrower region where v eje,fin exceeds v esc , i.e. a more gentle mass loss as seen in Sect. 3.2. From the figure the lower solid line of ∆t inj = 4.8 × 10 6 sec results in no unbound CSM, while the upper solid line with shorter ∆t inj = 10 6 sec results in about 2M of unbound CSM. These are consistent with the results obtained from numerical simulations within a factor of a few.
However, we note that this method may give inaccurate estimates of the unbound CSM mass when v eje,min is marginally greater than v esc and only a little mass close to the breakout radius is ejected. One reason is that the calculation of C(x 0 ) assumes adiabatic expansion. This might lead to an overestimation if radiation can efficiently escape from the shock downstream. Another issue may arise from a small but finite error in the polynomial fitting for obtaining v shock . In our fitting, r shock was found to have an error of a few ×10 10 cm, and the interval of the time sample for the fitting is ≈ 3 × 10 5 s, thus v shock contains an error of ∼ 1 km s −1 at each timestep. The escape velocity is somewhat insensitive to the mass measured from the surface, with a relative change of 1 km s −1 in ∼ 10 −2 M . This fitting error thus prevents this method from accurately obtaining the CSM mass if it is around this resolution.
Third, the amount of unbound matter due to the wave-driven mass loss is analytically calculated in Matzner & Ro (2021), neglecting energy loss due to radiation and assuming that the injection time is much shorter than the dynamical timescale. Matzner & Ro (2021) found that the amount of the unbound matter is proportional to f θ inj . Here, θ is expressed as where β ≈ 0.19, γ is the adiabatic index, and n is the polytropic index. Adopting (γ, n) = (5/3, 3/2) appropriate for RSG models, we obtain θ ≈ 2.09. For the case where ∆t inj t dyn , the combinations of injection energy and the amount of the unbound matter in our work (f inj , unbound mass) are (0.8, 0.99M ), (0.5, 0.31M ) and (0.3, 0.032M ). Our results at f inj = 0.5 (0.3) predict a factor of 1.2 (4.0) lower mass than that predicted from the dependence ∝ f 2.09 inj normalized at f inj = 0.8. We presume this to be because the smaller the injection energy, the more significant is the energy loss due to radiation.

Final Density and Velocity Profiles of the CSM
The density profiles of the final CSM for several parameter sets are shown in the top panel of Figure 6. When f inj is fixed, the CSM generated with a shorter ∆t inj extends to larger distances, due to the slightly larger shock speed around shock breakout.
The shape of the density profile appears to be a double power law with the inner region following ρ ∝ r −1.5 (solid lines in Figure 6), in agreement with that derived analytically in Tsuna et al. (2021) and numerically in Kuriyama & Shigeyama (2020) in the case of instantaneous energy injection. The inner power law originates from the fallback of the bound part of the CSM, which is rather insensitive to the energy injection history. The discontinuity in Figure 6 at r ∼ 10 14 cm is the shock created by the collision of the falling CSM and the progenitor. We thus conclude that this density profile of the CSM is a robust outcome independent of the mass eruption mechanism, given that a sufficient time passes for the bound CSM to start falling back.
The velocity profile of the final CSM for several parameter sets are shown in the bottom panel of Figure 6. The outermost velocities are around 100 km s −1 , slightly larger than the escape velocity of the original RSG progenitor. This is around the lower end of the observed CSM velocity from Hα lines of SNe IIn, of about 100-1000 km s −1 (Kiewe et al. 2012;. We note that the observed CSM velocity may be affected by acceleration due to radiation emitted in the previous phase of the SN. Alternatively, a more compact hydrogen-rich progenitor like blue supergiants may be responsible for a fraction of the explosions with highest CSM velocities. In this paper, we investigated how eruptive mass loss of a massive star's envelope is affected with the timescale of the energy injection at its base. We found that not only sufficient energy is required for eruptive mass loss, but also that the energy must be released within a certain time period, as shown in Figure 3. This is qualitatively consistent with Leung et al. (2021a), which put limitation on the energy injection for mass-loss to occur for a hydrogen-free progenitor. This limits the injection luminosity to be L min 8 × 10 40 erg s −1 to cause mass loss for our typical RSG stellar model. For an injection energy of ∼ 10 47 erg, this limits the injection timescale to be within a few weeks to reproduce a IIn-like CSM.

Application to Proposed Mass-loss Mechanisms
We apply our results to several proposed mass loss mechanisms. Recently Wu & Fuller (2021) examined pre-SN mass loss due to heating of the envelope by gravity waves excited at the inner core. Their calculations for progenitor models with initial masses of 15M show that during the neon/oxygen burning phase, an energy of 10 47 erg is supplied with a luminosity of 10 41 -10 42 erg s −1 , i.e. an energy injection timescale of 10 5 -10 6 s. Comparing this with Figure 3, we find that this energy injection can indeed result in eruptive mass loss. However, it should be noted that for wave-driven mass loss the energy dissipation occurs near the base of the star's outer envelope. Uncertainties in the position of the energy deposition not considered in our work can alter the binding energy, and hence can largely affect the conditions for mass eruption.
Mass ejection may also occur through binary interactions, e.g. when the massive star forms a common envelope with a compact object (CO). Such mechanism has been discussed as the origin of the massive CSM (e.g., Chevalier 2012), and may be followed by an explosion, with energy comparable to or even exceeding SNe, when the CO further spirals into the star's helium core (e.g., Fryer & Woosley 1998;Barkov & Komissarov 2011). When a CO spirals in the common envelope, the orbital energy between the CO and the star's core can be (partially) used to expel the envelope. The fraction of the dissipated energy that can be used for envelope ejection (the α parameter; Livio & Soker 1988;Ivanova et al. 2013) is expected to be α < 0.7, or even α ∼ 0.1 (Klencki et al. 2021, and references therein). Assuming that the CO spirals in at the base of the envelope and the orbital energy E orb is dissipated within roughly the period t Kep of the Kepler motion there, the energy injection rate E inj /∆t inj can be crudely expressed as The orbital quantities scale as Adopting the masses of the progenitor's He core and the CO to be 5 M and 1 M , respectively and the radius at the base of the envelope to be 10 12 cm, E orb /t Kep is on the order of 10 42 erg s −1 . Thus, for a range of α ≈ 0.1-1 mentioned in Klencki et al. (2021), we find from Figure 3 that this energy injection can indeed result in eruptive mass loss. However, our results are based on a one-dimensional study assuming spherical symmetry. The above constraint can therefore be somewhat modified since the common envelope process is clearly multi-dimensional.

Possible Caveats
We conclude by brefly discussing several possible caveats of our work.
We have calculated the eruptive mass loss from a RSG using 1D models under the assumption of spherical symmetry and found conditions for the energy injection to lead to mass eruption. First, these conditions are obtained for a particular RSG model. Though we have extended our results to more general contexts, in reality we should investigate other models to prove whether the conditions really work in general contexts. Nevertheless, we have obtained the values of E outer /t dyn in equation (17) for other sample stellar models (13 − 20 M ) available in CHIPS , and find that they are within 30% of the 15 M model used in this work. Thus we expect that our conditions will be applicable to RSGs in general.
Second, there is growing evidence for asymmetry in the CSM of some SNe IIn progenitors, e.g., from observations of polarization (e.g., Burrows et al. 1995;Plait et al. 1995;Larsson et al. 2016). For example, a disklike CSM was suggested from the optical spectropolarimetry of SN 1998S (Leonard et al. 2000). Such strongly aspherical CSM may be difficult to reproduce by massloss under the wave heating mechanism, even though multi-dimensional effects such as Rayleigh-Taylor instabilities operate upon energy deposition (Leung & Fuller 2020). This may hint a different energy injection mechanism at play for the CSM of SN 1998S, such as binary interactions as discussed above. Investigating the outcome of the CSM structure with a more model-agnostic aspherical energy injection would be an interesting future work.
Third, the energy injection was assumed to occur only once. However past observations of some SNe have found signs of multiple mass eruptions before the terminal explosion exemplified by SN 2009ip (e.g., Pastorello et al. 2013;Prieto et al. 2013;Ofek et al. 2013;Margutti et al. 2014). Energy injection can also significantly alter the density structure of the star on the thermal timescale, which might not be recovered at the onset of the next eruption. After the envelope inflates, eruption can become substantially easier due to the weakened gravitational binding of the envelope (Ouchi & Maeda 2019;. Therefore, if the energy injection is repeated, the condition for eruption can be different for each injection. Such multiple eruptions should be the subject of future research. Software: MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019, CHIPS , NumPy (Harris et al. 2020) version 1.22, Matplotlib (Hunter 2007) version 3.5.