SHOCK REVIVAL IN CORE-COLLAPSE SUPERNOVAE: A PHASE-DIAGRAM ANALYSIS

The physical mechanism that drives core-collapse supernovae remains an outstanding problem after several decades of research. While there exists clear evidence that massive stars do explode (Smartt 2009), a viable explosion mechanism has not yet been demonstrated robustly by theoretical means (for recent reviews, see Janka 2012; Burrows 2013; Foglizzo et al. 2015).

Since first proposed by Wilson (1985) and Bethe & Wilson (1985), "delayed neutrino heating" has generally been considered the most plausible mechanism for driving core-collapse supernovae (see, however, Kushnir 2015; Papish et al. 2015 for different views). The process envisioned for this mechanism is essentially two-staged. First, the iron core of the massive star collapses to a proto-neutron-star (PNS) that is stabilized by the strong nuclear force, creating an accretion shock that plows through the incoming material that flows inward. This shock stalls, owing to heavy energy losses in neutrino cooling of the shocked accretion layer (and also through dissociation of nuclei as they cross the shock front), but is eventually revived when heating of this layer by neutrinos emitted from the PNS overcomes the energy losses and the inertia of the incoming material. The fundamental issue regarding this mechanism is how (and if!) the competition between neutrino heating and cooling, as well as gravity and ram pressure, can revive the accretion shock and drive it into an outgoing expansion, eventually disrupting the entire envelope of the star.

After three decades of simulations of neutrino heating following core collapse, the overall picture is still a complicated one. There is a broad consensus that self-consistent, one-dimensional simulations generally fail to explode for all but the lowest-mass progenitors (Liebendorfer et al. 2001; Kitaura et al. 2006). This fact has motivated a shift toward two- and subsequently three-dimensional (3D) simulations. In these simulations, multi-dimensional effects, such as turbulence, convection, rotation, and instabilities, come into play and have been demonstrated to generally aid an explosion (see, e.g., Couch 2013; Murphy et al. 2013; Couch & Ott 2014; Fernández 2015; Melson et al. 2015). On the other hand, state-of-the-art 3D simulations have yet to resolve the long-standing problem: they tend to predict sub-energetic explosions (kinetic energies of a few times ${10}^{50}\;\mathrm{erg}$ or less), or no explosion at all (see the discussion in Burrows 2013 and references therein). However, the complexity of the neutrino-driven mechanism (diverse in both the physical processes and the wide range of distances and timescales involved) means that further advances may yet change this conclusion (see, e.g., Lentz et al. 2015 for a more favorable outlook).

This inherent intricacy has also led to an additional class of studies: effective simulations and calculations, which use simplified assumptions and physics. Such works compromise on quantitative accuracy in order to facilitate a qualitative understanding and a more straightforward parameter survey for identifying the underlying principles that are necessary to generate an explosion. A cornerstone of this line of research has been to invoke the neutrino luminosity as a free parameter in the simulations (a "neutrino light-bulb"), instead of generating it self-consistently in the simulations. Unsurprisingly, a fiducial increase of the neutrino luminosity does lead to an explosion. These analyses gave rise to the concept of a "critical neutrino luminosity," which is the minimal luminosity that drives a runaway explosion when all other parameters of the system are predetermined. The critical luminosity was first introduced by Burrows & Goshy (1993), who considered the problem in the stationary approximation, identifying as critical the minimum luminosity for which no steady-state solution exists.

Stationary models of similar nature were applied also in recent works, in search of the origin of the critical nature of the problem, as well as a clearer condition for an explosion (see Yamasaki & Yamada 2005, 2006; Keshet & Balberg 2012; Pejcha & Thompson 2012; Murphy & Dolence 2015 and references therein). For example, under simplifying assumptions such as a neutrinosphere of blackbody temperature ${T}_{\nu }=4{T}_{4}{\rm{}}\;{\rm{MeV}}$ and mass density ${\rho }_{\nu }={10}^{11}{\rho }_{11}{\text{g cm}}^{-3},$ the critical luminosity for which no stationary solution exists is (Keshet & Balberg 2012)

$$\begin{eqnarray}&&{L}_{c}\simeq 6.9\times {10}^{52}{\left[\displaystyle \frac{[{M}_{P,1.3}{T}_{4}{]}^{2}}{[1+0.12\mathrm{ln}\left(\displaystyle \frac{{M}_{P,1.3}^{2}{\rho }_{11}}{{\dot{M}}_{1}{T}_{4}^{3/2}},),{]}^{2}\right)}\right]}^{\mathrm{}}\;{\text{erg s}}^{-1},\end{eqnarray} \tag{ 1 }$$

where ${\dot{M}}_{1}=\dot{M}/(1\;{M}_{\odot }{{\rm{s}}}^{-1})$ and ${M}_{P,1.3}\equiv {M}_{P}/1.3{M}_{\odot }$ are the normalized mass accretion rate and mass of the PNS, respectively.

Treating the neutrino luminosity as a free parameter in simulations has also led to the key observation that the critical luminosity in multi-dimensional simulations is lower than in one-dimensional simulations under similar conditions (Murphy & Burrows 2008; Couch & O'Connor 2014), thus explicitly highlighting the importance of multi-dimensional effects. We note that there are conflicting results regarding whether the critical luminosity in three dimensions is higher than in two dimensions (Bruenn et al. 2009; Nordhaus et al. 2010; Hanke et al. 2012; Dolence et al. 2013; Couch & O'Connor 2014; Takiwaki et al. 2014).

In this work we revisit the critical nature of the transition from a steady accretion to a runaway explosion, with the aid of effective, one-dimensional simulations. While some important multi-dimensional features are necessarily discarded when using this approach, it still qualitatively describes much of the dynamics that dictates the evolution of the accretion shock and the shocked material. By nature, one-dimensional simulations are better suited for parameter surveys, being easier to evaluate as a quantitatively well-defined problem. We are specifically motivated by the fact that in simulations, explosions can occur after the shock goes through a series of increasingly strong oscillations, rather than directly accelerating from a standing shock to runaway expansion (Onishi et al. 2006; Murphy & Burrows 2008; Fernández 2012); naturally, this feature cannot be assessed by purely stationary models (in which the stability of the solution can only crudely be investigated; Nagakura et al. 2013). Our general goal is to examine which quantitative aspects of the flow determine the transition from an oscillating accretion shock to an explosion, as well as their relation to the critical neutrino luminosity. Furthermore, the critical neutrino luminosity for an explosion initiated by oscillations appears to be somewhat lower than that predicted by the stationary models, and we aim to uncover the reason for this trend.

The structure of this paper is as follows. In Section 2 we review the spherically symmetric physical model used in this work. Our simulation code, developed originally for this work, is described in Section 3. Typical results concerning the accretion flow that transitions from oscillations to an explosion are shown in Section 4. Our fundamental observation, considering the conditions for a positive shock acceleration, is presented in Section 5. Here, with the aid of Appendices A and B, we develop a quasi-stationary approximation that allows us to study the shock properties in the phase space of the shock radius and velocity. We compare our conclusions to a frequently suggested timescale criterion for an explosion in Section 6. In Section 7 we summarize our conclusions and discuss some of their implications. The oscillation period is estimated in Appendix C.

Simulations generally indicate that following core bounce, there is a transient phase of order 100 ms, after which the incoming mass accretion rate, ${\dot{M}}_{0},$ and the neutrino luminosity, L, settle on roughly fixed values during the evolution of the shock (Burrows et al. 2007; Marek & Janka 2009). Correspondingly, it is common practice in simplified models of the explosion process to set these two parameters as constants and study the dynamical dependence on their values. Here, we reduce the problem to an idealized spherically symmetric flow, as has been done in many similar studies (see, e.g., Fernández 2012). In spherical symmetry, the equations of motion, used to calculate the dynamics, are the conservation of mass, momentum, and energy:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\ \partial r}\left({r}^{2}\rho u\right)=0,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \left(\rho u\right)}{\partial t}+\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\ \partial r}\left({r}^{2}\rho {u}^{2}\right)=-\displaystyle \frac{\partial p}{\partial r}-\displaystyle \frac{{{GM}}_{P}\rho }{{r}^{2}},\end{eqnarray} \tag{ 3 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{\partial \left(\rho {e}_{\mathrm{tot}}\right)}{\partial t}+\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\ \partial r}\left[{r}^{2}\rho u\left({e}_{\mathrm{tot}}+\displaystyle \frac{p}{\rho }\right)\right]=\rho \dot{q},\end{eqnarray} \tag{ 4 }$$

where u, ρ, and p are, respectively, the fluid radial velocity (in the lab, i.e., the PNS, rest frame), mass density, and pressure,

$$\begin{eqnarray}&&{e}_{\mathrm{tot}}=\displaystyle \frac{1}{2}{u}^{2}+e-\displaystyle \frac{{{GM}}_{P}}{r}\end{eqnarray} \tag{ 5 }$$

is the specific total energy, and e is the specific internal energy. The net neutrino deposition rate (heating minus cooling) per unit fluid mass is denoted by $\dot{q}.$

In the gravitational terms (G being Newton's constant) M_P is the PNS mass. In the general case, M_P should be replaced with the mass enclosed within a radius r at time t, but since in realistic scenarios the mass of the PNS dominates over the mass of the accretion layer, we neglect this layer's self-gravity and use a fixed central mass.

2.1. The Basic Physical Model

Equations (2)–(4) are solved for an accretion layer located between the PNS and the shock radius, ${R}_{S}.$ The immediate downstream of the shock is thus used as an outer boundary. It is common practice in simplified simulations of neutrino-driven explosions to also specify an inner boundary at a radius identified with the PNS surface, ${R}_{P}.$ Physically speaking, the radius ${R}_{P}$ roughly corresponds to the neutrinosphere, above which the neutrino luminosity is approximately constant.

The outer boundary conditions at ${R}_{S}$ are determined by approximating the preshocked material as in a pressure-less free fall with a constant mass accretion rate ${\dot{M}}_{0}=4\pi {R}_{S}^{2}{\rho }_{0}{u}_{0}.$ The velocity (in the lab frame, i.e., the PNS frame) and density of infalling material at the shock are then

$$\begin{eqnarray}&&{u}_{0}=-\alpha {\left(\displaystyle \frac{{{GM}}_{P}}{{R}_{S}}\right)}^{1/2},\quad {\rho }_{0}=\displaystyle \frac{| {\dot{M}}_{0}| }{4\pi \alpha {({{GM}}_{P})}^{1/2}}{R}_{S}^{-3/2},\end{eqnarray} \tag{ 6 }$$

where α is a constant of order unity ($\sqrt{2}$ for perfect free fall). The properties of the material in the immediate downstream of the shock are then determined by the Rankine–Hugoniot jump conditions,

$$\begin{eqnarray}&&{\rho }_{1}{v}_{1}={\rho }_{0}{v}_{0},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{P}_{1}+{\rho }_{1}{v}_{1}^{2}={p}_{0}+{\rho }_{0}{v}_{0}^{2},\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}{v}_{1}^{2}+{e}_{1}+\displaystyle \frac{{P}_{1}}{{\rho }_{1}}=\displaystyle \frac{1}{2}{v}_{0}^{2}+{e}_{0}+\displaystyle \frac{{P}_{0}}{{\rho }_{0}}-{q}_{d}.\end{eqnarray} \tag{ 9 }$$

Here indices 0 and 1 denote quantities upstream and downstream of the shock, respectively, and $v=u-{V}_{S}$ is the fluid velocity relative to the shock, taking into account the (lab frame) shock velocity, ${V}_{S}.$ Hereafter we denote ${\rho }_{1}\equiv {\rho }_{S},$ ${u}_{1}\equiv {u}_{S},$ and ${P}_{1}\equiv {P}_{S},$ indicating the mass density, fluid (lab frame) velocity, and pressure just beneath the shock, respectively.

Finally, q_d in Equation (9) is the energy removed per unit mass through the dissociation of infalling ions by the shock. The value of q_d has several important consequences for the entire profile of the accretion layer. In terms of the boundary conditions, it determines the compression ratio across the shock, $\beta ={\rho }_{1}/{\rho }_{0},$ given by

$$\begin{eqnarray}&&\beta =\displaystyle \frac{\gamma +1}{\gamma -\sqrt{1+2({\gamma }^{2}-1)\theta }},\;\quad \mathrm{where}\quad \theta \equiv \displaystyle \frac{{q}_{d}}{{({u}_{0}-{V}_{S})}^{2}}.\end{eqnarray} \tag{ 10 }$$

Here $\gamma =P/(\rho e)+1$ is the local adiabatic index of the shocked material, which is typically radiation dominated and so $\gamma \sim 1.4$ (Janka 2001). In the limit of zero dissociation $({q}_{d}\to 0,\theta \to 0),$ the compression factor equals that of a simple strong shock, $\beta \to (\gamma +1)/(\gamma -1),$ but if dissociation is significant (compared to the kinetic energy of the infalling material in the shock reference frame), the compression factor will be larger.

2.2. Further Simplifications

As our goal in this work is a qualitative interpretation of the conditions for shock revival through neutrino heating, we apply some further simplifications to our model, which allow for a clearer insight into the numerical simulations. First, we use a simplified equation of state for the shocked material, describing radiation, nonrelativistic nucleons, and relativistic electrons of zero chemical potential and degeneracy (Yamasaki & Yamada 2005; Keshet & Balberg 2012):

$$\begin{eqnarray}&&p=\displaystyle \frac{11}{12}{{aT}}^{4}+\displaystyle \frac{{k}_{{\rm{B}}}\rho T}{{m}_{{\rm{n}}}},\quad e=\displaystyle \frac{11}{4}\displaystyle \frac{{{aT}}^{4}}{\rho }+\displaystyle \frac{3}{2}\displaystyle \frac{{k}_{{\rm{B}}}T}{{m}_{{\rm{n}}}},\end{eqnarray} \tag{ 11 }$$

where a, k_B, and m_n are the radiation constant, Boltzmann constant, and nucleon mass, respectively. While this equation of state neglects the finer aspects of the composition of the shocked material, it allows for more efficient simulations while still capturing the main essence of the problem, especially the transition from a radiation-dominated state $(\gamma \approx 4/3)$ near the shock to a matter-dominated one $(\gamma \approx 5/3)$ near the PNS, as the density increases by several orders of magnitude. Correspondingly, the dissociation energy, q_d, is treated as a free parameter, which can range between the non-dissociation limit q_d = 0 and full dissociation of iron ions into free nucleons with ${q}_{d}={q}_{\mathrm{Fe}}\equiv 8.5\times {10}^{18}\;\mathrm{erg}\;{{\rm{g}}}^{-1}.$ In reality, the actual dissociation across the shock is partial, owing to recombination processes of nucleons into α particles (Fernández & Thompson 2009), so that some of the dissociation energy is later added back to the accreted material. Notice that when dissociation is included, there is a physical upper limit on the radius of the stagnation shock for free-falling material, ${R}_{S}(t=0)\lesssim 200{M}_{P,1.3}{({q}_{d}/{q}_{\mathrm{Fe}})}^{-1}\;\mathrm{km}$.

We also use a simple recipe for neutrino heating and cooling, which is often applied in simplified analytical and numerical models (Janka 2001; Murphy & Burrows 2008). The simple formula, for the total neutrino heating, ${\dot{q}}_{{\rm{He}}}$ and cooling, ${\dot{q}}_{{\rm{C}}},$ are

$$\begin{eqnarray}&&{\dot{q}}_{{\rm{H}}}=1.54\times {10}^{20}{L}_{52}{\left(\displaystyle \frac{r}{100\;\mathrm{km}}\right)}^{-2}{\left(\displaystyle \frac{{T}_{{\nu }_{e}}}{4\;\mathrm{MeV}}\right)}^{2}\;\mathrm{erg}\;{{\rm{g}}}^{-1}\;{{\rm{s}}}^{-1}\end{eqnarray} \tag{ 12 }$$

and

$$\begin{eqnarray}&&{\dot{q}}_{{\rm{C}}}=1.40\times {10}^{20}{\left(\displaystyle \frac{T}{2\;\mathrm{MeV}}\right)}^{6}\;\mathrm{erg}\;{{\rm{g}}}^{-1}\;{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 13 }$$

where L₅₂ is the electron-neutrino luminosity in units of ${10}^{52}\;\mathrm{erg}\;{{\rm{s}}}^{-1},$ ${T}_{{\nu }_{e}}$ is the electron-neutrino temperature at the neutrinosphere, and T is the (radius-dependent) temperature of the matter and photons, assumed to be in equilibrium. In Equation (12) it is assumed that the electron antineutrino luminosity is equal to the electron neutrino luminosity, and that the contribution of the other neutrino types to heating can be neglected. The total energy deposition rate in Equation (4) is then $\dot{q}={\dot{q}}_{{\rm{H}}}-{\dot{q}}_{{\rm{C}}}.$

Finally, since we are only considering matter above the neutrinosphere, where the optical depth is small, we assume a fixed neutrino luminosity above radius ${R}_{P},$ so that L in Equation (12) is independent of radius. In principle, the heating term can be corrected by a factor of ${e}^{-\tau },$ where $\tau (r)$ is the optical depth between radius ${R}_{P}$ and r, but we find that applying this correction has a minor effect on the bulk properties of the accretion layer, which are the focus of our present work. Therefore, we do not apply this correction in the simulations we present in this work.

We solve the flow Equations (2)–(4) with the aforementioned boundary conditions and simplifications using a one-dimensional Lagrangian code. The code implements a standard von Neumann and Richtmyer staggered mesh method (von Neumann & Richtmyer 1950; Richtmyer & Morton 1967) for the equations of motion. The energy equation is solved implicitly with non-adiabatic contributions given by Equations (12)–(13).

The outer boundary is followed by continuously adding Lagrangian cells above the shock in free-fall velocity according to Equation (6). The shock dynamics are calculated with a quadratic artificial viscosity, σ, added to the pressure term of the flow, but this requires some additional caution. Energy losses by dissociation at the shock (when included) can amount to a significant fraction of the internal energy of the post-shock material, and as a result, numerical disturbances may arise within the artificial-viscosity scheme in which the shock is smoothed over a few grid cells. We circumvent this problem by applying the energy loss gradually as the cell passes through the shock region. For every cell entering the shock (with $\sigma \geqslant 0.5P$), the density and internal energy are saved, and the asymptotic post-shock compression factor, β, is calculated according to the density profile of the flow. We regulate the amount of energy lost by a cell, ${\tilde{q}}_{d},$ so that it reaches q_d only after the cell is compressed to the asymptotic value. Quantitatively, ${\tilde{q}}_{d}$ is calculated by

$$\begin{eqnarray}{\tilde{q}}_{d}\left(\rho \right)={q}_{d}\times \left\{\begin{array}{ll}1 & {\rm{}}\;{\rm{if}}\;{\rm{}}\rho \gt {\rho }_{1};\\ \displaystyle \frac{\rho -{\rho }_{0}}{{\rho }_{1}-{\rho }_{0}} & {\rm{}}\;{\rm{if}}\;{\rm{}}{\rho }_{0}\lt \rho \leqslant {\rho }_{1};\\ 0 & {\rm{}}\;{\rm{if}}\;{\rm{}}\rho \leqslant {\rho }_{0},\end{array}\right.\end{eqnarray} \tag{ 14 }$$

where ${\rho }_{0}$ is the pre-shock density of the element and ${\rho }_{1}=\beta {\rho }_{0}$ is the expected density after the shock. The cell then loses up to ${\tilde{q}}_{d}$ energy, but only as long as the internal energy does not drop below its value before entering the shock. This recipe guarantees that the internal energy in a cell cannot drop to negative values through rapid dissociation losses.

In order to ensure stability in the shock downstream, non-adiabatic processes (neutrino cooling and heating) are incorporated only after the cell has lost a total energy of q_d and has compressed to ρ₁. This gradual dissociation loss recipe over a typical shock width (a few cells) is a minor correction since neutrino heating and cooling are generally much weaker near the shock than farther downstream, and we find that applying this recipe changes the critical luminosity by no more than a few percent, while guaranteeing numerical stability.

At the inner boundary, mass elements that enter the PNS (having $u\lt 0$ at R_P) interact with a constant pressure P_P in a ghost cell just below R_P, until they drop completely below R_P (are absorbed in the PNS) or alternatively attain a positive velocity. The latter occurs only when the flow is well into a runaway expansion.

The simulations were typically calculated with about 500–1000 cells in the accretion layer, maintaining a decreasing cell width toward the PNS: the cell widths are adjusted so that each cell is thinner by a factor of ${(1+{\rm{\Delta }})}^{-1}$ with respect to its upper neighbor, with ${\rm{\Delta }}\approx {10}^{-2}\mathrm{to}{10}^{-3}.$ Rezoning is applied when necessary. We find that this resolution allows for numerical convergence, both in terms of the flow near the shock and in the steep density gradient near the PNS.

We initialize a simulation for a given combination of $L,{M}_{P},{R}_{P}$, and ${\dot{M}}_{0}$ by determining a stationary profile, i.e., solving the flow Equations (2)–(4) when all partial time derivatives are set to zero, as is the shock velocity, ${V}_{S}(t=0)=0.$ The outer boundary conditions follow by setting a specific shock radius, ${R}_{S}(t=0),$ so we may solve the entire profile by integrating the stationary flow equations from the shock radius inward to ${R}_{P}.$ The resulting pressure, P_P in the ghost simulation cell just below ${R}_{P},$ is then determined and subsequently serves as a time-independent inner boundary condition at ${R}_{P}.$ As mentioned above, this pressure is assumed to be a property of the (unsimulated) PNS and is thus kept fixed during the simulation. The dynamical evolution is then initiated with some small perturbation, either by intentionally shifting the shock radius (typically by a few hundred meters) or simply by numerical noise. We find that as long as the initial perturbation is small, the evolution that follows does not depend on the specifics.

We mostly varied L and ${\dot{M}}_{0}$ while keeping the PNS mass fixed at ${M}_{P}=1.3{M}_{\odot }$ and its radius at ${R}_{P}=40\;\mathrm{km}.$ We also varied the dissociation parameter (q_d) between the non-dissociation limit q_d = 0 and full dissociation of iron with ${q}_{d}={q}_{\mathrm{Fe}}.$ The choice of the initial radius of the stalled shock ${R}_{S}(t=0)$ warrants some discussion. Without enforcing some additional constraint, this radius has no unique value. In the stationary approximation (Burrows & Goshy 1993; Keshet & Balberg 2012; Pejcha & Thompson 2012), the shock radius is uniquely determined, usually by requiring that the optical depth between the PNS and the shock be equal to 2/3. As shown by Keshet & Balberg (2012), for such an optical depth the shock radius is a slowly varying function of L (when other parameters are kept fixed) and is typically 100–250 km. In the dynamic simulations we follow here, the accretion layer goes through a range of shock radii and optical depths, and so we opt to begin all simulations in a specific parameter survey with the same initial shock radius, typically ${R}_{S}(t=0)=100\;\mathrm{km}.$ In Section 5.4 we vary the initial shock radius, in order to demonstrate that it too affects the critical neutrino luminosity.

One of our goals is to identify the process of shock revival through the growth of oscillations, until the onset of runaway expansion, which—following Fernández (2012)—we define as the onset of an oscillatory explosive flow. Examples of the time-dependent evolution of the shock radius and velocity for ${M}_{P}=1.3\;{M}_{\odot },$ ${R}_{P}=40\;\mathrm{km}$, $| {\dot{M}}_{0}| =0.8\;{M}_{\odot }\ {{\rm{s}}}^{-1}$, and different neutrino luminosities are presented in Figure 1. Dissociation energy losses at the shock were neglected $({q}_{d}=0)$ in these simulations. All these simulations are initiated with the accretion shock set at ${R}_{S}(t=0)=100\;\mathrm{km},$ with a stationary profile that corresponds to this radius.

Zoom In Zoom Out Reset image size

As seen in the figure, most simulations show some initial oscillations of the shock radius and velocity, which have typical periods of tens of milliseconds (see Appendix C for the reason for this narrow range of periods). The distinction between an explosive and a nonexplosive case is clearly evident: when the neutrino luminosity is large enough, the oscillations grow in amplitude until eventually the shock velocity no longer reverts to a negative value, but rather continues to increase with a positive acceleration, culminating in a runaway expansion. The lowest neutrino luminosity that drives such a flow corresponds to the critical luminosity—in this case, about ${L}_{\mathrm{crit},52}\approx 8,$ while for higher luminosities the onset of runaway expansion occurs earlier (after fewer oscillations); for ${L}_{52}=10,$ barely one oscillation is completed before the shock velocity growth becomes exponential (this can be considered a "non-oscillatory mode"; Fernández 2012). The curve for ${L}_{52}=7.5$ represents results typical of $L\lt {L}_{\mathrm{crit}}:$ small oscillations are damped and the flow settles back to stationary accretion.

In Figure 2 we show the results of simulations similar to those of Figure 1, except that full dissociation losses at the shock are included with ${q}_{d}={q}_{\mathrm{Fe}}.$ Evidently, the inclusion of dissociation losses does not qualitatively change the character of the flow; we use this feature in Section 5.1 when we apply the quasi-stationary model. The inclusion of dissociation energy losses does cause a shift in the dependence of the flow on the neutrino luminosity, requiring higher luminosities to achieve a similar flow. This is to be expected, of course, since increased luminosity is required to compensate for the energy-loss rate $| {\dot{M}}_{0}| {q}_{d}\approx 1.4\times {10}^{52}\;\mathrm{erg}\;{{\rm{s}}}^{-1}.$ In practice, the critical luminosity is increased by a larger margin than $| {\dot{M}}_{0}| {q}_{d},$ to ${L}_{\mathrm{crit},52}\approx 13.75,$ since dissociation also changes the compression factor at the shock and hence the outer boundary conditions.

Zoom In Zoom Out Reset image size

In this section we develop the quasi-stationary approximation to derive an approximate expression for the shock acceleration, a_S, given the shock radius, ${R}_{S}$, and velocity, ${V}_{S}.$ This derivation allows us to identify regions in the ${R}_{S}-{V}_{S}$ plane in which both the velocity and acceleration remain invariably positive, hence indicating in a phase-space fashion when the conditions are favorable for a runaway expansion.

First, consider the spherical (with respect to the center of the PNS) moment of inertia of the accretion layer between ${R}_{P}$ and ${R}_{S},$

$$\begin{eqnarray}&&I(t)\equiv {\displaystyle \int }_{{R}_{P}}^{{R}_{S}}{r}^{2}{dm}={\displaystyle \int }_{{R}_{P}}^{{R}_{S}}4\pi {r}^{4}\rho \;{dr}.\end{eqnarray} \tag{ 15 }$$

Its first time derivative is

$$\begin{eqnarray}\displaystyle \frac{{dI}}{{dt}} & = & {\displaystyle \int }_{{R}_{P}}^{{R}_{S}}4\pi {r}^{4}\displaystyle \frac{\partial \rho }{\partial t}{dr}+{\left[4\pi {r}^{4}\rho \displaystyle \frac{{dR}}{{dt}}\right]}_{P}^{S}\\ & = & {\displaystyle \int }_{{R}_{P}}^{{R}_{S}}4\pi {r}^{4}\displaystyle \frac{\partial \rho }{\partial t}{dr}+4\pi {R}_{S}^{4}{\rho }_{S}{V}_{S}.\end{eqnarray} \tag{ 16 }$$

The notation ${[\cdots ]}_{P}^{S}$ stands for the difference between the expression in the immediate downstream of the shock and just outside the PNS, where the latter term is null in our model since we fix ${{dR}}_{P}/{dt}=0.$

We find that a general feature in the simulated flows is that to a very good approximation, the local mass accretion rate is uniform in the accretion layer,

$$\begin{eqnarray}&&\dot{M}(r)\simeq \dot{M}\equiv 4\pi {r}^{2}\rho u,\end{eqnarray} \tag{ 17 }$$

corresponding to slow changes in the density profile. This reflects the highly subsonic nature of the downstream flow. Note that $\dot{M}$ is not identical to the incoming mass accretion rate, ${\dot{M}}_{0},$ since it is modified by the shock velocity. Indeed, the mass accretion rate just below the shock is then

$$\begin{eqnarray}&&{\dot{M}}_{S}={\dot{M}}_{0}+4\pi {R}_{S}^{2}{V}_{S}{\rho }_{0}(\beta -1).\end{eqnarray} \tag{ 18 }$$

For a radius-independent mass accretion rate, the first term on the right-hand side of Equation (16) is approximately zero, and so

$$\begin{eqnarray}&&\displaystyle \frac{{dI}}{{dt}}\simeq 4\pi {R}_{S}^{4}{\rho }_{S}{V}_{S}.\end{eqnarray} \tag{ 19 }$$

Given that dI/dt is now approximated as a function of R_S and V_S, the second time derivative of the moment of inertia can be used to create an implicit relation for the shock acceleration,

$$\begin{eqnarray}&&\left(\displaystyle \frac{{d}^{2}I}{{{dt}}^{2}}\right)\simeq {V}_{S}\displaystyle \frac{\partial }{\partial {R}_{S}}\left(\displaystyle \frac{{dI}}{{dt}}\right)+{a}_{S}\displaystyle \frac{\partial }{\partial {V}_{S}}\left(\displaystyle \frac{{dI}}{{dt}}\right).\end{eqnarray} \tag{ 20 }$$

We emphasize that this is only an approximate equality due to the assumption of a uniform accretion rate in the entire accretion layer, resulting in dI/dt being a function of ${R}_{S}$ and ${V}_{S}$ alone.

Equations (19)–(20) relate bulk properties of the accretion layer, dI/dt and ${d}^{2}I/{{dt}}^{2},$ to the quantities at the accretion shock. We find that this relationship is reproduced to a very high accuracy in our simulations. In principle, this relation can be used to calculate the shock acceleration for a given profile, allowing us to predict the oscillatory movement of the shock, and whether or not runaway is to be expected. Generally, however, ${d}^{2}I/{{dt}}^{2}$ cannot be derived without a dynamical simulation, for which a prediction of a_S is redundant. Nevertheless, using some dedicated approximations, we derive in the following a quasi-stationary model, which does allow for the prediction of a_S, mapping it as a function of the instantaneous shock and flow properties.

5.1. A Quasi-stationary Approximation

We modify the simple, stationary approximation by including a non-zero shock velocity when setting up the conditions at the shock and then solving Equations (7)–(9) as described above. This approximation, which we refer to as "quasi-stationary," is applicable when the shock velocity is small with respect to the velocities in the accretion layer, which in turn are smaller than the typical sound speed in the flow. Note that the assumption of a uniform accretion rate implies that the velocity profile in the accretion layer adjusts quickly to changes in the shock radius, while a subsonic flow (also assumed in the stationary approximation) ensures that thermodynamic changes in the shock quickly advect throughout the shocked material and influence the inner regions of the flow.

Quantitatively, these conditions can be assessed through the typical timescales in the flow, defined by the oscillation period, t_osc,

$$\begin{eqnarray}&&{t}_{\mathrm{sc}}\equiv {\displaystyle \int }_{{R}_{P}}^{{R}_{S}}\displaystyle \frac{{dr}}{| {c}_{s}| }\end{eqnarray} \tag{ 21 }$$

and

$$\begin{eqnarray}&&{t}_{\mathrm{adv}}\equiv {\displaystyle \int }_{{R}_{P}}^{{R}_{S}}\displaystyle \frac{{dr}}{| u| }\approx \displaystyle \frac{{M}_{\mathrm{env}}}{\dot{| M| }},\end{eqnarray} \tag{ 22 }$$

which are the sound crossing time and advection time, respectively. In the second equality for the advection time, M_env is the mass of the envelope, i.e., the accretion layer. The second equality for t_adv is exact only if the mass accretion rate is indeed uniform in the entire layer. To be precise, t_osc is relevant when the flow goes through several oscillations before growing to runaway expansion (or damping out). Once the flow evolves to runaway expansion—or alternatively, if the expansion is initially non-oscillatory—the appropriate measure of the shock evolution should be its dynamical time, ${t}_{S}\equiv {R}_{S}/{V}_{S}.$

In the parameter range of interest the post, shock flow is very subsonic, with t_sc being a few milliseconds and the shortest of the three timescales. The competition between the oscillation and advection timescales is more complicated, since the latter is very sensitive to the compression factor at the shock and hence to dissociation losses. When neglecting dissociation losses, we find that for luminosities close to or exceeding L_crit, the advection timescale is ${t}_{\mathrm{adv}}\sim 10\;\mathrm{ms},$ while the oscillation period is about 50–60 ms. This typical value for the oscillation period can be quantitatively assessed with our quasi-stationary model; see Appendix C.

The hierarchy ${t}_{\mathrm{sc}}\lt {t}_{\mathrm{adv}}\lt {t}_{\mathrm{osc}}$ implies that the quasi-stationary approximation should be applicable. We demonstrate this explicitly in Figure 3, which compares four snapshots of the velocity, specific energy, and mass density profiles in the accretion layer, for the case $| {\dot{M}}_{0}| =0.8{M}_{\odot }\ {{\rm{s}}}^{-1},$ ${L}_{52}=8$ (the solid blue curve in Figure 1), as found in the simulation and in the quasi-stationary approximation (using the instantaneous shock radius and velocity from the simulation). We also show the stationary profiles (found for the instantaneous shock radius but setting ${V}_{S}=0$) at every snapshot. Clearly, the quasi-stationary approximation offers a significant improvement over the stationary case, offering an almost exact fit not only in specific energy and density but also in the velocity profile of the flow. Correspondingly, the quasi-stationary approximation we use below to analyze the shock acceleration is suitable for the regime of small dissociation losses—or gradual dissociation in the context of the microphysics of the shocked material.

Zoom In Zoom Out Reset image size

The quasi-stationary approximation is not quite as successful in the opposite limit of full dissociation into free nucleons across the shock. For ${q}_{d}={q}_{\mathrm{Fe}},$ the compression factor at the shock is higher, generally between 10 and 20, which results in a low fluid velocity; for luminosities close to critical we find that t_adv can reach 20 ms. On the other hand, t_osc is relatively insensitive to the level of dissociation, and so the quantitative deviations of the actual shock acceleration from that predicted by the quasi-stationary model become larger. Nonetheless, we find that the qualitative behavior of the flow and runaway expansion is similar even in the limit of full dissociation (since the oscillation period remains the largest timescale in the problem), and so we do assess that the quasi-stationary approximation is generally a useful starting point in the qualitative analysis of the transition to runaway expansion.

The advantage of the quasi-stationary approximation is that it allows us to estimate the second time derivative of the moment of inertia, ${({d}^{2}I/{{dt}}^{2})}_{\mathrm{QS}},$ which has units of energy. The full derivation is given in Appendix A; here we quote the final result,

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}\displaystyle \frac{{d}^{2}I}{{{dt}}^{2}}\simeq {{\mathbb{E}}}_{\mathrm{QS}}+{W}_{\mathrm{PNS}}+{W}_{B},\end{eqnarray} \tag{ 23 }$$

where the three components on the right-hand side are an effective energy, ${{\mathbb{E}}}_{\mathrm{QS}},$ a work term associated with the PNS, W_PNS, and the energy advected across both boundaries, W_B. These are defined as follows.

The effective energy, composed of kinetic, gravitational, and internal contributions, can be written in the form

$$\begin{eqnarray}&&{{\mathbb{E}}}_{\mathrm{QS}}\approx \tilde{K}+\tilde{{\rm{\Omega }}}+\tilde{U}+{\tilde{B}}_{S}\approx \tilde{K}+\tilde{{\rm{\Omega }}}+\tilde{U},\end{eqnarray} \tag{ 24 }$$

where the last expression can be estimated from the quasi-stationary model. Here we define

$$\begin{eqnarray}&&\tilde{K}\equiv {\displaystyle \int }_{{R}_{P}}^{{R}_{S}}\left(\displaystyle \frac{3-\gamma }{\gamma }\right)\displaystyle \frac{1}{2}{u}^{2}{dm},\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&\tilde{{\rm{\Omega }}}\equiv {\displaystyle \int }_{{R}_{P}}^{{R}_{S}}\left(\displaystyle \frac{3-2\gamma }{\gamma }\right)\left(-\displaystyle \frac{{{GM}}_{P}}{r}\right){dm},\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\tilde{U}\equiv {\displaystyle \int }_{{R}_{P}}^{{R}_{S}}\left(3\displaystyle \frac{\gamma -1}{\gamma }\right)\left[\displaystyle \frac{1}{| \dot{M}| }Q(r)\right]{dm},\end{eqnarray} \tag{ 27 }$$

and

$$\begin{eqnarray}&&{\tilde{B}}_{S}\equiv {\displaystyle \int }_{{R}_{p}}^{{R}_{s}}\left(3\displaystyle \frac{\gamma -1}{\gamma }\right){b}_{S}(m)\;{dm},\end{eqnarray} \tag{ 28 }$$

where γ is the local adiabatic index, Q(r) is the rate of change of the total internal energy due to heating and cooling in a layer extending between r and the shock,

$$\begin{eqnarray}&&Q(r)\equiv {\displaystyle \int }_{r}^{{R}_{s}}\dot{q}\;{dm}={\displaystyle \int }_{r}^{{R}_{s}}4\pi \rho {r}^{2}\dot{q}\;{dr},\end{eqnarray} \tag{ 29 }$$

and b_S is the Bernoulli function (specific energy),

$$\begin{eqnarray}&&b\equiv \displaystyle \frac{1}{2}{u}^{2}+e+\displaystyle \frac{p}{\rho }-\displaystyle \frac{{{GM}}_{p}}{r},\end{eqnarray} \tag{ 30 }$$

evaluated for every mass element according to its value when crossing the shock. The non-adiabatic (heating and cooling) contribution to the energy is contained in the double integral $\tilde{U},$ which tracks the total gained non-adiabatic energy in the accretion layer.

The non-standard coefficients of the energy integrands in Equations (25)–(27) arise from writing ${{\mathbb{E}}}_{\mathrm{QS}}$ in a form that is susceptible to the quasi-stationary approximation. In this approximation, terms associated with the instantaneous advection that is implied by assuming a steady state with a moving shock are isolated in ${\tilde{B}}_{s}$ and subsequently neglected; see Appendix A.

Next, there is work done by the PNS, arising from static and ram pressures at R_P. In general, static pressure dominates, so that W_PNS can be approximated by

$$\begin{eqnarray}&&{W}_{\mathrm{PNS}}\simeq 4\pi {R}_{P}^{3}{P}_{P}.\end{eqnarray} \tag{ 31 }$$

Finally, the energy advected across the boundaries is

$$\begin{eqnarray}&&{W}_{B}\simeq \alpha {({{GM}}_{P})}^{1/2}| {\dot{M}}_{0}| {R}_{S}^{1/2}\\ &&+\displaystyle \frac{| {\dot{M}}_{0}| }{7\alpha {({{GM}}_{P})}^{1/2}}\displaystyle \frac{d}{{dt}}\left\{\left[(\beta -1){\left(\displaystyle \frac{{R}_{P}}{{R}_{S}}\right)}^{2}+1\right]\displaystyle \frac{d}{{dt}}\left({R}_{S}^{7/2}\right)\right\},\end{eqnarray} \tag{ 32 }$$

where the second term is the correction that arises owing to the motion of the shock.

Substituting our result for ${({d}^{2}I/{{dt}}^{2})}_{\mathrm{QS}}$ in Equation (20), we arrive at a closed form expression for the shock acceleration in the quasi-stationary approximation:

$$\begin{eqnarray}&&{a}_{S}\approx \displaystyle \frac{{{\mathbb{E}}}_{\mathrm{QS}}+{W}_{\mathrm{PNS}}-| {\dot{M}}_{0}{R}_{S}{u}_{0}| -{V}_{S}^{2}(\partial \zeta /\partial {R}_{S})}{\chi \zeta },\end{eqnarray} \tag{ 33 }$$

where

$$\begin{eqnarray}&&\zeta \equiv \displaystyle \frac{\left(\beta -1\right)\left({R}_{S}^{2}-{R}_{P}^{2}\right)| {\dot{M}}_{0}| }{2| {u}_{0}| }\end{eqnarray} \tag{ 34 }$$

and

$$\begin{eqnarray}&&\chi \equiv \displaystyle \frac{1}{\beta -1}\displaystyle \frac{\partial }{\partial {V}_{S}}\left[{V}_{S}\left(\beta -1\right)\right].\end{eqnarray} \tag{ 35 }$$

All the variables in Equation (33) are determined by the set of external parameters of the flow (including P_P, which we fix at the beginning of the simulation) and the dynamical variables ${R}_{S}$ and ${V}_{S}$ (recall that $\beta =\beta ({R}_{S},{V}_{S},{M}_{P},{\dot{M}}_{0},{q}_{d})$).

5.2. The ${a}_{S}({R}_{S},{V}_{S})$ Phase Space

With the aid of the quasi-stationary approximation, we calculate the partition of the ${R}_{S}-{V}_{S}$ plane into regions of positive and negative acceleration of the shock. Given a set of external parameters, ${\dot{M}}_{0},$ L, M_P, and ${R}_{P},$ and the conditions of the initially stationary shock, determined by ${R}_{S}(t=0)$ and ${V}_{S}(t=0)=0$ (which also dictate the boundary conditions at the PNS), every point in the phase space is calculated with the relevant quasi-stationary profile, yielding an estimate of the shock acceleration. A set of examples with $| {\dot{M}}_{0}| =0.8\;{M}_{\odot }\ {{\rm{s}}}^{-1},$ ${M}_{P}=1.3\;{M}_{\odot },$ ${R}_{P}=40\;\mathrm{km}$, and ${R}_{S}(t=0)=100\;\mathrm{km}$ are shown in Figure 4. The figure shows the predicted values of a_S in the ${R}_{S}\mbox{--}{V}_{S}$ plane for three neutrino luminosities that are sufficient to drive a runaway expansion (see Figure 1). The actual evolutionary path of the shock radius and velocity, as found in each simulation, is superimposed on each of the plots.

The most prominent feature in Figure 4 is the distinct separation of the ${R}_{S}\mbox{--}{V}_{S}$ plane into regions of positive and negative acceleration. At small and large shock radii the acceleration is positive, separated by a trough of negative acceleration. In the region of small shock radii, the acceleration is dominated by the effect of the PNS surface pressure, which is why in this region the dependence on the shock velocity is relatively weak, while the spatial gradient of the acceleration is quite steep. In contrast, for larger shock radii, a positive acceleration depends on a more complicated combination of the total energy in the flow and on the boundary conditions at the shock itself, both of which depend on the shock velocity as well.

Zoom In Zoom Out Reset image size

A key result is the fact that the actual evolutionary paths of the shocks in the ${R}_{S}\mbox{--}{V}_{S}$ plane follow the predicted acceleration of the quasi-stationary model very closely, oscillating between negative and positive velocities in accordance with the sign of the acceleration. We find that the approximation for the shock acceleration is qualitatively robust but tends to overestimate its magnitude when compared to the full simulations (the term ${\tilde{B}}_{S}$ neglected above does tend to reduce $| {a}_{S}|$). Nonetheless, the contours of zero acceleration usually correspond very accurately to the points along the path where the shock acceleration in the simulations changes sign, especially for oscillating shocks. This result lends strong support to the validity of the quasi-stationary approach and the phase-space analysis.

The transition to a runaway expansion occurs when the oscillations grow sufficiently large such that a velocity of order $1000\;\mathrm{km}\;{{\rm{s}}}^{-1}$ pushes the shock through the negative acceleration trough into the region of positive acceleration on the right-hand side of the phase space. Typically, the path passes close to the saddle point of ${a}_{S}({R}_{S},{V}_{S})$ and arrives in the right-handed region of positive acceleration with a positive velocity, so the shock continues to expand. Hence, when the shock arrives in this region, the oscillations (if present) cease, and exponential expansion ensues. We hereon refer to the right curve (larger ${R}_{S}$) of zero acceleration as the "critical" a_S = 0 curve. In contrast, the left a_S = 0 curve (smaller ${R}_{S}$) corresponds to the change of sign of the shock acceleration in oscillatory motion and does not necessarily signify a transition to runaway expansion. We refer to this curve as the "oscillatory" a_S = 0 curve.

It is important to note that we do not find a case where the shock acceleration becomes negative at larger radii, to the right of the critical curve. The reason for this trend is that as the shock radius increases further, the competition between heating and other energetic effects maintains a positive acceleration (Equation (33)). Specifically, while $\tilde{{\rm{\Omega }}}$ and ($-| {\dot{M}}_{0}{R}_{S}{u}_{0}|$) become more negative with increasing shock radius, the effective heating term, $\tilde{U},$ becomes larger (more positive) and generally dominates the change in a_S (we note that W_PNS remains constant while $\tilde{K}$ and the last term in Equation (33) are essentially negligible; see Appendix B). Hence, once the evolving shock crosses the critical a_S curve from left to right, runaway expansion is inevitable.

The magnitude of the neutrino luminosity has a substantial effect on the acceleration phase space, as is readily seen when comparing the three plots in Figure 4. Clearly, as we raise L, the negative acceleration trough becomes narrower and shallower, and so the shock requires fewer oscillations to reach a point in the ${R}_{S}\mbox{--}{V}_{S}$ plane that will carry it over to the positive acceleration region. To further elucidate this point, we show in Figure 5 the phase space and evolutionary paths for a subcritical luminosity of ${L}_{52}=6$ (see Figure 1), for a transitional luminosity ${L}_{52}=9.7,$ and for a high luminosity of ${L}_{52}=11.$ In the subcritical case, small perturbations are damped owing to local stability (see Section 5.3), eventually converging onto a stationary solution, whereas larger perturbations are bounded owing to the wide and deep negative acceleration area in the phase space. As the luminosity increases, criticality is reached as the damping turns into anti-damping, and the shock continues to expand through a series of oscillations. Increasing the luminosity further narrows the negative acceleration trough until eventually the two a_S = 0 curves intersect. At even higher luminosities the curves detach again, and this detachment reflects a qualitative change in the explosive evolution. Now the initial profile lies on the critical a_S = 0 curve rather than on the oscillatory one (recall that we dictate initial profiles with ${V}_{S}({R}_{S}=100\;\mathrm{km})=0$ and ${a}_{S}({R}_{S}=100\;\mathrm{km})=0$). Such an initial profile will necessarily end in runaway expansion if the initial perturbation drives it into larger shock radii, but we find that it may also reach runaway expansion after one oscillation, if the perturbation is in the opposite, radially negative, direction. We conclude that the intersection of the a_S = 0 curves corresponds to a transition between oscillatory and direct runaway expansion.

Zoom In Zoom Out Reset image size

5.3. Paths to Runaway Expansion

The phase-space interpretation leads us to the following conclusion: even when all external parameters—$\{{\dot{M}}_{0},{M}_{P},{R}_{P},{P}_{P}$ and $L\}$—are set, both the initial shock radius and its velocity determine whether or not the flow will become unstable and reach runaway expansion. We demonstrate this by means of an example in Figure 6, in which we repeat the simulations in the case of ${L}_{52}=7.5,$ including the same pressure at the PNS, but start with static ${V}_{S}(t=0)=0$ shocks at different shock radii. In order to show the evolutionary paths corresponding to these different radii superimposed on a single ${a}_{S}({R}_{S},{V}_{S})$ phase-space map, we do not adjust the boundary PNS pressure according to each initial shock radius (i.e., P_P is set according to ${R}_{S}(t=0)=100\;\mathrm{km}$). Consequently, most profiles have a non-zero initial shock acceleration and evolve to cover a wide range of shock radii and velocities. The plot shows that there exist two regimes, at small and large initial radii, for which the flow results in a runaway expansion, while in an intermediate regime of shock radii it does not. Conversely, we can consider any combination of $\{{R}_{S},{V}_{S}\}$ as an initial profile for the simulation. There is a finite region in the ${R}_{S}\mbox{--}{V}_{S}$ plane (its border lying between the paths labeled R_S = 130 km and R_S = 140 km), for which oscillations are eventually damped, and the flow settles on the stationary solution. The physical reason for this limited range of $\{{R}_{S},{V}_{S}\}$ for which runaway expansion is denied (in this particular setup) is that at sufficiently small radii the pressure close to the PNS creates a large positive shock acceleration, which can drive the shock all the way beyond the critical a_S = 0 curve. Hence, runaway expansion will occur either if the flow is initiated in this region or if it is initiated in combinations of $\{{R}_{S},{V}_{S}\}$ that eventually evolve so that the shock penetrates deeply into the left-handed side region of positive acceleration. Of course, the critical a_S = 0 curve serves as an additional limit: a profile that is initiated in that region will accelerate exponentially if the initial shock velocity is positive (or even only slightly negative).

Zoom In Zoom Out Reset image size

In the general case of arbitrary initial $\{{R}_{S},{V}_{S}\}$ the outcome of the evolution can be found by calculating the entire evolutionary sequence (directly or by following it after mapping the appropriate phase space). Moreover, when the problem is restricted to initially stationary profiles $({V}_{S}(t=0)=0,\;{a}_{S}(t=0)=0),$ the situation is greatly simplified and we find that the outcome of the evolution can be assessed by the initial signs of the partial derivatives $\partial {a}_{S}/\partial {R}_{S}$ and $\partial {a}_{S}/\partial {V}_{S}.$ This can be seen by noting that here, perturbations are governed by

$$\begin{eqnarray}&&{\ddot{R}}_{S}={a}_{S}=\displaystyle \frac{\partial {a}_{S}}{\partial {R}_{S}}({R}_{S}-{R}_{S,0})+\displaystyle \frac{\partial {a}_{S}}{\partial {V}_{S}}{V}_{S},\end{eqnarray} \tag{ 36 }$$

which is simply the damped (or anti-damped) linear oscillator equation. We find that a linear stability analysis is generally applicable to the evolution of the shock, suggesting the following behavior:

• 1.  
  An initially stationary accretion layer with $\partial {a}_{S}/\partial {R}_{S}\gt 0$ necessarily corresponds to a point on the critical a_S = 0 curve. As is the case in the bottom panel of Figure 5, we find that such a profile will invariably evolve to a runaway expansion in a non-oscillatory fashion. In essence, this class of initially stationary profiles is inherently unstable.
• 2.  
  An initially stationary profile with $\partial {a}_{S}/\partial {R}_{S}\lt 0$ necessarily corresponds to an $\{{R}_{S},{V}_{S}=0\}$ point on the oscillatory a_S = 0 curve. Whether an initial perturbation diverges or damps out can be gauged locally by the sign of $\partial {a}_{S}/\partial {V}_{S}.$ If $\partial {a}_{S}/\partial {V}_{S}\lt 0,$ we expect a stable scenario: oscillations will tend to damp out, and the flow will settle onto the stationary solution. In contrast, if this derivative is positive, runaway expansion is expected, but its nature is more subtle. For small perturbations, we would expect that a further distinction will arise from the sign of ${\rm{\Delta }}\equiv {(\partial {a}_{S}/\partial {V}_{S})}^{2}+4\partial {a}_{S}/\partial {R}_{S}.$ If ${\rm{\Delta }}\lt 0,$ the shock should evolve through a series of increasing ("anti-damped") oscillations, finally resulting in runaway expansion. On the other hand, the combination ${\rm{\Delta }}\gt 0$ and $\partial {a}_{S}/\partial {V}_{S}\gt 0$ should lead to a non-oscillatory, runaway expansion (even though $\partial {a}_{S}/\partial {R}_{S}\lt 0$).

We test this insight in Figure 7. In the figure we show the results of a survey of simulations, presenting the critical luminosities as a function of mass accretion rate for ${M}_{P}=1.3\;{M}_{\odot }$ and ${R}_{P}=40\;\mathrm{km}$. Each simulation is initiated with a small perturbation around a stationary profile at ${R}_{S}=100\;\mathrm{km}$. The survey allows us to differentiate between stable configurations and runaway expansion and also to distinguish whether runaway occurs through an oscillatory or a non-oscillatory path. We define runaway expansion as non-oscillatory when there is no more than one cycle in which the shock radius drops below its initial value; for example, the ${L}_{52}=9$ case in Figure 1 is thus an example of a marginally non-oscillatory runaway expansion. The figure shows the critical luminosities for both oscillatory and non-oscillatory runaway expansion, which we denote ${L}_{c,\mathrm{osc}}$ and ${L}_{c,\mathrm{dir}},$ respectively. These data are based on our highest resolution, and through extrapolating our numerical convergence we estimate that the results are accurate to within 5%. Error bars delineate the range found in the simulations. The simulated critical luminosities are compared to curves that show the loci of $\partial {a}_{S}/\partial {R}_{S}=0,$ of $\partial {a}_{S}/\partial {V}_{S}=0,$ and of ${\rm{\Delta }}=0,$ as calculated with the quasi-stationary approximation (i.e., finding all combinations of L and $| {\dot{M}}_{0}|$ for which an initially stationary profile with ${R}_{S}(t=0)=100\;\mathrm{km}$ yields these derivatives).

As seen in the figure, there is a good qualitative consistency between the simulated critical luminosities and the theoretical curves, evaluated by the shock acceleration in the quasi-stationary approximation. The $\partial {a}_{S}/\partial {R}_{S}=0$ curve is especially indicative of the threshold for runaway expansion through the non-oscillatory path at higher mass accretion rates, for which $\partial {a}_{S}/\partial {V}_{S}\lt 0$ for all luminosities. Interestingly, the $\partial {a}_{S}/\partial {V}_{S}=0$ curve strongly depends on the accretion rate, requiring exceedingly higher L for increasing $| {\dot{M}}_{0}| ,$ eventually diverging at $| {\dot{M}}_{0}| \simeq 0.94\;{M}_{\odot }\;{{\rm{s}}}^{-1}.$ On the other hand, the $\partial {a}_{S}/\partial {R}_{S}=0$ curve has a rather weak dependence on the mass accretion rate, so only it survives at larger accretion rates. Furthermore, the $\partial {a}_{S}/\partial {R}_{S}=0$ and $\partial {a}_{S}/\partial {V}_{S}=0$ curves intersect at a slightly lower, critical, mass accretion rate of ${\dot{M}}_{c}\simeq 0.91\;{M}_{\odot }\;{{\rm{s}}}^{-1}.$ This intersection reflects the fact that for accretion rates of $| {\dot{M}}_{0}| \gt {\dot{M}}_{c},$ initial profiles with ${R}_{S}(t=0)=100\;\mathrm{km}$ that explode lie on the critical a_S = 0 curve, rather than the oscillatory one. As is indeed confirmed in the simulations, runaway expansion in this regime can occur only through a non-oscillatory evolution. We comment that at very high mass accretion rates, the a_S phase-space analysis becomes increasingly inaccurate owing to a large contribution from the Bernoulli function, which cannot be neglected in highly non-oscillatory expansions, when the shock dynamical time becomes comparable to the advection time (see Appendix A).

For lower mass accretion rates, $| {\dot{M}}_{0}| \leqslant {\dot{M}}_{c},$ the behavior of profiles with $\partial {a}_{S}/\partial {R}_{S}\lt 0$ is somewhat more complicated. Nevertheless, it too is qualitatively consistent with small perturbations in the quasi-stationary approximation. First, for profiles with $\partial {a}_{S}/\partial {V}_{S}\lt 0$ (low luminosities), runaway expansion is inherently denied: initial perturbations are damped, and the system relaxes to the stationary profile, as expected. The $\partial {a}_{S}/\partial {V}_{S}=0$ curve, above which oscillations are expected to anti-damp and lead to runaway expansion, indeed lies very close to the critical luminosity curve ${L}_{c,\mathrm{osc}},$ and hence provides a very good estimate of the critical condition for this type of evolution. The two curves do not quite coincide, and for most ${\dot{M}}_{0}$ there is a small range of luminosities for which $\partial {a}_{S}/\partial {V}_{S}\gt 0$ but runaway expansion does not yet occur. We find that in this range oscillations commence, but they tend to settle on a constant or almost-constant (increasing very slowly) amplitude. An example of such an evolution is presented in Figure 8, which shows the time-dependent radius of the shock for $| {\dot{M}}_{0}| =0.4\;{M}_{\odot }\;{{\rm{s}}}^{-1}$ and different neutrino luminosities. There is a general resemblance to Figure 1, but the case of ${L}_{52}=2.5$ is unique: initially the oscillations grow, but then settle on an approximately constant amplitude. Presumably, the existence of such an evolutionary path is due to the non-local nature of the ${a}_{S}({V}_{S})$ relation (in other words, the oscillator in Equation (36) develops larger amplitudes where it becomes nonlinear).

Figure 7 also shows the limiting luminosities for which there is a transition from oscillatory to non-oscillatory runaway expansion in the $\partial {a}_{S}/\partial {V}_{S}\gt 0$ region. According to linear theory, this curve should coincide with the ${\rm{\Delta }}=0$ curve (which is very close to the $\partial {a}_{S}/\partial {R}_{S}=0$ curve, and so is of little importance for the phase-space analysis). However, in this particular aspect we are strongly constrained by numerical resolution, which prevents us from initiating the simulation with truly infinitesimal perturbations. In general, it is impractical to initiate the simulations with perturbations that are smaller in amplitude than about 200 m, and the eventual dynamics of the shock depend on the changes of $\partial {a}_{S}/\partial {V}_{S}$ and $\partial {a}_{S}/\partial {R}_{S}$ over this scale. In the simulations, a perturbation of this magnitude can grow in radius by about 1 km during a single oscillation, which is significant in terms of the gradients of the acceleration. We show this by plotting in Figure 7 a curve that corresponds to a growth of the initial perturbation to 1 km in a single oscillation.³ The quantitative agreement between this curve and the critical non-oscillatory luminosity indeed suggests that the finite resolution causes us to identify oscillatory runaway as non-oscillatory when the linear analysis breaks down over a 1 km scale—obviously below the ${\rm{\Delta }}=0$ curve. In any case, this transitional luminosity is less important than the critical luminosities, which separate exploding and non-exploding profiles in the study of shock revival.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.4. Dependence on the Initial Shock Radius

So far we have set the initial shock radius arbitrarily at ${R}_{S}(t=0)=100\;\mathrm{km}$. Our phase-space analysis suggests that there should be a non-trivial dependence of the evolution, and therefore of the critical luminosity, on this radius, even when the set of parameters ${\dot{M}}_{0},{M}_{P},$ and R_P is fixed. In Figure 9 we confirm this expectation for the case of initial stationary profiles, where we compare the critical luminosity as a function of the external accretion rate, ${L}_{\mathrm{crit}}({\dot{M}}_{0}),$ when the profiles were constrained by two different choices of the initial shock radius: ${R}_{S}(t=0)=100\;\mathrm{km}$ and ${R}_{S}(t=0)=120\;\mathrm{km}$ (in all simulations, ${M}_{P}=1.3\;{M}_{\odot }$ and ${R}_{P}=40\;\mathrm{km}$). We also compare the results with the critical luminosity found when the initial conditions are determined by requiring a $\tau =2/3$ neutrino optical depth through the accretion layer, so that ${R}_{S}(t=0)$ is constrained for each simulation separately; for all but the lowest $| {\dot{M}}_{0}| ,$ this condition forces a larger initial shock radius, ${R}_{S}(t=0)=150\mbox{--}170\;\mathrm{km}$. It is evident from the figure that a dependence on ${R}_{S}(t=0)$ does indeed exist (weak at small $| {\dot{M}}_{0}| ,$ increasingly pronounced for larger accretion rates). Our results in this regard are consistent with the conclusions of Fernández (2012) and are distinct from implications of the stationary model, which generates a single physical solution when constrained by the additional condition of a definite optical depth.

Zoom In Zoom Out Reset image size

Also shown in Figure 9 is an analogous set of critical luminosities, found in simulations that incorporated full dissociation across the shock, ${q}_{d}={q}_{\mathrm{Fe}}.$ The qualitative dependence on ${R}_{S}(t=0)$ is retained when full dissociation is included, but all the critical luminosities are shifted upward, as expected. Notably, in this case the shock radius corresponding to $\tau =2/3$ generally lies between 100 and 120 km and so is the intermediate, rather than the lowest, critical luminosity curve. Once again, the general trend of the results suggests that our principal conclusions can be carried over to the more realistic case that incorporates dissociation across the shock and later during accretion.

The Burrows & Goshy (1993) critical luminosity, ${L}_{c,\mathrm{BG}}({\dot{M}}_{0}),$ is determined by the absence of a stationary solution. In cases where a stationary $\tau =2/3$ solution exists, but its instability leads to runaway expansion, the critical luminosity in simulations, ${L}_{\mathrm{crit}}({\dot{M}}_{0}),$ must be lower than ${L}_{c,\mathrm{BG}}({\dot{M}}_{0}).$ A subtle point here is that when using the $\tau =2/3$ constraint, ${L}_{c,\mathrm{BG}}$ is generally only a mild overestimate of L_crit even though the issue of stability is ignored (Fernández 2012). We qualitatively relate this to the phase-space analysis by noticing that in stationary models, as the neutrino luminosity is increased toward ${L}_{c,\mathrm{BG}},$ the initial shock radius changes significantly; see Figure 10. As a result, the phase-space structure of ${a}_{S}({R}_{S},{V}_{S})$ too changes in a substantial manner and toward instability, which is therefore reached at luminosities close to (but lower than) ${L}_{c,\mathrm{BG}}.$ Thus, for any given external mass accretion rate, the two critical luminosities, corresponding to the disappearance of a stationary solution and to its instability, lie in close proximity.

Zoom In Zoom Out Reset image size

Our phase-space interpretation, combined with the quasi-stationary approximation, enabled us to identify a simple criterion for determining if an initially stationary profile will explode, based on the signs of the partial derivatives of the shock acceleration. We emphasize that even though they are expressed by quantities pertaining to the shock, these are in fact global criteria, in the sense that they depend on the integral properties of the entire accretion layer. In this section we use the phase-space analysis and the quasi-stationary approximation to evaluate a commonly discussed timescale criterion for the onset of an explosion, comparing heating and advection in the so-called gain region of the flow (Murphy & Burrows 2008; Fernández 2012; Pejcha & Thompson 2012; Ott et al. 2013; Dolence et al. 2015).

The specific neutrino heating scales as the squared inverse of the radial position of the mass element in the accretion layer, while the specific neutrino cooling scales as the sixth power of the element's temperature (recall Equations (12)–(13)). As a result, at smaller radii where the temperatures are higher, cooling dominates over heating, $\dot{q}(r)\lt 0,$ whereas closer to the shock heating dominates and $\dot{q}(r)\gt 0.$ Such a structure is also found in supernova simulations, including studies that use more accurate modeling than our simplified equations (Fernández 2012; Ott et al. 2013). Correspondingly, it is common practice to identify the radial position in the flow where $\dot{q}=0$ as the "gain radius," R_G, and refer to the region extending between this radius and the shock as the "gain region."

The argument concerning the timescales is that if in the gain region the heating time,

$$\begin{eqnarray}&&{t}_{\mathrm{heat},G}=\displaystyle \frac{{E}_{G}}{{\displaystyle \int }_{{R}_{G}}^{{R}_{S}}\dot{q}{dm}},\end{eqnarray} \tag{ 37 }$$

is smaller than the corresponding advection time in this region,

$$\begin{eqnarray}&&{t}_{\mathrm{adv},G}={\displaystyle \int }_{{R}_{G}}^{{R}_{S}}\displaystyle \frac{{dr}}{u}\approx \displaystyle \frac{{M}_{G}}{| \dot{M}| },\end{eqnarray} \tag{ 38 }$$

then the flow should evolve toward runaway expansion. In these equations, E_G and M_G are the total internal energy of matter and the total mass, respectively, in the gain region, and the second equality for t_adv is exact if the mass accretion rate is uniform. The criterion ${t}_{\mathrm{heat},G}\lt {t}_{\mathrm{adv},G}$ reflects the expectation that if the material can heat significantly before it passes through the gain region, steady-state accretion could be disturbed. In the following, we demonstrate that the phase-space structure of a_S can be approximately related to this ratio of timescales.

First, note that in Equation (33), the factor ζ is invariably positive, and χ is always positive (because ${V}_{S}(\beta -1)$ tends to increase with shock velocity). Hence, the sign of a_S is determined by the numerator of this equation. This numerator includes the effective energy ${{\mathbb{E}}}_{\mathrm{QS}}\approx \tilde{K}+\tilde{U}+\tilde{{\rm{\Omega }}},$ which characterizes the entire accretion layer, and additional terms that include contributions from the PNS and the shock. In Appendix B we show that in the region of interest, the PNS term ${W}_{\mathrm{PNS}}\simeq 4\pi {R}_{P}^{3}{P}_{P}$ dominates over the other terms, and so we can generally claim that ${{\mathbb{E}}}_{\mathrm{QS}}\gt 0$ along with ${V}_{S}\gt 0$ provide a sufficient condition for runaway expansion. We show this explicitly in Figure 11, in which we plot again the shock acceleration in the ${R}_{S}\mbox{--}{V}_{S}$ plane, and superimpose on it the ${{\mathbb{E}}}_{\mathrm{QS}}=0$ curve. As expected, the ${{\mathbb{E}}}_{\mathrm{QS}}=0$ curve lies in the runaway expansion region to the right of the critical a_S = 0 curve, clarifying that a positive effective energy is indeed a sufficient—but not necessary—condition.

Zoom In Zoom Out Reset image size

The total energy ${{\mathbb{E}}}_{\mathrm{QS}}$ involves an integration over the entire accretion layer. Conversely, we can define a radius-dependent effective energy, ${\mathbb{E}}(r),$ corresponding to the integrals (25)–(27) taken from r to R_S. In general, we expect ${\mathbb{E}}(r)\gt {\mathbb{E}}({R}_{P})\equiv {{\mathbb{E}}}_{\mathrm{QS}},$ since matter close to the PNS experiences the strongest cooling and gravitational potential. As a result, there can exist combinations of ${R}_{S}$ and ${V}_{S}$ for which the profile has ${\mathbb{E}}(r)\gt 0$ for some range of radii ${R}_{P}\lt r\lt {R}_{S},$ while ${{\mathbb{E}}}_{\mathrm{QS}}\lt 0.$ Numerically, we find that for most initial configurations that eventually lead to runaway expansion, the function ${\mathbb{E}}(r)$ has a maximum at some inner radius in the accretion layer, which we denote as R_E. This radius is always close to but slightly lower than R_G: just below the gain radius the effective energy still increases inward, owing to the positive contribution of $\tilde{U}$ (the integrated heating, $Q(r),$ remains positive for some distance below R_G, so in this region $\partial \tilde{U}/\partial r\lt 0$). We demonstrate the locations of R_E and R_G in Figure 12, which shows ${\mathbb{E}}(r)$ in the accretion layer calculated for various initial shock radii and a given set of external parameters (${L}_{52}=8,$ $| {\dot{M}}_{0}| =0.8\ {M}_{\odot }\;{{\rm{s}}}^{-1},\;{M}_{P}=1.3\ {M}_{\odot },$ and q_d = 0).

Zoom In Zoom Out Reset image size

Repeating the process of identifying and calculating R_E and R_G as a function of R_S and V_S for a given set of external conditions, we find the loci of all the points in the ${R}_{S}\mbox{--}{V}_{S}$ plane for which ${\mathbb{E}}({R}_{G})=0$ and ${\mathbb{E}}({R}_{E})=0.$ The two corresponding curves are also shown in Figure 11 and can be seen to lie very close to each other, both to the left of the critical a_S = 0 curve. Our conclusion is that ${\mathbb{E}}({R}_{E})\gt 0$ and ${\mathbb{E}}({R}_{G})\gt 0$ are both thresholds that must be crossed when the flow evolves to runaway expansion, but they are not sufficient conditions to ensure that runaway will occur. Indeed, for a wide range of parameters, crossing the ${\mathbb{E}}({R}_{G})=0$ curve does not ensure runaway expansion.

The quasi-stationary approximation can now be used to relate the condition ${t}_{\mathrm{adv},G}/{t}_{\mathrm{heat},G}\gt 1$ to the condition ${\mathbb{E}}({R}_{G})\gt 0.$ With some manipulation of Equations (58)–(65), we derive that

$$\begin{eqnarray}&&{\mathbb{E}}({R}_{G})\approx {K}_{G}+{U}_{G}+\left(2{\hat{\gamma }}_{{E}_{G}}-3\right){E}_{G},\end{eqnarray} \tag{ 39 }$$

with K_G, E_G, and U_G being the total kinetic energy, total internal energy, and energy gain between R_G and R_S, respectively (the latter is the gain region equivalent to Equation (59)), and ${\hat{\gamma }}_{{E}_{G}}$ is the mean adiabatic index of the material in the gain region, weighted by the internal energy:

$$\begin{eqnarray}&&{\hat{\gamma }}_{{E}_{G}}=\displaystyle \frac{1}{{E}_{G}}{\displaystyle \int }_{{R}_{G}}^{{R}_{S}}\gamma {edm}.\end{eqnarray} \tag{ 40 }$$

Since E_G is positive by construction, the condition for ${\mathbb{E}}({R}_{G})\gt 0$ becomes

$$\begin{eqnarray}&&\displaystyle \frac{{K}_{G}}{{E}_{G}}+\displaystyle \frac{1}{{E}_{G}}{\displaystyle \int }_{{R}_{G}}^{{R}_{S}}\displaystyle \frac{1}{| \dot{M}| }Q(r){dm}\gt 3-2{\hat{\gamma }}_{{E}_{G}},\end{eqnarray} \tag{ 41 }$$

where Q(r) was defined in Equation (29). Since $\dot{q}\geqslant 0$ in the gain region, Q(r) monotonically decreases with the radius, and $Q({R}_{S})=0.$ The total non-adiabatic integral over the entire gain region can therefore be approximated by

$$\begin{eqnarray}&&{\displaystyle \int }_{{R}_{G}}^{{R}_{S}}Q(r){dm}=\eta Q({R}_{G}){M}_{G},\end{eqnarray} \tag{ 42 }$$

where η is a numerical factor reflecting the gradient of Q(r) ($\eta =1/2$ for a linear function). We find that for the typical configurations in Figure 4, $\eta \sim 0.6$ and ${\hat{\gamma }}_{{E}_{G}}\sim 1.4.$

Finally, since the accretion layer typically involves low velocities, the first term in Equation (41) is negligible, and so the condition for ${\mathbb{E}}({R}_{G})\gt 0$ can be summarized as

$$\begin{eqnarray}&&\displaystyle \frac{{\displaystyle \int }_{{R}_{G}}^{{R}_{S}}\dot{q}{dm}}{{E}_{G}}\displaystyle \frac{{M}_{G}}{| \dot{M}| }\gt \displaystyle \frac{1}{\eta }\left(3-2{\hat{\gamma }}_{{E}_{g}}\right)\end{eqnarray} \tag{ 43 }$$

or (using Equations (37) and (38))

$$\begin{eqnarray}&&\displaystyle \frac{{t}_{\mathrm{adv},G}}{{t}_{\mathrm{heat},G}}\gt \displaystyle \frac{1}{\eta }\left(3-2{\hat{\gamma }}_{{E}_{G}}\right).\end{eqnarray} \tag{ 44 }$$

For the typical values of $\eta =0.6,{\gamma }_{{E}_{G}}=1.4,$ and zero dissociation, the expression on the right-hand side of these equations is about 0.3. A comparison of this result to actual simulations is presented in Figure 13. Indeed, runaway expansion is seen to commence in conjunction with the ratio of timescales in the gain region exceeding this threshold.

Zoom In Zoom Out Reset image size

We note that dissociation losses can be included in the estimate for the required value of ${t}_{{\mathrm{adv}}_{G}}/{t}_{\mathrm{heat},G}$ by replacing $Q(r)/| \dot{M}|$ by $[Q(r)/| \dot{M}| -{q}_{d}]$ in Equation (41). While this approximation is somewhat crude (since the quasi-stationary approximation is less appropriate for large dissociation losses), it does provide some quantitative estimate of the effect. The actual derivation then yields a correction of $+({q}_{d}{M}_{G})/({E}_{G}\eta )$ to the right-hand sides of Equations (43) and (44), demonstrating that a larger dissociation energy requires a larger timescale ratio, since heating must compensate for the lost energy. With the aid of the quasi-stationary approximation we predict that in the case of full dissociation, the necessary value of ${t}_{{\mathrm{adv}}_{G}}/{t}_{\mathrm{heat},G}$ is about 1.4. Again, this value appears to be consistent with the simulations, as demonstrated in Figure 13.

Finally, we point out that the value of ${t}_{\mathrm{adv},G}/{t}_{\mathrm{heat},G}$ increases from left to right (i.e., with increasing radii) in the ${R}_{S}\mbox{--}{V}_{S}$ plane, so a ratio of unity should in general not be far from the critical a_S = 0 curve. However, the two are by no means identical criteria, and we find no distinct relation between ${t}_{\mathrm{adv},G}/{t}_{\mathrm{heat},G}$ and any specific feature in the phase space. Moreover, the value of unity does not appear to hold any particular significance: it occurs well into runaway in the absence of dissociation, while being insufficient (occurring during the oscillations prior to runaway) when full dissociation is assumed.

We summarize our analysis of the ratio between advection and heating timescales with two claims:

• 1.  
  Strictly speaking, ${t}_{\mathrm{adv},G}/{t}_{\mathrm{heat},G}\gt 1$ is not a special condition for an explosion. Rather, it is a value that apparently lies in the path of a shock that is on its way to a runaway expansion and occurs close to (before or after) the transition to a positive shock acceleration.⁴
• 2.  
  We also do not identify the gain region as being in any way unique in terms of a necessary condition for runaway expansion. The curves in Figure 11 imply that ${\mathbb{E}}({R}_{E})\gt 0,$ or some other alternative curve, could equally well be used as a necessary, but not sufficient, condition. Clearly, runaway expansion requires significant heating, but evolution toward a runaway expansion is a quality of the entire accretion layer, rather than just of the gain region.

In this work we investigated the driving of a stalled accretion shock into runaway expansion in the context of the neutrino heating mechanism in core-collapse supernovae. We used spherically symmetric simulations in conjunction with analytic derivations and examined the evolutionary paths the system may take toward runaway expansion, if it occurs. Our numerical results included a parameter survey of the neutrino luminosity and the incoming accretion rate, as well as of the shock parameters, namely, its instantaneous radius and velocity. Our main conclusions are as follows:

• 1.  
  The instantaneous shock acceleration, a_S, depends on both the shock radius, R_S, and velocity, V_S. We examine the shock acceleration in the ${R}_{S}\mbox{--}{V}_{S}$ plane and show how it is generally divided into regions of positive and negative shock acceleration (see Figures 4 and 5). We expanded the commonly used stationary approximation, by including a non-zero shock velocity and solving the resulting profile in the accretion layer with modified boundary conditions (Section 5.1). The quasi-stationary approximation developed here allows us to map this plane and nicely predicts the evolution of a given accretion layer. The quasi-stationary approximation is especially accurate when the shock dynamical time, ${t}_{S}={R}_{S}/{V}_{S},$ and oscillation period, t_osc, are much longer than the advection time through the accretion layer, which is the case if the immediate dissociation energy losses across the shock are not too large. The quasi-stationary approximation can actually be used to estimate the oscillation period, and we show why it is typically 25–50 ms, very weakly dependent on the mass accretion rate (see Appendix C). This result is consistent with simulations.
• 2.  
  Generally, the ${R}_{S}-{V}_{S}$ phase space shows positive shock acceleration at small and large shock radii, with an intermediate trough of negative acceleration. The trough is bounded by two ${a}_{S}=0$ curves. The curve corresponding to smaller radii is an "oscillatory" ${a}_{S}=0$ curve: if the shock moves from it inward, positive acceleration due to the high pressure close to the PNS can cause the shock to bounce back and oscillate. Such oscillations either damp out or grow, depending on the parameters of the flow. In contrast, the ${a}_{S}=0$ curve at large radii is "critical," since crossing it with a positive velocity leads to a region of positive acceleration, so runaway expansion is guaranteed (i.e., profiles that correspond to this curve are inherently unstable). The magnitude of the neutrino luminosity comes into play through quantitative aspects of this phase-space structure. For low luminosities, the trough of negative acceleration is both deep and wide, so when the initial shock radius and velocity correspond to this trough, the shock is generally stable and oscillations are damped. When the neutrino luminosity is high, the trough of negative acceleration becomes shallow and narrow, and larger regions of it correspond to profiles that are unstable, so the shock eventually evolves beyond the critical ${a}_{S}=0$ curve and toward runaway expansion.
• 3.  
  For arbitrary combinations of ${R}_{S}$ and ${V}_{S}$ serving as initial conditions, the evolution of a model can be tracked in the phase space estimated by the quasi-stationary approximation. These evolutionary paths are confirmed to accurately match the shock evolution in simulations. Our analysis demonstrates that the shock velocity must be accounted for when assessing whether or not the evolution will end in runaway expansion.
• 4.  
  In the special case of initially stationary profiles (${V}_{S}=0,\;{a}_{S}=0$), we find that the partial derivatives of the shock acceleration offer a satisfactory indication regarding the outcome. If $\partial {a}_{S}/\partial {R}_{S}\gt 0$ and $\partial {a}_{S}/\partial {V}_{S}\lt 0$, the initial profile lies on the critical ${a}_{S}=0$ curve, and the profile will evolve to runaway expansion in a non-oscillatory fashion. In contrast, if $\partial {a}_{S}/\partial {R}_{S}\lt 0$, the initial profile lies on the oscillatory ${a}_{S}=0$ curve and its evolution depends on the sign of $\partial {a}_{S}/\partial {V}_{S}$. If this derivative too is negative, the system is stable and any small initial perturbation will damp out and settle to the stationary profile. On the other hand, if $\partial {a}_{S}/\partial {V}_{S}\gt 0$, oscillations will tend to grow unstably (be "anti-damped"), leading to runaway expansion through an oscillatory path. This analysis bares out very well in the simulations, in the sense that we can predict the critical luminosity as a function of mass accretion rate, ${L}_{\mathrm{crit}}({\dot{M}}_{0})$, for both oscillatory and non-oscillatory explosions; see Figure 7. Quantitatively, $\partial {a}_{S}/\partial {V}_{S}=0$ is a very good indicator for finding the critical luminosity for oscillatory explosions, and $\partial {a}_{S}/\partial {R}_{S}=0$ does well in predicting the critical luminosity for non-oscillatory explosions. We note that we slightly underestimate ${L}_{\mathrm{crit}}({\dot{M}}_{0})$ for the oscillatory mode, because there is a small range of luminosities that generate stable oscillations.
• 5.  
  Both modes of runaway expansion are available for low mass accretion rates. For a given low $| {\dot{M}}_{0}|$, oscillatory runaway with $\partial {a}_{S}/\partial {V}_{S}\gt 0$ is met at lower luminosities than $\partial {a}_{S}/\partial {R}_{S}\gt 0$, so the actual critical luminosity is due to anti-damped oscillations. We find that the luminosity at which $\partial {a}_{S}/\partial {V}_{S}=0,{L}_{c,\mathrm{osc}}$, is strongly dependent on the mass accretion rate, and there exists some finite rate for which ${L}_{c,\mathrm{osc}}$ diverges. Thus, there exists a critical mass accretion rate, ${\dot{M}}_{c}$, so that for $| {\dot{M}}_{0}| \gt {\dot{M}}_{c}$, only non-oscillatory runaway expansion is possible, with a critical luminosity ${L}_{c,\mathrm{dir}}$, which is weakly dependent on ${\dot{M}}_{0}$. The value of ${\dot{M}}_{c}$ may depend on the specifics of the model (PNS properties, initial shock radius, dissociation losses, etc.), but for reasonable parameters it appears to be about $1\;{M}_{\odot }\;{{\rm{s}}}^{-1}$.
• 6.  
  We also applied the quasi-stationary approximation to examine the commonly used ratio between the advection timescale and the heating timescale in the gain region of the accretion layer (see Section 6). We find that for the gain region to have a positive effective energy, this ratio must be at least a few tenths, with some dependence on the specific loss of energy due to dissociation across the shock (Figure 13). This result implies that the ratio of timescales, while not a fundamental condition for an explosion, can serve as an indication if a given system will evolve to an explosion, as suggested in previous works. However, there are several qualifications to this observation. First, a positive energy in the gain region appears to be a threshold that must be crossed on the path to runaway expansion, yet it is not an actual condition for an explosion, since it does not exactly coincide with the region of positive shock acceleration (Figure 11). Second, we find that in many cases the advection time becomes longer than the heating time only after the transition to a runaway expansion has occurred, implying that an equality between the two timescales can be a consequence of a successful explosion, rather than a physical condition for initiating one. Third, we do not identify the gain region as being unique in some physical sense; in principle, an alternative criterion can be formulated by applying similar considerations on a different region in the accretion layer, such as that which holds the largest total energy. We therefore suggest that the properties of the gain region should not be used as a singular measure of the likelihood of an explosion.

Since our simulations are inherently limited to spherical symmetry (with additional assumptions and simplifications), quantitative differences are to be expected in the actual problem of the neutrino driven mechanism. Most notably, a multi-dimensional flow allows for unstable non-radial modes, which could have shorter periods than the radial oscillations that we focused on here. These modes include the advective-acoustic instabilities (SASI—Standing Accretion Shock Instability) and convective instabilities (see, e.g., Yamasaki & Yamada 2007; Fernández et al. 2014; Foglizzo et al. 2015). The latter typically evolve over a Brunt–Väisälä timescale, while the former may be even faster. The linear stability analysis by Yamasaki & Yamada (2007) suggests that while radial modes appear to dominate at larger neutrino luminosities, at lower luminosities the cycle time of the advective-acoustic instabilities could be as short as the advection time so that they can drive exponential growth and oscillations in the linear regime. The competition between radial and non-radial modes is likely to be a delicate one even at low luminosities, especially for nonlinear or saturated modes; note that the difference in timescales is not large, and the hydrodynamic evolution usually spans several oscillations.

In principle, a generalized version of the formalism we introduced here can be developed to include non-radial modes since the fundamental underlying principles should be similar. For example, the concepts of representing the problem through an "extended virial theorem" (see Appendix A) and including the shock motion in an analysis of nonlinear stability in an appropriate phase space are quite general and could possibly be applied to non-radial modes.

Other aspects of the qualitative nature of our results should be applicable to the full, multi-dimensional problem. In particular, we emphasize the dependence of the shock acceleration on its radius and velocity; this fact certainly facilitates dynamical spherical simulations to explode with a lower neutrino luminosity than predicted by the stationary approximation. We speculate that this behavior does carry over to the three-dimensional reality, where the shock can sample a wider range of combinations of ${R}_{S}$ and ${V}_{S}$ than in the one-dimensional case. This should result from initial aspherical conditions in the silicon–oxygen layer (Couch & Ott 2015a; Couch et al. 2015), and also from the onset of turbulence (of course, turbulent pressure should also contribute to facilitating an explosion; see Couch & Ott 2015b; Fernández 2015). This combination may be responsible—at least in part—for the observation that the critical luminosity in two- and three-dimensional simulations is lower than in the one-dimensional case.

A natural expansion of this work would be to examine two- and three-dimensional simulations and to quantitatively compare the actual acceleration of the shock with that predicted from a spherical model of identical parameters. Such a comparison will clarify the dependence of the shock characteristics on the physical processes inside the accretion layer, such as turbulence and convection.

We thank E. Livne, L. Metzker, and O. Pejcha for helpful discussions. U.K. is supported by the European Union Seventh Framework Programme, by an IAEC-UPBC joint research foundation grant, and by an ISF-UGC grant.

In this appendix we derive the various terms that appear in the quasi-stationary approximation for ${d}^{2}I/{{dt}}^{2}$ (Equation (23)). Recall that the goal of the approximation is to include the effect of a finite shock velocity on the profile of the accretion layer. To do so, we assume that the profile adjusts instantaneously to changes in the shock radius and velocity. As we show below, this requires some care, because the advection time is not negligible.

A.1. Virial Theorem for Spherical Accretion

In general form, the first and second time derivatives of the spherical moment of inertia in the accretion layer are

$$\begin{eqnarray}&&\displaystyle \frac{{dI}}{{dt}}={\displaystyle \int }_{{R}_{p}}^{{R}_{s}}8\pi {r}^{3}\rho {udr}-{\left[{r}^{2}\dot{M}(r)-4\pi {r}^{4}\rho V\right]}_{P}^{S},\end{eqnarray} \tag{ 45 }$$

where we used the continuity Equation (2), and

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}I}{{d}^{2}t}={\displaystyle \int }_{{R}_{P}}^{{R}_{S}}8\pi {r}^{3}\displaystyle \frac{\partial (\rho u)}{\partial t}{dr}-{\left[{r}^{2}\ddot{M}-\displaystyle \frac{d}{{dt}}\left(4\pi {r}^{4}\rho V\right)\right]}_{P}^{S},\end{eqnarray} \tag{ 46 }$$

where $\dot{M}(r)=4\pi {r}^{2}\rho$u is the local mass accretion rate, $\ddot{M}=d\dot{M}/{dt}$ is its time derivative, and $V$ is the velocity of the flow boundaries. As defined in Section 5, the notation ${[\cdots ]}_{P}^{S}$ stands for the difference between the expression in the immediate downstream of the shock and just outside the PNS. The last term in Equation (46) is a Lagrangian time derivative of the boundaries.

Using the momentum conservation Equation (3), we express ${d}^{2}I/{{dt}}^{2}$ in a "virial theorem" fashion with boundary conditions

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}\displaystyle \frac{{d}^{2}I}{{{dt}}^{2}}={\mathbb{E}}+{\rm{\Delta }}{{\mathbb{E}}}_{\mathrm{PNS}}+{\rm{\Delta }}{{\mathbb{E}}}_{S},\end{eqnarray} \tag{ 47 }$$

where

$$\begin{eqnarray}&&{\mathbb{E}}=2\;K+3{\displaystyle \int }_{{R}_{p}}^{{R}_{s}}\displaystyle \frac{p}{\rho }{dm}+{\rm{\Omega }},\end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{{\mathbb{E}}}_{\mathrm{PNS}}={\left[4\pi {r}^{3}\left(p+\rho {u}^{2}\right)+\displaystyle \frac{1}{2}{r}^{2}\ddot{M}-\displaystyle \frac{d}{{dt}}\left(2\pi {r}^{4}\rho V\right)\right]}_{r\to {R}_{P}},\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{{\mathbb{E}}}_{S}={\left[-4\pi {r}^{3}\left(p+\rho {u}^{2}\right)-\displaystyle \frac{1}{2}{r}^{2}\ddot{M}+\displaystyle \frac{d}{{dt}}\left(2\pi {r}^{4}\rho V\right)\right]}_{r\to {R}_{S}},\end{eqnarray} \tag{ 50 }$$

and K and Ω are the total kinetic and gravitational energies, defined as

$$\begin{eqnarray}&&K\equiv {\displaystyle \int }_{{R}_{p}}^{{R}_{s}}\displaystyle \frac{1}{2}{u}^{2}{dm}\end{eqnarray} \tag{ 51 }$$

and

$$\begin{eqnarray}&&{\rm{\Omega }}\equiv -{\displaystyle \int }_{{R}_{p}}^{{R}_{s}}\displaystyle \frac{{GM}(m)}{r}{dm}\approx -{\displaystyle \int }_{{R}_{p}}^{{R}_{s}}\displaystyle \frac{{{GM}}_{p}}{r}{dm}.\end{eqnarray} \tag{ 52 }$$

The physical interpretation of ${\mathbb{E}}$ and ${\rm{\Delta }}{{\mathbb{E}}}_{\mathrm{PNS}}+{\rm{\Delta }}{{\mathbb{E}}}_{S}$ is the energetic contributions of the accretion layer and of the boundaries (at the PNS and at the shock, respectively). We note again that the boundary terms are calculated within the accretion layer. The time derivative $d/{dt}$ in Equations (49)–(50) is a Lagrangian derivative of values in the post-shocked region, which is identically zero for a stationary PNS surface.

The ${{dI}}^{2}/{{dt}}^{2}$ expression in Equations (47)–(52) is exact but cannot be evaluated directly in a quasi-stationary approximation. This is due to the fact that the $P/\rho$ term in ${\mathbb{E}}$ is continuously modified by the motion of the shock, and this effect cannot be neglected when the advection time is not negligible with respect to the characteristic time for changes in the shock velocity, e.g., the oscillation time. Therefore, we next seek an alternative formulation that is more susceptible to the quasi-stationary approximation.

A.2. The Effective Energy Term

We derive the effective energy term ${{\mathbb{E}}}_{\mathrm{QS}}$, in Equation (23), by rearranging the terms in ${\mathbb{E}}$ and applying the quasi-stationary approximation. Integrating the energy-conservation Equation (4) over volume in a quasi-stationary profile ($\partial /\partial t=$ 0) yields

$$\begin{eqnarray}&&\dot{M}\left[{b}_{S}-b(r)\right]=Q(r),\end{eqnarray} \tag{ 53 }$$

where $b(r)$ is the Bernoulli function (specific energy) at $r$,

$$\begin{eqnarray}&&b(r)\equiv \displaystyle \frac{1}{2}{u}^{2}+e+\displaystyle \frac{p}{\rho }-\displaystyle \frac{{{GM}}_{p}}{r},\end{eqnarray} \tag{ 54 }$$

${b}_{S}\equiv b({R}_{S})$, and $Q(r)$, also defined in the main text (Equation (29)), is the total non-adiabatic energy deposition rate between r and ${R}_{S}$,

$$\begin{eqnarray}&&Q(r)\equiv {\displaystyle \int }_{r}^{{R}_{s}}4\pi \rho {r}^{2}\dot{q}{dr}={\displaystyle \int }_{r}^{{R}_{s}}\dot{q}\;{dm}.\end{eqnarray} \tag{ 55 }$$

Using Equation (9) for energy conservation across the shock yields

$$\begin{eqnarray}{b}_{S} & = & {b}_{0}-{q}_{d}+{V}_{S}({u}_{1}-{u}_{0})\\ & = & {-q}_{d}-\left(1-\displaystyle \frac{1}{\beta }\right)\left(1-\displaystyle \frac{{V}_{S}}{{u}_{0}}\right){u}_{0}{V}_{S},\end{eqnarray} \tag{ 56 }$$

where the Bernoulli function just outside the shock, ${b}_{0}=0$, vanishes in the pre-shocked material (Equation (6)). The Bernoulli function inside the post-shocked flow is then obtained by

$$\begin{eqnarray}&&b(r)=\displaystyle \frac{1}{2}{u}^{2}+e+\displaystyle \frac{p}{\rho }-\displaystyle \frac{{{GM}}_{P}}{r}=\displaystyle \frac{1}{| \dot{M}| }Q(r)+{b}_{S},\end{eqnarray} \tag{ 57 }$$

and the ratio $P/\rho$ in the shock can be rewritten as

$$\begin{eqnarray}&&\displaystyle \frac{p}{\rho }=\displaystyle \frac{\gamma -1}{\gamma }\left(\displaystyle \frac{1}{| \dot{M}| }Q(r)+{b}_{s}+\displaystyle \frac{{{GM}}_{p}}{r}-\displaystyle \frac{1}{2}{u}^{2}\right),\end{eqnarray} \tag{ 58 }$$

where $\gamma$ is the local adiabatic index of the flow. The total non-adiabatic energy gained by the material in the accretion layer is then

$$\begin{eqnarray}&&U={\displaystyle \int }_{{R}_{P}}^{{R}_{S}}\displaystyle \frac{1}{| \dot{M}| }Q(r){dm}.\end{eqnarray} \tag{ 59 }$$

Substituting these derived quantities into the integrals included in the effective energy ${\mathbb{E}}$ from Equation (23), we arrive at a form for the effective total energy in the quasi-stationary approximation, ${{\mathbb{E}}}_{\mathrm{QS}}$,

$$\begin{eqnarray}&&{{\mathbb{E}}}_{\mathrm{QS}}=\tilde{K}+\tilde{U}+\tilde{{\rm{\Omega }}}+{\tilde{B}}_{S}.\end{eqnarray} \tag{ 60 }$$

The terms in Equation (60) are effective manipulations of the original functions, including coefficients that depend on the local adiabatic index in the accretion layer. We define $\tilde{K},\tilde{U},\tilde{{\rm{\Omega }}}$, and ${\tilde{B}}_{S}$ as the effective kinetic energy, effective gained energy, effective gravitational potential, and effective Bernoulli function, respectively:

$$\begin{eqnarray}&&\tilde{K}\equiv {\displaystyle \int }_{{R}_{p}}^{{R}_{s}}\left(\displaystyle \frac{3-\gamma }{\gamma }\right)\displaystyle \frac{1}{2}{u}^{2}\;{dm},\end{eqnarray} \tag{ 61 }$$

$$\begin{eqnarray}&&\tilde{{\rm{\Omega }}}\equiv {\displaystyle \int }_{{R}_{p}}^{{R}_{s}}\left(\displaystyle \frac{3-2\gamma }{\gamma }\right)\left(-\displaystyle \frac{{{GM}}_{p}}{r}\right)\;{dm},\end{eqnarray} \tag{ 62 }$$

$$\begin{eqnarray}&&\tilde{U}\equiv {\displaystyle \int }_{{R}_{p}}^{{R}_{s}}\left(3\displaystyle \frac{\gamma -1}{\gamma }\right)\left(\displaystyle \frac{1}{| \dot{M}| }Q(r)\right){dm},\end{eqnarray} \tag{ 63 }$$

and

$$\begin{eqnarray}&&{\tilde{B}}_{S}\equiv {\displaystyle \int }_{{R}_{p}}^{{R}_{s}}\left(3\displaystyle \frac{\gamma -1}{\gamma }\right){b}_{S}(m)\;{dm}.\end{eqnarray} \tag{ 64 }$$

The last expression warrants some discussion. The quantity ${b}_{S}$ in the integral should reflect the Bernoulli function ${b}_{S}(m)$ of the mass element, $m$, when it crossed the shock, and so is not a constant across the accretion layer. For example, when dissociation is neglected, ${b}_{S}={V}_{s}({u}_{s}-{u}_{0})$, which vanishes in the stationary case and fluctuates between positive and negative values as the shock oscillates. In general, for small shock velocities, this causes the integral ${\tilde{B}}_{S}$ to be essentially negligible with respect to the other terms of the effective energy, which we confirm explicitly in the simulations. Hence, in the quasi-stationary approximation:

$$\begin{eqnarray}&&{{\mathbb{E}}}_{\mathrm{QS}}\approx \tilde{K}+\tilde{U}+\tilde{{\rm{\Omega }}}.\end{eqnarray} \tag{ 65 }$$

A.3. The Boundary Terms

The boundary terms in ${d}^{2}I/{{dt}}^{2}$ include distinct contributions from the PNS and from the shock. A key feature of our model is that the boundary contributions ${\rm{\Delta }}{{\mathbb{E}}}_{\mathrm{PNS}}$ and ${\rm{\Delta }}{{\mathbb{E}}}_{S}$ (Equations (49)–(50)) are invariant under the Rankine–Hugoniot relations (7)–(9). Therefore, we can determine the boundary terms using the known conditions at the PNS surface and at the immediate upstream of the shock.

At the PNS radius the inner boundary can be estimated with the quantities just below ${R}_{P}$, which include the predetermined pressure, ${P}_{P}$. The last term in Equation (49) is null since ${R}_{P}$ is fixed $({V}_{P}=0)$, and so

$$\begin{eqnarray}&&{\rm{\Delta }}{{\mathbb{E}}}_{\mathrm{PNS}}=4\pi {R}_{P}^{3}\left({P}_{P}+{\rho }_{P}{u}_{P}^{2}\right)+\displaystyle \frac{1}{2}{R}_{P}^{2}{\ddot{M}}_{P}.\end{eqnarray} \tag{ 66 }$$

The result in Equation (66) is general, but typically static pressure dominates over ram pressure at the PNS and we can approximate ${P}_{P}+{\rho }_{P}{u}_{P}^{2}\approx {P}_{P}$. In fact, in constructing the phase space, we found that accuracy can be increased further by imposing the fixed total pressure at the PNS as the sum of the actual thermal pressure at the PNS and the initial value of ${\rho }_{p}{u}_{P}^{2}$, which is a correction of a few percent. Thus, the approximation is only that ${\rho }_{p}{u}_{P}^{2}$ is not updated as the flow evolves, and this error is limited to less than 1%. In the quasi-stationary approximation $\dot{M}$ is uniform in space (although allowed to vary in time). Hence, we set ${\dot{M}}_{P}\simeq \dot{M}$ and ${\ddot{M}}_{P}\simeq \ddot{M}$ and use Equation (18).

At the shock, the quantity ${\rm{\Delta }}{{\mathbb{E}}}_{S}$ is calculated from the upstream values (index 0) instead of the shocked values (index $1$), and so

$$\begin{eqnarray}{\rm{\Delta }}{{\mathbb{E}}}_{S} & = & -\left[4\pi {R}_{S}^{3}\left({P}_{1}+{\rho }_{s}{u}_{1}^{2}\right)\right.\\ & & \left.+\displaystyle \frac{1}{2}{R}_{S}^{2}\ddot{M}({R}_{S})-\displaystyle \frac{d}{{dt}}\left(2\pi {R}_{S}^{4}{\rho }_{s}{V}_{S}\right)\right]\\ & = & -\left[4\pi {R}_{S}^{3}\left({P}_{0}+{\rho }_{0}{u}_{0}^{2}\right)\right.\\ & & \left.+\displaystyle \frac{1}{2}{R}_{S}^{2}{\ddot{M}}_{0}-\displaystyle \frac{d}{{dt}}\left(2\pi {R}_{S}^{4}{\rho }_{0}{V}_{S}\right)\right].\end{eqnarray} \tag{ 67 }$$

Finally, for accretion through the shock arriving as pressure-less free fall at a constant accretion rate, ${\dot{M}}_{0}$ (so ${\ddot{M}}_{0}=0$), and

$$\begin{eqnarray}&&{\rm{\Delta }}{{\mathbb{E}}}_{S}=-\alpha {({{GM}}_{P})}^{1/2}| {\dot{M}}_{0}| \left({R}_{s}^{1/2}-\displaystyle \frac{1}{7{\alpha }^{2}{{GM}}_{P}}\displaystyle \frac{{d}^{2}}{{{dt}}^{2}}\left({R}_{s}^{7/2}\right)\right).\end{eqnarray} \tag{ 68 }$$

In the quasi-stationary model laid out in the main text, we expressed the contributions from the boundaries through the work done by the PNS, denoted separately as ${W}_{\mathrm{PNS}}$, and the energy advected across both boundaries, combined in ${W}_{B}$. These two quantities are therefore

$$\begin{eqnarray}&&{W}_{\mathrm{PNS}}\simeq 4\pi {R}_{P}^{3}{P}_{P},\end{eqnarray} \tag{ 69 }$$

$$\begin{eqnarray}{W}_{B} & \simeq & -\alpha {({{GM}}_{P})}^{1/2}| {\dot{M}}_{0}| {R}_{S}^{1/2}\\ & & +\displaystyle \frac{| {\dot{M}}_{0}| }{7\alpha {({{GM}}_{P})}^{1/2}}\displaystyle \frac{d}{{dt}}\left[(\delta +1)\displaystyle \frac{d}{{dt}}\left({R}_{S}^{7/2}\right)\right],\end{eqnarray} \tag{ 70 }$$

where $\delta \equiv \left(\beta -1\right){({R}_{P}/{R}_{S})}^{2}$ is of the order of unity. The two equations are approximate owing to assuming a uniform accretion rate and neglecting the ram pressure at the PNS. Note that under the above assumptions we maintain the sum ${W}_{\mathrm{PNS}}+{W}_{B}\simeq {\rm{\Delta }}{{\mathbb{E}}}_{\mathrm{PNS}}+{\rm{\Delta }}{{\mathbb{E}}}_{S}$.

It is noteworthy that the energy advected across the boundaries includes an explicit dependence on the shock radius. The first term in Equation (70) is negative and restrains the shock, corresponding to the impulse of the accreting matter on the shock, while the second term is positive (assisting shock expansion) and accounts for the energy of the shocked matter and the work it does on the shock.

Summarizing the results, we conclude that the second time derivative of the spherical moment of inertia in the quasi-stationary approximation presented in the main text is given by

$$\begin{eqnarray}\displaystyle \frac{1}{2}\displaystyle \frac{{d}^{2}I}{{{dt}}^{2}} & \simeq & \tilde{K}+\tilde{U}+\tilde{{\rm{\Omega }}}+4\pi {R}_{P}^{3}{P}_{P}\\ & & -\alpha {({{GM}}_{P})}^{1/2}| {\dot{M}}_{0}| {R}_{S}^{1/2}\\ & & +\displaystyle \frac{| {\dot{M}}_{0}| }{7\alpha {({{GM}}_{P})}^{1/2}}\displaystyle \frac{d}{{dt}}\left[(\delta +1)\displaystyle \frac{d}{{dt}}\left({R}_{S}^{7/2}\right)\right].\end{eqnarray} \tag{ 71 }$$

Our final estimate for the shock acceleration in Equation (33) includes the effective energy of the accretion layer, ${{\mathbb{E}}}_{\mathrm{QS}}$, the work done by the PNS, ${W}_{\mathrm{PNS}}$, and two terms in ${W}_{B}$ that depend on the properties of the shock. Here we show that ${W}_{\mathrm{PNS}}$, which is invariably positive, generally dominates over the advected energy terms. Correspondingly, a positive ${{\mathbb{E}}}_{\mathrm{QS}}$ can serve as a sufficient (albeit not necessary) condition for a positive shock acceleration.

First, recall that ${P}_{P}$ is calculated through a stable, stationary accretion flow that corresponds to a point on the ${V}_{S}=0$ axis and is bound such that ${{\mathbb{E}}}_{\mathrm{QS}}\;\lt$ 0. Requiring ${d}^{2}I/{{dt}}^{2}=0$ along with ${{\mathbb{E}}}_{\mathrm{QS}}\;\lt$ 0 in Equation (23) implies that in the stationary solution, the pressure term at the PNS must dominate over the mass influx at the shock:

$$\begin{eqnarray}&&4\pi {R}_{P}^{3}{P}_{P}\gt \alpha {({{GM}}_{P})}^{1/2}| {\dot{M}}_{0}| {R}_{S}^{1/2}=| {\dot{M}}_{0}{R}_{S}{u}_{0}| .\end{eqnarray} \tag{ 72 }$$

This inequality holds as long as the profile is still in the oscillatory regime. It becomes invalid once the flow has reached runaway expansion (large shock radii), but by then the explosion is assumed to be well under way with a large shock velocity, and our quasi-stationary approximation breaks down in any case.

The second term originating from the flow across the boundaries in Equation (33) is $({-V}_{S}^{2}\partial \zeta /\partial {R}_{S})$, with $\zeta$ defined in Equation (34). Taking the partial derivative with respect to the shock radius yields

$$\begin{eqnarray}{V}_{S}^{2}\displaystyle \frac{\partial \zeta }{\partial {R}_{S}} & \simeq & {V}_{S}^{2}\displaystyle \frac{\left(\beta -1\right)| {\dot{M}}_{0}| }{2\alpha {({{GM}}_{P})}^{1/2}}\left[\left({R}_{S}^{2}-{R}_{P}^{2}\right)\displaystyle \frac{{R}_{S}^{-1/2}}{2}+2{R}_{S}^{3/2}\right]\\ & = & \displaystyle \frac{{V}_{S}^{2}}{{u}_{0}^{2}}\displaystyle \frac{\left(\beta -1\right)}{2}\left(\displaystyle \frac{5}{2}-\displaystyle \frac{{R}_{P}^{2}}{2{R}_{S}^{2}}\right)| {\dot{M}}_{0}{R}_{S}{u}_{0}| .\end{eqnarray} \tag{ 73 }$$

The equality is approximate because we neglected the weak dependence of $\beta$ on the shock radius. Now, for the small shock velocities we consider in the oscillatory phase, ${V}_{S}/{u}_{0}\lt 0.1$, the coefficient of $| {\dot{M}}_{0}{R}_{S}{u}_{0}|$ in the expression in Equation (73) is much smaller than unity ($\beta \approx 6-10$ and $\left(5/2-{R}_{P}^{2}/2{R}_{S}^{2}\right)\lt 5/2$). Correspondingly, ${V}_{S}^{2}\partial \zeta /\partial {R}_{S}$ is subdominant to the first shock-related term (Equation (72)) and does not change the hierarchy in which the PNS term dominates.

In this appendix we show that the extended virial theorem can be used to explain the typical oscillation period of tens of milliseconds we find in the simulations. Coming back to Equations (19), (47), and (66)–(70), we can write an oscillator-like equation for the shock radius when assuming a uniform accretion rate in the postshocked region,

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}\left(2\pi \phi {\rho }_{0}{R}_{S}^{4}{V}_{s}\right)\simeq {\mathbb{E}}+4\pi {R}_{P}^{3}{P}_{P}-| {\dot{M}}_{0}{R}_{S}{u}_{0}| ,\end{eqnarray} \tag{ 74 }$$

where $\phi \equiv (\beta -1)[1-{({R}_{P}/{R}_{S})}^{2}]$ and ${\mathbb{E}}$ is defined in Equation (48) (i.e., the quasi-stationary approximation for the energy is not necessary in this context). Equation (74) was found to be in good agreement with the shock motion for both zero and full dissociation parameters. We analyze the equation around a stationary solution for which ${V}_{S}=0$ at a shock radius ${R}_{S,0}={R}_{S}(t=0)$. For these values the right-hand side of the Equation (74) is zero.

We now assume that during the oscillations, the change in the effective energy, ${\mathbb{E}}$, is roughly proportional to the change of the inertia crossing the shock: ${\rm{\Delta }}{\mathbb{E}}=\mu | {\dot{M}}_{0}{R}_{S}{u}_{0}|$, with 0 $\lt \;\mu \;\lt$ 1. This is a lowest-order approximation, in which the change in ${\mathbb{E}}$ depends only on the shock radius (the dependence on the shock velocity is neglected). We confirmed this assumption quantitatively in the simulations. It conveys the fact that in small oscillations the accretion layer adjusts to include the material that is either added or lost as the shock moves. Correspondingly (recalling that ${W}_{\mathrm{PNS}}$ is constant in the quasi-stationary approximation),

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}\left[\displaystyle \frac{\phi }{2\alpha {({{GM}}_{P})}^{1/2}}| {\dot{M}}_{0}| {R}_{S}^{5/2}{V}_{S}\right]\\ &&\quad =-(1-\mu )\alpha {({{GM}}_{P})}^{1/2}| {\dot{M}}_{0}| \left({R}_{S}^{1/2}-{R}_{S,0}^{1/2}\right).\end{eqnarray} \tag{ 75 }$$

Considering that the shock radius is significantly larger than the PNS radius, $1-{({R}_{P}/{R}_{S})}^{2}\;\approx$ 1, and combined with the fact that the compression ratio depends weakly on the shock radius, we can treat $\phi$ as roughly constant during the oscillations. Equation (75) can now be rewritten as

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}\ {R}_{S}^{7/2}}{{{dt}}^{2}}=-\displaystyle \frac{7(1-\mu ){\alpha }^{2}{{GM}}_{P}}{\phi }\left({R}_{S}^{1/2}-{R}_{S,0}^{1/2}\right).\end{eqnarray} \tag{ 76 }$$

It is noteworthy that the mass accretion rate ${\dot{M}}_{0}$ has canceled out of the equation. For small oscillations, this is a harmonic oscillator equation with a time period of

$$\begin{eqnarray}&&{T}_{O}=2\pi \displaystyle \frac{{R}_{S,0}}{| {u}_{0}| }\sqrt{\displaystyle \frac{\phi }{1-\mu }.}\end{eqnarray} \tag{ 77 }$$

For a compression ratio in the range $\beta =6-10$ and ${R}_{P}/{R}_{S}\lt 0.5$, the value of $\sqrt{\phi }$ is confined to values of 2–3. Unless $\mu$ is very close to unity, we conclude that the oscillation period should be (10–20)${R}_{S,0}/| {u}_{0}|$, or 25–50 ms. This is indeed in good agreement with the results of the simulations. We also recover the general trend that the period should be weakly dependent on the mass accretion rate and on the neutrino luminosity, as their effect is limited to the finer details of $\sqrt{\phi /(1-\mu )}$.

Finally, we note that in the approximation above, the oscillations are unconditionally stable. This is due to the neglect of a $\partial ({\rm{\Delta }}{\mathbb{E}})/\partial {V}_{S}$ term in the derivation. This partial derivative is directly related to the $\partial {a}_{S}/\partial {V}_{S}$ derivative discussed in the main text, which determines whether the shock oscillations around the stationary solutions will damp or grow to a runaway expansion.