Gas Dynamics of the Nickel-56 Decay Heating in Pair-instability Supernovae

Very massive 140–260 M⊙ stars can die as highly energetic pair-instability supernovae (PI SNe) with energies of up to 100 times those of core-collapse SNe that can completely destroy the star, leaving no compact remnant behind. These explosions can synthesize 0.1–30 M⊙ of radioactive 56Ni, which can cause them to rebrighten at later times when photons due to 56Ni decay diffuse out of the ejecta. However, heat from the decay of such large masses of 56Ni could also drive important dynamical effects deep in the ejecta that are capable of mixing elements and affecting the observational signatures of these events. We have now investigated the dynamical effect of 56Ni heating on PI SN ejecta with high-resolution two-dimensional hydrodynamic simulations performed with the CASTRO code. We find that expansion of the hot 56Ni bubble forms a shell at the base of the silicon layer of the ejecta ∼200 days after the explosion but that no hydrodynamical instabilities develop that would mix 56Ni with the 28Si/16O-rich ejecta. However, while the dynamical effects of 56Ni heating may be weak they could affect the observational signatures of some PI SNe by diverting decay energy into internal expansion of the ejecta at the expense of rebrightening at later times.


Introduction
Observations by Humphreys & Davidson (1979), Davidson & Humphreys (1997), and Crowther et al. (2010) indicate that stars with masses >100 M e can form in the local universe. Cosmological simulations also suggest that the initial mass function (IMF) of primordial (or Population III) stars is topheavy and that many of them would have died with masses above 100 M e (e.g., Hirano et al. 2014). Very massive stars (VMS; 140-260 M e ) are thought to explode as pair-instability supernovae (PI SNe; Woosley et al. 2002). The original idea of the pair instability was introduced by Barkat et al. (1967) and Rakavy & Shaviv (1967) and further developed by Ober et al. (1983), Glatzel et al. (1985), Stringfellow & Woosley (1988), , and Heger & Woosley (2010). When the core of a VMS evolves to temperatures above 10 9 K, thermal photons in the tail of their Maxwellian distribution (hν1 MeV) begin to freeze out as electron-positron (e − -e + ) pairs through collisions with nuclei. Pair-production occurs at the expense of thermal pressure support of the core against gravity and it begins to contract and become hotter. Rising temperatures and densities in the core eventually ignite explosive oxygen and silicon burning that can completely disrupt the star. PI SNe can produce 10 52 -10 53 erg of energy and 0.1-30 M e of 56 Ni and may be the most energetic thermonuclear explosions in the universe. Unlike core-collapse (CC) SNe, whose central engines are not fully understood, PI SN explosions are much simpler and are emergent features of stellar evolution models. Stellar evolution models suggest that the progenitors of PI SNe can die with a variety of structures, ranging from red supergiants (RSGs) to blue supergiants (BSGs) whose radii differ by a factor of 100. Chatzopoulos et al. (2015) recently found that stellar rotation can even cause PI SN progenitors to lose their hydrogen envelopes and die as bare helium cores. They also found that rotation can shift the lower mass limit of PI SN progenitors from 140 M e down to 85 M e ( 56 Ni production can drop sharply in such events; Chen 2015). It was originally thought that the shells of elements built up in the interior of the star over its life would be severely disrupted and mixed by the passage of the PI SN shock through them prior to breakout from the surface. However, multidimensional models have since shown that the star expands homologously during the explosion, with little mixing between shells (Chen et al. 2011(Chen et al. , 2014Joggerst & Whalen 2011;Gilmer et al. 2017).
PI SN light curves and spectra depend heavily on the structure of the star at death and 56 Ni production and internal mixing during the explosion (e.g., Kasen et al. 2011;Dessart et al. 2013;Whalen et al. 2013aWhalen et al. , 2014aKozyreva et al. 2014;Kozyreva & Blinnikov 2015;Smidt et al. 2015;Mazzali et al. 2019). They are generally characterized by an intense, brief radiation pulse when the shock breaks out of the star followed by a decline in luminosity as the fireball expands and cools. The luminosity can then rise again after a few weeks to months when photons due to the decay of 56 Ni begin to diffuse out of the ejecta. This rebrightening phase can last from weeks to months depending on 56 Ni mass and the structure of the star.
PI SNe could be the ultimate cosmic lighthouses because they can be detected in the near-infrared (NIR) at cosmic dawn at z∼25 by the James Webb Space Telescope and at later epochs by the Nancy Grace Roman Space Telescope and the next generation of extremely large telescopes (Whalen et al. 2013b). They therefore could probe the masses of the first generation of stars (Chatzopoulos & Wheeler 2012;Hummel et al. 2012;Pan et al. 2012;de Souza et al. 2013de Souza et al. , 2014Meiksin & Whalen 2013;Whalen et al. 2013b;Mesler et al. 2014;Hartwig et al. 2018) and the origins of extremely metal-poor stars (Joggerst et al. 2010;Ishigaki et al. 2018;Takahashi et al. 2018). RSG PI SNe have the strongest NIR signals at high redshift while those of BSGs and stripped helium cores are much weaker (Smidt et al. , 2015Whalen et al. 2014b). Several PI SN candidates have now been found (e.g., Gal-Yam et al. 2009;Cooke et al. 2012;Gomez et al. 2019;Nicholl et al. 2020).
A unique aspect of PI SNe are the large masses of 56 Ni they produce, and the effects of its decay on their light curves, spectra, and ejecta dynamics are not fully understood. 5 M e of 56 Ni releases ∼1×10 51 erg (about the energy of a CC SN) as γ-rays. They downscatter in energy as they diffuse out of the ejecta on timescales that depend on the mass between the 56 Ni and the surface of the star, and thus its structure at death. Many light curve and spectrum models take nearly all of these photons to escape and cause the SN to rebrighten at later times when in reality only some do because the remainder heat material deep in the ejecta and create a hot bubble that expands and does PdV work on its surroundings. In principle, the expansion of this bubble could trigger the formation of hydrodynamical instabilities and mix 56 Ni with other elements. Consequently, 56 Ni heating in PI SNe could reduce rebrightening and change the order in which lines appear in the spectra at later times.
These issues have never been resolved because no multidimensional explosion model has ever been run for long enough times to determine the dynamics of the hot bubble (>100 days). We have now evolved PI SNe out to 300 days in two-dimensional (2D) simulations to evaluate the effects 56 Ni heating on the dynamics of the ejecta. Our numerical method and explosion models are described in Section 2. The dynamics of PI SNe at later times (and of the 56 Ni bubble in particular) are examined in Section 3. We discuss the implications 56 Ni heating in Section 4 and conclude in Section 5.

Numerical Method
Our three PI SN progenitor models were first exploded in 1D in KEPLER. The blast profiles were then mapped onto a 2D grid in CASTRO and evolved out to 300 days.

KEPLER Models
Rotating VMS can produce bare He cores, as discussed earlier, while mixing in the cores of even Population III stars can create RSGs with radii 100 times greater than those of the main-sequence star (e.g., Meakin & Arnett 2007;Arnett et al. 2009;Woodward et al. 2014). After the main sequence, carbon and oxygen in the central convective cores of massive stars become closer to the hydrogen-burning shell. If any convective overshooting or other convective boundary mixing occurs, these heavier elements can be mixed into the hydrogen layer, dramatically increasing energy production in the core and expanding the star into an RSG. We therefore chose three models that bracket the likely range of structures of PI SN progenitors: 105 M e and 110 M e He cores (models He105 and He110) and a Z=10 −4 Z e 225 M e RSG that has retained its hydrogen envelope (model U225). These stars explode with energies of 48.3 B, 55.3 B, and 46.6 B (1 B=10 51 erg) and produce 8.53 M e , 13.13 M e , and 16.52 M e of 56 Ni, respectively. The physical properties of all three stars and their explosions are listed in Table 1.
We evolve all three stars from the onset of the PI through core contraction, explosion, and the end of all nuclear burning in the KEPLER stellar evolution code (Weaver et al. 1978). Our simulations were initialized with profiles from Heger & Woosley (2010) and evolved until there was no further change in the energy or composition of the ejecta, ∼10 3 -10 4 s after the explosion. At this time the forward shock in He105 and He110 has just broken through the surface of the star while it has only entered the hydrogen envelope in U225. We perform the explosion in 1D in KEPLER because it naturally emerges from the stellar evolution calculation and we can follow nucleosynthesis in much more detail than would be practical in 2D or 3D.
There are also no significant departures from spherical symmetry in the ejecta over the short times we evolve the explosion in KEPLER. Chen et al. (2011Chen et al. ( , 2014 found that explosive burning did drive the formation of some instabilities in the O shell but they could not grow to large amplitudes because burning ended after a few tens of seconds. They are much smaller than those in CC SNe. We approximate these structures in our 2D CASTRO models by seeding the grid with random velocity perturbations of the order of ∼1% of the local sound speed.
We show elemental abundances at the end of all three KEPLER runs in Figure 1, which show the onion-like structure of elements in the ejecta. No 56 Ni overlaps with the 12 C or 16 O layers. Corresponding density and velocity profiles are shown in Figure 2. The velocity of the forward shock in He105 and He110 is ∼2×10 9 cm s −1 . The helium stars have become unbound and their central densities have dropped to ∼200 g cm −3 . In the U225 model the shock has just entered the hydrogen envelope and its velocity is also ∼2×10 9 cm s −1 . Its central density is ∼1 g cm −3 . The absence of any contact discontinuities in these profiles indicates that no fluid instabilities would have formed in these three explosions by these times in multidimensional simulations. Therefore, if any mixing occurs it is due to the propagation of the shock or the expansion of the 56 Ni bubble at later times.

2D CASTRO Simulations
We then port our 1D KEPLER blast profiles onto a 2D cylindrical grid in CASTRO with the conservative mapping scheme developed by Chen et al. (2013). This method preserves the mass, energy, and momentum of the ejecta on  the new grid over a large range of spatial scales from features in the O and Si shells at ∼10 9 cm to the radius of the star itself, which is four to six magnitudes larger. CASTRO is a multidimensional adaptive mesh refinement (AMR) hydrodynamics code for astrophysical simulations Zhang et al. 2011). CASTRO uses an unsplit piecewise parabolic method (PPM) hydro scheme (Colella & Woodward 1984) with multispecies advection and has several different equations of state (EOS).
We advect the 13 species that constitute the PI SN ejecta: 1 H, 4 He, 12 C, 16 O, 20 Ne, 24 Mg, 28 Si, 36 Ar, 40 Ca, 44 Ti, 48 Cr, 56 Fe, and 56 Ni. We use the Helmholtz EOS (Timmes & Swesty 2000), which includes contributions by both degenerate and nondegenerate relativistic and nonrelativistic electrons, positron-electron pairs, ions, and radiation, during the early phase of the explosion and switch to an ideal gas EOS later when densities fall below 10 −12 g cm −3 . The gravity solver applies the monopole approximation by constructing a spherically symmetric gravitational potential from the radial average of the density and then calculating the corresponding gravitational force everywhere in the AMR hierarchy. It is a reasonable approximation because departures from spherical symmetry in the ejecta are small.
Our CASTRO root grid is 1.875×10 14 cm in r and z with a resolution of 256 2 . We center eight nested grids on the star for a maximum spatial resolution of 256×2 8 , or 2.86×10 9 cm. This zone size was found by Chen et al. (2014) to be sufficient to resolve the onset of fluid instabilities driven by explosive Si and O burning in multidimensional PI SN simulations. Up to eight levels of AMR refinement are also allowed during the initial mapping of the KEPLER profiles into CASTRO and the run thereafter. Refinement is done on gradients in density, velocity, and pressure. We use both static nested grids and AMR because the nested grids ensure that the 56 Ni-rich region always receives the highest resolution while the AMR properly resolves important flows that are outside the most deeply nested grids. Reflecting and outflow boundary conditions are applied to the simulation box on the inner and outer boundaries in r and z, respectively, and we simulate one octant of the explosion.
Because we evolve the explosion for 300 days the SN shock soon reaches the outer boundaries of the original simulation domain. When the shock reaches a grid boundary we quadruple the size of the mesh while holding the number of mesh points constant and conservatively map the flows onto this new grid using the approach by Chen et al. (2013). In our runs the original box is quadrupled three times to a final size of ∼1.2×10 16 cm, which is ∼100 and 10,000 times larger than the radii of the RSG and helium stars, respectively (note that the original mesh fully enclosed the helium stars but not the RSG).
We surround all three stars with a circumstellar medium (CSM) density profile of the form , where ρ s and r s are the surface density and radius of the star. This is done to prevent gas densities from becoming negative after the shock crashes out of the star. The total masses of these envelopes out to the final, largest mesh boundaries are ∼4×10 −5 M e for the helium stars and ∼0.05 M e for the RSG, which are negligible in comparison to the masses of the stars themselves. We chose diffuse envelopes to prevent the formation of reverse shocks that could drive the growth of Rayleigh-Taylor (RT) instabilities at later times as the SN shock plows up the ambient medium. Although the CSM of actual stars like those in our study are not known, the ones we chose are consistent with the weak winds and mild mass loss that are usually associated with lowmetallicity stars (ρ∝r −2 ).

56 Ni
Decay 56 Ni is synthesized during PI SNe at the center of the star by explosive Si burning. It decays with a half-life of 6.1 days to 56 Co, which then decays with a half-life of 77.1 days to 56 Fe. The γ-rays heat the ejecta as they downscatter in energy through it. We do not transport γ-rays in CASTRO so we deposit their energy locally as internal energy. This approximation holds as long as the surrounding ejecta is optically thick to γ-rays. The energy generation rate per unit volume for 56 Ni decay,   Ni 56 , is: Ni (Nadyozhin 1994). The 56 Co mass fraction, X Co 56 (t), is a function of X Ni 56 : , is:  (Nadyozhin 1994).
We plot mass fractions, energy release rates, and total energy released for one M e of 56 Ni in Figure 3. A total of 1.8×10 50 erg is released from its decay over 250 days, 5.92×10 49 erg from 56 Ni   56 Co and 1.28×10 50 erg from 56 Co   56 Fe. Assuming that all the decay energy exits the ejecta as radiation at later times (100 days after explosion), 5 M e  56 Ni is required for a PI SN become a superluminous supernova (SLSN) with a total radiation energy budget 10 51 erg.

Ejecta Dynamics
Gas densities and 56 Ni mass fractions for all three explosions are shown at 20, 100, and 300 days in Figure 4. Since we do not trace the advection of 56 Co and 56 Fe they are included in the 56 Ni (their mass fractions can be extracted from Figure 3). All three stars are completely destroyed, with no compact remnants left behind. Most of the energy from 56 Ni decay has been released by 20 days, and Rayleigh-Taylor (RT) fingers have appeared behind the forward shock in the He105 and He110 models. They are caused by a reflection wave rebounding inward from the surface of the star when the forward shock breaks through it. Their amplitudes remain small and never reach the 56 Ni at the center.
More prominent RT instabilities appear in the U225 explosion that are due to the formation of a reverse shock when the forward shock begins to plow up the H layer of the star. When the shock enters this extended envelope it decelerates and a reverse shock forms and then detaches from it. The two shocks become separated by a contact discontinuity that is prone to RT instabilities if where ρ is the gas density and P is the gas pressure. Since the density of the ejecta decreases with radius, the reverse shock creates a pressure inversion that is in the opposite direction of its density gradient. The contact discontinuity destabilizes and RT fingers appear (such features have been found in earlier simulations of RSG explosions; Chen et al. 2014). They form much closer to the 56 Ni-rich core than in the He105 and He110 runs and have already begun to perturb its outer layers. 56 Co becomes the dominant source of decay energy 60 days after the explosion and small RT fingers appear in the carbon shell at 100 days in He105 and He110 but they do not affect the 56 Ni, as shown in Figure 5. In contrast, the RT instabilities visible at 20 days in the U225 model have now disrupted the outer 56 Ni layers. At 300 days more than 95% of all the decay energy has been released and the 56 Ni bubble has grown to r∼10 16 cm in all three SNe. Heat due to radioactive decay has caused the bubble to expand into the surrounding ejecta, plowing it up into a shell. The thickness of the shell, δr, is ∼10 15 cm and it remains mostly spherical in He105 and He110 but has become heavily disrupted in U225.
We show the 56 Ni, 28 Si, 16 O, and 12 C layers in the He110 and U225 runs at 100 and 200 days in Figure 5. These times capture instabilities and mixing when they are most prominent in these regions. Trace amounts of mixing appear in 12 C in He110 but not in the other three layers, which remain essentially frozen in mass coordinate out to 300 days. The small RT fingers in the outer layer of 12 C form when a weak reverse shock steps back through the ejecta as the SN shock plows up the CSM. It is weak because the CSM is diffuse. Mixing is much more extensive in U225, with severe disruptions of the 12 C, 16 O, and 28 Si shells and the outer layers of the 56 Ni core caused by instabilities driven by the reverse shock.

Mixing Due to 56 Ni Heating
To better quantify mixing in our runs we plot 1D angleaveraged mass fractions for 4 He, 12 C, 16 O,28 Si,and 56 Ni at 20 and 200 days after explosion in Figure 6. At 20 days 56 Ni overlaps with 28 Si in all three SNe because it forms from explosive 28 Si burning, but it also extends slightly into the lower layers of the 16 O shell. In contrast, 56 Ni and 12 C never mix in the helium star explosions. In the U225 run 12 C, 16 O, and 28 Si all appear in the outer layers of the 56 Ni core, indicating that RT instabilities have already caused mixing in the innermost regions of the ejecta by 20 days. The velocity of the forward shock at this time is 1-2×10 9 cm s −1 in all three models. The 4 He at the center of the ejecta is due to photodisintegration of 56 Ni during explosive burning. 56 Ni is formed during the first 20 s of the explosion in 28 Si but then begins to be destroyed at 15-20 s by thermal photons created in the extreme temperatures of rapid burning (see Figure 5 of Chen et al. 2014). Note that the CSM of the two helium stars has H, 4 He, 12 C, and 16 O mass fractions of 0.7, 0.25, 0.03, and 0.02, consistent with the stripping away of the H envelope prior to explosion (these mass fractions were also applied to U225 for consistency).  , and U225 at 20, 100, and 300 days. They remain mostly spherical in He105 and He110 but are disrupted by RT instabilities in U225, which dredge up some 56 Ni by 100 days. A dense shell plowed up by the expansion of the hot 56 Ni bubble is clearly visible in all three models. No mixing due to the expansion of the bubble occurs in any of the models (the mixing in U225 is due to instabilities formed by the reverse shock, not the expansion of 56 Ni).
By 200 days the forward shock has reached r∼1.2×10 16 cm, ∼90% of the energy due to radioactive decay has been released, and elemental abundances have begun to freeze in the mass coordinate. The 56 Ni and 16 O mass fractions still cross each other at 10 −2 , indicating that little 56 Ni has been dredged up in He105 and He110. 56 Ni overlaps with 16 O and 12 C at r∼6-8×10 15 cm in U225, but only small amounts of it appear in 12 C. We show angle-averaged density profiles for all three models at 200 days in Figure 7. The dense shell plowed up by the expansion of the hot 56 Ni bubble within the ejecta is at r∼4-6×10 15 cm. It has a width Δr∼10 15 cm and a density that is three to four times higher than its surroundings. The absence of instabilities in the shell in the helium star runs at this time indicates that 56 Ni heating does not cause much mixing (the disruption of the shell in U225 is due to instabilities driven by the reverse shock, not radioactive decay).
Mixing at these depths in the ejecta could be greater in actual PI SNe than in our models because convection in the core of the star could produce hot spots, off-center ignitions, and asymmetric eruptions that cannot occur in our 1D burn calculations. These can cause prompt 56 Ni mixing and produce chemical yields that are different from those in our simulations. Our models also do not include radiation transport so they cannot capture the decoupling of radiation and gas at later times, which is subject to a variety of radiation-hydrodynamical instabilities.
The total 56 Ni decay energy in our models is 1.58-3.07× 10 51 erg, about 20 times smaller than the explosion energy, so it plays only a minor role in the dynamics of the flows. If all the decay energy is converted into kinetic energy at the base of the 40-50 M e  28 Si shell its net velocity gain would be 2-3× 10 8 cm s −1 , which is consistent with our simulations. We find that this energy is only sufficient to plow up a shell, not accelerate it to velocities capable of driving the formation of fluid instabilities.

Comparison to Previous Models
Joggerst & Whalen (2011) performed the first multidimensional simulations of PI SNe and found little mixing but only evolved them for short times, ending them just after shock breakout. Chen et al. (2011) later studied both core contraction and explosion in 2D PI SN models and saw minor mixing in the oxygen burning shell but they too only evolved the blast for short times. Chatzopoulos & Wheeler (2012) and Chen et al. (2014) performed a suite of 2D PI SN runs with nuclear burning and, like Joggerst & Whalen (2011), found that mixing depends heavily on the structure of the progenitor and is stronger in RSGs than in BSGs. Most recently, Gilmer et al. (2017) carried out the first 3D PI SN simulations with nuclear burning and found minor mixing that was consistent with earlier work. The longest that any of these studies evolved the explosion was less than a month (Chen et al. 2014) so they could not evaluate the effects of 56 Ni decay, but our runs exhibit similar degrees of mixing out to the times these earlier models were run.

Light Curves/Spectra
Mixing deep in the ejecta could alter the spectra of PI SNe at late times by changing the order in which some lines appear. Prominent 28 Si, 24 Mg, and 16 O lines appear at different times with maximum intensities that depend on mixing at the 28 Si/ 16 O interface. However, Jerkstrand et al. (2016) and Chatzopoulos et al. (2019) recently calculated light curves and spectra from multidimensional PI SN models and found that mixing does not affect PI SN light curves and that their spectra and color evolution does not change with viewing angle. However, the degree to which the expansion of the hot 56 Ni bubble diverts energy from rebrightening at later times remains unknown, but it will clearly produce transients that are less luminous than those predicted by most light curve calculations.

PI SNe as SLSNe
If even half of the energy of the 56 Ni decay in our models escapes the ejecta as radiation it could easily power an SLSN with a peak bolometric luminosity that exceeds 10 44 erg for 100 days. If so, how could it be distinguished from magnetarpowered SNe, which have also been proposed as SLSNe? The latter are powered by the spin-down energy of millisecond pulsars with extreme magnetic fields (Blondin et al. 2001;Chevalier & Irwin 2011;Chen et al. 2016Chen et al. , 2017aChen et al. , 2020. Although both PI and magnetar-powered SNe can reproduce SLSN light curves, their central engines can be distinguished by their spectra. Unlike PI SNe, mixing deep in the ejecta of magnetar-powered explosions is rampant and produces rapidly evolving spectra (Chen et al. 2016). In contrast, the decay of large masses of 56 Ni should produce strong 56 Fe lines in the nebular phases of PI SNe.

Rotation
Stellar rotation, although not included in our models, would likely reduce 56 Ni production, decay heating, hot bubble expansion, and therefore mixing deep in the ejecta. As mentioned earlier, rotation can cause VMS to encounter the PI at lower masses because rotational mixing over the life of the star builds up more massive He cores for a given progenitor mass. However, Chen (2015) found that centrifugal forces due to rotation can partially counteract core contraction and decrease explosion energies and 56 Ni synthesis in PI SNe. Consequently, the PI SNe of rotating stars would exhibit less mixing at early times and less radioactive heating and internal expansion of the ejecta at later times.

Mixing in 3D
How might mixing change in our models in 3D? Numerical simulations by Chen et al. (2017b) suggest that mixing is weaker in 3D because of the nature of turbulence. Fluid instabilities are more violent in 2D than in 3D in a given model because turbulent cascades are stronger, so they produce more mixing. Indeed, simulations are often run in 2D before being run at much higher cost in 3D to determine the upper limit to the mixing that would be expected in the problem. Gilmer et al. (2017) found that differences in mixing in 2D and 3D were most pronounced at the Si/O interface in PI SNe. The RT fingers were similar but their overdensities were larger in 3D. However, mixing is sensitive to grid artifacts and initial perturbations so it is not yet clear if mixing in PI SNe is stronger in 2D or 3D.

Conclusion
We find that 56 Ni bubble dynamics does not affect the spectra of PI SNe but can reduce bolometric luminosities during rebrightening. How 56 Ni decay energy is partitioned between the internal and kinetic energies of the ejecta and the radiation that escapes the ejecta in PI SNe remains uncertain. Previous calculations of PI SN light curves often ignored the transformation of decay energy into ejecta dynamics. Kozyreva & Blinnikov (2015) and Kozyreva et al. (2017) calculated 1D PI SN light curves with STELLA, which includes radiation hydrodynamics, so they were able to track deviations from the homologous expansion assumed in some earlier calculations (e.g., Dessart et al. 2013). They found density structures created by 56 Ni bubble expansion that were consistent with those in our models but few effects on rebrightening. However, their models were limited to very low mass resolution (∼150 Lagrangian zones) that may have produced spurious radiation transport. The formation of dense shells due to 56 Ni expansion in our models would clearly come at the expense of rebrightening at later times.
Mixing driven by processes other than 56 Ni expansion can alter the spectra of some PI SNe but not others. For example, RT instabilities driven by the formation of reverse shocks in RSG explosions can jumble together the 12 C, 16 O, and 28 Si shells at early times and affect the order in which their spectral lines emerge at later times when photons due to 56 Ni decay finally diffuse out of the ejecta. Such mixing does not occur in the explosions of bare He cores so there are few if any changes to the order in which spectral lines later appear.
We note that even if some of the energy of radioactive decay is diverted from rebrightening into the internal expansion of the 56 Ni bubble, it does not disqualify PI SNe as SLSN candidates. Only ∼5 M e of 56 Ni is required to produce a superluminous event if most of the energy of decay escapes the ejecta as photons. Our models produce 8.5-16.5 M e of 56 Ni, so even if the majority of the decay energy is lost to work enough could still escape to create an extremely bright transient. We are now developing high-resolution 1D radiation-hydrodynamical simulations of PI SNe with CASTRO with γ-ray transport to determine how energy due to 56 Ni decay is partitioned in PI SNe and its effects on rebrightening.