Pair-Instability Supernova Simulations: Progenitor Evolution, Explosion, and Light Curves

TL;DR

This paper presents 1D, 2D, and 3D hydrodynamic simulations of pair-instability supernovae, analyzing their explosion mechanisms, ejecta mixing, and synthetic light curves to compare with observations and explore spectral effects.

Contribution

It introduces the first 3D hydrodynamic simulations of PISNe and investigates how ejecta mixing influences supernova light curves and spectra.

Findings

Slight deviations from spherical symmetry increase with dimensionality.

No significant impact on bolometric light curves from mixing.

Ejecta mixing may affect spectral features, not light curves.

Abstract

In recent years, the viability of the pair-instability supernova (PISN) scenario for explaining superluminous supernovae has all but disappeared except for a few slowly-evolving examples. However, PISN are not predicted to be superluminous throughout the bulk of their mass range. In fact, it is more likely that the first PISN we see (if we have not seen one already) will not be superluminous. Here, we present hydrodynamic simulations of PISNe for four stellar models with unique envelope properties spanning the PISN mass range. In addition, we compute synthetic light curves for comparison with current and future observations. We also investigate, in the context of our most massive model, the prospect of mixing in the supernova ejecta alleviating discrepancies between current PISN models and the remaining superluminous candidate events. To this end, we present the first published 3D hydrodynamic simulations of PISNe. After achieving convergence between 1D, 2D, and 3D simulations we examine mixing in the supernova ejecta and its affect on the bolometric light curve. We observe slight deviations from spherical symmetry which increase with the number of dimensions. We find no significant effects on the bolometric light curve, however we conclude that mixing between the silicon and oxygen rich layers caused by the Rayleigh-Taylor instability may affect spectra.

Tables

Table 1. Table 1: Pre-SN properties for models P150, P175, P200, and P250.

Model	$M_{\text{tot}}$	$M_{\text{CO}}$	Radius	Surface Composition
	(M⊙)	(M⊙)	(R⊙)
P150	90.8	65.7	1267	20% H; 80% He
P175	102.8	81.4	1107	18% H; 82% He
P200	109.9	100.9	80.1	6% H; 94% He
P250	126.7	126.7	2.4	34% He; 39% C; 27% O

Table 2. Table 2: 1D and 3D Explosion Properties for P150, P175, P200, and P250 at the fiducial resolution ( 6.5 × 10 7 6.5 superscript 10 7 6.5\times 10^{7}  cm).

Model	M${}_{\text{CO}}$ (M⊙)	E${}_{\text{exp}}$ (B)		${\rho }_{\text{c}}$ (${10}^6$ g/cm3)		T${}_{\text{c}}$ (${10}^9$ K)		M${}_{\text{Ni}}$ (M⊙)		M${}_{\text{Si}}$ (M⊙)		M${}_{\text{O}}$ (M⊙)
		1D	3D	1D	3D	1D	3D	1D	3D	1D	3D	1D	3D
P150	65.7	5.58	5.68	1.42	1.39	3.41	3.39	$2.86\times {10}^{-3}$	$2.59\times {10}^{-3}$	5.85	5.82	48.5	48.6
P175	81.4	17.7	17.2	2.27	2.21	3.92	3.89	$3.16\times {10}^{-1}$	$3.00\times {10}^{-1}$	12.9	12.8	51.3	51.5
P200	100.9	49.4	48.7	3.92	3.79	4.89	4.85	$1.21\times {10}^1$	$1.16\times {10}^1$	22.4	22.4	42.2	42.6
P250	126.7	82.1	81.7	7.03	6.64	5.83	5.75	$3.43\times {10}^1$	$3.35\times {10}^1$	24.5	24.6	35.9	35.7

Table 3. Table 3: 1D, 2D, and 3D Explosion Properties for P250 at high resolution ( 3.25 × 10 7 3.25 superscript 10 7 3.25\times 10^{7}  cm).

	E${}_{\text{exp}}$	${\rho }_{\text{c}}$	T${}_{\text{c}}$	M${}_{\text{Ni}}$	M${}_{\text{Si}}$	M${}_{\text{O}}$
	(B)	(${10}^6$ g/cm3)	(${10}^9$ K)	(M⊙)	(M⊙)	(M⊙)
1D	81.9	6.90	5.80	34.0	24.5	35.6
2D	81.9	6.85	5.79	34.0	24.5	35.6
3D	81.8	6.77	5.77	33.8	24.6	35.7

Equations

$$E_{\exp }=E_{\mathrm{tot},\mathrm{f}}-E_{\mathrm{tot},\mathrm{i}}+E_{\mathrm{tot},\mathrm{p}},$$

Full text

hydrodynamics, radiative transfer, stars:evolution, stars:interiors, stars:massive, (stars:) supernovae: general

††software: GENEC (Ekström et al., 2012), FLASH (v4.3; Fryxell et al. (2000) & Dubey et al. (2009)), STELLA (Blinnikov et al., 2006)

1 Introduction

Pair-instability supernovae (PISNe) are the explosive deaths of very massive stars (VMS; defined by Vink (2014) as stars with initial masses greater than 100 M*⊙) that produce carbon-oxygen (CO) cores in the mass range 60 M${}_{\odot}\lesssim$ MCO $\lesssim$ 130 M⊙. Stellar models predict that, for non-rotating stars with zero metallicity, this corresponds to a zero-age main sequence (ZAMS) mass range of 140 M⊙* $<$ MZAMS $<$ 260 M*⊙* (Heger & Woosley, 2002). For stars in this mass range, life is cut short when the pair-instability (PI) triggers an implosion of the core shortly after core carbon burning. The implosion is reversed by explosive nuclear burning (of primarily oxygen) which releases enough energy to totally unbind the star. The PI occurs when the radiation pressure in the stellar core is reduced by the reaction $\gamma+\gamma\rightarrow e^{-}+e^{+}$ and was first shown to cause explosions in simulations of isentropic oxygen cores (Rakavy & Shaviv, 1967; Barkat et al., 1967). At the time, detection of such an event was thought to be highly unlikely since massive enough progenitors within the range of detectability were thought to be extremely rare, if they existed at all. Since then, observations of VMS along with advances in stellar physics have made the search much more promising.

A few VMS, have been detected in the Large Magellanic Cloud, specifically in the cluster R136 (Crowther et al., 2010b). In addition, there are several good VMS candidates near the Galactic center (Martins, 2014). However, it is unclear whether or not these stars will be able to retain enough mass to explode as PISNe. They exist in regions where the metallicity is near solar which is thought to drive very high mass-loss rates. Langer et al. (2007) find a limiting metallicity of Z*⊙*/3 above which no star will be able to explode as a PISN. In any case, their existence suggests that similarly massive stars have formed in regions of lower metallicity and in the early Universe.

Including rotation leads to a more chemically homogeneous stellar evolution allowing more of the initial mass to be converted into heavier elements. This effect facilitates the formation of larger carbon-oxygen cores for a given ZAMS mass, shifting the minimum ZAMS mass down from 140 M*⊙* to 65 M*⊙* for stars with initial rotation rates of 80% of Keplerian velocity (Chatzopoulos & Wheeler, 2012a). In other words, in a model for a 65 M*⊙* star rotating at this velocity, all of its initial mass was converted into carbon and oxygen. Additionally, the presence of a magnetic field at the surface of a VMS can quell mass loss rates allowing for the possibility of PISNe even from solar-metallicity progenitors (Georgy et al., 2017). These results substantially lower the bar for finding PISN progenitor candidates.

Many recent numerical simulations agree that the Pop III Initial Mass Function (IMF) is dominated by stars around 100 M*⊙* (Abel et al., 2002; Bromm & Loeb, 2004; Yoshida et al., 2008). However, some very recent simulations find that fragmentation may lead to stellar populations that extend below even 1 M*⊙* (Stacy & Bromm, 2014; Stacy et al., 2016) which means VMS would be significantly less abundant than previously thought. Additionally, it is unclear how massive an individual Pop III star can get. The maximum mass may be limited by ionization feedback (Krumholz, 2014). If the maximum mass is as high as it appears to be in the local Universe, then VMS are likely to exist in the early Universe. Furthermore, such stars would experience lower mass-loss rates due to the absence of metals in their atmospheres, allowing them to retain enough mass to explode as PISNe. There is also numerical evidence that suggests pockets of pristine star-forming gas exist even at relatively low redshift ($2<z<5$) (Tornatore et al., 2007) which increases the prospects of PISN progenitor stars existing in the local Universe.

PISNe are expected to produce a wide variety of SN types as well as span a large range in peak luminosity. On the lower end of the PISN mass range, the explosion energies can be only a few Bethe (1 Bethe equals $10^{51}$ erg) and the nickel yields may be less than those of ordinary core-collapse supernovae. Consequently, the light curves (LCs) and spectra may resemble those resulting from other supernova mechanisms. For example, red supergiant and stripped core PISN progenitors would likely look like long-duration luminous Type II-P and Type Ib/Ic SNe, respectively (Kasen et al., 2011; Kozyreva et al., 2014a). Low-mass PISNe may even explain some “.Ia supernovae” (named as such because their explosion strength is 1/10th that of ordinary supernovae) (Whalen et al., 2014).

Conversely, near the upper end of the PISN mass range the explosions can be extremely energetic. The yield of radioactive nickel in such explosions can approach 55 M*⊙* (Heger & Woosley, 2002) causing a very luminous long-duration SN. It is this potential for high luminosity that made high-mass PISNe an attractive model for Superluminous Supernovae (SLSNe) when they (SLSNe) were first observed about a decade ago.

SLSNe are defined as any SN with a peak absolute magnitude brighter than -21 (Gal-Yam, 2012). They are classified as either SLSN-I (for events without hydrogen) or SLSN-II (for events with hydrogen) just as in the classification scheme for normal SNe. As is the case for normal SNe, SLSNe display a large amount of diversity within the two main types. The PISN model is most well-suited to explain the slowest evolving SLSNe-I that exhibit post-peak decline rates consistent with the radioactive decay of 56Ni to 56Co and 56Fe (Inserra et al., 2017; Jerkstrand et al., 2017).

Gal-Yam et al. (2009) proposed a PISN explanation for such a SLSN, namerly SN 2007bi (although Woosley et al. (2007) had made the case that SN 2006gy was produced by circumstellar interaction with shells originating from the pulsational pair-instability). SN 2007bi was discovered in a relatively nearby dwarf galaxy ($z=0.1279$) which means that, if it was a PISN, VMS can form and retain enough mass to explode as PISNe in the local Universe. The PISN interpretation of SN 2007bi was both supported Kasen et al. (2011); Kozyreva et al. (2014a) and critiqued Dessart et al. (2012, 2013); Chatzopoulos et al. (2015); Jerkstrand et al. (2016) by subsequent works. PISN models were able to sufficiently match key observables such as the bolometric light curve and photospheric velocity, however they were not able to explain the blue nebular spectra of SN 2007bi.

More recently, PTF12dam has captured the interest of PISN enthusiasts. PTF12dam has late-time (Nicholl et al., 2013) and host galaxy (Chen et al., 2015) properties that are very similar to those of SN 2007bi but, unlike SN 2007bi, it was caught before peak luminosity (Quimby et al., 2012). The relatively fast rise to peak luminosity together with the spectral evolution over this period pose serious problems for PISN models (Nicholl et al., 2013). However, as shown by Kozyreva et al. (2017) with two of the models used here (P200 and P250), stripped-envelope PISN models at relatively high metallicity ($Z=0.001$) predict shorter rise times and higher color temperatures than their Pop III cousins. Here, we extend the mass range to include two lower mass models (P150 and P175) to cover the PISN mass range. In addition, we extend the dimensionality of our simulations to 2D and 3D in order to examine the effects of mixing of 56Ni in the SN ejecta and its observational consequences.

We describe our methods in Section 2 including the stellar evolution, hydrodynamic (in 1D, 2D, and 3D) and radiation-hydrodynamic simulations. In Section 3, we present the explosion properties from the hydrodynamic simulations and in Section 4 we discuss the lightcurves. Finally, we summarize the main points of the paper in Section 5.

2 Inputs and Numerical Setup

2.1 Stellar Models

Four VMS models were computed with the GENEC stellar evolution code (Ekström et al., 2012) and with the same input physics as in Yusof et al. (2013). GENEC uses adaptive spatial and temporal resolutions. Spatial resolution is set to resolve gradients of key quantities like temperature and hydrogen content. Around 200 zones are used on the ZAMS and 300-800 zones are used towards the end of the evolution (more zones are used for more extended/cooler envelopes). Models took between 25,000 and 30,000 time steps from start to end. Mass loss strongly affects the evolution of very massive stars. We therefore list the prescriptions used to calculate the models in this study. For main-sequence stars, we used the prescription for radiative line driven winds from Vink et al. (2001), which compares rather well with observations (Crowther et al., 2010a; Muijres et al., 2011). For stars in a domain not covered by the Vink et al. (2001) prescription ($\log(T_{\text{eff}})<3.9$), we applied the de Jager et al. (1988) prescription to models with $\log(T_{\text{eff}})>3.7$. For $\log(T_{\text{eff}})\leq 3.7$, we performed a linear fit to the data from Sylvester et al. (1998) and van Loon et al. (1999) (see Crowther, 2001). The formula used is given in Eq. 2.1 in Bennett et al. (2012). In the calculations, we consider a transition from O-type or giant to Wolf-Rayet (WR) star when the surface hydrogen mass fraction, $X_{s}<0.3$ and the effective temperature, $\log(T_{\text{eff}})>4$. The mass loss rate used during the WR phase depends on the WR sub-type. For the eWNL phase (when $0.3>X_{s}>0.05$), the Gräfener & Hamann (2008) recipe is used (in the validity domain of this prescription, which usually covers most of the eWNL phase). In many cases, the WR mass-loss rate of Gräfener & Hamann (2008) is lower than the rate of Vink et al. (2001), in which case, we used the latter. For the eWNE phase – when $0.05>X_{s}$ and the ratio of the mass fractions of $({{}^{12}\text{C}}+{{}^{16}\text{O}})/{{}^{4}\text{He}}<0.03$ – and WC/WO phases – when $({{}^{12}\text{C}}+{{}^{16}\text{O}})/{{}^{4}\text{He}}>0.03$ – we used the corresponding prescriptions of Nugis & Lamers (2000). Note also that both the Nugis & Lamers (2000) and Gräfener & Hamann (2008) mass-loss rates account for clumping effects (Muijres et al., 2011).

The metallicity dependence of mass loss rates is included in the following way. The mass loss rate used at a given metallicity, $\dot{M}(Z)$, is the mass loss rate at solar metallicity, $\dot{M}(Z_{\odot})$, multiplied by the ratio of the metallicities to the power of $\alpha$: $\dot{M}(Z)=\dot{M}(Z_{\odot})(Z/Z_{\odot})^{\alpha}$, where $\alpha$ was set to 0.85 for the O-type phase and WN phase and 0.66 for the WC and WO phases following Eldridge & Vink (2006). Note that for WR stars the initial metallicity rather than the actual metallicity was used in the equation above. The parameter $\alpha$ was set to 0.5 for the de Jager et al. (1988) prescription. Finally, $\alpha$ was set to 0 (no dependence) if $\log(T_{\text{eff}})\leq 3.7$ (note that none of the models presented in this study reach such low effective temperatures).

All models are non-rotating and have an initial metallicity of $Z=0.001$. Considering solar composition to be $Z=0.014$, this means that the initial metallicity $Z\simeq 0.07\,Z_{\odot}$. Their ZAMS masses are given in the names of the models (P150, P175, P200, and P250) and their pre-SN properties are listed in Table 1. The masses were chosen to span the PISN mass range given in Heger & Woosley (2002), that is 140 M*⊙* $<$ MZAMS $<$ 260 M*⊙*.

The evolution of the models is presented in Figures 1 (structure), 2 (HRD), 3 (left: mass loss; right: Eddington parameter), and 4 (central conditions). The models have very large convective cores and are very luminous, which is typical for very massive stars. This leads to very strong mass loss ranging between $10^{-6}$ and $10^{-1}$ solar masses per year. The peak in mass loss, seen around log(time left)$\sim 5.4$ in Figure 3 (left) corresponds to when the model reach cool parts of the HRD, for which the de Jager et al. (1988) prescription is used. This is an empirical prescription, which mimics the strong mass loss experienced by luminous blue variable stars. The zigzag pattern (repeated spikes, best seen in model P150 at the end of its evolution) is due to the star getting cooler than $\log(T_{\text{eff}})<3.9$, the mass loss rates switches from the Vink et al. (2001) to the de Jager et al. (1988). This leads to a sharp increase in mass loss, which causes a contraction of the surface back to hotter temperature. Mass loss decreasing by to the Vink et al. (2001) prescription, the star expands again and the cycle continues creating a zigzag pattern. The very strong mass loss explains why model P250 loses not only most of its hydrogen-rich envelope but also most of its helium envelope, ending its life as a compact CO core. Model P250 is very similar to more metal-rich models as those presented in Yusof et al. (2013), to which we refer the reader for more details about the evolution of the structure and mass loss in very massive stars. The absence of an extended envelope surrounding the exploding core should lead to a much faster rise to peak luminosity, closer to that of some of the slowest rising SLSNe (Kozyreva et al., 2017).

As the models in this study have a sub-solar metallicity ($Z=0.001$), models P150, P175, and P200 manage to retain a fraction of their hydrogen-rich envelope. They thus are located in cooler parts of the HR-diagram and have larger radii at the end of their evolution than model P250. This is more typical of very low and metal-free models of VMS (see Hirschi, 2007; Ekström et al., 2008; Yoshida & Umeda, 2011; Yoon et al., 2012; Chatzopoulos & Wheeler, 2012b; Kozyreva et al., 2017, and references therein)

As expected the fraction of the initial mass lost due to stellar winds increases with initial mass due to the dependence of the mass loss on luminosity. Mass loss prescriptions for VMS are still uncertain though and VMS get close to the Eddington limit towards the end of their evolution. Gräfener et al. (2011) suggested enhanced mass-loss rates (with respect to Vink et al., 2001) for stars with high Eddington parameters ($\Gamma_{\rm Edd}\geq$ 0.7) that they attribute to the Wolf-Rayet stage. In the present work, we did not use an increased mass loss rate close to the Eddington limit. In order to know whether it would have had an impact on the models, we discuss here the proximity of our models to the Eddington limit. Figure 3 (right) shows the evolution of the Eddington parameter, $\Gamma_{\text{Edd}}=L/L_{\text{Edd}}=\kappa L/(4\pi cGM)$. The initial values for $\Gamma_{\text{Edd}}$ range between $0.4-0.6$, so well below the Eddington limit, $\Gamma_{\text{Edd}}=1$, and below the limiting value of 0.7 where enhanced mass-loss rates are expected according to Gräfener et al. (2011). As the evolution proceed, however, the Eddington parameter increases to values above 0.7. Additional mass loss may thus be able to remove the rest of the hydrogen-rich envelope during the late stages, even at very low metallicities. The reader is referred to the recent book on VMS for more detail (Vink, 2015).

The evolution of the central conditions is shown in Figure 4. The models were evolved until at least carbon burning. The end point of the track is the point at which the models were mapped into FLASH.

2.2 Hydrodynamic Simulations in 1D

The stellar cores from the models described in Table 1 were mapped into hydrodynamics code FLASH (version 4.3) (Fryxell et al., 2000; Dubey et al., 2009) during core carbon burning, at the end of the evolutionary tracks shown in Fig. 4. Note that Figure 4 shows the time evolution of the central properties of the models (not the profiles within the star). Figure 5 shows the profiles of the adiabatic index for the GENEC model mapped into FLASH. As seen in Figure 5, the instability develops first off-center and not in the very center. This explains why the tracks of the central properties (Figure 4) do not reach the unstable (grey) region. This is due to neutrino cooling being stronger in the very center during the contraction of the core after He-burning. With the exception of P150, the mapping to FLASH was done at the point during core carbon burning at which the cores have become slightly unstable due to the PI. P150 constitutes a special case in which the GENEC model crashed before reaching the instability (see Figure 5). The difficulty in evolving this particular model is due to the envelope of the model expanding to large radii following the core contraction at the end of core He-burning, whereas the more massive models (with little or no H-rich envelope) remain compact at the end of their evolution. The choice of the evolutionary stage for the mapping is critical to generating an explosion with minimum error in stellar evolution. GENEC does not include the contribution to the pressure from electron-positron pairs so the mapping must be done before this contribution becomes too great. Conversely, FLASH is not a stellar evolution code so the input model must be sufficiently evolved for collapse to occur. We were able to achieve collapse in FLASH for models in which the initial pair pressure barely exceeded 1% of the total pressure in the core. For P150, we had to use a somewhat different method for mapping. We will first explain the method for exploding the three fully evolved models: P175, P200, and P250.

Only the inner cores of radii $3.334\times 10^{10}$ cm (P175), $4.167\times 10^{10}$ cm (P200), and $5.000\times 10^{10}$ cm (P250) respectively, comprising the CO cores plus a small part of the helium shells were mapped into FLASH. For brevity, we will call these inner spheres the ‘cores’. The radii of the ‘cores’ ($R_{\text{core}}$) were chosen so that the ratio $R_{\text{core}}/R_{\text{CO}}$ was about 1.3. We mapped the initial models using the same scheme as in Chatzopoulos et al. (2014)

We use the new directionally-unsplit hydrodynamics solver (Lee & Deane, 2009) coupled with the Aprox19 nuclear burning network (Timmes & Swesty, 2000) which are both included in FLASH. We employ the Helmholtz equation of state (Timmes & Swesty, 2000) and the block-structured grid implementation Paramesh4dev (which includes Adaptive Mesh Refinement; AMR). Self-gravity was computed with the new Multipole solver in FLASH (Couch et al., 2013). The grids were set to an initial minimum resolution of $1.3\times 10^{8}$ cm with the freedom to refine once, in response to high enough density and/or temperature gradients, to $6.5\times 10^{7}$ cm. We employ the ‘reflect’ (‘diode’) boundary condition for the inner (outer) boundary.

The minimum resolution specified is higher than would normally be required for a PISN simulation. Here, however, we had noticed that mass would slowly flow in through the outer zone in FLASH during the time leading up to the collapse (even though the ‘diode’ boundary condition is meant to prevent this). Since we had to evolve the star in FLASH for a considerable amount of time before the collapse occurred, a high minimum resolution was required so that, by the time collapse occurred, only a few solar masses of material had been added to the simulation. The refinement criteria were set to half of their default values in FLASH which allowed the maximum refinement level to be achieved in the inner cores during the explosive burning phases. This was important for determining the energetics and nucleosynthetic yields of the explosion. During explosive burning, the shock wave is launched from just outside the exploding core. The simulation is halted just before the shock wave exits the domain (at this point nuclear burning rates are negligible).

To explode P150, we used the same resolution and refinement settings as we did for the rest of the models. The difference is an extra step in between the stellar evolution and the explosion simulation where we use FLASH, but with a different hydrodynamics implementation, to prime the model for collapse. We took the final P150 time step from the stellar evolution (shown in Figure 5) and mapped the inner $1.67\times 10^{11}$ cm (including the carbon-oxygen core plus part of the helium shell) into FLASH to evolve towards the instability using the directionally-split hydrodynamics solver. This solver was better suited for the slow initial approach to the PI. We then re-mapped the unstable core into FLASH, this time evolving with the new unsplit hydrodynamics solver as for the other three models to follow the collapse and explosion. If the transition is made too early the core will not collapse. If the transition is made too close to collapse the small numerical effects from switching solvers do not have time to dissipate. Thus, we transition at the earliest time that yields a collapse and explosion. This corresponds to a time when the minimum adiabatic index in the core is 1.302. The subsequent simulation then takes 650.1 s to reach maximum compression, allowing enough time for the core to relax out the numerical effects from the transition. With the additional mapping required for model P150, the uncertainties are more difficult to quantify. However, we stress that the carbon-oxygen core, the mass of which is an excellent predictor of the explosion properties, was fully formed during the stellar evolution with GENEC. In addition, the collapse and explosion phases were fully simulated using the unsplit hydrodynamics solver, which is consistent with the three more massive models.

2.3 Hydrodynamic Simulations in 2D and 3D

Simulations in 1D, 2D, and 3D with each of the four models were done beginning from profiles taken from the 1D simulations described in Section 2.2 just before collapse (about 20s prior to maximum compression). We used the recommended geometry settings in FLASH for each dimensionality: “spherical” for 1D, “cylindrical” for 2D, and “Cartesian” geometry for 3D. Since we were mapping a spherical grid onto non-spherical geometries in the 2D and 3D cases, a slightly smaller domain was evolved (a cube (3D) or half of a square (2D) that could fit inside a sphere (3D) or semicircle (2D) of radius equal to the domain size of the 1D simulation). It was important that the original 1D simulations included a small part of the envelope in their ‘cores’ so that the domains in 2D and 3D could still contain the entire CO cores. The dimensions of our new grids are in 1D: $0<r<\text{xmax}$ (1D), in 2D: $0<r<\text{xmax}$, $-\text{xmax}<z<\text{xmax}$, and in 3D: $-\text{xmax}<x,y,z<\text{xmax}$, where for each model xmax had to be less than $1/\sqrt{3}$ ($1/\sqrt{2}$) times the radius of the 1D simulation so that the initial state of the 3D cube (2D half square) could be completely specified by the spherical input model. For P150, P175, P200, and P250 xmax was set to $2.500\times 10^{10}$ cm, $1.670\times 10^{10}$ cm, $2.083\times 10^{10}$ cm, and $2.500\times 10^{10}$ cm, respectively.

We elected to use a fixed grid for the multidimensional simulations after noticing that derefinement led to spurious mixing at shell interfaces. This fixed grid consisted of nested cubes (3D) or half squares (2D) of different, but fixed, resolution settings. In 3D, the innermost cube (with edge length of $1.0\times 10^{10}$ cm) had the maximum resolution used. This innermost cube was surrounded by a second, larger cube (edge length of $1.5\times 10^{10}$ cm) which has the next lower refinement level (one factor of two) in the volume not occupied by the innermost cube. The second cube is again surrounded by another (third) cube (edge length of xmax) which has again the next lower refinement level in the volume not occupied by the second cube. In 2D, we used a similar hierarchy of half squares in which the longer edges of the half squares corresponded to the cube edge lengths given above. We use the ‘diode’ boundary condition for all boundaries in 3D and for all boundaries in 2D except the inner boundary for which we use ‘reflect’. We ran each simulation until the shock reached the edge of the grid which was long after all significant nucleosynthesis had occurred.

To facilitate the comparison of our multidimensional results with our 1D results, we further processed the 3D data by computing angular averages for all the variables. The angular average is a two-step process. First, the state variables are averaged over a ‘block’ ($4^{n}$ cells, where $n$ is the number of dimensions). Then, we perform another mass-weighted average on to a coarser grid in which the bin widths are $1.5\times 10^{8}$ cm. During this second average, the $5^{th}$ and $95^{th}$ percentile values for each variable in each bin are also computed. This last step allows us to see how wide the range of values for a particular variable can be within a single radial bin and is used for the shaded regions in Figure 8 in Section 3.2.

2.4 Re-appending the Envelope

For the purpose of comparing to current and future observations, we computed synthetic LCs for all our models. For accurate light curves, it was necessary to include the entire star in the simulations. On the timescale of our hydrodynamic simulations (several $\times 10^{4}$ s) we did not expect the envelope to change during the collapse of the core, so we re-appended the pristine envelope from the stellar evolution simulation to the exploding ‘core’ for each model. During the hydrodynamic simulation, the outer edge of the ‘core’ decreased in density and temperature causing a discontinuity between the final ‘core’ profile and the envelope profile. However, we needed to join the two profiles in a smooth manner while still preserving the outer envelope structure which is of the greatest importance to the light curve.

For the 1D simulations, we chose two points between which the density and temperature would be artificially set to follow a straight line connecting the two points in the $\log(\rho)$-radius or $\log(T)$-radius plane. For this modified region, the composition is uniform so the mass fractions are all set according to this uniform composition. One point was chosen to be in front of the shock. The other point was chosen at a point in the envelope and had to fulfill two criteria: (i) the local slope at this point has to be comparable to the slope of the connecting line and (ii) it has to be far enough from the surface so that the structure of the outer envelope was preserved. This method is illustrated using the density profile of P250 as an example in Figure 6. The mass and internal energy lost in this process is comparable to but slightly lower than the mass (and accompanied internal energy) that got added through the outer boundary prior to collapse in the FLASH simulations. The total mass (energy) changes by less than 3% (2%) by this procedure.

In addition to the method described above, we investigated two additional limiting cases. For this test we used a slightly earlier time step than described above because it was necessary in order to apply the first method. Method A involved flattening the density and temperature profiles in the region between the shock and the edge of the core to match the values at the inner edge of the envelope. This method added a small amount of material ($\sim 1M_{\odot}$) and a small amount of internal energy ahead of the shock. Method B was as described in the paragraph above, except that the exploding core remained untouched. Instead, the density and temperature in the inner region of the envelope were replaced by a straight line connecting the outermost point in the core to the same fixed point in the envelope from before. This method caused a loss of a similarly small amount of material (again, $\sim 1M_{\odot}$) and internal energy ahead of the shock as method A. For both methods, the resulting profiles were mapped back into FLASH for further evolution and then mapped into STELLA (Blinnikov et al., 2006) using the standard procedure described below and in Section 4) to compute light curves. The resulting light curves were almost identical, indicating that both methods are suitable for this study. Our standard method for re-appending the envelope represents an intermediate method between method A and B, in which matter and energy are removed ahead of the shock (shown by the red line in Figure 6).

In 3D, the simulations were computationally too expensive to include the entire star. In this case, we conjoined an angular average of the ‘core’ with the pristine envelope using the same method as in the case of the 1D simulations. The only difference is that the angular-averaged ‘core’ profile from the 3D simulation represents a sphere of $2.5\times 10^{10}$ cm, whereas the pristine envelope begins at $5\times 10^{10}$ cm due to mapping from a sphere (on a spherical 1D grid to a Cartesian cube in 3D, as described in Section 2.3). We artificially set the densities and temperatures as described above (linear slope in the $\log{(\rho)}-r$ and $\log(T)-r$ plane, respectively). The velocities in this region are set to zero as during stellar evolution. For the composition we choose the composition of the envelope (which is almost perfectly uniform throughout the envelope). While this step erases information of this intermediate region, it preserves the information from the multi-dimensional simulation of the ‘core’, such as the amount of outward mixing of nickel.

After re-appending the stellar envelope, the entire star was further evolved in FLASH before mapping to the radiation-hydrodyanmics code STELLA to compute light curves (see Section 4). For this, we use a maximum resolution of $6.5\times 10^{7}$ cm and relax the refinement criteria back to their default values. We also allow a lower minimum resolution of $5.2\times 10^{8}$ cm (which is lower than in the simulations of the ‘cores’ only) since the low densities at the outer edge of the envelopes effectively remove the problem of mass inflow at the outer boundary of the computational domain. As before, we employ the ‘reflect’ (’diode’) boundary condition for the inner (outer) boundary. We choose to map to STELLA after the shock has traversed half the radius of the star. The light curve calculations with STELLA described in Section 4 use these final profiles.

3 Explosion Properties

3.1 PISN Explosions in 1D

We present explosion properties of all four models considered here in Table 2. The Table also includes yields of Ni, Si, and O as calculated within FLASH. Between the lowest mass model (P150) and the highest mass model (P250) the nickel yield for the 1D simulations increases by more than four orders of magnitude. This is a consequence of the steep dependence on density and temperature of the associated nuclear reactions together with the dependence of the strength of the collapse on the mass of the CO core. The dependence on CO core mass of the results are in good agreement with previous work (Heger & Woosley, 2002; Kasen et al., 2011; Dessart et al., 2013; Whalen et al., 2013, 2014; Kozyreva et al., 2014b; Chatzopoulos et al., 2015).

Silicon is produced from primarily oxygen (and some carbon). Then, if the temperatures and densities get high enough (as in P200 and P250), nickel is produced in substantial quantities primarily from fusion of silicon. To first order, the explosion converts some number of solar masses of oxygen into silicon and nickel. Thus, the final oxygen masses are substantially lower than their pre-explosion values while both the silicon and nickel yields are substantially higher than their pre-explosion values. Note however that a substantial amount of oxygen remains unburnt in each case and thus constitutes a large fraction of the ejected mass.

Figure 7 shows the compositional profiles of the ejecta after evolution in FLASH is complete. The wide range in nickel yield between the models is evident from the Figure as the higher mass models reach core conditions sufficient for explosive silicon burning in addition to the explosive oxygen burning which is achieved in all the models. Also note how the helium appears in the core for the higher mass models reflecting the dissociation of nickel into helium that occurs due to increasingly high core temperatures. This is the mechanism by which cores of higher mass will collapse directly to black holes. The unburnt oxygen is confined to a shell in between two regions where nucleosynthesis occurred. The region interior to the unburnt oxygen shell experiences the bulk of the nuclear burning, however, there is also a small region (in terms of mass coordinate) outside of this shell in which oxygen captures alpha particles from the inner edge of the stellar envelope. Such burning produces mainly silicon and is triggered by shock heating.

3.2 PISN Explosions in Multi-D

In addition to the simulations in 1D, we have also performed simulations at the same fiducial resolution ($6.5\times 10^{7}$ cm) in 2D and 3D. For the rest of the paragraph, we focus on the 3D simulations. The 2D simulations yield very similar results for all the properties in Table 2. In 3D (and also in 2D) we see a slight decrease in nickel yields (and overall explosion strength) for all models (see Table 2) at the fiducial resolution of $6.5\times 10^{7}$ cm. In 3D, the collapse and explosive burning phases occur with very little asphericity which means that the difference in explosion strengths is likely a resolution effect. The 3D simulations utilized a fixed grid but had the same maximum resolution as the 1D simulations. Although the simulations used the same maximum resolution, the computational cells have different shapes. In 1D the cells are simply radial bins. In 3D, they are cubes whose orientation with respect to the radial direction depends on their position in space. The effective radial resolution is thus spatially dependent and is at best equal to the 1D resolution. This leads to the expectation that higher resolution is required to achieve convergence in 3D.

To test this, we increased the resolution everywhere in the grid by a factor of two and ran the simulation again in 1D, 2D, and 3D. Because of limitations in computing time, we restrict ourselves to the most massive model, P250. Table 3 shows the explosion properties in 1D, 2D, and 3D at the new resolution of $3.25\times 10^{7}$ cm. At this resolution we are fully converged in explosion energy, silicon yield, and nickel yield between the 1D, 2D, and 3D simulations. We calculate explosion energies similarly to Chatzopoulos et al. (2015) (see their Table 1):

[TABLE]

where the difference between $E_{\mathrm{tot,f}}$ and $E_{\mathrm{tot,i}}$ is in essence the energy released by nuclear burning and $E_{\mathrm{tot,p}}$ is the initial negative binding energy of the progenitor model. The total energies include contributions from the kinetic, internal, and gravitational potential energies, where the internal energy receives all energy released by nuclear reactions.

We now turn our attention to the multidimensional effect of mixing in the supernova ejecta. The mixing is caused by the growth of the Rayleigh-Taylor (RT) instabilities at compositional boundaries. Chatzopoulos et al. (2014) and Chen et al. (2014) have seen this effect in their 2D simulations. In our simulations, mixing sets in after maximum compression in the expanding ejecta. The instabilities grow for about 15s at which point the shock reaches the edge of the grid and our simulations are halted. Figure 8 shows the 2D (Left) and 3D (Right) angular averaged profiles at the end of the simulation for model P250 at the increased resolution. The shaded vertical width for each nuclear species is computed for each bin in the angular averaging scheme described in Section 2.3 from the distribution of mass fractions involved in the bin average. The lower (upper) edge of the shaded region corresponds to the mass fraction that is greater than 5% (95%) of the values occurring in the bin average. Thus, the vertical width can be thought of as a measure of the degree of asphericity for each mass coordinate. Comparing the left and right sides of the figure shows that mixing is stronger in 3D than in 2D. This result is not so surprising as it is generally agreed that the RT instability grows faster in 3D than in 2D, at least initially (Kuchugov et al., 2014). Furthermore, Calder et al. (2002) find that this is true in FLASH for single-mode perturbations.

On the right panel of Figure 8, the most apparent feature is the mixing around the Si-O interface. In addition, significant smearing of the interface is apparent in the relatively shallow mass fraction gradients present. We also note weak mixing at the Ni-Si interface, the effect of which is much too small to have an appreciable effect on the LC rise time (see Section 4 and Figure 15). It should be noted that the multidimensional phases of the simulations only followed the evolution until just before the shock exits the grid along the axes (when the shock radius passes 2.5 $\times$ 1010 cm) while the total radius of the model is 1.67 $\times$ 1011 cm. Thus, the amount of mixing seen should be regarded as a lower limit to the amount we would expect from the simulation at shock breakout.

In addition to the strength, the character of the mixing differs between 2D and 3D at the Si-O interface. Figure 9 shows the 28Si mass fraction color coded on a log-scale for the 2D (left) and a slice of the 3D (right) simulations at the end of the simulation (the same data that was used to generate the 28Si data for Figure 8). Each panel is centered on the most unstable layer (the Si-O interface). Also visible in the top-right corner of each panel is the inner edge of the outer silicon shell produced via shock burning (see Section 3.1). Examining the left panel of Figure 9 we note that, in 2D, very thin RT fingers exist on both the inside (the Ni-Si interface) and outside (the Si-O interface) of the Si-rich shell. Focusing now on the right panel, we see that the 3D simulation exhibits similar behavior to the 2D simulation at the Ni-Si interface but very different behavior at the Si-O interface where we see RT plumes whose growth results in much stronger mixing.

Also evident from Figure 9 is the angular dependence of the RT features. There is a distinct lack of RT plumes directed along the axes. The Cartesian grid has certainly imposed some numerical artifacts on the data. If mixing at this interface occurs in nature, as we expect, there would be no such angular dependence. However, it is unclear if there is an artificial suppression of RT mixing along the axes or if RT mixing is enhanced off axis. Simulations using a spherical geometry may help to shed light on this issue. Additionally, there is evident substructure in the RT plumes. We provide Figure 10 for a closer look at the plume at 45*∘* in the right panel of Figure 9. At this scale one can even make out the computational cells. Even finer features may develop in a similar simulation with higher resolution that could further change the character (and perhaps also the strength) of the mixing. We leave this prospect for future work.

The contrasting character of the RT mixing between 2D and 3D is also evident from the densities at the Si-O interface. Figure 11 gives a zoomed-in view (the same view as in Figure 10) but with density color coded on a log-scale for the 2D (left) and the same slice of the 3D (right) simulations. We see similar thin fingers and plumes in 2D and 3D, respectively. Note that the underdense regions correspond to the outgoing ${}^{28}Si$ fingers and plumes. The overdensities of these features are on the order of half of a percent in 2D and 5% in 3D. Unlike our result that the mixing is stronger in 3D than in 2D, the larger scale of the RT features in 3D compared to 2D is quite surprising and needs further investigation.

4 Light Curves

We used the 1D radiation-hydrodynamics code STELLA (see Blinnikov et al. (1998, 2000); Kozyreva et al. (2017) for details) to follow the post-explosion evolution of the PISN ejecta and for calculating light curves. For this, we mapped the output from the FLASH simulations into STELLA when the shock was halfway in the hydrogen-rich envelope in P150/P175 models, and just before shock breakout in P200/P250 models. When velocity gets close to the speed of light in P200 and P250, we cut the FLASH profiles at about 30% of speed of light, because STELLA does not include relativistic corrections for the radiative transfer equations. Using this procedure, we calculated bolometric and broad-band light curves for all four models (P150, P175, P200, and P250). We present synthetic LCs from the 1D explosions for all four models. In addition, we show the LC from the fixed-grid, 3D simulation of model P250 using the highest maximum resolution. Since STELLA is a 1D code, an angular average of the 3D FLASH simulation of model P250 was used for this LC (see Section 2.4).

The first light indicating the explosion is radiation emitted as the shock breaks out on the surface of the progenitor. The luminosity at shock breakout depends mostly on the energy of the explosion, the duration depends on the radius of the progenitor. In Figure 12, we present shock breakouts computed with STELLA111We estimate the duration of the shock breakout event the same way as in Kasen et al. (2011), i.e. the full width at half-maximum.. The shock breakout lasts for a fraction of a minute for the compact model P250 (0.5 s). At redshift z=10 it will last about 5 s (according to the cosmological time dilation) and will appear in visual or infrared. As for the P200 model, shock breakout lasts 0.7 min (i.e. 7 min at z=10). It would be quite challenging to detect such events. Models P150 and P175 are extended (1107 R⊙ and 1267 R⊙) and have longer shock breakout durations: 4.2 h and 1.7 h for P150 and P175, respectively. The color temperature of the shock breakouts are: $7.8\times 10^{5}$ K (67 eV) for P150 and P175, $6.5\times 10^{6}$ K (560 eV) for P200, and $2.2\times 10^{6}$ K (200 eV) for P250.

In Figure 13, we present bolometric light curves of all four models222The light curve data are available via http://www.astro.keele.ac.uk/~kozyreva/LCindex.html.. The peak luminosity varies according to the amount of radioactive nickel produced in the PI explosion, i.e. from $10^{41}$ erg s*-1* for P150 (M${}_{\mathrm{Ni}}=0.003$ M⊙) to $3\times 10^{44}$ erg s*-1* for P250 (M${}_{\mathrm{Ni}}=33.5$ M⊙). According to the well-established criterion, only P250 is bright enough to be considered a candidate to explain some SLSNe. To generalize, only PISN models with CO-cores above about 110 M⊙, i.e. which produce more than 15 M⊙ of 56Ni, may result in SLSN-like events. Our P200 model looks like a luminous (but not superluminous) SN, while P175 stays around a luminosity of $10^{42}$ erg s*-1*, which is typical for type II SNe. One of the main properties of our PISN models is the width of the light curves. The light curves rise and slowly decay during 200 days, because of high ejecta mass, about 100 M⊙. The nickel-powered maximum phase resembles the plateau phase of a type IIP SN where the luminosity changes less than 1 magnitude over hundreds of days. Light curves of P175, P200, and P250 peak around day 110, i.e. relatively soon after the explosion compared to the previously published PISN light curves. We refer the reader to Kozyreva et al. (2017) for a discussion of the reason for the 100-day rise-time, which is shorter by about 50 days than previously published PISN light curves (Kasen et al., 2011; Dessart et al., 2013). Our P200 and P250 models are hydrogen-free, therefore, once exploded they appear as hydrogen and helium-free SNe (Type Ic), as photosphere is located deep in the hydrogen and helium-free layers of SN ejecta. The P150 light curve is governed by hydrogen recombination during the first 60 days. The photosphere recedes through the extended hydrogen-helium envelope, therefore the luminosity is relatively high. Later, the P150 light curve flattens because it is supported by oxygen recombination, and luminosity drops as oxygen is located deeper in the ejecta.

In Figure 14, we present light curves for our PISN models in U, B, V, and R broad bands. P150 is an outlier compared to other three models, as P150 produces only 0.003 M⊙ of 56Ni and has no Ni-powered maximum. Instead, P150 has a 50-day earlier maximum supported by a recombination wave receding through the hydrogen-helium envelope. When compared to the $150$ M*⊙* red supergiant model 150M by Kozyreva et al. (2014a), the recombination phase is half as long because our P150 model is relatively less extended (1200 R⊙ for vs. 3500 R⊙ for 150M) and contains 3.8 M⊙ of hydrogen versus 4.9 M⊙ in their 150M model (helium yields are 20 M⊙ and 24 M⊙ for P150 and 150M, respectively). The peak magnitude is still high in U band reaching –19 mags (–18 mags in B, V, and R), which is similar to plateau magnitudes of bright SNe IIP (Schmidt et al., 1994). The other three models (P175, P200, P250) have similarly shaped LCs with a broad peak. P200 and P250 are very bright in all bands, and even brighter in U band rather than in B and V bands. P175 has a short recombination phase similar to P150 (15 days above $-16$ mags in BVR bands; 30 days in U band) with the peak magnitude $-19$ mags in U band and $-18$ mags in other bands. If discovered during the first 30 days P175 will appear as hydrogen-rich type II SN, while it will be classified as type Ib if discovered at the day 110, i.e. around the main Ni-powered maximum.

Finally, in Figure 15, we compare the bolometric light curves for model P250 from the 1D and 3D explosions. We find that even though there is some degree of mixing at the Ni-Si and interfaces seen in the right panel of Figure 8, this does not translate to a significant difference in the bolometric light curves.

5 Summary

PISNe, first theorized in 1967 (Rakavy & Shaviv, 1967; Barkat et al., 1967), are now a much more intriguing observational prospect than they were initially. It is not known exactly how rare these spectacular explosions are and also whether or not they can occur in the local Universe (many uncertainties still exist in stellar evolution, especially in mass loss rates). However, the physics robustly gives rise to explosions for stars of high enough mass without tuning. The observation of VMS in the local Universe together with simulations of star formation at lower metallicities strongly suggest that PISNe exist. If such an event has already been detected, it is imperative to bridge the gap between modeling and observations. If the first PISN detection is yet to come, models across the mass range at different metallicities and different envelope properties will be useful in identifying such a detection. For these purposes we have further investigated PISN modeling.

We have simulated the life of four VMS models that span the PISN mass range from the main sequence to explosive death and computed the expected LCs. In the process we have used three different astrophysical codes: GENEC (stellar evolution), FLASH (PISN explosion), and STELLA (LC calculation). The mass loss prescriptions used in GENEC resulted in progenitor models with very little or no hydrogen at the surface and compact radii compared to published very low or zero metallicity models. Mass loss is the dominant uncertainty in VMS models. The models approach the Eddington limit towards the end of their evolution. This mean that mass loss could be more important in real stars, even at (very) low metallicities. On the other hand, magnetic fields may reduce the efficiency of mass loss (Georgy et al., 2017). STELLA’s coupled radiation-hydrodynamics and opacity tables gave relatively fast rise times, helping to alleviate the general discrepancies between PISN model LCs and observations. The versatility of FLASH allowed us to simulate the explosions in different dimensionalities and grid settings.

We find that our 1D explosion properties are in good agreement with those of other groups for similar CO core masses. This agreement with other groups, many of which used different stellar evolution and hydrodynamics codes, strengthens the validity of this work as well as previous works which we have compared to.

We extended our study to 2D and 3D to examine the effects of mixing on the ejecta structure and light curves. This extension was strongly motivated by the prospect of shortening the LC rise time by means of outward mixing of 56Ni (Kozyreva & Blinnikov, 2015). Such an effect would put our PISN models more in line with the PISN candidates among the already detected SLSNe. Focusing on our most massive model, P250 (the only model that meets the SLSN criterion), we achieved convergence in the explosion properties between all three dimensionalities. However, no significant mixing of 56Ni was apparent. Indeed, Figure 15 shows nearly identical LCs between the 1D and the angle-averaged 3D profiles.

Aside from the lack of a strong 56Ni mixing affect, we found that mixing was stronger in 3D than in 2D. The compositional interface that experienced the strongest mixing was the Si-O interface. Additionally, the character of mixing at this interface differs largely between 2D and 3D. In 2D we see the formation of thin RT fingers and in 3D we see larger-scale plumes. Mixing at this interface will likely have spectral consequences since the spatial distribution in the ejecta directly translates to the velocity distribution. This prospect needs to be explored further as RT growth is strongly resolution dependent.

Finally, we computed synthetic LCs of our four models in 1D, and also for the converged 1D and 3D simulations of P250. Our set of four 1D LCs (shown in Figure 13) exhibit very luminous shock break-outs and, at peak, span a range in bolometric luminosity greater than three orders of magnitude, from $10^{41}$ erg s*-1* for P150 to $3\times 10^{44}$ erg s*-1* for P250, according to the amount of radioactive nickel. All exhibit broad plateaus lasting 100–200 days (except P150).

Acknowledgments. MSG and CF are supported by the United States Department of Energy under an Early CAREER Award (DOE grant no. SC0010263). AK was partially supported by the I-Core center of excellence of the CHE-ISF. NY acknowledges support from University of Malaya Research Grant Programme (UMRG : RP020B-14AFR). RH acknowledges support from the World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan and from the “ChETEC” COST Action (CA16117), supported by COST (European Cooperation in Science and Technology). The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n. 306901. The FLASH and STELLA simulations were carried out on the DIRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment is funded by BIS National EInfrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1. DiRAC is part of the National EInfrastructure.