DYNAMICS OF A SPHERICAL ACCRETION SHOCK WITH NEUTRINO HEATING AND ALPHA-PARTICLE RECOMBINATION

Although tremendous progress has been made on the mechanism of core-collapse supernovae in recent years, we still do not have a clear picture of a robust path to an explosion in stars that form iron cores. Heating by the absorption of electron-type neutrinos significantly modifies the settling flow below the bounce shock but—in spite of early positive results (Bethe & Wilson 1985)—explosions are obtained in spherical collapse calculations only if the progenitor star is lighter than about 10–12 M_☉ (Kitaura et al. 2006). More massive stars fail to explode in spherical symmetry (Liebendörfer et al. 2001).

Two-dimensional collapse calculations show strong deformations of the shock and convective motions below it (Burrows et al. 1995, 2006, 2007; Janka & Mueller 1996; Buras et al. 2006a, 2006b; Marek & Janka 2009). It has long been recognized that convection increases the residency time of settling material in the zone of strong neutrino heating (Herant et al. 1992). It is also becoming clear that multidimensional explosions require the assembly of a smaller amount of material with positive energy, but the details of how this happens remain murky.

An early treatment of shock breakout by Bethe (1997) focused on the strong gradient in the ram pressure of the infalling material, but implicitly assumed that the shocked material had already gained positive energy. If large-scale density inhomogeneities are present below the shock, they will trigger a finite-amplitude, dipolar instability, thereby allowing accretion to continue simultaneously with the expansion of positive-energy fluid (Thompson 2000). The accretion shock is also capable of a dipolar oscillation which leads, above a critical amplitude, to a bifurcation between freshly infalling material and material shocked at earlier times (Blondin et al. 2003). Such an oscillation is easily excited in a spherical flow composed of a zero-energy, polytropic fluid (Blondin & Mezzacappa 2006) via a linear feedback between ingoing vortex and entropy waves and outgoing sound waves (Foglizzo et al. 2007). It can also be excited indirectly by neutrino heating, which if sufficiently strong will drive large-scale buoyant motions below the shock (Herant et al. 1994; Foglizzo et al. 2006). For relatively weak heating, the dipolar oscillation can decrease the damping effect of advection and trigger convective motions that would otherwise be suppressed (Foglizzo et al. 2006; Scheck et al. 2008).

A significant sink of thermal energy in the accretion flow arises from the dissociation of heavy nuclei. The heavy elements that flow through the shock are broken up into α-particles and nucleons when exposed to the high temperature (>1 MeV) of the post-shock region. The Bernoulli parameter b of the shocked fluid then becomes substantially negative. A significant fraction of this dissociation energy can be recovered if nucleons recombine into α-particles (Bethe 1996). But for this to occur, a decrease in the temperature is required and thus the shock must expand significantly beyond the radius at which it typically stalls (∼100–150 km).

One of the primary goals here, and in a previous paper (Fernández & Thompson 2009; hereafter Paper I), is to gauge the relative importance of these effects in setting the stage for a successful explosion. The persistence and amplitude of a dipolar oscillation fully can only be reliably measured in fully three-dimensional simulations (there are preliminary indications that it is less prominent in three dimensions than in axial symmetry; Iwakami et al. 2008). On the other hand, the interplay between α-particle recombination and hydrodynamical instabilities has received little attention. Although recombination is certainly present in previous numerical studies which employ finite-temperature equations of state (EOSs), it should be kept in mind that considerable uncertainties in the EOS remain at supranuclear densities. A softening or hardening of the EOS feeds back on the position of the shock for a given pre-collapse stellar model (Marek & Janka 2009). The parametric study of the critical neutrino luminosity by Murphy & Burrows (2008) is based on a single EOS; they obtain explosions in which the shock seems to break out from nearly the same radial position at ∼250–300 km. Variations in the density profile of the progenitor star will similarly modify the position of the shock, the concentration of α-particles below it, and the critical neutrino luminosity for an explosion.

In this paper, we study the interplay between nonspherical shock oscillations, neutrino heating, and α-particle recombination, when the heating rate is pushed high enough for an explosion to occur. Our focus is on the stalled shock phase, between ∼100 ms and 1 s after bounce. In a one-dimensional calculation, the accretion flow reaches a quasi-steady state during this interval, and the shock gradually recedes (e.g., Liebendörfer et al. 2001; Buras et al. 2006a).

Our approach is to introduce EOSs of increasing complexity into one- and two-dimensional, time-dependent hydrodynamic simulations. To this end, we use the code FLASH2.5 (Fryxell et al. 2000), which is well tested in problems involving nuclear energy release in compressible fluids (Calder et al. 2002). We adopt a steady-state model as our initial condition, and a constant mass accretion rate, neutrino luminosity, and fixed inner boundary. The steady-state approximation to the stalled shock phase was first introduced by Burrows & Goshy (1993), and has recently been used by Ohnishi et al. (2006) to study the nonlinear development of the shocked flow with a semi-realistic equation of state and neutrino heating. In Paper I, we modeled the accretion flow as a polytropic fluid, from which a fixed dissociation energy is removed immediately below the shock, and allowed for neutrino cooling but not heating.

Here we generalize this model to allow both for heating, and for nuclear statistical equilibrium (NSE) between neutrons, protons, and α-particles. Heating is introduced in a simple, parameterized way, without any attempt at simulating neutrino transport. The nuclear abundances are calculated as a function of pressure p and density ρ, using a complete finite-temperature, partially degenerate EOS. Our model for the shocked material retains one significant simplification from Paper I: we do not allow the electron fraction Y_e to vary with position below the accretion shock. Here there are two competing effects: electron captures tend to reduce Y_e, whereas absorption of ν_e and $\bar{\nu }_e$ tends to drive high-entropy material below the shock toward Y_e ≃ 0.5. Since we are interested especially in the dynamics of this high-entropy material, we set Y_e = 0.5 when evaluating the α-particle abundance. To obtain a realistic density profile, we continue to approximate the internal energy of the fluid as that of a polytropic fluid with a fixed index $\gamma = {4\over 3}$. In reality, the equation of state between the neutrinosphere and the shock depends in a complicated way on the degeneracy of the electrons and the effects of electron captures. The consequences of introducing these additional degrees of freedom will be examined in future work.

In spite of these simplifications, our results already show many similarities with more elaborate collapse calculations. Spherical explosions are due to a global instability resembling the one-dimensional standing accretion shock instability (SASI), but modified by heating. As in Paper I, we find that the period of the ℓ = 0 mode remains close to twice the post-shock advection time. Strong deformations of the shock in two-dimensional simulations are driven by material with positive Bernoulli parameter, which generally resides outside the radius r_α where the gravitational binding energy of an α-particle is equal to its nuclear binding energy. The recombination of α-particles can play a major role in creating this positive-energy material, but for this to happen the shock must be pushed beyond ∼200 km from the neutronized core.

In this paper, we consider only neutrino heating as the impetus for the initial expansion of the shock, rather than more exotic effects such as rotation or magnetic fields. We find that the critical heating rate is a strong function of the initial position of the shock with respect to r_α, which implies that a much higher neutrino luminosity is needed to revive a shock that stalls well inside r_α. The difference in the critical heating rate between one-dimensional and two-dimensional simulations also depends on the size of the shock, and thence on the structure of the forming neutron star.

We find that the shock develops a dipolar oscillation with a large amplitude only when the heating rate is also high enough to trigger a strong convective instability. We therefore surmise that buoyant motions driven by neutrino heating play a major role in driving the dipolar oscillations that are seen in more complete simulations of core collapse. Some evidence is found that acoustic wave emission by the convective motions also plays a role. We investigate the possibility of a heat engine within the gain layer (the region with a net excess of neutrino heating over cooling). We find that most of the heat deposition by neutrinos is concentrated in lateral flows at the base of the prominent convective cells. At the threshold for an explosion, neutrino heating plays a key role in pushing material to positive Bernoulli parameter, but only if the shock starts well inside r_α.

The plan of the paper is as follows. Section 2 presents our numerical setup, treatment of heating, cooling, and nuclear dissociation, and outlines the sequences of models. Sections 3 and 4 show results from one-dimensional and two-dimensional simulations, respectively. We focus on the relative effectiveness of neutrino heating and α-particle recombination in driving an explosion, and the relation between Bernoulli parameter and large-scale deformations of the shock. The critical heating rate for an explosion is analyzed in Section 5, and the competition between advective-acoustic feedback and convective instability is discussed in Section 6. We summarize our findings in Section 7. The appendices contain details about our EOS and the numerical setup.

As in Paper I, the initial configuration is a steady, spherically symmetric flow onto a gravitating point mass M. The flow contains a standing shock wave, and the settling flow below the shock cools radiatively in a narrow layer outside the inner boundary of the simulation volume. The space of such models is labeled basically by three parameters: accretion rate $\dot{M}$, luminosity L_ν in electron neutrinos and anti-neutrinos, and the radius r_* of the base of the settling flow, which corresponds roughly to the neutrinosphere radius. The mass M of the collapsed material represents a fourth parameter, but it covers a narrower range than the other three. The infalling material is significantly de-leptonized before hitting the shock only in the first 50 ms or so of the collapse (e.g., Liebendörfer et al. 2001).

We explore a two-dimensional surface through this three-dimensional parameter space by (1) fixing the ratio of r_* to the initial shock radius r_s0 in the absence of heating (r_*/r_s0 = 0.4); (2) allowing r_s0 to vary with respect to an appropriately chosen physical radius; and then (3) increasing the level of heating until an explosion is uncovered. In the full problem, the shock radius at zero heating is a unique function of $\dot{M}$ and r_*, with a small additional dependence on M and the composition of the flow outside the shock (Houck & Chevalier 1992). The secular cooling of the collapsed core forces a gradual decrease in r_*, and $\dot{M}$ also varies with time and with the progenitor model.

Given the important role that α-particle recombination plays in the final stages of an explosion, we implement (2) by referencing r_s0 to the radius where the gravitational binding energy of an α-particle equals its nuclear binding energy,

$$\begin{equation} r_\alpha = \frac{GMm_\alpha }{Q_\alpha } \simeq 254 M_{1.3}\;{\rm km}. \end{equation} \tag{ 1 }$$

Here Q_α ≃ 28.30 MeV is the energy needed to break up an α-particle into 2n and 2p, m_α the mass of an α-particle, and M_1.3 = M/(1.3 M_☉). Choice (1) allows us to consider models that have, implicitly, both a range of physical values of r_* and a range of $\dot{M}$. It is, of course, made partly for computational simplicity (the limited size of the computational domain) and also to facilitate a comparison between models that have different values of r_s0/r_α. Nuclear dissociation is taken into account either by removing a fixed specific energy ε right below the shock, or by enforcing NSE between n, p, and α throughout the settling flow. Once this choice is made, the normalization of the cooling function is adjusted to give r_*/r_s0 = 0.4. The heating rate remains freely adjustable thereafter.

We adopt this simplification because we do not intend to find the precise value of the critical neutrino luminosity, but instead to probe the behavior of the system around this critical point, whatever its absolute value.

We now describe the key components of this model in more detail and explain the setup of the hydrodynamic calculations. As in Paper I, the time evolution is carried out using the second-order, Godunov-type, adaptive-mesh-refinement code FLASH2.5 (Fryxell et al. 2000).

2.1. Initial Conditions

The introduction of heating causes a change in the structure of the initial flow configuration. The radius r_s of the shock in the time-independent solution to the flow equations increases with heating rate; that is, r_s ⩾ r_s0. The material above the shock is weakly bound to the protoneutron star, and in practice can be taken to have a zero Bernoulli parameter

$$\begin{equation} b = \frac{1}{2}v^2 + \frac{\gamma }{\gamma -1}\frac{p}{\rho }- \frac{GM}{r}. \end{equation} \tag{ 2 }$$

Here v is the total fluid velocity, which is radial in the initial condition, p is the pressure, ρ is the mass density, and G is Newton's constant. The flow upstream of the shock is adiabatic and has a Mach number $\mathcal {M}_1=5$ at a radius r = r_s0.

The composition of the fluid is very different upstream and downstream of the shock. Changes in internal energy due to nuclear dissociation and recombination are taken into account using the two models described in Section 2.1.1. For the internal energy density of the fluid e, we continue to use the polytropic relation e = p/(γ − 1); hence, the second term on the right-hand side of Equation (2). Because we are not explicitly including changes in electron fraction due to weak processes, we keep γ = 4/3 for both models of nuclear dissociation. This largely determines the density profile inside a radius ${\sim}{1\over 2}r_\alpha$, where the α-particle abundance is very low.

The upstream and downstream flow profiles are connected through the Rankine–Hugoniot jump conditions, which are modified so as to allow a decrement ε in b across the shock. The resulting compression factor is (Paper I)

$$\begin{eqnarray} \kappa \equiv \frac{\rho _2}{\rho _1} & = & (\gamma +1)\Bigg [ \left(\gamma + \mathcal {M}_1^{-2}\right) \nonumber \\ & &\left.-\,\sqrt{ \left(1-\mathcal {M}_1^{-2}\right)^2 + (\gamma ^2-1)\frac{2\varepsilon }{v_1^2}} \quad \right]^{-1}, \end{eqnarray} \tag{ 3 }$$

which reduces to κ → (γ + 1)/(γ − 1) for $\mathcal {M}_1\rightarrow \infty$ and ε = 0. Throughout this paper, the specific nuclear dissociation energy ε is defined to be positive. The subscripts 1 and 2 denote upstream and downstream variables, respectively.

All flow variables are made dimensionless by scaling radii to r_s0, velocities to v_ff 0 = (2GM/r_s0)^1/2, timescales to t_ff 0 = r_s0/v_ff 0, and densities to the upstream density at r = r_s0, ρ₁(r_s0) (Equation (B13)). See Paper I for further details. Throughout the paper, we denote the average of a function F(X, ...) over some variable X by 〈F〉_X.

2.1.1. Nuclear Dissociation

We model nuclear dissociation in two ways. First, we remove a fixed specific energy ε right below the shock, as done in Paper I. This represents the prompt and complete breakup of whatever heavy nuclei are present in the upstream flow. The main limitation of this approximation is that the dissociation energy does not change with the radius (or inclination) of the shock. The main advantage is simplicity: ε is independent of any dimensional parameters and can be expressed as a fraction of v²_ff 0.

We also use a more accurate dissociation model which allows for NSE between α-particles and nucleons.³ During the stalled shock phase of core-collapse supernovae, the shock sits at r ∼ 100–200 km, with a post-shock temperature T>1 MeV and density ρ ≳ 10⁹ g cm⁻³. In these conditions, the heavy nuclei flowing through the shock are broken up into α, p, and n.

A range of isotopes are present in the iron core of a massive star as well as in nuclear burning shells (Woosley et al. 2002), but since the binding energy per nucleon varies only by ∼10% we simply assume a single type of nucleus in the upstream flow. We focus here on the later stages of the stalled shock phase, during which the oxygen shell is accreted. An energy Q_O ≃ 14.44 MeV must be injected to dissociate an ¹⁶O nucleus into four α-particles (Audi et al. 2003), which corresponds to the specific dissociation energy

$$\begin{equation} \varepsilon _{\rm O} = \frac{Q_{\rm O}}{m_{\rm O}} \simeq 0.038 M_{1.3}^{-1} \left({r\over 150\;{\rm km}}\right)\,v_{\rm ff}^2(r). \end{equation} \tag{ 4 }$$

Here m_O ≃ 16m_u is the mass of an oxygen nucleus, with m_u being the atomic mass unit. The smallness of this number indicates that little oxygen survives in the post-shock flow, and so we set the equilibrium mass fraction of oxygen to zero below the shock, X^eq_O = 0. The binding energy of an α-particle is of course much larger, giving

$$\begin{equation} \varepsilon _\alpha = \frac{Q_\alpha }{m_\alpha } \simeq 0.297 M_{1.3}^{-1}\left({r\over 150\;{\rm km}}\right) \,v_{\rm ff}^2(r). \end{equation} \tag{ 5 }$$

We find that α-particles appear in significant numbers only at relatively large radii (≳0.5 r_α) and in material that has either (1) been significantly heated by electron neutrinos closer to the neutrinosphere or (2) been freshly shocked outside r_α. The electrons are only mildly degenerate in material that has a high entropy and α-particle content, so that neutrino heating drives Y_e close to ∼0.5 (or even slightly above: see; e.g., Buras et al. 2006b). We therefore set Y_e = 0.5 in the Saha equation that determines the equilibrium mass fractions X^eq_n, X^eq_p, and X^eq_α = 1 − X^eq_n − X^eq_p. These quantities are tabulated as functions of p and ρ using an ideal, finite-temperature and partially degenerate equation of state for electrons and nucleons; see Appendix A for details. Specific choices must then be made for the parameters r_s0, M, and $\dot{M}$; we generally take M = 1.3 M_☉ and $\dot{M} = 0.3\,M_\odot$ s⁻¹, but allow r_s0 to vary. An investigation of how changes in Y_e feed back onto the formation of α-particles is left for future work.

A specific energy

$$\begin{equation} e_\mathrm{nuc} = - X_{\rm O} (\varepsilon _{\rm O} + \varepsilon _\alpha) - \left(X_\alpha - X_\alpha ^\mathrm{eq}[\rho,p] \right) \varepsilon _\alpha \end{equation} \tag{ 6 }$$

is either released or absorbed within a single time step (it can be of either sign). Here X_O is nonvanishing only for fluid elements that have just passed across the shock, and we have set X^eq_O = 0. The quantity (6) is introduced as an energy source term in FLASH, and from it one readily obtains a rate of release of nuclear binding energy per unit mass,

$$\begin{equation} {{{d}}e_{\rm nuc}\over {{d}}t} \equiv {e_{\rm nuc}\over \Delta t}, \end{equation} \tag{ 7 }$$

where Δt is the simulation time step.

In the initial condition, the dissociation energy at the shock is obtained from Equation (6) using X_O = 1 and X_α = 0 upstream of the shock:

$$\begin{equation} \varepsilon (t=0) = \varepsilon _{\rm O} + \left(1 - X_\alpha ^\mathrm{eq}[\rho _2,p_2]\right)\varepsilon _\alpha. \end{equation} \tag{ 8 }$$

Figure 1 shows how ε(t = 0) and X^eq_α depend on the shock radius r_s0, for upstream flows composed⁴ of pure ¹⁶O and ⁵⁶Fe, and for different values of $\dot{M}$. The dissociation energy is approximately constant inside ∼75 km, where the downstream flow is composed of free nucleons, but decreases at greater distances, remaining ∼40% of the gravitational binding energy at the shock. The mass fraction of α-particles reaches 50% at r = 150–175 km, with a weak dependence on $\dot{M}$.

Zoom In Zoom Out Reset image size

2.1.2. Heating and Cooling in the Post-shock Flow

To allow direct comparison with our previous results, we employ a cooling rate per unit volume of the form

$$\begin{equation} {\LM\!}_{C} = C p^a \rho ^{b-a}, \end{equation} \tag{ 9 }$$

with a = 1.5, b = 2.5, and C a normalization constant. As in Paper I, we include a Gaussian entropy cutoff to prevent runaway cooling. The exponents in Equation (9) represent cooling dominated by the capture of relativistic, nondegenerate electrons and positrons on free nucleons (e.g., Bethe 1990). Inside the radius where the electrons become strongly degenerate, and α-particles are largely absent, one has ${\LM\!}_C \propto p_e^{3/2} n_p \propto (Y_e \rho)^3$. This gives essentially the same dependence of ${\LM\!}_C$ on r as Equation (9) when Y_e= constant and $\gamma = {4\over 3}$ (corresponding to ρ ∝ r⁻³ in a nearly adiabatic settling flow). In more realistic collapse calculations, Y_e grows with radius between the neutrinosphere and the shock, but ρ tends to decrease more rapidly than ∼r⁻³ (e.g., Buras et al. 2006a). Our chosen form for the cooling function results in a slightly wider gain region and, therefore, a slightly lower critical heating rate for an explosion. The bulk of the cooling occurs in a narrow layer close to the accretor at r = r_*, and the accreted material accumulates in the first few computational cells adjacent to the inner boundary without a major effect on the rest of the flow.

We model neutrino heating as a local energy generation rate per unit volume of the form

$$\begin{equation} {\LM\!}_H = H (1-X_\alpha) \rho /r^2. \end{equation} \tag{ 10 }$$

The normalization constant H measures the strength of the heating. The factor (1 − X_α) accounts for the fact that the cross section for neutrino absorption by α-particles is much smaller than that for free nucleons (Bethe 1990). For simplicity, we do not include the flux factor due to the transition between diffusion and free-streaming. Our focus here is on the nature of the instabilities occurring in the flow near the threshold for an explosion, and we do not attempt a numerical evaluation of the critical heating rate.

An additional energy source term arises from the change in the equilibrium fraction of α-particles as they are advected in the steady-state initial solution. The instantaneous adjustment of X_α to its equilibrium value, combined with Equation (6), yields an energy generation rate per unit volume

$$\begin{equation} {\LM\!}_\alpha = \rho v \varepsilon _\alpha \frac{{{d}}X^\mathrm{eq}_\alpha }{{{d}}r} = \rho v \varepsilon _\alpha \left[ \frac{\partial {X_\alpha ^{\rm eq}}}{\partial \rho }\frac{{{d}}\rho }{{{d}}r} + \frac{\partial {X_\alpha ^{\rm eq}}}{\partial p}\frac{{{d}}p}{{{d}}r}\right]. \end{equation} \tag{ 11 }$$

This energy generation rate is negative, as the temperature increases inward and thus the α-particle fraction decreases with decreasing radius (v is negative).

2.1.3. Numerical Setup

In our time-dependent calculations, we use one-dimensional and two-dimensional axisymmetric spherical coordinates with baseline resolution Δr_base = r_s0/320 and Δθ_base = π/192, with one extra level of mesh refinement inside r = r_* +  0.1(r_s0 − r_*) to better resolve the steep density gradient that arises in the cooling layer. We do not employ a hybrid Riemann solver because we do not see the appearance of the odd–even decoupling instability (Quirk 1994).

We employ a reflecting inner boundary condition at r = r_* for the sake of simplicity; we do not attempt to model the protoneutron star (as done by Murphy & Burrows 2008) or its contraction through a moving inner boundary (as done by Scheck et al. (2006, 2008)). The outer boundary condition is kept fixed at r = 7r_s0, and is set by the upstream flow at that position.

To trigger convection below the shock, we introduce random cell-to-cell velocity perturbations in v_r and v_θ at t = 0, with an amplitude 1% of the steady-state radial velocity. To study the interplay between shock oscillations and convection, we also drop overdense shells with a given Legendre index ℓ, as done in Paper I, without random velocity perturbations.

In order to track the residency time of the fluid in the gain region, we assign a scalar to each spherically symmetric mass shell in the upstream flow. This scalar is passively advected by FLASH2.5. Through this technique, we are able to assign a "fluid" time to each element in the domain, corresponding to the time at which the mass shell would cross the instantaneous angle-averaged shock position if advected from the outer boundary at the upstream velocity:

$$\begin{equation} t_F = t_{\rm OB} + \int _{\langle r_s(t)\rangle _\theta }^{r_{\rm OB}} \frac{{{d}}r}{|v_r|}. \end{equation} \tag{ 12 }$$

Here t_OB is the time at which the fluid enters through the outer radial boundary at r = r_OB, and 〈r_s(t)〉_θ is the angle-averaged shock position. Initially, t_OB = 0 and all the fluid below the shock is set to t_F = 0. This prescription works well for statistical studies (Section 4.2), tracing large-scale fluid patches, despite some inevitable turbulent mixing on small scales.

An explosion is defined as either (1) a collision between the shock and the outer boundary of the simulation volume (r = 7r_s0) within 1000t_ff0 of the start of the simulation; or (2) in the special case of the one-dimensional constant-ε models, a transient expansion that breaks a quasi-steady pattern within the same time frame. Even in the one-dimensional simulations, very small changes in heating rate can lead to dramatic changes in shock behavior, and so this definition of explosion is good enough for our purposes.

2.2. Model Sequences

We choose six sequences of models, each with a range of heating parameters H ⩾ 0, and each evolved both in spherical (one-dimensional) and axial (two-dimensional) symmetry. Their parameters are summarized in Table 1. In each sequence, the normalization of the cooling function is chosen so that r_*/r_s = 0.4 at zero heating. Three sequences have a constant dissociation energy, which take the values ε/v²_ff 0 = {0.1, 0.15, 0.2}. The other three sequences assume NSE below the shock, and have shock radii r_s0 = {50, 75, 125} km at zero heating. This means that the physical value of the cooling radius also takes on different values, namely {20, 30, 50} km. In effect, our models are probing different sizes for the neutrinosphere, and different times following the collapse. The other parameters in the NSE models are M = 1.3 M_☉ and $\dot{M}=0.3\;M_\odot$ s⁻¹.

Table 1. Sample Configurations

ε/v2ff 0	Hv−3ff 0rs0−1	rs/rs0	ε/v21	κ	χ
0.1	0	1.00	0.10	7.3	0
	8.00 E−3	1.27	0.13	7.7	4.5
	1.48E−2a	2.57	0.31	11.0	22
0.15	0	1.00	0.15	8.6	0
	7.00E−3	1.29	0.20	9.6	9.0
	1.17E−2a	2.34	0.41	18.9	40
0.2	0	1.00	0.20	10.1	0
	5.50E−3	1.30	0.27	12.5	19
	8.38E−3a	2.08	0.47	34.9	74
rs0 (km)	Hv−3ff 0rs0−1	rs/rs0	ε(t = 0)/v21	κ	χ	Xeqα(rs)
50	0	1.00	0.11	7.6	0	5.5E-6
	8.00E−3	1.29	0.15	8.1	5.5	6.1E-5
	1.43E−2a	3.01	0.26	8.6	27	0.43
75	0	1.00	0.17	9.0	0	4.3E-4
	6.50E−3	1.30	0.22	10.1	11	4.5E-2
	1.15E−2a	3.61	0.21	6.8	51	0.83
125	0	1.00	0.21	10.7	0	0.26
	3.50E−3	1.33	0.21	9.9	22	0.51
	7.28E−3a	3.99	0.18	6.1	120	0.99

Note. ^aMaximum heating rate for a steady flow solution (Burrows & Goshy 1993).

Download table as:  ASCIITypeset image

Table 1 samples some properties of a few models from each sequence: one with zero heating, another with H close to the critical value for an explosion, and a third with the largest heating parameter that will allow a steady solution. Note that the shock starts out at ∼1.3 r_s0 in the time-independent, spherical flow solution, and quickly saturates at ∼(1.8–2)r_s0 in the two-dimensional models with heating just below threshold for an explosion. The quantity ε/v²₁ references the dissociation energy to (twice) the kinetic energy of the upstream flow, and is the key free parameter determining the compression rate κ across the shock (Equation (3)).

When examining how the prescription for nuclear dissociation influences the results, we will focus on the ε = 0.15v²_ff 0 sequence and the NSE sequence with r_s0 = 75 km, which have similar initial density profiles (due to the low initial α-particle abundance in the NSE model).

The six initial models at zero heating are shown in Figure 2(a). Panel (b) shows the sequence of initial models with r_s0 = 75 km and a range of heating parameters. The model with H = 0.007v³_ff0r_s0 is close to the threshold for an explosion, while the one with H = 0.009v³_ff0r_s0 is well above threshold. At higher values of H, cooling by α-particle dissociation (Equation (11)) can be significant in a layer below the shock, causing the density profile to steepen slightly.

Zoom In Zoom Out Reset image size

Figure 2(c) shows how our constant-γ, ideal gas approximation to the internal energy of the flow compares with the full EOS containing finite-temperature and partially degenerate electrons (see Appendix B for details). The curves labeled "+α" include our prescription for heating/cooling by α-particle recombination/dissociation, and those labeled "−α" do not. We show the sequence with the largest shock radius (r_s0 = 125 km) so that NSE allows some α's to be present. The neglect of electron captures below the shock results in an adiabatic index between 4/3 and 5/3 in the zone where α-particles are absent. This causes the EOS to stiffen, so that the density profile is well approximated by an ideal gas with γ ≃ 1.48 at zero heating. Adding in heating tends to flatten the density profile even more, and with γ = 1.48 it would be much flatter than is typically seen in a realistic core collapse model. Hence we choose an EOS with γ = 4/3.

3.1. Shock Oscillations and Transition to Explosion

An explosion in spherical symmetry involves the development of an unstable ℓ = 0 SASI mode. We showed in Paper I that, in the absence of neutrino heating, the period of this mode is essentially twice the post-shock advection time. As heating is introduced into the flow, we find that this relation is maintained. The ℓ = 0 mode is damped until the heating rate is pushed above a critical value, which we now discuss.

It should be emphasized that this critical heating rate is generally lower than that defined by Burrows & Goshy (1993), which marked the disappearance of a steady, spherically symmetric solution to the flow equations. Large amplitude one-dimensional shock oscillations have been witnessed near the threshold for explosion in calculations by Ohnishi et al. (2006) and Murphy & Burrows (2008). Both calculations employed a realistic EOS, but like us included neutrino heating as a local source term in the energy equation. Oscillations have also been seen by Buras et al. (2006b) in more elaborate calculations with Boltzmann neutrino transport.

The origin of the one-dimensional SASI oscillation can be briefly summarized as follows. An initial outward shock displacement generates an entropy perturbation, which is negative for γ ≲ 5/3. This entropy perturbation is advected down to the cooling layer, where it causes, at constant ambient pressure, an increase in the cooling rate, $\delta {\LM\!}_C/{\LM\!}_C = -[(\gamma -1)/\gamma ](b-a) \delta S > 0$. The resulting negative pressure perturbation is rapidly communicated to the shock, which recedes and generates an entropy perturbation of the opposing sign. One more iteration results in a shock displacement of the same sign as the initial displacement, and allows the cycle to close. The duration of the ℓ = 0 mode is therefore nearly twice the advection time from the shock to the cooling layer,

$$\begin{equation} {2\pi \over \omega _{\rm osc}} \simeq 2\int _{r_*}^{r_s} {{{d}}r\over |v_r|}. \end{equation} \tag{ 13 }$$

The cycle is stable for γ = 4/3 and r_*/r_s0 = 0.4.

When heating is added, the density profile flattens. Increasing γ has the same effect, and has been found to push up the growth rate of linear SASI modes (e.g., Paper I). There is a critical heating rate for which the damping effect of the one-dimensional SASI is neutralized and there is no net growth. We find that, once the heating rate exceeds this critical value, the system always explodes.

We therefore define the critical heating rate in our one-dimensional simulations to be the minimum heating rate for growing shock oscillations.⁵ Figure 3 shows real and imaginary eigenfrequencies as a function of heating rate for our one-dimensional initial configurations with constant ε. The curves were obtained by solving the differential system of Foglizzo et al. (2007), modified to account for a constant rate of nuclear dissociation (Paper I) as well as incorporating the heating function in Equation (10). The runs marked by stars explode within a time 1000t_ff0, and so require a small, but finite, positive growth rate.

Zoom In Zoom Out Reset image size

In an exploding run, the expansions become longer and contractions shorter as the shock oscillation develops a large amplitude. Eventually the accretion flow is halted during a contraction. This marks the point of explosion, beyond which the feedback between the shock and the cooling layer is broken. Material then tends to pile up in the gain region, is further heated, and more material reverses direction. The net effect is to push the shock outward. Movie 1 in the online material illustrates this chain of events.

3.2. Effects of Alpha-particle Recombination

Shock breakout is controlled by the build-up of positive energy fluid downstream of the shock, and therefore is sensitive to the density profile below the shock. Heating by neutrinos is concentrated fairly close to the protoneutron star, inside a distance ∼(2–3)r_*. Heating by α-particle recombination is concentrated at a greater distance ∼r_α (Equation (1)), but still can reach a comparable amplitude.

The dependence of shock breakout on heating rate is displayed in Figure 4 for two accretion models and several values of H close to H_cr (see Table 1). The initial expansion of the shock during the explosion phase is very similar for models with constant ε and with NSE in the shocked fluid. However, the time evolution bifurcates near the radius r_α.

Zoom In Zoom Out Reset image size

Figure 5 shows successive profiles of the shocked flow in the exploding run with H = 1.09H_cr and r_s0 = 75 km. The α-particle fraction approaches unity as the shock reaches the radius r_α. The second panel shows the specific nuclear energy generation rate (Equation (7)) normalized to the adiabatic rate of change of the enthalpy,

$$\begin{equation} w_{\rm ad} = \frac{1}{\rho }\frac{{{d}}p}{{{d}}t} = -c_s^2 \nabla \cdot \mathbf {v}. \end{equation} \tag{ 14 }$$

Here c_s = (γp/ρ)^1/2 is the sound speed. The third panel compares the amplitude and distribution of neutrino and recombination heating, and the bottom panel plots the radial velocity in the post-shock region.

Zoom In Zoom Out Reset image size

We can summarize this behavior as follows: during the initial expansion phase, fluid below the shock continues to move inward, and the dissociation of α-particles removes energy from the flow (as expected from Equation (11)). Some fluid behind the shock begins to move outward around 300t_ff0, but nuclear dissociation still causes a net loss of internal energy. However, the recombination of α-particles sets in above r_α, especially in regions where X_α ≲ 0.5. By the time the shock hits the outer boundary, de_nuc/dt exceeds one-half of |w_ad|.

The dependence of the density contrast κ (Equation (3)) on radius also has an influence on the details of breakout. When the dissociation energy ε is held fixed, κ increases toward larger radius. This has the effect of creating a dense layer of fluid below the shock when the shock has reached a radius where ε ∼ v²_ff/2. In spherical symmetry, the breakout of the shock is then impeded by this dense layer, which cannot exchange position with the lighter material below it. It can happen that the energy in the expanding region is no longer able to sustain the heavier material above, and the shock collapses, as shown in Figure 4 for the constant-ε run with H = 1.08H_cr. This obstruction is avoided when statistical equilibrium between n, p, and α is maintained below the shock, because ε/v²₁ and κ both decrease gradually as the shock expands to distances much larger than r_s0 (Figure 1). This limit to the shock expansion does not occur in two dimensions, as the superposition of dense fluid over lighter fluid is Rayleigh–Taylor unstable on the dynamical time t_ff0.

Extending the flow calculation to two dimensions reveals some subtle patterns of behavior. The time evolution of the shock is shown in the left panel of Figure 6 for a range of heating rates near the threshold for explosion. In contrast with the one-dimensional runs, the breakout of the shock looks similar in models with constant dissociation energy and with NSE between n, p, and α below the shock. Both types of models are subject to buoyancy-driven instabilities, which allow cold material below the shock to interchange position with hotter material within the gain region. As a result, the shock is highly asymmetric at breakout in both cases. In Section 5, we compare the critical heating rate for explosion in one- and two- dimensional simulations, and examine how it is influenced by α-particle recombination.

Zoom In Zoom Out Reset image size

Around the threshold for explosion, all of our runs develop vigorous convective motions before the SASI has a chance to undergo even a few oscillations. At high heating rates, we find that the convective instability is driven by the negative entropy gradient within the layer of maximal neutrino heating. In nonexploding runs, the shock settles to a quasi-equilibrium state with oscillations taking place over a range of angular (Legendre) index ℓ, as previously seen by Ohnishi et al. (2006), Scheck et al. (2008), and Murphy & Burrows (2008). The amplitude of the ℓ = 1 and 2 modes remains small until the heating parameter H has begun to exceed about one-half the critical value for an explosion. The competition between SASI growth and convective instability is examined in detail in Section 6.

One gains considerable insight into the mechanism driving shock breakout by examining the distribution of Bernoulli parameter (Equation (2)) in the shocked fluid. We first consider the NSE runs with r_s0 = 50 km and 125 km, with the heating parameter H just above the threshold for an explosion. Two snapshots from each of these runs are shown in Figure 7. In the first case, the initial equilibrium shock radius is only ∼r_α/4 km, and α-particles are essentially absent below the shock. In the second, the shock starts at ∼2r_α/3 and X_α ∼ 0.5 initially in the post-shock flow.

Zoom In Zoom Out Reset image size

Large deformations of the outer shock are caused by convective plumes that carry positive energy. Strong neutrino heating is generally concentrated within an inner zone where the material is gravitationally bound (b < 0). The degree of symmetry of this bound material depends on the α-particle abundance. In the r_s0 = 125 km run, it is spherically stratified and the material with b>0 is generally excluded from it. Strong recombination heating is present both below and above the surface where b ≃ 0, indicating that it is mainly responsible for imparting positive energy to the shocked material. The mean shock radius expands by a factor ∼2.5 between the upper two frames in Figure 7, but the growth in the volume of positive-energy material is not accompanied by a significant expansion of the inner bound region, whose outer radius remains fixed at r ≃ r_α.

This segregation of bound from unbound material is broken when the shock is more compact initially, as is seen in the lower two panels of Figure 7. A single dominant accretion plume is continuously present, which funnels cold and dense material into the zone of strong neutrino heating. Alpha-particles are present only well outside the boundary between b < 0 and b>0. The competition between α-particle and neutrino heating is discussed in more detail in Section 4.1, and the influence of α-particles on the threshold heating rate for an explosion is examined in Section 5.

The accumulation of a bubble of hot material right behind the shock is a consequence of the balance of the buoyancy force acting within the bubble, and the ram pressure of the pre-shock material. The ratio of force densities is (Thompson 2000)

$$\begin{equation} \frac{F_{\rm buoy}}{F_{\rm ram}} \simeq \left(\frac{\rho - \rho _{\rm bubble}}{\rho }\right)\left(\frac{2GM/r_s}{v_r^2}\right) \Delta \Omega _{\rm bubble}, \end{equation} \tag{ 15 }$$

where r_s is the shock radius, v_r is the ambient radial flow speed, ρ is the ambient density, ρ_bubble the density of the bubble, and ΔΩ_bubble is its angular size. A low-density bubble (ρ − ρ_bubble ∼ ρ) can resist being entrained by the convective flow once it grows to a size $\Delta \Omega _{\rm bubble}\sim {\cal M}_{\rm con}^2$ Sr, where ${\cal M}_{\rm con}$ is the convective Mach number. On the other hand, the bubble must attain a much larger angular size ΔΩ_bubble ∼ 1 Sr if the buoyancy force is to overcome the upstream ram (|v_r| ∼ v_ff) and force a significant expansion of the shock surface. Figure 7 shows that the extent of the shock expansion is indeed correlated with the angular width of the region where hot material accumulates.

Another interesting feature of Figure 7 is the presence of secondary shocks, which are triggered once the outer shock becomes significantly nonspherical. Their locations are marked by discrete jumps in the rate of recombination heating. Secondary shocks are also prevalent throughout the nonlinear phase in the constant-ε models. Figure 8 shows the normalized pressure gradient (r/p)|∇p| for collapse models of both types, when H is just above threshold for an explosion (right before the shock hits the outer boundary of the simulation volume). The online version of this paper contains an animated version of Figure 8 showing the complete evolution. In both cases, secondary shocks extend over the whole post-shock domain, signaling the dissipation of supersonic turbulence which is stirred by accretion plumes that penetrate into the gain region.

Zoom In Zoom Out Reset image size

4.1. Distribution of Alpha-particle Recombination Heating

Heat input by neutrino absorption and by α-particle recombination have very different distributions within the shocked fluid: strong neutrino heating is concentrated inside r_s0, whereas recombination heating of a comparable amplitude is distributed throughout the settling flow. Strong recombination heating quite naturally extends below the zone where α-particles are present in significant numbers, as is seen in Figure 9. The first and third panels of this figure depict the pre-explosion steady state of the r_s0 = 75 km model with H = 1.02H_cr, while the second and fourth panels show the last time before the shock hits the outer boundary. At the latter time, one sees that the strongest recombination heating is concentrated in a layer where X_α ≲ 0.5, at the base of the extended α-rich plumes. Just as in the one-dimensional simulations (e.g., Figure 5), X_α approaches unity during shock breakout.

Zoom In Zoom Out Reset image size

The relative strength of neutrino heating and recombination heating depends on the initial radius of the shock, and on the Bernoulli parameter of the post-shock material. Figure 10 separates out cooling by α-particle dissociation from heating by recombination and neutrino radiation during the pre-explosion quasi-steady state (leftmost panels), at the onset of explosion (second panel left to right), and during breakout (two rightmost panels). See Figure 6 for comparison. The colored curves show the positive, negative, and net contributions from nuclear energy generation. The sharp negative spike near b = 0 represents α-particle dissociation in fresh, cold downflows. The formation of material with b>0 is primarily due to α-particle recombination in the r_s0 = 125 km run. As the initial radius of the shock is reduced with respect to r_α, neutrino heating makes a proportionately larger contribution near breakout.

Zoom In Zoom Out Reset image size

The strength of the boost given to the shock by recombination heating can be gauged by comparing de_nuc/dt to the adiabatic rate of change w_ad of the enthalpy of the flow (Equation (14)). Figure 11 shows the result for all three NSE sequences with H just above H_cr. In all cases, de_nuc/dt ≃ w_ad in various parts of the shocked fluid once the shock extends beyond a radius ≃r_α. Most of the heat input by recombination is concentrated where X_α ∼ 0–0.5, just as in spherical symmetry.

Zoom In Zoom Out Reset image size

Histograms of de_nuc/dt versus b and X_α are shown in the right panel of Figure 12. The rapid dissociation of α-particles in fresh downflows is represented by the long tail toward large negative values of de_nuc/dt, showing that the overall contribution of nuclear energy generation is negative. The α-particle concentration is very stratified, with higher X_α occurring at larger radius. Most of the mass with positive Bernoulli parameter is located at large radii.

Zoom In Zoom Out Reset image size

4.2. Residency Time

A long residency time of material in the gain region is commonly viewed as a key ingredient in a successful neutrino-driven explosion. To calculate t_res, we use the method described in Section 2.1.3: we first assign a unique "fluid time" t_F (Equation (12)) to each infalling radial mass shell in the simulation, which is effectively the time at which it passes the shock. We then define⁶ the residency time of the fluid as

$$\begin{equation} t_{\rm res} = t - t_F, \end{equation} \tag{ 16 }$$

where t is the present time. A related method (tracer particles) is used by Murphy & Burrows (2008) to calculate the residency time in collapse simulations with a more realistic EOS.

As material with positive Bernoulli parameter accumulates below the shock, we indeed find that its t_res grows larger. The shock starts running outward if the energy of this unbound material grows on a timescale shorter than the convective time. The right panel of Figure 6 shows the characteristic expansion time of the shock t_exp ∼ r_s,max/|dr_s,max/dt| (as measured at its outermost radius), alongside 〈t_res〉_vol (as measured within the material comprising the upper part of the residency time distribution). The final breakout of the shock seen in the left panel of Figure 6 corresponds to the time when t_exp ∼ 〈t_res〉_vol. This lengthening of the mean residency time can largely be ascribed to the increased dynamical time of the expanding shock. What changes most dramatically during breakout is the ratio of the expansion time to the dynamical time.

The breakout is a bit more gradual in the r_s0 =  125 km model with heating just above threshold for an explosion (H = 1.04H_cr; see animated version of Figures 7(a) and 7(b) in the online material). In this case, the expansion time of the shock remains somewhat longer than the residency time of material below the shock, which implies that the breakout depends on the continuing release of nuclear binding energy.

Note that large changes in the distribution of t_res are concentrated in regions of positive b. In the left panel of Figure12, we plot the distribution of b and t_res in the ε = 0.15v²_ff 0 run that is just above the threshold for an explosion. Regions with small or negative residency time represent freshly injected fluid. The distribution is stratified in b around t_res ∼ 40t_ff0 (which is approximately an overturn period of a convective cell, see Section 4.3 and Figure 13). Material with more negative b resides on average at a smaller radius. Fluid with a longer residency time has mostly positive b, corresponding to material transported upwards by convective cells.

Zoom In Zoom Out Reset image size

It is also apparent from the lower panel of Figure 12 that material with a longer residency time tends to have lower X_α, as is expected because it also tends to have a higher temperature.

4.3. Heat Engine in a Two-dimensional Explosion

Fresh material that is accreted through an oblique shock has a relatively low entropy, but once it reaches the base of the gain region it is exposed to an intense flux of electron-type neutrinos. Some of this heated material rises buoyantly, and forces an overturn of the fluid below the shock. Material with a longer residency time may therefore undergo multiple episodes of heating. On this basis, Herant et al. (1992, 1994) suggested that a convective flow would mediate a heat engine below the shock that would drive a secular increase in the energy of the shocked fluid.

We now investigate whether a heat engine operates in our simulations, and how it depends on the heating parameter H. We focus on a model with a constant nuclear dissociation energy, ε = 0.15v²_ff 0. In this class of models, the infalling heavy nuclei are completely broken up below the shock, and no heating by the reassembly of α-particles is allowed. As a first step, we average the convective flow over windows of width 50 t_ff 0, which de-emphasizes short-term fluctuations in the averaged velocity field 〈v〉_t. Figure 13 shows 〈v〉_t and 〈b〉_t (Equation (2)) at four different times in the run with H just above the threshold for explosion (H = 1.02H_cr). The radius of maximum heating (r ≃ 0.66r_s0) coincides with the lower boundary of the convective cells, across which material flows horizontally. The overturn period in these large-scale cells is ∼40–50t_ff0, as found by integrating streamlines of the mean flow (an example is shown on the lower left panel of Figure 13). Heated fluid accumulates in the region in between the top of convective cells and the shock. A strong deformation of the shock allows a plume of fresh material to descend diagonally between the convective cells. The tilt of this cold downflow intermittently flips in sign, and the averaged circulation pattern typically has an "∞" shape. The heating of fluid parcels in the two hemispheres is also intermittent, and sometimes two circulation flows are established simultaneously, thereby causing a bipolar expansion of the shock.

We have identified a useful figure of merit which connects a secular increase in the shock radius to the strength of neutrino heating at the base of the gain region. Figure 14 shows the absolute value of the Bernoulli parameter b at the radius of maximum heating (r_H,max ≃ 0.66r_s0) as a function of polar angle θ. In the top panel, the four sets of thin solid lines correspond to the four snapshots of Figure 13, and the bottom panel shows the analogous results for a nonexploding run. Overplotted as thick solid lines is the quantity

$$\begin{equation} \Theta = \bigg \langle \frac{{\LM\!}_H - {\LM\!}_C}{\rho }\bigg \rangle _{t,\theta *,r=r_{\rm H,max}} \frac{r_{\rm H,max}}{|\langle v_\theta \rangle _{t,\theta *}|}. \end{equation} \tag{ 17 }$$

This measures the specific energy that is absorbed from neutrinos by the material that flows laterally along the lower boundary of the convective cells. In Equation (17), the angular average of the heating rate and meridional velocity is restricted to a single convective cell.⁷

Zoom In Zoom Out Reset image size

In nonexploding models, the circulation in the gain region settles to a quasi-steady state, with no net amplification of the mass in material with positive b. The heat absorbed at the base of the convective cells is of the same order of the Bernoulli parameter of the fluid, that is, Θ ≲ |b|. In an exploding model, Θ will often exceed |b| by a factor 2–3. As is shown in the upper panel of Figure 14, Θ grows with time as the system approaches the explosion.

The stability of the averaged flow pattern appears to be, in part, an artifact of the axisymmetry of the flow. This imposes strong restrictions on the motion of convective cells, causing vorticity to accumulate on the largest spatial scales. Our observation that the bulk of the neutrino heating takes place within horizontal flows suggests that the ratio of heating timescale to radial advection time in the gain layer may be a less precise diagnostic of the conditions for explosion in axisymmetry: the horizontal convective velocity is typically low compared with the downward velocity of the main accretion plume. We do observe that the main accretion plume becomes strongly distorted near the threshold for an explosion, so that a significant fraction of the plume material enters one of the convection cells. This effectively decreases the amount of material that accretes to the protoneutron star and thus increases the overall advection timescale across the gain region.

An explosion occurs when the heating parameter H is raised above a critical value⁸ H_cr. We now explore how H_cr depends on the details of the EOS and the initial radius of the shock. One can express H simply in terms of the ratio of the heating rate ($4\pi r^3 {{\LM\!}}_H$) to the accretion luminosity ($GM\dot{M}/r$), in the idealized (but unrealistic) case where the flow is composed only of free nucleons and moves hypersonically. Then this ratio depends on H but not on the accretion rate $\dot{M} = 4\pi r^2 \rho (r) |v_{\rm ff}(r)|$. The precise value of the reference radius is unimportant; we choose r_s0, the shock radius in the time-independent, spherical flow solution at H = 0. Then

$$\begin{equation} \frac{4\pi {r_\mathrm{s0}}^3{\LM\!}_H[\rho _1(r_{s0})]}{GM\dot{M}/r_{s0} }{\Bigg |}_{X_\alpha = 0;\;{\cal M}=\infty } \;= \; {2H\over {r_\mathrm{s0}} v_{\rm ff\,0}^3}. \end{equation} \tag{ 18 }$$

This quantity is ∼10⁻³ to 10⁻² in the models we examine, which are below or near the threshold for explosion.

Note that the cooling in our model is concentrated at the base of the settling flow. As a result, the width of the gain region (relative to the shock radius) does not change significantly between different models. The critical heating parameter is therefore only indirectly related to the amplitude of the cooling function through the structure of the settling flow below the shock. Our purpose here is to explore how the critical heating rate depends on the strength of the gravitational binding of the shocked fluid to the collapsed core, and on the abundance of alpha particles.

Figure 15 displays H_cr for all of our model sequences. The abscissa is ε/v²₁, where ε is the nuclear dissociation energy and v₁ is the flow speed upstream of the shock in the initial configuration (that is, in the time-independent, spherical flow solution). In the case of the NSE equation of state, this quantity can be translated into an initial value of the shock radius using Figure 1. (Note that ε/v²₁ has a weak dependence on r_s in the NSE sequence.)

Zoom In Zoom Out Reset image size

A few interesting features of Figure 15 deserve comment. First, a comparison with Table 1 shows that the critical heating rate for explosion is ∼50%–70% of the maximum heating rate for which a steady-state flow solution can be found. The maximal heating parameter H_steady for a steady flow corresponds directly to the one first determined by Burrows & Goshy (1993) using a more realistic EOS. Note also that the values of H_cr in the one-dimensional and two-dimensional models are much closer to each other than they are to H_steady. This result is perhaps not surprising, given that the explosion is not immediate, but is approached through a series of transient fluid motions.

Second, H_cr is lower when NSE between n, p, and α is maintained below the shock. In this case, the dissociation energy at the shock is not fixed, but is (roughly) inversely proportional to radius. However, the difference between the NSE models and the constant-ε models is only ∼10% in H_cr when the shock starts out well below r_α (Equation (1)). An explosion is significantly easier when the fluid below the shock starts out with a significant population of α particles, as in the models with r_s0 = 125 km.

Third, H_cr tends to decrease with increasing ε/v²₁: a slightly lower heating rate per unit mass is required to explode a flow with a larger density contrast κ across the shock. Because almost all the gravitating mass is in the collapsed core, the gravitational binding energy of the gain region is approximately proportional to κ, whereas the net heat absorbed over the advection time is a stronger function of density, $t_{\rm adv} \int ({\LM\!}_H-{\LM\!}_C) {{d}}^3 r \propto \kappa ^2$. (One factor of κ comes from the advection time t_adv as given by Equation (13), and the other from the density dependence of ${{\LM\!}}_H$.) For example, Table 1 shows that κ is ∼1.6 times larger for ε/v²_ff 0 = 0.2 than for ε/v²_ff 0 = 0.1, and that H_cr is smaller by the inverse of the same factor.

Fourth, the two-dimensional runs all require less heating than their one-dimensional counterparts to explode. A major reason for this is that all two-dimensional configurations explode along one or both poles (see Figures 11 and 13), so that less material must be lifted through the gravitational field than in a fully spherical explosion. We have found that the precise value of the difference between the critical heating rate in the one-dimensional and two-dimensional explosions depends on the choice of r_*/r_s0, and therefore on the normalization C of the cooling function. The fact that we find a smaller difference than Murphy & Burrows (2008) may be a consequence of our simpler cooling function and equation of state.

The critical heating rate depends in an interesting way on the starting radius of the shock, in a way that points to the recombination of α-particles as an important last step in the transition to an explosion. Figure 16 shows that H_cr in the NSE models grows rapidly as the initial shock radius⁹r_s is pushed inside r_α. Here we normalize the heating parameter at a fixed physical radius, namely r_α. Translated into the context of a realistic core collapse, this means that the critical neutrino luminosity for an explosion decreases with increasing shock radius. The radius of the stalled shock depends, in turn, on the EOS above nuclear matter density: Marek & Janka (2009) find that a softer EOS corresponds to a larger shock radius, mainly due to the higher accretion luminosity onto the neutronized core. Here we have subsumed this uncertainty in the high-density EOS into a single free parameter, the ratio r_s0/r_α. Hydrodynamic instabilities are effective at driving an explosion to the extent that they push the shock radius close to r_α; beyond this point, the remainder of the work on the flow is done largely by α-particle recombination.

Zoom In Zoom Out Reset image size

One also notices from Figure 16 that the difference between H_cr in spherical and axial symmetry depends on the starting radius of the shock. The closer r_s0 is to r_α, the weaker the dependence of the critical heating rate on the dimensionality of the flow.

Overturns of the fluid below the shock can be triggered in two distinct ways: through the development of Ledoux convection in the presence of a strong negative entropy gradient, or via the nonlinear development of the SASI (a linear feedback between ingoing entropy and vortex waves, and an outgoing sound wave). We now show that the amplitude of the dipolar mode that is excited in the shock is strongly tied to the level of neutrino heating, and so thermal forcing plays a crucial role in maintaining the oscillation. To a certain extent, this distinction is of secondary importance, in the sense that memory of the linear phase of the instability is lost once the inflow of fresh material below the shock bifurcates from older shocked fluid. Nonetheless, the origin of the convective motions does have implications for the stability of ℓ = 1 and 2 modes in three-dimensional simulations: one expects that large-scale oscillations will change shape and direction more rapidly if they are triggered primarily by neutrino heating.

We can ask whether a heating parameter H that yields an explosion will also form an unstable entropy gradient below the shock. Convection develops through a competition between inward advection and neutrino heating. A detailed analysis by Foglizzo et al. (2006) shows that the parameter

$$\begin{equation} \chi \equiv \int \left|{\omega _{\rm BV}\over v_r}\right|\, dr \end{equation} \tag{ 19 }$$

must exceed a critical value ≃3 for unstable convective plumes to grow before being advected downward through the gain region below the shock. Here ω_BV is the Brunt–Väisälä frequency.

Using our initial flow models, we can translate H into a value for χ, and find (Table 1) that typically χ ∼ 5–20 at the threshold for a neutrino-driven explosion. The implication for convection below the shock is illustrated in Figure 17, which shows two snapshots for models with χ = 2.1 and 5.5, neither of which explodes. At the lower heating rate, the time required for convection to develop depends on the strength of the seed perturbation, whereas at the higher heating rate convective overturns develop rapidly within the layer of strong neutrino heating and spread throughout the post-shock region over a few dozen dynamical times. (The figure shows the result in the case where the seed perturbation is dominated by a small spherical startup error in the initial model.)

Zoom In Zoom Out Reset image size

We conclude that, near the threshold for a neutrino-driven explosion and for our given set of physical assumptions, convection is driven primarily through the development of a strong, negative entropy gradient within the gain region, rather than through the nonlinear development of SASI modes. The growth of the SASI requires at least a few oscillations, each with a period comparable to the advection time. The SASI is therefore subdominant when χ ≳ 3. It is worth comparing this with the results of Scheck et al. (2008), who forced the inner boundary of the simulation volume to move inward to model core contraction, thereby generating large advection velocities. The net result was that the flow barely reached χ ≃ 3 in exploding configurations. While this effect may be important for relatively prompt explosions, it should be kept in mind that the rate of contraction of the neutrinosphere has slowed subtantially a few hundred milliseconds after core bounce. Marek & Janka (2009) found evidence for the delayed explosion of a 15 M_☉ progenitor around 600 ms after core bounce, for which the effect of core contraction is not likely to dominate the dynamics.

In Paper I we considered the nonlinear, saturated state of the SASI in the absence of neutrino heating, and showed that the amplitude of the shock oscillations drops significantly as the dissociation energy ε is increased. We now explore how the rms amplitude of the shock oscillations correlates with the strength of heating. To eliminate the effect of secular shock motions around or above the threshold for explosion, we first calculate the running average 〈a_ℓ〉_50t of the shock Legendre coefficients a_ℓ over an interval 50t_ff0, and then calculate the rms of $\hat{a}_\ell \equiv a_\ell - \langle a_\ell \rangle _{50t}$ over the duration of each simulation. The result is plotted in Figure 18 as a function of H for two model sequences (ε = 0.15v²_ff 0 and NSE with r_s0 = 75 km).

Zoom In Zoom Out Reset image size

There is a clear trend of increasing ℓ = 1, 2 mode amplitude with increasing heating rate. This confirms our previous suggestion that large-amplitude shock oscillations require strong heating when the dissociation energy behind the shock exceeds ∼0.15v²_ff 0. The models with H = 0 reveal a slight exception to the overall trend: the rms amplitude of the shock oscillations appears larger than it does in models with small but finite heating rate, because the oscillations are strongly intermittent at H = 0 (see Paper I). The shock oscillations grow much stronger just below the threshold for explosion (exploding runs are marked by stars), above which they seem to saturate. Note also that their amplitudes do not vary much with the choice of the dissociation model. The rms amplitudes relative to the running average of a₀ at the threshold for explosion are {5%, 12%, 8%} for the ℓ = 1, 2, 3 modes in the ε/v²_ff 0 = 0.15 sequence, and {6%, 12%, 7%} in the NSE r_s0 = 75 km sequence.

We have performed an additional sequence of runs in which we drop an overdense shell with a given Legendre index ℓ through the shock (see also Paper I). This has the effect of selectively triggering individual SASI modes. Figures 19 and 20 display the Legendre coefficients of the shock alongside the angle-averaged entropy gradient. We find that the amplitude of the ℓ = 1 and 2 modes is strongly tied to the strength of convection. For both types of dissociation models, convection is quenched by the accretion flow when χ < 3: it grows intermittently in strength, but never reaches large enough amplitudes to significantly distort the shock surface. As a consequence, the entropy gradient remains shallow and negative most of the time. Coherent shock oscillations are also seen, but they have a low amplitude due to the large dissociation energy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Convection grows much more rapidly when χ>3, distorting the shock surface before the SASI has the chance to execute a few oscillations. In this case, the entropy gradient is initially more negative, but quickly flattens. Indeed, the ℓ = 1, 2 amplitudes only become large in models where neutrino-driven convection is strong enough to flatten the entropy gradient.

Another way of seeing that convection is the forcing agent behind shock oscillations when χ>3 is to analyze the distribution of vorticity in the gain region. Figure 21 shows a histogram of vorticity versus Mach number at three different instants in the evolution of the run with ε = 0.15v²_ff 0, ℓ = 1 seed perturbation, and χ = 3.9 (upper two panels in Figure 19). At t = 20t_ff0, convection is just getting started and the vortical motions are restricted to Mach numbers ≲ 0.3. However, by t = 35t_ff0 the Mach number distribution extends up to ${\cal M} \gtrsim 0.5$ and has almost reached its asymptotic form (t = 60t_ff0), at the same time that convection has filled the region below the shock. The dipolar mode of the shock develops a large amplitude only after this fully developed convective state has been reached. The convective rolls are a source of acoustic radiation (e.g., Goldreich & Kumar 1988), which will drive a dipolar oscillation of the shock if the overturn frequency is comparable to the frequency of the ℓ = 1 mode, |∇ × v| ∼ 2 × 2π/t_adv. This zone is marked by the vertical solid lines in Figure 21, and indeed encompasses most of the mass.

Zoom In Zoom Out Reset image size

We have investigated the effects of α-particle recombination and neutrino heating on the hydrodynamics of core-collapse supernovae, when the heating rate is pushed high enough to reach the threshold for an explosion. The effect of dimensionality has been probed by comparing one-dimensional (spherical) and two-dimensional (axisymmetric) time-dependent hydrodynamic calculations. Our main results can be summarized as follows:

• 1.  
  The critical heating parameter that yields an explosion depends sensitively on the starting position of the shock relative to r_α. This means that the critical neutrino luminosity depends sensitively on the stall radius of the shock and, in turn, on the core structure of the progenitor star and the density profile in the forming neutron star. Within the framework explored in this paper, we find two extreme types of explosion. In the first, neutrino heating does most of the work, with a significant final boost from α-particle recombination. In the second, neutrino heating is generally less important at promoting material below the shock to positive energies.
• 2.  
  During the final stages of an explosion, the heat released by α-particle recombination is comparable to the work done by adiabatic expansion. This heat is concentrated in material that has previously been heated by neutrinos. Significantly more energy is lost through α-particle dissociation in fresh downflows, so that nuclear dissociation remains on balance an energy sink within the accretion flow.
• 3.  
  The large-amplitude oscillations that are seen in one-dimensional runs near an explosion are the consequence of the ℓ = 0 SASI as modified by heating. In contrast with the ℓ = 1, 2 modes of a laminar accretion flow, the period of these oscillations is close to twice the post-shock advection time. The critical heating rate for an explosion (assuming constant mass accretion rate, neutrino luminosity, and inner boundary) corresponds to neutral stability for the ℓ = 0 mode.
• 4.  
  The critical heating parameter H_cr for an explosion is generally lower in axial than in spherical symmetry, but the difference becomes smaller as the starting radius of the shock approaches r_α. The difference depends somewhat on the ratio r_*/r_s0 and thus on the cooling efficiency and equation of state.
• 5.  
  Nonspherical deformations of the shock are tied to the formation of large-scale plumes of material with positive energy. Our two-dimensional explosions with a super-critical heating rate involve a large-scale convective instability that relies on the accumulation of vorticity on the largest spatial scales. Volume-filling convective cells are apparent in a time-averaged sense. Transient heating events create positive-energy material that accumulates in between the convective cells and the shock. A significant fraction of the heating occurs in horizontal flows at the base of the convective cells, which are fed by a dominant equatorial accretion plume. If the heating parameter is large enough, this results in an amplifying cycle and explosion.
• 6.  
  The amplitude of the ℓ = 1 and 2 modes correlates strongly with the value of the heating parameter, and is coupled to the appearance of vigorous neutrino-driven convection below the shock. In agreement with the work of Foglizzo et al. (2006) and Scheck et al. (2008), we find that χ ≈ 3 marks the transition from a strong linear instability in a nearly laminar flow below the shock, to a volume-filling convective instability. In all of our simulations, the threshold for explosion lies well within the latter regime. This highlights a basic difference between one-dimensional and two-dimensional explosions: the mechanism is fundamentally nonlinear in two dimensions.
• 7.  
  We have explored essentially one ratio of cooling radius to shock radius, namely r_*/r_s0 = 0.4 at zero heating (corresponding to r_*/r_s ∼ 0.2 near the threshold for an explosion). The growth of the ℓ = 1 SASI mode is strongest for this particular aspect ratio when ε = 0 (see Figure 12 of Paper I). As dissociation is introduced into the flow, we found that the peak growth rate moves to larger values of r_*/r_s0. On the other hand, detailed collapse calculations indicate ratios of neutrinosphere radius to shock radius that are even smaller than ∼0.2 following ∼100 ms after collapse (e.g., Marek & Janka 2009). We conclude that the l = 1 SASI mode is not being artificially suppressed by our choice of initial shock size.
• 8.  
  Vortical motions with a Mach number ∼0.3–0.5 first appear at the onset of convective instability around the radius of maximal neutrino heating, but before the dipolar mode of the shock reaches its limiting amplitude. These vortices are a source of acoustic waves, which have a similar period to the large-scale oscillation of the shock. Near the threshold for an explosion, the turbulence in the gain region becomes supersonic, as the existence of widespread secondary shocks attests. These shocks convert turbulent kinetic energy to internal energy, increasing the effective heating rate.

There are at least two reasons why explosions by the mechanism investigated here may be more difficult in fully three-dimensional simulations. First, the existence of more degrees of freedom for the low-order modes of the shock in three dimensions implies that the amplitude of individual shock oscillations is lower. As a result, it is more difficult for the shock to extend out to the radius where α-particle recombination gives it the final push. Second, an axisymmetric explosion that is driven by neutrino heating involves the accumulation of vorticity on the largest spatial scales, an effect that is special to two dimensions. A full resolution of these issues is possible only with high-resolution three-dimensional simulations.

We are grateful to Jonathan Dursi for scientific discussions and help with FLASH. We also thank Adam Burrows, Christian Ott, and Jeremiah Murphy for stimulating discussions. Careful and constructive comments by an anonymous referee helped to improve the presentation of this paper. The software used in this work was in part developed by the DOE-supported ASC/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago. Computations were performed at the CITA Sunnyvale cluster, which was funded by the Canada Foundation for Innovation. This research was supported by NSERC of Canada. R.F. is supported in part by the Ontario Graduate Scholarship Program.

We calculate the α-particle mass fraction X_α in nuclear statistical equilibrium by limiting the nuclear species to α-particles and free nucleons, and fixing the electron fraction Y_e = 0.5. We tabulate X_α and temperature T as a function of pressure p and density ρ, and then use these tables to calculate the rate of release of nuclear binding energy by the method described in Section 2.1.1. The temperature does not appear explicitly in the FLASH hydrodynamic solver, and only enters the flow equations indirectly through X_α.

We include the contributions to p from radiation, relativistic and partially degenerate electron–positron pairs, and nonrelativistic α-particles and nucleons. When k_BT>m_ec²/2, it can be written as (Bethe et al. 1980)

$$\begin{equation} p = \frac{1}{12}\frac{(k_B T)^4}{(\hbar c)^3}\left[ \frac{11\pi ^2}{15} + 2\eta ^2 + \frac{1}{\pi ^2}\eta ^4\right] + \left(1 - \frac{3}{4}X_\alpha \right) \frac{\rho }{m_u} k_B T, \end{equation} \tag{ A1 }$$

where η = μ_e/(k_BT) the normalized electron chemical potential, also known as degeneracy parameter, and ℏ, c, and m_u are Planck's constant, the speed of light, and the atomic mass unit, respectively. The density and degeneracy parameter are further related by

$$\begin{equation} \rho = \frac{m_u}{3\pi ^2 Y_e}\left(\frac{k_B T}{\hbar c}\right)^3\eta (\pi ^2+\eta ^2), \end{equation} \tag{ A2 }$$

where Y_e is the electron fraction. The equilibrium fraction of α-particles is given by the nuclear Saha equation,

$$\begin{eqnarray} X_n^2 X_p^2 &=& \frac{1}{2} X_\alpha \left[\frac{m_u n_Q(T)}{\rho }\right]^3 \exp {\left(-\frac{Q_\alpha }{k_B T} \right)}; \nonumber\\ && n_Q(T) = \left(\frac{m_u k_B T}{2\pi \hbar ^2} \right)^{3/2} \end{eqnarray} \tag{ A3 }$$

as supplemented by the conditions of mass and charge conservation,

$$\begin{equation} X_n + X_p + X_\alpha = 1 \end{equation} \tag{ A4 }$$

$$\begin{equation} X_p + \case{1}{2}X_\alpha = Y_e. \end{equation} \tag{ A5 }$$

In Equations (A3)–(A5), X_n and X_p are the mass fractions of free neutrons and protons, respectively, and Q_α = 28.3 MeV is the binding energy of an α-particle. Combining Equations (A1) and (A2) gives η and T in terms of p and ρ. The equilibrium mass fraction X^eq_α is calculated from ρ and T. For numerical calculations, we tabulate X^eq_α, ∂X^eq_α/∂ln ρ, and ∂X^eq_α/∂ln p for a grid of density and pressure. In addition, we tabulate partial derivatives of T to substitute into Equations (B10) and (B11).

Figure 22 shows contours of constant X^eq_α and constant entropy for different variables as a function of density. The entropy per nucleon is obtained by adding the contributions from the different components (e.g., Bethe et al. 1980),

$$\begin{eqnarray} S &=& \pi ^2 Y_e \frac{(11\pi ^2/15 + \eta ^2)}{\eta (\pi ^2+\eta ^2)} + \left(1-\frac{3}{4}X_\alpha \right)\nonumber\\ &&\times\,\left[\frac{5}{2} + \ln {\left\lbrace \frac{m_u n_Q(T)}{\rho } \right\rbrace }\right] -\,X_p\ln {X_p} - X_n\ln {X_n}\nonumber\\ && -\frac{1}{4}X_\alpha \ln {\left(X_\alpha /32\right)}. \end{eqnarray} \tag{ A6 }$$

The post-shock density in the initial configuration is typically ρ₂ ∼ 10⁹ g cm⁻³ with an entropy ∼10–15k_B/nucleon. The formation of α-particles that is seen in Figure 1 results from an expansion of the shock into the part of the thermodynamic plane in Figures 22(a) and 22(b) where ρ₂ < 10⁹ g cm⁻³ and T ≲ 1 MeV. In this regime, the electrons are nondegenerate and the pressure in photons and pairs begins to exceed the nucleon pressure. The dip in the adiabatic index seen in Figure 22(f) results from α-particle dissociation/recombination, which partially compensates the change in internal energy due to compression/expansion.

Zoom In Zoom Out Reset image size

We write down the ordinary differential equations that are used to compute the initial flow, and the density profiles in Figures 2(a)–2(c). The steady-state Euler equations in spherical symmetry are

$$\begin{equation} \frac{1}{v_r}\frac{{{d}}v_r}{{{d}}r} + \frac{1}{\rho }\frac{{{d}}\rho }{{{d}}r} + \frac{2}{r} = 0, \end{equation} \tag{ B1 }$$

$$\begin{equation} v_r\frac{{{d}}v_r}{{{d}}r} + \frac{1}{\rho }\frac{{{d}}p}{{{d}}r} + g = 0, \end{equation} \tag{ B2 }$$

$$\begin{equation} \rho v_r \frac{{{d}}e_\mathrm{int}}{{{d}}r} - \frac{p v_r}{\rho } \frac{{{d}}\rho }{{{d}}r} = {\LM\!}_H - {\LM\!}_C + {\LM\!}_\alpha, \end{equation} \tag{ B3 }$$

where e_int is the internal energy per unit mass, ${\LM\!}_H$, ${\LM\!}_C$, and ${\LM\!}_\alpha$ the source terms described in Equations (9)–(11), and g = GM/r². Since two variables suffice to describe the thermodynamic state of a system, we write

$$\begin{equation} \frac{{{d}}e_\mathrm{int}}{{{d}}r} \equiv E_p \frac{{{d}}p}{{{d}}r} + E_\rho \frac{{{d}}\rho }{{{d}}r} \end{equation} \tag{ B4 }$$

and

$$\begin{eqnarray} &&{\LM\!}_\alpha \equiv A_p \frac{{{d}}p}{{{d}}r} + A_\rho \frac{{{d}}\rho }{{{d}}r}. \end{eqnarray} \tag{ B5 }$$

The coefficients E_i and A_i encode the dependence on the equation of state. Replacing Equations (B4) and (B5) in (B3), and using Equations (B1) and (B2) to eliminate the pressure derivative, we obtain

$$\begin{equation} \frac{{{d}}\rho }{{{d}}r} = \frac{(\rho v_r E_p - A_p)\left(\rho g - 2\rho v_r^2/r\right) + \left({\LM\!}_H -{\LM\!}_C\right)}{(\rho v_r E_\rho - p v_r/\rho - A_\rho) + v_r^2(\rho v_r E_p - A_p)}. \end{equation} \tag{ B6 }$$

The coefficients in Equations (B4) and (B5) work out to

$$\begin{eqnarray} E_p & = & \frac{1}{(\gamma -1)\rho }\qquad \qquad \quad ({\rm constant }\; \gamma) \end{eqnarray} \tag{ B7 }$$

$$\begin{eqnarray} E_\rho & = & -\frac{p}{(\gamma -1)\rho ^2}\qquad \qquad ({\rm constant }\; \gamma). \end{eqnarray} \tag{ B8 }$$

For a constant-γ equation of state, e_int = p/[(γ − 1)ρ]. The pressure (Equation (A1)) in the NSE model described in Appendix A can be decomposed into contributions from relativistic particles and from nucleons, p = p_rel + p_mat, and the specific internal energy is

$$\begin{equation} e_\mathrm{int} = \frac{1}{\rho }\left(3p_\mathrm{rel} + \frac{3}{2}p_\mathrm{mat} \right) = 3\frac{p}{\rho } - \frac{3}{2}\left(1 - \frac{3}{4}X_\alpha \right)\frac{k_B T}{m_u}. \end{equation} \tag{ B9 }$$

One therefore finds

$$\begin{equation} E_p = \frac{3}{\rho } + \frac{9}{8}\frac{k_B T}{m_u}\frac{\partial X_\alpha }{\partial p} - \frac{3}{2}\left(1 - \frac{3}{4}X_\alpha \right)\frac{1}{m_u}\frac{\partial (k_B T)}{\partial p}, \quad(\hbox{NSE}) \end{equation} \tag{ B10 }$$

$$\begin{eqnarray} E_\rho &=& -\frac{3p}{\rho ^2} + \frac{9}{8}\frac{k_B T}{m_u}\frac{\partial X_\alpha }{\partial \rho } - \frac{3}{2}\left(1-\frac{3}{4}X_\alpha \right)\frac{1}{m_u}\frac{\partial (k_B T)}{\partial \rho }.\nonumber\\ &&\qquad\qquad\qquad\qquad\qquad\qquad (\hbox{NSE})\qquad \end{eqnarray} \tag{ B11 }$$

The initial post-shock solution is obtained by integrating the above equations from r_s to an inner radius r_* at which the flow stagnates. We iterate the normalization of the cooling function in Equation (9) so that r_* = 0.4r_s0 in the absence of heating. When adding heating, the cooling normalization and r_* are kept fixed, which results in an expansion of the shock from its initial position to r_s>r_s0 (Figures 2(a) and 2(b)). The initial Mach number at the inner boundary is chosen so as to satisfy

$$\begin{eqnarray} &&\hspace{-48pt} \bigg | \sum _i \left({\LM\!}_{H,\,i} - {\LM\!}_{C,\,i} + {\LM\!}_{\alpha,\,i}\right) V_i \bigg |\nonumber\\ && \simeq 0.995\left[ \frac{GM}{r_*} - \varepsilon (t=0) \right] \big | \dot{M} \big |, \end{eqnarray} \tag{ B12 }$$

where the sum is taken over the computational cells below the shock at our fixed resolution (see Section 2.1.3), V_i is the volume of each computational cell, and the source terms are evaluated at the inner radial cell face. The numerical coefficient on the right-hand side depends on the radial resolution, and is chosen empirically to prevent runaway cooling due to discreteness effects in time-dependent calculations. The resulting inner Mach number is typically 10⁻³ to 10⁻².

When including α-particles in the EOS, one needs to calculate self-consistently the value of X^eq_α below the shock, the corresponding dissociation energy ε(t = 0) (Equation (8)), and compression factor κ (Equation (3)). The density upstream of the shock is obtained from

$$\begin{equation} \rho _1(r_s) = \frac{\dot{M}}{4\pi r_s^2 |v_1(r_s)|}, \end{equation} \tag{ B13 }$$

where

$$\begin{equation} v_1(r_s) = -\frac{v_\mathrm{ff}(r_s)}{\sqrt{1 + 2\mathcal {M}_1^{-2}(r_s)/(\gamma -1)}} \end{equation} \tag{ B14 }$$

is the upstream velocity at r = r_s, while the upstream pressure satisfies

$$\begin{equation} p_1(r_s) = \frac{\rho _1(r_s) [v_1(r_s)]^2}{\gamma \mathcal {M}_1^2(r_s)}. \end{equation} \tag{ B15 }$$

Equations (B13) and (B15) are transformed to physical units for input to the NSE model by adopting ${\cal M}_1({r_\mathrm{s0}}) = 5$, $\dot{M} = 0.3\,M_\odot$ s⁻¹, M = 1.3 M_☉, and a particular value for the shock radius r_s0 in the absence of heating.

Here we give some further details of the time evolution of our initial models using FLASH2.5. Heating and cooling are applied in an operator split way in between hydrodynamic sweeps. Nuclear dissociation in the constant-ε model is implemented through the fuel+ash module in FLASH, with the modifications described in Paper I. The rate of change of specific internal energy is computed using the current hydrodynamic variables and timestep, after which the EOS subroutine is called to ensure that the variables are thermodynamically consistent.

In the NSE nuclear dissociation module, numerical stability is maintained using an implicit update of the pressure in between hydro sweeps,

$$\begin{equation} p_\mathrm{new} = p_\mathrm{cur} + (\gamma -1)\rho e_\mathrm{nuc}(\rho,p_\mathrm{new}), \end{equation} \tag{ C1 }$$

where e_nuc is the energy generation per timestep in Equation (6), and the subscripts cur and new refer to the current and new value of the pressure, respectively. The density is kept constant across this step, so as to be consistent with the other source terms. Equation (C1) usually converges in 3–4 Newton iterations, adding a negligible overhead to our execution time. We restrict the timestep of the simulation so that, in addition to the standard Courant–Friedrichs–Levy condition, it enforces |e_nuc| < 0.8(p/ρ)/(γ − 1). To prevent α-particle recombination in the cooling layer (due to the decrease of internal energy), we adopt a cutoff in density, so that X_α = 0 if ρ>3 × 10¹⁰ g cm⁻³.

Both the constant-ε and the NSE dissociation modules require that the Mach number remain below a fixed value $\mathcal {M}_{\rm burn}$ for burning, which has the effect of preventing dissociation or recombination upstream of the shock (Paper I). This results in a small amount of incomplete burning in the presence of strong shock deformations, a phenomenon which is also encountered in the full collapse problem. We set the threshold to $\mathcal {M}_{\rm burn} = 2$ in most of our simulations. The critical heating parameter H_cr depends weakly on ${\cal M}_{\rm burn}$: changes in ${\cal M}_{\rm burn}$ cause small changes in the amount of unburnt material with zero Bernoulli parameter, and only slightly alters the net energy of the gain region. At the outer boundary of the simulation volume, $\mathcal {M}_{\rm burn}$ is just below the Mach number of the upstream flow. We have tried to expand the outer boundary to r = 9r_s0 (with a somewhat smaller $\mathcal {M}_{\rm burn}$) and found that runs that did hit r = 7r_s0 still hit the new outer boundary.

Aside from the inclusion of heating and NSE dissociation, the numerical setup for our runs is identical to that of Paper I. In that paper, the numerical output was verified by comparing the measured growth rates of linearly unstable modes of the shock with the solution to the eigenvalue problem. We have tested the implementation of our NSE dissociation model by verifying that, in the absence of initial perturbations, our steady-state initial conditions remain steady. Spherical transients present in the initial data die out in a few ℓ = 0 oscillation cycles, and are present even when nuclear burning is omitted. (See Paper I for a more extended description.)