Very Low Energy Supernovae: Light Curves and Spectra of Shock Breakout

TL;DR

This paper uses numerical simulations to study the light curves and spectra of shock breakouts in very low-energy supernovae, revealing their faint, red, and long-lasting nature, and discusses observational prospects.

Contribution

First application of the CASTRO code to shock breakout, extending analysis to very low-energy supernovae with detailed simulations and comparison to analytic models.

Findings

VLE SNe light curves are faint, red, and long-lived.

Shock breakouts are brighter but shorter, aiding detection.

Simulations cover energies from 10^{47} to 10^{50} ergs for red supergiants.

Abstract

The brief transient emitted as a shock wave erupts through the surface of a presupernova star carries information about the stellar radius and explosion energy. Here the CASTRO code, which treats radiation transport using multigroup flux-limited diffusion, is used to simulate the light curves and spectra of shock breakout in very low-energy supernovae (VLE SNe), explosions in giant stars with final kinetic energy much less than 10$^{51}$ erg. VLE SNe light curves, computed here with the KEPLER code, are distinctively faint, red, and long-lived, making them challenging to find with transient surveys. The accompanying shock breakouts are brighter, though briefer, and potentially easier to detect. Previous analytic work provides general guidance, but numerical simulations are challenging due to the range of conditions and lack of equilibration between color and effective temperatures. We consider previous analytic work and extend discussions of color temperature and opacity to the lower energy range explored by these events. Since this is the first application of the CASTRO code to shock breakout, test simulations of normal energy shock breakout of SN1987A are carried out and compared with the literature. A set of breakout light curves and spectra are then calculated for VLE SNe with final kinetic energies in the range $10^{47} - 10^{50}$ ergs for red supergiants with main sequence masses 15 Msun and 25 Msun. The importance of uncertainties in stellar atmosphere model, opacity, and ambient medium is discussed, as are observational prospects with current and forthcoming missions.

Figures

Tables

Table 1. Table 1 SN1987A Breakout Model Results

Model	Star	$v$aaVelocity at breakout. (km/s)	KEbbKinetic energy at breakout. (ergs)	$L_{peak}$ (erg/s)	$\Delta t$ccFull-width half-max of travel-time-corrected light curve. (s)	Max $T_{ef}$ (K)	Max $T_c$ (K)	${\nu }_{peak}$ (Hz/Å)
1.0 B	SN1987A	2.0$\times {10}^4$	1.08$\times {10}^{51}$	6.77$\times {10}^{44}$	70.7	5.13$\times {10}^5$	1.1$\times {10}^6$	6.47e16
2.3 B	SN1987A	3.2$\times {10}^4$	2.01$\times {10}^{51}$	1.53$\times {10}^{45}$	45.0	6.29$\times {10}^5$	$\cdots$	$\cdots$

Table 2. Table 2 Presupernova Star Parameters

Star	Final Mass ($M_{\odot }$​)	He Core Mass ($M_{\odot }$​)	Radius (cm)
SN1987A	15	8	3.2$\times {10}^{13}$
RSG15	12.79	4.27	6$\times {10}^{13}$
RSG25	15.84	8.20	1.07$\times {10}^{14}$

Table 3. Table 3 VLE Breakout Model Results

Model	Star	$v$aaVelocity at breakout. (km/s)	KEbbbKinetic energy at breakout. (ergs)	KEfccFinal kinetic energy of supernova as calculated by KEPLER. (ergs)	$L_p$ (erg/s)	$L_{cr}$ddCorrected for light travel time. (erg/s)	$\Delta t$eeFull-width half-max of travel-time-corrected light curve. (h)	Max $T_{ef}$ (K)	Max $T_c$ (K)	${\lambda }_{max}$ (Å)
A15	RSG15	80	3.86$\times {10}^{46}$	6.58$\times {10}^{46}$	9.50$\times {10}^{39}$	9.57$\times {10}^{39}$	68.4	8.15$\times {10}^3$	$\cdots$	$\cdots$
B15	RSG15	150	6.43$\times {10}^{47}$	1.54$\times {10}^{48}$	3.89$\times {10}^{41}$	3.82$\times {10}^{41}$	35.2	2.06$\times {10}^4$	$\cdots$	$\cdots$
C15	RSG15	320	4.98$\times {10}^{48}$	1.21$\times {10}^{49}$	8.39$\times {10}^{42}$	8.33$\times {10}^{42}$	8.1	4.44$\times {10}^4$	6.40$\times {10}^4$	797
D15	RSG15	650	2.13$\times {10}^{49}$	5.04$\times {10}^{49}$	5.43$\times {10}^{43}$	5.38$\times {10}^{43}$	5.3	7.08$\times {10}^4$	1.15$\times {10}^5$	443
E15	RSG15	1200	5.43$\times {10}^{49}$	1.23$\times {10}^{50}$	2.13$\times {10}^{44}$	2.02$\times {10}^{44}$	3.1	9.97$\times {10}^4$	2.24$\times {10}^5$	228
F15	RSG15	1920	2.16$\times {10}^{50}$	5.07$\times {10}^{50}$	8.25$\times {10}^{44}$	8.07$\times {10}^{44}$	1.83	1.40$\times {10}^5$	2.35$\times {10}^5$	217
G15	RSG15	2400	2.54$\times {10}^{50}$	1.20$\times {10}^{51}$	1.68$\times {10}^{45}$	1.66$\times {10}^{45}$	1.42	1.67$\times {10}^5$	3.18$\times {10}^5$	160
A25	RSG25	150	1.38$\times {10}^{48}$	$\cdots$	9.07$\times {10}^{41}$	9.03$\times {10}^{41}$	67.0	1.57$\times {10}^4$	$\cdots$	$\cdots$
B25	RSG25	179	1.53$\times {10}^{48}$	$\cdots$	1.23$\times {10}^{42}$	1.23$\times {10}^{42}$	57.9	1.70$\times {10}^4$	$\cdots$	$\cdots$
C25	RSG25	210	2.27$\times {10}^{48}$	6.12$\times {10}^{48}$	2.76$\times {10}^{42}$	2.77$\times {10}^{42}$	37.0	2.08$\times {10}^4$	$\cdots$	$\cdots$
D25	RSG25	280	3.52$\times {10}^{48}$	$\cdots$	5.85$\times {10}^{42}$	5.83$\times {10}^{42}$	25.9	2.51$\times {10}^4$	$\cdots$	$\cdots$
E25	RSG25	480	1.10$\times {10}^{49}$	$\cdots$	3.67$\times {10}^{43}$	3.61$\times {10}^{43}$	9.3	3.97$\times {10}^4$	$\cdots$	$\cdots$

Table 4. Table 4 Comparison to Other Models

Model	KEfaaFinal kinetic energy as measured in KEPLER (ergs)	$L_{peak}$ (erg/s)	$L_K$bbPeak luminosity of light curve in KEPLER. (erg/s)	Pred. $L$ (erg/s)	Max $T_{ef}$ (K)	$T_{ef}$ ${}_{,K}$ccPeak effective temperature in KEPLER.(K)	Pred. $T_{ef}$ (K)
A15	6.58$\times {10}^{46}$	9.50$\times {10}^{39}$	5.09$\times {10}^{39}$	4.25$\times {10}^{39}$	8.15$\times {10}^3$	6.93$\times {10}^3$	6.29$\times {10}^3$
B15	1.54$\times {10}^{48}$	3.89$\times {10}^{41}$	3.94$\times {10}^{41}$	3.10$\times {10}^{41}$	2.06$\times {10}^4$	1.99$\times {10}^4$	1.84$\times {10}^4$
C15	1.21$\times {10}^{49}$	8.31$\times {10}^{42}$	6.32$\times {10}^{42}$	5.11$\times {10}^{42}$	4.44$\times {10}^4$	4.02$\times {10}^4$	3.71$\times {10}^4$
D15	5.04$\times {10}^{49}$	5.43$\times {10}^{43}$	3.95$\times {10}^{43}$	3.56$\times {10}^{43}$	7.08$\times {10}^4$	6.36$\times {10}^4$	6.02$\times {10}^4$
E15	1.23$\times {10}^{50}$	2.13$\times {10}^{44}$	1.17$\times {10}^{44}$	1.20$\times {10}^{44}$	9.97$\times {10}^4$	8.34$\times {10}^4$	8.15$\times {10}^4$
F15	5.07$\times {10}^{50}$	8.25$\times {10}^{44}$	6.17$\times {10}^{44}$	8.06$\times {10}^{44}$	1.40$\times {10}^5$	1.30$\times {10}^5$	1.31$\times {10}^5$
G15	1.20$\times {10}^{51}$	1.68$\times {10}^{45}$	1.76$\times {10}^{45}$	2.65$\times {10}^{45}$	1.67$\times {10}^5$	1.69$\times {10}^5$	1.77$\times {10}^5$

Equations

$$\Delta E$$

$${\kappa }_P{}_b$$

$${\kappa }_P{}_p$$

$$\eta$$

$$\frac{\partial E_g}{\partial t}+\vec{\nabla }\cdot \left(\frac{3-f_g}{2}{E}_g\vec{v}\right)-\vec{v}\cdot \vec{\nabla }\left(\frac{1-f_g}{2}{E}_g\right)$$

$$c{\kappa }_g(\alpha T^4-E_g)+\vec{\nabla }\cdot \left(\frac{c\lambda }{{\chi }_R}\vec{\nabla }{E}_g\right)$$

$$+{\int }_g\frac{\partial }{\partial \nu }\left[\left(\frac{1-f}{2}\vec{\nabla }\cdot \vec{v}+\frac{3f-1}{2}\hat{n}\hat{n}:\vec{\nabla }\vec{v}\right)\nu E_{\nu }\right]$$

$$-\frac{3f-1}{2}{E}_g\hat{n}\hat{n}:\vec{\nabla }\vec{v}$$

$$\frac{\partial E_g}{\partial t}+\frac{4}{3}\vec{\nabla }\cdot E_g\vec{v}-\frac{1}{3}\vec{v}\cdot \vec{\nabla }\left(E_g\right)$$

$$c{\kappa }_g(\alpha T^4-E_g)+\vec{\nabla }\cdot \left(\frac{c}{3{\chi }_R}\vec{\nabla }{E}_g\right)$$

$$+\frac{1}{3}{\int }_g\frac{\partial }{\partial \nu }\left[\left(\vec{\nabla }\cdot \vec{v}\right)\nu E_{\nu }\right]$$

$${\tau }_x^{-1}$$

$${\tau }_v^{-1}=\frac{\vec{v}}{3L}$$

$$\frac{{\tau }_v}{{\tau }_x}=\theta =\left(c{\kappa }_P\right)\left(\frac{3c}{{\chi }_R{\vec{v}}^2}\right)$$

$$L_{obs}(t)$$

$$L_{bo}$$

$$T_{obs}$$

Full text

1. Introduction

Supernova shock breakout occurs when the leading edge of the explosion erupts through a star’s surface. As the outgoing shock propagates through the stellar envelope, radiation builds up behind an optically thick leading edge. When the surrounding matter reaches a sufficiently small optical depth, this radiation diffuses out in a relatively short times to produce a bright flash. This flash is the second indication, after the neutrino pulse, that a core-collapse supernova has occurred. Because their properties are determined by the star’s structure in a thin layer near the surface, shock breakouts convey unique information about the progenitor’s surface gravity, radius, composition, and explosion energy.

In a recent survey, Sukhbold et al. (2015) found that neutrino-powered supernovae from 9 to 120 $\mathrm{M}_{\odot}$ ​ with central engines calibrated to reproduce the known properties of the Crab supernova and SN 1987A, fell into two categories: failures, which did not blow the star up at all and made black holes, and robust explosions with energies above 1050 erg. While the treatment of neutrino transport was approximate, this behavior suggests a sensitivity to presupernova properties and an explosion mechanism that is “gated”, i.e. it either works well or not at all. Achieving a neutrino-powered explosion of e.g., $\sim 10^{49}$ erg would require fine tuning. The same study also noted, however, that the observed mass spectrum of black holes is inconsistent with the frequent collapse of the entire star (see also Kochanek 2014a). Some uncertain phenomenon is needed to frequently remove the envelope of failed supernovae, either during black hole formation or before. Given that black holes masses are only determined in binary systems, the envelope may have been frequently lost to a close companion star, but there should also be many cases where detached red supergiants, with the same helium core structure, collapse to black holes. Various possibilities for envelope ejection in these cases have been discussed including the Nadyozhin-Lovegrove effect (Lovegrove & Woosley 2013), envelope ejection via acoustic energy input (Quataert & Shiode 2012; Shiode & Quataert 2014), thermonuclear instability during silicon burning (Woosley & Heger 2015b) pulsational-pair instability supernovae (Woosley et al. 2007; Woosley & Heger 2015a), and the outbursts of luminous blue variables (Smith et al. 2011). These events have in common the ejection of substantial envelope matter, by a shock wave in most cases, with energy much less than 1051 erg. The term “very-low-energy supernova” (VLE SN) is therefore defined in this work as a shock-powered transient from a massive, evolved star with a final kinetic energy significantly less than the $\sim$ 0.6 B expected from core collapse (1 B = $10^{51}$ erg of final kinetic energy).

Gerke et al. (2015) have reviewed the current status of searches for failed supernovae as well as the faint transients arising from the Nadyozhin-Lovegrove effect. Four candidates were initially selected, but followup observations ruled out three of these as they later reappeared. The fourth candidate event satisfied the criteria for a VLE SN and continued to be observed. Recently Adams et al. (2016) reported that the source had dimmed significantly below the progenitor luminosity. Modeling of possible dust effects compared to optical and IR source data are consistent with the event indeed being terminal. The transient’s cool, dim properties cannot be explained by ordinary dust. This event might therefore be the first observed example of the Nadyozhin-Lovegrove effect and an excellent example of a real VLE SN though further observations are needed to make sure the star has indeed disappeared.

VLE SNe are faint and challenging to observe in the general supernova population. Other quantities, like mass, radius, and opacity, being equal, Popov (1993) and Kasen & Woosley (2009) predict that the luminosity of a Type IIp supernova during its plateau stage scales as $E_{\rm exp}^{5/6}$. Thus a supernova with explosion energy, $E_{\rm exp}\sim 10^{48}$ erg would have a luminosity 300 times fainter than a typical Type IIp supernova with $E_{\rm exp}\sim 10^{51}$ erg, or about 1040 erg s*-1*. The duration of the plateau, which goes as $E_{\rm exp}^{-1/6}$ would be 4.6 times longer, or about 400 days.

Shock breakout on the other hand is briefer, brighter, and potentially easier to detect by large field-of-view transient surveys. The breakout from a 5$\times 10^{48}$ erg explosion would have a typical luminosity of 8$\times 10^{42}$ erg s*-1* and last several hours. Unlike breakout in common supernovae, the temperature is relatively low, about 1.15$\times 10^{5}$ K, and a larger fraction of the energy would be emitted longward of the Lyman limit. Observing the shock breakouts of VLE SNe can give occurrence rates for this unusual sort of supernovae as well as constrain the properties of the presupernova star and the explosion energy. As we shall see, the properties of these low energy transients are well defined, but a major uncertainty is the event rate.

Shock breakout has previously been simulated in models for SN1987A (Ensman & Burrows 1992; Tolstov et al. 2013), Type Ib and Ic SNe (Tolstov et al. 2013), Type IIp supernovae (Tominaga et al. 2011), and pair-instability supernovae at cosmic distances (Kasen et al. 2011; Whalen et al. 2013). Although analytic estimates exist (Piro 2013), no numerical simulations have yet been carried out for VLE SNe. In particular, the color temperature, which is distinct from the effective temperature and essential to determining the bolometric correction, has not been calculated because its determination requires a multigroup treatment of radiation. In the coming age of large-scale transient surveys accurate light curves and spectra will be vital for mission planning and analysis. Without such models, these transients might easily be confused with other phenomena having similar time scales and luminosities - for instance, novae, failed Type Ia SNe, tidal disruption events, or common envelope ejection from binary mergers - and the light curve automated selection criteria can inadvertently misclassify interesting transients or even ignore them altogether.

The CASTRO code (Almgren et al. 2010) is a multi-dimensional compressible hydrodynamics code that uses adaptive-mesh refinement. The code incorporates multigroup flux-limited diffusion (MGFLD) transport of radiation (Zhang et al. 2013). Here, CASTRO is used with a constant mesh and only in one dimension. The coordinates are Eulerian and spherical. The diffusion approximation for radiation closes the radiation transport equations by assuming a Fick’s Law diffusion relation between radiation energy density $E_{r}$ and flux $\vec{F}$, $\vec{F}=\vec{\nabla}{E_{r}}$. This approximation is only appropriate in an optically-thick regime, however. When the material becomes optically thin, the diffusion treatment leads to superluminal velocities and produces unphysical behavior. To avoid this, the equations can be modified to incorporate a flux limiter, $\lambda$, that gives the correct limiting behavior in both the diffusive and free-streaming regimes and a stable transition in between (Levermore & Pomraning 1981). Propagation velocities will not exceed the speed of light and the pressure tensor will have the correct limiting behavior. This approximation allows a stable modeling of the crucial transition between optically-thick and optically-thin media that shapes the breakout flash.

Grey or two-temperature (“2T”) radiation transport models radiation as a fluid with a separate temperature $T_{rad}$ that can vary from the gas temperature, but assumes the radiation spectrum is always a blackbody with temperature $T_{rad}$. CASTRO’s multigroup capability, unlike 2T transport, allows the radiation to have a non-thermalized spectrum. In the case of shock breakout, the radiation spectrum is expected to be a dilute blackbody, i.e. it has a blackbody form but peaks at a different wavelength than would be predicted from the radiation energy density. In this paper the temperature computed from a simple $E_{r}=aT^{4}$ relation is called the “effective temperature” $T_{ef}$ ​, while the Wien’s law temperature corresponding to the emitted spectrum’s peak wavelength is the color temperature $T_{c}$ ​. For further details on MGFLD and its specific implementation in CASTRO, see Krumholz et al. (2007) (derivation) and Zhang et al. (2011) (implementation).

Opacity is considered in § 2 and applied to color temperature behavior in VLE SNe in § 3. Results from our SN1987A test problem are first given in § 4. In § 5 we describe the code setup and the starting models for the RSG runs. In § 6 we describe the methods for sampling and post-processing the data. The red supergiant results are described in detail and compared to analytic and numerical predictions in § 7. Prospects for observing are considered in § 8.4.

2. Opacity

The properties of shock breakout are determined by the optically thick-thin transition layer in the star’s atmosphere and modeling it requires careful attention to the opacity. Unfortunately, detailed opacity calculations in the low density, hot regime of shock breakout are not straightforward, especially when bound atoms play a role. Existing studies of shock breakout have often assumed the dominance of electron scattering opacities during shock breakout. This is a reasonable assumption at standard supernova energies where temperatures behind the shock will be in the range 1$\times 10^{5}$ - 1$\times 10^{7}$ K. This assumption begins to break down, however, in the regime probed by low energy events. For VLE SNe breakouts, we expect temperatures in the range 1$\times 10^{4}$ - 5$\times 10^{5}$ K and densities in the range 1$\times 10^{-9}$ - 1$\times 10^{-12}$ g/cm3. At these conditions different opacity processes begin to play a role, and the tabulated opacities used by many supernova and stellar evolution codes do not extend to the low-density regime, requiring the code to do its own opacity calculations.

Nakar & Sari (2010) define a parameter, $\eta$, that predicts whether a shell can thermalize its trapped radiation. This parameter depends strongly on the strength of local absorption. Therefore the absorption opacity, quantified in CASTRO by the Planck mean $\kappa_{P}$ ​, must be considered in detail to ensure accurate color temperature calculations. This section discusses four major opacity processes, their effects on both regular and low-energy shock breakout, and their treatment in these models.

2.1. Free-Free Processes (Compton scattering and bremsstrahlung)

2.1.1 Analytic Representation

At the temperatures, densities, and metallicities explored in this paper, a major contribution to scattering opacity comes from photons colliding with free electrons (Compton scattering), which in these low-energy models can be considered to be in the Thomson limit. This process also contributes to absorption opacity, however, because the collisions are not perfectly elastic and a small amount of energy is exchanged between photon and electron. The energy exchanged per collision is:

[TABLE]

where $T_{g}$ is the gas temperature. If the photon energy $h\nu$ is less than the gas/electron energy $4kT_{g}$ then the photons will gain energy (the inverse Compton process); otherwise they will lose energy to the gas. As the radiation temperature drops, the average photon energy $h\nu$ does as well, and each scattering exchanges less energy. If the only absorption source considered is Compton scattering, then as the breakout energy drops absorption will become inefficient at equilibrating the radiation. The simulation will therefore show the chromosphere retreating within the star to higher temperatures and densities, increasing the $T_{c}$ /$T_{ef}$ ratio. This will be discussed further in § 3.

Bremsstrahlung emission occurs when an electron passes close to another charged particle, generally an atomic nucleus, and the change in its energy causes the emission of a photon. There is an associated inverse bremsstrahlung process wherein a photon strikes an electron moving near a nucleus and causes an increase in its energy. This results in a transfer of energy from radiation to gas and therefore functions as an absorptive opacity. Inverse bremsstrahlung in regions with fully ionized hydrogen and helium follows a Kramer’s Law opacity and can be calculated analytically as:

[TABLE]

where $C_{b}$ is a numerical constant and $l$ is the ionization level of the nucleus.

While Compton scattering depends only on the electron number density $n_{e}$, since it considers only photons scattering off electrons, inverse bremsstrahlung depends on both $n_{e}$ and $n_{i}$, as both an electron and an ion are involved in each interaction. For a completely ionized atmosphere consisting of mostly hydrogen and helium, $n_{i}\approx n_{e}$. This gives $\kappa_{P}\,{}_{b}\propto n_{i}^{2}$, and in fully ionized regions $\rho\approx n_{i}\mu$ where $\mu$ is the mean molecular weight, resulting in a $\rho^{2}$ dependency. Inverse bremsstrahlung therefore drops off more rapidly with $\rho$ than Compton scattering, relevant since shock breakout takes place in a region of sharply decreasing density. The steepness of this density profile is a property of the progenitor star and therefore does not differ between standard CCSNe and VLE SNe.

The gas temperature at the shock front, however, does vary significantly between standard CCSNe and VLE SNe. Compton scattering depends linearly on gas temperature, but inverse bremsstrahlung goes as $\propto T^{-3.5}$. Thus a small drop in temperature causes the inverse bremsstrahlung opacity to increase much faster than Compton scattering, fast enough to overcome the suppression from the $\rho^{2}$ term, making inverse bremsstrahlung much more significant in models at lower temperatures. In the VLE case temperatures during breakout are on the order of 5$\times 10^{4}$ - 3$\times 10^{5}$ K, and in this regime inverse bremsstrahlung plays a significant role and must be considered in simulation. The different contributions of various opacity processes at different energies are illustrated at low energy in Fig. 1 and at a more average energy in Fig. 2.

2.1.2 Implementation

Energy exchange via photon-electron collisions occurs at the same rate as electron scattering. The CASTRO opacity network used for breakout represents absorptive opacity from Compton scattering as a fraction of electron scattering opacity $\kappa_{C}=\epsilon\chi_{e}$, where $\chi_{e}$ is the Thomson electron scattering opacity and $\epsilon$ is some factor $<1.0$. In our simulations $\epsilon$ is generally in the range 1$\times 10^{-3}$ - 1$\times 10^{-4}$.

For inverse bremsstrahlung calculations the CASTRO opacity network implements eq. (2). Note that this approximation becomes invalid if hydrogen and helium are not fully ionized. When the gas temperature drops below approximately 4.5$\times 10^{4}$ K, helium and hydrogen begin to recombine, and this equation gives progressively more unreasonable answers as inverse bremsstrahlung no longer follows a Kramer’s Law form. The opacity network therefore stops calculating bremsstrahlung opacity below this temperature.

2.2. Bound-Free Processes (photoionization)

Photoionization is a complex process. Because cross-sections for ionization depend strongly on the energy of the incoming photon, precise formulas for this opacity are heavily frequency-dependent. However, above 5000 K hydrogen may be assumed to be ionized, and above 4.5$\times 10^{4}$ K helium may also be assumed to be completely ionized. Therefore when $T>4.5\times 10^{4}$ K only metals will contribute to photoionization opacity. In a star of solar metallicity the fraction of the atmosphere that consists of metals is small, but their large cross-sections and low ionization energies make their opacities significant. In this regime the photoionization can be treated as a grey opacity with a Kramer’s Law form.

[TABLE]

where $X_{H}$ is the hydrogen mass fraction and $Z$ the metal mass fraction (not the proton number). However, this raises the question of how to treat photoionization when $T<4.5\times 10^{4}$ K. Here the Kramer’s Law opacity breaks down and truly accurate calculations are quite time-intensive. An approximate grey opacity may be found by doing a number of frequency-dependent calculations and taking their Planck mean, but this is again quite time-intensive. It is therefore worth considering whether photoionization opacities below the helium ionization limit are relevant; in these models the opacity network ceases to calculate bound-free opacities below this temperature.

2.3. Bound-Bound Processes (line opacities)

Line opacities are not accounted for in these simulations due to the complexity and computational cost. Incorporating them would require either generating tables for the appropriate temperature and densities or manually implementing a much larger opacity network that would be frequency-dependent, negating the work done to place other opacities in simpler grey forms and requiring lengthy multigroup simulations for even the bolometric light curves.

The effects of line broadening due to velocity shear in the region of the photosphere are also not included in these calculations, but fortunately the shocks considered here are low-velocity as a rule, and at the explosion energies where electron scattering does not dominate envelope velocities are $\sim 500$ km/s. Thus velocity shear effects can be ignored.

To examine this analytic opacity model, tabulated values from the OPAL database were used for comparison. OPAL tables incorporate the effects of lines, as well as many other opacity processes, to a high degree of precision. OPAL can only give total opacity and not separate absorption and scattering opacities, which these simulations require. Still, it is a useful point of comparison. Fig. 3 presents opacities as calculated by the analytic formulas discussed here compared to calculations using OPAL for a low-energy explosion near shock breakout. The OPAL opacities are systematically higher; this is expected to be the case as OPAL includes more absorption processes. Thus the opacities in these simulations should be taken as lower bounds. This implies that the peak luminosity and color temperature results from these simulations are upper bounds.

3. Color Temperature Behavior in VLE SNe

Shock breakout theory allows for a difference between $T_{ef}$ and $T_{c}$ because the spectral shape of the radiation may be set at a chromosphere depth $R_{c}$ while the overall energy is set at the photosphere depth $R_{p}$ ​, and during breakout these two radii are not necessarily equal. The ratio between $T_{c}$ and $T_{ef}$ is not set directly by the physics of shock breakout. Nakar & Sari (2010) predicts a ratio 1.8 in red supergiants and 2.1 in blue supergiants. Rabinak & Waxman (2011), who consider a low-temperature phase of the supernova akin to conditions during a VLE SN breakout, calculate a ratio in the range 1.1 - 1.8. Many (but not all) numerical simulations of breakout show a ratio of 2-3 between $T_{c}$ and $T_{ef}$ (Ensman & Burrows 1992; Tolstov et al. 2013; Klein & Chevalier 1978). It is reasonable, therefore, to ask what kind of ratio can be expected in the VLE case.

In the following analysis it is important to distinguish between the optical depth as measured from the stellar surface, denoted as $\tau^{\prime}$, as opposed to the optical depth of a single shell, denoted simply $\tau$. In a simulation context a “shell” is assumed to refer to one cell of a 1D spherical model, each cell having a single value for both opacity and density. $\tau$ is therefore equal to $\kappa\rho\Delta r$. In all following discussions of opacity we also assume both electrons and ions to be at the single gas temperature $T_{g}$ tracked by CASTRO.

As Nakar & Sari (2010) describe, if a shell of gas can produce sufficient photons to create a blackbody spectrum at its local gas temperature before the breakout radiation diffuses through it, it will set the radiation’s color temperature to its gas temperature. The emitted peak color temperature therefore reflects the gas temperature of the outermost shell both a) still coupled to the radiation and b) capable of producing sufficient photons to meet this thermalization criterion. This is not necessarily the radius predicted from optical depth arguments. Nakar & Sari (2010) use this criterion to create the parameter $\eta$, which is the ratio of the number of photons a shell must produce in order to set the radiation to its local gas temperature, divided by the number of photons that can be produced in the appropriate time:

[TABLE]

where $t_{s}$ is the time radiation spends in the shell, which is at minimum $\Delta r/c$. While the photons are propagating in the diffusive regime, the time they spend in each shell is $\tau\Delta r/c$. If $\eta<1$, the shell can produce sufficient photons and it will thermalize the radiation. The emitted peak color temperature therefore reflects the gas temperature of the outermost shell still coupled to the radiation and having $\eta<1$, which is here called the chromosphere $R_{c}$. In the limit that radiation spends a full diffusion time in each cell, Nakar & Sari (2010) notes that the condition $\eta=1$ can be alternately expressed as each photon experiencing on average one absorption in the shell, placing the chromosphere at the location where $\tau_{a}\tau\approx 1$.

Shock breakout begins when $\tau^{\prime}\approx c/v_{s}$ and the radiation can escape the star ahead of the shock. Opacity drops steeply in the pre-shock region and once radiation can travel ahead of the shock it will escape from the rest of the star without further interaction. Thus the condition $\tau^{\prime}=c/v_{s}$ is also the condition $\tau=c/v_{s}$ and the corresponding shell is the photosphere of breakout.111Assuming that the medium around the star is transparent to the breakout flash. This is assumed to be true in the models considered here, but in the case of stars with dense CSM/strong stellar winds, it may not be, and the photosphere may be external to the star entirely. The modeling of breakout in stars with complex CSMs is of great interest, but well beyond the scope of this work. If $\tau\approx\tau_{a}$, then $\eta=1$ is $\tau_{a}^{2}\approx 1$ and $R_{c}$ should now be at $\tau_{a}\approx 1$. However, the breakout criterion is now $\tau_{a}\approx c/v_{s}$. As $c/v_{s}>1$, this would place $R_{p}$ at a higher optical depth, i.e. a smaller radius, than $R_{c}$; but $R_{c}$ cannot exceed $R_{p}$ since past $R_{p}$ the radiation no longer spends a diffusion time in each cell and on average is not expected to interact. Thus as long as $\tau_{a}\approx\tau$, $R_{c}\approx R_{p}$ and $T_{c}\approx T_{ef}$.

A typical shock velocity near the surface in a blue supergiant progenitor with a standard-energy ($>$ 0.6 B) shock model is 1.5$\times 10^{4}$ km/s. At this velocity breakout would occur at the radius where $\tau=c/v_{s}=20$. If $\tau_{a}\geq 0.05$ at the same radius, the spectrum will remain in thermal equilibrium and $R_{c}=R_{p}$; otherwise the radiation will drop out of thermal equilibrium before reaching the photosphere and $T_{c}\,>T_{ef}\,$. A typical value for $\tau_{a}$ at the shock front in the SN87A case is of order $10^{-3}$. Note that since the temperature is changing rapidly just behind the shock, $R_{c}$ does not have to be much smaller than $R_{p}$ to give a significant $T_{c}>T_{ef}$.

The primary differences in the VLE SNe case are the significantly lower gas temperatures and shock velocities. Absorptive opacity is higher due to the low-T effects on opacity as discussed in § 2. At the greater optical depth a typical photon experiences many more collisions (essentially $\tau^{2}$) and will thus have more opportunity to thermalize even in the presence of a small absorptive opacity. At $v=1\times 10^{3}$ km/s breakout would occur at $\tau\approx 300$ and $R_{c}=R_{p}$ if $\tau_{a}\geq 0.003$ at that location. This is an easier condition to meet, especially in a regime with significant contributions from inverse bremsstrahlung, and values of $\tau_{a}$ in this region are of order $10^{-2}-10^{-1}$.

To further illustrate the difference with energy, consider a simplified model implementing only electron scattering for both $\chi_{R}$ and $\kappa_{P}$ , with the absorptive component represented by a prescribed ratio $\epsilon=\kappa_{P}\,/\chi_{R}\,$ (see § 2.1). Thus $\tau_{a}=\epsilon\tau$ and the criterion for color temperature becomes $\epsilon\tau^{2}=1$, so the color temperature is set at the depth $\tau=1/\sqrt{\epsilon}$. Then, to achieve $T_{c}$ = $T_{ef}$ , the standard-energy example requires $\epsilon>2.5\times 10^{-3}$ while the VLE example would only require $\epsilon>1\times 10^{-5}$.

The color temperature can therefore be expected to converge with the effective temperature as the shock temperature goes down. This results in a breakout that is brighter in the IR and optical windows than would be expected because the bolometric correction is smaller. The observational impacts of this effect are discussed in § 8.1.

3.1. Velocity Effects

Initial attempts to simulate shock breakout in SN 1987A produced unexpectedly high radiation color temperatures, particularly at the shock front; gas temperature remained steady, but $T_{c}$ increased by an order of magnitude. This discrepancy was ultimately resolved by implementing higher, more accurate absorption opacities. It is worth considering why the simulation results were so sensitive to changes in absorption opacity, as an illustration of the necessity of considering opacity details. Examining the different terms affecting the spectrum in the flux-limited diffusion approximation provides the solution. During breakout, two terms in the radiation-hydro equations dominate the spectrum: an advection term in frequency space and the matter-radiation coupling term. The multi-group treatment of flux-limited diffusion takes the form:

[TABLE]

Consider the optically thick limit where the gas has the most influence on the radiation, in which case $f_{g}=1/3$ and all $\hat{n}\hat{n}$ terms (representing free streaming radiation) disappear.

[TABLE]

On the left-hand side, the first term in the equation represents the quantity to be solved for, the change in group energy with time. The second and third terms represent energy transfer by bulk fluid motion. The three terms on the right-hand side respectively represent energy exchange via absorption/emission, diffusion of radiation, and energy transfer between frequency groups.

This equation can be interpreted as an advection equation in frequency space (see Zhang et al. (2013) for details). The “sound speed” in this medium includes a term dependent on the divergence of velocity, $\nabla\cdot\vec{\vec{v}}$. In practice, this means that as the gas is compressed by the shock, a velocity divergence arises that shifts the spectrum towards the blue. This term can cause the radiation spectral temperature to diverge from the gas temperature even in optically thick regions. In particular, if this term is the only one considered, the spectrum will harden significantly just before breakout as the shock reaches the stellar atmosphere and increases in velocity. The magnitude of this effect is directly linked to the velocity at the shock front - higher velocity leads to greater divergence and more hardening of the spectrum.

Conversely, the absorption-emission source exchange term, which depends on absorptive opacity, will tend to equilibrate the radiation and the matter, so in the case of shock breakout it will cool the radiation. The relative strength of these two terms determines the ultimate color temperature. If the source exchange term dominates, $T_{c}$ will follow the gas temperature. If the advection term dominates, however, $T_{c}$ can be significantly higher. Therefore the relative timescales of these terms must be considered in order to ensure the resulting color temperature has the correct qualitative behavior.

Let the scale length $L=c/\chi_{R}\,\vec{v}$ where $\chi_{R}$ is the scattering opacity. The terms that will alter the spectrum are the frequency group coupling term $(1/3)(\vec{\nabla}\cdot\vec{v})\nu E_{g}$, which depends on velocity divergence, and the absorption/emission energy exchange term $c\kappa_{P}\,(\alpha T^{4}-E_{r})$. The time scale of the exchange term is the inverse of $c\kappa_{P}\,$​. The time scale of the velocity divergence term is the inverse of $\vec{v}/3L$.

[TABLE]

If $\theta>1$, $\tau_{v}>\tau_{x}$, the exchange term dominates, and the radiation spectral temperature will follow the gas. If $\theta<1$, $\tau_{v}<\tau_{x}$, the velocity term dominates, and the radiation spectral temperature will increase above the gas temperature in regions of high velocity divergence, e.g. shock fronts.

For the SN1987A model, $v\sim 15000$ km s*-1* at breakout. If $\kappa_{P}\,/\chi_{R}\,$ = $10^{-6}$, $\theta=0.0012$; if $\kappa_{P}\,/\chi_{R}\,$ increases to $10^{-3}$, then $\theta=1.2$, implying a significant shift in the spectral behavior of the simulation. Fig. 4 shows the temperatures of three SN1987A simulations run up to the moment of breakout in order to test this hypothesis. Since at these temperatures electron scattering dominates both $\kappa_{P}$ and $\chi_{R}$ , the opacity network was simplified to only calculate this process, which allows easy variation by changing the parameter $\epsilon=\kappa_{P}\,/\chi_{R}\,$. Solid lines represent gas temperature and dashed lines represent color temperature. Where the two lines overlap, the radiation and gas are in equilibrium. 222The stair-step appearance of the $T_{c}$ curves is an artifact of the multigroup approximation, where radiation is approximated as a set of groups each corresponding to a range of frequencies. A color temperature calculation done by selecting the group with the highest energy density and applying Wien’s Law to its central frequency will therefore show discrete changes in value corresponding to the boundaries of the frequency groups. These simulations used 64 frequency groups spaced logarithmically in the range 1$\times 10^{15}$ - 1$\times 10^{18}$ Hz.

The red and green dashed lines, representing color temperature in simulations with $\epsilon=10^{-6}$ and $10^{-4}$ respectively, both spike above the gas temperature at the shock front

• the region of highest velocity divergence - and then diffuse forward. By contrast in the $\epsilon=10^{-3}$ simulation (blue) the color temperature follows the gas temperature through the hydrodynamic shock. Thus, the anomalously high color temperatures in the initial simulations, which were run with $\epsilon=10^{-6}$, are due to velocity effects. Changing this ratio to a more realistic $\epsilon\sim 10^{-3}$ brought $T_{c}$ in line with previous simulations, as will be discussed in § 4.

Shock front velocities in the red supergiant models, especially the VLE explosions modeled here, are much slower than in SN 1987A. Velocities at breakout range from 80 to 1200 km s*-1* for RSG15. and 150 to 500 km s*-1* for RSG25. Even for $\kappa_{P}\,/\chi_{R}\,$ = $10^{-6}$ the exchange term dominates in most RSG cases. Velocity effects are therefore unlikely to play a role in VLE SNe breakouts.

In cases with high velocity, low metallicity, or both, this effect may still come into play. Supernovae of more standard energies exploding in very metal poor blue supergiants, as theorized for some models of first-generation stars, might also have unusually high color temperatures. While a full simulated exploration of this space is beyond the scope of this paper, we can qualitatively predict some behavior. Stars with high envelope metallicities and/or low shock velocities will likely show the same 2 - 3 ratio of $T_{c}$ to $T_{ef}$ as SN1987A, but stars with very low envelope metallicities, such as Pop III stars, might exhibit spectral temperatures dramatically higher than predicted. Kasen et al. (2011) estimated the shock velocity at breakout in pair SNe in a red supergiant model at $1.5\times 10^{4}$ km/s, similar to the SN1987A model discussed in this section; in blue supergiants this number can be substantially higher. At the temperatures of a standard-energy breakout inverse bremsstrahlung is highly suppressed and metals contribute most of the photoionization opacity; in a low-metal star total absorption opacity may therefore be quite low. At high velocities and low absorptive opacity the velocity divergence discussed here can easily dominate the exchange term and make the spectrum much hotter than opacity arguments would predict.

4. Verification and Validation of CASTRO: SN1987A

Because CASTRO’s MGFLD module is new, it was important to first verify its capabilities for this sort of simulation. A frequently modeled shock breakout event is that of SN1987A, which has been previously simulated by Ensman & Burrows (1992) (EB92) with the VISPHOT code and Tolstov et al. (2013) (T12) in 2012 with the STELLA code. SN1987A also has upper limits on the temperature and luminosity of its shock breakout that were derived after the fact from observations of the ionization of the surrounding gas (Fransson et al. 1989; Lundqvist & Fransson 1996). The progenitor chosen was an 18 $\mathrm{M}_{\odot}$ blue supergiant with a radius of $3.2\times 10^{12}$ cm, shown in Fig. 5, the same model studied in EB92. This model was also studied in Sukhbold et al. (2015). Two different explosion energies were sampled. One model at 1.0 B was close to the observed value, near $1.3\times 10^{51}$ erg (Arnett et al. 1989) and the other model, at 2.3 B, was chosen to match the second simulation in EB92. The high temperatures in SN1987A’s breakout ensure that electron scattering opacity is dominant.

The 1.0 B explosion was simulated in both single-group (grey radiation) and multigroup mode; the 2.3 B explosion was studied only in single-group. The results are shown in Fig. 6 which gives the bolometric light curve, and Fig. 7 which gives the spectrum of the 1.0 B model at peak $T_{c}$. Peak $T_{c}$ will occur near the very beginning of the breakout, while the bolometric luminosity will rise on a longer timescale; in SN1987A the peak luminosity occurs about 100 seconds after peak color temperature. Bremsstrahlung and photoionization opacities are negligible in this energy regime and are not included. Comptonization is included and the ratio $\kappa_{P}$ /$\chi_{R}$ has been set to $5\times 10^{-3}$ in order to ensure the source exchange term is dominant (§ 3.1). The bolometric light curve is similar to past studies with a luminosity that peaks around $10^{45}$ erg/s, indicating a high-energy explosion. The breakout light curve has a width of $\sim 100$ s, indicating a small stellar radius. Peak effective temperature occurs at peak luminosity, and the ratio of peak temperatures $T_{c}$ /$T_{ef}$ $=1.1\times 10^{6}$ K$/5.1\times 10^{5}$ K = 2.2, matching previous calculations. The 2.3 B explosion has a higher and sharper peak than the 1.0 B. A summary of the results is given in Table 1. CASTRO’s results therefore match both prior simulations and theoretical predictions for the test case of SN1987A.

5. Simulation Setup For the Main Study

5.1. Progenitors

Based on a number of recent surveys (e.g., O’Connor & Ott 2011; Pejcha & Thompson 2015; Sukhbold et al. 2015) the presupernova stars that most commonly produce black holes in stars that still retain their envelopes had main sequence masses in the range 15 to 35 $\mathrm{M}_{\odot}$ ​. Especially prolific in black hole production may be the stars 18 to 26 $\mathrm{M}_{\odot}$ (Horiuchi et al. 2011; Ugliano et al. 2012; Brown & Woosley 2013; Sukhbold & Woosley 2014; Kochanek 2014b; Clausen et al. 2015). In order to sample this range, two progenitor stars are used here, both red supergiants, with main sequence masses 15 $\mathrm{M}_{\odot}$ (RSG15) and 25 $\mathrm{M}_{\odot}$ (RSG25). Their presupernova structures are shown in Fig. 8, Fig. 9, and Table 2. Since the radiation produced by breakout depends chiefly on the shock wave energy, which will be varied, and the presupernova radius, which changes only by about a factor of two in the mass range of interest, these two models should suffice. The models here are taken from Woosley & Heger (2007) and are very similar to those more recently explored by Sukhbold & Woosley (2014) and Sukhbold et al. (2015).

These stars were evolved using the KEPLER code (Weaver et al. 1978; Woosley et al. 2002) from ignition on the main sequence to the time of collapse. An artificial shock of variable energy was initiated by first extracting the iron core and then dumping energy in the bottom 10 zones. Since the calculations here focus on the early light curve and its plateau and do not depend upon nucleosynthesis or details of how the explosion was powered, this simple approach should suffice.

It is important to note that the energy deposited, the shock energy at the time of breakout, and the final kinetic energy at infinity of the ejecta are three different quantities. Not all energy deposited becomes kinetic, and the shock energy at breakout is not the same as the total energy, which includes both positive internal and negative gravitational energy. Both the kinetic energy at breakout, which is relevant for breakout, and the final kinetic energy, which is relevant for the later light curve, are given in Table 3.

In most core-collapse simulations the properties of the pre-collapse stellar atmosphere are irrelevant, since by the time the explosion is seen, the layers that governed breakout have already become transparent. Care must be taken to treat it accurately in shock breakout studies. Additional modifications were therefore made to the progenitor models’ atmospheres. These modifications are discussed in § 5.1.2.

Three calculations were done for each case: a precise breakout simulation using CASTRO (§ 7); a simplified breakout simulation using KEPLER for comparison purposes; and a precise calculation of the later plateau-stage light curve using KEPLER (§ 7.3). KEPLER uses only single temperature flux-limited radiation transport, but as long as a KEPLER progenitor model is mapped into CASTRO at a time when the material is still very optically thick and the radiation is fully thermalized, no loss of precision occurs. Breakout results from KEPLER give information on the bolometric luminosity and effective emission temperature that is useful to compare with CASTRO, and with adequate zoning, give reasonable estimates. These results are discussed further in § 7.3. Only the CASTRO multigroup calculation gives information on the color temperature and spectrum.

Because the energy in the shock at breakout depends only on the local velocity structure in the hydrogen envelope, models for the 25 $\mathrm{M}_{\odot}$  case were generated by simply multiplying the velocities in KEPLER model C25 by a constant factor before inputting them into CASTRO. All RSG15 models were calculated individually using KEPLER.

5.1.1 Choice of Energies

Explosion energies in the range $10^{47}-10^{50}$ ergs were explored for red supergiants. The lower bound on this range is set by the approximate binding energy of the hydrogen envelope. Even the recombination of the hydrogen and helium in the envelope gives this much energy and it is difficult to imagine the ejection of most of the envelope with less energy. This also gives very low velocity and a long transient that might be difficult to detect or be easily confused with other sources. This sort of event might result from nuclear instabilities in the core during oxygen and silicon burning. For example, the stronger silicon flashes studied in stars of near 10 $\mathrm{M}_{\odot}$ by Woosley & Heger (2015b) imparted a kinetic energy to the envelope of $\sim 5\times 10^{49}$ erg, but the weaker flashes imparted far less, down to a few $\times 10^{47}$ erg. Unfortunately those very weak shocks did not make it to the surface before the core collapsed and a new shock was launched. A shock of $3\times 10^{48}$ erg did reach the surface though. Shiode & Quataert (2014) estimate that convection in presupernova stars can drive shocks with $10^{46}-10^{48}$ ergs of energy. Energies of 1047 - 1048 erg were also seen in the failed supernovae studied by Nadezhin (1980) and Lovegrove & Woosley (2013).

“Low mass” pulsational pair instability supernovae (main sequence masses near 70 - 80 $\mathrm{M}_{\odot}$ ) can also give very weak shock waves due to the late onset and weak nature of the pulsations (Woosley 2016, in prep.) During the last few hours of the star’s life, pulsations in the oxygen burning shell drive strong sounds waves through the helium core that steepen into shocks at the base of the hydrogen envelope. By the time the shock reaches the stellar surface, at 1 - 2 $\times 10^{14}$ cm, the shock energy has degraded and ejection speeds are only a few hundred km s*-1* or less, corresponding to kinetic energies $\sim 10^{47}-10^{49}$ erg. Slightly more massive stars give pulses from 10${}^{49}-10^{50}$ erg.

The upper bound to the energy range studied is set by the lowest energy supernovae that have already been detected. The Crab supernova for example, is thought to have had an energy of near 1050 erg (Yang & Chevalier 2015). Although light curves and breakout above $\sim 1\times 10^{50}$ erg have already been studied (Tominaga et al. 2011; Tolstov et al. 2013; Kasen et al. 2011; Whalen et al. 2013; Ensman & Burrows 1992), models F15 (0.5 B) and G15 (1.2 B) are included here to study the transition between standard-energy and low-energy behavior.

5.1.2 Stellar Atmospheres

The KEPLER code is designed to study the internal structure and nucleosynthesis of stars and supernovae and does not treat the outer stellar atmosphere very carefully. However shock breakout is fundamentally an atmospheric phenomenon, since its properties are governed by the critical thick-to-thin transition region near the photosphere. In the SN1987A models, this atmosphere is on the order of $1\times 10^{12}$ cm thick; in the RSG models it is $\sim 6\times 10^{12}$. The effect of variations in the stellar atmosphere was tested by running two different simulations of shock breakout in SN 1987A. The first model had its atmosphere replaced by a power law fit of the form log${}_{10}\,\rho=\alpha_{1}\mathrm{log}_{10}\,r+\beta_{1},\ \mathrm{log}_{10}\,T=\alpha_{2}\mathrm{log}_{10}\,r+\beta_{2}$, where $\alpha,\beta$ were determined by fitting the original KEPLER data. The fits gave $\alpha_{1}=-42.1,\ \beta_{1}=518.0,\ \alpha_{2}=-26.4,\ \beta_{2}=334.5$. The second model then had its atmosphere replaced by the same power law form using $\alpha_{1}/2,\alpha_{2}/2$, generating a significantly less steep gradient. The resulting atmospheric gradients are shown in Fig. 10. As can be seen in Fig. 11, the different atmospheres result in quite different breakout profiles. The two breakouts have the same integrated energy, since that energy is being released from the shock wave itself, but the more slowly varying atmosphere influences the time scale and temperature of the observed transient. A shallower atmospheric gradient creates a wider transition region and thus a breakout that is longer lasting and correspondingly less luminous at peak, while a steeper gradient creates a harder, faster transient.

The true atmosphere of SN 1987A’s progenitor is not expected to vary this much between presupernova models, but the comparison does demonstrate that differences in the atmospheric gradient can produce corresponding and significant differences in the resulting light curves. During KEPLER’s modeling of the propagation of the shock waves from the core to the CASTRO link point, physics changes designed to improve simulation of the shock wave structure can lead to secondary responses in the atmosphere. Since explosions at different energies take different amounts of time to reach the surface, this can also lead to variations between different models in the same pre-supernova progenitor. These divergences are clearly unphysical and modifications to the CASTRO input models are therefore required.

An alternative to a simple power-law extrapolation of the density structure in the KEPLER presupernova star’s outer zones is to use a model stellar atmosphere. For the 15$\mathrm{M}_{\odot}$ red supergiant models, realistic model atmospheres were available from the MARCS database (Gustafsson, B. et al. 2008). The MARCS atmospheres were calculated using a specialized code designed to help observers fit spectral data for red supergiants. The atmospheres include LTE and NLTE effects not simulated in either KEPLER or CASTRO, as well as the effects of line blanketing. The 15 $\mathrm{M}_{\odot}$ progenitor has $T_{ef}$ $\sim$ 3550 K, solar metallicity, and log surface gravity of -0.32. The two MARCS atmospheres closest to the 15 $\mathrm{M}_{\odot}$  progenitor were selected based upon $T_{ef}$ , metallicity, and log specific surface gravity, then interpolated to fit the KEPLER progenitor’s properties. The MARCS atmosphere data extend approximately $5\times 10^{12}$ cm below and $1\times 10^{12}$ cm above the photosphere. In order to ensure a smooth transition a power law is fit to the combined MARCS atmosphere and KEPLER model and used to replace all RSG15 atmospheres. This replacement provided a much more accurate atmospheric gradient. The final progenitor atmospheres are shown in Fig. 12.

The atmosphere of a 25 $\mathrm{M}_{\odot}$ red supergiant is more complex and has been less well-studied in simulation. This progenitor’s envelope is very extended ($\sim 10^{14}$ cm) and loosely bound. Once the star reaches this stage of its life it will have likely lost part or all of the envelope to winds or other instabilities; uncertainties in mass loss make it difficult to reliably estimate the atmosphere’s final structure. The section of parameter space explored by MARCS did not extend near enough to our progenitor star to reliably extrapolate the data. Thus for RSG25 a power law of the same form as the 87A fit described above was fit to the existing KEPLER model and used to replace the progenitor atmospheres.

5.2. Equation of State and Network

The CASTRO simulations used an 18-species network, terminating at nickel, with two auxiliary variables, electron fraction $Y_{e}$ and mean molecular weight $1/\mu$. Nuclear reactions were turned off for the present study and all species were passively advected. In the regions of interest the composition was dominated by hydrogen and helium at relatively low temperatures, so an ideal gas law plus radiation provided an accurate equation of state. KEPLER has its own general equation of state that is good under all conditions except extremely high density. There was good agreement between the KEPLER pressure, internal energy, temperature and density and the equivalent quantities in CASTRO after the remap.

5.3. Opacity

As discussed in § 2, shock breakout is a transition of the shock wave from an optically-thick to optically-thin region, and an accurate treatment of opacity is key to obtaining realistic results. The physics of shock breakout occur at and behind the shock front where the gas temperature is in the range 4.5$\times 10^{4}$ - 1$\times 10^{6}$ K when breakout occurs, depending on the model’s energy. At these temperatures and the relevant mass fractions the primary species of interest are H and He. A Saha solver is used to calculate the ionization fraction of H and He. The electron abundance is then used to calculate both a scattering opacity $\chi_{R}$ and a small contribution to absorptive opacity $\kappa_{P}\,{}_{c}$, computed by assuming a fixed ratio $\epsilon=\kappa_{P}\,{}_{c}/\chi_{R}\,$​. Both Thomson scattering and its contribution to absorptive opacity are grey i.e. insensitive to frequency. Inverse bremsstrahlung absorption is calculated according to the equations discussed in § 2.1. Photoionization absorption contributed by metals is calculated according to the equations discussed in § 2.2. Line opacities, bremsstrahlung processes below 4.5$\times 10^{4}$ K, and photoionization processes below 4.5$\times 10^{4}$ K are not accounted-for, as the former requires much more detailed calculations (§ 2.3) and the grey opacity laws for the latter two processes break down once their assumptions of complete H and He ionization are significantly violated. For details on both the theory and implementation of these opacities, see § 2.

CASTRO is an Eulerian code and thus requires a certain amount of mass in every cell on the grid, or severe instabilities will result. A very thin ambient medium is therefore placed around the progenitor star to keep it stable while the processes of interest run. In the case of shock breakout this ambient medium must be made optically thin to ensure that its presence does not distort the resulting lightcurves. The ambient medium is generally made as cold and thin as possible while still maintaining numerical stability, to avoid affecting simulation results. If the opacity does not have the correct low-temperature behavior, this medium can become opaque and influence the results in an unphysical way. The H-He Saha solver becomes (and stays) correctly transparent in the ambient medium. Photoionization opacity is neglected (see previous section). If it were not, this medium might have some slight absorptive opacity. The medium used in this study had a density of order $10^{-12}$ g cm*-3*. Most of the shock breakout transients themselves, and thus our simulations, are completed by the time the hydrodynamic shock itself reaches the ambient medium, and if the medium is appropriately transparent, it will not interact with the radiation either. Our assumed density for the medium therefore does not affect the breakout transient hydrodynamically in our simulations.

In reality stars near the ends of their lives – especially red supergiants – are prone to shedding large and uncertain amounts of mass. Narrow-line Type II supernovae demonstrate clear signs of CSM interactions. Modeling the full scope of breakout-CSM interactions, however, is beyond the scope of this paper.

6. Light Curves and Spectra

Light curves and spectra were evaluated by sampling the flux in a single distant cell, representing the observer. The measured flux is corrected from the comoving frame back to the lab frame, and then the entire light curve corrected for light travel time. Because of the star’s curvature, a distant observer will not see the breakout front erupting uniformly across the star; rather, more distant portions of the disk will light up at later times since the light must travel slightly farther to reach the observer. The comparatively small (3\times 10^{12}$$\mathrm{cm}) SN1987A progenitor has a light travel time of only 100 seconds, but RSG15 ($\sim$6\times 10^{13}$$\mathrm{cm}) and RSG25 ($\sim$1\times 10^{14}$$\mathrm{cm}) have light travel times of 2000 and nearly 10000 seconds, respectively. The overall effect of this correction is to smear the light curve out in time, increasing the peak width and lowering the peak brightness. We use the same simple light travel correction formula as T12:

[TABLE]

This formula makes three assumptions: that the distance to the observer is large, that the radiation is isotropic, and that the photospheric radius $R_{p}$ remains stationary. The first two assumptions are easily satisfied. The third is less accurate at later times as the envelope begins to expand, but the speed of the photosphere is much smaller than the speed of light and $R_{p}$ can be effectively taken as constant during the light crossing time.

Spectra are calculated by sampling the individual group fluxes in the same cell as the light curve. All multi-group models in this paper were calculated using 64 logarithmically-spaced groups. CASTRO automatically places any energy at frequencies below the the lower limit in the lowest group, and any at frequencies higher than the upper limit in the highest group; thus a failure to resolve the correct spectral range can be detected by checking for anomalously high energies in the lowest or highest group. $T_{c}$ is calculated by selecting the frequency group with the highest energy density and applying Wien’s Law to its central frequency. $T_{ef}$ is calculated from the bolometric radiation flux using the standard $F=\sigma_{SB}T_{ef}\,^{4}$ blackbody relation. We considered breakout to be complete when the bolometric light curve had declined to at least half of peak brightness.

7. Breakout in Red Supergiants

7.1. Bolometric Light Curve

Due to their large radii, breakout transients in red supergiants have a much longer time scale than SN 1987A, further lengthened by the low energies considered. Typical durations of the breakout peak are hours to days. For some of the lower energies, the supernova fails almost completely. For others the envelope is ejected, but not the helium core. The wide range of kinetic energies examined results in a diverse set of peak luminosities. Perhaps unfortunate for their detection, duration and peak brightness are inversely related - that is, the brighter the breakout, the shorter it will likely be.

In RSG15, the final kinetic energies for the VLE models ranged from $6.6\times 10^{46}$ erg to $1.23\times 10^{50}$ erg as listed in Table 3. Their peak luminosities ranged from 9.57$\times 10^{39}$ - 2.02$\times 10^{44}$ erg. The two standard-energy models, F15 (0.5 B) and G15 (1.2 B), had peak luminosities of 8.07$\times 10^{44}$ and 1.66$\times 10^{45}$ respectively. Bolometric light curves for RSG15 are shown in Fig. 13. Full-width half-maximum durations range from 3 - 35 h, with an outlying 70 h for the lowest energy model, A15. Though mass loss leaves RSG15 and RSG25 with similar presupernova masses (12.79 $\mathrm{M}_{\odot}$ vs. 15.84 $\mathrm{M}_{\odot}$ ​), the radius of RSG25 is nearly double the radius of RSG15. Breakouts with similar energies are thus expected to have longer light curves and lower peak energies in RSG25 versus RSG15. Bolometric light curves for RSG25 (Fig. 14) show peak luminosities around $10^{42}$ erg/s and durations in the 25 - 70 h range, excepting E25, whose kinetic energy at breakout was significantly higher.

Models B15, C15, E15, and G15 form a sequence where each model’s final kinetic energy increases by approximately an order of magnitude, spanning $10^{48}-10^{51}$ erg. Their peak luminosities span five orders of magnitude ($10^{41}-10^{45}$ erg/s), with the largest relative increases occurring between B15 and C15, and C15 and E15, about a factor of 20. The duration of B15 (1.54$\times 10^{48}$) is about 35 hours while the durations of C15 (1.12$\times 10^{49}$ erg), E15 (1.23$\times 10^{50}$ erg), and G15 (1.2$\times 10^{51}$ erg) are 8.1h, 3.1h, and 1.4h, respectively, meaning a significant shift in breakout duration occurs between $\sim 10^{48}$ and $10^{49}$ erg of final kinetic energy.

Models D25 and C15 have similar kinetic energies at breakout (3.52$\times 10^{48}$ vs. 4.98$\times 10^{48}$ erg); their durations, however, are significantly different. C15 has a duration of 8.1 hours vs. 25.9 for D25. Models E25 and D15 have similar energies (1.10$\times 10^{49}$ vs. 2.13$\times 10^{49}$ erg) and similar durations of 9.3 vs. 5.3 hours. Again a strong shift in behavior is observed with energy, in this case the influence of progenitor structure on breakout. The much larger radius of RSG25 means that a shock will spend much more time traveling through the envelope than it would in RSG15. This longer travel time and increased interaction with the envelope may amplify the effect of progenitor structure on the resulting breakout. Shock waves energetic enough to traverse the envelope relatively quickly may have less time to interact and therefore show less variation with progenitor structure. Kinetic energy therefore emerges again as the most important parameter governing breakout behavior. This agrees with the analytic equations derived by Piro (2013), discussed further in § 7.4, which give a dependence $L\propto E_{48}^{1.36}$.

Both RSG15 and RSG25 show minor anomalies in their light curves post-peak. RSG15’s shock breakouts show a curious blip in their in decline rate post-peak - initially declining quite steeply, but briefly changing to a shallower slope. Four of RSG15’s breakouts display this anomaly and it is correspondingly extended in the longer-duration light curves. RSG25 shows a similar variation post-peak, although in this case it shows multiple peaks/changes in the rate of decline. None of these small features represent significant variations relative to the light curve of each breakout, but the fact that they appear in all light curves of a single progenitor merits further investigation. It is unclear what physical behavior causes these anomalies, but their appearance seems to be related to the formation of a high-density spike at the photosphere as the hydrodynamic shock begins to expand the envelope into the ambient medium. This spike is likely an artifact of 1D simulation. As noted, none of these features are significant relative to the overall light curves, so they can be neglected.

Only a small fraction of the envelope achieves escape speed in Model A15, as shown in Fig. 15, and the dynamics of this model are not as well determined as the others.

7.2. Color Temperature

CASTRO’s multigroup transport was used to generate the spectra of RSG15 breakouts with various shock energies. These are shown in Fig. 16. All have a dilute blackbody form. Models A15 and B15 were excluded from the multigroup simulation because their large optical depth at breakout posed computational difficulty and the opacity was not realistic anyway (see below). All multigroup models were calculated using 64 logarithmically-spaced groups. Model C15 was run with a frequency range $1.5\times 10^{14}-5\times 10^{16}$ Hz. Models D15 and E15 were run with the range 5.0$\times 10^{14}$ - 1.0$\times 10^{17}$ Hz. Models F15 and G15 were run with range 1.5$\times 10^{14}$ - 1.0$\times 10^{17}$ Hz.

Near peak luminosity in Models C15 - E15 $T_{c}$ ranges from $6.4\times 10^{4}$ to $3.18\times 10^{5}$ K, moving the blackbody emission peak down from the soft X-ray bands into the hard UV. Previous works have found a color temperature higher than the effective temperature by a factor of roughly 2 - 3 (Ensman & Burrows 1992; Tolstov et al. 2013; Klein & Chevalier 1978). Ratios of $T_{c}$ to $T_{ef}$ in Models C15 - G15 in Table 3 range from 1.4 to 2.2, similar to the values suggested by Rabinak & Waxman (2011).

However, as noted in § 3, these previous works considered only high energy explosions where electron scattering grossly dominated the opacity. In the low-energy breakouts here, a larger absorptive component helps to thermalize the diffusing photons. Given the small fraction of absorptive opacity needed for thermalization (§ 3) and the opacities shown in Figs. 1 - 3, we expect the near equality of color and effective temperature at peak for Models A15, B15, and C15. Higher energies, approaching those of SN 1987A, should be dominated by electron scattering with $T_{c}$ close to the maximum values given in Table 3 (Models F15 and G15). Models D15 and E15 should have $T_{c}$ somewhere between the electron scattering value and $T_{ef}$. This is actually good news for the low energy cases since it implies that breakout there can be calculated without recourse to multigroup transport, i.e., by codes using a single temperature.

Further work on these very low energy cases would benefit from using an atomic opacity database such as OPAL, but the implementation of this poses difficulty. OPAL tables give the total opacity and do not distinguish the absorptive and scattering components, as is needed for CASTRO. In the past EB92 overcame a similar limitation by assuming scattering opacity to come entirely from Compton scattering, computing this value, then subtracting it from the OPAL value to get an absorptive opacity. In cases where the scattering opacity is dominant, this technique can lead to errors as one large number is being subtracted from another similar number, giving an uncertain small difference.

7.3. Comparison to KEPLER Results

As noted in Section § 5.1, the model atmospheres were replaced in the move from KEPLER to CASTRO. Fig. 17 shows bolometric luminosity of two breakouts calculated in CASTRO, one using the MARCS atmosphere (blue) and one using a fitted version of the original KEPLER atmosphere (green). The differences are slight but present. KEPLER itself calculates a bolometric light curve that differs from CASTRO, but the divergence appears to be resolution-related, as finer zoning in KEPLER significantly reduced the difference. Full results are shown in Table 4. KEPLER slightly underestimates the peak luminosity of the CASTRO results and falls within a factor of 1 - 3 of the analytic predictions. KEPLER can therefore produce a useful estimate of the bolometric light curve, even with its simpler radiation transport. Of the two codes, CASTRO is the only one that can compute a color temperature, so that result cannot be compared directly. However, if the theoretical work considered in § 2 can be used to estimate the location of the chromosphere, a color temperature may still be calculated from KEPLER results.

7.4. Comparison to Analytic Predictions

Analytic predictions for shock breakout in “normal” Type IIp supernovae already exist in the literature (Tominaga et al. 2011; Matzner & McKee 1999; Katz et al. 2010). Piro (2013) extended these formulas to consider the specific case of low-energy supernovae. These formulas predict the bolometric luminosity and timescale of breakouts based on the properties of the progenitor star and the explosion and can be tested against our numerical results. Piro (2013) gives:

[TABLE]

where $E_{48}=E_{kin}/10^{48},M_{10}=M_{ej}/10\mathrm{M}_{\odot}\,\!,R_{1000}=R_{*}/1000R_{\odot},$ and $\kappa_{0.34}=\kappa/0.34$. Stellar radii are given in Table 2. It is assumed in both cases that the ejecta mass $M_{ej}$ is equal to the size of the hydrogen envelope, that the adiabatic index $\gamma=5/3$, and that the factor $(\rho_{1}/\rho_{*})\sim 1$. The results are shown in Table 4.

There is some ambiguity in the equations as to when the quantity $E_{48}$ should be measured; the kinetic energy at breakout differs from the ejecta’s final kinetic energy at infinity. For RSG15, the peak luminosity predictions fell much closer to the KEPLER and CASTRO results when $E_{48}$ was assumed to be kinetic energy at infinity as measured in KEPLER. The analytic predictions and the KEPLER results are very close and both slightly underestimate the CASTRO luminosities. RSG25 shows much greater variance; the analytic predictions underestimate the numerical results by a factor of 4 - 10, with the inaccuracies increasing with kinetic energy. The analytic formulas therefore give reasonable if not precise estimates for a breakout’s brightness.

8. Observing Prospects

The most realistic and refined simulations are still of little use without observations to test them. This section discusses current and near-future surveys as well as recovered candidates for VLE SNe and CCSNe failures, and provides considerations and guidance for future observations of the transients discussed in this paper, especially with regard to the expected color temperatures in VLE SNe as opposed to standard-energy breakouts.

Fig. 18 shows the evolution of the RSG15 models post-shock breakout as simulated by KEPLER. Fig. 19 shows the same for the RSG25 models. The RSG15 plateau durations vary significantly with energy, as might be expected. Plateau magnitudes tend to be some 1.5 - 2 orders of magnitude lower than the breakout peak.

8.1. VLE SNe Breakout in Optical & IR

The practical upshot of the opacity effects discussed in § 3 is that the color temperatures of these faint breakouts are significantly cooler than those of more energetic events, both because of their lower shock energies and because of the convergence of $T_{c}$  and $T_{ef}$ ​. Although the bolometric luminosity of these events is much lower than normal CCSNe, more of the energy will be emitted at low frequencies, and in an IR band the low-energy breakout may actually appear brighter than its more bolometrically-energetic counterpart. The reported detections of shock breakout by the Kepler satellite (Garnavich et al. 2016), which observes primarily in the optical and near-infrared, emphasize the advantage of this cooler spectrum.

A simple estimate of the brightness of breakout transients in any given band can be found by assuming the spectrum is a blackbody at $T_{c}$  and calculating the fraction of blackbody energy inside that bandpass. More accurate calculations would be done using a measured filter curve. The Kepler satellite has an IR bandpass of 0.4 - 0.9 microns. Fig. 20 shows the results of calculating the fraction of blackbody energy emitted within that bandpass, assuming that $T_{c}$ is equal to the maximum $T_{ef}$  in low-energy events. Color temperature is unlikely to go lower than this value, meaning these curves represent upper bounds. The dynamic range of peak luminosities is significantly compressed, as the brighter breakouts also have higher color temperatures. However the duration of the transient is unaffected, and thus breakout events of similar luminosity can still be distinguished in energy by measuring the duration. The peak luminosities for these filtered curves range from 4$\times 10^{41}$ - 4$\times 10^{39}$ erg/s. These are still dim events, but not out of reach of current and future surveys; notably, they are significantly brighter for longer than a standard-energy breakout in these bands.

Fig. 21 shows simulated KEPLER light curves for the same 0.4 - 0.9 micron range. The KEPLER code cannot directly compute a color temperature, but based on analytic and numerical results we can make some assumptions about its value. The black curve shows the predicted light curve in the 0.4 - 0.9 micron range of an explosion with final KE 2.4 B, calculated by applying formulas from Nakar & Sari (2010) to locate the chromosphere and using the gas temperature at that layer as a color temperature. The red curve shows a similar prediction for a VLE SNe at 1.2$\times 10^{49}$ erg, corresponding to model C15. At this energy color temperature is predicted to converge with effective temperature and a light curve can be made with a 1-T code that assumes $T_{c}$ = $T_{ef}$ . The low-energy supernova, though nearly a thousand times less energetic, peaks at a higher luminosity and stays bright for longer in the IR. Again, these are still dim events, but not out of reach.

The cooler temperatures of RSG breakouts make studying VLE SNe at larger distances a more attractive prospect, as redshifting will move the light into optical and IR windows. But the majority of the supernova’s energy is still emitted at wavelengths shorter than the Lyman $\alpha$ cutoff ($>95$% at a color temperature $1\times 10^{5}$ K), and if an extragalactic supernova is assumed to be embedded in a UV-absorbing ISM within its host galaxy, this energy may be absorbed at the source. Since these transients are already faint, nearby (z $\ll$ 1) events are the primary target, and the majority of absorption will occur either at the source or in the Galaxy. The UV attenuation will therefore depend on the viewing angle through the Galaxy and the unknown circumstellar medium at the source. Nearby supernovae will still offer a better target.

8.2. Searches for Failed Supernovae

Kochanek et al. (2008) proposed to begin a novel search for completely failed CCSNe by looking not for the presence but for the absence of sources. The “Survey About Nothing” monitors 1$\times 10^{6}$ red supergiants with the Large Binocular Telescope looking for the abrupt disappearance of any of these stars. In addition to potentially capturing a core-collapse failure, this survey could also detect VLE SNe coming from one of these sources. The supernovae themselves would be visible as a sudden brightening of the “star” for of order a year, followed by a gradual but complete disappearance. After 4 years of observations Gerke et al. (2015) reviewed the survey data searching for both complete failures as well as the neutrino-mediated transients created by the Nadyozhin-Lovegrove effect. Four candidates were initially recovered, but followup observations ruled out three sources as they later reappeared. The final candidate event satisfied the criteria for a very low energy supernova and continued to be observed. Recently Adams et al. (2016) reported 3 more years of observation on this candidate, showing that the source had dimmed significantly below the progenitor luminosity. Modeling of possible dust effects compared to optical and IR source data suggest that the observed event was terminal and that the transient’s cool, dim properties cannot be explained by dust. This event may therefore be the first observed example of the Nadyozhin-Lovegrove effect and an excellent example of a real VLE SNe. Unfortunately this survey’s cadence is not frequent enough to catch a breakout event.

Reynolds et al. (2015) conducted a search through HST archival data looking for collapse events that were not flagged by survey selection rules at the time and recovered one candidate in the range 25 - 35 $\mathrm{M}_{\odot}$ that may have undergone an optically dark collapse.

8.3. Candidate Shock Breakout Events

Several candidate shock breakout events have been published in the literature, but most are high-energy events suspected to come from compact progenitors. Soft X-ray breakout bursts are easier to detect because of the large number of existing space-based X-ray transient satellites designed specifically for the wide-field coverage and rapid slew time needed to capture breakout. In 2008 Soderberg et al. (2008) serendipitously captured an X-ray transient when a supernova went off during a Swift observation of its host galaxy. Soderberg et al. (2008) attribute this event to a Type Ib/c CCSN breaking out from a dense stellar wind surrounding its progenitor, a scenario consistent with the high mass loss rates of Type Ib/c progenitors near the end of their lives. Unfortunately this rapid high-energy event, while of great interest on its own and as a proof of concept for shock breakout observations, bears little relation to the transients explored in this work.

Closer to the VLE SNe regime, UV observations using the GALEX satellite in 2008 detected two CCSNe very close to the time of explosion (Gezari et al. 2008): one with fading and one with rapidly rising UV emission, suggesting that the latter had been caught during its breakout phase. KEPLER hydrodynamic models combined with the CMFGEN radiation transport code were used to model the observed UV light curve of these events as breakout in a 15 $\mathrm{M}_{\odot}$  red supergiant exploding with a final kinetic energy 1.2 B. The calculated effective temperature was similar to the high-energy RSG15 models presented here. The authors noted an effect that will also be important in the case of VLE SNe, namely that as the bolometric light curve fades the spectral temperature declines, bringing more luminosity into the UV (or optical) observing window. The net effect is an apparent plateau phase that is actually the result of competing processes.

8.3.1 Kepler Satellite Observations

In early 2016 Garnavich et al. (2016) announced the observation of two CCSNe with the Kepler satellite333Not to be confused with the KEPLER stellar evolution code! and provided data suggesting the telescope had also captured the associated shock breakouts. The imaging cadence of Kepler is still insufficient to resolve the breakout itself, but subtracting models of the expected light curve from the data shows a systematic excess consistent with a breakout event producing additional luminosity at the beginning of the transient. While the observed supernovae are much more energetic than the events modeled here, the detections serve as an interesting proof of concept for shock breakout observation, particularly for VLE SNe transients that would have significantly longer timescales.

8.4. Current & Upcoming Observing Programs

The key to capturing these breakouts is high survey cadence, preferably hourly. Even a daily measurement can miss the more energetic breakouts entirely. Followup will be simpler than for compact star breakouts since observers will have a response window measured in hours rather than seconds. Spectroscopic followup is strongly recommended to measure $T_{c}$ ​. The spectrum of a red supergiant breakout will be dominated by the blackbody continuum and hydrogen-helium lines. It may also show absorption features from the surrounding nebula. Unlike the later supernova light curve, the breakout will show little to no nickel or iron emission. The light curve will transition to photometric and spectroscopic behavior typical of a normal (albeit dim) Type II SN. Some planned missions could take good observations of breakout transients. The WFIRST mission would provide wide-field observations in the near-IR and the ULTRASAT program currently in design would launch a rapid-cadence UV satellite that would be ideal for detecting shock breakouts. Among ground-based observatories, Pan-STARRS, PTF/ZTF, LCOGT, and eventually LSST could all make useful observations, although none of these networks are a perfect fit.

9. Conclusions

Supernova shock breakout is a promising tool for detecting otherwise dim supernova and retrieving information about their progenitor stars. Observations - or null detections - of VLE SNe can help understand the full range of CCSNe outcomes and place realistic, observational constraints on the failure and partial failure rate of supernovae. Shock breakouts in VLE SNe, in particular, may be easier to observe than those in regular CCSNe, despite the larger bolometric luminosity of the latter, because of their extended duration and cooler spectral temperature. At the same time they are easier to observe than their associated VLE SNe because of their higher luminosity.

The shock breakout of SN1987A is modeled first as a test of the CASTRO multigroup radiation transport module, giving peak luminosities, durations, and color temperatures similar to other published studies. Careful examination of the spectrum of this model motivates a discussion of opacity, and of velocity terms in the radiation transport equations. The impact of velocity on color temperature in low-opacity high-energy models is illustrated.

Two red supergiants, RSG15 and RSG25, are selected for detailed study. These masses approximately bracket the suspected mass range of failed supernovae, though, in fact, it is the stellar radius and shock energy that matter most. The choice of stellar atmosphere is considered and found to have a potentially significant impact on breakout behavior. A range of low-energy explosions is then modeled in KEPLER and transferred to CASTRO. Two additional events with more standard energies are also simulated in RSG15. All these events give light curves and spectra that show clear variations with both explosion energy and progenitor radius. Peak bolometric luminosities for VLE SNe range from $10^{39}-10^{44}$ erg/s depending on explosion energy and progenitor mass. Peak values are compared to results from the KEPLER code and from analytic predictions and found to be in reasonable agreement.

RSG15 models are then further simulated using CASTRO’s multigroup radiation transport module in order to determine the departure of their spectrum from that calculated assuming radiation and matter are in equilibrium at the photosphere. The relative roles of absorption and scattering processes is found to be critical. Analytic formulas for different absorption processes that become significant at the lower temperatures common in VLE SNe are presented and discussed. Analytic arguments suggest that $T_{c}$ in VLE SNe will behave differently from that in high-energy events. At lower energies, the higher ratio of absorptive to scattering opacity processes will cause $T_{c}$ and $T_{ef}$ to converge, resulting in a cooler spectrum than would otherwise be predicted. From a consideration of sources of absorptive opacity (Figs. 1 -3) and rough estimates of the thermalization criteria (§ 3), we estimate that breakout in RSG with kinetic energies less than about 1049 erg will have nearly equal effective and color temperatures while those close to 1051 erg and above will have values like those determined with the neglect of absorption (Table 3). Intermediate values of explosion energy will have color temperatures between $T_{ef}$ and $T_{c}$ . Further work using more realistic opacities is needed to give accurate quantitative results for these intermediate cases.

VLE SNe breakout is expected to produce a blue ($>$ 1$\times 10^{4}$ K) transient with an approximate duration 3 - 70 h and a bolometric luminosity $10^{39}-10^{44}$ erg. A slim but reasonable fraction of this energy will be emitted at observable optical and IR wavelengths because of the significantly lower color temperatures. Light curves are presented for RSG15 models observed in IR, demonstrating that at these wavelengths shock breakout in VLE SNe can in fact outshine its more energetic counterparts. Current and upcoming observation programs are assessed for suitability in detecting VLE SNe breakout. Cadence is the limiting factor for both existing and future surveys; breakout observations require cadence less than a day and preferably hourly.

10. Acknowledgements

We thank Ann Almgren and John Bell for developing the BOXLIB framework and the CASTRO code as well as Mike Zingale and Max Katz for their substantial contributions to it. Weiqun Zhang provided tireless and invaluable code support throughout this project. EL also thanks Jimmy Fung, Mack Kenamond, and Rob Lowrie for assistance in understanding flux-limited diffusion, and George Fuller for project support. Dan Kasen provided important discussions on the topic of opacity. Elizabeth Baumel provided debugging support; Katie Hamren pointed us towards the SVO Filter Profile Service, which provided clear, helpful, and thorough data on observing filters. This research has been supported by the NASA Theory Program (NNX14AH34G), the DOE High Energy Physics Program (DE-FC02-09ER41438) and a UC Lab Fees Research Award (12-LR-237070).