Three-dimensional hydrodynamic simulations of supernova ejecta with a central energy source

TL;DR

This study uses 3D relativistic hydrodynamic simulations to explore how different levels of central energy injection influence supernova ejecta structure and observable features.

Contribution

It demonstrates how varying injected energies alter supernova ejecta morphology, highlighting potential observational differences linked to central engine activity.

Findings

High energy injection causes hot bubble breakout and clumpy ejecta.

Lower energy injection results in well-stratified outer envelopes.

Different energy levels may explain diverse supernova spectral features.

Abstract

We present the results of three-dimensional special relativistic hydrodynamic simulations of supernova ejecta with a powerful central energy source. We assume spherical supernova ejecta freely expanding with the initial kinetic energy of $10^{51}$ erg. We performed two simulations with different total injected energies of $10^{51}$ and $10^{52}$ erg to see how the total injected energy affects the subsequent evolution of the supernova ejecta. When the injected energy well exceeds the initial kinetic energy of the supernova ejecta, the hot bubble produced by the additional energy injection overwhelms and penetrates the whole supernova ejecta, resulting in clumpy density structure. For the smaller injected energy, on the other hand, the energy deposition stops before the hot bubble breakout occurs, leaving the outer envelope well-stratified. This qualitative difference may indicate that central engine powered supernovae could be observed as two different populations, such as supernovae with and without broad-line spectral features, depending on the amount of the total injected energy with respect to the initial kinetic energy.

Equations

$$\rho(r,t_{0})=\left\{\begin{array}[]{ccl}\rho_{0}\left(\frac{v}{v_{\mathrm{br}}}\right)^{-\delta}&\mathrm{for}&v\leq v_{\mathrm{br}},\\ \rho_{0}\left(\frac{v}{v_{\mathrm{br}}}\right)^{-m}&\mathrm{for}&v_{\mathrm{br}}<v\leq v_{\mathrm{max}},\end{array}\right.$$

$$v_{\mathrm{br}}=\sqrt{\frac{2f_5{E}_{\mathrm{sn}}}{f_3{M}_{\mathrm{ej}}}},$$

$$f_l=\frac{(m-l)(l-\delta )}{m-\delta -(l-\delta )(v_{\mathrm{br}}/v_{\max }{)}^{m-l}}.$$

$${\rho }_0=\frac{f_3{M}_{\mathrm{ej}}}{4\pi v_{\mathrm{br}}^3{t}_0^3}.$$

$$f_l=\frac{(m-l)(l-\delta )}{m-\delta }.$$

$$v_{\mathrm{br}}=\sqrt{\frac{2(m-5)(5-\delta )E_{\mathrm{sn}}}{(m-3)(3-\delta )M_{\mathrm{ej}}}}.$$

$$R_{\mathrm{in}}(t)=R_{\mathrm{in},0}\left(\frac{t}{t_0}\right),$$

$$R_{\mathrm{c}}={\xi }_{\mathrm{c}}A{t}^{\alpha },$$

$$A={\left\{\frac{3(\gamma -1)(2-\gamma )L}{{\alpha }^2\gamma \left[1+3\alpha (\gamma -1)\right]{\xi }_{\mathrm{c}}^3{\eta }_{\mathrm{c}}(v_{\mathrm{br}}{)}^{\delta -3}{f}_3{M}_{\mathrm{ej}}}\right\}}^{1/(5-\delta )},$$

$$E_{\mathrm{kin}}=E_{\mathrm{sn}}+\frac{(1+3\alpha )(\gamma -1)}{1+3\alpha (\gamma -1)}Lt,$$

$$E_{\mathrm{int}}=\frac{(2-\gamma )}{1+3\alpha (\gamma -1)}Lt,$$

$$f_{\mathrm{br}}=\frac{2{\alpha }^2\gamma \left[1+3\alpha (\gamma -1)\right]{\xi }_{\mathrm{c}}^3{\eta }_{\mathrm{c}}{f}_5}{3(\gamma -1)(2-\gamma )}=5.1,$$

$$q(r)=\frac{1}{\Delta V}{\int }_{r-\Delta r/2}^{r+\Delta r/2}{r}^{2}dr{\int }_{4\pi }d\Omega q(x,y,z),$$

$$M_{\mathrm{sw}}\simeq \frac{m-3}{m-\delta }{M}_{\mathrm{ej}},$$

$$v_{\mathrm{fs}}\simeq {\left[\frac{2(m-\delta )E_{\mathrm{inj}}}{(m-3)M_{\mathrm{ej}}}\right]}^{1/2},$$

$$E_{\mathrm{int}}(t)=\left\{\begin{array}[]{cl}\frac{(2-\gamma)}{1+3\alpha(\gamma-1)}L(t-t_{0})&\mathrm{for}\ t\leq t_{\mathrm{term}},\\ \frac{(2-\gamma)}{1+3\alpha(\gamma-1)}E_{\mathrm{inj}}\left(\frac{t}{t_{\mathrm{term}}}\right)^{-5/4}&\ \mathrm{for}\ t_{\mathrm{term}}<t.\end{array}\right.$$

$$E_{\mathrm{kin}}(t)=E_{\mathrm{sn}}+L(t-t_0)-E_{\mathrm{int}}(t).$$

$$R_{\mathrm{fs},\mathrm{f}}{\left(\frac{t}{20t_{\mathrm{c}}}\right)}^{5/4}=v_{\max }t,$$

$$t=\frac{v_{\max }^4(20t_{\mathrm{c}}{)}^5}{R_{\mathrm{fs},\mathrm{f}}^4}\simeq 2\times 10^2\mathrm{yrs}.$$

$$\rho (t,r)=B{t}^{m-3}{r}^{-m}.$$

$$R_{\mathrm{fs}}(t)=C{t}^{\lambda },$$

$$V_{\mathrm{b}}=\frac{4\pi }{3}{R}_{\mathrm{fs}}^3=\frac{4\pi }{3}{C}^3{t}^{3\lambda }.$$

$$\rho (t,R_{\mathrm{fs}})=B{C}^{-m}{t}^{-\lambda m+m-3}.$$

$$p_{\mathrm{f}}=\frac{2}{\gamma +1}\rho (t,R_{\mathrm{fs}}){\left(\frac{{dR}_{\mathrm{fs}}}{dt}-\frac{R_{\mathrm{fs}}}{t}\right)}^2=\frac{2(\lambda -1{)}^2}{\gamma +1}B{C}^{2-m}{t}^{(2-m)\lambda +m-5}.$$

$$E_{\mathrm{int}}\propto V_{\mathrm{b}}^{1-\gamma }\propto t^{3\lambda (1-\gamma )}.$$

$$p_{\mathrm{b}}=(\gamma -1)\frac{E_{\mathrm{int}}}{V_{\mathrm{b}}}\propto V_{\mathrm{b}}^{-\gamma }\propto R_{\mathrm{fs}}^{-3\gamma }\propto t^{-3\lambda \gamma }.$$

$$(2-m)\lambda +m-5=-3\lambda \gamma ,$$

$$\lambda =\frac{m-5}{m-3\gamma -2}.$$

Full text

supernova: general – gamma-ray burst: general – shock waves

1 INTRODUCTION

There is growing evidence that several classes of unusual core-collapse supernovae (CCSNe) harbor a powerful “engine” at the center. Broad-lined Ic SNe (SNe Ic-BL) are characterized by broad absorption features in their spectra, indicating large explosion energies. It is widely known that some SNe Ic-BL are associated with long-duration gamma-ray bursts (GRBs) and thus require an ultra-relativistic jet penetrating the star (Woosley & Bloom, 2006; Hjorth & Bloom, 2012; Cano et al., 2017). These energetic CCSNe are also recognized as bright radio emitters, which strongly indicates the presence of a fast ejecta component and a (mildly) relativistic blast wave giving rise to radio synchrotron emission (e.g., Kulkarni et al., 1998; Soderberg et al., 2010). In addition, one of the most plausible scenarios for the emerging new class of hydrogen-poor superluminous SNe (also known as type I superluminous SNe, hereafter SLSNe-I; Quimby et al. 2007; Barbary et al. 2009; Pastorello et al. 2010; Quimby et al. 2011; Chomiuk et al. 2011) is the so-called central engine scenario, in which the central compact object, a highly rotating magnetized neutron star (Kasen & Bildsten, 2010; Woosley, 2010) or a black hole accretion disk (Dexter & Kasen, 2013), injects an additional energy into the surrounding supernova ejecta to give rise to bright thermal emission (see, e.g., Quimby et al., 2011; Gal-Yam, 2012; Moriya et al., 2018, for review). Despite a lot of investigations on the effects of the putative central engine left in SN ejecta, there are still many problems remained unanswered, such as, how it is produced, how exactly it deposits an additional energy into the surrounding SN ejecta, and what kind of conditions should be met for their progenitor systems.

Recently, SLSNe and their origin have been paid great attention. Despite their small occurrence rate (e.g. Quimby et al., 2013; McCrum et al., 2015; Prajs et al., 2017), their extreme brightness in optical and UV bands makes it feasible to find them even in distant galaxies (Cooke et al., 2012; Tanaka et al., 2012, 2013; Inserra & Smartt, 2014; Inserra et al., 2018a). The recent discoveries of spectroscopically confirmed SLSNe at $z\sim 2$ by Dark Energy Survey Supernova Program (Pan et al., 2017; Smith et al., 2018) and Subaru/Hyper Sprime-Cam (HSC) (Moriya et al., 2019; Curtin et al., 2019) indeed demonstrated their potential as a probe of star-forming activity in high-z universe.

However, the power source(s) of the bright emission of SLSNe-I is still unclear. In addition to the central engine scenario introduced above, there are two competing scenarios; the interaction between the circumstellar medium and the SN ejecta (CSM scenario; Chevalier & Irwin 2011; Ginzburg & Balberg 2012; Moriya et al. 2013), and pair-instability SNe (PISN scenario; Barkat et al. 1967; Rakavy & Shaviv 1967; Heger & Woosley 2002; Woosley et al. 2007; Gal-Yam et al. 2009). These three scenarios have both advantages and disadvantages in explaining the whole populations of SLSNe-I. Therefore, whether the observed population of SLSNe-I can further be divided into two or more subclasses and whether SLSNe-I share some common properties with other types, are important questions ever since SLSNe were recognized as a distinguishing population. Gal-Yam (2012) proposed the presence of the so-called SLSNe-R subclass; SLSNe-I whose late-time light curves can be accounted for by radioactive 56Ni. The required amount of 56Ni is huge, $1$–$10\ M_{\odot}$, and therefore is unlikely synthesized in normal CCSNe, but may be in line with very massive bare carbon-oxygen cores invoked in the PISN scenario. However, the PISN scenario is also known to have difficulty in explaining some rapidly evolving SLSNe-I because theoretical PISN models predict much longer photon diffusion times (Nicholl et al., 2013). Furthermore, they do not succeed in reproducing blue spectra of SLSNe-I (e.g., Dessart et al., 2012; Jerkstrand et al., 2016). In addition, there are some observational studies suggesting the presence of fast- and slowly-evolving populations of SLSNe-I (e.g., Inserra et al. 2013, 2017, see also Nicholl et al. 2015, 2017; Kangas et al. 2017), which might indicate two distinct channels or progenitors for SLSNe-I.

Many authors have provided model light curves of SLSNe-I (Inserra et al., 2013; Chatzopoulos et al., 2013; Nicholl et al., 2015; Wang et al., 2015a, b, 2016, 2017; Nicholl et al., 2017), mostly based on the standard treatment of thermal photon diffusion in spherically expanding ejecta (Arnett’s solution; see, Arnett 1980, 1982, 1996). However, photometric data alone appear to be difficult to distinguish one scenario from the others. Recently, thanks to the currently ongoing transient survey programs, such as the Zwicky Transient Facility (ZTF; Bellm, 2014), Pan-STARRS (Kaiser et al., 2002), and HSC Subaru Strategic Program, there are plenty of observational data of SLSNe-I accumulated (e.g., De Cia et al., 2018; Lunnan et al., 2018). With these rich photometric and spectroscopic data sets, some authors have performed analyses with more sophisticated statistical approaches (e.g., Liu et al., 2017; Nicholl et al., 2017). Most recently, Inserra et al. (2018b) developed a statistical method to distinguish and classify SLSNe-I by using photometric datasets combined with Gaussian process regression. They claim that there are two distinct populations characterized by fast and slow evolutions. Quimby et al. (2018) have also developed their own method for distinguishing SLSNe-I from other types by systematically analyzing light curves and spectra of SNe from the Palomar Transient Factory (PTF) archival data. They suggest that their SLSN-I samples can be divided into two groups, PTF12dam-like and SN2011ke-like objects, according to their spectral similarities. The latter objects show smoother spectra around the optical maximum than the former objects. Nicholl et al. (2019) performed a principal component analysis of the nebular spectra of SLSN-I samples and found no firm evidence of sub-populations.

On the other hand, SLSNe-I are known to share some common properties with SNe Ic-BL and GRB-SNe in terms of spectroscopic features (Pastorello et al., 2010; Nicholl et al., 2016; Jerkstrand et al., 2017; Liu et al., 2017) and host galaxies (Neill et al., 2011; Lunnan et al., 2014; Leloudas et al., 2015a; Thöne et al., 2015; Angus et al., 2016; Perley et al., 2016; Chen et al., 2017; Schulze et al., 2018; Hatsukade et al., 2018). Moreover, the luminous SN component superposed on the afterglow emission of the ultra-long GRB 111209A, dubbed SN 2011kl (Greiner et al., 2015; Kann et al., 2019), was brighter than any other SNe associated with GRBs and was remarkably similar to SLSNe-I. These similarities may further support the idea that at least some SLSNe-I are actually powered by a central energy source and motivate some authors to put forward unified scenarios for energetic and bright CCSNe (e.g., Metzger et al., 2015; Kashiyama et al., 2016). If both SLSNe-I and SNe Ic-BL are powered by a central engine, the engine should play multiple roles. In other words, it sometimes powers bright optical emission lasting up to hundreds of days, while it produces energetic ejecta with large kinetic energies or even launches an ultra-relativistic jet. It is still unclear whether a single scenario covers these distinct time scales. On one hand, the two different energy deposition mechanisms are realized in a single object with different settings. For example, Margalit et al. (2018b) suggest that the inclination of the dipole magnetic field with respect to the rotational axis of a rotating magnetized proto-neutron star may lead to both GRB-like jet and quasi-spherical wind. On the other hand, this diversity may be explained by physically distinct central engines, neutron stars for a population and black holes for the other. Currently, what the central engine is and how exactly it deposits energy into the supernova ejecta are still unclear.

How the deposited energy affects the dynamical evolution of supernova ejecta is of critical importance in observations of CCSNe. Properties of supernova ejecta with a central energy source have been investigated mainly in the context of Galactic pulsar wind nebula embedded in supernova remnants (e.g., Chevalier 1984; Chevalier & Fransson 1992; Jun 1998, see Gaensler & Slane 2006, for a review). Recently, this issue attracts further attention in the context of the central engine model for SLSNe-I. However, it is only recently that multi-dimensional numerical simulations are performed in this context (Chen et al., 2016; Suzuki & Maeda, 2017; Blondin & Chevalier, 2017). In our previous study (Suzuki & Maeda 2017; hereafter, SM17), we performed a long-term hydrodynamic simulation of SN ejecta with a relativistic wind injected from the center in 2D cylindrical geometry. The simulation revealed that hydrodynamic instabilities caused by the ejecta-wind interaction play a vital role in determining the resultant ejecta structure. However, the assumed axisymmetry resulted in an artificial bipolar structure in the density distribution of the ejecta in the free expansion stage. Therefore, in this work, we perform simulations with similar numerical settings in 3D cartesian coordinates. We find that the dynamical evolution of the ejecta in the 3D simulation is qualitatively similar to the previous 2D simulation and find that the analytical consideration in SM17 can also apply to the ejecta in 3D. Furthermore, we also consider the case in which the total amount of the injected energy is comparable to the original kinetic energy of the ejecta.

This paper is organized as follows. In Section 2, we describe numerical setups. We perform a couple of numerical simulations with different amounts of the injected energy. Results of the numerical simulations are presented and compared in Section 3. We discuss some implications of the numerical results on properties of SN ejecta with a central energy source in Section 4. Finally, Section 5 concludes this paper.

2 Numerical Setups

The numerical setups of the 3D simulation are similar to SM17 for the purpose of comparison. In the following, we briefly describe the assumptions and conditions of our 3D simulations.

2.1 Supernova ejecta

We assume spherical SN ejecta with the mass $M_{\mathrm{ej}}=10M_{\odot}$ and the kinetic energy $E_{\mathrm{sn}}=10^{51}$ erg. The radial velocity $v$ of a layer is given by its radius $r$ divided by the elapsed time $t$, $v=r/t$ for $v<v_{\mathrm{max}}$. The evolution of the ejecta starts at $t=t_{0}$. The radial density profile $\rho(r,t)$ at $t=t_{0}$ is described by a broken power-law function,

[TABLE]

(Chevalier & Soker, 1989; Matzner & McKee, 1999). The velocity at the break $v_{\mathrm{br}}$ separates the inner and outer parts of the ejecta. The non-dimensional parameters $\delta$ and $m$ specify the inner and outer density slopes. The break velocity is expressed in terms of the ejecta mass and the kinetic energy as follows;

[TABLE]

with

[TABLE]

The density normalization constant $\rho_{0}$ is given by

[TABLE]

For centrally concentrated ejecta with a steep outer slope (i.e., a small $v_{\mathrm{br}}/v_{\mathrm{max}}$ and a large $m$) as assumed in our simulations, the numerical factor $f_{l}$ can be approximated as,

[TABLE]

Accordingly, Equation (2) yields,

[TABLE]

We assume the same inner and outer slopes as SM17, $\delta=1$ and $m=10$. The break velocity yields $v_{\mathrm{br}}=3.8\times 10^{8}$ cm s*-1* with this parameter set. We further assume that the maximum velocity is 10 times larger than the break velocity, $v_{\mathrm{max}}=3.8\times 10^{9}$ cm s*-1*.

2.2 Energy injection

The energy injection at the center of the ejecta is realized as follows. We specify a small spherical region at the center of the numerical domain, in which we deposit an additional energy. In the energy injection region, thermal gas is injected at a constant energy injection rate $L$ and then the injection is terminated at $t=t_{\mathrm{term}}$, when the total injected energy reaches $E_{\mathrm{inj}}=L(t_{\mathrm{term}}-t_{0})$. The corresponding mass injection rate is denoted by $\dot{M}c^{2}=L/\Gamma_{\mathrm{cr}}$, where $\Gamma_{\mathrm{cr}}$ is a parameter governing the baryon richness of the injected gas and set to $\Gamma_{\mathrm{cr}}=20$. This value indicates that the rest mass energy of the injected gas is smaller than the total injected energy by a factor of $20$. This baryon richness is so small that the injected gas is relativistic. In the previous 2D simulation (SM17), we adopted a fixed energy injection region. In 3D simulations, however, it is numerically expensive to keep resolving the fixed energy injection region as its radius progressively becomes small compared with the physical scale of the ejecta. Therefore, we adopt an expanding energy injection region. The energy injection radius linearly increases with time according to the free expansion,

[TABLE]

where $R_{\mathrm{in,0}}=10^{11}$ cm is the initial values of the energy injection radius. The energy injection rate is set to $L=10^{46}$ erg s*-1*, following SM17.

We carry out a couple of simulations with different amounts of the injected energy. In model E52, we continue the energy injection until the total amount of the injection energy reaches $E_{\mathrm{inj}}=10^{52}$ erg at $t=t_{0}+10^{6}$ s. This model is the 3D counterpart of our previous 2D cylindrical simulation. We perform another simulation with the total injected energy of $E_{\mathrm{inj}}=10^{51}$ erg at $t=t_{0}+10^{5}$ s (referred to as model E51), while the energy injection rate is same as model E52. In this model, the injected energy is comparable to the initial kinetic energy of the supernova ejecta. In the following, we compare these two models to see how the dynamical evolution of the supernova ejecta depends on the amount of the injected energy.

It is convenient to normalize the elapsed time by the characteristic time scale of the energy injection, $t_{\mathrm{c}}=E_{\mathrm{sn}}/L=10^{5}$ s. We begin our simulations at $t=t_{\mathrm{0}}=0.02t_{\mathrm{c}}$ and evolve the system up to $t=20t_{\mathrm{c}}$.

2.3 Adaptive mesh refinement and coarsening

Our special relativistic hydrodynamics code is equipped with an adaptive mesh refinement technique (AMR; Berger & Colella 1989), which automatically covers regions requiring finer resolutions by smaller grids. In particular, shock waves driven by the outermost layer of the ejecta and the hot bubble created by the energy injection are captured by numerical cells with relatively fine resolutions. However, tracking the propagation of the shock waves with a fixed resolution requires an impractically large number of numerical cells to cover their structures, especially, in 3D simulations. Therefore, in order to ensure realistic computational costs, we eventually relax the numerical resolution by gradually reducing the maximum level of the AMR grid. As a result, the minimum resolved length increases with time. Nevertheless, it is always kept small compared to the physical scale of the entire ejecta. The base grid, i.e., the AMR level $l=0$, covers the numerical domain of $-4.8\times 10^{16}\mathrm{cm}\leq x,y,z\leq 4.8\times 10^{16}\mathrm{cm}$ by $64^{3}$ numerical cells. The maximum AMR level is initially set to $l=16$, corresponding to the numerical domain effectively covered by $\sim(4.2\times 10^{6})^{3}$ cells. Despite the slightly different treatments of the energy injection and the coarsening of the numerical grid, we obtain qualitatively similar results to the 2D counterpart as we shall see below.

3 Results

In the following, we present the results of the numerical simulations.

3.1 Model E52

3.1.1 Hot bubble expansion and breakout

We first focus on the results of model E52. The spatial distributions of the 4-velocity $\Gamma\beta$, the density $\rho$, and the pressure $p$, in the 2D slices around $y=0$ are shown in Figure 1. The panels in Figure 1 represent snapshots of the simulation at $t/t_{\mathrm{c}}=3.0$, $4.0$. $4.5$, and $5.0$, during which the forward shock emerges from the surface of the ejecta. The dynamical evolution of the ejecta in model E52 is qualitatively similar to that of our previous 2D cylindrical simulation (SM17) and calculations by other authors (Chen et al., 2016; Blondin & Chevalier, 2017).

At the beginning of the evolution, a quasi-spherical region with high pressure (hot bubble) appears around the center of the ejecta due to the energy injection. The forward shock driven by the pressure of the hot bubble propagates in the supernova ejecta. When the forward shock front is still in the inner part of the supernova ejecta, the ram pressure exerting on the post-shock gas is sufficiently large to contain the thermal pressure of the hot bubble, keeping the hot bubble nearly spherical. However, once the forward shock front passes through the interface separating the inner and outer parts of the ejecta, $r=v_{\mathrm{br}}t$, the ram pressure of the outer ejecta can no longer contain the hot bubble, leading to an efficient acceleration of the forward shock. On the course of the hot bubble expansion, the Rayleigh-Taylor instability develops around the contact surface, mixing the injected gas with the supernova ejecta. The forward shock modified by the development of the Rayleigh-Taylor instability efficiently expands in the outermost layer of the supernova ejecta, resulting in the violent breakout of the shock into the surrounding space. This hot bubble breakout plays an important role in producing clumpy ejecta, which is distinguished from the 1D spherical picture of supernova ejecta with central energy injection.

In our 3D simulation, the presence of the hot bubble and its breakout are successfully reproduced, as is seen in Figure 1. In our previous 2D cylindrical simulation, shocked region and the density distribution exhibit global bipolar structure despite the energy injection in a spherical manner. As we have noted in SM17, this is an artificial effect due to the assumed axisymmetry in the 2D cylindrical simulation. The 3D results presented in Figure 1 do not suffer from such an artificial effect and exhibit globally spherical ejecta after the hot bubble breakout. However, one caveat should be noted. In the initial stage of the development of the Rayleigh-Taylor instability, the shock front exhibits an artificial symmetry. In this simulation, numerical errors arising from the 3D cartesian grid structure predominantly contribute to the seed of the instability. Therefore, the angular distributions of hydrodynamic variables should be treated carefully. This issue could be remedied by simulations with much higher resolution. As we shall see below, however, the radial distributions of hydrodynamic variables after the hot bubble breakout can be understood in the same analytical way as the 2D cylindrical study, and thus our conclusions are unlikely affected by this numerical issue.

Figure 2 shows the growth of the hot bubble in a more quantitative way. We identify the numerical cell with the highest density (referred to as the density peak) in the computational domain and plot the distance $R_{\mathrm{peak}}$ of the density peak from the center as a function of time in the upper panel of Figure 2. The growth of $R_{\mathrm{peak}}$ is expected to follow that of the contact surface separating the shocked injected gas and the surrounding ejecta. Its temporal evolution, which is shown as a solid line, is compared with a self-similar solution, which has been developed for SN remnants with a pulsar wind nebula in its center (Chevalier, 1984; Chevalier & Fransson, 1992; Jun, 1998). The radius $R_{\mathrm{c}}$ of the contact surface is given by

[TABLE]

with $\alpha=5/4$ and $\xi_{\mathrm{c}}=0.9849$ for the parameter set adopted in this calculation. The constant $A$ is given by

[TABLE]

with $\eta_{\mathrm{c}}=0.1706$ (see, SM17 for the derivation of the self-similar solution). In the upper panel of Figure 2, we also plot the temporal evolution of the minimum resolved length $\Delta x_{\mathrm{min}}$. In the lower panel of Figure 2, we plot the kinetic and internal energies, $E_{\mathrm{kin}}$ and $E_{\mathrm{int}}$, as a function of time. From the self-similar solution, they are expected to grow as follows,

[TABLE]

for the kinetic energy, and

[TABLE]

for the internal energy.

The temporal evolution of $R_{\mathrm{peak}}$ agrees with that of $R_{\mathrm{c}}$, Equation (8), up to $t/t_{\mathrm{c}}\simeq 4.0$. The growth of the kinetic and internal energies are also well described by Equations (10) and (11) in the same time interval. This means that the dynamical evolution of the ejecta immediately approaches the self-similar solution, thereby suggesting that the results of the simulations are less dependent on the value of the initial time $t_{\mathrm{0}}$ and the initial condition. The radius $R_{\mathrm{peak}}$ starts deviating from the self-similar solution after $t/t_{\mathrm{c}}\simeq 3.0$–$4.0$. This time corresponds to the breakout of the hot bubble from the outermost layer of the ejecta (see, Figure 1). From the analytical considerations, the breakout is expected to occur when the forward shock driven by the hot bubble reaches the interface separating the inner and the outer parts of the ejecta, at $r=v_{\mathrm{br}}t$. The self-similar solution predicts that this occurs at $t=f_{\mathrm{br}}t_{\mathrm{c}}$, where the numerical factor is given by

[TABLE]

for the adopted parameter set. The hot bubble breakout occurs slightly earlier in the numerical simulation than the prediction by the 1D self-similar solution. This is probably due to the shock front modified by the Rayleigh-Taylor instability. Because of the non-linear development of the Rayleigh-Taylor instability, shocks propagating into some radial directions can “overshoot” with respect to the average radius of the shock front. Such overshooting components emerge from the inner ejecta earlier than shocks along other directions, which makes the hot bubble breakout happening earlier. This qualitatively explains the earlier hot bubble breakout. However, the development of the Rayleigh-Taylor instability before the hot bubble breakout is highly non-linear, which makes it difficult to predict the exact onset time of the hot bubble breakout.

3.1.2 Density structure

Figure 3 shows the density distribution at the end of the simulation, $t/t_{\mathrm{c}}=10$. The hot bubble breakout and a bunch of relativistic outflows penetrating the ejecta realize patchy density structure. The outermost layer of the supernova ejecta is accelerated by shocks driven by the hot bubble and their velocities are close to the speed of light. Since the central energy injection is terminated at $t/t_{\mathrm{c}}=10.02$, the outer part of the ejecta is already expanding almost freely.

Despite the patchy density structure, physical variables exhibit simple global distributions. We obtain the radial profile of a physical variable $q$ by dividing the numerical domain into successive concentric shells with a small width $\Delta r$ and averaging the variable over the shell,

[TABLE]

where $\Delta V=4\pi r^{2}\Delta r$ is the volume of the concentric shell and $\Omega$ is the solid angle. In Figure 4, we plot the radial profiles of the density, the radial velocity $v_{r}$, and the kinetic luminosity at $t/t_{\mathrm{c}}=20$. The kinetic luminosity at $r$ is defined as the energy going through the surface area $4\pi r^{2}$ at $r$ per unit time and is calculated by the product of the non-relativistic kinetic energy flux $\rho v_{r}^{3}/2$ and the surface area $4\pi r^{2}$.

The radial velocity profile follows $v_{r}\propto r$, which is naturally expected for freely expanding ejecta. The kinetic luminosity shown in the bottom panel of Figure 4 exhibits a flat profile extending out to layers traveling at relativistic speeds, $\Gamma v_{r}\sim c$. This is understood as a result of the central energy injection at a constant rate and the efficient energy transport throughout the ejecta. As discussed in SM17, the channel flows emanating the ejecta play a role in transporting the injected energy and depositing it into the outer layers of the ejecta, resulting in a flat kinetic energy distribution. This is a distinguishing feature of the multi-dimensional picture of supernova ejecta with a long-lived central energy source. In spherically symmetric models, the transport of the injected energy is realized by the spherical forward shock driven by the pressure of the hot bubble, ending up as a spherical shell made up of the shocked ejecta (Kasen & Bildsten, 2010).

SM17 have considered the radial density profile realized when freely expanding spherical ejecta follow a flat kinetic energy distribution, and analytically derived a single power-law density profile, $\rho\propto t^{-3}v_{r}^{-n}$ or $\rho\propto t^{n-3}r^{-n}$, with an exponent in the range of $5\leq n\leq 6$. In the top panel of Figure 4, we compare the radial density profile obtained by our hydrodynamic simulation with single power-law functions with exponents $-5$ and $-6$. The comparison shows a good agreement, demonstrating the validity of the analytical consideration.

3.2 Model E51

3.2.1 Dynamical evolution of ejecta with moderate energy injection

The dynamical evolution of model E51 shows a marked difference from its more energetic counterpart. Figure 5 shows spatial distributions of physical variables realized in model E51 after the central energy injection is terminated. Although the hot bubble can be seen in the pressure distribution in Figure 5, it is still confined by the unshocked ejecta in a quasi-spherical manner. Even at the end of the simulation, $t/t_{\mathrm{c}}=20.0$ (the bottom panels of Figure 5), the forward shock front is far from the outermost layer of the ejecta. This is clearly due to the smaller injected energy than model E52.

This outcome can be understood again by the self-similar solution described by Equations (8) and (9). In fact, the total amount of the injected energy $E_{\mathrm{inj}}$ relative to the initial kinetic energy $E_{\mathrm{sn}}$ of the supernova ejecta plays a critical role in determining whether the forward shock reaches the interface before the energy injection is terminated. When the energy injection still continues after the forward shock reaches the interface, the continuously supplied energy drives the violent breakout of the hot bubble as we have described in Section 3.1. However, without the on-going energy supply, the hot bubble appears to be stuck in the ejecta even after the forward shock reaches the interface. For the adopted density profile, the threshold energy of $5.1E_{\mathrm{sn}}\simeq 5\times 10^{51}$ erg (see, SM17) is between $10^{51}$ and $10^{52}$ erg, which reasonably explains the different dynamical evolutions of models E52 and E51. Figure 6 presents the temporal evolutions of $R_{\mathrm{peak}}$, $\Delta x_{\mathrm{min}}$, and the kinetic and the internal energies. In contrast to model E52, the termination of the central energy injection at $t/t_{\mathrm{c}}=1.0$ makes the internal energy suddenly decline at this time.

3.2.2 Self-similar expansion of the hot bubble

When the injected energy is insufficient to realize the hot bubble breakout, the hot bubble is confined in the ejecta even after the forward shock front reaches the outer part of the ejecta with the steep density gradient. Figure 7 shows the radial distributions of the density, the radial four-velocity and the pressure at $t/t_{\mathrm{c}}=1.0$, $2.0$, $5.0$, $10.0$, and $20.0$ for model E51. The radial profiles are composed of the unshocked ejecta and the embedded hot bubble. As is seen in the radial profiles of the pressure, the hot bubble is clearly separated from the unshocked ejecta by the forward shock front. The density profiles exhibit peaks around the interface between the hot bubble and the surrounding ejecta. One of the remarkable features is the almost uniform pressure distribution within the hot bubble. This suggests that the time scale of the shock propagation into the surrounding ejecta is sufficiently longer than the sound crossing time within the hot bubble, i.e., the nearly entire region of the hot bubble is causally connected. Therefore, the pressure equilibrium is achieved inside the hot bubble.

Furthermore, the profiles of physical variables at different epochs are similar to each other, indicating that the hot bubble and the surrounding ejecta evolve in a self-similar way, although the contact surface has already suffered from the development of the Rayleigh-Taylor instability. This suggests that the evolution is described by another self-similar solution, which corresponds to a point-like explosion in a spherical expanding medium. In Appendix A, we consider the self-similar expansion of a spherical hot region in a freely expanding medium with a power-law radial density profile, by assuming that the pressure equilibrium is realized in the entire shocked region. In deriving the self-similar expansion law, we have assumed that the hot bubble expands in an adiabatic way (equation A7). We found that; (i) in a medium with a shallow radial density slope, $-d\ln\rho/d\ln r\leq 3\gamma+2$, the hot bubble is confined in the ejecta and passively expands according to the freely expanding ambient gas, and (ii) in a medium with a steep radial density slope, $-d\ln\rho/d\ln r>3\gamma+2$, the forward shock expands outwardly with respective to the expanding ambient gas. This finding suggests that the hot bubble is stuck in the inner ejecta when its outer radius does not reach the interface separating the inner and outer ejecta, $r=v_{\mathrm{br}}t$. On the other hand, once the hot bubble emerges from the interface, the latter self-similar expansion law for steep radial density slopes can apply.

In fact, the injected energy relative to the initial ejecta kinetic energy, $E_{\mathrm{inj}}/E_{\mathrm{sn}}$, again plays a role in determining which case should be realized. In order for the outer radius of the hot bubble to catch up with the interface, the forward shock velocity $v_{\mathrm{fs}}$ relative to the velocity of the ambient ejecta should be larger than or the same order of magnitude as the break velocity $v_{\mathrm{br}}$. When the hot bubble sweeps a mass $M_{\mathrm{sw}}$, the forward shock velocity can roughly be estimated to be $v_{\mathrm{fs}}\sim(2E_{\mathrm{inj}}/M_{\mathrm{sw}})^{1/2}$. Setting the swept mass to be the mass of the inner ejecta,

[TABLE]

the forward shock velocity is found to be

[TABLE]

when it reaches the interface. For $E_{\mathrm{inj}}\simeq E_{\mathrm{sn}}$, which is of particular interest here, this forward shock velocity is the same order of magnitude as $v_{\mathrm{br}}$. Therefore, in this particular model E51, the hot bubble can reach the transition layer $r=v_{\mathrm{br}}t$. As is shown in Figures 5 and 7, the hot bubble certainly reaches the outer ejecta at later epochs. In addition, in the bottom panel of Figure 7, we plot the dashed line representing the scaling relation $p\propto r^{-3\gamma}$ expected for the adiabatic expansion, and confirm the agreement with the declining hot bubble pressure. Therefore, the late-time behavior of the hot bubble in the numerical simulation is well explained by the latter case. On the other hand, for a sufficiently small injected energy $E_{\mathrm{inj}}\ll E_{\mathrm{sn}}$, the hot bubble would never reach the interface and end up being dissolved in the inner ejecta.

To summarize, the dynamical evolution of the hot bubble in model E51 first follows the self-similar solution discussed in Section 3.1.1, and then makes a transition to the other self-similar evolution discussed in Appendix A.

Due to the terminated energy injection, the internal energy of the ejecta decreases after the linear increase described by Equation (11). While the forward shock is propagating in the outer part of the ejecta, the self-similar expansion law expects the internal energy decreasing as $E_{\mathrm{int}}\propto t^{-5/4}$. Thus, we obtain the following temporal evolution of the internal energy by connecting the two different self-similar laws,

[TABLE]

Accordingly, the kinetic energy evolves as,

[TABLE]

In the bottom panel of Figure 6, we compare the temporal evolutions of the kinetic and internal energies in the numerical simulation with the analytical evaluation, which shows an overall agreement.

We estimate the time when the forward shock front catches up with the outermost layer of the ejecta. At the end of the simulation, $t/t_{\mathrm{c}}=20$, the forward shock is located around $R_{\mathrm{fs,f}}=10^{15}$ cm, which continues to grow as $R_{\mathrm{fs}}\propto t^{5/4}$. Thus, the forward shock front reaches the surface of the ejecta when the following condition,

[TABLE]

is satisfied, where we have neglected the interaction between the outermost layer and the surrounding medium for simplicity. For the maximum velocity of $v_{\mathrm{max}}\simeq 3.8\times 10^{9}$ cm s*-1*, the expected “touch-down” time is found to be

[TABLE]

The outermost layer of the ejecta is eventually decelerated by the collision with the ambient medium, which makes the touch-down happen earlier to some extent. Nevertheless, this time scale is clearly much longer than the diffusion time scale for thermal photons in the ejecta. We also note that the diffusion time scale shorter than the touch-down time violates the assumption of the adiabatic expansion in deriving the self-similar expansion law. When approximately taking into account the effective radiative loss, where the adiabatic exponent of the post-shock gas is approximately close to unity, the hot bubble is stuck in the ejecta as we discuss in Appendix A.

Finally, in Figure 8, we compare the radial distributions of the density, the radial 4-velocity, and the pressure at $t=20t_{\mathrm{c}}$ for models E52 and E51. We again see the differences between the ejecta with and without the hot bubble breakout. We also plot the radial distributions of these variables for the freely expanding ejecta without the central energy injection nor the interaction with the ambient medium. For model E52, the ejecta is clearly accelerated and thus the outer layer of the ejecta is well ahead of the free expansion case. For model E51, while the inner part of the ejecta is significantly affected by the energy injection, the outer part exhibits similar radial profiles to the freely expanding case.

4 Discussion

Through the 3D hydrodynamical simulations, we have shown that supernova ejecta with a long-lived central energy source follow different dynamical evolutions depending on the total amount of the injected energy. In the following, we summarize the dynamical evolution of supernova ejecta with a central energy source and discuss their observational implications.

4.1 Supernova ejecta with hot bubble breakout

First, we consider the evolution of supernova ejecta with a sufficiently large injection energy for the hot bubble breakout. Supernova ejecta with relatively large kinetic energies would have an advantage in explaining some energetic and/or bright supernovae. Several energetic SNe Ic-BL associated with GRBs, e.g., SNe 1998bw (Kulkarni et al., 1998; Galama et al., 1998), 2006aj (Campana et al., 2006; Mazzali et al., 2006; Modjaz et al., 2006; Pian et al., 2006; Soderberg et al., 2006), and 2010bh (Starling et al., 2011; Cano et al., 2011; Bufano et al., 2012; Olivares E. et al., 2012) are believed to be driven by a powerful central engine. Even without any associated gamma-ray signal explicitly manifesting central engine activities, radio-bright SNe Ic-BL sharing similar properties with GRB-SNe, such as SN 2009bb and 2012ap, are also suspected to be an explosion driven by a central engine (Soderberg et al., 2010; Milisavljevic et al., 2015; Chakraborti et al., 2015). In addition, some authors pointed out that the energetic type Ib SN 2012au also share common emission properties with energetic SNe Ic-BL (Milisavljevic et al., 2013; Takaki et al., 2013; Milisavljevic et al., 2018). The observational properties of these energetic SNe are worth comparing with the multi-dimensional picture of supernova ejecta with a powerful central energy source revealed by the present simulation.

4.1.1 Photometric and spectral evolution

Our 3D hydrodynamic simulation turns out to be qualitatively similar to the 2D cylindrical counterpart, confirming the multi-dimensional picture of supernova ejecta with a powerful central energy source. As has already been suggested by several work (Chen et al. 2016; SM17), one of the important outcomes of the hot bubble breakout is the efficient mixing of materials out to the surface of the ejecta. The supernova ejecta penetrated by channel flows likely show spatially inhomogeneous or clumpy distributions of elements having synthesized in the steady and explosive burning stages prior to the energy injection. Therefore, heavy elements supposed to be embedded in deeper layers for normal supernova ejecta could be found in both inner and outer layers for supernova ejecta with hot bubble breakout. These mixed heavy elements would contribute to opacities in outer layers of the ejecta, affecting the photometric and spectroscopic properties of SNe. Radioactive nickel mixed into outer layers of SN ejecta is also known to affect light curves of SNe. Recent systematic analysis of SNe Ic-BL from the PTF (Taddia et al., 2019) found that all of them require nickel mixing to some extent in order to explain their light curves. The spectroscopic observations of GRB 171205A/SN 2017iuk at unprecedentedly early epochs (Izzo et al., 2019) recently revealed the presence of heavy metals in the fast ejecta component preceding the SN ejecta. These are in agreement with the elemental mixing due to a long-lived central engine.

Another important consequence of the hot bubble breakout is relatively flat radial density distributions. The radial density structure is of fundamental importance in spectral formation. As we have seen in Section 3.1, the radial density profile of the ejecta is well described by a power-law function with an exponent between $-5$ and $-6$. This density profile is remarkably shallow compared with outer envelopes of normal SNe (Chevalier & Soker, 1989; Matzner & McKee, 1999). For the energetic SN Ic-BL 1997ef, Mazzali et al. (2000) introduced a flat outer density slope, $\rho\propto r^{-4}$, for better reproducing its spectral evolution, rather than the steeper one, $\rho\propto r^{-8}$, in the hypernova model of Iwamoto et al. (2000). This modification requires the total explosion energy of $1.75\times 10^{52}$ erg for SN 1997ef, well exceeding the typical value of $10^{51}$ erg. Although the slope with the exponent $-4$ is even shallower than that found in our simulation, the ejecta structure with the flat density slope and the sufficiently large kinetic energy is in an overall agreement with the hydrodynamic model explored in this work. More recently, the spectral modelings of several SLSNe-I by Mazzali et al. (2016) have found that a steeper outer density slope, $\rho\propto r^{-7}$, is more appropriate, which is roughly consistent with the original density structure expected in a standard SN explosion. The different outer density structure may indicate the difference in the energy redistribution process in the two energetic populations of SNe Ic-BL and SLSNe-I.

In addition, the clumpy SN ejecta make it easier for high-energy photons and particles from the central energy source to penetrate through the ejecta (Arons, 2003). As a result, the energy deposition from the central compact object to the surrounding gas is realized in a more extended manner than well-structured spherical ejecta. Recently, Dessart (2019) carried out non-LTE spectral synthesis calculations of SNe with an extended energy injection region. He confirmed that the way of the energy injection surely affects light curves. The more extended the energy deposition is, the earlier thermal emission can escape from the ejecta, making them brighter in early epochs.

The ionization states of different layers in the ejecta and the associated emission lines would also be affected. Indeed, in the recent modelings of the nebular spectra of some SLSNe-I by Jerkstrand et al. (2017), they introduced clumpy ejecta in order to explain O I recombination lines and other spectral features. Dessart (2019) also investigated the effect of clumpy ejecta (see, also, Dessart et al., 2018). The comparison of spectral synthesis calculations with some SLSNe-I spectra suggests that SLSNe-I require clumpy ejecta to some extent even at maximum light. The hot bubble breakout and the subsequent efficient mixing of the ejecta give one possible explanation for the clumpy ejecta.

4.1.2 Non-thermal emission

Another important feature of supernova ejecta with the hot bubble breakout is the presence of a mildly relativistic component in the ejecta. The channel flows deposit a relatively large amount of energies into a small fraction of the ejecta at outer layers, realizing a high energy-to-mass ratio at the outermost layer. It is widely recognized that supernova ejecta colliding with an ambient medium create the forward and reverse shocks, producing non-thermal particles (see Chevalier & Fransson, 2017, for a review). The fast ejecta can efficiently dissipate the kinetic energy in the presence of a dense ambient medium and give rise to bright non-thermal emission via synchrotron and inverse Compton processes. Some energetic SNe Ic-BL are indeed known as bright radio sources. In the context of radio-bright SNe, radial density profiles with a flat slope are likely to explain their bright radio emission. For example, Chakraborti et al. (2015) suggested $\rho\propto r^{-6}$ for SN 20012ap. Recently, we have calculated radio and X-ray emission from supernova ejecta with power-law density profiles implied by the hot bubble breakout (Suzuki & Maeda, 2018). We found that ejecta with radial profiles of $\rho\propto r^{-5}$ or $\rho\propto r^{-6}$, could produce bright radio and X-ray emission with luminosities similar to those of radio-bright SNe Ic-BL, SN 1998bw and 2009bb.

On the other hand, in spite of dedicated efforts to detect radio, X-ray, and gamma-ray emission from SLSNe-I, any plausible signal have not been detected up to a few years after the explosion (Levan et al. 2013; Coppejans et al. 2018; Margutti et al. 2018; Renault-Tinacci et al. 2018; Hatsukade et al. 2018, but see the case of SCP 06F6; Gänsicke et al. 2009). This may indicate that SLSNe-I do not accompany a high-velocity ejecta component, although it is a natural outcome of the violent hot bubble breakout. As we suggested in Suzuki & Maeda (2018), in SLSNe-I, most of the injected energy may be used for bright thermal emission rather than accelerating the forward shock, resulting in less violent or even the absence of the hot bubble breakout. However, the sample size is still limited and thus it is probably too early to conclude. Future survey and follow-up observations of SLSNe-I across wide wavelength ranges would figure out if most SNSNe-I lack radio and X-ray signals indicating central engine activities.

The recent detection of a persistent radio source at the position consistent with the SLSN-I PTF10hgi (Eftekhari et al., 2019) might indicate that the SLSN-I was actually powered by a central compact object and the radio emission from the wind nebula driven by the compact source is now unveiled $\sim 9$ years after its discovery. If the spatial coincidence would turn out to be true, this detection along with tight radio upper limits for other SLSNe-I in a few years might imply that radio emission can penetrate the surrounding gas only after it has become transparent to radio waves in $\sim 10$ years, which is in fact consistent with theoretical expectations (Omand et al., 2018; Margalit et al., 2018a). The wind nebula embedded in the SN ejecta is also expected to produce high-energy photons, which ionize the SN ejecta and eventually escape into interstellar space (i.e., ionization breakout; Metzger et al., 2014; Margalit et al., 2018a). Since the clumpy density structure of the SN ejecta would make it easier for such high-energy photons and/or radio waves to escape, such electromagnetic signals may be observed earlier than the theoretical expectations based on spherical SN ejecta.

4.2 Supernova ejecta without hot bubble breakout

In contrast, the density structure of the supernova ejecta with a moderately powerful energy source would be well stratified in the absence of the hot bubble breakout. The hot bubble is well confined by the ram pressure of the ejecta balancing the thermal pressure of the hot gas. Therefore, the internal energy of the hot bubble is eventually lost due to adiabatic expansion and radiative diffusion. It is the diffusion time throughout the ejecta that determines the dominant energy loss process. When the diffusion time scale is much longer than the time scale of the adiabatic expansion, the hot bubble eventually cools by exerting pressure on the surrounding gas and most energy would be shared among the ejecta as the kinetic energy. On the other hand, with the diffusion time scale shorter than the expansion time scale, the internal energy of the hot bubble can diffuse out through layers stratified on the hot bubble rather than being lost by adiabatic cooling. Specifically, the radiation front can overtake the shock front and then emerge from the outermost layer. As we have estimated in Section 3.2.2, the time scale required for the forward shock to reach the outer edge of the ejecta is of the order of $100$ years. Thus, the diffusion time scale is clearly much shorter than this time scale, suggesting that a considerable fraction of the internal energy of the hot bubble would be lost via radiative diffusion.

In this case, the resultant ejecta are composed of a high-density shell confining the hot bubble and well-stratified outer layers of the unshocked ejecta. The outer layers are kept hot due to diffusing photons from the hot bubble. These components would keep spherical symmetry to some extent as long as the energy injection and the original SN ejecta are roughly spherical. In other words, situations expected in 1D spherical calculations (e.g., Kasen & Bildsten, 2010) can apply. Therefore, bright thermal emission in optical wavelengths is expected, depending on the time scales of adiabatic cooling and radiative loss. The photospheric temperature is high and elements in the outer layers of the ejecta are highly ionized because of the power source. However, the expected density slope of the outermost layer is similar to normal stripped-envelope CCSNe, leading to spectral lines with moderate widths. As we have discussed in Section 4.1, SNe with bright thermal emission but without bright radio and X-ray emission could be explained by this scenario.

4.3 Two types of supernovae powered by central engines?

The considerations in Sections 4.1 and 4.2 suggest that there may be two types of supernovae with central power sources, supernovae with relatively powerful and moderate additional energy sources. The former population would show bright optical emission and broad-line spectral features, and potentially be a bright source in radio, X-ray, and even gamma-ray. The latter population could also give rise to bright optical emission powered by a central energy source but without bright non-thermal emission.

These findings can be related to the potential diversity of SLSNe-I. Whether or not SLSNe-I are made up of two or more sub-classes has been paid great attention. Investigating their observational properties, some authors have claimed that there are fast- and slow-evolving populations (e.g., Inserra et al., 2013, 2017). Quimby et al. (2018) classified SLSNe-I in their PTF samples into PTF12dam-like and 2011ke-like objects, where the latter show smoother spectra around their optical maxima than the former. One simple interpretation of the systematic difference in their spectral appearance is that it reflects different ejecta structures. For a shallow radial density profile, a specific range of density corresponds to wider radius and velocity ranges than for a steep radial density profile, making absorption troughs broad. Therefore, the population of SLSNe-I with smooth spectra may indicate relatively flat density profiles.

The different light curve evolution of SLSNe-I can also be a key to unveiling the potential diversity of SLSNe. Fast-evolving SLSNe-I appear to have smaller ejecta mass than the slowly-evolving population (Nicholl et al., 2015). Given that the ejecta mass is positively correlated with the explosion energy for normal stripped-envelop SNe (e.g., Lyman et al., 2016), ejecta with a smaller total mass are more likely to be affected by the central energy injection. In other words, even when the same injection energy is assumed, less massive ejecta (supposedly larger $E_{\mathrm{inj}}/E_{\mathrm{sn}}$) more likely experience the efficient mixing caused by the central energy injection. The 2011ke-like objects analyzed by Quimby et al. (2018) indeed exhibit fast evolving light curves, indicating smaller ejecta masses. Thus, the different spectral behavior may be explained by the different impacts of the central energy injection into the SN ejecta.

Unified scenarios for SLSNe-I and SN Ic-BL (e.g., Metzger et al., 2015; Kashiyama et al., 2016) indicate that they share several common photometric and spectroscopic properties. In fact, there are some observational implications on the connection, such as, spectroscopic properties (Pastorello et al., 2010; Nicholl et al., 2016; Jerkstrand et al., 2017; Liu et al., 2017), and host galaxies (Neill et al., 2011; Lunnan et al., 2014; Leloudas et al., 2015a; Thöne et al., 2015; Angus et al., 2016; Perley et al., 2016; Chen et al., 2017; Schulze et al., 2018). More recently, Blanchard et al. (2019) reported a transitional event, SN 2017dwh, which outshined as brightly in optical bands as SLSNe-I at early epochs, but showed a dramatic transition from an SLSNe-like blue spectrum to a red and broad-lined one similar to SN Ic-BL spectra. They also identified an absorption feature around $3200\mathrm{\AA}$, which they attributed to Co II. This requires highly efficient mixing of iron-peak elements, which is usually embedded in the deep interior of the ejecta, to outer layers. A sufficiently powerful engine can potentially explain the highly bright emission and efficient mixing at the same time, as demonstrated by our simulations.

4.4 Final remarks

Finally, we make some remarks on the simulation results and the scenario introduced above.

The first remark is the resolution dependence of the numerical simulations. Even with the high spatial resolution realized by the AMR technique, the development of the hydrodynamic instability in our simulations show some artificial structure along the numerical grid. For example, as seen in Figure 1, the Rayleigh-Taylor fingers are elongated along the $x$- and $z-$axes before the breakout. This is probably because of the grid structure. Therefore, care must be taken for the angular dependence of the column density or other hydrodynamic variables important for leakage of photons and high energy particles, and viewing angle effects. Nevertheless, the angle-averaged distributions, which are mainly discussed in this study, can be understood in analytical ways and therefore is robust.

We also note that our simulations do not take into account radiative transfer effects. For short diffusion time scales for photons in the ejecta, the energy loss via radiative diffusion must decrease the pressure of the hot bubble, probably making the breakout less violent or even suppressing it. This is expected when the diffusion time scale is comparable to the characteristic time scale of the energy injection, $E_{\mathrm{inj}}/L$. This effect should be investigated by radiation-hydrodynamic simulations appropriately incorporating radiative transfer effects, which we leave as a future study.

5 Conclusions

In this work, we performed 3D hydrodynamic simulations of supernova ejecta with a central energy source and presented their results. We performed simulations with the injected energy much larger than and comparable to the initial kinetic energy of the ejecta to reveal how the amount of the injected energy affects the dynamical evolution of the ejecta. The simulations with two different setups show a remarkable contrast, confirming the expectation that the total amount of the injected energy relative to the original explosion energy of the supernova is one of the most important parameters. The two qualitatively different models lead to different density structures, which can be probed by spectroscopic observations of extraordinary SNe. We point out that the engine-driven SN ejecta with different density structure can explain the diversity of energetic SNe.

We appreciate the anonymous referee for his/her constructive comments, which greatly helped us improve the manuscript. Numerical computations were carried out in part on the Cray XC50 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan. A.S. acknowledges support by Japan Society for the Promotion of Science (JSPS) KAKENHI Grand Number JP19K14770. K.M. acknowledges support by JSPS KAKENHI Grant Numbers JP18H05223, JP18H04585, and JP17H02864.

Appendix A Self-similar expansion of a hot bubble in a power-law ejecta

In this section, we consider the expansion of a spherical hot region in freely expanding ejecta. We especially focus on the self-similar behavior of the forward shock and the pressure of the hot bubble so that the solution can be applied to the numerical result in Section 3.

We assume the following power-law density profile for the freely expanding ejecta corresponding to the outer envelope,

[TABLE]

The spherical hot region is located in the central region of the ejecta and is surrounded by the forward shock front propagating in the ejecta. The forward shock is located at $r=R_{\mathrm{fs}}(t)$ and evolves in a self-similar way,

[TABLE]

where the exponent $\lambda$ is determined later. In order for the shock front to travel outward with respect to the freely expanding medium, the exponent $\lambda$ is required to be larger than unity, $\lambda>1$. The volume $V_{\mathrm{b}}$ of the hot bubble is thus given by,

[TABLE]

The temporal evolution of the pre-shock density is obtained by inserting the forward shock radius into Equation (A1),

[TABLE]

The jump condition for a strong shock leads to the following post-shock pressure $p_{\mathrm{f}}$,

[TABLE]

On the other hand, the internal energy $E_{\mathrm{int}}$ of the hot bubble evolves in an adiabatic way in the absence of the energy injection,

[TABLE]

From the assumed self-similarity, the pressure $p_{\mathrm{b}}$ of the hot bubble is proportional to the post-shock pressure $\propto t^{(2-m)\lambda+m-3}$. In the following, we assume that the pressure in the hot bubble is uniform and identical with the post-shock pressure. In other words, the entire region in the hot bubble is causally connected, which is justified by the numerical result in Section 3. Therefore, the internal energy of the hot bubble is simply given by the product of the volume of the hot bubble and the internal energy density $p_{\mathrm{b}}/(\gamma-1)$,

[TABLE]

By comparing Equations (A5) and (A7), the exponent $\lambda$ should satisfy the following relation,

[TABLE]

or equivalently,

[TABLE]

The requirement $\lambda>1$ leads to $m>3\gamma+2$. In other words, the density slope should be sufficiently steep, $m>6$ for $\gamma=4/3$.

We apply this scaling relation to the late-time evolution of model E51, where the hot bubble is expanding into the spherical supernova ejecta with a steep density profile $m=10$ and $\gamma=4/3$. In this case, the self-similar expansion law derived above can apply and the exponent $\lambda$ is found to be $\lambda=5/4$. Accordingly, the internal energy decays as $E_{\mathrm{int}}\propto t^{3\lambda(1-\gamma)}=t^{-5/4}$.

Another case worth considering is the hot bubble expansion in a shallow density slope, $m<3\gamma+2$, as in the inner ejecta, $-d\ln\rho/d\ln r=1<3\gamma+2$, considered in this study. In this case, the exponent $\lambda$ is below unity and apparently violates the assumption of the shock propagating outward with respect to the freely expanding medium. Namely, the shock front is receding into deeper layers of the ejecta, which is unphysical. In fact, under this condition, the pressure $p_{\mathrm{b}}$ of the hot bubble (equation A7), declines faster than the post-shock pressure $p_{\mathrm{f}}$ required for the outward shock propagation (equation A5). Therefore, the pressure balance will never be reached. As a result, the hot bubble is confined in deeper region and passively expands according to the expansion of the ambient ejecta, $\lambda\sim 1$.

We also note that the assumption of the adiabatic expansion is no longer valid when the radiative diffusion timescale is much shorter than the dynamical timescale of the hot bubble expansion. Significant radiative loss at the shock front effectively reduces the adiabatic exponent of the post-shock gas down to unity, $\gamma\rightarrow 1$. Accrodingly, the exponent $\lambda$ also approaches unity, $\lambda\rightarrow 1$. In other words, the hot bubble can also be stuck in the expanding ejecta when the energy loss is considerable.