Document Sample
ms Powered By Docstoc
  Preprint typeset using LTEX style emulateapj v. 6/22/04

                                       DYNAMICS OF A SPHERICAL ACCRETION SHOCK WITH
                                                  RODRIGO F ERNÁNDEZ 1 AND C HRISTOPHER T HOMPSON 2
                                                                    Submitted to ApJ

              We investigate the effects of neutrino heating and α-particle recombination on the hydrodynamics of core-
           collapse supernovae. Our focus is on the non-linear dynamics of the shock wave that forms in the collapse, and
           the assembly of positive energy material below it. To this end, we perform one- and two-dimensional, time-
           dependent hydrodynamic simulations with FLASH2.5. These generalize our previous calculations by allowing
           for bulk neutrino heating and for nuclear statistical equilibrium between n, p and α. The heating rate is freely
           tunable, as is the starting radius of the shock relative to the recombination radius of α-particles. An explosion
           in spherical symmetry involves the excitation of an overstable mode, which may be viewed as the = 0 version
           of the ‘Standing Accretion Shock Instability’. In 2D simulations, non-spherical deformations of the shock are
           driven by plumes of material with positive Bernoulli parameter, which are concentrated well outside the zone of
           strong neutrino heating. The non-spherical modes of the shock reach a large amplitude only when the heating
           rate is also high enough to excite convection below the shock. The critical heating rate that causes an explosion
           depends sensitively on the initial position of the shock relative to the recombination radius. Weaker heating
           is required to drive an explosion in 2D than in 1D, but the difference also depends on the size of the shock.
           Forcing the infalling heavy nuclei to break up into n and p below the shock only causes a slight increase in the
           critical heating rate, except when the shock starts out at a large radius. This shows that heating by neutrinos
           (or some other mechanism) must play a significant role in pushing the shock far enough out that recombination
           heating takes over.
           Subject headings: hydrodynamics — instabilities — nuclear reactions — shock waves — supernovae: general

                             1. INTRODUCTION                                   critical amplitude, to a bifurcation between freshly infalling
   Although tremendous progress has been made on the mech-                     material and material shocked at earlier times (Blondin et al.
anism of core-collapse supernovae in recent years, we still do                 2003). Such an oscillation is easily excited in a spherical flow
not have a clear picture of a robust path to an explosion in stars             composed of a zero-energy, polytropic fluid (Blondin & Mez-
that form iron cores. Heating by the absorption of electron-                   zacappa 2006) via a linear feedback between ingoing vor-
type neutrinos significantly modifies the settling flow below                     tex and entropy waves and outgoing sound waves (Foglizzo
the bounce shock but – in spite of early positive results (Bethe               et al. 2007). It can also be excited indirectly by neutrino heat-
& Wilson 1985) – explosions are obtained in spherical col-                     ing, which if sufficiently strong will drive large-scale buoyant
lapse calculations only if the progenitor star is lighter than                 motions below the shock (Herant et al. 1994; Foglizzo et al.
about 10-12 M (Kitaura et al. 2006). More massive stars fail                   2006).
to explode in spherical symmetry (Liebendörfer et al. 2001).                      A significant sink of thermal energy in the accretion flow
   Two-dimensional collapse calculations show strong defor-                    arises from the dissociation of heavy nuclei. The iron, silicon,
mations of the shock and convective motions below it (Bur-                     and oxygen that flows through the shock is broken up into
rows et al. 1995; Janka & Mueller 1996; Buras et al. 2006a,b;                  α-particles and nucleons when it is exposed to the high tem-
Burrows et al. 2006, 2007; Marek & Janka 2007). It has                         perature (> 1 MeV) of the postshock region. The Bernoulli
long been recognized that convection increases the residency                   parameter b of the shocked fluid then becomes substantially
time of settling material in the zone of strong neutrino heat-                 negative. A significant fraction of this dissociation energy can
ing (Herant et al. 1992). It is also becoming clear that mul-                  be recovered if nucleons recombine into α-particles (Bethe
tidimensional explosions require the assembly of a smaller                     1996). But for this to occur, a decrease in the temperature is
amount of material with positive energy, but the details of how                required and thus the shock must expand significantly beyond
this happens remain murky.                                                     the radius at which it typically stalls (∼ 100 − 150 km).
   An early treatment of shock breakout by Bethe (1997) fo-                       One of the primary goals here, and in a previous paper (Fer-
cused on the strong gradient in the ram pressure of the in-                    nández & Thompson 2008) [hereafter Paper I], is to gauge
falling material, but implicitly assumed that the shocked ma-                  the relative importance of these effects in setting the stage
terial had already gained positive energy. If large-scale den-                 for a successful explosion. The persistence and amplitude of
sity inhomogeneities are present below the shock, they will                    a dipolar oscillation can only be reliably measured in three-
trigger a finite-amplitude, dipolar instability, thereby allowing               dimensional (3D) simulations (there are preliminary indica-
accretion to continue simultaneously with the expansion of                     tions that it is less prominent in 3D than in 2D; Iwakami et al.
positive-energy fluid (Thompson 2000). The accretion shock                      2008). On the other hand, the interplay between α-particle
is also capable of a dipolar oscillation which leads, above a                  recombination and hydrodynamical instabilities has received
                                                                               little attention. Although recombination is certainly present in
  1 Department of Astronomy and Astrophysics, University of Toronto.           previous numerical studies which employ finite-temperature
Toronto, Ontario M5S 3H4, Canada.                                              equations of state (EOSs), it should be kept in mind that con-
  2 CITA, 60 St. George St., Toronto, Ontario M5S 3H8, Canada.
                                                                               siderable uncertainties in the EOS remain at supranuclear den-
2                                                   FERNÁNDEZ & THOMPSON

sities. A softening or hardening of the EOS feeds back on the       shock in 2D are driven by material with positive Bernoulli pa-
position of the shock for a given pre-collapse stellar model        rameter, which generally resides outside the radius rα where
(Marek & Janka 2007). The parametric study of the critical          the gravitational binding energy of an α-particle is equal to
neutrino luminosity by Murphy & Burrows (2008) is based             its nuclear binding energy. The recombination of α-particles
on a single EOS; they obtain explosions in which the shock          plays a major role in creating this positive-energy material,
seems to break out from nearly the same radial position at          but for this to happen the shock must be pushed beyond ∼ 200
∼ 250 − 300 km. Variations in the density profile of the pro-        km from the neutronized core.
genitor star will similarly modify the position of the shock, the      In this paper, we consider only neutrino heating as the im-
concentration of α-particles below it, and the critical neutrino    petus for the initial expansion of the shock, rather than more
luminosity for an explosion.                                        exotic effects such as rotation or magnetic fields. We find that
   In this paper, we study the interplay between non-spherical      the critical heating rate is a strong function of the initial po-
shock oscillations, neutrino heating, and α-particle recombi-       sition of the shock with respect to rα , which implies that a
nation, when the heating rate is pushed high enough for an          much higher neutrino luminosity is needed to revive a shock
explosion. Our focus is on the stalled shock phase, between         that stalls well inside rα . The difference in the critical heat-
∼ 100 ms and 1 s after bounce. In a one-dimensional calcu-          ing rate between 1D and 2D simulations also depends on the
lation, the accretion flow reaches a quasi-steady state during       size of the shock, and thence on the structure of the forming
this interval, and the shock gradually recedes (e.g. Liebendör-     neutron star.
fer et al. 2001; Buras et al. 2006a).                                  We find that the shock develops a dipolar oscillation with a
   Our approach is to introduce equations of state of increasing    large amplitude only when the heating rate is also high enough
complexity into one- and two-dimensional, time-dependent            to trigger a strong convective instability. We therefore sur-
hydrodynamic simulations. To this end, we use the code              mise that buoyant motions driven by neutrino heating play a
FLASH2.5 (Fryxell et al. 2000), which is well tested in prob-       major role in driving the dipolar oscillations that are seen in
lems involving nuclear energy release in compressible fluids         more complete simulations of core collapse. Some evidence
(Calder et al. 2002). We adopt a steady state model as our          is found that acoustic wave emission by the convective mo-
initial condition, and a constant mass accretion rate, neutrino     tions also plays a role. We investigate the possibility of a heat
luminosity, and fixed inner boundary. The steady-state ap-           engine within the gain layer (the region with a net excess of
proximation to the stalled shock phase was first introduced          neutrino heating over cooling). We find that most of the heat
by Burrows & Goshy (1993), and has recently been used by            deposition by neutrinos is concentrated in lateral flows at the
Ohnishi et al. (2006) to study the non-linear development of        base of the prominent convective cells. At the threshold for an
the shocked flow with a semi-realistic equation of state and         explosion, neutrino heating plays a key role in pushing mate-
neutrino heating. In Paper I we modeled the accretion flow as        rial to positive Bernoulli parameter, but only if the shock starts
a polytropic fluid, from which a fixed dissociation energy is         well inside rα .
removed immediately below the shock, and allowed for neu-              The plan of the paper is as follows. Section 2 presents our
trino cooling but not heating.                                      numerical setup, treatment of heating, cooling, and nuclear
   Here we generalize this model to allow both for heating,         dissociation, and outlines the sequences of models. Sections 3
and for nuclear statistical equilibrium (NSE) between neu-          and 4 show results from 1D and 2D simulations, respectively.
trons, protons, and α-particles. Heating is introduced in a         We focus on the relative effectiveness of neutrino heating and
simple, parametrized way, without any attempt at simulating         α-particle recombination in driving an explosion, and the re-
neutrino transport. The nuclear abundances are calculated as        lation between Bernoulli parameter and large-scale deforma-
a function of pressure p and density ρ, using a complete finite-     tions of the shock. The critical heating rate for an explosion
temperature, partially degenerate EOS. Our model for the            is analyzed in §5, and the competition between advective-
shocked material retains one significant simplification from          acoustic feedback and convective instability is discussed in
Paper I: we do not allow the electron fraction Ye to vary with      §6. We summarize our findings in §7. The appendices con-
position below the accretion shock. Here there are two com-         tain details about our EOS and the numerical setup.
peting effects: electron captures tend to reduce Ye , whereas
absorption of νe and νe tends to drive high-entropy material
                        ¯                                                                2. NUMERICAL MODEL
below the shock toward Ye 0.5. Since we are interested es-             As in Paper I, the initial configuration is a steady, spheri-
pecially in the dynamics of this high-entropy material, we set      cally symmetric flow onto a gravitating point mass M. The
Ye = 0.5 when evaluating the α-particle abundance. To obtain        flow contains a standing shock wave, and the settling flow be-
a realistic density profile, we continue to approximate the in-      low the shock cools radiatively in a narrow layer outside the
ternal energy of the fluid as that of a polytropic fluid with a       inner boundary of the simulation volume. The space of such
fixed index γ = 3 . In reality, the equation of state between the    models is labeled basically by three parameters: accretion rate
neutrinosphere and the shock depends in a complicated way            ˙
                                                                    M, luminosity Lν in electron neutrinos and anti-neutrinos, and
on the degeneracy of the electrons and the effects of electron      the radius r∗ of the base of the settling flow, which corre-
captures. The consequences of introducing these additional          sponds roughly to the neutrinosphere radius. The mass M
degrees of freedom will be examined in future work.                 of the collapsed material represents a fourth parameter, but
   In spite of these simplifications, our results already show       it covers a narrower range than the other three. The infalling
many similarities with more elaborate collapse calculations.        material is significantly de-leptonized before hitting the shock
One-dimensional explosions are due to a global instability re-      only in the first 50 ms or so of the collapse (e.g. Liebendörfer
sembling the one-dimensional Standing Accretion Shock In-           et al. 2001).
stability (SASI), but modified by heating. As in Paper I, we            We explore a two-dimensional surface through this three-
find that the period of the = 0 mode remains close to twice          dimensional parameter space by i) fixing the ratio of r∗ to the
the post-shock advection time. Strong deformations of the           initial shock radius rs0 in the absence of heating (r∗ /rs0 = 0.4);
                                                                    ii) allowing rs0 to vary with respect to an appropriately chosen
                                   ALPHA PARTICLE RECOMBINATION IN SUPERNOVAE                                                                    3

physical radius; and then iii) increasing the level of heating       modified so as to allow a decrement ε in b across the shock.
until an explosion is uncovered. In the full problem, the shock      The resulting compression factor is (Paper I)
radius at zero heating is a unique function of M and r∗ , with
a small additional dependence on M and the composition of                             ρ2
the flow outside the shock (Houck & Chevalier 1992). The                         κ≡       = (γ + 1)      γ + M−2 −
secular cooling of the collapsed core forces a gradual decrease
in r∗ , and M also varies with time and with progenitor model.                                                      2ε
   Given the important role that α-particle recombination                                        1 − M−2 + (γ 2 − 1) 2
                                                                                                      1                                 ,      (3)
plays in the final stages of an explosion, we implement ii) by
referencing rs0 to the radius where the gravitational binding
energy of an α-particle equals its nuclear binding energy,           which reduces to κ → (γ + 1)/(γ − 1) for M1 → ∞ and ε =
                                                                     0. Throughout this paper, the specific nuclear dissociation
                        GMmα                                         energy ε is defined to be positive. The subscripts 1 and 2
                  rα =            254M1.3 km.                  (1)
                          Qα                                         denote upstream and downstream variables, respectively.
Here Qα 28.30 MeV is the energy needed to break up an                   All flow variables are made dimensionless by scaling radii
α-particle into 2n and 2p, mα the mass of an α-particle, and         to rs0 , velocities to vff 0 = (2GM/rs0 )1/2 , timescales to tff 0 =
M1.3 = M/(1.3M ). Choice i) allows us to consider models             rs0 /vff 0 , and densities to the upstream density at r = rs0 , ρ1 (rs0 )
that have, implicitly, both a range of physical values of r∗ and     [eq. B13]. See Paper I for further details. Throughout the
a range of M. It is, of course, made partly for computational        paper we denote the average of a quantity F(X, ...) over some
simplicity (the limited size of the computational domain) and        variable X by F X .
also to facilitate a comparison between models that have dif-
ferent values of rs0 /rα . Nuclear dissociation is taken into ac-                         2.1.1. Nuclear Dissociation
count either by removing a fixed specific energy ε right below
                                                                        We model nuclear dissociation in two ways. First, we re-
the shock, or by enforcing NSE between n, p, and α through-
                                                                     move a fixed specific energy ε right below the shock, as done
out the settling flow. Once this choice is made, the normal-
                                                                     in Paper I. This represents the prompt and complete breakup
ization of the cooling function is adjusted to give r∗ /rs0 = 0.4.
                                                                     of whatever heavy nuclei are present in the upstream flow.
The heating rate remains freely adjustable thereafter.
                                                                     The main limitation of this approximation is that the dissoci-
   We adopt this simplification because we do not intend to
                                                                     ation energy does not change with the radius (or inclination)
find the precise value of the critical neutrino luminosity, but
                                                                     of the shock. The main advantage is simplicity: ε is indepen-
instead to probe the behavior of the system around this critical
                                                                     dent of any dimensional parameters and can be expressed as a
point, whatever its absolute value.
                                                                     fraction of v2 0 .
   We now describe the key components of this model in more
                                                                        We also use a more accurate dissociation model which
detail, and explain the setup of the hydrodynamic calcula-
                                                                     allows for NSE between α-particles and nucleons.3 Dur-
tions. As in Paper I, the time evolution is carried out using the
                                                                     ing the stalled shock phase of core-collapse supernovae, the
second-order, Godunov-type, adaptive-mesh-refinement code
                                                                     shock sits at r ∼ 100 − 200 km, with a postshock temperature
FLASH2.5 (Fryxell et al. 2000).
                                                                     T > 1 MeV and density ρ 109 g cm−3 . In these conditions,
                     2.1. Initial Conditions                         the heavy nuclei flowing through the shock are broken up into
                                                                     α, p, and n.
  The introduction of heating causes a change in the structure          A range of isotopes are present in the iron core of a mas-
of the initial flow configuration. The radius rs of the shock in       sive star as well as in nuclear burning shells (Woosley et al.
the time-independent solution to the flow equations increases         2002), but since the binding energy per nucleon varies only
with heating rate; that is, rs ≥ rs0 . The material above the        by ∼ 10% we simply assume a single type of nucleus in the
shock is weakly bound to the protoneutron star, and in practice      upstream flow. We focus here on the later stages of the stalled
can be taken to have a zero Bernoulli parameter                      shock phase, during which the oxygen shell is accreted. An
                       1        γ p GM                               energy QO 14.44 MeV must be injected to dissociate an
                   b = v2 +           −     .                (2)     16
                       2      γ −1 ρ    r                               O nucleus into 4 α-particles (Audi et al. 2003), which cor-
Here v is the total fluid velocity, which is radial in the initial    responds to the specific dissociation energy
condition, p is the pressure, ρ is the mass density, and G is                             QO                r
Newton’s constant. The flow upstream of the shock is adia-                          εO =      = 0.038M1.3        v2 (r).                        (4)
batic and has a Mach number M1 = 5 at a radius r = rs0 .                                  mO             150 km ff
  The composition of the fluid is very different upstream and         Here mO 16mu is the mass of an oxygen nucleus, with mu
downstream of the shock. Changes in internal energy due to           the atomic mass unit. The smallness of this number indicates
nuclear dissociation and recombination are taken into account        that little oxygen survives in the post-shock flow, and so we
using the two models described in §2.1.1. For the internal           set the equilibrium mass fraction of oxygen to zero below
energy density of the fluid e, we continue to use the polytropic                     eq
                                                                     the shock, XO = 0. The binding energy of an α-particle is
relation e = p/(γ − 1); hence the second term on the right-          of course much larger, giving
hand side of eq. (2). Because we are not explicitly including
changes in electron fraction due to weak processes, we keep                               Qα        −1     r
                                                                                   εα =      = 0.32M1.3        v2 (r).                         (5)
γ = 4/3 for both models of nuclear dissociation. This largely                             mα            150 km ff
determines the density profile inside a radius ∼ 2 rα , where            3 Although heavier nuclei can begin to recombine once the shock moves
the α-particle abundance is very low.                                significantly beyond rα , this generally occurs only after the threshold for an
  The upstream and downstream flow profiles are connected              explosion has been reached, and makes a modest additional contribution to
through the Rankine-Hugoniot jump conditions, which are              the recombination energy.
4                                                             FERNÁNDEZ & THOMPSON
   We find that α-particles appear in significant numbers only
at relatively large radii ( 0.5 rα ) and in material that has ei-                            .             -1                    eq
                                                                                                                           Xα [ 56Fe ]
                                                                                             M = 1.0 M . s
ther i) been significantly heated by electron neutrinos closer                                .             -1
                                                                                 0.8         M = 0.3 M . s
to the neutrinosphere; or ii) been freshly shocked outside rα .                              .             -1
                                                                                             M = 0.1 M . s
The electrons are only mildly degenerate in material that has
a high entropy and α-particle content, so that neutrino heat-                                                                                  eq
                                                                                                                                             Xα [ 16O ]
ing drives Ye close to ∼ 0.5 (or even slightly above: see, e.g.,
Buras et al. 2006b). We therefore set Ye = 0.5 in the Saha
equation that determines the equilibrium mass fractions Xn ,                     0.4
  eq        eq         eq     eq
Xp and Xα = 1 − Xn − Xp . These quantities are tabulated
as functions of p and ρ using an ideal, finite-temperature and                                                                            2
                                                                                                                          ε(t=0) / vff0 [ 56Fe ]
partially degenerate equation of state for electrons and nucle-                  0.2
ons; see Appendix A for details. Specific choices must then
be made for the parameters rs0 , M, and M; we generally take
                                                                                                                 ε(t=0) / vff0 [ 16O ]
M = 1.3 M and M    ˙ = 0.3 M s−1 , but allow rs0 to vary. An in-                   0
                                                                                    50            100           150               200                250    300
vestigation of how changes in Ye feed back onto the formation                                                         rs0 [km]
of α-particles is left for future work.
   A specific energy                                                                F IG . 1.— Equilibrium mass fraction of α-particles Xα , and ratio of initial
                                                                                 dissociation energy ε(t = 0) [eq. 8] to v2 behind a spherical shock positioned
            enuc = −XO (εO + εα ) − Xα − Xα [ρ, p] εα ,                   (6)    at radius rs0 . Curves of different shadings correspond to different mass ac-
                                                                                 cretion rates. The Rankine-Hugoniot shock jump conditions and dissociation
is either released or absorbed within a single time step (it can                 energy are calculated self-consistently, as described in Appendix B. Square
be of either sign). Here XO is non-vanishing only for fluid ele-                  brackets refer to the upstream composition of the accretion flow, which for
ments that have just passed across the shock, and we have set                    simplicity is taken to be pure 56 Fe or 16 O. The Mach number upstream of the
  eq                                                                             shock is M1 = 5, and the central mass is M = 1.3M .
XO = 0. The quantity (6) is introduced as an energy source
term in FLASH, and from it one readily obtains a rate of re-
lease of nuclear binding energy per unit mass,
                           denuc    enuc
                                 ≡       ,                     (7)               in a nearly adiabatic settling flow). In more realistic collapse
                             dt      ∆t                                          calculations, Ye grows with radius between the neutrinosphere
where ∆t is the simulation time step.                                            and the shock, but ρ tends to decrease more rapidly than ∼ r−3
   In the initial condition, the dissociation energy at the shock                (e.g. Buras et al. 2006a). Our choice of cooling function ap-
is obtained from equation (6) using XO = 1 and Xα = 0 up-                        pears to widen the gain region slightly compared with these
stream of the shock:                                                             calculations, and therefore to reduce the critical heating rate
                                     eq                                          for an explosion. The bulk of the cooling occurs in a narrow
                ε(t = 0) = εO + 1 − Xα [ρ2 , p2 ] εα .                    (8)
                                                                                 layer close to the accretor at r = r∗ , and the accreted material
Figure 1 shows how ε(t = 0) and Xα depend on the shock ra-                       accumulates in the first few computational cells adjacent to
dius rs0 , for upstream flows composed4 of pure 16 O and 56 Fe,                   the inner boundary without a major effect on the rest of the
and for different values of M. The dissociation energy is ap-                    flow.
proximately constant inside ∼ 75 km, where the downstream                           We model neutrino heating as a local energy generation rate
flow is composed of free nucleons, but decreases at greater                       per unit volume of the form
distances, remaining ∼ 40% of the gravitational binding en-                                                LH = H(1 − Xα )ρ/r2 .                           (10)
ergy at the shock. The mass fraction of α-particles reaches
50% at r = 150 − 175 km, with a weak dependence on M. ˙                          The normalization constant H measures the strength of the
                                                                                 heating. The factor (1 − Xα ) accounts for the fact that the cross
      2.1.2. Heating and Cooling in the Post-shock Flow                          section for neutrino absorption by α-particles is much smaller
                                                                                 than that for free nucleons (e.g., Bethe 1993). For simplicity,
  To allow direct comparison with our previous results, we
                                                                                 we do not include the flux factor due to the transition between
employ a cooling rate per unit volume of the form
                                                                                 diffusion and free-streaming. Our focus here is on the nature
                             LC = Cpa ρb−a ,                              (9)    of the instabilities occurring in the flow near the threshold for
                                                                                 an explosion, and we do not attempt a numerical evaluation of
with a = 1.5, b = 2.5, and C a normalization constant. As in Pa-                 the critical heating rate.
per I, we include a gaussian entropy cutoff to prevent runaway                     An additional energy source term arises from the change in
cooling. The exponents in eq. (9) represent cooling domi-                        the equilibrium fraction of α-particles as they are advected in
nated by the capture of relativistic, non-degenerate electrons                   the steady state initial solution. The instantaneous adjustment
and positrons on free nucleons (e.g., Bethe 1990). Inside the                    of Xα to its equilibrium value, combined with eq. (6), yields
radius where the electrons become strongly degenerate, and                       an energy generation rate per unit volume
α-particles are largely absent, one has LC ∝ pe n p ∝ (Ye ρ)3 .                                            eq          eq       eq
This gives essentially the same dependence of LC on r as eq.                                            dXα          ∂Xα dρ ∂Xα dp
                                                                                         Lα = ρvεα            = ρvεα        +        .                     (11)
(9) when Ye = constant and γ = 3 (corresponding to ρ ∝ r−3                                               dr           ∂ρ dr   ∂ p dr
   4 In the case where the upstream flow is pure 56 Fe, we replace ε in eq. (8)   This energy generation rate is negative, as the temperature
with εFe = QFe /mFe 0.093M1.3 (r/150 km) v2 (r), and set the electron frac-
                                −1                                               increases inwards and thus the α-particle fraction decreases
tion to Ye = 26/56 in the NSE calculation behind the shock.                      with decreasing radius (v is negative).
                                    ALPHA PARTICLE RECOMBINATION IN SUPERNOVAE                                                                    5

                    2.1.3. Numerical Setup                                       10000
                                                                                           (a)                                         2
                                                                                                                                ε / vff0 = 0.1
   In our time dependent calculations, we use 1D and 2D                                                                                2
axisymmetric spherical coordinates with baseline resolution                       1000                                          ε / vff0 = 0.15
∆rbase = rs0 /320 and ∆θbase = π/192, with one extra level of                                                                          2
                                                                                                                                ε / vff0 = 0.2
mesh refinement inside r = r∗ + 0.1(rs0 − r∗ ) to better resolve

                                                                     ρ / ρ1
                                                                                   100                                          rs0 = 50 km
the steep density gradient that arises in the cooling layer. We
                                                                                                                                rs0 = 75 km
do not employ a hybrid Riemann solver because we do not see                        10
the appearance of the odd-even decoupling instability (Quirk                                      γ = 4/3                       rs0 = 125 km
1994).                                                                                             H=0
   We employ a reflecting inner boundary condition at r = r∗                          1
for the sake of simplicity; we do not attempt to model the
protoneutron star (as done by Murphy & Burrows 2008) or                                    (b)                                     H=0
its contraction through a moving inner boundary [as done by                      1e+11                                             H = 0.003
Scheck et al. (2006) and Scheck et al. (2008)]. The outer                                                                          H = 0.005
boundary condition is kept fixed at r = 7rs0 , and is set by the                                                                    H = 0.007

                                                                     ρ [g cm ]
                                                                                 1e+10                                             H = 0.009
upstream flow at that position.
   To trigger convection below the shock, we introduce ran-
dom cell-to-cell velocity perturbations in vr and vθ at t = 0,                   1e+09
                                                                                                 rs0 = 75 km
with an amplitude 1% of the steady state radial velocity. To                                       γ = 4/3
study the interplay between shock oscillations and convection,
we also drop overdense shells with a given Legendre index ,                      1e+08
as done in Paper I, without random velocity perturbations.
   In order to track the residency time of the fluid in the gain                  1e+11
                                                                                            (c)                                 γ = 4/3 + α
region, we assign a mass-scalar to each spherically symmet-                                                                     γ = 4/3 − α
ric mass shell in the upstream flow. This scalar is passively                                                                    full + α
                                                                     ρ [g cm ]

advected by FLASH2.5. Through this technique, we are able                        1e+10                                          γ = 1.48 + α
to assign a “fluid" time to each element in the domain, corre-                                                                   full - α
sponding to the time at which the mass shell would cross the                                                                    γ = 1.46 − α
                                                                                 1e+09       rs0 = 125 km
instantaneous angle averaged shock position if advected from
the outer boundary at the upstream velocity:                                                       H=0
                     tF = tOB +                      .       (12)                    0.4               0.6     0.8         1         1.5          2
                                  rs (t)   θ
                                               |vr |
                                                                                                                 r / rs0
Here tOB is the time at which the fluid enters through the
outer radial boundary at r = rOB , and rs (t) θ is the angle av-       F IG . 2.— Sample initial density profiles, which are solutions to the spher-
eraged shock position. Initially, tOB = 0 and all the fluid be-       ically symmetric and time-independent flow equations. Panel (a) shows the
                                                                     zero-heating configurations for all the sequences shown in Table 1. Other pa-
low the shock is set to tF = 0. This prescription works well                                                                        ˙
                                                                     rameters are γ = 4/3, M1 (rs0 ) = 5 for all configurations, and M = 0.3M s−1 ,
for statistical studies (§4.2), tracing large scale fluid patches,    M = 1.3M for the NSE models. Panel (b) shows a sequence with a fixed
despite some inevitable turbulent mixing among neighboring           cooling function and range of heating rates (H is given in units of rs0 v3 0 ).
fluid parcels.                                                        The dashed line shows the upstream flow. Panel (c) shows a sequence with
   An explosion is defined as either i) a collision between           different equations of state, rs0 = 125 km, and H = 0. The labels “+α” and
                                                                     “−α” mean with and without α-particles included in the EOS, while “full”
the shock and the outer boundary of the simulation volume            means that the EOS explicitly includes finite-temperature and partially de-
(r = 7rs0 ) within 1000tff0 of the start of the simulation; or ii)   generate electrons, black body photons, and ideal-gas ions. All other param-
in the special case of the 1D constant-ε models, a transient         eters are the same as in (a). Only the γ = 4/3 upstream flow is shown. See
expansion that breaks a quasi-steady pattern within the same         Appendix B for further details.
timeframe. Even in the 1D simulations, very small changes
in heating rate can lead to dramatic changes in shock behav-
ior, and so this definition of explosion is good enough for our       critical value for an explosion, and a third with the largest
purposes.                                                            heating parameter that will allow a steady solution. Note that
                                                                     the shock starts out at ∼ 1.3 rs0 in the time-independent, spher-
                     2.2. Model Sequences
                                                                     ical flow solution, and quickly saturates at ∼ (1.8 − 2)rs0 in
   We choose six sequences of models, each with a range of           the 2D models with heating just below threshold for an explo-
heating parameters H ≥ 0, and each evolved both in spherical         sion. The quantity ε/v2 references the dissociation energy to
(1D) and axial (2D) symmetry. Their parameters are summa-            (twice) the kinetic energy of the upstream flow, and is the key
rized in Table 1. In each sequence, the normalization of the         free parameter determining the compression rate κ across the
cooling function is chosen so that r∗ /rs = 0.4 at zero heating.     shock (eq. [3]).
Three sequences have a constant dissociation energy, which              When examining how the prescription for nuclear dissoci-
take the values ε/v2 0 = {0.1, 0.15, 0.2}. The other three se-
                    ff                                               ation influences the results, we will focus on the ε = 0.15v2 0  ff
quences assume NSE below the shock, and have shock radii             sequence and the NSE sequence with rs0 = 75 km, which have
rs0 = {50, 75, 125} km at zero heating. The other parameters         similar initial density profiles (due to the low initial α-particle
in the NSE models are M = 1.3M and M = 0.3M s−1 .                    abundance in the NSE model).
   Table 1 samples some properties of a few models from each            The six initial models at zero heating are shown in Fig. 2a.
sequence: one with zero heating, another with H close to the         Panel (b) shows the sequence of initial models with rs0 =
6                                                                  FERNÁNDEZ & THOMPSON

                                  TABLE 1
                           S AMPLE C ONFIGURATIONS                                           0.01

     ε/v2 0
        ff        Hv−3 rs0 −1
                    ff          rs /rs0      ε/v2
                                                1        κ     χ

                                                                                ωgrow tff0
      0.1             0         1.00         0.10        7.3    0
                  8.00E-3       1.27         0.13        7.7   4.5
                  1.48E-2*      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-2*      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                            -0.01
                  8.38E-3*      2.08         0.47       34.9   74
    rs0 [km]      Hv−3 rs0 −1
                    ff          rs /rs0   ε(t = 0)/v2
                                                    1    κ     χ     Xα (rs )                                        2
                                                                                                       ε = 0.1vff0
      50              0         1.00         0.11        7.6     0   5.5E-6                                              2
                  8.00E-3       1.29         0.15        8.1   5.5   6.1E-5                  0.15      ε = 0.15vff0
                  1.43E-2*      3.01         0.26       8.6     27    0.43                             ε = 0.2vff0
      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

                                                                                ωosc tff0
                  1.15E-2*      3.61         0.21       6.8     51    0.83
      125             0         1.00         0.21       10.7     0    0.26                     0.1
                  3.50E-3       1.33         0.21        9.9    22    0.51
                  7.28E-3*      3.99         0.18       6.1    120    0.99

    * No-steady   state heating rate (e.g., Burrows & Goshy 1993).

75km and a range of heating parameters. The model with                                         0.004             0.006                           0.008    0.01
H = 0.007v3 rs0 is close to the threshold for an explosion,
             ff0                                                                                                     H / (rs0 vff0 )
while the one with H = 0.009v3 rs0 is well above thresh-
old. At higher values of H, cooling by α-particle dissocia-                         F IG . 3.— Linear growth rates (top) and oscillation frequencies (bottom)
                                                                                  of 1D models with constant dissociation energy ε, as a function of heating
tion (eq. [11]) can be significant in a layer below the shock,                     parameter H around the threshold for explosion. Stars denote configurations
causing the density profile to steepen slightly.                                   that explode within 1000tff0 . Increased heating makes the system more un-
   Fig 2c shows how our constant-γ, ideal gas approximation                       stable because the density profile flattens, akin to an increase in γ. Dotted
to the internal energy of the flow compares with the full EOS                      lines show the frequency ωosc = 2π/(2tadv ). Oscillation frequencies decrease
                                                                                  with increasing heating rate because rs moves out relative to r∗ , so that the
containing finite-temperature and partially degenerate elec-                       advection time tadv (eq. [13]) increases.
trons (see Appendix B for details). The curves labeled “+α”
include our prescription for heating/cooling by α-particle re-
combination/dissociation, and those labeled “−α” do not. We
show the sequence with the largest shock radius (rs0 = 125                        explosion in calculations by Ohnishi et al. (2006) and Mur-
km) so that NSE allows some α’s to be present. The neglect                        phy & Burrows (2008). Both calculations employed a realis-
of electron captures below the shock results in an adiabatic                      tic EOS, but like us included neutrino heating as a local source
index between 4/3 and 5/3 in the zone where α-particles are                       term in the energy equation. Oscillations have also been seen
absent. This causes the EOS to stiffen, so that the density pro-                  by Buras et al. (2006b) in more elaborate calculations with
file is well approximated by an ideal gas with γ 1.48 at zero                      Boltzmann neutrino transport.
heating. Adding in heating tends to flatten the density profile                        The origin of the 1D SASI oscillation can be briefly summa-
even more, and with γ = 1.48 it would be much flatter than                         rized as follows. An initial outward shock displacement gen-
is typically seen in a realistic core collapse model. Hence we                    erates an entropy perturbation, which is negative for γ 5/3.
choose an EOS with γ = 4/3.                                                       This entropy perturbation is advected down to the cooling
                                                                                  layer, where it causes, at constant ambient pressure, an in-
                   3. ONE DIMENSIONAL SIMULATIONS                                 crease in the cooling rate, δLC /LC = −[(γ − 1)/γ](b − a)δS >
       3.1. Shock Oscillations and Transition to Explosion                        0. The resulting negative pressure perturbation is rapidly
                                                                                  communicated to the shock, which recedes and generates an
   An explosion in spherical symmetry involves the develop-                       entropy perturbation of the opposing sign. One more iteration
ment of an unstable = 0 SASI mode. We showed in Paper                             results in a shock displacement of the same sign as the initial
I that, in the absence of neutrino heating, the period of this                    displacement, and allows the cycle to close. The duration of
mode is essentially twice the post-shock advection time. As                       the = 0 mode is therefore nearly twice the advection time
heating is introduced into the flow, we find that this relation is                  from the shock to the cooling layer,
maintained. The = 0 mode is damped until the heating rate
is pushed above a critical value, which we now discuss.                                                        2π                       dr
   It should be emphasized that this critical heating rate is gen-                                                           2               .            (13)
                                                                                                               ωosc              r∗    |vr |
erally lower than that defined by Burrows & Goshy (1993),
which marked the disappearance of a steady, spherically sym-                      The cycle is stable for γ = 4/3 and r∗ /rs0 = 0.4.
metric solution to the flow equations. Large amplitude 1D                            When heating is added, the density profile flattens. This is
shock oscillations have been witnessed near the threshold for                     analog to an increase in γ, which increases the linear growth
                                          ALPHA PARTICLE RECOMBINATION IN SUPERNOVAE                                                                              7

           6        0.92 Hcr                                              2
                                                          ε = 0.15 vff0
                    0.95 Hcr
                    0.98 Hcr
                    1.02 Hcr
           4        1.05 Hcr
                    1.08 Hcr

rs / rs0

                   0.91 Hcr
           6                                                   NSE
                   0.94 Hcr
                   0.98 Hcr                              rs0= 75 km
                   1.02 Hcr
                   1.06 Hcr
           4       1.09 Hcr
                                                              rα Μ1.3


               0    200          400          600          800                1000
                                      t / tff0
   F IG . 4.— Shock radius as a function of time for two sequences of 1D
 simulations. The upper panel shows runs with constant dissociation en-
 ergy ε = 0.15v2 0 , and a range of heating coefficients H near the critical
 value Hcr = (0.006625 ± 0.000125)v3 0 rs0 . The lower panel shows runs with
 rs0 = 75 km and an α-particle contribution to the EOS. In this case, Xα is
 initially negligible everywhere below the shock (see Table 1), but grows as
 the shock expands. The horizontal dotted line labels the radius rα at which
 the nuclear binding energy Qα of an α-particle equals its gravitational bind-
 ing energy (eq. 1). The critical heating for this second sequence is lower,
 Hcr = (0.006125 ± 0.000125)v3 0 rs0 .

 rates and thus pushes the system towards instability (e.g., Pa-
 per I). There is a critical heating rate for which the damping                        F IG . 5.— Radial profiles of various quantities during shock breakout in the
 effect of the 1D SASI is neutralized and there is no net growth.                    NSE model with H = 1.09Hcr and rs0 = 75 km (Fig. 4). Top panel: mass frac-
                                                                                     tion of α-particles. Second panel: rate of release of specific nuclear binding
 We find that, once the heating rate exceeds this critical value,                     energy denuc /dt compared with the (adiabatic) rate of change of enthalpy wad
 the system always explodes.                                                         [eq. 14]. Third panel: net neutrino heating rate per unit volume LH − LC
    We therefore define the critical heating rate in our 1D simu-                     (thin curves) and denuc /dt (thick curves), both normalized to the local value
 lations to be the minimum heating rate for growing shock os-                        of c2 = γ p/ρ. Bottom panel: radial velocity normalized to vff0 at radius rs0 .
                                                                                     Both denuc /dt and wad are smoothed in radius for clarity.
 cillations.5 Figure 3 shows linear eigenfrequencies as a func-
 tion of heating rate for our 1D initial configurations with con-
 stant dissociation. These values were obtained by solving the
 differential system of Foglizzo et al. (2007), modified to ac-                       feedback between the shock and the cooling layer is broken.
 count for a constant rate of nuclear dissociation (Paper I) as                      Material then tends to pile up in the gain region, is further
 well as incorporating the heating function in eq. (10). The                         heated, and more material reverses direction. The net effect
 runs marked by stars explode within a time 1000tff0 , and so                        is to push the shock outward. Movie 1 in the online material
 require a small, but finite, positive growth rate.                                   illustrates this chain of events.
    In an exploding run, the expansions become longer and con-
 tractions shorter as the shock oscillation develops a large am-                               3.2. Effects of Alpha-Particle Recombination
 plitude. Eventually the accretion flow is halted during a con-                          Shock breakout is controlled by the build-up of positive en-
 traction. This marks the point of explosion, beyond which the                       ergy fluid downstream of the shock, and therefore is sensitive
    5 We define our critical heating parameter H to be the average of the val-        to the density profile below the shock. Heating by neutrinos
 ues in the exploding and non-exploding runs that are closest to the threshold       is concentrated fairly close to the protoneutron star, inside a
 for explosion, within our fiducial 1000tff0 cutoff.                                  distance ∼ (2 − 3)r∗ . Heating by α-particle recombination is
8                                                                    FERNÁNDEZ & THOMPSON
           6                                ε = 0.15 vff0
                                                                       1.06 Hcr                                         ε = 0.15 vff0
                                                                       1.02 Hcr
                                                                       0.98 Hcr                                      rs,max / |drs,max/ dt|
                                                                       0.91 Hcr                                      〈 tres 〉 vol (upper 10%)
           4                                                                                                         〈 tres 〉 vol (upper 50%)

           2                                                                                              50

                                             rs0 = 75 km                1.06 Hcr
                                                                        1.02 Hcr
                                                                        0.98 Hcr                         150             rs0 = 75 km
                                                                        0.90 Hcr
                                                                     rα M1.3

           2                                                                                              50
                                                                                           time [tff0]
rs / rs0

                                              rs0 = 50 km
           6                                                                   -1
                                                                      rα M1.3                                            rs0 = 50 km
                                                                       1.08 Hcr
           4                                                           1.02 Hcr
                                                                       0.98 Hcr                          100
                                                                       0.95 Hcr



           6                              rs0 = 125 km                 1.19 Hcr
                                                                       1.04 Hcr
                                                                       0.96 Hcr
                                                                                                                          rs0 = 125 km
                                                                       0.89 Hcr                          300

                                                                -1                                       200
                                                    rα M1.3

           0                                                                                               0
               0      200           400          600                 800            1000                       0        200                        400
                                         t / tff0                                                                             t / tff0
  F IG . 6.— Left panel: Angle-averaged shock radius (solid lines) and maximum shock radius (dotted lines) for various 2D models around the threshold for
explosion. Upper panel shows runs with constant dissociation energy ε = 0.15v2 0 , while lower panels displays NSE runs with rs0 as labeled. Critical heating
rates Hcr are different for each configuration, and can be found in Figure 15. Right panel: Black lines show expansion timescale of maximum shock radius
texp ∼ rs,max /|drs,max /dt|, computed using a polynomial fit for the runs just above the threshold for explosion (corresponding to the grey solid lines on the right
panels). Green and red lines show the average residency time over the 50% and 10% of the gain region volume with highest tres , respectively (see §4.2 for the
definition of this timescale). Shock breakout occurs whenever texp ∼ tres vol , except in the model where recombination heating is dominant (rs0 = 125 km).
                                  ALPHA PARTICLE RECOMBINATION IN SUPERNOVAE                                                       9

concentrated at a greater distance ∼ rα (eq. [1]), but still can    that the convective instability is driven by the negative entropy
reach a comparable amplitude.                                       gradient within the layer of maximal neutrino heating. In non-
  The dependence of shock breakout on heating rate is dis-          exploding runs, the shock settles to a quasi-equilibrium state
played in Fig. 4 for two accretion models and several values        with oscillations taking place over a range of angular (Legen-
of H close to Hcr (see Table 1). The initial expansion of the       dre) index and of temporal frequencies, as previously seen
shock during the explosion phase is very similar for models         by Ohnishi et al. (2006), Scheck et al. (2008), and Murphy
with constant ε and with NSE in the shocked fluid. However,          & Burrows (2008). The amplitude of the = 1 and 2 modes
the time evolution bifurcates near the radius rα .                  remains small until the heating parameter H has begun to ex-
  Figure 5 shows successive profiles of the shocked flow in           ceed about one half the critical value for an explosion. The
the exploding run with H = 1.09Hcr and rs0 = 75 km. The             competition between SASI growth and convective instability
α-particle fraction approaches unity as the shock reaches the       is examined in detail in §6.
radius rα . The second panel shows the specific nuclear energy          One gains considerable insight into the mechanism driving
generation rate [eq. (7)] normalized to the adiabatic rate of       shock breakout by examining the distribution of Bernoulli pa-
change of the enthalpy,                                             rameter (eq. [2]) in the shocked fluid. We first consider the
                          1 dp                                      NSE runs with rs0 = 50 km and 125 km, with the heating pa-
                    wad =       = −c2 · v.
                                     s                      (14)    rameter H just above the threshold for an explosion. Two
                          ρ dt                                      snapshots from each of these runs are shown in Fig. 7. In the
Here cs = (γ p/ρ)1/2 is the sound speed. The third panel com-       first case, the initial equilibrium shock radius is only ∼ rα /4
pares the amplitude and distribution of neutrino and recombi-       km, and α-particles are essentially absent below the shock. In
nation heating, and the bottom panel plots the radial velocity      the second, the shock starts at ∼ 2rα /3 and Xα ∼ 0.5 initially
in the postshock region.                                            in the postshock flow.
   We can summarize this behavior as follows: during the ini-          Large deformations of the outer shock are generally caused
tial expansion phase, fluid below the shock continues to move        by convective plumes that carry positive energy. Strong neu-
inward, and the dissociation of α-particles removes energy          trino heating is generally concentrated inside the inner bound
from the flow (as expected from eq. [11]). Some fluid behind          (b < 0) zone. The degree of symmetry of the bound material
the shock begins to move outward around 300tff0 , but nuclear       depends on the α-particle abundance. In the rs0 = 125 km run,
dissociation still causes a net loss of internal energy. However,   the bound region is spherically stratified and the material with
the recombination of α-particles sets in above rα , especially      b > 0 is generally excluded from it. Strong recombination
in regions where Xα 0.5. By the time the shock hits the             heating is present both below and above the surface where
outer boundary, denuc /dt exceeds one-half of |wad |.               b 0, indicating that it is mainly responsible for imparting
   The dependence of the density contrast κ (eq. [3]) on radius     positive energy to the shocked material. The mean shock ra-
also has an influence on the details of breakout. When the           dius expands by a factor ∼ 2.5 between the upper two frames
dissociation energy ε is held fixed, κ increases toward larger       in Fig. 7, but the growth in the volume of positive-energy ma-
radius. This has the effect of creating a dense layer of fluid       terial is not accompanied by a significant expansion of the in-
below the shock when rs is such that ε ∼ v2 /2. In spheri-
                                                                    ner bound region, whose outer radius remains fixed at r rα .
cal symmetry, the breakout of the shock is then impeded by             This segregation of bound from unbound material is bro-
the accumulation of a dense layer of material that cannot ex-       ken when the shock is more compact initially, as is seen in the
change position with the lighter material below it. It can hap-     lower two panels of Fig. 7. A single dominant accretion plume
pen that the energy in the expanding region is no longer able to    is continuously present, which funnels cold and dense mate-
sustain the heavier material above, and the shock collapses, as     rial into the zone of strong neutrino heating. Alpha-particles
shown in Fig. 4 for the constant ε run with H = 1.08Hcr . This      are present only well outside the boundary between b < 0 and
obstruction is avoided when statistical equilibrium between n,      b > 0. The competition between α-particle and neutrino heat-
p, and α is maintained below the shock, because ε/v2 and   ff
                                                                    ing is discussed in more detail in §4.1, and the influence of
κ both decrease gradually with radius (Fig. 1). Obviously,          α-particles on the threshold heating rate for an explosion is
this limit to the shock expansion does not occur in 2D, as this     examined in §5.
situation would be Rayleigh-Taylor unstable.                           The accumulation of a bubble of hot material right behind
                                                                    the shock is a consequence of the balance of the buoyancy
             4. TWO-DIMENSIONAL SIMULATIONS                         force acting within the bubble, and the ram pressure of the
  Extending the flow calculation to two dimensions reveals           preshock material. The ratio of force densities is (Thompson
some subtle patterns of behavior. The time evolution of the         2000)
shock is shown in Fig. 6a for a range of heating rates near               Fbuoy      ρ − ρbubble    2GM/rs
the threshold for explosion. In contrast with the 1D runs, the                                                  ∆Ωbubble ,      (15)
breakout of the shock looks similar in models with constant               Fram            ρ           v2
dissociation energy and with NSE between n, p, and α below          where rs is the shock radius, vr is the ambient radial flow
the shock. Both types of models are subject to buoyancy-            speed, ρ is the ambient density, ρbubble the density of the bub-
driven instabilities, which allow cold material below the shock     ble, and ∆Ωbubble is its angular size. A low-density bubble
to interchange position with hotter material within the gain        (ρ − ρbubble ∼ ρ) can resist being entrained by the convective
region. As a result, the shock is highly asymmetric at breakout     flow once it grows to a size ∆Ωbubble ∼ M2 Sr, where Mcon
in both cases. In §5 we compare the critical heating rate for       is the convective Mach number. On the other hand, the bubble
explosion in 1D and 2D, and examine how it is influenced by          must attain a much larger angular size ∆Ωbubble ∼ 1 Sr if the
α-particle recombination.                                           buoyancy force is to overcome the upstream ram (|vr | ∼ vff )
  Around the threshold for explosion, all of our runs develop       and force a significant expansion of the shock surface. Fig-
vigorous convective motions before the SASI has a chance to         ure 7 shows that the extent of the shock expansion is indeed
undergo even a few oscillations. At high heating rates, we find      correlated with the angular width of the region where hot ma-
10                                                            FERNÁNDEZ & THOMPSON

  F IG . 7.— Snapshots of two separate models, with heating parameter H just above the threshold for an explosion, and initial shock radius either well inside rα
(rs0 = 50 km, without heating) or close to rα (rs0 = 125 km, without heating). Within each panel, the top figure displays Bernoulli parameter b; the middle figure
the rate of change of nuclear binding energy per unit mass; and the bottom figure the net neutrino heating rate per unit mass. Top left: the early development of
an asymmetric plume with positive b; top right: the same run just before the shock hits the outer boundary. In this rs0 = 125 km run, the heating by α-particle
recombination is enhanced with respect to neutrino heating due to the large Xα in the initial stationary model. Recombination heating straddles the central zone
with b < 0, which maintains a nearly spherical boundary near the radius rα ( 2.0 rs0 ). Bottom left: α-particles begin to form as the shock approaches rα in
the rs0 = 50 km run, but neutrino heating remains much stronger than recombination heating. Bottom right: the same run just before the shock hits the outer
boundary. When the shock starts off well inside rα , neutrino heating dominates the initial expansion, and material with b > 0 forms well inside rα (see §4.3).
Animations showing the evolution of these two configurations are available in the online version of the article.
                                        ALPHA PARTICLE RECOMBINATION IN SUPERNOVAE                                                            11

                                                                                terial accumulates.
                                                                                   Another interesting feature of Fig. 7 is the presence of sec-
                                                                                ondary shocks, which are triggered once the outer shock be-
                                                                                comes significantly non-spherical. Their locations are marked
                                                                                by discrete jumps in the rate of recombination heating. Sec-
                                                                                ondary shocks are also prevalent throughout the nonlinear
                                                                                phase in the constant-ε models. Figure 8 shows the nor-
                                                                                malized 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 simu-
                                                                                lation volume). The online version of the article contains an
                                                                                animated version of Fig. 8 showing the complete evolution.
                                                                                In both cases, secondary shocks extend over the whole post-
                                                                                shock domain, signaling the dissipation of supersonic turbu-
                                                                                lence which is stirred by accretion plumes that penetrate into
                                                                                the gain region.

                                                                                 4.1. Distribution of Alpha Particle Recombination Heating
                                                                                   Heat input by neutrino absorption and by α-particle recom-
                                                                                bination have very different distributions within the shocked
 F IG . 8.— Normalized pressure gradient (r/p)| p| showing the secondary        fluid: strong neutrino heating is concentrated inside rs0 ,
shock structure during breakout. The top model is ε = 0.15v2 0 and the bottom
                                                           ff                   whereas recombination heating of a comparable amplitude
NSE with rs0 = 75 km. Both have heating rates just above the threshold for
an explosion. An animation showing the time evolution is available in the       is distributed throughout the settling flow. Strong recombi-
online version of the article.                                                  nation heating quite naturally extends below the zone where
                                                                                α-particles are present in significant numbers, as is seen in
                                                                                Fig. 9. The first and third panels of this figure depict the
                                                                                pre-explosion steady state of the rs0 = 75 km model with
                                                                                H = 1.02Hcr , 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 1D simulations (e.g.
                                                                                Fig. 5), Xα approaches unity during shock breakout.
                                                                                   The relative strength of neutrino heating and recombina-
                                                                                tion heating depends on the initial radius of the shock, and on
                                                                                the Bernoulli parameter of the postshock material. Figure 10
                                                                                separates out cooling by α-particle dissociation from heating
                                                                                by recombination. The sharp negative spike near b = 0 rep-
                                                                                resents α-particle dissociation in fresh, cold downflows. The
                                                                                formation of material with b > 0 is primarily due to α-particle
                                                                                recombination in the rs0 = 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.
                                                                                   The strength of the boost given to the shock by recombina-
                                                                                tion heating can be gauged by comparing denuc /dt to the adia-
                                                                                batic rate of change wad of the enthalpy of the flow (eq. [14]).
                                                                                Figure 11 shows the result for all three NSE sequences with
                                                                                H just above Hcr . In all cases, denuc /dt wad in various parts
                                                                                of the shocked fluid once the shock extends beyond a radius
                                                                                   rα , just as in 1D. Most of the heat input by recombination
                                                                                is concentrated where Xα ∼ 0 − 0.5.
                                                                                   Histograms of denuc /dt versus b and Xα are shown in
                                                                                Fig. 12. The rapid dissociation of α-particles in fresh down-
                                                                                flows is represented by the long tail toward large negative val-
                                                                                ues of denuc /dt. The α-particle concentration is very stratified,
                                                                                with higher Xα occurring at larger radius.

                                                                                                     4.2. Residency Time
  F IG . 9.— Top panels: rate of release of specific nuclear binding energy
denuc /dt. Bottom panels: mass fraction of α-particles Xα . We show two
                                                                                  A long residency time of material in the gain region is com-
instants in the exploding NSE run with rs0 = 75 km and H = 1.02Hcr . The        monly viewed as a key ingredient in a successful neutrino-
shock contour is approximated by the white line which marks XO = 90%.           driven explosion. To calculate tres , we use the method de-
                                                                                scribed in §2.1.3: we first assign a unique “fluid time" tF
                                                                                (eq. [12]) to each infalling radial mass shell in the simulation,
12                                                                            FERNÁNDEZ & THOMPSON

                              0.5                    rs0 = 50 km


                                      t = 250 tff                       t = 350 tff                 t = 415 tff                    t = 455 tff
    dE / db [ρ1 vff,0 rs0 ]

                               0.5                   rs0 = 75 km


.                                    t = 200 tff                       t = 290 tff                  t = 370 tff                    t = 434 tff

                               0.5                  rs0 = 125 km


                                      t = 100 tff                      t = 200 tff                   t = 270 tff                   t = 333 tff

                                     -0.4   -0.2        0        0.2   -0.4   -0.2      0    0.2    -0.4   -0.2     0      0.2    -0.4    -0.2       0      0.2
                                                                                             b [vff,0 ]
 F IG . 10.— Distribution of volume-integrated heating rate of material, with respect to Bernoulli parameter b. We restrict attention to material in the gain region
(defined by LH > LC ) in the three 2D NSE runs just above the threshold for explosion. Black curves: net heating rate resulting from neutrino absorption and
emission. Red/green curves: heating/cooling rate by α-particle recombination and dissociation in material with denuc /dt > 0 and denuc /dt < 0, respectively.
Blue curves: net heating/cooling rate due to changing α-particle abundance. 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 rs0 = 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.

which is effectively the time at which it passes the shock. We                                material comprising the upper part of the residency time dis-
then define6 the residency time of the fluid as                                                 tribution). The final breakout of the shock seen in Fig. 6a
                                                                                              corresponds to the time when texp ∼ tres vol . The breakout is
                                               tres = t − tF ,                        (16)
                                                                                              a bit more gradual in the rs0 = 125 km model with heating just
where t is the present time. A related method (tracer parti-                                  above threshold for an explosion (H = 1.04Hcr ; see animated
cles) is used by Murphy & Burrows (2008) to calculate the                                     version of Figure 7a,b in the online material). In this case, the
residency time in collapse simulations with a more realistic                                  expansion time of the shock remains somewhat longer than
EOS.                                                                                          the residency time of material below the shock, which implies
  As material with positive Bernoulli parameter accumulates                                   that the breakout depends on the continuing release of nuclear
below the shock, we indeed find that its tres grows larger.                                    binding energy as well as on the separation of positive-energy
The shock starts running outward if the energy of this un-                                    fluid from the accretion flow.
bound material grows on a timescale shorter than the con-                                        This lengthening of the mean residency time can largely
vective time. Figure 6b shows the characteristic expansion                                    be ascribed to the increased dynamical time of the expanding
time of the shock texp ∼ rs,max /|drs,max /dt| (as measured at its                            shock. What changes most dramatically during breakout is
outermost radius), alongside tres vol (as measured within the                                 the ratio of the expansion time to the dynamical time. Note
                                                                                              that large changes in the distribution of tres are concentrated in
   6 Since t is defined in terms of the angle-averaged radius of the shock,
            F                                                                                 regions of positive b. In Fig.12, we plot the distribution of b
there is a modest error in tF due to non-radial deformations of the shock.
Given the lack of substantial large-scale mixing between the single accre-
                                                                                              and tres in the ε = 0.15v2 0 run that is just above the threshold
tion funnel and convective cells, this prescription serves well as a tracer of                for an explosion (see Fig. 13).
different fluid populations.
                                          ALPHA PARTICLE RECOMBINATION IN SUPERNOVAE                                                                               13

                                                                                                                Mass-weighted Radius [rs0]

                                                                                                        0.41     0.60      1.00    1.50       2.00    3.50
                                                                                                 0.2                       _
                                                                                                                  t = (119 + 25)                           _
                                                                                                                                                  t = (189 + 25)

                                                                                   b [vff02]

                                                                                                 0.2                       _
                                                                                                                  t = (264 + 25)                           _
                                                                                                                                                  t = (364 + 25)

                                                                                   b [vff02]
                                                                                                    -20 0      20 40 60 80                0   20 40 60 80 100
                                                                                                                 tres [tff0]                    tres [tff0]

                                                                                                                Mass-weighted Radius [rs0]

  F IG . 11.— Ratio of denuc /dt (rate of release of specific nuclear binding                            0.41     0.60      1.00    1.50       2.00    3.50
energy, eq. [7]) to wad (adiabatic rate of change of the enthalpy, eq. [14]) for
the three NSE models with H just above Hcr .                                                     0.6
                                                                                                        rs0 = 75 km, t = 434 tff
                                                                                   b [vff02]

   It is also apparent from Fig. 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. Before the shock
gets close to a radius rα , most of the shocked material has
negative Bernoulli parameter and tres 50tff0 , pointing to ma-
terial lingering for about one convective overturn.
       4.3. Heat Engine in a Two-Dimensional Explosion                                            80
                                                                                   tres [tff0]

   Fresh material that is accreted through an oblique shock has                                   60
a relatively low entropy, but once it reaches the base of the                                     40
gain region it is exposed to an intense flux of electron-type                                      20
neutrinos and is heated. A buoyancy force will push some                                           0
of this heated material upward, and force an overturn of the                                     −20
fluid below the shock. Material with a longer residency time                                        0.0 0.2 0.4 0.6 0.8                       -0.2       0.0        0.2
may therefore undergo multiple episodes of heating. On this                                                   Xα                          denuc / dt [vff03 / rs0]
basis, Herant et al. (1992, 1994) suggested that a convective
flow would mediate a heat engine below the shock that would                           F IG . 12.— Top panel: Histogram of Bernoulli parameter b and residency
drive a secular increase in the energy of the shocked fluid.                        time tres in the exploding run with ε = 0.15v2 0 (see also Fig. 13). The colors
   We now investigate whether a heat engine operates in our                        label the mass weighted radius, and we include all material experiencing a
                                                                                   net excess of neutrino heating over cooling. Bottom panel: Histogram of
simulations, and how it depends on the heating parameter H.                        Bernoulli parameter b and residency time tres versus α-particle mass fraction
We focus on a model with a constant nuclear dissociation en-                       Xα and rate of release of specific nuclear binding energy denuc /dt, in the
ergy, ε = 0.15v2 0 . In this class of models, the infalling heavy
                                                                                   exploding run with rs0 = 75 km (Fig. 9).
nuclei are completely broken up below the shock, and no heat-
ing by the reassembly of α-particles is allowed. As a first step,
we average the convective flow over windows of width 50tff 0 ,
which de-emphasizes short term fluctuations in the averaged                         fluid is lifted by buoyancy forces and accumulates in the re-
velocity field v t . Figure 13 shows v t and b t (eq. [2]) at                       gion in between the top of convective cells and the shock.
four different times in the run with H just above the threshold                    A strong deformation of the shock allows a plume of fresh
for explosion (H = 1.02Hcr ). The radius of maximum heating                        material to descend diagonally between the convective cells.
(r 0.66rs0 ) coincides with the lower boundary of the convec-                      The tilt of this cold downflow intermittently flips in sign, and
tive cells, across which material flows horizontally. Heated                        the averaged circulation pattern typically has an “∞" shape.
14                                                                               FERNÁNDEZ & THOMPSON

 F IG . 13.— Four snapshots of the exploding model with ε = 0.15v2 0 and H = 1.02Hcr . The Bernoulli parameter (color map) and the velocity field (white arrows)
are averaged over intervals of duration 50tff 0 . The thick white contours show the surface with 50% mass fraction in heavy nuclei (time-averaged).

                                             |〈 b〉 t |               Θ                      The heating of fluid parcels in the two hemispheres is also
                           0.4                                                              intermittent, and sometimes two circulation flows are estab-
                                     H = 1.02 Hcr                t = (119± 25)   tff0       lished simultaneously, thereby causing a bipolar expansion of
                                                                 t = (189± 25)   tff0       the shock.
 specific energy [vff0 ]

                           0.3                                   t = (264± 25)   tff0          We have identified a useful figure of merit which connects a
                                                                 t = (364± 25)   tff0       secular increase in the shock radius to the strength of neutrino
                                                                                            heating at the base of the gain region. Figure 14 shows the
                           0.2                                                              absolute value of the Bernoulli parameter b at the radius of
                                                                                            maximum heating (rH,max 0.66rs0 ) as a function of polar
                                                                                            angle θ. In the top panel, the four sets of lines correspond
                           0.1                                                              to the four snapshots of Fig. 13, and the bottom panel shows
                                                                                            the analogous results for a non-exploding run. Overplotted as
                                                                                            dashed lines is the quantity

                                                                                                                    LH − LC                       rH,max
                                                                                                           Θ=                                               .         (17)
                                     H = 0.91 Hcr                t = (200± 25) tff0                                    ρ          t,θ∗,r=rH,max | vθ t,θ∗ |
 specific energy [vff0 ]

                                                                 t = (500± 25) tff0

                                                                 t = (975± 25) tff0         This measures the specific energy that is absorbed from neu-
                                                                                            trinos by the material that flows laterally along the lower
                           0.2                                                              boundary of the convective cells. In eq. (17), the angular av-
                                                                                            erage of the heating rate and meridional velocity is restricted
                                                                                            to a single convective cell.7
                           0.1                                                                 In non-exploding 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 pa-
                                 0                  1        2                          3   rameter 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
  F IG . 14.— Bernoulli parameter (solid lines) and specific energy absorbed
during lateral advection Θ (eq. [17], dashed line) at the radius of maxi-                      7 To identify the range of angles comprising the lower boundary of a con-
mum heating rH,max . These quantities are averaged over intervals of dura-                  vective cell (r = rH,max ), we first average vθ t over all polar angles, and then
tion 50tff 0 . The upper panel shows the exploding run with ε = 0.15v2 0 and
                                                                                            define a single cell as a zone where | vr t | < | vθ t,θ | and vθ has the same
H = 1.02Hcr (the same as in Fig. 13), and the lower panel shows the non-                    sign. Once the convective cells have been so identified, the angular average
                              2 and H = 0.91H .
exploding run with ε = 0.15vff 0                cr
                                                                                            is repeated within each cell. The quantity | vθ t,θ∗ | appearing in eq. (17)
                                                                                            represents this more restricted average, which typically covers ∼ 1 rad in the
                                                                                            polar direction (e.g., Fig. 13).
                                                        ALPHA PARTICLE RECOMBINATION IN SUPERNOVAE                                                                         15
                      0.01                                                                                                   0.02

                                                                              1D constant ε
                                                                              2D constant ε                                                                    1D
                     0.008                                                    1D NSE                                                                           2D

                                                                                                 M1.3 Hcr / [rα vff(rα) ]
                                                                              2D NSE                                        0.015
 Hcr / (rs0 vff0 )




                         0.1              0.15            0.2               0.25                                                0.2   0.3   0.4          0.5    0.6
                                                        ε / v1                                                                               M1.3 (rs / rα)
  F IG . 15.— Critical heating parameter Hcr that yields an explosion, for                        F IG . 16.— Critical heating parameter Hcr that yields an explosion, for
all the model sequences explored in this paper (Table 1). The abscissa is                       the runs that include α-particles in the EOS. The abscissa is the ratio of the
the ratio of ε to v2 in the initial flow configuration (v1 being the radial flow
                   1                                                                            initial shock radius rs to rα . Error bars have the same meaning as in Fig. 15.
velocity just upstream of the shock). Error bars show the separation between                    The critical heating parameter (a close analog of Lν ) decreases substantially
exploding and non-exploding models, with the points marking the average.                        with increasing shock radius. The differences in Hcr between the 1D and 2D
                                                                                                models also decreases.

panel of Fig. 14, Θ grows with time as the system approaches
the explosion.
   The stability of the averaged flow pattern appears to be, in                                  and v1 is the flow speed upstream of the shock in the initial
part, an artifact of the axisymmetry of the flow. This im-                                       configuration (that is, in the time-independent, spherical flow
poses strong restrictions on the motion of convective cells,                                    solution). In the case of the NSE equation of state, this quan-
causing vorticity to accumulate on the largest spatial scales.                                  tity can be translated into an initial value of the shock radius
Our observation that the bulk of the neutrino heating takes                                     using Fig. 1. (Note that ε/v2 has a weak dependence on rs .)
place within horizontal flows suggests that the ratio of heat-                                      A few interesting features of Fig. 15 deserve comment.
ing timescale to radial advection time in the gain layer may be                                 First, a comparison with Table 1 shows that the critical heat-
a less precise diagnostic of the conditions for explosion in 2D:                                ing rate for explosion is ∼ 50 − 70% of the maximum heat-
the horizontal convective velocity is typically low compared                                    ing rate for which a steady-state flow solution can be found.
with the downward velocity of the main accretion plume. We                                      The maximal heating parameter Hsteady for a steady flow cor-
do observe that the main accretion plume becomes strongly                                       responds directly to the one first determined by Burrows &
distorted near the threshold for an explosion, so that a signif-                                Goshy (1993) using a more realistic EOS. Note also that the
icant fraction of the plume material enters one of the convec-                                  values of Hcr in the 1D and 2D models are much closer to
tion cells. This effectively decreases the amount of material                                   each other than they are to Hsteady . This result is perhaps not
that accretes to the protoneutron star and thus increases the                                   surprising, given that the explosion is not immediate, but is
overall advection timescale across the gain region.                                             approached through a series of transient fluid motions.
                                                                                                   Second, Hcr is lower when NSE between n, p and α is
                         5. CRITICAL HEATING RATE FOR EXPLOSION
                                                                                                maintained below the shock. The difference between the NSE
   An explosion occurs when the heating parameter H is raised                                   models and the constant-ε models is only ∼ 10% in Hcr when
above a critical value8 Hcr . We now explore how Hcr de-                                        the shock starts out well below rα (eq. [1]). In spherical sym-
pends on the details of the EOS and the initial radius of                                       metry, this effect can be ascribed to the slow decline of ε with
the shock. One can express H simply in terms of the ra-                                         increasing shock radius: this makes it easier for the shock to
tio of the heating rate (4πr3 LH ) to the accretion luminosity                                  expand through the range of radii where ε/v2 0 approaches 2 .
(GM M/r), in the idealized (but unrealistic) case where the                                     An explosion is significantly easier when the fluid below the
flow is composed only of free nucleons and moves hyperson-                                       shock starts out with a significant population of α particles, as
ically. Then this ratio depends on H but not on the accretion                                   in the models with rs0 = 125 km.
rate M = 4πr2 ρ(r)|vff (r)|. The precise value of the reference                                    Third, Hcr tends to decrease with increasing ε/v2 : a slightly
radius is unimportant; we choose rs0 , the shock radius in the                                  lower heating rate per unit mass is required to explode a
time-independent, spherical flow solution at H = 0. Then                                         flow with a larger density contrast κ across the shock. Be-
                                                                                                cause almost all the gravitating mass is in the collapsed core,
                             4πrs0 3 LH [ρ1 (rs0 )]                        2H                   the gravitational binding energy of the gain region is ap-
                                                                     =            .      (18)
                                 GM M/rs0                                rs0 v3 0
                                                                                                proximately proportional to κ, whereas the net heat absorbed
                                                      Xα =0; M=∞                                over the advection time is a stronger function of density,
This quantity is ∼ 10−3 − 10−2 in the models we examine,                                        tadv (LH − LC )d3 r ∝ κ2 . (One factor of κ comes from the
which are below or near the threshold for explosion.                                            advection time tadv as given by eq. [13], and the other from the
  Figure 15 displays Hcr for all of our model sequences. The                                    density dependence of LH .) For example, Table 1 shows that
abscissa is ε/v2 , where ε is the nuclear dissociation energy
                                                                                                κ is ∼ 1.6 times larger for ε/v2 0 = 0.2 than for ε/v2 0 = 0.1,
                                                                                                                                 ff                     ff
                                                                                                and that Hcr is smaller by the inverse of the same factor.
       8        Our method for determining Hcr is discussed in §3.1.                               Fourth, the 2D runs all require less heating than their 1D
16                                                           FERNÁNDEZ & THOMPSON

counterparts to explode. A major reason for this is that all
two-dimensional configurations explode along one or both
poles (see Figs. 11 and 13), so that less material must be lifted
through the gravitational field than in a fully spherical explo-
   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 Hcr in the
NSE models grows rapidly as the initial shock radius9 rs 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 neu-
trino 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
(2007) 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
rs0 /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.
   One also notices from Fig. 16 that the difference between
Hcr in 1D and 2D depends on the starting radius of the shock.
The closer rs0 is to rα , the weaker the dependence of the crit-
ical heating rate on the dimensionality of the flow.
                   6. CONVECTION AND THE SASI
   Overturns of the fluid below the shock can be triggered in
two distinct ways: through the development of Ledoux con-                        F IG . 17.— The development of a convective instability is strongly lim-
vection in the presence of a strong negative entropy gradient,                 ited when the parameter χ 3. These panels show snapshots of entropy
or via the non-linear development of the SASI (a linear feed-                  (normalized to initial postshock value) in a NSE run (rs0 = 75 km) with two
                                                                               different heating rates. When χ = 2.1, convective cells of a limited extent are
back between ingoing entropy and vortex waves, and an out-                     triggered in the layer where the net heating rate is strongest, but they do not
going sound wave). We now show that the amplitude of the                       propagate into the upper parts of the gain region. Convection becomes much
dipolar mode that is excited in the shock is strongly tied to                  more vigorous and widespread when χ = 5.5. Note that both of these models
the level of neutrino heating, and so depends on a balance be-                 are non-exploding. An animation showing the time evolution of these two
                                                                               configurations is available in the online version of the article.
tween non-linear excitation and damping mechanisms. 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 bi-                   the threshold for a neutrino-driven explosion. The implica-
furcates from older shocked fluid. Nonetheless, the origin of                   tion for convection below the shock is illustrated in Fig. 17,
the convective motions does have implications for the stabil-                  which shows two snapshots for models with χ = 2.1 and 5.5,
ity of = 1 and 2 modes in 3D simulations: one expects that                     neither of which explodes. At the lower heating rate, the time
large scale oscillations will change shape and direction more                  required for convection to develop depends on the strength of
rapidly if they are triggered primarily by neutrino heating.                   the seed perturbation, whereas at the higher heating rate con-
   We can ask whether a heating parameter H that yields an ex-                 vective overturns develop rapidly within the layer of strong
plosion will also form an unstable entropy gradient below the                  neutrino heating and spread throughout the post-shock region
shock. Convection develops through a competition between                       over a few dozen dynamical times. (The figure shows the
inward advection and neutrino heating. A detailed analysis                     result in the case where the seed perturbation is dominated
by Foglizzo et al. (2006) shows that the parameter                             by a small spherical startup error in the initial model.) We
                                 ωBV                                           conclude that, near the threshold for a neutrino-driven explo-
                        χ≡              dr,                (19)                sion and for our given set of physical assumptions, convec-
                                                                               tion is driven primarily through the development of a strong,
must exceed a critical value         3 for unstable convective                 negative entropy gradient within the gain region, rather than
plumes to grow before being advected downward through the                      through the non-linear development of SASI modes.
gain region below the shock. Here ωBV is the Brunt-Väisälä                        In Paper I we considered the non-linear, saturated state of
frequency.                                                                     the SASI in the absence of neutrino heating, and showed that
   Using our initial flow models, we can translate H into a                     the amplitude of the shock oscillations drops significantly as
value for χ, and find (Table 1) that typically χ ∼ 5 − 20 at                    the dissociation energy ε is increased. We now explore how
  9 Note that r is the shock radius in the time-independent flow solution.      the r.m.s. amplitude of the shock oscillations correlates with
For a fixed cooling function, rs is a monotonically increasing function of H,   the strength of heating. To eliminate the effect of secular
and equals rs0 at H = 0.                                                       shock motions around or above the threshold for explosion,
                                                        ALPHA PARTICLE RECOMBINATION IN SUPERNOVAE                                                                                      17
∆ (a0 - <a0>) / rs0
                                   2                                                                                                  ε = 0.15vff0              χ = 1.5       χ = 3.9
                       0.1   ε / vff0 = 0.15

                        0                                                                                                  0.1

                                                                                                           a1 / rs0

                                                                                     ∆ (a1-<a1>)/rs0
                      -0.2                                                     0.5

                                                                               0                                           -0.1
∆ (a2-<a2>)/rs0

                      0.5                                                      -0.5                                        0.5

                                                                                                       rs0〈 ds/dr〉 θ,med
                        0                                                                                                    0

                      -0.5                                                                                                 -0.5

                                                                                     ∆ (a0-<a0>)/rs0
                             rs0 = 75 km (NSE)                                 0.2

∆ (a1-<a1>)/rs0

                                                                                                              a2 / rs0

                        0                                                                                                    0
                                                                                     ∆ (a2-<a2>)/rs0
                      -0.5                                                     0.5
                                                                                                       rs0〈 ds/dr〉 θ,med
                               0         0.002     0.004       0.006   0.008
                                                           3                                                               -0.5
                                               H / (rs0 vff0 )
    F IG . 18.— Amplitude (r.m.s.) of = 0, 1, 2 modes of the shock in two model                                             -1
                                                                                                                                  0        20            40              60   80        100
   sequences with varying heating parameter H. Stars indicate exploding runs.
   We show the r.m.s. fluctuation of the difference between the instantaneous                                                                                  t / tff0
   Legendre coefficient a and a running average a 50t that is computed over
   a window of width 50tff0 (see text). This subtracts the secular movement of                            F IG . 19.— Principal shock Legendre coefficient and the radial median of
   the shock in runs that are close to or above the threshold for explosion. Note                       the angle-averaged entropy gradient for the ε = 0.15v2 0 model and two dif-
   that the amplitude is measured in absolute units (rs0 ).                                             ferent heating rates, corresponding to χ = 1.5 (H = 0.002 rs0 v3 0 , blue curves)
                                                                                                        and χ = 3.9 (H = 0.004 rs0 v3 0 , red curves). Both runs are below the threshold
                                                                                                        for an explosion, but vigorous convection is established throughout the gain
                                                                                                        region in the run with the higher heating rate. Top (bottom) two panels: seed
                                                                                                        perturbation is a shell with = 1 ( = 2) density profile. A running average of
   we first calculate the running average a 50t of the shock Leg-                                        the entropy gradient (temporal width 20tff0 ) appears as thick solid lines.
   endre coefficients a over an interval 50tff0 , and then calculate
   the r.m.s. fluctuation of a ≡ a − a 50t over the duration of
   each simulation. The result is plotted in Fig. 18 as a function
   of H for two model sequences (ε = 0.15v2 0 and NSE with
                                                 ff                                                     20 display the Legendre coefficients of the shock alongside
   rs0 = 75 km).                                                                                        the radial median (over the gain region) of the angle-averaged
      There is a clear trend of increasing = 1, 2 mode ampli-                                           entropy gradient. We find that the amplitude of the = 1 and 2
   tude with increasing heating rate. This confirms our previ-                                           modes is strongly tied to the strength of convection. For both
   ous suggestion that large-amplitude shock oscillations require                                       types of dissociation models, convection is quenched by the
   strong heating when the dissociation energy behind the shock                                         accretion flow when χ < 3: it grows intermittently in strength,
   exceeds ∼ 0.15v2 0 . The models with H = 0 reveal a slight ex-
                     ff                                                                                 but never reaches large enough amplitudes to significantly dis-
   ception to the overall trend: the r.m.s. amplitude of the shock                                      tort the shock surface. As a consequence, the entropy gradient
   oscillations appears larger than it does in models with small                                        remains shallow and negative most of the time, with coherent
   but finite heating rate, because the oscillations are strongly in-                                    shock oscillations taking place simultaneously. The shock os-
   termittent at H = 0 (see Paper I). The shock oscillations grow                                       cillations have a low amplitude due to the large dissociation
   much stronger just below the threshold for explosion (explod-                                        energy.
   ing runs are marked by stars), above which they seem to satu-                                           Convection grows much more rapidly when χ > 3, dis-
   rate. Note also that their amplitudes do not vary much with the                                      torting the shock surface before the SASI has the chance to
   choice of dissociation model. The r.m.s. amplitudes relative                                         execute a few oscillations. In this case, the entropy gradi-
   to the running average of a0 at the threshold for explosion are                                      ent is initially more negative, but quickly flattens. Indeed,
   {5%, 12%, 8%} for the = 1, 2, 3 modes in the ε/v2 0 = 0.15
                                                          ff                                            the = 1, 2 amplitudes only become large in models where
   sequence, and {6%, 12%, 7%} in the NSE rs0 = 75 km se-                                               neutrino-driven convection is strong enough to flatten the en-
   quence.                                                                                              tropy gradient.
      We have performed an additional sequence of runs in which                                            Another way of seeing that convection is the forcing agent
   we drop an overdense shell with a given Legendre index                                               behind shock oscillations when χ > 3 is to analyze the dis-
   through the shock (see also Paper I). This has the effect of                                         tribution of vorticity in the gain region. Figure 21 shows a
   selectively triggering individual SASI modes. Figures 19 and                                         histogram of vorticity vs. Mach number at three different in-
    18                                                               FERNÁNDEZ & THOMPSON
                               rs0 = 75 km          χ = 2.1         χ = 5.5
                    0.1                                                                                           Normalized Mass Histogram
    a1 / rs0

                      0                                                                                          0.01 0.10 0.30 0.50 0.70 0.90
                    -0.1                                                                                    1
                                                                                                                                           t = 20 tff0
                    0.5                                                                                    0.8
rs0〈 ds/dr〉 θ,med


                    0.1                                                                                    1.0
                                                                                                                                           t = 35 tff0
       a2 / rs0

                                                                                             Mach number
rs0〈 ds/dr〉 θ,med


                                                                                                                                           t = 60 tff0
                           0         20        40              60   80        100
                                                    t / tff0                                               0.6
   F IG . 20.— Same as Fig. 19, but for NSE models with α-particles and                                    0.4
 rs0 = 75 km.
 stants in the evolution of the run with ε = 0.15v2 0 , = 1 seed
                                                                                                            -10       -5          0             5        10
 perturbation, and χ = 3.9 (upper two panels in Figure 19). At                                                             φ-Vorticity [tff0-1]
 t = 20tff0 , convection is just getting started and the vortical
 motions are restricted to Mach numbers 0.3. However, by
                                                                                      F IG . 21.— Histogram of vorticity vs. Mach number in the gain region,
 t = 35tff0 the Mach number distribution extends up to M 0.5                        weighted by mass, at three different instants in the evolution of the run with
 and has almost reached its asymptotic form (t = 60tff0 ), at                       ε = 0.15v2 0 , = 1 seed perturbation, and χ = 3.9 (corresponding to the upper
 the same time that convection has filled the region below the                       two panels in Fig. 19). The distribution broadens once convection is fully
 shock. The dipolar mode of the shock develops a large ampli-                       developed, but before the dipolar shock mode shows significant growth. The
                                                                                    vertical lines show the vorticity of a convective flow with period equal to
 tude only after this fully developed convective state has been                     (solid) the mean radial advection time and (dashed) the period of a lateral
 reached. The convective rolls are a source of acoustic radia-                      sound wave at the midpoint between r∗ and rs . An animated version of this
 tion (e.g. Goldreich & Kumar 1988), which will drive a dipo-                       figure is available in the online version of the article.
 lar oscillation of the shock if the overturn frequency is compa-
 rable to the frequency of the = 1 mode, | ×v| ∼ 2×2π/tadv .
 This zone is marked by the vertical solid lines in Fig. 21, and
 indeed encompasses most of the mass.                                               turn, on the core structure of the progenitor star and the den-
                                                                                    sity profile in the forming neutron star. Within the framework
                                             7. SUMMARY                             explored in this paper, we find two extreme types of explo-
   We have investigated the effects of α-particle recombi-                          sion. In the first, neutrino heating does most of the work,
 nation and neutrino heating on the hydrodynamics of core-                          with a significant final boost from α-particle recombination.
 collapse supernovae, when the heating rate is pushed high                          In the second, neutrino heating is generally less important at
 enough to reach the threshold for an explosion. The effect                         promoting material below the shock to positive energies. By
 of dimensionality has been probed by comparing 1D (spher-                          way of comparison, the shock stalls at 100 − 140 km and no
 ical) and 2D (axisymmetric) time-dependent hydrodynamic                            explosion results in the 15 M collapse calculation of Buras
 calculations. Our main results can be summarized as follows:                       et al. (2006b); whereas the shock stalls at ∼ 200 km and ex-
                                                                                    plosion is (just barely) found in the 11.2 M model of Buras
 1. – The critical heating parameter that yields an explosion                       et al. (2006a).
 depends sensitively on the starting position of the shock rel-                     2. – During the final stages of an explosion, the heat released
 ative to rα . This means that the critical neutrino luminosity                     by α-particle recombination is comparable to the work done
 depends sensitively on the stall radius of the shock and, in                       by adiabatic expansion. This heat is concentrated in material
                                  ALPHA PARTICLE RECOMBINATION IN SUPERNOVAE                                                         19

that has previously been heated by neutrinos. Significantly            ence between 1D and 2D explosions: the mechanism is fun-
more energy is lost through α-particle dissociation in fresh          damentally non-linear in two dimensions.
downflows, so that nuclear dissociation remains on balance             7. – Vortical motions with a Mach number ∼ 0.3-0.5 first ap-
an energy sink within the accretion flow.                              pear at the onset of convective instability around the radius
3. – The large-amplitude oscillations that are seen in 1D runs        of maximal neutrino heating, but before the dipolar mode of
near an explosion are the consequence of the = 0 SASI as              the shock reaches its limiting amplitude. These vortices are
modified by heating. In contrast with the = 1, 2 modes of a            a source of acoustic waves, which have a similar period to
laminar accretion flow, the period of these oscillations is close      the large-scale oscillation of the shock. Near the threshold
to twice the post-shock advection time. The critical heating          for an explosion, the turbulence in the gain region becomes
rate for an explosion (assuming constant mass accretion rate,         supersonic, as the existence of widespread secondary shocks
neutrino luminosity, and inner boundary) corresponds to neu-          attests. These shocks convert turbulent kinetic energy to in-
tral stability for the = 0 mode.                                      ternal energy, increasing the effective heating rate.
4. – The critical heating parameter Hcr for an explosion is
generally lower in 2D than in 1D, but the difference becomes             There are at least two reasons why explosions by the mech-
smaller as the starting radius of the shock approaches rα .           anism investigated here may be more difficult in 3D than in
5. – Non-spherical deformations of the shock are tied to the          2D. First, the existence of more degrees of freedom for the
formation of large-scale plumes of material with positive en-         low-order modes of the shock in 3D implies that the ampli-
ergy. Our two-dimensional explosions with a super-critical            tude of individual shock oscillations is lower. As a result, it is
heating rate involve a large-scale convective instability that        more difficult for the shock to extend out to the radius where
relies on the accumulation of vorticity on the largest spa-           α-particle recombination gives it the final push. Second, an
tial scales. Volume-filling convective cells are apparent in a         explosion that is driven by neutrino heating in 2D involves
time-averaged sense. Transient heating events create positive-        the accumulation of vorticity on the largest spatial scales, an
energy material that accumulates in between the convective            effect that is special to two dimensions. A full resolution of
cells and the shock. A significant fraction of the heating oc-         these issues is possible only with high-resolution 3D simula-
curs in horizontal flows at the base of the convective cells,          tions.
which are fed by a dominant equatorial accretion plume. If
the heating parameter is large enough, this results in an am-
plifying cycle and explosion.
6. – The amplitude of the = 1 and 2 modes correlates                     We are grateful to Jonathan Dursi for scientific discussions
strongly with the value of the heating parameter, and is cou-         and help with FLASH. The software used in this work was
pled to the appearance of vigorous neutrino-driven convection         in part developed by the DOE-supported ASC / Alliance Cen-
below the shock. In agreement with the work of Foglizzo et al.        ter for Astrophysical Thermonuclear Flashes at the University
(2006) and Scheck et al. (2008), we find that χ ≈ 3 marks the          of Chicago. Computations were performed at the CITA Sun-
transition from a strong linear instability in a nearly laminar       nyvale cluster, which was funded by the Canada Foundation
flow below the shock, to a volume-filling convective instabil-          for Innovation. This research was supported by NSERC of
ity. In all of our simulations, the threshold for explosion lies      Canada. R. F. is supported in part by the Ontario Graduate
well within the latter regime. This highlights a basic differ-        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 Ye = 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 §2.1.1. The
temperature does not appear explicitly in the FLASH hydro 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 kB T > me c2 /2, it can be written (Bethe et al. 1980):
                                         1 (kB T )4 11π 2           1                3   ρ
                                     p=                   + 2η 2 + 2 η 4 + 1 − Xα           kB T,                         (A1)
                                        12 (¯c)3
                                              h       15           π                 4   mu
where η = µe /(kB T ) the normalized electron chemical potential, also known as degeneracy parameter, and ¯, c, and mu are
Planck’s constant, the speed of light, and the atomic mass unit, respectively. The density and degeneracy parameter are further
related by
                                                        mu     kB T
                                                   ρ= 2                 η(π 2 + η 2 ),                                    (A2)
                                                      3π Ye ¯c  h
where Ye is the electron fraction. The equilibrium fraction of α-particles is given by the nuclear Saha equation,
                                                          3                                          3/2
                               2 2   1     mu nQ (T )          Qα                          mu kB T
                             Xn Xp = Xα                exp −         ;         nQ (T ) =                                          (A3)
                                     2         ρ              kB T                         2π¯2h
as supplemented by the conditions of mass and charge conservation,
                                                      Xn + Xp + Xα = 1                                                            (A4)
                                                         Xp + Xα = Ye .                                                           (A5)
20                                                                                            FERNÁNDEZ & THOMPSON
                                                                                                                                     eq                                      eq
                                       S = 5 kB/nucleon                                   S = 25 kB/nucleon                     Xα = 0.5                                 ∂ Xα /∂ ln ρ
                                                                                                   eq                                eq                                      eq
                                       S = 15 kB/nucleon                                  Xα = 0.9                              Xα = 0.1                                 ∂ Xα /∂ ln p

                     30                                                                       10                                                                   100
 log10 p [erg cm ]

                     29       (a)                                                                                                                                                                          (e)

                                                                                                                                             Prel / Pmat
                     28                                                                                                                                            10

                     26                                                                                                           (c)                                1

                                                                 ∂ Xα /∂ ln p, ∂ Xα /∂ ln ρ

                                                                                                                                            γ = (∂ ln p/∂ ln ρ)S
                                                                                                          (d)                                                      1.6            S = 10 kB/nucleon        (f)
                              (b)                                                              1
   kBT [MeV]


                      2                                                                        0

                      1                                                                       -1                                                                   1.2

                      0                                                                       -2                                                                     1
                          6     7      8         9        10                                       6         7      8       9         10                                 6        7       8       9          10
                                                     -3                                                                         -3                                                                    -3
                                log10 ρ [g cm ]                                                             log10 ρ [g cm ]                                                       log10 ρ [g cm ]
 F IG . A1.— Equation of state of a fluid containing n, p, α, photons, and finite-temperature and partially degenerate electrons, in nuclear statistical equilibrium.
We solve eqns. (A1)-(A5) and tabulate all quantities on a grid of density and pressure. Panels (a), (b) and (c): pressure, temperature, and degeneracy parameter as
a function of density, for fixed values of Xα and entropy. Dotted line: η = π, the approximate boundary between degenerate and non-degenerate electrons. The
                                                                      eq                                         eq
gray area shows the region of the thermodynamic plane where XO = 0.5. Panel (d): partial derivatives of Xα with respect to density (solid lines, positive) and
pressure (dashed lines, negative). Panel (e): ratio of relativistic pressure (photons and pairs) to material pressure (α-particles and nucleons). Panel (f): adiabatic
index γ for different adiabats (dotted line: γ = 4/3).

In eqs. (A3)-(A5), Xn and Xp are the mass fractions of free neutrons and protons, respectively, and Qα = 28.3 MeV is the binding
energy of an α-particle. Combining eqs. (A1) and (A2) gives η and T in terms of p and ρ. The equilibrium mass fraction Xα
                                                                           eq     eq               eq
is calculated from ρ and T . For numerical calculations, we tabulate Xα , ∂Xα /∂ ln ρ, and ∂Xα /∂ ln p for a grid of density and
pressure. In addition, we tabulate partial derivatives of T to substitute in eqs. (B10) and (B11).
   Figure A1 shows contours of constant Xα 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),
                                           (11π 2 /15 + η 2 )      3                                    5      mu nQ (T )                              1
                               S = π 2Ye                      + 1 − Xα                                    + ln                  − Xp ln Xp − Xn ln Xn − Xα ln Xα /32 .                                       (A6)
                                             η(π 2 + η 2 )         4                                    2          ρ                                   4
The postshock density in the initial configuration is typically ρ2 ∼ 109 g cm−3 with an entropy ∼ 10 − 15kB /nucleon. The for-
mation of α-particles that is seen in Fig. 1 results from an expansion of the shock into the part of the thermodynamic plane in
Figs. A1a,b where ρ2 < 109 g cm−3 and T 1 MeV. In this regime, the electrons are non-degenerate and the pressure in pho-
tons and pairs begins to exceed the nucleon pressure. The dip in the adiabatic index seen in Fig. A1f results from α-particle
dissociation/recombination, which partially compensates the change in internal energy due to compression/expansion.

  We write down the ordinary differential equations that are used to compute the initial flow, and the density profiles in Fig. 2a,b,c.
The steady state Euler equations in spherical symmetry are
                                                1 dvr 1 dρ 2
                                                       +      + =0                                                              (B1)
                                                vr dr ρ dr r
                                                    dvr 1 dp
                                                 vr     +     +g=0                                                              (B2)
                                                    dr ρ dr
                                                    deint pvr dρ
                                               ρvr       −       = LH − LC + Lα ,                                               (B3)
                                                     dr    ρ dr
where eint is the internal energy per unit mass, LH , LC , and Lα the source terms described in eqns. (9)-(11), and g = GM/r2 .
Since two variables suffice to describe the thermodynamic state of a system, we write
                                                      deint      dp     dρ
                                                            ≡ Ep    + Eρ .                                                (B4)
                                                       dr        dr     dr
                                 ALPHA PARTICLE RECOMBINATION IN SUPERNOVAE                                                       21

                                                               dp      dρ
                                                                   + Aρ .
                                                            Lα ≡ A p                                                  (B5)
                                                               dr      dr
The coefficients Ei and Ai encode the dependence on the equation of state. Replacing eqns. (B4) and (B5) in (B3), and using
eqns. (B1) and (B2) to eliminate the pressure derivative, we obtain
                                             dρ   ρvr E p − A p ρg − 2ρv2 /r + (LH − LC )
                                                =                                          .                                    (B6)
                                             dr    ρvr Eρ − pvr /ρ − Aρ + v2 ρvr E p − A p

  The coefficients in eqs. (B4) and (B5) work out to
                                                   Ep =                              (constant γ)                               (B7)
                                                        (γ − 1)ρ
                                                   Eρ = −                            (constant γ).                              (B8)
                                                          (γ − 1)ρ2
for a constant-γ equation of state, eint = p/[(γ − 1)ρ]. The pressure (eq. [A1]) in the NSE model described in Appendix A can be
decomposed into contributions from relativistic particles and from nucleons, p = prel + pmat , and the specific internal energy is
                                                   1        3               p 3    3                 kB T
                                          eint =     3prel + pmat         =3 −  1 − Xα                    .                     (B9)
                                                   ρ        2               ρ 2    4                  mu
One therefore finds
                                    3 9 kB T ∂Xα 3       3     1 ∂(kB T )
                               Ep =  +           −   1 − Xα               ,                                   (NSE)           (B10)
                                    ρ 8 mu ∂ p 2         4    mu ∂ p
                                     3p 9 kB T ∂Xα 3       3     1 ∂(kB T )
                               Eρ = − 2 +          −    1 − Xα              .                                 (NSE)           (B11)
                                     ρ    8 mu ∂ρ 2        4     mu ∂ρ
   The initial postshock solution is obtained by integrating the above equations from rs to an inner radius r∗ at which the flow
stagnates. We iterate the normalization of the cooling function in eq. (9) so that r∗ = 0.4rs0 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
rs > rs0 (Fig. 2a,b). The initial Mach number at the inner boundary is chosen so as to satisfy
                                                                                          GM             ˙
                                            LH, i − LC, i + Lα, i Vi           0.995          − ε(t = 0) M ,                  (B12)
where the sum is taken over the computational cells below the shock at our fixed resolution (see §2.1.3), Vi 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−3 − 10−2 .
  When including α-particles in the EOS, one needs to calculate self-consistently the value of Xα below the shock, the corre-
sponding dissociation energy ε(t = 0) [eq. (8)], and compression factor κ [eq. (3)]. The density upstream of the shock is obtained
                                                       ρ1 (rs ) =    2 |v (r )|
                                                                                ,                                           (B13)
                                                                  4πrs 1 s
                                                                              vff (rs )
                                                     v1 (rs ) = −                                                             (B14)
                                                                    1 + 2M−2 (rs )/(γ − 1)

is the upstream velocity at r = rs , while the upstream pressure satisfies
                                                                        ρ1 (rs )[v1 (rs )]2
                                                           p1 (rs ) =                       .                                 (B15)
                                                                           γM2 (rs )
Eqns. (B13) and (B15) are transformed to physical units for input to the NSE model by adopting M1 (rs0 ) = 5, M = 0.3 M s−1 ,
M = 1.3 M , and a particular value for the shock radius rs0 in the absence of heating.

                                                          C. TIME EVOLUTION
  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 hydro 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.
22                                                             FERNÁNDEZ & THOMPSON

  In the NSE nuclear dissociation module, numerical stability is maintained using an implicit update of the pressure in between
hydro sweeps,
                                             pnew = pcur + (γ − 1)ρenuc (ρ, pnew ),                                        (C1)
where enuc 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 |enuc | <
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 × 1010 g cm−3 .
   Both the constant-ε and the NSE dissociation modules require that the Mach number stays below a fixed value Mburn 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 Mburn = 2 in most of our simulations. The critical heating parameter Hcr depends
weakly on Mburn : changes in Mburn 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, Mburn is just below the
Mach number of the upstream flow. We have tried expanded the outer boundary to r = 9rs0 (with a somewhat smaller Mburn ) and
found that runs that did hit r = 7rs0 still hit the new outer boundary.
   Aside from the inclusion of heating and the NSE dissociation module, 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.)

Audi, G., Wapstra, A. H., & Thibault, C. 2003, Nuclear Physics A, 729, 337       Goldreich, P., & Kumar, P. 1988, ApJ, 326, 462
Bethe, H. A. 1990, Reviews of Modern Physics, 62, 801                            Herant, M., Benz, W., & Colgate, S. 1992, ApJ, 395, 642
—. 1993, ApJ, 412, 192                                                           Herant, M., Benz, W., Hix, W. R., Fryer, C. L., & Colgate, S. A. 1994, ApJ,
—. 1996, ApJ, 469, 737                                                             435, 339
Bethe, H. A., Applegate, J. H., & Brown, G. E. 1980, ApJ, 241, 343               Houck, J. C., & Chevalier, R. A. 1992, ApJ, 395, 592
Bethe, H. A., & Wilson, J. R. 1985, ApJ, 295, 14                                 Iwakami, W., Kotake, K., Ohnishi, N., Yamada, S., & Sawada, K. 2008, ApJ,
Blondin, J. M., & Mezzacappa, A. 2006, ApJ, 642, 401                               678, 1207
Blondin, J. M., Mezzacappa, A., & DeMarino, C. 2003, ApJ, 584, 971               Janka, H.-T., & Mueller, E. 1996, A&A, 306, 167
Buras, R., Janka, H.-T., Rampp, M., & Kifonidis, K. 2006a, A&A, 457, 281         Kitaura, F. S., Janka, H.-T., & Hillebrandt, W. 2006, A&A, 450, 345
Buras, R., Rampp, M., Janka, H.-T., & Kifonidis, K. 2006b, A&A, 447, 1049        Liebendörfer, M., Mezzacappa, A., Thielemann, F.-K., Messer, O. E., Hix,
Burrows, A., & Goshy, J. 1993, ApJ, 416, L75                                       W. R., & Bruenn, S. W. 2001, Phys. Rev. D, 63, 103004
Burrows, A., Hayes, J., & Fryxell, B. A. 1995, ApJ, 450, 830                     Marek, A., & Janka, H. . 2007, ApJ submitted, astro-ph/0708.3372
Burrows, A., Livne, E., Dessart, L., Ott, C. D., & Murphy, J. 2006, ApJ, 640,    Murphy, J. W., & Burrows, A. 2008, ApJ, 688, 1159
  878                                                                            Ohnishi, N., Kotake, K., & Yamada, S. 2006, ApJ, 641, 1018
—. 2007, ApJ, 655, 416                                                           Quirk, J. J. 1994, Int. Jour. Num. Meth. Fluids, 18, 555
Calder, A. C., Fryxell, B., Plewa, T., Rosner, R., Dursi, L. J., Weirs, V. G.,   Scheck, L., Janka, H.-T., Foglizzo, T., & Kifonidis, K. 2008, A&A, 477, 931
  Dupont, T., Robey, H. F., Kane, J. O., Remington, B. A., Drake, R. P.,         Scheck, L., Kifonidis, K., Janka, H.-T., & Müller, E. 2006, A&A, 457, 963
  Dimonte, G., Zingale, M., Timmes, F. X., Olson, K., Ricker, P., MacNeice,      Thompson, C. 2000, ApJ, 534, 915
  P., & Tufo, H. M. 2002, ApJS, 143, 201                                         Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Reviews of Modern
Fernández, R., & Thompson, C. 2008, ApJ submitted, astro-ph/0811.2795              Physics, 74, 1015
Foglizzo, T., Galletti, P., Scheck, L., & Janka, H.-T. 2007, ApJ, 654, 1006
Foglizzo, T., Scheck, L., & Janka, H.-T. 2006, ApJ, 652, 1436
Fryxell, B., Olson, K., Ricker, P., Timmes, F. X., Zingale, M., Lamb, D. Q.,
  MacNeice, P., Rosner, R., Truran, J. W., & Tufo, H. 2000, ApJS, 131, 273

Shared By: