Fast Blue Optical Transients Due to Circumstellar Interaction and the Mysterious Supernova SN 2018gep

Interaction between the supernova (SN) ejecta and dense circumstellar material (CSM) (Chevalier 1982) is one of the models (e.g., Woosley et al. 2007; Chatzopoulos et al. 2013; Moriya et al. 2013; Morozova et al. 2015; Sorokina et al. 2016; Blinnikov 2017; Kasen 2017; Smith 2017; Moriya et al. 2018) proposed to explain the diversity of highly luminous SNe, in parallel with magnetar-powered SNe (Maeda et al. 2007; Kasen & Bildsten 2010; Woosley 2010; Kasen et al. 2016), accretion-powered SNe (Dexter & Kasen 2013; Wang et al. 2018), and pair-instability SNe (Kasen et al. 2011; Gilmer et al. 2017). CSM interaction occurs when the rapidly expanding stellar ejecta collides with the quasi-static CSM. The ejecta drives a shock through the CSM, which converts the ejecta kinetic energy into thermal energy, leading to the bright event observed. The model has been applied to recent superluminous SNe such as SN 2008gy (Woosley et al. 2007), SN 2007bi, SN 2010gx and PTF 09cnd (Sorokina et al. 2016), iPTF 14hls (Woosley 2018), PTF 12dam (Tolstov et al. 2017), SN 2016iet (Lunnan et al. 2018), iPTF 16eh (Gomez et al. 2019), AT 2018cow (Leung et al. 2020b), and PS 15dpn (Wang & Li 2020).

To explain the origin of the CSM, a number of mechanisms can trigger greatly enhanced mass loss prior to the final stellar explosion, including common-envelope-triggered mass loss (Chevalier 2012; Schrøder et al. 2020), pulsation-induced mass loss in pulsational pair-instability SNe (PPISNe; Umeda & Nomoto 2003; Woosley 2017; Leung et al. 2019; Marchant et al. 2019; Leung et al. 2020a; Woosley 2019; Renzo et al. 2020a), enhanced stellar wind in super-asymptotic giant branch stars (Jones et al. 2013; Moriya et al. 2014; Nomoto & Leung 2017; Leung & Nomoto 2019; Tolstov et al. 2019), and wave-driven mass loss (Quataert & Shiode 2012; Shiode & Quataert 2014; Fuller 2017; Fuller & Ro 2018; Ouchi & Maeda 2019; Leung & Fuller 2020; Kuriyama & Shigeyama 2021).

Pulsation-induced mass loss relies on the electron–positron pair-creation catastrophe (Barkat et al. 1967), which happens in very massive stars (∼80–140 M_⊙; Heger & Woosley 2002; Umeda & Nomoto 2002; Ohkubo et al. 2009; Hirschi 2017; Yoshida et al. 2016). These radiation-dominated stars lose core pressure support from photons during their conversion into electron–positron pairs. The contraction of the star's core becomes dynamical, triggering explosive burning of C and O. The excess energy generation makes the core bounce, driving a shock through the envelope that ejects mass from the surface. The star can experience several mass-loss events, depending on the available carbon and oxygen in the core (Woosley 2017; Leung et al. 2019; Marchant et al. 2019; Woosley 2019; Renzo et al. 2020a).

In Leung et al. (2020b), we explored the possibility of applying the pulsation-induced mass-loss model to explain the rapid transient AT 2018cow. Like SN 2018gep, AT 2018cow is also classified as a "Fast Blue Optical Transient" (FBOT). That work, together with Tolstov et al. (2017), suggests that PPISNe provides the flexibility to span the wide diversity of transient objects from FBOT to superluminous SNe. However, the unusually rapid transient in SN 2018gep leads to speculation whether other mass-loss mechanisms are necessary to explain the optical signals of this object.

Wave-driven mass loss (Quataert & Shiode 2012; Shiode & Quataert 2014; Fuller 2017; Fuller & Ro 2018; Wu & Fuller 2021) relies on the vigorous convective motions in the massive star's core during its late-phase nuclear burning (e.g., carbon-burning and later advanced burning stages). This in turn excites internal gravity waves that propagate through the radiative core, where some of the wave energy is transmitted into the envelope via acoustic waves. When the waves reach the surface where the density gradient is large, they develop into weak shocks that dissipate and deposit their energy in the surroundings. Even though only a small fraction of the wave energy is leaked to the envelope, it can be sufficient to eject a substantial amount of mass.

SN 2018gep (ZTF 18abukavn) is an SN discovered by the Zwicky Transient Factory (ZTF; Bellm et al. 2019; Graham et al. 2019), first analyzed by Ho et al. (2019) and later by Pritchard et al. (2020). This object has photometric and spectroscopic features similar to some other recent rapid transients, such as AT 2018cow (Prentice et al. 2018; Smartt et al. 2018; Margutti et al. 2019; Perley et al. 2019) and iPTF 16asu (Wang et al. 2019; Whitesides et al. 2017). SN 2018gep is remarkable for its very early detection by ZTF and its proximity at ∼143 Mpc from Earth, which allowed for detailed follow-up observations.

Among all rapid transients, SN 2018gep has a fast rise time of 0.5–3 days (Ho et al. 2019) until it reached a large peak luminosity of 3 × 10⁴⁴ erg s⁻¹. Its low-mass and compact host galaxy has a metallicity of only one-fifth of solar metallicity. The surface temperature of the transient at its peak of ∼40,000 K is among the highest observed in stripped-envelope SNe. The upper limits of detection by high-energy bands (X-ray and gamma-ray) and the radio band have led to speculations that this explosion arises from interaction with a compact and dense CSM. Future polarization and X-ray observations of this object can further pin down the possibility of dense CSM (see, e.g., Huang et al. 2019; Rivera Sandoval et al. 2018, for applications in AT 2018cow.)

In Section 2, we perform hydrodynamical modeling of mass loss from helium star SN progenitors due to wave heating, predicting the resulting CSM structure at the time of explosion. In Section 3, we compare the CSM properties from our wave-driven mass-loss models to radiative transfer models of this work and those in Ho et al. (2019). We also include PPISN models from the literature to contrast the corresponding CSM properties of this class of stars. In Section 4, we present radiative transfer modeling and an explosion parameter search to fit the bolometric light curve of SN 2018gep. We then discuss the robustness of our results in Section 5, examining possible degeneracies, comparing to other models of SN 2018gep, and discussing the physical origins of this event. Finally we summarize the findings of this work and highlight our interpretation of SN 2018gep.

To investigate wave-driven mass loss and to construct precollapse models for radiative transfer modeling, we use the one-dimensional stellar evolution code MESA (Module for Experiments in Stellar Astrophysics; Paxton et al. 2011, 2013, 2015, 2018), version 8118. We model stars with initial masses M_ini = 20–80 M_⊙, which correspond to He cores with masses M_He ∼ 5–40 M_⊙ (Hirschi 2017; Woosley 2017). In addition to the default massive star settings provided in MESA, we set the star to be nonrotating, with a Dutch mass-loss ³ scaling factor η = 0.8 and an overshoot parameter f_ov = 0.001. The final mass before collapse, even without dynamical mass loss, depends strongly on the initial metallicity due to mass loss from line-driven winds. We follow the entire stellar evolution from the onset of core H burning until the onset of gravitational collapse. ⁴

We consider only stripped-envelope stars in this work because this class of stars produces hydrogen-free Type Ib/c SNe like SN 2018gep. We first model a zero-age main-sequence (ZAMS) star until hydrogen is exhausted in the core. Then we remove the H envelope by relaxing the mass to exclude the H-rich matter and continue the evolution. The loss of the majority of the H envelope is most likely to arise from a binary interaction (either stable or unstable mass transfer). While binary stripping may not remove all of the hydrogen in reality, the remaining hydrogen is quickly lost via winds at the masses and metallicities we consider. The final mass is thus primarily controlled by (uncertain) wind mass loss rather than the details of the binary interaction.

In Table 1, we list relevant masses for the progenitor models in this study. Although the ZAMS mass is the model parameter we control, we will focus on the pre-explosion mass ${M}_{\exp }$, as it determines the SN ejecta mass, along with the core evolution and hence the wave-driven mass loss. The core mass is defined by the outer mass coordinate where the abundances of those elements drop below threshold values, where we use the default value 0.01.

Table 1. The Initial, He-core, C-core, and Pre-explosion Masses of the Models Studied in This Work

Mini (M⊙)	MHe(M⊙)	MC (M⊙)	${M}_{\exp }$ (M⊙)
25	6.56	3.00	5.44
40	12.97	5.58	9.20
50	17.45	8.39	10.46
60	21.78	9.58	12.22

Note. We Assume Z = 0.02 And A Dutch wind1 coefficient n = 0.8. The He-core mass is taken at the end of core H-burning. The C-core mass is taken from the values at the end of simulations (The definition of "the core mass" is given in the text of Section 2.1).

Download table as:  ASCIITypeset image

To investigate the process of wave-driven mass loss, we follow a scheme similar to that described in Fuller & Ro (2018). Once an oxygen-burning core has developed in the star, we terminate the simulation. We set the inner boundary at 1 M_⊙, interior to the carbon core outer boundary. We keep the outer carbon layer as a buffer, but in general, the mass loss is very small and the carbon layer remains unchanged.

To calculate how the wave energy is deposited in the envelope, we follow the method of Fuller & Ro (2018), in which the outgoing acoustic wave luminosity L_wave decreases as

$$\begin{eqnarray}&&\displaystyle \frac{{{dL}}_{\mathrm{wave}}}{{dM}}=-\displaystyle \frac{{L}_{\mathrm{wave}}}{{M}_{\mathrm{wave},\mathrm{shock}}}-\displaystyle \frac{{L}_{\mathrm{wave}}}{{M}_{\mathrm{wave},\mathrm{diff}}},\end{eqnarray} \tag{ 1 }$$

The local value of dL_wave/dM is the deposited energy per unit mass within the star. Above, M_wave,shock and M_wave,diff are the expected damping mass by shock heating and radiative diffusion, given by

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

$$\begin{eqnarray}&&{M}_{\mathrm{wave},\mathrm{diff}}=\displaystyle \frac{2{L}_{\max }}{{\omega }^{2}K},\end{eqnarray} \tag{ 3 }$$

where γ, c_s ², and K are the adiabatic index, sound speed squared, and thermal diffusivity (Fuller 2017). ${L}_{\max }=2\pi \rho {r}^{2}{c}_{s}^{3}$ is the maximum energy transportable by linear acoustic waves, and ω is the wave frequency. The wave frequencies excited by convection are typically in the range ω = 10⁻²–10⁻³ s⁻¹ and are a model parameter in this work. Additionally, the outflowing matter can expand with supersonic velocity, stretching the acoustic wavelength and increasing the value of M_wave,diff and M_wave,shock as described in Fuller & Ro (2018).

Similar to our previous work (Leung & Fuller 2020), when we model the envelope dynamics, we treat the duration of energy deposition as a parameter. This allows us to model the resulting interaction-powered light curve as a function of both the wave-energy deposition and its duration, but we note that the two parameters are linked when modeling the inner core evolution of the star, as investigated in Shiode & Quataert (2014) and Wu & Fuller (2021). We retain all of the ejected mass shells, which later form the CSM, on our model grids because the typical timescale from the onset of O burning to the final collapse is ∼0.1 yr, and the corresponding CSM radius is less than 10⁴ R_⊙.

Below, we find that the mass-loss rates via wave-energy deposition are much larger than those due to line-driven winds, which are typically ${\dot{M}}_{\mathrm{wind}}\sim {10}^{-4}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ near the end of the progenitor's life. Hence, while winds remove much more mass than wave heating during prior phases of evolution, wave-driven mass loss dominates over wind mass loss in the final months of the star's life, so we ignore wind mass loss when constructing CSM density profiles around the progenitor star.

We first study how the CSM properties depend on the energy deposition duration. We consider a 5.44 M_⊙ He star with Z = 0.02. In the energy deposition phase, we deposit a total of 6 × 10⁴⁷ erg energy in the envelope, a value near the upper end of the envelope wave-energy deposition range found in Wu & Fuller (2021). The duration, which depends on the exact duration of O burning, is treated as a free parameter from 0.05 to 1.35 yr, which covers the typical O-burning duration for a wide range of stellar mass. In Figure 1, we plot the hydrodynamical profiles of these models at the end of the simulations.

Zoom In Zoom Out Reset image size

The density profile shows a generally smooth structure with radial dependence similar to (though slightly steeper than) that for a steady wind, ρ ∝ r⁻². The density profile is also punctuated by spikes associated with internal shocks, which arise due to the uneven energy deposition (and hence uneven outflow velocity) at early times. Near its outer edge, the CSM has a typical density of 10⁻¹⁰–10⁻¹³ g cm⁻³. The temperature also falls rapidly with radius, except in the outermost optically thin layers where it is constant. The ejecta moves away from the core with a typical velocity of ∼3 − 5 × 10⁷ cm s⁻¹, but with an exception for the longest energy deposition model (t_dep = 1.35 yr), in which the expansion of the star is quasi-hydrostatic. The escape velocity at large radii is about 3 × 10⁶ cm s⁻¹, confirming that the ejected matter can successfully escape from the gravity of the star.

The luminosity profile demonstrates how the star distributes the deposited wave energy. In the bound portion of the star, the luminosity is constant at ≈10⁵ L_⊙. Across the energy deposition zone, the stellar luminosity quickly rises by one to two orders of magnitude but is smaller than the total wave-heating rate. This shows that much of the wave energy is used to unbind material (i.e., increase its gravitational potential energy) and convert it to kinetic energy of the outflow. An exception is the model with t_dep = 1.35 yr, in which matter is not efficiently ejected and most of the wave energy is radiated from the star. A long t_dep corresponds to a low L_dep = E_dep/t_dep for a given deposited energy E_dep. The heated layer has sufficient time to quasi-hydrostatically expand. This prevents the formation of fast winds and/or shocks and hence reduces mass ejection. Note also that the star's luminosity can decrease or increase within the outflowing material, the latter occurring as kinetic energy from colliding shells is converted into heat that is radiated outward.

The external mass coordinate profile of M−m shows the radial dependence of the ejected mass, with most of the mass located in the outer part of the CSM, where the M−m profiles turn over between roughly 300–3000 R_⊙. We see that the ejected mass is a few hundredths of a solar mass for all but the longest energy deposition timescale.

We next examine the energy deposition profile when the simulation ends, recalling that the energy deposition zone is time dependent. The energy deposition zone moves toward the stellar core as the heated zone expands. We also see that the waves damp at a smaller radius when t_dep becomes small, because the corresponding L_wave = E_dep/t_dep becomes large, so that the acoustic waves develop into weak shocks at a higher density, and thus at a smaller radius. The peak of the heat deposition also occurs near the sonic point of the outflow, i.e., near the boundary between the nearly hydrostatic star and the outflowing CSM.

We next compare how the total energy deposition changes the CSM profile of the ${M}_{\exp }$ = 5.44 M_⊙ model, fixing the deposition duration at 0.15 yr. In general, the net energy deposition is determined by the core structure and evolution, and Wu & Fuller (2021) find typical energies of a few × 10⁴⁷ erg for a wide number of models, with many massive star models in the range 10⁴⁷–10⁴⁸ erg. Hence, we consider total energy depositions of 2 × 10⁴⁷ erg, 6 × 10⁴⁷ erg, and 1.8 × 10⁴⁸.

Figure 2 shows that the total energy deposition is another primary factor affecting the CSM. While the density and temperature profiles are fairly similar between the models, the ejecta mass varies by about an order of magnitude, with the model with E_dep = 1.8 × 10⁴⁸ erg ejecting about 0.063 M_⊙, whereas the model with 2 × 10⁴⁷ erg ejects about 0.009 M_⊙. The ejected mass scales approximately linearly with the injected energy. The CSM also extends to slightly larger radii for larger energy deposition. The ejecta velocity shows only minor variations between the models, with the low-energy model exhibiting slightly larger velocity at the end of the simulation. The luminosity also does not monotonically change with E_dep, with a much higher luminosity for the largest E_dep model.

Zoom In Zoom Out Reset image size

We proceed to examine how the wave-heating process affects helium stars of different mass. A higher progenitor mass leads to a more compact core with a smaller radius during oxygen burning, which results in an envelope of higher binding energy. In Figure 3 we compare models with a pre-explosion mass of ${M}_{\exp }$ = 5.44, 9.20, and 12.22 M_⊙ (corresponding to ZAMS models of M_ini = 25, 40, and 60 M_⊙, respectively; see Table 1). For ease of interpretation, we assume a fixed energy deposition of 6 × 10⁴⁷ erg for a duration of 0.15 yr, similar to the O-burning duration of an M_ini = 40 M_⊙ star.

Zoom In Zoom Out Reset image size

In Figure 3, we plot the hydrodynamical profiles of these models at the end of these simulations. Higher stellar mass leads to a much lower amount of mass being ejected, from more than 10⁻² M_⊙ in the ${M}_{\exp }$ = 5.44 M_⊙ model, down to about 10⁻⁴ M_⊙ in the ${M}_{\exp }$ = 12.22 M_⊙ model. The outermost layers of the CSM extend to a few thousand solar radii in each case. Again the density profile falls off slightly steeper than r⁻², and the temperature profile flattens at a temperature of T ≈ 6000 K in the optically thin outer regions of the ejecta. The luminosity of the models is about 10 to a hundred times the original luminosity and increases with ${M}_{\exp }$.

In this section, we use radiative transfer modeling to derive constraints for the CSM based on the observational data of SN 2018gep. We use the SN explosion code (SNEC; Morozova et al. 2015) to compute the post-explosion light curve. This code follows the prototype presented in Bersten et al. (2011, 2013). The code has been widely used in the literature so we refer the reader to the instrumentation paper about the detailed implementation and its test cases. In this section, we choose ${M}_{\mathrm{ini}}=25\,{M}_{\odot }({M}_{\exp }=5.44\,{M}_{\odot })$ as our default pre-explosion model. We also included a model of M_ini = 50 M_⊙, whose He core of M_He = 17.45 M_⊙ is reduced to a pre-explosion mass of ${M}_{\exp }$ = 10.46 M_⊙ by winds (see Table 1). The model is used for comparison in Figure 10. The ejecta mass is chosen by selecting the mass cut. Here we only briefly review the our input.

After constructing ordinary stellar evolutionary models (i.e., models without any wave heating), we extract the density and chemical composition profile from MESA. Then we determine some input parameters necessary for SNEC. This includes (1) the inner mass cut M_cut, (2) the explosion energy ${E}_{\exp }$, and (3) the ⁵⁶Ni distribution inside the star X(⁵⁶ Ni). The minimum mass cut is the mass coordinate of the Fe core, assuming the Fe core to be promptly collapsing due to gravitational instability. A larger mass cut corresponds to an explosion with fallback or an aspherical explosion.

For ease of model comparison, we do not use the exact progenitor models from Section 2. Instead, we add parameterized CSM to our models, with properties motivated by wave-driven and pulsation-driven mass loss. Outside the original stellar surface, we add a CSM that extends to radius R_CSM with mass M_CSM. The CSM is initially isothermal with a density dependence ρ ∝ 1/r². For PPISN models, the extended and smooth mass ejection history reported in Leung et al. (2019, Figures (23)–(26)) also suggests a similar CSM profile. These parameters should not be over-interpreted since the real CSM density profile is not well predicted by stellar evolution models. Similar to Leung et al. (2020b), we assume the ejecta is uniformly mixed. In Figure 5, we plot the chemical abundance of the default model for illustration.

Zoom In Zoom Out Reset image size

We emphasize that the light-curve models are most sensitive to the ejecta–CSM interaction and are not very sensitive to the progenitor model. Therefore, for a given ejecta mass and composition, the results and constraints on CSM are applicable to both wave-driven mass loss and PPISN models. In Section 5.5, we further discuss how the aspherical explosion of massive stars is connected to these properties.

First we present our default model, which fits many of the features in the bolometric light curve of SN 2018gep. To find a good match, we need to search over the progenitor mass M_ini, ejecta mass M_ej, CSM mass M_CSM, CSM radius R_CSM, explosion energy ${E}_{\exp }$, and ⁵⁶Ni mass in the ejecta M(⁵⁶ Ni). In such a high-dimensional parameter space, we do not attempt to fit the light curve perfectly. Instead, we try to understand how the light-curve shape depends on each of the parameters. Table 2 lists the parameters of our default model.

Table 2.  Parameters of Our Default Progenitor and Explosion Model for SN 2018gep

Parameter	Value
Pre-explosion mass (${M}_{\exp }$)	5.44 M⊙
MHe, MC	6.56, 3.00 M⊙
Explosion energy (${E}_{\exp }$)	1 × 1052 erg
Ejecta mass (Mej)	2.0 M⊙
CSM mass (MCSM)	0.3 M⊙
CSM radius (RCSM)	1100 R⊙
Nickel mass (M(56 Ni))	0.33 M⊙
CSM slope	−2

Download table as:  ASCIITypeset image

In Figure 6, we plot the bolometric light curve (top left), photospheric radius (top right), effective temperature (bottom left), and photospheric velocity (bottom right) of our default model. For the bolometric light curves, we also include control experiments to demonstrate how the ⁵⁶Ni decay (green dotted line) and the pure CSM interaction (purple dashed line) contribute individually to the formation of the full light curve (blue solid line). SN 2018gep shows a two-component structure in the light curve. In the first several days, there is a plateau in the light curve where the luminosity is ∼10⁴⁴ erg s⁻¹. During this time, the luminosity is almost entirely created by thermal diffusion out of the shock-heated CSM. Around day 5, the CSM contribution decreases rapidly as the CSM becomes optically thin. The light curve falls more steeply until it flattens again when ⁵⁶Ni decay dominates the luminosity.

Zoom In Zoom Out Reset image size

In this work, we define t = 0 as the moment of explosion of the model, which takes place at UT 2018 September 8 20:20. This is about 0.3 days before the definition of t = 0 in Ho et al. (2019) based on the g-band light-curve brightening. Indeed, shock breakout occurs roughly 0.3 days after explosion in our default model, roughly consistent with the observed time of rapid brightening.

The observed photospheric radius R_ph of SN 2018gep shows a steady rise from an initially high value of ∼3 × 10¹⁴ to ∼3 × 10¹⁵ cm. Our models show similar evolution, albeit with a more rapid initial increase in R_ph and a less rapid late-time increase in R_ph. After day 5, there is a small but sudden drop in R_ph in our models when the shock-heated CSM becomes transparent. For the effective temperature T_eff, the default model follows the trend of SN 2018gep fairly well. The model drops from ∼10⁵ K on day 1 to ∼10⁴ K at day 6, falling more gradually at later times. The largest discrepancy appears to be at the points between days 3–6, where our model underpredicts the temperature and overpredicts the photospheric radius. The imperfect fit could also indicate a CSM with a slightly larger radial extent or different density structure could better match the data.

We note the small bumps in the light curve around day 1.5 and day 4.5. These bumps correspond to short temperature plateaus at the photosphere that occur at the first and second recombination of helium. The nearly constant photospheric temperature but increasing photospheric radius creates small and brief increases in the bolometric luminosity.

Lastly, we examine the evolution of the photospheric velocity, which is defined in the models as the Lagrangian velocity at the photosphere. The observed velocities of Ho et al. (2019) are measured using the C/O lines at early times (which are difficult to measure) and the Fe ii line at late times, which typically forms at smaller optical depth (and hence larger velocity) than the blackbody photosphere (Morozova et al. 2020). At most times, the model velocities are within ∼20% of those observed. The model photosphere expands at a very high velocity of ∼3−4 × 10⁹ cm s⁻¹ at early times when it is located within the CSM, but it falls sharply to about 2 × 10⁹ cm s⁻¹ when the CSM becomes optically thin. This sharp drop is not apparent in the data, though the sparseness and large uncertainties of the data prevent a detailed comparison.

We now study how the CSM shock cooling takes place in the default model by examining the evolving hydrodynamical structure. Figure 7 shows the density (top left), temperature (top right), velocity (middle left), the opacity (middle right), the free electron fraction (bottom left), and optical depth (bottom right) at several times after the explosion. Note that the innermost ejecta layers have a mass coordinate of 3.4 M_⊙ due to the mass cut used in this model. The initial density profile is that of the stellar model and the constant-density CSM. When the shock hits the CSM, it creates a density bump at the contact discontinuity, which gradually vanishes as the star expands. Afterward, the density decreases steadily due to nearly homologous expansion, apart from another density bump that develops near the middle of the ejecta after day 3 at M(r) ∼ 4.3 M_⊙. This is produced by a reverse shock propagating inwards from the contact discontinuity with the CSM. The compression associated with the reverse shock creates a density bump. The reverse shock evidently dies out at M(r) ∼ 4.3 M_⊙, but the density bump is sustained long afterward.

Zoom In Zoom Out Reset image size

As expected, the velocity increases with mass coordinate, and the star reaches locally homologous expansion by day 1. Before that, at days 0.01 and 0.1, we capture the moment when the shock hits the CSM and propagates through it. The majority of matter in the progenitor star (up to 5.44 M_⊙) expands with a velocity ∼2 × 10⁹ cm s⁻¹, while the low-density CSM obtains much larger velocities of 3−4 × 10⁹ cm s⁻¹. The extremely thin layer of outer CSM can reach as high as 10¹⁰ cm s⁻¹.

Like the density, temperature generally decreases with time due to (initially adiabatic) expansion. The early shock–CSM collision leads to a hot outer layer of T > 10⁵ K at early times. A temperature bump propagates inward due to the reverse shock propagating through the inner ejecta, which weakens as the star expands. At late times, the ejecta approaches isothermality as it becomes optically thin.

The opacity profile is closely linked to the temperature profile. After shock heating, the opacity rises by a factor of 2–3 where the shock-heated CSM is the most opaque. The opacity is typically largest where the temperature is ∼10⁵ K and free–free absorption dominates. At much larger temperatures, electron scattering dominates the opacity and κ ≈ 0.2 cm² g⁻ for the hydrogen-free ejecta. At much smaller temperatures (T ≲ 10⁴ K), the ejecta recombines and the opacity plummets. This occurs starting around day 5, after which the free-electron fraction and opacity decrease sharply.

Finally we examine the optical depth evolution in the bottom right panel of Figure 7, where the horizontal line at τ = 2/3 defines the photosphere. In the first 3 days, the photosphere lies in the outermost layers because the shocked-heated CSM is opaque. After day 3, the shocked CSM cools and begins to recombine. This causes the optical depth to decrease rapidly such that the photosphere recedes, consistent with the sudden drop of photospheric radius on day 7 in Figure 6. By day 10, the photosphere has receded into the stellar ejecta below the CSM, and by day 15 nearly all of the ejecta is optically thin.

In this section, we examine how parameter variations away from the default model affect the shape of the light curve.

4.4.1. Dependence on Explosion Energy

The first parameter we examine is the explosion energy. Because SN 2018gep was a broad-lined Ic SN, the explosion was likely powered by energy sources such as a rapidly rotating magnetar or accretion onto a black hole, which can power an explosion with an energy of order 10⁵² erg. In Figure 8, we plot the SN light curve for models with different explosion energies. The explosion energy has greater importance in the shock-cooling phase (before day 6), where we expect the emergent luminosity to scale linearly with the explosion energy (see Section 5.2), ${L}_{p}\propto {E}_{\exp }$. It is less important in the radioactive decay phase (after day 6) where the luminosity is determined primarily by the ⁵⁶Ni mass. The duration of the shock cooling phase also matches expectations, scaling approximately as ${t}_{p}\propto {E}_{\exp }^{-1/4}$ (Section 5.2), due to the faster expansion of the ejecta that shortens the photon diffusion time.

Zoom In Zoom Out Reset image size

4.4.2. Dependence on CSM Mass

4.4.3. Dependence on CSM Radius

4.4.4. Dependence on ⁵⁶Ni Mass

We next examine how the ⁵⁶Ni mass affects the light-curve structure. In Figure 9, we compare the default light curve with models containing no ⁵⁶Ni or two times more. The zero ⁵⁶Ni model is identical to that in Figure 6. Again, the early time light curves are essentially identical in all three cases, demonstrating that the ⁵⁶Ni and ⁵⁶Co decay do not contribute appreciably to the early time evolution of the light curve. On the other hand, the amount of ⁵⁶Ni determines the luminosity of the late-time light curve, similar to normal Type Ib/c SNe. Because the CSM mass and ⁵⁶Ni mass are both ∼0.3 M_⊙, it is unlikely that ⁵⁶Ni comes from pulsation and is more likely to be generated during core collapse.

Zoom In Zoom Out Reset image size

4.4.5. Dependence on Ejecta Mass

4.4.6. Dependence on CSM Density Profile

Only a couple of detailed models have been published for SN 2018gep. Based on the sharp rise and high surface temperature, Ho et al. (2019) argued that an extended shock breakout in a circumstellar medium (CSM) powers the initial peak. Their gray-opacity radiative transfer calculation based on analytic CSM profiles found the light curve could be approximately reproduced with M_ej = 8 M_⊙, explosion energy E_ej = 2 × 10⁵² erg, along with a CSM in a shell with mass M_CSM = 0.02 M_⊙ and radius R_CSM = 4000 R_⊙ with a small width 400 R_⊙. These numbers indicate the possibility that this compact SN Ic-BL has experienced mass loss very shortly before the explosion. The model of Ho et al. (2019) is distinctive from our default model above, as they predict a higher R_CSM but a lower M_CSM. While their explosion energy is also larger, the energy per unit mass is similar, as they find ${E}_{\exp }/M=2.5\times {10}^{51}$ erg M_⊙ ⁻¹, and we find ${E}_{\exp }/M\,=3.3\times {10}^{51}$ erg/M_⊙ ⁻¹.

To attempt to distinguish between the competing models, we repeat the light-curve calculation with SNEC, using their values as the model input. We thus choose a higher progenitor mass model (M_ini = 50 M_⊙, with a pre-explosion mass ${M}_{\exp }\sim 10\,{M}_{\odot }$) such that we can allocate 8 M_⊙ of ejecta mass. In Figure 10, we plot the light curve using their derived values as numerical inputs, along with our default model but with a M_CSM similar to theirs. After sharp early light-curve peaks, all of these models fade much more rapidly than the data because the low-mass CSM becomes optically thin very quickly. The luminosity plateaus after day 10, when ⁵⁶Ni decay dominates the energy source. We also repeat the experiment with a larger CSM mass M_CSM = 0.04 M_⊙, but the light curve still fades too fast.

Zoom In Zoom Out Reset image size

We suspect the main difference (see discussion below) is that the constant opacity used in Ho et al. (2019) is larger than the temperature-dependent opacity in our models. The use of a constant opacity allows the CSM to remain optically thick even after recombination, which slows down the luminosity drop in an unphysical manner. Another difference between their modeling and ours is the use of a thin shell of CSM in Ho et al. (2019), which creates a higher CSM density for a given CSM mass, thus allowing the high luminosity to be sustained for a longer time. These model differences also change the interpretation of the data between days 1–4. In our models, shock cooling of the envelope creates a plateau in the first few days, while the extended shock breakout models of Ho et al. (2019) create a large spike in between the data points at 1 and 3 days. Available photometry during that time seems to indicate a luminosity plateau rather than a spike, favoring our shock-cooling model. However, the lack of UV data prevents an accurate bolometric luminosity measurement, so we cannot confidently exclude the extended shock breakout model.

To get a rough estimate of the explosion properties of SN 2018gep, Pritchard et al. (2020) compared multicolor light curves to semianalytic MOSFiT models (Guillochon et al. 2018). In their best-fit CSM interaction model, they derive an ejecta mass 0.49 M_⊙, ⁵⁶Ni mass of 0.13 M_⊙, CSM mass of 0.11 M_⊙, and CSM radius of 3000 R_⊙. Their inferred CSM mass is in between our shock-cooling model (0.3 M_⊙) and the shock breakout model (0.02 M_⊙) from Ho et al. (2019), and their CSM radius is closer to that of Ho et al. (2019). However, the best-fit constant opacity of κ = 0.58 cm²g⁻¹ for their CSM models is likely unphysically large based on the opacities from our more realistic models (Figure 7). This problem is even more pronounced for their best-fit magnetar models, so magnetar powering of SN 2018gep is unlikely. For these reasons, we believe the parameters of our more detailed CSM interaction models are closer to the actual characteristics of SN 2018gep. However, a somewhat larger CSM radius and smaller CSM mass than our default model is certainly possible given degeneracies in light-curve fitting (see Section 5.3).

For an early light curve is dominated by shock cooling of the extended CSM, our models reveal clear trends in terms of light-curve plateau luminosity L_p , duration t_p , and slope. Here we compare with the scalings predicted by analytic models (Piro et al. 2021). The asymptotic ejecta speed (Matzner & McKee 1999) is

$$\begin{eqnarray}&&{v}_{t}\sim {\left(\displaystyle \frac{{E}_{\exp }}{m(r)}\right)}^{1/2}{\left(\displaystyle \frac{m(r)}{\rho {r}^{3}}\right)}^{\beta },\end{eqnarray} \tag{ 4 }$$

with m(r) the enclosed mass, and β ≃ 0.19. Our models are characterized by a radially extended but low-mass CSM with M_CSM ≪ M_ej. Hence, within the CSM, m(r) ≃ M_ej and ρr³ ∼ M_CSM, such that the CSM ejecta velocity is roughly

$$\begin{eqnarray}&&{v}_{t}\sim {E}_{\exp }^{1/2}{M}_{\mathrm{ej}}^{\beta -1/2}{M}_{\mathrm{CSM}}^{-\beta }.\end{eqnarray} \tag{ 5 }$$

The duration of the shock-cooling light curve is determined by the diffusion timescale t_d of photons out of the shocked CSM,

$$\begin{eqnarray}&&{t}_{d}\sim {\left[\displaystyle \frac{\kappa {M}_{\mathrm{CSM}}}{{v}_{t}c}\right]}^{1/2},\end{eqnarray} \tag{ 6 }$$

where κ is the opacity. At early times (t ≲ t_d ), the luminosity evolves as

$$\begin{eqnarray}&&L(t)\sim \displaystyle \frac{{{cR}}_{\mathrm{CSM}}{v}_{t}^{2}}{\kappa }{\left(\displaystyle \frac{{t}_{d}}{t}\right)}^{4/(n-2)},\end{eqnarray} \tag{ 7 }$$

The luminosity and photospheric radius fall sharply after a timescale

$$\begin{eqnarray}&&{t}_{\mathrm{ph}}\sim {\left[\displaystyle \frac{\kappa {M}_{\mathrm{CSM}}}{{v}_{t}^{2}}\right]}^{1/2}.\end{eqnarray} \tag{ 8 }$$

In the constant opacity approximation, the plateau luminosity is therefore expected to scale approximately as ${L}_{p}\,\propto {R}_{\mathrm{CSM}}{E}_{\exp }{M}_{\mathrm{CSM}}^{-2\beta }{M}_{\mathrm{ej}}^{2\beta -1}$. The light curve will drop sharply after a plateau duration ${t}_{\mathrm{ph}}\propto {E}_{\exp }^{-1/2}{M}_{\mathrm{CSM}}^{1/2+\beta }{M}_{\mathrm{ej}}^{1/2-\beta }$. Hence, for a constant plateau luminosity and duration, we expect ${R}_{\mathrm{CSM}}\,\propto {M}_{\mathrm{CSM}}^{-1}$ and ${E}_{\exp }\propto {M}_{\mathrm{CSM}}^{1+2\beta }{M}_{\mathrm{ej}}^{1-2\beta }$. The scalings provided above are approximately consistent with the light-curve modeling results of the previous sections.

The above analytic relations for the light-curve luminosity and plateau duration predict degeneracies between the CSM mass, radius, and explosion energy, as are well known for Type II-P SNe (Goldberg & Bildsten 2020). To understand whether our inferred CSM parameters for SN 2018gep are subject to such degeneracies, we computed additional numerical models, which are analytically predicted to have the same plateau luminosity and duration. In Figure 11, we plot the light curves for the default model (solid blue line) and a contrasting model with 2.60 × E_exp, 2 × M_CSM, and 0.5 × R_CSM (Model 1, purple dashed line), and another with 0.38 × E_exp, 0.5 × M_CSM, and 2 × R_CSM (Model 2, green dashed line).

Zoom In Zoom Out Reset image size

We find that the analytic scaling relations work fairly well (though not perfectly) in the parameter range considered, as the numerically computed luminosity, photospheric radius, and effective temperature for the three models are similar. The more compact model (Model 1) has a slightly lower peak luminosity and slower fall in luminosity during the shock-cooling phase. Minor deviations from the analytic predictions occur because the differing density and temperature change the opacity, which affects how the photosphere recedes. These results echo the findings of Dessart & Hillier (2019); Goldberg & Bildsten (2020) for Type II-P SN light curves, where they showed that the light-curve shape allows for degeneracy between the mass and radius of the progenitor's hydrogen envelope.

However, in contrast to Type II-P SNe where the envelope mass dominates the ejecta mass, the inferred CSM mass for SN 2018gep is much smaller than the total ejecta mass. This means the photospheric velocity scales approximately as v ∝ E_exp/M_ej and is not plagued by the same degeneracy noted by Goldberg & Bildsten (2020) for Type II-P SNe, so the photospheric velocity can potentially be used to break the degeneracy. In this case, a higher E_exp (Model 1) leads to a higher photosphere velocity, which predicts photospheric velocities substantially larger than those observed at early times (Figure 11). A better match is obtained for our default model or lower E_exp (Model 2). We caution that the measured Fe II line velocity is highly uncertain and is likely to be slightly larger than the continuum photospheric velocity because it is formed farther out in the ejecta (e.g., Paxton et al. 2018). For example, in iPTF 14hls (Arcavi et al. 2017), the photosphere radius derived from the Fe II line exceeds that of the blackbody radius after 100 days, exemplifying that the two radii may be very different.

We conclude that the CSM properties likely lie somewhere in between the default model and Model 2, i.e., in the range 0.15 M_⊙ ≲ M_CSM ≲ 0.3 M_⊙, 1100 R_⊙ ≲ R_CSM ≲ 2200 R_⊙, and $5\times {10}^{51}\,\mathrm{erg}\lesssim {E}_{\exp }\lesssim {10}^{52}\,\mathrm{erg}$. More detailed modeling should be performed to identify other possible degeneracies (e.g., between E_exp and M_ej, whether these can be broken by the Ni-powered portion of the light curve and to robustly estimate uncertainties in the inferred explosion parameters). In the future, it will be also beneficial to use realistic radiative transfer models during the mass ejection process so that the exact density profile of the CSM can be accurately determined. Spectral synthesis that distinguishes velocities of individual elements will also provide an alternative method to lift the degeneracy.

Dessart et al. (2020) showed that a star with no H-rich envelope explodes as a Type Ib (Ic) SN depending on its X(He) at the photosphere. According to Table 3 and Figure 3 in that work, a star with a photospheric He mass fraction X(He) ≲ 0.2 explodes as a Type Ic and X(He) ≳ 0.5 as a Type Ib SN. They also mentioned that ⁵⁶Ni mixing is less important because of the long mean free path of gamma rays. Our default model of M_ini = 25 M_⊙ contains X(He) ∼0.36 (Figure 5) and thus lies at the transition between Ib and Ic SNe. In higher-mass models in Section 4.4.6, X(He) = 0.16, 0.037, and 0.035 after mixing for M_ini = 40, 50, and 60 M_⊙, respectively. These models would appear as Type Ic SNe.

It is necessary to conduct detailed radiation hydrodynamical simulations to reliably determine the SN spectral type. Additionally, the mixing is likely nonspherical, so that ⁵⁶Ni and He may not be microscopically mixed. For example, for a high-energy jet-induced explosion (Nomoto 2017), much of the ejected mass that powers the light curve may originate near the CO core where energy deposition takes place. Such jet-associated ejecta and fallback materials effectively reproduce mixing (Figure 1 of Tominaga et al. 2007), and in this case we expect Type Ic spectra.

The relatively low ejecta mass of our default model is very small compared to the progenitor mass for PPISNe (i.e., progenitor masses M_ini of roughly 80–140 with pre-explosion masses M_exp of roughly 35–50 M_⊙). Hence, adopting the PPISN interpretation likely requires that the ejecta be highly aspherical (e.g., a bipolar jet-driven outflow) or that nearly all of the star collapses into a black hole and only a small fraction is ejected (see, e.g., Wongwathanarat et al. 2013 for aspherical mass ejection in the neutron star case.) Powell et al. (2021) examines core collapse and the (partial) explosion of PPISN models. For their low-mass PPISN model, a successful explosion of ∼3 × 10⁵¹ erg is observed but the high binding energy of the remnant implies black hole formation with partial mass ejection. This phenomenon could account for the low ejecta mass found by our light curve modeling. The escape of ⁵⁶Ni into the ejecta could be produced by this partial explosion, or in a jet or BH accretion disk wind.

Both of these possibilities may be realized in the collapsar scenario for driving the explosion, in which the explosion energy arises from a jet and/or disk wind from an accreting black hole (Woosley & Bloom 2006; Woosley & Heger 2006). Numerical simulations of collapsars show that the conservation of angular momentum leads to an accretion disk with a standing shock (e.g., Molteni et al. 1996; Stone et al. 1999). The disk can drive an outflow containing a couple of solar masses, a few tenths of a solar mass of ⁵⁶Ni, and an energy of ∼10⁵² erg (see, e.g., MacFadyen & Woosley 1999; Kohri et al. 2005; Zenati et al. 2020), similar to the properties of our default model. Alternatively, the explosion could be driven by a rapidly rotating magnetar (e.g., Obergaulinger & Aloy 2021). The jet (e.g., Maeda et al. 2002; Nagataki et al. 2003; Nomoto et al. 2013) can cause shock heating in the envelope, which creates a distinctive nucleosynthetic pattern including lots of ⁵⁶Ni (Tominaga 2009). The breakout can lead to a long gamma-ray burst (Tominaga et al. 2007; Nagataki 2011). However, whether the jet breaks out depends strongly on its energetics (Zhang et al. 2008; Barnes et al. 2018), with a smothered jet resulting in a Type Ic-BL event like SN 2018gep.

One-dimension models of massive stellar explosions (e.g., Heger & Woosley 2010; Nomoto 2017; Limongi & Chieffi 2020) predict that an explosion energy of roughly 2 × 10⁵¹ erg is necessary to produce ∼0.3 M_⊙ of ⁵⁶Ni in the ejecta, for progenitor models with masses M_ini ≲ 20 M_⊙. Higher explosion energies of roughly 5 × 10⁵¹ erg are required for higher-mass stars (Umeda & Nomoto 2008). The low-mass progenitor possibility is roughly consistent with our wave-driven models (though the explosion energy is lower than that inferred from light-curve modeling), while the high-mass possibility is more consistent with our PPISN models. In the low-mass progenitor scenario, the small ejecta mass inferred from light-curve modeling is expected from the ∼4−5 M_⊙ helium core of the progenitor star.

In the high-mass PPISN progenitor scenario, a small ejecta mass requires asymmetric explosion models in which most of the ⁵⁶Ni and explosion energy is contained within a small amount of mass, possibly ejected within bipolar jets. In this case, the ⁵⁶Ni considered in our light-curve models comes from the final explosion (and not from prior pulses). It is possible that, for strong pulsations, a significant amount of ⁵⁶Ni is synthesized prior to the final explosion. The radioactive energy can thermalize in the core and lead to further mass loss. However, this usually happens for high-mass PPISN near the pair-instability SN limit, whereas our light-curve modeling suggests lower-mass PPISN progenitors.

It will be interesting future work to self-consistently model the pre-SN evolution of a rapidly rotating helium star with pre-SN mass of ≈40 M_⊙. If it can produce pre-SN ejecta of ∼0.3 M_⊙ due to pulsational pair instability before collapse and generate a ∼10⁵² erg explosion with a few solar masses of high-velocity ejecta, it could likely produce an event like SN 2018gep. Recent work (Marchant & Moriya 2020) shows that whether pulsation occurs depends on how fast the star is rotating. Their rigidly rotating models show pulsations, which release less than half of the stellar angular momentum. This favors the formation of a collapsar, which allows jet-energy deposition for aspherical mass ejection. Future work should further examine the interplay of rotation, stellar evolution, mass ejection, and the SN explosion.

In this work we have relied on the stellar evolution code MESA for modeling the mass ejection and the formation of CSM before the final explosion. MESA assumes diffusive radiation transport, and it assumes that matter is in local thermal equilibrium with radiation. These approximations will gradually break down as the density of the ejected matter drops as it expands. The transition from optically thick to optically thin matter allows the radiation to leak out more efficiently, hence lowering the radiation pressure driving expansion. The exact temperature profile and opacity profile will also depend on this radiation transport, further complicating the problem and potentially altering the density profile away from that of our models.

However, we do not expect the CSM mass or radius to be strongly altered by more accurate radiative transport calculations. As shown in the wave-driven mass-loss model in Figure 4 of Fuller & Ro (2018), the ejected material has nearly reached its asymptotic velocity by the time that radiation diffuses effectively, so the expansion velocity calculated from MESA is a reasonable estimate. As long as the ejecta mass is optically thick where it is being accelerated (i.e., near the sonic point and escape point), the wave-driven ejecta speed tends to be near or somewhat larger than the escape speed, as described in Quataert et al. (2016). We suspect PPISN mass loss to also be driven in the optically thick limit so that the calculation of Renzo et al. (2020a) is a reasonable approximation. For accurate CSM modeling out to lower optical depths, multiband radiative transfer alongside hydrodynamics will be necessary. Nonradial instabilities between colliding shells will also affect the density profile (Chen et al. 2014), so multidimensional calculations should also be examined. However, this will greatly increase the required computational resource, so we leave these investigations for future work.

In our radiative transfer modeling, we approximated the CSM density profile with a 1/r² dependence. As shown in Figures 1–2, this scaling captures the global CSM features. However, a fully consistent model will require a more realistic density profile which also contains density jumps at shell-shell collisions. Our exploratory model of a thin and dense CSM shown in Section 4.4.6 suggests that some quantitative features of the CSM will change slightly with the CSM internal structure. It will be an interesting follow-up project to understand how these features affect the SN light curve.