Supernova Ejecta Interacting with a Circumstellar Disk. I. Two-dimensional Radiation-hydrodynamic Simulations

The left-hand sides of Equations (1)–(5) represent the transport of mass, momentum, and energy in the physical space. Therefore, these equations with the source terms being zero can be integrated in standard numerical ways. We adopt the same steps as the previous work (Suzuki et al. 2016). For Equations (1) and (2), we employ the so-called M1 closure scheme (Levermore 1984), where the Eddington tensor ${D}^{{ij}}={P}_{{\rm{r}}}^{{ij}}/{E}_{{\rm{r}}}$ is given by the following analytic function of E_r and ${F}_{{\rm{r}}}^{i}:$

$$\begin{eqnarray}&&{D}^{{ij}}=\displaystyle \frac{1-\chi }{2}{\delta }^{{ij}}+\displaystyle \frac{3\chi -1}{2}{n}^{i}{n}^{j},\end{eqnarray} \tag{ 9 }$$

where ${n}^{i}={F}_{{\rm{r}}}^{i}/| {F}_{{\rm{r}}}|$ and the parameter χ is calculated as

$$\begin{eqnarray}&&\chi =\displaystyle \frac{3+4{f}^{i}{f}_{i}}{5+2\sqrt{4-3{f}^{i}{f}_{i}}},\end{eqnarray} \tag{ 10 }$$

with ${f}^{i}={F}_{{\rm{r}}}^{i}/{E}_{{\rm{r}}}$. For the spatial reconstruction of E_r and ${F}_{{\rm{r}}}^{i}$, we employ the third-order weighted, essentially nonoscillatory scheme (Liu et al. 1994; Jiang & Shu 1996).

For the hydrodynamics equations, we also use the commonly adopted Harten–Lax–van Leer–Einfeldt Riemann solver (Mignone & Bodo 2005) combined with the third-order MUSCL reconstruction.

The equations for radiation field and hydrodynamic variables are related by the coupling terms G⁰ and Gⁱ, whose functional forms are as follows:

$$\begin{eqnarray}&&{G}^{0}={\rm{\Gamma }}\bar{\rho }\bar{{\kappa }_{{\rm{a}}}}({a}_{{\rm{r}}}{\bar{T}}_{{\rm{g}}}^{4}-{\bar{E}}_{{\rm{r}}})-{\rm{\Gamma }}\bar{\rho }({\bar{\kappa }}_{{\rm{a}}}+{\bar{\kappa }}_{{\rm{s}}}){\beta }_{i}{\bar{F}}_{{\rm{r}}}^{i}\end{eqnarray} \tag{ 11 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{G}^{i} & = & -\bar{\rho }({\bar{\kappa }}_{{\rm{a}}}+{\bar{\kappa }}_{{\rm{s}}}){\bar{F}}_{{\rm{r}}}^{i}+{\rm{\Gamma }}\bar{\rho }{\bar{\kappa }}_{{\rm{a}}}\left({a}_{{\rm{r}}}{\bar{T}}_{{\rm{g}}}^{4}-{\bar{E}}_{{\rm{r}}}\right){\beta }^{i}\\ & & -\,\displaystyle \frac{{{\rm{\Gamma }}}^{2}}{{\rm{\Gamma }}+1}\bar{\rho }({\bar{\kappa }}_{{\rm{a}}}+{\bar{\kappa }}_{{\rm{s}}}){\beta }^{i}{\beta }_{j}{\bar{F}}_{{\rm{r}}}^{j},\end{array}\end{eqnarray} \tag{ 12 }$$

where ${\bar{\kappa }}_{{\rm{a}}}$ and ${\bar{\kappa }}_{{\rm{s}}}$ are the absorption and scattering opacities, respectively. We have assumed that the scattering process is isotropic in the comoving frame. The radiation energy density and the flux appearing in these expressions are defined in the comoving frame and related to those in the laboratory frame through the Lorentz transformation,

$$\begin{eqnarray}&&{E}_{{\rm{r}}}={{\rm{\Gamma }}}^{2}\left({\bar{E}}_{{\rm{r}}}+2{\beta }_{i}{\bar{F}}_{{\rm{r}}}^{i}+{\beta }_{i}{\beta }_{j}{\bar{P}}_{{\rm{r}}}^{{jk}}\right)\end{eqnarray} \tag{ 13 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{F}_{{\rm{r}}}^{i} & = & {\rm{\Gamma }}{\bar{F}}_{{\rm{r}}}^{i}+{\rm{\Gamma }}{\beta }_{j}{\bar{P}}_{{\rm{r}}}^{{ij}}\\ & & +\,{{\rm{\Gamma }}}^{2}{\beta }^{i}\left({\bar{E}}_{{\rm{r}}}+\displaystyle \frac{2{\rm{\Gamma }}+1}{{\rm{\Gamma }}+1}{\beta }_{j}{\bar{F}}_{{\rm{r}}}^{j}+\displaystyle \frac{{\rm{\Gamma }}}{{\rm{\Gamma }}+1}{\beta }_{j}{\beta }_{k}{\bar{P}}_{{\rm{r}}}^{{jk}}\right),\end{array}\end{eqnarray} \tag{ 14 }$$

(see, e.g., Mihalas & Mihalas 1984), which close the equations. The numerical methods to solve these equations are presented in Appendix A.1.

Our treatment of electron scattering has some limitations. While the free–free process immediately realizes gas–radiation equilibrium in sufficiently dense material, it is not always effective. For shocks propagating in dilute gas, for example, post-shock gas density could be so small that the free–free process is not effective for creating enough photons for the equilibrium. Alternatively, Compton scattering can be a dominant energy transfer process for gas and radiation. With an insufficient number of photons, the photon spectrum becomes a Wien function rather than a Planck function. This circumstance is expected for SN shock breakout in a dilute stellar wind (Weaver 1976). In this way, our two-temperature treatment sometimes does not capture the gas–radiation coupling correctly. In this study, however, we focus on shock breakout from a dense CSM, where free–free processes produce enough photons to achieve bright thermal emission.

We consider free–free absorption and electron scattering as the dominant radiative processes in this setting. We assume the commonly adopted opacity formulae,

$$\begin{eqnarray}&&{\bar{\kappa }}_{{\rm{a}}}=3.7\times {10}^{22}(1+{X}_{{\rm{h}}})({X}_{{\rm{h}}}+{X}_{\mathrm{he}})\rho {\bar{T}}_{{\rm{g}}}^{-7/2}\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1},\end{eqnarray} \tag{ 15 }$$

for free–free absorption and

$$\begin{eqnarray}&&{\bar{\kappa }}_{{\rm{s}}}=0.2(1+{X}_{{\rm{h}}})\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}\end{eqnarray} \tag{ 16 }$$

for electron scattering (in cgs units; see, e.g., Rybicki & Lightman 1979). Here we assume a fully ionized gas. We note that this is not always a valid assumption. As we shall see below, the photospheric temperature decreases down to the recombination temperature of hydrogen, <6000 K, at later epochs. Then, the assumption of fully ionized gas is no longer justified.

We initially assume freely expanding spherical SN ejecta. We start our simulation at t₀ = 1000 s. Thus, the radial velocity of a layer at a radius R is initially given by

$$\begin{eqnarray}&&v=\displaystyle \frac{R}{{t}_{0}}\end{eqnarray} \tag{ 17 }$$

for $v\lt {v}_{\max }$ and v = 0 otherwise. Here we denote the 3D radius by $R\equiv {({r}^{2}+{z}^{2})}^{1/2}$ in order to distinguish it from the radial coordinate r. We also introduce the inclination angle $\theta ={\cos }^{-1}(r/R)$ to specify an angle measured from the symmetry axis z = 0. The initial density structure is assumed to be the commonly adopted broken power-law distribution with the break velocity v_br (Chevalier & Soker 1989; Matzner & McKee 1999),

$$\begin{eqnarray}&&{\rho }_{\mathrm{ej}}(r,z)=\displaystyle \frac{{f}_{3}{M}_{\mathrm{ej}}}{4\pi {v}_{\mathrm{br}}^{3}{t}_{0}^{3}}g(R/{t}_{0}),\end{eqnarray} \tag{ 18 }$$

for R < v_maxt₀, where M_ej is the total mass of the ejecta and

$$\begin{eqnarray}&&{f}_{l}=\displaystyle \frac{\left(m-l\right)\left(l-\delta \right)}{m-\delta -\left(l-\delta \right){\left({v}_{\mathrm{br}}/{v}_{\max }\right)}^{m-l}}.\end{eqnarray} \tag{ 19 }$$

For a sufficiently steep outer density slope and a large ${v}_{\max }/{v}_{\mathrm{br}}$, this expression is approximated as follows:

$$\begin{eqnarray}&&{f}_{l}\simeq \displaystyle \frac{(m-l)(l-\delta )}{m-\delta }.\end{eqnarray} \tag{ 20 }$$

The nondimensional function g(v) is given by

$$\begin{eqnarray}g(v)=\left\{\begin{array}{ccc}{\left(\tfrac{v}{{v}_{\mathrm{br}}}\right)}^{-\delta } & \mathrm{for} & v\leqslant {v}_{\mathrm{br}},\\ {\left(\tfrac{v}{{v}_{\mathrm{br}}}\right)}^{-m} & \mathrm{for} & {v}_{\mathrm{br}}\lt v.\end{array}\right.\end{eqnarray} \tag{ 21 }$$

The exponents δ and m are fixed to be δ = 1 and m = 10 throughout this study. The break velocity v_br is determined by specifying the total mass and kinetic energy, M_ej and E_sn, of the ejecta,

$$\begin{eqnarray}&&{v}_{\mathrm{br}}={\left(\displaystyle \frac{2{f}_{5}{E}_{\mathrm{sn}}}{{f}_{3}{M}_{\mathrm{ej}}}\right)}^{1/2}\simeq {\left[\displaystyle \frac{2(m-5)(5-\delta ){E}_{\mathrm{sn}}}{(m-3)(3-\delta ){M}_{\mathrm{ej}}}\right]}^{1/2}.\end{eqnarray} \tag{ 22 }$$

We assume that the initial gas internal energy density of the ejecta is proportional to the initial kinetic energy distribution,

$$\begin{eqnarray}&&{\bar{E}}_{{\rm{g}}}(r,z)=\epsilon \displaystyle \frac{{\rho }_{\mathrm{ej}}(r,z)}{2}{\left(\displaystyle \frac{R}{{t}_{0}}\right)}^{2},\end{eqnarray} \tag{ 23 }$$

with  = 0.05. The internal-to-kinetic energy ratio is expected to be around unity at the very beginning of the expansion of the SN ejecta. Then, the internal energy rapidly decreases in the absence of any heating source and immediately becomes negligible compared to the kinetic counterpart. The expanding outer SN envelope with a power-law radial density profile is obtained as an asymptotic behavior after the internal energy becomes negligible to the kinetic energy (Chevalier & Soker 1989). Therefore, from the beginning, we set the internal energy of the ejecta to only 5% of the kinetic one so that the radial density structure mimics the cold ejecta realized well after the explosion. This treatment is justified when the initial thermal energy loaded on the SN ejecta is less likely to contribute to the bolometric luminosity as in interacting SNe. The radiation energy density and radiative flux are initially zero. Although there is no radiation field in the ejecta at the beginning of the simulation, the thermal equilibrium between gas and radiation, where the gas and radiation temperatures are identical, is immediately achieved in the SN ejecta.

In this study, we fix the total mass and initial kinetic energy of the ejecta to be M_ej = 10 M_⊙ and E_sn = 10⁵¹ erg. Thus, the corresponding break velocity is v_br = 3.8 × 10³ km s⁻¹. The outer ejecta initially extend out to v_maxt₀ at time t₀, where the outermost layer is adjacent to the inner edge of the CSM. In the following, we set v_max = 1. 26 × 10⁹ cm s⁻¹ or v_maxt₀ = 1.26 × 10¹² cm, which corresponds to v_br/v_max = 0.3.

We consider a centrally concentrated CSM embedded in a steady wind. Such CSMs are thought to be an outcome of still unclear mass-loss activities of massive stars 10–1000 yr before their core collapse, and therefore their density and velocity structures are highly uncertain. We assume that both the CSM and wind densities follow an inverse square law. The dense CSM is truncated at a radius of $5\times {10}^{15}$ cm, which corresponds to an enhanced mass loss 10–100 yr prior to the core collapse for a wind velocity of 10–100 km s⁻¹ (carbon- and oxygen-burning stages). Dense CSMs confined within a few 10¹⁵ cm are often adopted in light-curve modelings of Type IIn SNe (e.g., Ginzburg & Balberg 2012; Moriya et al. 2013a).

Thus, for spherical CSMs, the density profile is described by

$$\begin{eqnarray}&&\rho (r,z)={\rho }_{\mathrm{wind}}(R)+{\rho }_{\mathrm{sph}}(R),\end{eqnarray} \tag{ 24 }$$

where the wind and CSM components are given by

$$\begin{eqnarray}&&{\rho }_{\mathrm{wind}}(R)={A}_{\star }{R}^{-2}\end{eqnarray} \tag{ 25 }$$

and

$$\begin{eqnarray}&&{\rho }_{\mathrm{sph}}(R)=\displaystyle \frac{{A}_{\mathrm{sph}}}{{R}^{2}}\exp \left[-{\left(\displaystyle \frac{R}{{R}_{\mathrm{csm}}}\right)}^{p}\right],\end{eqnarray} \tag{ 26 }$$

with p = 10. In the following, the wind density parameter is set to A_⋆ = 5 × 10¹¹ g cm⁻¹, which corresponds to a steady mass loss at a rate of 10⁻⁶ M_⊙ yr⁻¹ for a constant wind velocity of 100 km s⁻¹. This value leads to a sufficiently dilute wind component so that it does not affect the expansion of the SN ejecta. The exponential factor in this spherical CSM density profile realizes a smooth cutoff around R = R_csm, beyond which the wind component dominates. The CSM mass M_csm is expressed as a function of the normalization constant A_sph and the outer radius R_csm,

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{csm}} & = & 4\pi {\int }_{0}^{\infty }{\rho }_{\mathrm{sph}}(R){R}^{2}{dR}\\ & = & 4\pi {A}_{\mathrm{sph}}{R}_{\mathrm{csm}}{\rm{\Gamma }}(1+1/p),\end{array}\end{eqnarray} \tag{ 27 }$$

where Γ(x) is a gamma function and Γ(1 + 1/p) ≃ 0.951 for p = 10.

For a CSM disk, we assume the following functional form:

$$\begin{eqnarray}&&\rho (r,z)={\rho }_{\mathrm{wind}}(R)+{\rho }_{\mathrm{disk}}(r,z),\end{eqnarray} \tag{ 28 }$$

with

$$\begin{eqnarray}&&{\rho }_{\mathrm{disk}}(r,z)=\displaystyle \frac{{A}_{\mathrm{disk}}}{{R}^{2}}\exp \left[-{\left(\displaystyle \frac{R}{{R}_{\mathrm{csm}}}\right)}^{p}-{\left(\displaystyle \frac{\vartheta }{{\vartheta }_{\mathrm{disk}}}\right)}^{q}\right],\end{eqnarray} \tag{ 29 }$$

where ϑ is an angle measured from the equator,

$$\begin{eqnarray}&&\vartheta ={\cos }^{-1}\left(\displaystyle \frac{r}{R}\right).\end{eqnarray} \tag{ 30 }$$

Throughout this work, we assume q = 4. Therefore, the density is enhanced around the region with a half-opening angle of ϑ_disk. The parameter q determines the slope of the upper and lower edges of the CSM disk. Although we adopt a very simplified CSM disk model, there would be room for improvement for future studies. For example, by assuming a steady disk, the scale height of the upper and lower edges of the disk is determined by the disk temperature. The geometry and structure of CSM disks would offer us a hint toward how they are produced and thus should be studied in more detail. The CSM mass M_csm is expressed as a function of R_csm and ϑ_disk,

$$\begin{eqnarray}&&{M}_{\mathrm{csm}}=4\pi {A}_{\mathrm{disk}}{\rm{\Gamma }}(1+1/p){F}_{q}({\vartheta }_{\mathrm{disk}}){R}_{\mathrm{csm}},\end{eqnarray} \tag{ 31 }$$

where F_q(ϑ_disk) is the covering fraction of the disk with respect to the solid angle,

$$\begin{eqnarray}&&{F}_{q}({\vartheta }_{\mathrm{disk}})={\int }_{0}^{\pi /2}\exp \left[-{\left(\displaystyle \frac{\vartheta }{{\vartheta }_{\mathrm{disk}}}\right)}^{q}\right]\cos \vartheta d\vartheta .\end{eqnarray} \tag{ 32 }$$

The covering fraction yields ${F}_{4}({\vartheta }_{\mathrm{disk}})=0.157$ and 0.310 for ϑ_disk = 10° and 20°. In Figure 1, we show the spatial distribution of the density, the gas and radiation energy densities, and the radial velocity shortly after t = t₀ for model D10_M10 (see Table 1 for the model parameters).

Zoom In Zoom Out Reset image size

In Figures 2 and 3, we plot the radial profiles of several hydrodynamic variables obtained by the simulations with M_csm = 0.1 and 10 M_⊙, respectively. The simulations qualitatively reproduce early studies on SN ejecta and the associated blast wave evolution in a dense CSM (e.g., Ginzburg & Balberg 2012; Moriya et al. 2013a). The shock front driven by the expanding ejecta sweeps the surrounding medium. Initially, the shock is almost adiabatic because of the short mean free path of photons in the dense medium. In this stage, the time required for the energy exchange between gas and radiation is much shorter than the dynamical timescale, resulting in the almost identical gas and radiation temperature distributions. At later epochs, however, the pre-shock gas density and corresponding optical depth gradually decrease. This allows the post-shock radiation to diffuse into the pre-shock region to create the so-called "precursor," where the pre-shock gas is heated and accelerated by radiation. The precursor is followed by a layer with a very high gas temperature (known as a Zel'dovich spike; see, e.g., Zel'dovich & Raizer 1967). This layer is immediately behind the shock front, at which the shock kinetic energy is directly converted into the internal energy of the post-shock gas. The internal energy of the post-shock gas is gradually shared by gas and radiation around the shock front, achieving gas–radiation equilibrium behind the shock front. When the precursor breaks out from the photosphere in the CSM, photons produced in the post-shock region can escape into the surrounding interstellar space. The hot layer behind the shock front has become wide because of the inefficient coupling between gas and radiation in the wind region, R > R_csm.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 4 shows the structure of the forward shock front. We plot the density and temperature profiles for the models with M_csm = 0.1 M_⊙ at t = 10⁵ s and M_csm = 10 M_⊙ at t = 10⁷ s. The high-temperature spike followed by a cooled, piled-up gas layer is clearly recognized. The width of the high-temperature layer is roughly estimated as follows. The shock velocity at later epochs shown in Figures 6 and 7 is typically v_shock ∼ 0.01c. Therefore, the post-shock gas temperature can be as high as

$$\begin{eqnarray}&&{T}_{{\rm{g}},\mathrm{shock}}\simeq \displaystyle \frac{\mu {m}_{{\rm{u}}}{v}_{\mathrm{shock}}^{2}}{2{k}_{{\rm{B}}}}\simeq 3\times {10}^{9}\,{\rm{K}}.\end{eqnarray} \tag{ 37 }$$

Neglecting the radiation energy and relativistic effects, the cooling rate of the gas energy density is given by

$$\begin{eqnarray}&&\displaystyle \frac{d{\bar{E}}_{{\rm{g}}}}{{dt}}\simeq -c\bar{\rho }{\bar{\kappa }}_{{\rm{a}}}{a}_{{\rm{r}}}{\bar{T}}_{{\rm{g}}}^{4},\end{eqnarray} \tag{ 38 }$$

which leads to the following cooling timescale:

$$\begin{eqnarray}&&{t}_{\mathrm{cool}}=\displaystyle \frac{{\bar{E}}_{{\rm{g}}}}{d{\bar{E}}_{{\rm{g}}}/{dt}}=\displaystyle \frac{{k}_{{\rm{B}}}}{(\gamma -1)\mu {m}_{{\rm{u}}}{\bar{\kappa }}_{{\rm{a}}}{a}_{{\rm{r}}}{\bar{T}}_{{\rm{g}}}^{3}}.\end{eqnarray} \tag{ 39 }$$

Accordingly, the width l_cool of the relaxation layer is roughly estimated to be

$$\begin{eqnarray}\begin{array}{rcl}{l}_{\mathrm{cool}} & \simeq & {v}_{\mathrm{shock}}{t}_{\mathrm{cool}}\\ & \simeq & 8\times {10}^{12}{\left(\displaystyle \frac{\rho }{{10}^{-11}{\rm{g}}{\mathrm{cm}}^{-3}}\right)}^{-1}{\left(\displaystyle \frac{{v}_{\mathrm{shock}}}{0.01c}\right)}^{2}\ \mathrm{cm}.\end{array}\end{eqnarray} \tag{ 40 }$$

Taking ρ = 10⁻¹¹ g cm⁻³ and v_shock = 0.01c, for example, this order-of-magnitude estimation gives l_cool ∼ 10¹³ cm, which explains the width of the high-temperature spike in the upper panel of Figure 4. The spike in the lower panel is wider than that in the upper panel, probably reflecting its lower post-shock density. In reality, the density of the gas flow following the spike widely varies due to the piling up of material, which makes it difficult to correctly estimate the width of the high-temperature layer. However, the above order-of-magnitude estimation seems to work well. As seen in Figure 4, the high-temperature spikes are covered by less than 10 numerical cells. Although we avoid a spike represented by only a single numerical cell, the shock structure is likely smeared out to some extent.

Zoom In Zoom Out Reset image size

For the model with M_csm = 0.1 M_⊙, the CSM mass is much smaller than the total ejecta mass of M_ej = 10 M_⊙. Therefore, only a small fraction of the outer ejecta is swept by the reverse shock front. The maximum velocity of the ejecta before the shock emergence from the outer edge of the CSM (t = 10⁶ s) is ≃6 × 10³ km s⁻¹. Then, after the emergence, the shock front accelerates again due to the steep density gradient around the outer edge of the CSM.

For the model with ${M}_{\mathrm{csm}}=10{M}_{\odot }$, on the other hand, the total mass of the CSM is comparable to the ejecta. Thus, the ejecta significantly decelerate and dissipate a considerable fraction of the kinetic energy. As seen in Figure 3, the maximum velocity of the ejecta decreases down to 2 × 10³ km s⁻¹ at t = 10⁷ s. The forward shock front has not reached the outer edge of the CSM even at t = 10⁷ s and thus still contributes to the bright thermal emission. In this model, a dense and geometrically thin shell is formed. The density of the shell is higher than that of the surrounding unshocked gas by about 2 orders of magnitude, which cannot be achieved by the adiabatic shock jump condition alone. This is the so-called cooling shell, in which the effective energy loss makes material piling up in the layer between the forward and reverse shock fronts.

The bolometric light curves of the 1D spherical models are presented in Figure 5. In the upper panel, we plot the cumulative radiated energy (Equation (36)). In the lower panel, the light curves are compared with each other and the radioactive energy deposition rate ${\dot{E}}_{\mathrm{Ni}+\mathrm{Co}}$, with a nickel mass of M_Ni = 0.1 M_⊙ (Nadyozhin 1994). We note that some spikes in the light curves are artificially produced due to numerical treatments, especially the nonuniform AMR grid structure (this is alleviated in 2D simulations; see below). Nevertheless, the general trends of the light curves around the peak luminosity are nicely captured.

Zoom In Zoom Out Reset image size

The CSM mass governs the peak and characteristic timescale of the light-curve evolution. For the lower CSM mass model, M_csm = 0.1 M_⊙, the contribution of the CSM interaction is very limited. It certainly contributes to the light curve within ∼10 days around the peak, at which the forward shock is still in the CSM. At later epochs, the luminosity continuously decays down to values below the radioactive energy deposition rate. The contribution of the CSM interaction becomes more significant for increasing CSM mass. The models with larger CSM masses are more capable of dissipating the kinetic energy of the SN ejecta, but it takes longer. The peak luminosity and total radiated energy plotted in Figure 5 indeed increase for increasing CSM mass. At the same time, the elevated CSM densities make the rise and decay timescales longer. At the beginning, the CSM interaction occurs at a deep layer of the CSM, which is initially hidden from the observer. The higher CSM density and thus the larger optical depth lead to a longer diffusion time. Therefore, it takes a longer time for the photons produced in the interaction layer to emerge from the outer edge of the CSM. As a result, for the model with the highest CSM mass, bright emission with ${L}_{\mathrm{bol},\mathrm{iso}}\gtrsim {10}^{43}$ erg s⁻¹ lasts for more than 100 days. A more quantitative discussion of the rise and decline timescales is found in Section 3.4. The total radiated energy reaches terminal values of ${E}_{\mathrm{rad},\mathrm{iso}}=1.3\times {10}^{49}$, 1.2 × 10⁵⁰, and 4.4 × 10⁵⁰ erg for the models with M_csm = 0.1, 1.0, and 10 M_⊙, corresponding to conversion efficiencies of 1.3%, 12%, and 44%.

While the outer part of the SN ejecta (ρ ∝ r^−m) is interacting with a CSM with a power-law radial density profile with an index −s, ρ ∝ R^−s, the kinetic energy dissipation rate is proportional to ${t}^{(6s-15+2m-{ms})/(m-s)}$ (Moriya et al. 2013b). For the adopted parameter set, the exponent is found to be −3/8. When the dissipated kinetic energy is assumed to be converted to radiation at a constant rate, the light curve of the interaction-powered emission roughly follows this time dependence at timescales longer than the photon diffusion time in the CSM. Therefore, the bolometric luminosity behaves as L_bol ∝ t^{−3/ 8}. In Figure 5, we compare this power-law time dependence with the 1D spherical model with M_csm = 1.0 M_⊙. After the bump around the maximum light, at which the timescale of the light-curve evolution is determined by the photon diffusion timescale, the light curve flattens and decays at a similar rate to the power-law function. The numerical light curve is declining more rapidly than the power-law time dependence. This may be caused by a time-dependent conversion efficiency. Recently, Tsuna et al. (2019) developed a light-curve model for interacting SNe by including the effect of a time-dependent conversion efficiency. They claim that the time-dependent conversion efficiency can make light curves decline more rapidly. They also found that the forward and reverse shock components behave differently due to the different post-shock densities, resulting in a double power-law light curve. In the earlier part, where the forward shock dominates, the luminosity decreases more rapidly than t^−3/8 due to the time-dependent conversion efficiency.

Figures 6 and 7 show the temporal evolution of the SN ejecta interacting with the CSM with M_csm = 0.1 and 10 M_⊙ (models S_M01 and S_M10). Their dynamical evolution is similar to that of the 1D models, except for the ejecta–CSM interface. In 2D simulations, the contact surface is subject to hydrodynamic instabilities. As we have noted in Section 3.1, for M_csm = 0.1 M_⊙, the radiative precursor develops even at several 10⁴ s after the explosion. This can also be seen in the spatial distributions of the radiation energy density (lower panels of Figure 6). The radiation front propagates well ahead of the shock front that is still deeply embedded in the CSM, resulting in the pre-shocked region with a high radiation energy density E_r > 1 erg cm⁻³. On the other hand, for the model with the highest CSM mass (lower panels of Figure 7), the region with a high radiation energy density is well confined in the massive CSM at early epochs (up to several 10⁶ s). The radiation front then breaks out from the outer edge of the CSM and fills the surrounding space with radiation. We note that the shock front is still in the CSM even at t = 10⁷ s.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One of the important hydrodynamic features in these 2D simulations is the Rayleigh–Taylor instability. The forward shock front accumulates increasing amounts of the surrounding medium and thus decelerates as it propagates in the CSM. Because of the deceleration, the shocked SN ejecta feel the inertial force toward the direction of the shock propagation, which tries to replace dense gas in the shocked ejecta with relatively dilute gas in the shocked CSM. As a result, Rayleigh–Taylor fingers appear at the ejecta–CSM interface and travel into the shocked CSM. The shortest unstable wavelength of the Rayleigh–Taylor instability is determined by several physical conditions, including diffusion length (e.g., Chevalier & Klein 1978). In these simulations, however, the instability develops from numerical errors mainly caused by AMR grid structures, and the physically determined shortest wavelength is not resolved even with AMR, especially at early epochs.

Another instability potentially playing a role in this setting is the Vishniac instability (also known as the nonlinear thin shell instability). Vishniac (1983) and Ryu & Vishniac (1987) analytically considered the stability of a blast wave traveling in a uniform medium and found that it is overstable when the effective adiabatic index is close to unity (γ_eff < 1.2). In the Vishniac instability, the wave front or shell oscillates with an amplitude increasing with time in a power-law fashion. In a radiative shock, a part of the dissipated shock kinetic energy is lost as escaping radiation, effectively reducing the adiabatic index (Bertschinger 1986; Draine & McKee 1993). Later numerical studies confirm that this kind of instability certainly occurs in radiative shocks in various astrophysical phenomena, including SN remnants in the radiative phase (e.g., Blondin et al. 1998; Michaut et al. 2012; Badjin et al. 2016; Minière et al. 2018).

These two instabilities both make the shock front deviate from its initial spherical shape. While the Rayleigh–Taylor instability is a purely hydrodynamic instability that can develop in any stage of the dynamical evolution, the Vishniac instability plays a role only when the radiative loss in the post-shock region is significant and the effective adiabatic index of the post-shock gas is below a threshold value. When the optical depth from the outer edge of the CSM to the shock front is sufficiently large, the total pressure is dominated by the radiative one, and thus the adiabatic index of the post-shock gas becomes close to that of photon gas, i.e., γ = 4/3. As the forward shock approaches the outer edge of the CSM, the radiative loss becomes increasingly significant. Then, the mixture of the gas and radiation behaves as a gas with an adiabatic index close to unity. This is when the Vishniac instability possibly develops. For lower CSM densities, however, insignificant radiative loss prevents the Vishniac instability from developing even at later epochs. In our simulations, the free–free process is the only radiative process for direct gas cooling. Therefore, for lower CSM densities with insignificant radiative cooling, gas should behave as an ideal gas with an adiabatic index of γ = 5/3, as assumed, and the Vishniac instability is not effective within the timescale covered by the simulations.

In Figure 8, we compare the shock structure at t = 5 × 10⁶ s in the three 2D spherical CSM models. In the spatial distribution shown in the top panel of Figure 8, the ejecta–CSM interface indeed exhibits filamentary structure, which is likely created as a result of the Rayleigh–Taylor instability. In contrast to the models with higher CSM densities, the forward shock front and contact discontinuity are well separated from each other. This suggests that the post-shock gas almost behaves as an ideal gas without significant radiative loss. Under these circumstances, the motion of gas around the ejecta–CSM interface can be turbulent, leading to additional dissipation of the kinetic energy due to collisions of filaments. In most parts, however, the filaments and perturbed shell are confined in the narrow shocked layer between the forward and reverse shock fronts. Even at several 10⁷ s after the explosion, the development of the hydrodynamic instability is not strong enough to modify the global shape of the ejecta. We have also checked that the radial profiles of the hydrodynamic variables are similar to those in 1D spherical models. Thus, the effect of the additional energy dissipation would be quite limited as far as these particular simulations are concerned. As we shall see below, the light curves of 2D spherical models are not significantly different from the corresponding 1D spherical models. However, one caveat is that our simulations assume no initial density and velocity perturbations, and the instabilities develop from small numerical errors. We cannot exclude a possibility that large density and velocity perturbations are present in pre-explosion stars and/or CSMs and enhance the development of hydrodynamics instabilities after the explosion. In such cases, its influence on the energy dissipation and light curves is expected to be more significant.

Zoom In Zoom Out Reset image size

In the bottom panel of Figure 8, the shocked gas in model S_M10 seems to suffer from significant radiative loss. Unlike the models with the smaller CSM masses, the forward shock front and contact discontinuity are not well separated. The density distribution shows a geometrically thin shell, as in the 1D spherical models, although its shape is far from spherical. The creation of a cooling shell suggests that radiative cooling is very effective in this particular model. The nonspherical cooling shell is probably due to the development of the Vishniac instability.

In Figure 9, we plot isotropic equivalent bolometric light curves with different viewing angles (Θ_obs = 0°, 30°, 45°, 60°, and 90°) and compare them with those of 1D spherical simulations. For most parts, the light curves with different viewing angles are similar to each other and the 1D spherical counterparts at early epochs, although small differences are recognized at later epochs. At around the maximum light, photons produced in the ejecta–CSM interface predominantly contribute to the emission. Such photons would have experienced multiple absorption and scattering episodes, which make the radiation field nearly isotropic. This explains why the light curves are similar to each other at early epochs. In contrast, at later epochs, photons from the nonspherical interaction layer escape into the surrounding space without significant absorption and scattering. Then, light curves with different viewing angles start deviating from each other. Nevertheless, the difference between the light curves shown in each panel of Figure 9 is not significant.

Zoom In Zoom Out Reset image size

First, we note that the light curves with Θ_obs = 0° are exceptionally different from those with different viewing angles. A caveat is that this feature is probably due to an artificial effect associated with the assumed axisymmetry. We impose a reflecting boundary condition at the symmetry axis. Therefore, incoming gas motion toward the reflecting boundary encounters the opposite gas motion with the same amplitude, enhancing existing perturbations and numerical errors. We indeed recognize artificial structure around the symmetry axis in the density distributions shown in Figures 6 and 7. Thus, this deviation would probably be resolved in 3D simulations.

Next, we mention the systematic difference in light curves with other viewing angles (e.g., the light curves viewed from 30° are always above those from 45°). In 2D cylindrical geometry, the development of the Raleigh–Taylor instability along a particular radial direction systematically depends on the inclination angle, which causes the differences in the light curves, because the perturbations that developed in the simulation are not actually clumps but rings. This makes the development of perturbations in different orientations systematically different. Since the differences in late-time light curves are probably due to different ways of energy dissipation at the interaction layer at each direction, the late-time light curves are certainly influenced by this systematic effect. These late-time differences would behave correctly in 3D simulations. However, the deviations of the late-time light curves from each other are within a factor of a few, except for Θ_obs = 0°. The peak bolometric luminosities of 2D spherical models that are larger than the corresponding 1D spherical models are at most 37%, 22%, and 7% for models with M_csm = 0.1, 1.0, and 10 M_⊙, respectively. This indicates that the effect of the nonspherical interaction layer due to the Rayleigh–Taylor instability is limited. At least, the instability does not change the luminosity by an order of magnitude. As we have noted in Section 3.2.1, however, preexisting large density and velocity perturbations may enhance the development of the instability to produce global nonspherical structure. In such cases, the luminosity could vary more widely because of the additional energy dissipation by ejecta–clump interaction.

In Figure 9, we also plot light curves of some Type IIn SNe, 1998S, 2010jl, and 2006gy, in the literature. Although it is not our purpose here to precisely reproduce the light curves of particular interacting SNe by adjusting some model parameters, this comparison clearly demonstrates that there are some Type IIn SNe with light-curve shapes similar to the simulation results.

In summary, deviations from 1D spherical models are not significant. Therefore, we can conclude that the effects of multidimensional hydrodynamic instabilities on light curves are negligible, at least for spherical SN ejecta interacting with a spherical CSM.

Figures 10 and 11 present the dynamical evolution of the ejecta–CSM interaction in models D10_M1 and D20_M1. The ejecta start interacting with the surrounding media immediately after the beginning of the simulations. In the presence of a CSM disk, a part of the ejecta traveling around the equator decelerates efficiently, resulting in highly aspherical ejecta structure (top panels of Figures 10 and 11). The polar part of the ejecta is not covered by the CSM; thus, the ejecta can expand almost freely. As a result, the ejecta are squeezed in the presence of the CSM disk, and a ringlike region with no radially expanding ejecta appears behind the CSM disk (hereafter referred to as the void region).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One of the important consequences of this ejecta–disk interaction is that the kinetic energy dissipation is most efficiently happening in the equator, where the ejecta are adjacent to the CSM disk. This is clearly seen in the spatial distributions of the radiation energy density, particularly in the bottom panels of Figures 10 and 11. The region with the highest radiation energy density is found at an off-center location on the equatorial plane.

Another important issue for the emission from SN ejecta is the impact of the viewing angle. As we will see below, the light curves of an SN interacting with a CSM disk are significantly affected by whether they are seen directly or through the CSM. Therefore, what fraction of the ejecta is affected by the CSM disk plays a critical role in determining the expected emission powered by the ejecta–disk interaction. From the density distributions in Figures 10 and 11, there is a general trend that SN ejecta colliding with CSM disks with larger opening angles ϑ_disk suffer more significantly from the squeezing effect. As a result, the solid angle corresponding to an almost freely expanding part of the SN ejecta appears to be a decreasing function of the disk half-opening angle. It is expected that the corresponding solid angle is roughly given by 4π[1 − F_q(ϑ_disk)], where F_q(ϑ_disk) is the covering fraction of the CSM disk (Equation (32)).

To be more precise, however, the boundary between the SN ejecta and the void region is determined by a more complicated hydrodynamics process. As seen in the top panels, the density distributions, of Figure 11, the inner part of the CSM disk is swept by the forward shock immediately after it starts interacting with the SN ejecta. As a result, the inner disk is compressed, and the associated pressure increase should also affect the ejecta that should have avoided the effect of the CSM disk in the ballistic collision case. In this way, a part of the kinetic energy of the ejecta dissipated around the equator is redistributed into ejecta components close to the void region. In the density distributions shown in Figures 10 and 11, the ejecta components traveling along the boundary between the ejecta and the void region precede those along the symmetry axis. This can be understood as a consequence of the complex hydrodynamic interaction between SN ejecta and a CSM disk, which cannot be achieved by considering a ballistic collision alone. The opening angle of the void region indeed seems to be slightly larger than the assumed half-opening angle of the CSM disk.

In addition, the radially traveling ejecta adjacent to the disk surface could suffer from Kelvin–Helmholtz instability, leading to further dissipation of the kinetic energy. Figure 12 shows the density and velocity structures of the interface between the ejecta and the disk for model D20_M1. The radial and angular components of the velocity at coordinates (r, z) are defined as

$$\begin{eqnarray}&&{v}_{R}={v}_{r}\displaystyle \frac{r}{R}+{v}_{z}\displaystyle \frac{z}{R}\end{eqnarray} \tag{ 41 }$$

and

$$\begin{eqnarray}&&{v}_{\theta }={v}_{r}\displaystyle \frac{r}{R}-{v}_{z}\displaystyle \frac{z}{R}.\end{eqnarray} \tag{ 42 }$$

The density distribution around the ejecta–disk interface (2 × 10¹⁵ cm ≤ r ≤ 4 × 10¹⁵ cm) exhibits nearly evenly spaced perturbations with nonzero angular velocities. Since the ejecta–disk interface is a shear layer, where the radially expanding gas is adjacent to the CSM nearly at rest (middle panel of Figure 12), this perturbed interface is likely produced by Kelvin–Helmholtz instability. The development of the instability makes the ejecta around the interface clumpy, which would affect the emission traveling around the ejecta–disk interface. Although it seems to be the case in this particular simulation, it is not clear whether the instability always develops. The growth of the Kelvin–Helmholtz instability is sensitive to the density and velocity gradients at the shear layer, namely the vertical structure of the CSM disk. As long as the pre-SN mass-loss process responsible for creating confined CSM is unclear, the role played by shear layers is also uncertain.

Zoom In Zoom Out Reset image size

As we have discussed, some hydrodynamic instabilities, the Rayleigh–Taylor, Vishniac, and possibly Kelvin–Helmholtz, produce small-scale structures in the ejecta–CSM interface. These instabilities develop even from numerical errors and thus destroy the equatorial symmetry. As seen in Figures 10 and 11, the spatial distributions of physical variables in the upper and lower hemispheres certainly exhibit small equatorial asymmetries, although their global structures are similar. This leads to slight differences in the light curves seen from the upper and lower hemispheres, as we shall see below.

We also note that the treatment of the upper and lower edge of the CSM disk is also important for the opening angle of the void region. In our simulations, we assume the CSM density structure in Equation (28) with q = 4. However, the realistic density structure of the envelope of a CSM disk would likely depend on its formation process and therefore is highly uncertain.

In Figures 13 and 14, we plot the isotropic equivalent bolometric light curves for the disk CSM models. The light curves with different viewing angles show marked differences from those of the spherical CSM models with the same CSM mass (Figure 9). In addition, some viewing angle effects are clearly recognized. We can divide the presented light curves into two classes: fast and slow risers. The light curves with smaller viewing angles are characterized by a fast rise. On the other hand, an observer around the equator (Θ_obs = 90°) would see less luminous and slowly evolving emission lasting for ∼1 yr.

Zoom In Zoom Out Reset image size

The steep rise seen in light curves with smaller viewing angles can be explained by emission escaping into the regions that are not covered by the CSM disk. Without high-density material preventing radiation from escaping, the observer can see the photosphere located in the SN ejecta directly. The radiation energy initially attached to the SN ejecta is simply released into interstellar space as the ejecta expand. However, the initial radiation energy in the SN ejecta alone cannot explain the luminosity at the initial peak. Therefore, the ejecta–CSM interaction happening around the equator also contributes to the emission along small viewing angles. In other words, the energy dissipated at the ejecta–CSM interface diffuses throughout the SN ejecta and is then released into regions that are not covered by the CSM disk.

Among the light curves in each panel of Figures 13 and 14, those with larger viewing angles (Θ_obs = 90° for Figure 13 and Θ_obs = 60° and 90° for Figure 14) are classified into the slowly evolving case. These light curves exhibit a slow rise followed by a slow decay. As seen in the density distributions in Figures 10 and 11, the lines of sight corresponding to these viewing angles are intercepted by the CSM. Therefore, the observer sees the emission diffusing through the CSM disk.

Zoom In Zoom Out Reset image size

As seen in Figures 10 and 11, the shocked ejecta and CSM do not show strict equatorial symmetry. As a result, light curves with a viewing angle Θ_obs = Θ and its symmetric counterpart Θ_obs = 180° − Θ are not identical to each other. In Figure 15, we compare light curves with Θ_obs = 30° and 150°, 45° and 135°, and 60° and 120° for models with ϑ_disk = 10°. Although each pair of light curves shows slight differences, they mostly behave in a similar way. We have also checked the light curves of the models with ϑ_disk = 20° and confirmed that the effect of the equatorial asymmetry on light curves is less than a factor of 2.

Zoom In Zoom Out Reset image size

Some light curves in Figures 13 and 14 show nonmonotonic temporal evolution. The light curves appear to be bumpier for a larger CSM mass. This is a result of the complex hydrodynamic interaction between the SN ejecta and the inner disk. We consider model D20_M10, which shows the widest variety of light curves among the disk models. Figure 16 shows the spatial distributions of the density and the outgoing flux at several epochs for model D20_M10. The outgoing flux is calculated by Equation (33) with ${l}_{{\rm{v}}}^{i}=(r/R,z/R)$ at each numerical cell. These outgoing flux distributions reflect the distributions of the density and the radiation energy source at the corresponding epoch. At early epochs, the ejecta are interacting with an inner and denser part of the disk; thus, the disk region is not penetrated by radiation, leading to small outgoing fluxes. The outgoing flux is larger around the region with lower inclination angles, which is not covered by the disk. Furthermore, at the earliest epoch (t = 5 × 10⁵ s), even regions close to the symmetry axis, z = 0, show high outgoing fluxes, which indicates that the radiation produced by the interaction region can easily escape into directions with smaller viewing angles. At the early stages of the ejecta–disk interaction, the energy dissipation happens at outer layers of the ejecta. Therefore, the relatively small optical depth of the outer layers makes it easy for the dissipated energy to radiate away into wider solid angles, including inclination angles close to 0°. Then, the interaction region digs through the inner regions of the ejecta (in mass coordinates) and serves as an energy source nearly at the center. Therefore, for ejecta with smaller inclination angles, the energy dissipated at the deeply embedded interaction region should diffuse through the stratified layers. These two effects explain the light curves of D20_M10 with Θ_obs = 0° and 30° (Figure 14). They exhibit an initial spike followed by a slowly evolving second peak. The initial spike is created by the impact of the outer ejecta colliding with the inner disk, which can easily escape into smaller viewing angles. The second peak is powered by photons diffusing from the deeply embedded energy dissipation region all the way toward the outermost layer.

Zoom In Zoom Out Reset image size

The light curve with Θ_obs = 45° behaves in a more complicated way because its viewing angle is close to the orientation of the ejecta–disk interface. As seen in Figure 16, the region with the highest outgoing flux is aligned with the upper and lower disk surfaces. This is naturally expected for the ejecta–disk interaction where the energy dissipation happens most efficiently in the inner disk edge. Above (below) the upper (lower) disk surface with a large outgoing flux, however, regions with lower outgoing fluxes can be found. These regions are associated with the ejecta dynamically influenced by the collision with the inner disk. In the ejecta–disk collision, the ejecta initially expanding in the equatorial direction are pushed away by the inner disk digging through the inner part of the ejecta and accumulate above (below) the upper (lower) disk surface. The resultant region with an enhanced density around the disk surface can be seen in the density distributions in the top panels of Figure 16. In these regions, the boosted optical depth along the radial direction creates "shadowed regions" (the regions with relatively lower F_out extending along θ = 30°–45° and 135°–150°) in the outgoing flux distribution, leading to striped outgoing flux distributions, as seen in the bottom panels of Figure 16. In addition, the orientation of the region with enhanced density can change with time due to the compression of the disk by the ejecta. Accordingly, the shadowed region is closer to the equatorial plane at late epochs. It is how a particular line of sight intersects with striped flux distributions that determines the temporal behavior of the corresponding light curve. As a result of these effects, the light curves exhibit nonmonotonic evolution.

We also note that the outgoing flux distribution at t = 5 × 10⁶ s shows some small-scale structure in the shadowed region. This region corresponds to clumpy ejecta along the ejecta–disk interface. Therefore, this structure may be related to the shear layer affected by the possible growth of the Kelvin–Helmholtz instability, although we cannot exclude the possibility of numerical oscillations. Since the typical length of this small-scale perturbation is less than 10¹⁵ cm, which yields a light-crossing time of less than 1 day, this behavior does not produce the 10 day long bumps observed in light curves, even when it is really a physical phenomenon.