© ESO 2021

1. Introduction

Growing observational evidence suggests that massive stars sometimes experience episodic and intense mass loss accompanied by temporal brightening. Eta Carinae is one of the most well-studied and well-known objects that experienced such an event (e.g., Davidson & Humphreys 1997). This intense mass loss or brightening event has been considered to be related with the activity of luminous blue variables (LBVs), which were introduced by Conti et al. (1984). On the other hand, some recent observations suggest that Wolf-Rayet stars (WR stars) may also experience such events (Foley et al. 2007; Pastorello et al. 2008; Smith et al. 2020).

These mass loss events could affect the observational features of subsequent core-collapse supernovae (SN). After a massive star abruptly loses a significant portion of the envelope, dense circumstellar matter (CSM) is formed around the star. If a core-collapse SN (CCSN) takes place in this dense CSM, the ejecta collide with the CSM. Then the kinetic energy of the ejecta is dissipated at shocks and becomes the main energy source (e.g., Chugai 1992; Smith et al. 2017). These kinds of SNe are classified as Type IIn supernovae (SNe IIn) in case of hydrogen-rich CSM (Schlegel 1990) or Type Ibn supernovae (SNe Ibn) in case of helium-rich CSM (Pastorello et al. 2007).

Recent observations have revealed that some progenitors of SNe IIn experienced temporary brightening phase a few years or a few decades before the core collapse. For example, SN 2018cnf (Pastorello et al. 2019a), SN 2016bdu (Pastorello et al. 2018), SN 2013gc (Reguitti et al. 2019), PTF12cxj (Ofek et al. 2014), and SN 2011ht (Fraser et al. 2013a) were reported to exhibit such brightening. Since the peak luminosity exceeds the Eddington luminosity of a massive star, this brightening must lead to an eruptive mass loss. The sparsely observed light curves prior to these SNe show that the brightening phase typically lasted for several years and indicate that the progenitors repeated episodic mass loss events during this period. SN 2009ip is one of the most famous SNe IIn whose progenitor star experienced multiple brightening phases likely associated with episodic mass loss events. The progenitor star of SN 2009ip in a brightening phase was first detected in 2009, and repeatedly exhibited brightening with short intervals less than 50 days in 2011 (Pastorello et al. 2013). Eventually this object experienced the most luminous outburst “2012b” in 2012. While some authors have argued that it was a genuine CCSN (e.g., Mauerhan et al. 2013; Smith et al. 2014; Graham et al. 2017), some have suggested other scenarios (e.g., merger-burst event, Soker & Kashi 2013; pulsational pair-instability event, Pastorello et al. 2013). The progenitor mass of SN 2009ip is estimated as 50 − 80 M_⊙ (Smith et al. 2010a).

The mass loss rates from progenitors of SNe IIn are estimated at 0.026 − 0.12 M_⊙ yr⁻¹ (Kiewe et al. 2012), 10⁻⁴ − 10⁻² M_⊙ yr⁻¹ (Taddia et al. 2013), or more than 10⁻³ M_⊙ (Moriya et al. 2014). These values are so high that they cannot be reconciled with a steady wind mass loss model (Vink et al. 2001; van Loon et al. 2005; Smith 2014). In addition, some SNe IIn show bumps in their light curves that are thought to be related to bumpy density structures of CSM (Reguitti et al. 2019; Nyholm et al. 2017; Stritzinger et al. 2012). Although these bumps seem to be rare cases among SNe IIn (Nyholm et al. 2020), this bumpy CSM structure also implies not steady but episodic mass loss events from the progenitor star. Therefore, a dynamical phenomenon, which is not included in most of current stellar evolution models, would occur during the temporary brightening phase.

However, the extra energy sources that trigger these mass loss events and dynamical mechanism of the mass eruption reproducing observational features have not been fully understood. There are several works that investigate the dynamics of the stellar envelope of massive stars into which an extra energy is deposited (Quataert et al. 2016; Fuller 2017; Fuller & Ro 2018; Ouchi & Maeda 2019; Dessart et al. 2010; Owocki et al. 2019; Kuriyama & Shigeyama 2020). In Kuriyama & Shigeyama (2020), we investigated the observational properties of mass eruption from the progenitors of SNe IIn/Ibn assuming an extra energy is supplied into the envelope of single massive stars. We deposited energy into the bottom of the stellar envelope with a short timescale and calculated the dynamics of the subsequent mass eruption using 1D radiation hydrodynamical simulation code.

Although energy is deposited only once and the corresponding single dynamical eruption event is discussed in these studies including Kuriyama & Shigeyama (2020), an intense mass loss event is often repeated in real situations, as discussed above, and this could become a key factor for understanding the mass loss mechanism. Once a dynamical eruption takes place, the expanding envelope keeps its altered density profile for the thermal timescale. If another eruption event occurs again before the envelope returns to a hydrostatic equilibrium state, the property of the second eruption may be completely different because of a different density profile of the envelope.

To deal with this problem, in this work we study the dynamics of eruptive mass loss that is repeated twice, as a supplemental work of Kuriyama & Shigeyama (2020). Of course eruption can repeat more than twice, but we just focus on a comparison between eruptions from the original (hydrostatic) envelope and the expanding envelope. We assume that eruptive mass loss is related to the local peak of the energy generation rate of nuclear burning shown in Fig. 1, although the specific mechanism of energy transportation from the burning region to the envelope is not assumed as was the case in Kuriyama & Shigeyama (2020). We deposit extra energy at the bottom of the hydrogen envelope twice and investigate the properties of each eruption and difference between the first and second eruptions. The second energy injection is conducted before the stellar envelope relaxes from the first injection. Each extra energy is injected for a short period of time and the resulting dynamical eruption is investigated by radiation hydrodynamical simulations. In Sect. 2, we introduce the progenitor models, which we use as initial models in our simulation and the method of radiation hydrodynamical calculation. In Sect. 3 we present the results of calculation and find out the differences between the first and second eruptions. We clarify the implications to the observations in Sect. 4 and present our conclusions in Sect. 5.

Fig. 1. Evolution of the energy generation rate by nuclear burning for each model indicated by labels. The rate becomes higher toward the core collapse. There are some local peaks for each model around ten years, five years, and one year before core collapse. We adopted the peaks at 11.2 years before core collapse (model RSG) and at 7.2 years before core collapse (BSG) as initial models for our calculations. These peaks correspond to the core neon burning phase.

2. Setup and methods

2.1. Progenitor models

Although it has been often considered that eruptive mass loss and SNe IIn are related to activities of LBVs (e.g., Kotak & Vink 2006; Gal-Yam et al. 2007; Langer 2012) as described in Sect. 1, observations suggest that red supergiants (RSGs) may also be progenitors of SNe IIn (Smith et al. 2009; Bilinski et al. 2015). The progenitor mass is also ambiguous and varies with the subclass of SNe IIn (Smith 2014). In this work, we adopted two types of 15 M_⊙ progenitors, namely the blue supergiant model (model BSG) and the red supergiant model (model RSG). We made these progenitor models using the stellar evolution code MESA (release 10398; Paxton et al. 2011, 2013, 2015, 2018, 2019). We did not include the effect of rotation in our calculations. We stopped the calculations of MESA 11.2 yr before core collapse for model RSG and 7.9 yr before core collapse for model BSG, respectively (Fig. 1). We remeshed these models to resolve the structure of radiation mediated shocks in the outer envelope and used them as initial models for our hydrodynamical simulations presented in Sect. 2.2. We set the numbers of grid points to 10 000 cells for model RSG and 5000 cells for model BSG. We tested other mass resolutions for model RSG (8000, 12 000, 14 000, and 16 000 cells), and we find that as the number of grid points becomes larger the amount of ejected mass by energy injection tends to increase. The increase converges above 14 000 cells, and there is a roughly ten percent difference in ejected mass between in case of 10 000 cells (the resolution adopted in this paper) and 16 000 cells. Therefore, model RSG has an ejected mass accuracy of ten percent at worst. Model BSG has more cells in the ejected matter than model RSG. Thus we think model BSG has an ejected mass accuracy comparable to or better than model RSG.

These two models have the same zero-age main sequence mass (15 M_⊙) but different metallicities (Table 1). The difference in metallicity causes different opacities in the stellar envelope. As a consequence, a star with higher metallicity evolves to a RSG and another star with lower metallicity becomes a BSG (different density distributions in Fig. 2; radii and effective temperatures in Table 1). The time evolution of the energy generation rate of nuclear burning is also slightly different between the two models. The timescale of each nuclear burning stage in model RSG is longer than that in model BSG (Fig. 1) because model RSG has a lower core mass (Table 1). Model RSG has experienced more intense steady mass loss because of higher metallicity and thus has a smaller core. We started hydrodynamical calculations at 11.2 yr before core collapse for models RSG and 7.9 yr for models BSG corresponding to the local peaks of nuclear burning luminosity, as mentioned in Sect. 2.2. Although we can identify other higher peaks within one year before core collapse, we only focus on these peaks because it could take more than one year for the ejected envelope to extending far enough to reproduce the observed typical distance from the progenitor star to a CSM (e.g., 160 AU, Smith & McCray 2007). Detailed methods of constructing the two progenitor models are described in Appendix A.

Fig. 2. Density distribution of each progenitor model. Model RSG has a more extended envelope than model BSG because of higher metallicity (Table 1) and opacity. In addition to that, model RSG has experienced stronger mass loss due to higher metallicity and thus has a lower total mass.

2.2. 1D radiation hydrodynamical simulation

As described in Sect. 1, the main purpose of this paper is to reveal differences in the properties between the first and second mass eruptions from SNe IIn progenitors. To investigate this topic, we used the same 1D radiation hydrodynamical simulation method as in Kuriyama & Shigeyama (2020). This code was originally developed by Shigeyama & Nomoto (1990) and was revised to solve the Riemann problem by the piecewise parabolic method (Colella & Woodward 1984); radiation transfer is treated by flux-limited diffusion approximation (Levermore & Pomraning 1981). We used Kramer’s opacity and did not include the effect of dust formation. Therefore, a calculation of a low temperature region in which the dust opacity becomes important should contain some uncertainty. We injected extra energy into the bottom of the hydrogen envelope twice without specifying from where and how the extra energy is supplied; the mechanism for supplying the extra energy has not been fully understood so far (see Sect. 1). The interval between the first and the second energy injection Δt was scaled with the dynamical timescales of the envelope t_dyn (98 days for RSG and 3.5 days for BSG). This is because we focus on the dynamics of the envelope, although the interval depends on the physics of energy source in reality. The eruption of the envelope proceeds on the timescale of t_dyn and after that it shrinks on the timescale of t_KH. In our model, t_KH (∼200 yr for model RSG and BSG) is much longer than the remaining time to core collapse, and thus only t_dyn seems to play a role in determining the physics of the multiple eruption. Therefore we used t_dyn as a fiducial interval between the two injections. As shown in Table 2, the first energy E_inj1 was injected at the time t = 0 and the second energy E_inj2 at t = Δt_inj. For each of the progenitor models RSG and BSG, five patterns of calculation with different amounts of the injected energy and Δt_inj were conducted (ten patterns of calculations in total). We call these models RSG1-f, s, m, l, RSG2-f, BSG1-f, s, m, l, and BSG2-f, respectively (Table 2). The letter f stands for “fiducial” and the other letters s, m, and l stand for short, medium, and long intervals, respectively. The quantity E_inj1 was roughly set to one-third or one-quarter of |E_envelope| (E_envelope represents the total energy of the envelope) for every model. This value is enough to expel 10⁻² − 10⁻¹ M_⊙ material (Kuriyama & Shigeyama 2020), which is a typical amount of CSM in SNe IIn.

3. Results

3.1. First energy injection and corresponding mass eruption

Two injections of extra energies at different epochs cause two distinct mass eruption events (Table 3) and two distinct peaks of luminosity (Figs. 3 and 4) and effective temperature (Figs. 5 and 6). About one- quarter (model RSG1-f, s, m, l, and RSG2-f) or one-third (BSG1-f, s, m, l, and BSG2-f) of |E_envelope| was injected into the bottom of the hydrogen-rich envelope in the first energy injection E_inj1. A relatively strong shock wave propagates toward the surface and breaks out with a rapidly rising light curve in all the models. The first energy injection ejects matter with masses of ∼0.01 M_⊙ (model RSG1-f, s, m, l and RSG2-f) and ∼0.2 M_⊙ (model BSG1-f, s, m, l and RSG2-f). The ejected mass in model RSG1-f, s, m, l and RSG2-f (model BSG1-f, s, m, l and BSG2-f) are equal in the first eruption because the same amount of energy E_inj1 is injected. Although the ratio of E_inj1/E_envelope for model RSG1-f and RSG2-f is similar to that for model BSG1-f and BSG2-f, there is an order of magnitude difference in the ejected mass above. This is because the envelope of model BSG is denser than that of model RSG and therefore a more energetic shock wave propagates outward and expels a considerable amount of the envelope. For every model, the hydrogen-rich envelope remains inflated after the first eruption as indicated from the significantly altered radius and density profiles in Figs. 7 and 8.

Fig. 3. Light curves for models RSG1-f, s, m, l, and RSG2-f. The first eruption triggered by the injected energy Einj1 produces a peak at day ∼80. The second peak in each model corresponds to the second eruption triggered by the injected energy Einj2.

Fig. 4. Same as Fig. 3 but for models BSG1-f, s, m, l, and BSG2-f.

Fig. 5. Evolution of effective temperature for each RSG progenitor model.

Fig. 6. Evolution of effective temperature for each BSG progenitor model.

Fig. 7. Radii as functions of the enclosed mass for the fiducial models. The three lines in each panel represent the profiles before the first energy injection (dotted line), between the first and the second energy injections (dashed line), and at the end of the hydrodynamical simulations (solid line).

Fig. 8. Same as Fig. 7 but for the density profiles.

3.2. Second energy injection and corresponding mass eruption

Before the expanding envelope shrinks and returns to a hydrostatic equilibrium state, the second energy injection was conducted. The amount of injected energy is the same as in the first injection for models RSG1-f, s, m, l and BSG1-f, s, m, l. On the other hand, for models RSG2-f and BSG2-f, twice (model RSG2-f) or 1.5 times (model BSG2-f) more energy is injected. These second energy injections trigger the second mass eruption events accompanied with the second peaks in luminosity. The matter ejected by the second energy injection collides and interacts with the expanding envelope or matter ejected by the first energy injection (dashed lines in Fig. 9; discontinuities around M_r ≃ 14.5 M_⊙ in each panel correspond to shock). After that the gravitationally bound part of expanding envelope begins to fall back, and this fall-back matter undergoes damped oscillations (Fig. 10). The trajectory of a fluid element follows that of a free-fall particle as seen from the red bold line. This also indicates that the gravity dominates the motion of the fall-back matter. Different periods of oscillation between each radius cause collisions (e.g., collision between M_r = 12.5 M_⊙ and 12.0 M_⊙ at t = 2.3 × 10⁶ s in Fig. 10). These kinds of collisions result in additional shock waves (e.g., solid line at M_r = 12.3 M_⊙ and 13.4 M_⊙ in bottom panel of Fig. 9) and discontinuity in density profiles (solid lines in each panel of Fig. 8).

Fig. 9. Velocity profiles as functions of the enclosed mass for the fiducial models. The three lines in each panel represent the profiles before the second energy injection (dotted line), after the second energy injection (dashed line), and at the end of the hydrodynamical simulations (solid line).

Fig. 10. Trajectories of some fluid elements for model BSG2-f. The red bold line represents the trajectory of a matter with a mass coordinate Mr = 12.5 M⊙ assuming free-fall motion from the maximum radius.

There are three key differences between the first and second mass eruptions. First, the brightening phase associated with the second mass eruption lasts more than ten times longer than that in the first mass eruption (Figs. 3 and 4). The luminosity gradually increases because photons can easily diffuse out from the shock wave in the inflated low-density envelope. After the peak, the luminosity slowly declines and the brightening phase lasts ∼100 days in model RSG1-f, s, m, l, 2 − f and ∼ten days in model BSG1-f, s, m, l, 2 − f. We describe what makes the different duration of brightening phases between two eruptions in Sect. 3.3.

Second, there is a significant difference in the color between the first and second mass eruption events (Figs. 5 and 6). For every model, when the first mass eruption takes place, its effective temperature rises by a factor of ∼ten. On the other hand, in the second mass eruption, the effective temperature rises but does not reach 10 000 K. The temperature at the local maximum in the second mass eruption is even lower than the effective temperature of the original state before the first energy injection (namely, the output value from MESA). Thus a progenitor star observed a few years before the SN does not necessarily follow the standard core mass luminosity relation constructed by the standard stellar evolution models that assume quasi-hydrostatic evolution.

Third, the amount of ejected mass is different between the first and second mass eruption even with the same injected energy (E_inj1 = E_inj2 for RSG1-f, s, m, l and BSG1-f, s, m, l). Interestingly, RSG1-f and BSG1-f show the opposite results. While the second energy injection ejects ∼2.5 times higher mass than the first injection in model BSG1-f, it ejects only one-tenth mass of the first injection in model RSG1-f. These results come from two antagonism effects, which are the weakly bound envelope after the first eruption (this effect enhances eruption) and shock wave attenuation due to diffusion of photons (this effect weakens eruption). While the latter effect works in model RSG1-f, the former effect works more effectively in model BSG1-f.

The amount of ejected mass is sensitive to energy input. In model RSG2-f, the second energy injection ejects ∼5 × 10² times higher mass than that in model RSG1-f, s, m, l even though the difference of E_inj2 between them is only a factor of two. This is because the velocity of matter where the shock wave locates (v_matter) and the escape velocity (v_esc) are very close to each other in cases in which a part of the envelope is ejected. A certain part of the envelope is marginally bound (v_matter ≃ v_esc) and a slight difference in energy input and resulting shock velocity determine whether it is bound or unbound. For example, when we compare second eruptions in RSG1-f and RSG2-f at t = 2.5 × 10⁷ s, v_matter ≃ 4 × 10⁶ cm s⁻¹ for RSG1-f and v_matter ≃ 6 × 10⁶ cm s⁻¹ for RSG2-f (dashed lines at top and second top panel in Fig. 9). On the other hand v_escape is around 5 − 6 × 10⁶ cm s⁻¹ in this region. Thus, the magnitude relations between v_matter and v_esc are different between RSG1-f and RSG2-f and it makes a large difference in ejected mass (Table 3). This kind of sensitivity implies that, although there are a wide diversity of observed SNe IIn or II-P in terms of the presence or the amount of CSM, this diversity may be caused by a smaller difference in the late-phase progenitor evolution and extra energy supply between different progenitors than has previously been thought.

We have seen the three key differences between the first and second mass eruptions. These differences hold at least for an interval in the range 0.5t_dyn <  Δt_inj <  4.0t_dyn and would hold over the wider range of Δt_inj because the envelope contracts on the Kelvin-Helmholtz timescale (∼200 yr for model RSG and BSG), which is much longer than the dynamical timescale. Meanwhile, when we compare the properties of the second eruptions between models that have different intervals of energy injection Δt_inj (−f, −s, −m and −l), a model with a longer interval tends to have a higher peak luminosity (Fig. 11). Correspondingly, the total radiated energy also has the same trend. On the other hand, models that have different intervals eject a similar amount of matter with a similar kinetic energy (Table 3). These behaviors originate from the velocity of decelerated matter once ejected in the first eruption. As the interval becomes longer, the ejected matter in the first eruption is gradually decelerated by the gravity of the star, while matter velocity in the downstream region of the outward-shock front is almost constant because the same amount of energy is injected. Thus, the difference in velocities between the upstream and downstream regions increases for longer intervals (e.g., for models RSG1-s, f, m, l, 1.9 × 10⁶, 2.7 × 10⁶, 4.3 × 10⁶, and 7.9 × 10⁶  km s⁻¹, respectively, at the luminosity peak), leads to higher pressure and temperature in the shocked region, and results in a higher peak luminosity and radiated energy. The peak luminosity and radiated energy in Fig. 11 seem to increase more gradually with increasing Δt_inj and would converge when Δt_inj is large enough that the deceleration and resultant fallback of matter once ejected in the first eruption event settles in.

Fig. 11. Peak luminosity (top panel) and total radiated energy (bottom panel) in the second eruption as functions of the interval Δtinj. The total radiated energy is evaluated by the cumulative luminosity from the time when the luminosity begins to rise to 200 days (20 days) after the peak for models RSG (BSG).

Each ejected fluid element that has a positive total energy (summation of thermal, gravitational, and kinetic energy) keeps expanding at an almost constant velocity until the core collapse of the star. Erupted material reaches ∼2 × 10¹⁵ cm (RSG) or ∼3 × 10¹⁶ cm (BSG) from the progenitor when the star undergoes core collapse. Their velocities are several 10⁶ cm s⁻¹ (RSG) or several 10⁷ cm s⁻¹ (BSG) (Fig. 12). These values do not contradict observations (e.g., Kiewe et al. 2012). When we express the velocity-density relation in Fig. 12 as ρ ∝ v^−s, the CSM formed by the first eruption (broken line) have larger s values compared with that by the second ejection (solid line), which roughly follow a line with s ∼ 1.5. This is because a part of photons diffusing out from the shock wave formed by the second ejection can be absorbed and accelerate the inner part of the CSM formed by the first eruption. This “steepening” of s value from ∼1.5 (and resulting CSM) can affect the light curve of SNe powered by CSM interaction modeled by Chevalier (1982).

Fig. 12. Density profiles of CSM as functions of velocity at the core collapse for the fiducial models. The broken lines (solid lines) correspond to the CSM erupted by the first injection (second injection). The gray line represents the relation ρ ∝ v−1.5, which we adopt as a fiducial slope in the case of eruptive mass loss (Kuriyama & Shigeyama 2020; Tsuna et al. 2020).

3.3. Compare light curves between the first and second eruptions

The light curves of first and second eruptions show different features (Figs. 3 and 4). In this section, we discuss the difference between the first and second eruptions using model RSG1-f as an example.

In the first eruption, the outward shock wave transports energy without leaking photons until it approaches the surface of the star, and thus we can observe strong shock breakout (Fig. 3). In the second eruption, photons diffuse out from the shock front long before it reaches the photosphere, where τ = 1, because the first eruption lowers the density from the original (MESA output) state. Thus photons travel longer distances to the photosphere, which leads to a broader peak in the light curve. We can see these diffusing out photons in Fig. 13. When t = 2.20 × 10⁷ s, the shock wave is located at M_r = 14.3 M_⊙, but the pre-shocked region (M_r >  14.3 M_⊙) has already heated compared to t = 1.85 × 10⁷ s owing to leaking photons from the shock. Correspondingly, the luminosity reaches peak at t = 2.20 × 10⁷ s (Fig. 3).

Fig. 13. Temperature and optical depth for model RSG1-f during the first eruption (left) and the second eruption (right). Broken lines (6.4E6 in the left panel and 1.85E7 in the right panel) correspond to profiles just before eruptions, long-broken lines correspond to profiles at the time of peaks in the luminosity due to eruptions, and solid lines show the profiles after the luminosity peaks.

The duration of the brightening phase is determined by the photon diffusion timescale. Therefore, firstly we consider the region whose radiation energy can contribute to the light curves and then estimate the diffusion timescale of the region. When the diffusion velocity v_diffusion in a certain region becomes greater than the wave velocity v_wave, photons can diffuse out from the wave. Since v_diffusion is written as

$$\begin{array}{c}\hfill v_{\mathrm{diffusion}}\simeq c/(3\tau ),\end{array}$$(1)

the radius R_τ, where v_diffusion ≃ v_wave holds at the time of shock breakout, has the optical depth

$$\begin{array}{c}\hfill \tau \simeq \frac{c}{3v_{\mathrm{wave}}}.\end{array}$$(2)

The radiation energy density at the place where τ <  c/(3v_wave) is transferred at the diffusion velocity and contributes the tail part of the light curve. Thus, we can obtain a rough estimate of the emergent radiation flux during the light curve peak. As a result, the luminosity L with the duration Δt is estimated as

$$\begin{array}{cc}\hfill L& =4\pi R_{\tau }^{2}a{T}_{\tau }^4\frac{c}{3\tau }=4\pi R_{\tau }^{2}a{T}_{\tau }^4{v}_{\mathrm{wave}},\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill \mathrm{\Delta }t& =\frac{R_{\mathrm{ph}}-R_{\tau }}{v_{\mathrm{diff}}}=\frac{R_{\mathrm{ph}}-R_{\tau }}{c/(3\tau )}.\hfill \end{array}$$(4)

We obtain v_wave ≃ 53.4 km s⁻¹, R_τ ≃ 4.91 × 10¹³ cm and R_ph ≃ 4.96 × 10¹³ cm in the first eruption from the results of model RSG1-f. The above equations give τ ≃ 1870, L ≃ 4.1 × 10³⁸ erg s⁻¹, and Δt ≃ 1 days. In the second eruption, v_wave ≃ 53.5 km s⁻¹, R_τ ≃ 8.43 × 10¹³ cm, and R_ph ≃ 1.26 × 10¹⁴ cm. Thus, τ ≃ 1870, L ≃ 3.2 × 10³⁸ erg s⁻¹, and Δt ≃ 89 days. A large difference in R_ph − R_τ between the two eruptions causes different behaviors of the two brightening events. After the first eruption, the expansion of the envelope significantly increases R_ph. Thus, photons need to travel longer distances from the shock front to the photosphere. This results in a broader peak of the light curve owing to the second eruption. The mass coordinate, where τ ≃ 1870, at shock breakout is M_r ∼ 14.6 M_⊙ for the first eruption and ∼14.3 M_⊙ for the second eruption (Fig. 13), and thus, the enclosed mass range that contributes to the light curve is also wider in the second eruption compared with the first eruption.

4. Discussion

Our results show unambiguous differences between the first and second mass eruption in terms of light curves, colors, and amounts of erupted unbound mass. Therefore, when we discuss the observation of pre-CCSN mass eruption events and compare them with some physical models or simulations, we may have to take into consideration the effect of repeated mass eruptions seen in our simulations and the altered density profile of the envelope. It should be noted that our results might be restricted to stars with an initial mass of 15 M_⊙. It is important to explore multiple eruption events from stars with different masses, although this is beyond the scope of the present paper.

The motivation of this work is to interpret the nature of multiple mass eruption events accompanied by rapid and repeated luminosity variance. In case of erratic pre-CCSN phase of SN 2009ip during 2011 March-November, its magnitude oscillated with an amplitude of about three mag (Pastorello et al. 2013). The period of this luminosity oscillation is less than 50 days and this phase lasts about ten months. The peak magnitude and the time interval between each local peak are almost unchanged during this phase. Although the progenitor mass of SN 2009ip (50 − 80M_⊙; Smith et al. 2010a) and our models (15 M_⊙) are quite different, these observed properties seem to be incompatible with our results. In our work, once the first mass eruption occurs, the density structure of the envelope is altered. Thus the local peak luminosity and color are significantly different between the first and second mass eruption (Figs. 3–6), although the same amount of energy is deposited on the same timescale (model RSG1-f, s, m, l and BSG1-f, s, m, l). In addition, our models predict a slower decline of the luminosity after the second local peak because of the extended envelope as shown in Fig. 14, where the peak luminosity of each eruption in R band is plotted as a function of half-decay time of the luminosity and some observations are plotted for comparison. A typical half-decay time of the luminosity in models BSG is ∼5 × 10³  s for the first eruption but ∼ several 10⁵  s for the second eruption, which is longer than that of SN 2009ip progenitor outbursts in 2009 and 2011. If these behaviors in our model hold for a more massive star including the progenitor of SN 2009ip (though it is not obvious), our results may suggest that the nature of mass eruption from SN 2009ip in 2011 cannot be explained by the spherically symmetric eruption of a single star triggered by extra energy supply into the envelope.

Fig. 14. Peak luminosity in R band as a function of the half-decay time of the luminosity for each eruption in our simulation. The bolometric luminosity in our simulation is filtered by Johnson-Cousins R band under the assumption of blackbody radiation. Filter transparency data are taken from the website of Leibniz Institute for Astrophysics Potsdam (AIP). 3064 Jy is adopted as the zero magnitude flux in R band (Bessell et al. 1998). Some observational data of SN progenitor outbursts (SN 2009ip; Pastorello et al. 2013; Fraser et al. 2013b, LSQ13zm Tartaglia et al. 2016, and SN 2015zh; Ofek et al. 2016) and gap transients (intermediate luminosity red transient (ILRT) AT 2017be; Cai et al. 2018, luminous red nova (LRN) SN Hunt248; Mauerhan et al. 2015; Kankare et al. 2015, and SN impostor SN 2002bu; Kato et al. 2004) in R band or r band are also plotted for comparison. Half-decay times of the luminosities for these transients are estimated by linear interpolation of observed magnitudes as a function of time. SN 2009ip shows a lot of peaks (events) in the light curve during 2009–2012. Four representative and well-observed pre-SN outbursts at MJD ∼ 55098, 55826, and the two events (2012a and 2012b) introduced in Sect. 1, are chosen. SN Hunt248, which is thought to be a merger-burst event (Mauerhan et al. 2018) and classified as LRN (e.g., Metzger & Pejcha 2017; Pastorello et al. 2019b), shows triple peaks (outburst events) referred to as 2014a, b, and c (Kankare et al. 2015). The 2014b event had begun before the luminosity of the 2014a event declined by a factor of 2 from the peak. Accordingly, only the 2014b and c events are plotted. It should be noted that objects in this figure are representative and that this figure does not show all of the known pre-SN outbursts or gap transients.

Mass eruptions treated in this paper show a wide variety of timescale and peak luminosity depending on the progenitor type (RSG or BSG) and “history of eruptions” (first or second eruption) (Fig. 14). They tend to be a bit darker than the known “gap transients” (Kasliwal 2012) but brighter than classical novae. If ∼15 M_⊙ stars experience these kinds of eruptions prior to core collapse, upcoming deep and high cadence transient surveys will fill in the new gap region between the known gap transients and classical novae.

Observations have revealed the diversity of SNe IIn (Taddia et al. 2013; Smith et al. 2017) in the amount, composition, position, and morphology of CSM. In particular, there is a wide range of variety in the amount of CSM. While some progenitors are associated with more than 30 − 50 M_⊙ (SN 2010jl; Zhang et al. 2012) or 20 M_⊙ (SN 2006gy; Smith et al. 2010b) of CSM, others only have 0.003 M_⊙ (SN 1998S; Fassia et al. 2001) or less than 0.1 M_⊙ (SN 2005gl; Gal-Yam & Leonard 2009). This variety could originate from different mechanisms of eruption to form CSM, of course, from a different progenitor mass. As introduced in Sect. 1, there are a wide variety of scenarios that may explain extra energy supply and resulting mass loss and CSM formation. On the other hand, our results show that only a factor of two difference in the amount of deposited energy results in a difference of ejected mass by more than two orders of magnitude (the second mass eruption in RSG1-f and RSG2-f, Table 3). Thus a wide range of the amount of CSM might also reflect a small difference in the amount of extra energy supplied.

When we compare the results of models RSG and BSG, we find that model BSG requires more than ten times larger amount of energy than model RSG to eject the same mass; both models RSG2-f and BSG1-f expel ∼0.7 M_⊙ throughout double eruption events. From this aspect, model RSG (i.e., stars in a relatively metal-rich environment) is preferred as a progenitor of SNe IIn if we only focus on the envelope density stratification and ignore the metallicity dependence of the physics involved in the extra energy injection. At a glance, our results seem to be consistent with Graham (2019), who suggests that SNe IIn prefer relatively higher metallicity environment compared with the other subtypes of SNe II based on a database analysis work. However, there should be various other physical factors correlated with metallicity and affecting the classification of a (sub-)type of an SN. Therefore, our result that metal-rich stars seem to easily suffer from mass eruption is just one of the possible factors SNe IIn prefer.

5. Conclusions

We have carried out radiation hydrodynamical simulations to investigate the properties of the repeated dynamical mass eruption events from a single 15 M_⊙ star prior to the core collapse, in the context of the formation of CSM around SNe IIn. The key point of this work is that we caused mass eruption events not only once but twice by depositing energy two times. The dynamical properties of single mass eruption event have already been studied by Kuriyama & Shigeyama (2020), Owocki et al. (2019), and Dessart et al. (2010). However, recent observations of SN IIn progenitors indicate that mass eruption phase often occurs more than once just before core collapse. A mass eruption event can alter the density structure of the envelope and therefore the subsequent mass eruption event can show completely different observational features. From this view point, we deposited the extra energy into the bottom of the progenitor envelope twice and simulated the time evolution of both eruption events. The progenitor models were made by MESA. The amount of deposited energy is chosen to reproduce observed CSM, and roughly set to one-third or one-quarter of |E_envelope| (binding energy of the envelope). We did not deal with the origin of extra energy, which has not been completely understood so far, although it should be studied in our future work.

Results show non-negligible differences between the first and second mass eruption events. After the first eruption event occurs, the envelope is temporarily inflated on the Kelvin-Helmholtz timescale. Thus, when the second mass eruption takes place consecutively in the inflated envelope, it shows a fainter, redder, and long-term eruption. This result may conflict with observations for SN 2009ip during the outburst phase in 2011, although the progenitor mass and properties of the envelope in SN 2009ip are quite different from those of our models. The symmetric mass eruption model from a single star could be ruled out for this event in SN 2009ip if our claims about the properties of second eruption hold for more massive stars. Erupted material from the second eruption collides with the prior CSM from the first eruption and alter the density profile of the prior CSM as described in Sect. 3.2. This interaction can affect the light curve of the subsequent CSM interacting SN. Another suggestion of our work is that only a few factors of difference in extra energy makes a larger (by a few orders of magnitude) difference in the amount of erupted mass. In other words, a little difference in the progenitor stellar evolution can be amplified and emerges as a significant difference in the observational features during the pre-CCSN phase and SN itself. From this view point, the difference of the progenitor star evolution between SNe IIn and normal SNe II could be smaller than we think. The peak luminosities of eruptions from 15 M_⊙ stars in this paper lie in between that of known gap transients and classical novae. Higher cadence and deeper transient survey in the future will provide us a larger amount of detailed pre-SN activity data and some objects will show multi-peaked luminosity fluctuations. Interpreting these would require the inclusion of the effect of multiple eruptions and the corresponding altered density structure discussed in this work.

Acknowledgments

We thank the anonymous referee for giving us a lot of valuable comments and Amy Mednick for language editing. This work is partially supported by JSPS KAKENHI NOs. 16H06341, 20H05639, MEXT, Japan.

References

Appendix A: Detailed methods of making the two progenitor models using MESA

The two progenitor models, RSG and BSG, which were used in our simulation as the initial models, were made via a stellar evolution code MESA release 10398 (Paxton et al. 2011, 2013, 2015, 2018, 2019).

According to Table 1, initial mass and metallicity were set to the following values and we started to make these models evolve.

initial_mass = 15.0d0
initial_z = 0.02d0 (for model RSG)
Zbase = 0.02d0 (for model RSG)
initial_z = 0.0002d0 (for model BSG)
Zbase = 0.0002d0 (for model RSG)

After the termination of main sequence, we set the following options which designate the method of opacity and mass loss.

use_Type2_opacities = .true.
cool_wind_RGB_scheme = 'Dutch '
cool_wind_AGB_scheme = 'Dutch '
RGB_to_AGB_wind_switch = 1d-4
Dutch_scaling_factor = 0.8

Parameters other than those above were set to default values throughout the entire calculation. We stopped the calculation 11.2 yr before core collapse (model RSG) and 7.2 yr before core collapse (model BSG) and adopted as initial models for our calculations.

All Tables

All Figures

Fig. 1. Evolution of the energy generation rate by nuclear burning for each model indicated by labels. The rate becomes higher toward the core collapse. There are some local peaks for each model around ten years, five years, and one year before core collapse. We adopted the peaks at 11.2 years before core collapse (model RSG) and at 7.2 years before core collapse (BSG) as initial models for our calculations. These peaks correspond to the core neon burning phase.

In the text

	Fig. 2. Density distribution of each progenitor model. Model RSG has a more extended envelope than model BSG because of higher metallicity (Table 1) and opacity. In addition to that, model RSG has experienced stronger mass loss due to higher metallicity and thus has a lower total mass.
In the text

	Fig. 3. Light curves for models RSG1-f, s, m, l, and RSG2-f. The first eruption triggered by the injected energy Einj1 produces a peak at day ∼80. The second peak in each model corresponds to the second eruption triggered by the injected energy Einj2.
In the text

	Fig. 4. Same as Fig. 3 but for models BSG1-f, s, m, l, and BSG2-f.
In the text

	Fig. 5. Evolution of effective temperature for each RSG progenitor model.
In the text

	Fig. 6. Evolution of effective temperature for each BSG progenitor model.
In the text

	Fig. 7. Radii as functions of the enclosed mass for the fiducial models. The three lines in each panel represent the profiles before the first energy injection (dotted line), between the first and the second energy injections (dashed line), and at the end of the hydrodynamical simulations (solid line).
In the text

	Fig. 8. Same as Fig. 7 but for the density profiles.
In the text

	Fig. 9. Velocity profiles as functions of the enclosed mass for the fiducial models. The three lines in each panel represent the profiles before the second energy injection (dotted line), after the second energy injection (dashed line), and at the end of the hydrodynamical simulations (solid line).
In the text

	Fig. 10. Trajectories of some fluid elements for model BSG2-f. The red bold line represents the trajectory of a matter with a mass coordinate Mr = 12.5 M⊙ assuming free-fall motion from the maximum radius.
In the text

	Fig. 11. Peak luminosity (top panel) and total radiated energy (bottom panel) in the second eruption as functions of the interval Δtinj. The total radiated energy is evaluated by the cumulative luminosity from the time when the luminosity begins to rise to 200 days (20 days) after the peak for models RSG (BSG).
In the text

	Fig. 12. Density profiles of CSM as functions of velocity at the core collapse for the fiducial models. The broken lines (solid lines) correspond to the CSM erupted by the first injection (second injection). The gray line represents the relation ρ ∝ v−1.5, which we adopt as a fiducial slope in the case of eruptive mass loss (Kuriyama & Shigeyama 2020; Tsuna et al. 2020).
In the text

	Fig. 13. Temperature and optical depth for model RSG1-f during the first eruption (left) and the second eruption (right). Broken lines (6.4E6 in the left panel and 1.85E7 in the right panel) correspond to profiles just before eruptions, long-broken lines correspond to profiles at the time of peaks in the luminosity due to eruptions, and solid lines show the profiles after the luminosity peaks.
In the text

Fig. 14. Peak luminosity in R band as a function of the half-decay time of the luminosity for each eruption in our simulation. The bolometric luminosity in our simulation is filtered by Johnson-Cousins R band under the assumption of blackbody radiation. Filter transparency data are taken from the website of Leibniz Institute for Astrophysics Potsdam (AIP). 3064 Jy is adopted as the zero magnitude flux in R band (Bessell et al. 1998). Some observational data of SN progenitor outbursts (SN 2009ip; Pastorello et al. 2013; Fraser et al. 2013b, LSQ13zm Tartaglia et al. 2016, and SN 2015zh; Ofek et al. 2016) and gap transients (intermediate luminosity red transient (ILRT) AT 2017be; Cai et al. 2018, luminous red nova (LRN) SN Hunt248; Mauerhan et al. 2015; Kankare et al. 2015, and SN impostor SN 2002bu; Kato et al. 2004) in R band or r band are also plotted for comparison. Half-decay times of the luminosities for these transients are estimated by linear interpolation of observed magnitudes as a function of time. SN 2009ip shows a lot of peaks (events) in the light curve during 2009–2012. Four representative and well-observed pre-SN outbursts at MJD ∼ 55098, 55826, and the two events (2012a and 2012b) introduced in Sect. 1, are chosen. SN Hunt248, which is thought to be a merger-burst event (Mauerhan et al. 2018) and classified as LRN (e.g., Metzger & Pejcha 2017; Pastorello et al. 2019b), shows triple peaks (outburst events) referred to as 2014a, b, and c (Kankare et al. 2015). The 2014b event had begun before the luminosity of the 2014a event declined by a factor of 2 from the peak. Accordingly, only the 2014b and c events are plotted. It should be noted that objects in this figure are representative and that this figure does not show all of the known pre-SN outbursts or gap transients.

In the text