Hydrodynamic Simulations of Pre-supernova Outbursts in Red Supergiants: Asphericity and Mass Loss

Type II supernovae (SNe) are the explosions of massive hydrogen-rich stars, and it is widely accepted that the progenitors of Type II-P SNe are red supergiants. Observations indicate that progenitor stars of interacting Type IIn SNe usually produce pre-SN outbursts (Ofek et al. 2014), specific examples including SN 2009ip (Mauerhan et al. 2013), SN 2010mc (Ofek et al. 2013), and SN 2015bh (Elias-Rosa et al. 2016; Ofek et al. 2016). However, even fairly typical Type II-P SNe show evidence for circumstellar material (CSM) that may arise from pre-SN mass loss (Khazov et al. 2016; Yaron et al. 2017). Pre-SN mass eruption events have also been observed or hinted at in other types of SNe, including Type IIb such as SN 2013cu (Gal-Yam et al. 2014; Groh 2014; Gräfener & Vink 2016), Type Ibn such as SN 2006jc (Pastorello et al. 2007), SN 2015G (Shivvers et al. 2017), and SN 2015U (Shivvers et al. 2016), and even broad-lined Type Ic such as PTF11qcj (Corsi et al. 2014) and SN 2018gep (Ho et al. 2019).

Pre-SN mass ejection has been proposed to explain many features of both ordinary and superluminous SNe. The interaction between the ejecta from the final core-collapse explosion and the previously ejected CSM allows the efficient conversion of the ejecta kinetic energy to thermal energy by shock heating. A classical example is pulsation pair-instability SNe (Woosley et al. 2007; Yoshida et al. 2016; Woosley 2018; Leung et al. 2019a). Shock waves driven by unstable oxygen burning eject large amounts of mass prior to the star's collapse, creating dense CSM that allows for a superluminous SN (Sorokina et al. 2016; Morozova et al. 2017; Tolstov et al. 2017).

For less massive stars, wave energy transport within the progenitor may be able to trigger smaller amounts of pre-SN mass loss (see, e.g., Quataert & Shiode 2012). Due to the large rate of energy generation by nuclear reactions, convective motions near the star's core become very energetic during advanced burning phases such as carbon, neon, oxygen, and silicon burning (e.g., Woosley et al. 2002). The convective motion generates gravity waves that carry a power L_wave that scales with the convective luminosity L_con as L_wave ∼ M_conL_con (Goldreich & Kumar 1990; Alvan et al. 2015), where M_con is the mean Mach number of the convective motion. The gravity waves propagate outward and are partially reflected by overlying convective shells, though some wave energy is still transferred outward into acoustic waves after tunneling through these evanescent regions (see, e.g., Fuller et al. 2015; Takata 2016 for applications in low-mass red giant stars).

The acoustic waves dissipate via weak shock formation when the rate of energy transport by waves exceeds the maximum possible wave flux (Ro & Matzner 2017), which drops rapidly at the transition from the helium shell to the hydrogen envelope. The thermalized wave energy, depending on how fast it is injected, can trigger expansion of the envelope and also mass ejection (Fuller 2017). Recent studies of observational data have shown that pre-SN activity may be important for explaining the shape of the light curve (Moriya et al. 2017, 2018; Morozova et al. 2017, 2018), though too much wave heat (or heat deposited over too long a duration) is inconsistent with typical Type II-P SNe (Ouchi & Maeda 2019). Instead, more impulsive wave heating is favored based on light-curve modeling of SN 2017eaw, confirming that such energy deposition can better reproduce SN light curves (Morozova et al. 2020).

We note that wave-driven outbursts may not occur in many SN progenitors, as demonstrated by SN progenitor monitoring (Kochanek et al. 2017; Johnson et al. 2018) that rules out luminous and long-lasting outbursts in several nearby Type II-P SNe. The purpose of this paper is not to determine the frequency or likelihood of wave-driven outbursts, but rather their effect on the progenitor's structure if they do occur.

Fuller (2017) and Fuller & Ro (2018) compute one-dimensional stellar evolutionary models for a 15 M_⊙ star model with and without a hydrogen envelope. They demonstrate that wave heat has very different effects in hydrogen-rich and hydrogen-poor stars. In the latter, the energy can eject ∼10⁻²–10⁻¹ M_⊙ in an optically thick super-Eddington wind from the surface of the star. In hydrogen-rich stars, however, wave energy is thermalized at the base of the hydrogen-rich envelope and is likely insufficient to unbind the entire envelope. While small amounts of material may be unbound by a weak shock excited by wave heating, the one-dimensional models exhibit an inflation of the envelope as the extra thermal pressure blows a low-density "bubble" at the base of the hydrogen envelope.

However, the energy deposition can give rise to non-radial hydrodynamical instabilities that cannot be captured in any of the previous one-dimensional simulations. As mentioned above, wave energy deposition in a low-density bubble causes it to have a high thermal pressure but a low density. As it expands into the low-pressure and high-density envelope, Rayleigh–Taylor instability will occur at the interface between the bubble and the overlying envelope, strongly modifying the subsequent dynamics of the envelope's response. It therefore becomes interesting to examine multi-dimensional simulations of the star's response, which have never been performed in a systematic way for this process. In this work, we explore the general behavior of the envelope in response to such energy deposition and examine the aspherical effects on the pre-SN mass loss and the envelope structure.

The organization of the paper is as follows. In Section 2, we describe the hydrodynamical formalism for modeling the wave energy deposition in the stellar envelope and how we construct the initial model. Then in Section 3 we present a benchmark model that demonstrates the typical effects of wave heating in multiple dimensions. In Section 4, we examine how the mass loss process and the envelope structure change with the energy deposition and its duration. In Section 5, we discuss implications for supernovae and progenitor stars with different rates and durations of wave heating. In the Appendix, we further show how our simulations depend some of the input parameters, including the mesh resolution, the initial perturbation, and the boundaries.

To construct the initial models for the hydrodynamics run, we first use the one-dimensional stellar evolutionary code MESA (Modules for the Experiments in Stellar Astrophysics, Paxton et al. 2011, 2013, 2015, 2017, 2019) version 8118. We evolve a 15 M_⊙ star with solar metallicity, starting at the zero-age main sequence. The simulation continues until the central density approaches 10¹⁰ g cm⁻³. Then we extract the stellar profile when the star is carrying out core O-burning and map it to our two-dimensional mesh.

In the top left panel of Figure 1 we plot a Hertzsprung–Russell diagram of the model's evolution, and in the top right panel we plot the time evolution of the nuclear, neutrino, and surface luminosities of the stellar model. The surface luminosity does not vary during the late-stage vigorous nuclear burning in the core, but the neutrino luminosity in general follows the nuclear luminosity. Note that the shell burning of C and Ne triggers a brief jump in the nuclear luminosity, while O-burning provides a more steady energy production. In the bottom right panel of Figure 1, we plot the density and temperature profile of the MESA profile. We model only the outer part beyond 10 R_⊙ in the simulations, corresponding to m(r)/M ∼ 0.5 at the base of the H-envelope where the wave heating is expected to take place.

Zoom In Zoom Out Reset image size

We use a two-dimensional hydrodynamics code that we previously used extensively for SN simulations. The code is designed for multiple purposes and has been used in simulations of Type Ia SNe (Leung et al. 2015; Leung & Nomoto 2018, 2020), electron capture SNe (Leung & Nomoto 2017, 2019; Zha et al. 2019b; Leung et al. 2020) and accretion-induced collapse (Leung et al. 2019b; Zha et al. 2019a). The code utilizes the fifth-order shock-capturing weighted essentially non-oscillatory (WENO) scheme (Barth & Deconinck 1999) and the five-step third-order non-strong-stability-preserving Runge–Kutta scheme for time discretization (Wang & Spiteri 2007) of the Euler equations. In this work, we consider the equatorial plane of the star in polar coordinates (ρ, θ). This allows the use of a reflective boundary in the radial direction and a periodic boundary in the θ-direction. To perform the simulations in a computationally feasible time range, we run the simulations in a single quadrant.

To close the equations, we use the Helmholtz equation of state (Timmes et al. 2000). This consists of contributions from electron gas with arbitrary degeneracy, ions in the form of an ideal gas, photons in the Planck distribution, and electron–positron pairs. We use a reflective boundary for the radial inner boundary, and a free-flow boundary for the radial outer one. When the energy deposition is strong, the fluid near the inner boundary can expand significantly, causing mass flow outward. In that case, we provide additional mass to the innermost grid cells such that the outflow mass is balanced by the additional mass, as it would be in a real star by upwelling material. Then, we inject the corresponding gravitational energy to the system to maintain energy conservation and adjust the mass cut value. This ensures that the total mass and energy, including both the hydrodynamical grids and the mass cut, are conserved as a whole.

Below we describe the changes to the code in order to accommodate the wave heating physics. We consider a wave luminosity emerging from the core, ${L}_{\mathrm{wave},0}$, that propagates into our computational domain and then dissipates into heat. We model only the heat dissipation and not the waves themselves. In each time step, we compute how much energy is deposited at each of the computational mesh points by comparing the maximal luminosity ${L}_{\max }=2\pi {r}^{2}\rho {c}_{s}^{3}$ and the local L_wave, which decreases with radius as wave energy is converted into heat. When L_max > L_wave, no change in L_wave is made. Otherwise, the new L_wave is computed by solving (see, e.g., Landau & Lifshitz 1959; Ulmschneider 1970; Mihalas & Mihalas 1984; Fuller & Ro 2018 for derivations)

$$\begin{eqnarray}&&\displaystyle \frac{{{dL}}_{\mathrm{wave}}}{{dm}}=\displaystyle \frac{\gamma +1}{3\pi }\omega {c}_{s}^{2}{\left(\displaystyle \frac{{L}_{\mathrm{wave}}}{{L}_{\max }}\right)}^{3/2}.\end{eqnarray} \tag{ 1 }$$

This formula is appropriate for waves that gradually steepen to form weak shocks (i.e., the post-shock pressure only slightly exceeds the pre-shock pressure), which occurs when L_wave is slightly larger than L_max. As discussed in Fuller & Ro (2018), the weak shock damping prevents the waves from steepening into strong shocks with L_wave ≫ L_max, so the approximation of Equation (1) is good throughout most of the wave energy deposition region.

The relative change ΔL_wave/L_wave can be very steep for a high value of L_wave. We smooth the energy deposition by hand to a mass with a minimum of ${m}_{\mathrm{dep},\min }=0.1\,{M}_{\odot }$, which ensures that the process of energy deposition is resolved at the current resolution. To do this, the energy along each ray with a given θ is smoothed over a mass ${m}_{\mathrm{dep},\theta }=0.1\,{M}_{\odot }/{N}_{\theta }$, where N_θ is the number of grid points along the θ-direction. In each step, we search along a constant θ-direction for the transition position (r₀, θ₀) where L_max = L_wave. Then we integrate outward to look for the outer radius r₁ so that the accumulated mass ${\int }_{{r}_{0}}^{{r}_{1}}\mathrm{vol}(r,{\theta }_{0})\rho (r,{\theta }_{0})={m}_{\mathrm{dep},\theta }$. The local energy deposition in each grid cell is the mass of the cell divided by ${m}_{\mathrm{dep},\theta }$.

After evolving the stellar model, we take one snapshot during core O-burning, and map the density and pressure into our 2D simulation. The profile from the stellar evolution code obeys hydrostatic equilibrium. However, during mapping from the Lagrangian MESA mesh to the Eulerian simulation mesh, slight departures from hydrostatic equilibrium can be introduced. For the inner boundary, we chose a mass cut of mass ∼5.5 M_⊙ with a radius $\sim 7.5\times {10}^{11}$ cm, near the base of the hydrogen envelope and below the initial wave heating region. We use a uniform 1500 × 30 grid with a radial grid size of 1.5 × 10¹¹ cm and angular grid size of 3°.

In Appendix A, we further study how our model achieves hydrostatic equilibrium after it is mapped to our hydrodynamics model. While radial motion and changes to the density profile do occur, they are smaller than those caused by the wave heating that we examine below. Additionally, convection is not present in the initial model, which may underestimate the energy transport at early times if convective flows can help channel the deposited energy away.

The initial profile is spherically symmetric, so we add density perturbations in the form $\delta \rho (r,\theta )=\alpha \rho \sin (12\theta )\sin ({br}/R)$. We choose α = 0.01 in this work. In Appendix C we demonstrate that the exact value of α does not change our results. The density perturbation is added to mimic the inherent density fluctuations due to convection, and aims at seeding the growth of non-radial instabilities. The value b is chosen such that the total mass is conserved globally and along each radial direction. This prevents global motion of the star in any specific direction.

When the energy deposition starts, the heated region is usually within a narrow band in radius near the bottom of the simulation domain. As outlined in the previous section, the deposition region is smoothed along 0.1 M_⊙ starting from the innermost grid point with energy deposition. For our benchmark model the heating rate is set to 3 × 10⁶ L_⊙ and lasts for 120 days. In Figure 2, we show density color plots of the benchmark model at different time slices to show how the star expands and contorts due to energy deposition. The energy deposition creates a layered structure, due to weak shocks propagating through the envelope, as seen on day 114. Like the 1D models, the 2D models exhibit a "bubble" of lower density within the envelope, but with a weaker density inversion. At day 114, the bubble surface is near 6 × 10¹³ cm but expands outward, and by day 228, the bubble surface is wavy due to the density perturbations introduced in our initial conditions. The heating ceases after day 120, but the star does not expand significantly until after day ∼200, when the heat-induced pressure wave reaches the stellar surface. The core also develops non-radial structures due to non-radial instabilities and flows driven by the heating.

Zoom In Zoom Out Reset image size

We plot in Figure 3 the angle-averaged density profiles of the benchmark model at selected times to further analyze the density and velocity evolution of the ejecta. The energy deposition takes place near the core at ∼10¹³ cm, driving a density inversion that propagates outward. The density difference can be a factor of ∼3 between the minimum value in the trough and the peak. The angle-averaged density inversion gradually smooths out as the star expands, due to the growth of the Rayleigh–Taylor instability from the seed perturbations added to the model, which create the wavy patterns seen in Figure 2.

Zoom In Zoom Out Reset image size

In the left panel of Figure 4, we plot the angle-averaged radial velocity profile of the benchmark model. The energy deposition occurring in the first 120 days drives expansion of the envelope, with internal velocities ∼15 km s⁻¹. When this velocity/pressure wave reaches the surface, where the density is small, it steepens into a shock to trigger a rapid expansion of the envelope. The velocity can be as high as ∼40 km s⁻¹, near the escape speed of the star's surface. As the envelope expands, it slows down due to gravity, approaching ∼30 km s⁻¹ at 570 days after the start of energy deposition.

Zoom In Zoom Out Reset image size

To further understand the asphericity of the ejecta, we also plot in the bottom panel of Figure 4 the angle-averaged angular velocity snapshots. Most of the star has a smaller angular velocity than the radial velocity, but the angular velocity reaches ∼20 km s⁻¹ within the innermost 2 × 10¹³ cm. The direction of the angular flow changes with both radius and time due to its turbulent nature. Beyond ∼5 × 10¹³ cm, there is no observable motion in the angular velocity.

In this model, there is no direct mass ejection. The outer boundary at 2.3 × 10¹⁴ cm is sufficiently large to contain all matter during the simulation. Only the outer layers of the star achieve velocities close to the escape velocity. While pressure gradients may still provide extra acceleration, the sub-escape velocities of inner layers suggest that the net mass ejection will be low. Besides, after the O-burning has started, it takes ∼1 yr before the final collapse takes place, so there might not be sufficient time for any matter to escape before the final collapse.

Finally, we examine the energy budget of the benchmark model. In Figure 5, we plot the different components of energy as a function of time. The total energy grows linearly in the first 120 days due to the deposited energy in the model. Note also the small negative energy inflow due to the gravitational potential energy of mass added to the system as described in Section 2. Examining the trends in gravitational energy and internal energy, we see that the star does not significantly expand in the first 100 days, but then expands noticeably, lowering the gravitational energy. The early growth and subsequent decline of internal energy show that the system first stores the deposited energy in thermal energy, and later converts this into kinetic energy, which is then converted into gravitational energy as the star expands outward. We also note that the internal energy eventually decreases below its initial value, reflecting the nearly adiabatic losses during envelope expansion.

Zoom In Zoom Out Reset image size

We first study how models with the same total deposited energy behave, from which we can extract some general observations about the effects of wave energy deposition.

4.1.1. Hydrodynamics of LEMs

In Figure 6, we plot in the upper (lower) panel the density profiles of Model 1 (3) taken from the LEMs. In Model 1, which has a lower L_wave, the star shows a regulated expansion. A density inversion (like that shown in the benchmark model) develops at 8 × 10¹³ cm at day 114 and slowly propagates to 2 × 10¹⁴ cm by day 570. The density contrast across the inversion nearly stays fixed at roughly a few times the minimum density, while the sharp density gradient at the surface remains. In Model 3, the expansion proceeds steadily without forming a density inversion. A strong shock does not form in either case, leaving a sharp density gradient at the surface, and all material in the star remains gravitationally bound.

Zoom In Zoom Out Reset image size

Figure 7 shows the radial velocity profiles of the LEMs. In Model 1, the initial energy deposition triggers a pulse-like velocity profile at day 114 with a peak velocity of ∼30 km s⁻¹ at day 228. The whole star is gravitationally bound, including the shock-heated surface. Model 3 shows an erratic velocity profile characterized by shocks at early times, which can be seen at 3, 5, 7, and 9 × 10¹³ cm at day 114. After the shocks reach the surface, the star exhibits nearly "homologous" expansion in the sense that the expansion velocity is roughly proportional to radius. Again the velocity is too low to escape the star directly. The surface velocity is only 15 km s⁻¹ at the end of the simulation.

Zoom In Zoom Out Reset image size

These two models suggest that the actual changes in the density and velocity profiles are almost negligible in the low-energy regime. The heat triggers moderate expansion but no major restructuring of the density profile.

4.1.2. Hydrodynamics of MEMs

In Figure 8 we plot the density profiles of two models with the same energy input as the benchmark model, with E_dep = 1.195 × 10⁴⁷ erg, but at a low heating rate of 10⁶ L_⊙ for 360 days (Model 2) and at a high heating rate of 10⁷ L_⊙ for 40 days (Model 6). The lower energy deposition in Model 2 creates a density inversion that mildly accelerates with time as the density contrast approaches the surface. The density difference can be a few times higher than the minimum density at the trough. The star smoothly expands to ∼1.9 × 10¹⁴ cm. In Model 6, no significant density inversion develops in the whole simulation, just a small one near the surface. The stellar density profile is mostly monotonically decreasing and expanding, signifying a smooth expansion. The outer layers stagnate at the end of the simulation, suggesting that the star's expansion slows down, with all of its matter remaining bound.

Zoom In Zoom Out Reset image size

In Figure 9, we plot the velocity profiles of the MEMs. Similar to the density profiles, the velocity profiles also show interesting differences between the two models. Model 2 with a low heating rate does not show deceleration, and the star reaches nearly homologous expansion at the end of simulation. Model 6 shows almost the opposite behavior. The energy injection creates a two-peak shock structure at day 114, with the inner shock traveling faster than the outer. The shocks merge into a stronger shock at day 228, which reaches above 60 km s⁻¹, almost 50% higher than Model 2, despite the same total energy deposited. However, the expansion is not sustained by further heating, and the surface velocity drops by half and reaches ∼30 km s⁻¹ at the end of the simulation.

Zoom In Zoom Out Reset image size

4.1.3. Hydrodynamics of HEMs

We next examine the HEMs, Models 5 and 7, which have a higher deposited energy of 3.585 × 10⁴⁷ erg. Figure 10 shows that the density profiles are drastically altered by the large energy deposition. Both models develop large density inversions as material in the inner envelope is lifted outward. In Model 5, which has a lower L_wave, the density profile is flattened, with a density peak near 6 × 10¹³ cm at day 114 that propagates to nearly 2 × 10¹⁴ cm by day 570. In Model 7, the higher heating rate drives a stronger shock that expels material to larger radii. At day 456, the outer layers of the star reach the outer boundary of the simulation at ∼2.3 × 10¹⁴ cm.

Zoom In Zoom Out Reset image size

In Figure 11, we plot the velocity profiles of Models 5 and 7 from the HEM series. Model 5 shows that the star can still expand rapidly with a moderate L_dep but with a long duration. The velocity profile features a strong two-peak structure at intermediate times but approaches homologous expansion at late times. The outer layers asymptotically approach a velocity ∼40 km s⁻¹ larger than the escape speed. The inner layers slow down with time, developing a complex turbulent structure in the core.

Zoom In Zoom Out Reset image size

Dynamical features accompany the expansion in the velocity profiles of Model 7. When the strong shock breaks out from the surface of the star around day 228, the velocity peaks at 70 km s⁻¹, compared to the local escape velocity of ∼40 km s⁻¹. Afterwards, the ejecta slows down but can still escape. Turbulent motion persists in the inner part of the star. The outer layers of the ejecta reach the outer boundary of the simulation as shown in the profile at day 456. At the end of the run, the star approaches homologous expansion.

In this section, we discuss general differences between models with the same wave heating rate but different durations, or vice versa. The duration is determined by the lifetime of the core burning. The duration of the advanced burning, such as O-burning, generally decreases with stellar mass. The exact heating rate is also unclear because there is significant uncertainty in the wave luminosity generated by the convective core, but it is expected to increase for more massive stellar cores. Here, we explore the consequences of the variation of these values for the mass ejection process.

We examine how ${L}_{\mathrm{wave}}$ affects the evolution of the velocity profile by analyzing the three models with a heating duration of 120 days: Model 1 (LEM), Model 4 (MEM), and Model 7 (HEM). The energy deposition rate significantly affects how the star expands and how the star ejects mass. A lower L_wave makes the star expand but with no mass ejection. Higher values of L_wave can drive stronger shocks that lead to mass ejection from the surface. Only Model 7 with L_wave = 10⁷ L_⊙ exhibits significant mass ejection due to wave heating. The larger energy deposition of the HEM not only makes the star expand faster, but it also creates a prolonged trough in the density profile inside the star. In contrast, the LEM exhibits very little internal motion or mass ejection in its velocity profile.

We identify significant differences in the density and velocity structures between models with different heating rates. Models with low but sustained heating rates exhibit more steady expansion, with a small density inversion forming as the outer layers are gradually lifted outward by the expanding inner layers. Due to the low heating rate, no strong shock wave is excited, and matter is not shock-accelerated near the surface. We refer to this as a "bubble" phenomenon. In contrast, a larger heating rate excites a stronger shock that accelerates the surface layers outward, creating a low-density tail of material above the previous location of the photosphere, with a much smaller density inversion interior to this tail. We refer to this as a "bomb" phenomenon. Thus, low heating rates create a "bubble" phenomenon, while high heating rates (i.e., a heating timescale shorter than the local dynamical timescale) create a "bomb" phenomenon.

In Figure 12, we plot the total ejecta mass against time for the models presented. The ejecta mass is defined by the integrated mass elements when the sum of the local kinetic energy, internal energy, and potential energy is positive. This quantity is therefore time-dependent as demonstrated in previous sections. The seven models can be characterized into three groups based on the total energy deposition: LEMs exhibit little to no mass loss, MEMs exhibit moderate amounts of mass loss, and HEMs exhibit large amounts of mass loss. For LEMs (Models 1 and 3), there is effectively no ejected mass because the escape mass is as low as ∼10⁻² M_⊙. The MEMs (Models 2, 4, and 6) show an ejected mass that is roughly an order of magnitude higher at 0.5–1.5 M_⊙. Finally, in HEMs (Models 5 and 7), the ejected mass reaches asymptotic values of 4–6 M_⊙, where the whole envelope mass is also ∼6 M_⊙.

Zoom In Zoom Out Reset image size

The heating rate also affects the evolution of the ejected mass. Unsurprisingly, models with the largest heating rates (Models 6 and 7) show the fastest initial increase in unbound mass, while models with the lower heating rates require a longer time for any mass to accumulate enough energy to escape. Between models with the same energy deposition, lower heating rates with longer durations generally produce larger ejecta masses.

Here we compare the positions where the acoustic waves deposit energy. In Figure 13, we plot the position of the energy deposition against time for the models presented. The sharp drop in all curves occurs when heating abruptly shuts off. We see that models with a lower L_dep have their energy deposition at larger radii, which vary from ∼5 × 10¹³ cm in Models 1 and 2 down to 2 × 10¹² cm for Model 7. As the acoustic waves deposit their energy only when the maximum wave energy flux is exceeded (Equation (1)), a lower energy deposition rate requires a lower density and sound speed to match this criterion, so the heating occurs at larger radii. The roughly constant heating radius in Models 1–5 before day 100 suggests that the wave energy does not strongly alter the stellar structure in the wave heating region during this time. The trends in the heating radius are correlated with L_dep, because the wave heat causes the inner envelope to expand and decrease in density, causing the heating radius to move inward. This happens almost immediately for the largest heating rates. As seen in Models 5 and 7, the recession of the deposition sphere can reduce its size by almost an order of magnitude.

Zoom In Zoom Out Reset image size

At last, we compare the total energy as a function of time among all models and plot this in Figure 14. Unsurprisingly, the final energy is mostly determined by the total acoustic wave energy deposition. Note that only Models 5 and 7 finish with positive total energy, implying large amounts of unbound material, as verified by the results in Figure 12. However, as the star expands, extra contributions of energy emerge due to the energy content of inflowing/outflowing matter across the simulation boundaries, so the final energy in the simulation is not merely the sum of the initial energy and the heating energy.

Zoom In Zoom Out Reset image size

This is not just an artifact of the simulations, because mass in a real star will generally flow from below the wave heating region to above it, changing the net energy content of the envelope. This effect is negligible for models with lower L_dep, due to a smaller contribution from material flowing across the boundaries. Models with larger heating rates drive material away from the inner boundary, causing more matter to flow into the simulation domain from the inner boundary, which advects negative gravitational energy into the simulation domain.

In Model 5, there is a mild increase of energy between days 400 and 450. This is the period when the rapid expansion creates an energy discontinuity near the isothermal outer layer, which is set to be at the minimum temperature allowed in the simulation (∼3000 K). When the discontinuity passes across this region, fluid elements can be cooled below the minimum temperature. In this code we choose to stick with thermodynamic consistency instead of energy consistency, so the missing energy between the actual temperature and the temperature floor is added into the system. Since this only alters the total energy budget by about 1%, we do not believe it qualitatively affects our conclusions.

In general, the outburst process depends strongly on the wave excitation in the core, in addition to the properties of the envelope, and both core and envelope should be evolved together for robust estimates. Fortunately, most red supergiant envelopes have very similar, nearly polytropic structures. The envelope masses are also similar, ∼7 M_⊙, since more massive stars have more of their mass in their cores, and lose more envelope mass by winds. In contrast, stars of different mass have very different wave heating rates and heating durations (Shiode & Quataert 2014). Hence, in our parameter study, lower values of L_wave for longer durations are more representative of low-mass stars, while higher values of L_wave for shorter durations are more representative of high-mass stars. Different burning phases behave similarly: carbon shell burning can be much longer (∼10² yr) but the wave luminosity is correspondingly smaller. Silicon burning generates much larger wave fluxes but only lasts days, so the envelope does not have time to respond (Fuller et al. 2015).

It is possible that other processes occur in some stars, such as convective shell mergers. These would spur intense nuclear burning and vigorous convection (Meakin & Arnett 2006; Collins et al. 2018; Yadav et al. 2020), driving very large wave fluxes into the envelope. Another scenario is that the core dynamo amplifies the magnetic field in the convective shells near the core, which could also generate waves and lead to energy deposition near the surface (Soker & Gilkis 2017). While we do not investigate such processes here, an uncommon subset of stars could generate wave fluxes much greater than our fiducial estimates.

Our simulations do not include radiative transfer and cannot predict the detailed appearance of wave-driven pre-SN outbursts. From the 1D models of Fuller (2017), the photometric signal of a pre-SN outburst in a red supergiant is likely dominated by shocks propagating up to the photosphere after short but intense periods of wave heating, such as the core neon burning phase. Our 2D simulations with large heating rates also launch these shocks, which are not strongly affected by non-radial instabilities. Hence, we expect that the observable features of pre-SN outbursts are captured reasonably well by 1D models. When the outburst takes place in a binary system, the accretion of the ejected matter by the companion star might also produce a bright event (Mcley & Soker 2014; Danieli & Soker 2019).

In this work, we have extensively computed two-dimensional models showing how the wave heating generates asphericity and drives mass loss. In Figure 15, we further demonstrate how the asphericity develops at late times. We examine Models 2 and 6, finding that both models are spherical to a good approximation. However, notable turbulence is driven within the middle of the envelope, which is especially evident at radii around 5 × 10¹³ cm. These structures are much more pronounced in Model 6, likely because the larger heating rate (but shorter duration) fuels more expansion of the inner envelope, generating stronger Rayleigh–Taylor instability. Furthermore, the density profile in Model 2 shows observable "wrinkles," though it is not yet clear if these features are real or numerical artifacts. For LEMs and MEMs, the asphericity is small, and the wave heating primarily contributes to a quasi-spherical expansion of the stellar envelope.

Zoom In Zoom Out Reset image size

For comparison, we also plot in Figure 16 the density color plots of two higher energy models, Models 5 and 7. The most striking features of these models is the large-scale structure with three wavelengths per quadrant at radii around 10¹⁴ cm. This is clearly a consequence of the initial configuration of the simulations, where a sinusoidal perturbation with three periods in the radial and angular directions (and an amplitude 1%) was embedded in the density profile. However, by the end of the simulation, this 1% perturbation has grown to a large amplitude, producing a variation by as much as a factor of 10 in the horizontal density fluctuations. We attribute this large amplification to rapidly growing Rayleigh–Taylor instability. Interestingly, the radial location where the instability is most prominent depends on the wave heating rate, but it is clear that higher total heat deposition allows for more growth of density perturbations. Of course, a real star would have initial density fluctuations at a variety of scales due to turbulent convection in the envelope, but these models indicate that those perturbations can be magnified by wave heating effects, potentially producing significant asphericity in the exploding star.

Zoom In Zoom Out Reset image size

To quantify the level of asphericity in the models, we compute the mass-weighted asphericity via the Fourier modes

$$\begin{eqnarray}&&{S}_{n}=\displaystyle \frac{1}{M}{\int }_{{R}_{\mathrm{cut}}}^{R}{\int }_{0}^{\pi /2}8{r}^{2}{drd}\theta \,\rho \sin (n\theta ).\end{eqnarray} \tag{ 2 }$$

Here, the integer n is the Fourier harmonic and M is the total mass of the star. We also compute a similar expression with C_n using $\cos (n\theta )$ as the Fourier component, and we thus define the asphericity of nth order, ${A}_{n}=\sqrt{{C}_{n}^{2}+{S}_{n}^{2}}$.

In Figure 17 we plot the time evolution of the asphericity modes A₁ and A₃ of all seven models presented in previous sections. The n = 3 component is often dominant as discussed above. Comparing the two plots, A₁ is indeed smaller than A₃ by a factor of a few. The LEMs show a mild increase in asphericity, but their low values show that the stars remain nearly spherical. The MEMs show a larger asphericity that grows nearly linearly in time. The HEMs have the most significant asphericity, with the A₃ component of Model 7 reaching 10%. Despite its lower total energy deposition, Model 5 exhibits the largest value of A₁ and hence has the most elongated structure, though the elongation only reaches ≈3%. The angular dependence of the envelope density may give rise to interesting effects for the subsequent SN and its appearance.

Zoom In Zoom Out Reset image size

In Appendix B we further compare models with different initial aspherical seeds to see how this affects the star's final asphericity, finding that the results are not strongly dependent on this choice.

First, our simulations use a moderate resolution that is almost constant throughout the star with Δx = 1.5 × 10¹⁰ cm. This Δx is rather large compared to the mass cut, which is ∼7.5 × 10¹⁰ cm. Thus, the inner features of the star might not be fully captured and may be under-resolved. This value of Δx can provide us the necessary resolution for the outer part to keep track of the smaller features due to hydrodynamical instabilities. Most such features should be produced in the outer envelope because that is the place where the density inversion is present.

Second, our simulation does not resolve the acoustic waves emitted from the core, but instead only includes a heating term to account for their energy deposition. Studying how these waves deposit heat and momentum is numerically challenging because it takes place on a short length scale and a short timescale, i.e., the wavelength and wave period. The huge number of acoustic waves leaked into the envelope provides the possibility of modeling this process statistically, namely by the recognition of L_wave as a time-averaged quantity. Interestingly, our results imply that the net effect is for acoustic waves to damp via weak shocks, and to convert their energy into pressure pulses with greater length scale and timescale, which then continue to propagate out until they shock again. It would be an interesting extension to resolve both the original acoustic waves and the secondary pressure pulses to see if the results differ from the treatment adopted here.

Third, we have used the Helmholtz equation of state for our computation. But for modeling the envelope, we are approaching a density close to 10⁻¹¹ g cm⁻³, which is the density floor provided by the equation-of-state table. As discussed above, the temperature would also decrease below 10⁴ K if we did not impose a temperature floor. This limits us in trying to explore further how the ejecta behaves beyond day 570. In some of our more energetic models, vacuum zones develop inside the star, with a few grid cells along the same r reaching the density floor. These are not evident in the figures because the density profiles shown are angle-averaged. The large associated pressure differences begin to trigger unphysical flow. Further extrapolation to lower densities and temperatures will require an updated equation of state. Additionally, at these low densities, the matter might not be tightly coupled to the radiation field, so non-diffusive radiative transfer may be important. Furthermore, the ionization state can be out of thermodynamic equilibrium. These factors might largely complicate the computations and this extra physics will be interesting follow-up work.

The density floor is set to maintain code stability when calling the equation of state. However, the presence of such a floor also sets a minimum mass that can be resolved in our simulation. To estimate that, we consider the outermost grid cells used in our simulation, located at a radius ∼1.5 × 10¹⁴ cm. Using the density floor, the minimum mass traceable by this grid is ∼6 × 10⁻⁴ M_⊙. The corresponding masses for inner grids are smaller because the grid volume scales as ∼r². This means that the ejected mass ranging from ∼10⁻² to 6 M_⊙ in our models is likely well resolved by our code. For ejecta masses ≲0.01 M_⊙, a suitable extension to even lower density below 10⁻¹¹ g cm⁻³ will be necessary. However, extending to lower densities will also require additional physical considerations, such as the validity of thermodynamic equilibrium and the inclusion of radiative processes.

Fourth, our two-dimensional simulations might underestimate small-scale turbulent features. Hydrodynamical instabilities, such as the Kelvin–Helmholtz instabilities and Rayleigh–Taylor instabilities, can extend to the Kolmogorov length scale, where the viscous nature of the fluid becomes prominent. In fact, for a three-dimensional model, the development of Rayleigh–Taylor instability can be different because the extra dimension allows an extra degree of freedom for the development of the smaller features. This phenomenon is frequently seen in large-scale eddy modeling, where two-dimensional and three-dimensional turbulent energy cascades have different directions of energy flow. However, the large-scale features reported here should be similar in the two types of simulations.

Fifth, we do not consider the effects of recombination and radiative transfer in these simulations. In our progenitor model, the outer ≈0.1 M_⊙ has an optical depth τ ≲ 10, hence the assumption of high optical depth in those layers is not very good. While the approximation of radiation being in equilibrium with matter is well maintained in the rest of the star, the true ejecta mass of our less energetic models (with M_ej ≲ 0.1 M_⊙) may be reduced if radiation can leak out before their outer layers are accelerated toward the escape speed.

The effects of recombination can also potentially affect the expansion of the envelope. Our equation of state does not include recombination energy, and the ionization temperatures of ⁴He (∼29,000 K) and ¹H (∼10,000 K) overlap with the temperature range used in the simulation, meaning that partial recombination should take place during its expansion. Because of the envelope's low binding energy, this energy contributes substantially to the energy budget, and it may drive greater envelope expansion and mass loss than predicted by these simulations. On the other hand, recombination reduces the envelope opacity, allowing thermal energy to leak out, as described above. Resolving this competition between energy deposition by recombination and energy loss by radiation will be necessary for better estimates of the mass loss and pre-SN density profile.

Sixth, in the simulations we have made approximations for the initial conditions and energy deposition. In fact, how the exact energy deposition takes place in such a thin mass slice, and how the initial convective structure can be coupled to a hydrostatic model, should be further explored.

To clarify how these limitations affect the final results, we provide more comparison tests in the appendices. In Appendix A, we study how the initial density perturbation affects the later evolution without energy deposition. In Appendix B, we examine how the choice of density perturbation affects the later evolution. In Appendix C, we further compare how our two-dimensional energy deposition differs from the one-dimensional energy deposition model previously reported in Fuller (2017). In Appendix D, we study how our choice of smearing mass for the energy deposition affects the morphology and dynamics of the ejecta. In Appendix E, we provide further details on our coordinate system. The above tests aim at clarifying the physics necessary for a comprehensive parameter survey for our future work.

In this work, we study the multi-dimensional hydrodynamical response of a red supergiant envelope heated by wave energy transport from nuclear burning in the core. We extract the stellar envelope from a 15 M_⊙ stellar model, and we inject energy similar to that expected from the wave luminosity and duration of the core oxygen burning phase. To understand sensitivity to the uncertain wave luminosity and duration, we treat the heat deposition rate and its duration as model parameters. Unsurprisingly, we find that heat deposition in excess of the envelope's binding energy can drive large amounts of mass loss, while heat deposition lower than this threshold has a much smaller impact.

We also observe two classes of behaviors that depend on the wave heating rate. For a lower wave energy flux, the wave energy causes gradual expansion of the stellar envelope without strong mass loss, and the expansion is mostly spherical. With a high wave energy flux, the wave energy can amplify initial density asphericities and trigger uneven heating. This results in more asymmetric motion in the envelope, where the Rayleigh–Taylor instability arises, smoothing radial density inversions that occur in one-dimensional models. Larger heating rates also drive stronger shocks, expelling material above the original photosphere of the star. These low-density tails above the star's photosphere may have important effects on the early light curves of Type II-P SNe. Such events would also likely lead to observable pre-SN outbursts.

Comparison between our two-dimensional models and spherically symmetric models shows only modest differences, mostly arising from Rayleigh–Taylor instabilities that smooth out density inversions created by wave heating. Hence, one-dimensional models are likely to yield reasonably good estimates of the unbound mass and the photometric signals produced by pre-SN outbursts in red supergiants. However, multi-dimensional models or one-dimensional models incorporating the effects of Rayleigh–Taylor instability (Duffell 2016) should be used for estimates of the pre-SN density profile.

We further discuss how the results of this work can be applied to stars of different mass. More massive stars with more vigorous convection during oxygen burning can potentially result in large enough heating rates to launch shocks that trigger observable outbursts and mass loss from the stellar envelope before the star collapses. Less massive stars with lower wave heating rates are unlikely to suffer such outbursts, except in cases where degenerate burning flashes, convective shell mergers, or a different core structure permit larger wave fluxes into the envelope. The corresponding optical signal and the implications for SNe light curves will be interesting follow-up work.