Document Sample

S UBMITTED TO A P J Preprint typeset using LTEX style emulateapj v. 6/22/04 A DYNAMICS OF A SPHERICAL ACCRETION SHOCK WITH NEUTRINO HEATING AND ALPHA-PARTICLE RECOMBINATION RODRIGO F ERNÁNDEZ 1 AND C HRISTOPHER T HOMPSON 2 Submitted to ApJ ABSTRACT 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 signiﬁcant 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 ﬂow not have a clear picture of a robust path to an explosion in stars composed of a zero-energy, polytropic ﬂuid (Blondin & Mez- that form iron cores. Heating by the absorption of electron- zacappa 2006) via a linear feedback between ingoing vor- type neutrinos signiﬁcantly modiﬁes the settling ﬂow 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 sufﬁciently 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 signiﬁcant sink of thermal energy in the accretion ﬂow 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 ﬂows 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 ﬂuid then becomes substantially time of settling material in the zone of strong neutrino heat- negative. A signiﬁcant 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 signiﬁcantly 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 ﬁnite-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 ﬂuid (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 ﬁnite-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 proﬁle 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 ﬁelds. We ﬁnd 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 ﬂow 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 ﬁnd 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 ﬂuids 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 ﬁxed 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 ﬁrst introduced neutrino heating over cooling). We ﬁnd that most of the heat by Burrows & Goshy (1993), and has recently been used by deposition by neutrinos is concentrated in lateral ﬂows 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 ﬂow 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 ﬂow as rial to positive Bernoulli parameter, but only if the shock starts a polytropic ﬂuid, from which a ﬁxed 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 ﬁnite- 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 signiﬁcant simpliﬁcation 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 ﬁndings 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 conﬁguration is a steady, spheri- pecially in the dynamics of this high-entropy material, we set cally symmetric ﬂow onto a gravitating point mass M. The Ye = 0.5 when evaluating the α-particle abundance. To obtain ﬂow contains a standing shock wave, and the settling ﬂow be- a realistic density proﬁle, we continue to approximate the in- low the shock cools radiatively in a narrow layer outside the ternal energy of the ﬂuid as that of a polytropic ﬂuid with a inner boundary of the simulation volume. The space of such 4 ﬁxed 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 ﬂow, 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 simpliﬁcations, our results already show it covers a narrower range than the other three. The infalling many similarities with more elaborate collapse calculations. material is signiﬁcantly de-leptonized before hitting the shock One-dimensional explosions are due to a global instability re- only in the ﬁrst 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 modiﬁed by heating. As in Paper I, we We explore a two-dimensional surface through this three- ﬁnd that the period of the = 0 mode remains close to twice dimensional parameter space by i) ﬁxing 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 modiﬁed 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 ﬂow outside the shock (Houck & Chevalier 1992). The κ≡ = (γ + 1) γ + M−2 − 1 ρ1 secular cooling of the collapsed core forces a gradual decrease −1 ˙ in r∗ , and M also varies with time and with progenitor model. 2ε 2 Given the important role that α-particle recombination 1 − M−2 + (γ 2 − 1) 2 1 , (3) v1 plays in the ﬁnal 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 speciﬁc nuclear dissociation GMmα energy ε is deﬁned 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 ﬂow 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 ﬁxed speciﬁc 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 ﬁxed speciﬁc energy ε right below the shock, as done out the settling ﬂow. 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 ﬂow. The heating rate remains freely adjustable thereafter. The main limitation of this approximation is that the dissoci- We adopt this simpliﬁcation because we do not intend to ation energy does not change with the radius (or inclination) ﬁnd 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 . ff 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-reﬁnement 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 ﬂowing 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 ﬂow conﬁguration. 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 ﬂow 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 ﬂow. 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 ﬂuid velocity, which is radial in the initial responds to the speciﬁc dissociation energy condition, p is the pressure, ρ is the mass density, and G is QO r −1 Newton’s constant. The ﬂow 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 ﬂuid 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 ﬂow, 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 ﬂuid 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 1 determines the density proﬁle inside a radius ∼ 2 rα , where 3 Although heavier nuclei can begin to recombine once the shock moves the α-particle abundance is very low. signiﬁcantly beyond rα , this generally occurs only after the threshold for an The upstream and downstream ﬂow proﬁles 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 1 We ﬁnd that α-particles appear in signiﬁcant 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 signiﬁcantly 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 ] 0.6 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 eq 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, ﬁnite-temperature and 2 ε(t=0) / vff0 [ 56Fe ] partially degenerate equation of state for electrons and nucle- 0.2 ons; see Appendix A for details. Speciﬁc choices must then ˙ be made for the parameters rs0 , M, and M; we generally take 2 ε(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. eq A speciﬁc energy F IG . 1.— Equilibrium mass fraction of α-particles Xα , and ratio of initial eq dissociation energy ε(t = 0) [eq. 8] to v2 behind a spherical shock positioned ff0 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 ﬂuid ele- brackets refer to the upstream composition of the accretion ﬂow, 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 ﬂow). 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 eq Figure 1 shows how ε(t = 0) and Xα depend on the shock ra- accumulates in the ﬁrst few computational cells adjacent to dius rs0 , for upstream ﬂows 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- ﬂow. proximately constant inside ∼ 75 km, where the downstream We model neutrino heating as a local energy generation rate ﬂow 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 ﬂux 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 ﬂow 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 3/2 α-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) 4 (9) when Ye = constant and γ = 3 (corresponding to ρ ∝ r−3 dr ∂ρ dr ∂ p dr 4 In the case where the upstream ﬂow is pure 56 Fe, we replace ε in eq. (8) This energy generation rate is negative, as the temperature O 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 ff 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 reﬁnement 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 reﬂecting 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 ﬁxed at r = 7rs0 , and is set by the H = 0.007 ρ [g cm ] -3 1e+10 H = 0.009 upstream ﬂow 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 ﬂuid 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 ﬂow. This scalar is passively full + α ρ [g cm ] -3 advected by FLASH2.5. Through this technique, we are able 1e+10 γ = 1.48 + α to assign a “ﬂuid" 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 1e+08 rOB dr 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 ﬂuid enters through the outer radial boundary at r = rOB , and rs (t) θ is the angle av- F IG . 2.— Sample initial density proﬁles, which are solutions to the spher- eraged shock position. Initially, tOB = 0 and all the ﬂuid be- ically symmetric and time-independent ﬂow equations. Panel (a) shows the zero-heating conﬁgurations 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 conﬁgurations, and M = 0.3M s−1 , for statistical studies (§4.2), tracing large scale ﬂuid patches, M = 1.3M for the NSE models. Panel (b) shows a sequence with a ﬁxed despite some inevitable turbulent mixing among neighboring cooling function and range of heating rates (H is given in units of rs0 v3 0 ). ff ﬂuid parcels. The dashed line shows the upstream ﬂow. Panel (c) shows a sequence with An explosion is deﬁned 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 ﬁnite-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 ﬂow 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 deﬁnition 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 ﬂow 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 1 (1D) and axial (2D) symmetry. Their parameters are summa- (twice) the kinetic energy of the upstream ﬂow, 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 inﬂuences 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 proﬁles (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 0 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 eq 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 2 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). 0.05 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 ) 3 while the one with H = 0.009v3 rs0 is well above thresh- ff0 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 signiﬁcant in a layer below the shock, parameter H around the threshold for explosion. Stars denote conﬁgurations causing the density proﬁle 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 proﬁle ﬂattens, akin to an increase in γ. Dotted to the internal energy of the ﬂow 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 ﬁnite-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 ﬁle is well approximated by an ideal gas with γ 1.48 at zero Boltzmann neutrino transport. heating. Adding in heating tends to ﬂatten the density proﬁle The origin of the 1D SASI oscillation can be brieﬂy summa- even more, and with γ = 1.48 it would be much ﬂatter 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 ﬂow, we ﬁnd that this relation is from the shock to the cooling layer, maintained. The = 0 mode is damped until the heating rate rs 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 deﬁned 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 ﬂow equations. Large amplitude 1D When heating is added, the density proﬁle ﬂattens. 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 2 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 −1 rα Μ1.3 2 0 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 coefﬁcients H near the critical ff value Hcr = (0.006625 ± 0.000125)v3 0 rs0 . The lower panel shows runs with ff eq 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 . ff 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 proﬁles 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 speciﬁc nuclear binding We ﬁnd 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 deﬁne 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 . s 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 conﬁgurations with con- stant dissociation. These values were obtained by solving the differential system of Foglizzo et al. (2007), modiﬁed 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 ﬁnite, 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 ﬂow 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 ﬂuid downstream of the shock, and therefore is sensitive 5 We deﬁne our critical heating parameter H to be the average of the val- to the density proﬁle below the shock. Heating by neutrinos cr 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 ﬁducial 1000tff0 cutoff. distance ∼ (2 − 3)r∗ . Heating by α-particle recombination is 8 FERNÁNDEZ & THOMPSON 200 2 6 ε = 0.15 vff0 2 1.06 Hcr ε = 0.15 vff0 1.02 Hcr 150 0.98 Hcr rs,max / |drs,max/ dt| 0.91 Hcr 〈 tres 〉 vol (upper 10%) 4 〈 tres 〉 vol (upper 50%) 100 2 50 rs0 = 75 km 1.06 Hcr 6 1.02 Hcr 0.98 Hcr 150 rs0 = 75 km 0.90 Hcr 4 100 -1 rα M1.3 2 50 time [tff0] rs / rs0 0 rs0 = 50 km 6 -1 rα M1.3 rs0 = 50 km 150 1.08 Hcr 4 1.02 Hcr 0.98 Hcr 100 0.95 Hcr 2 50 0 6 rs0 = 125 km 1.19 Hcr 400 1.04 Hcr 0.96 Hcr rs0 = 125 km 4 0.89 Hcr 300 -1 200 rα M1.3 2 100 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 ff rates Hcr are different for each conﬁguration, 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 ﬁt 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 deﬁnition 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 ﬂuid. 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 proﬁles of the shocked ﬂow 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 speciﬁc 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 ﬂuid. We ﬁrst 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- ﬁrst 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 ﬂow. We can summarize this behavior as follows: during the ini- Large deformations of the outer shock are generally caused tial expansion phase, ﬂuid 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 ﬂow (as expected from eq. [11]). Some ﬂuid 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 stratiﬁed 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 inﬂuence on the details of breakout. When the dius expands by a factor ∼ 2.5 between the upper two frames dissociation energy ε is held ﬁxed, κ 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 ﬂuid terial is not accompanied by a signiﬁcant expansion of the in- below the shock when rs is such that ε ∼ v2 /2. In spheri- ff ner bound region, whose outer radius remains ﬁxed 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 inﬂuence 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 ﬂow 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 r dissociation energy and with NSE between n, p, and α below where rs is the shock radius, vr is the ambient radial ﬂow 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 ﬂow once it grows to a size ∆Ωbubble ∼ M2 Sr, where Mcon con 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 inﬂuenced 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 signiﬁcant 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 ﬁnd 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 ﬁgure displays Bernoulli parameter b; the middle ﬁgure the rate of change of nuclear binding energy per unit mass; and the bottom ﬁgure 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 conﬁgurations 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 signiﬁcantly 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 ﬂuid: 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 ﬂow. Strong recombi- online version of the article. nation heating quite naturally extends below the zone where α-particles are present in signiﬁcant numbers, as is seen in Fig. 9. The ﬁrst and third panels of this ﬁgure 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 downﬂows. 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 ﬂow (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 ﬂuid 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- ﬂows is represented by the long tail toward large negative val- ues of denuc /dt. The α-particle concentration is very stratiﬁed, with higher Xα occurring at larger radius. 4.2. Residency Time F IG . 9.— Top panels: rate of release of speciﬁc 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 ﬁrst assign a unique “ﬂuid time" tF (eq. [12]) to each infalling radial mass shell in the simulation, 12 FERNÁNDEZ & THOMPSON 0.5 rs0 = 50 km 0 -0.5 t = 250 tff t = 350 tff t = 415 tff t = 455 tff -1 dE / db [ρ1 vff,0 rs0 ] 2 0.5 rs0 = 75 km 0 -0.5 . t = 200 tff t = 290 tff t = 370 tff t = 434 tff -1 0.5 rs0 = 125 km 0 -0.5 t = 100 tff t = 200 tff t = 270 tff t = 333 tff -1 -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 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 (deﬁned 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 downﬂows. 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 deﬁne6 the residency time of the ﬂuid as tribution). The ﬁnal 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 ﬁnd 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- ﬂuid from the accretion ﬂow. 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 deﬁned 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 ff tion funnel and convective cells, this prescription serves well as a tracer of for an explosion (see Fig. 13). different ﬂuid populations. ALPHA PARTICLE RECOMBINATION IN SUPERNOVAE 13 Mass-weighted Radius [rs0] 0.41 0.60 1.00 1.50 2.00 3.50 0.3 0.2 _ t = (119 + 25) _ t = (189 + 25) 0.1 b [vff02] 0 -0.1 -0.2 0.2 _ t = (264 + 25) _ t = (364 + 25) 0.1 b [vff02] 0 -0.1 -0.2 -0.3 -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 speciﬁc 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 0.3 b [vff02] It is also apparent from Fig. 12 that material with a longer 0.0 residency time tends to have lower Xα , as is expected because it also tends to have a higher temperature. Before the shock -0.3 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. 100 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 ﬂux 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 ﬂuid 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 ﬂow 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 ﬂuid. time tres in the exploding run with ε = 0.15v2 0 (see also Fig. 13). The colors ff 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 speciﬁc nuclear binding energy denuc /dt, in the ergy, ε = 0.15v2 0 . In this class of models, the infalling heavy ff 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 ﬁrst step, we average the convective ﬂow over windows of width 50tff 0 , which de-emphasizes short term ﬂuctuations in the averaged ﬂuid is lifted by buoyancy forces and accumulates in the re- velocity ﬁeld 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 downﬂow intermittently ﬂips in sign, and tive cells, across which material ﬂows 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 ﬁeld (white arrows) ff 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 ﬂuid parcels in the two hemispheres is also 0.4 intermittent, and sometimes two circulation ﬂows 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 ] 2 0.3 t = (264± 25) tff0 We have identiﬁed a useful ﬁgure 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 2 0.3 t = (975± 25) tff0 This measures the speciﬁc energy that is absorbed from neu- trinos by the material that ﬂows 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 ampliﬁcation of the mass in material with positive b. The heat absorbed at the base 0 of the convective cells is of the same order of the Bernoulli pa- 0 1 2 3 rameter of the ﬂuid, 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 speciﬁc 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 ﬁrst average vθ t over all polar angles, and then tion 50tff 0 . The upper panel shows the exploding run with ε = 0.15v2 0 and ff deﬁne 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 identiﬁed, 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α) ] 3 2D NSE 0.015 Hcr / (rs0 vff0 ) 3 0.006 0.01 0.004 0.005 0.002 0.1 0.15 0.2 0.25 0.2 0.3 0.4 0.5 0.6 2 ε / 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 ﬂow conﬁguration (v1 being the radial ﬂow 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 ﬂow pattern appears to be, in and v1 is the ﬂow speed upstream of the shock in the initial part, an artifact of the axisymmetry of the ﬂow. This im- conﬁguration (that is, in the time-independent, spherical ﬂow 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 .) 1 place within horizontal ﬂows 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 ﬂow solution can be found. with the downward velocity of the main accretion plume. We The maximal heating parameter Hsteady for a steady ﬂow cor- do observe that the main accretion plume becomes strongly responds directly to the one ﬁrst 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 ﬂuid 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 . ff 1 ˙ (GM M/r), in the idealized (but unrealistic) case where the An explosion is signiﬁcantly easier when the ﬂuid below the ﬂow is composed only of free nucleons and moves hyperson- shock starts out with a signiﬁcant 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 1 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 ﬂow solution at H = 0. Then ﬂow 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 ff 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 1 κ 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 conﬁgurations explode along one or both poles (see Figs. 11 and 13), so that less material must be lifted through the gravitational ﬁeld than in a fully spherical explo- sion. 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 ﬁxed 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) ﬁnd 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 ﬂow 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 ﬂow. 6. CONVECTION AND THE SASI Overturns of the ﬂuid 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 conﬁgurations 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 inﬂow of fresh material below the shock bi- the threshold for a neutrino-driven explosion. The implica- furcates from older shocked ﬂuid. 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 ﬁgure 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- vr 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 ﬂow models, we can translate H into a the amplitude of the shock oscillations drops signiﬁcantly as value for χ, and ﬁnd (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 ﬂow solution. the r.m.s. amplitude of the shock oscillations correlates with s For a ﬁxed 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 2 ∆ (a0 - <a0>) / rs0 0.2 2 ε = 0.15vff0 χ = 1.5 χ = 3.9 0.1 ε / vff0 = 0.15 0 0.1 a1 / rs0 -0.1 0 ∆ (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 0 0.1 -0.2 ∆ (a1-<a1>)/rs0 a2 / rs0 0.5 0 0 ∆ (a2-<a2>)/rs0 -0.5 0.5 0.5 rs0〈 ds/dr〉 θ,med 0 0 -0.5 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. ﬂuctuation of the difference between the instantaneous t / tff0 Legendre coefﬁcient 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 coefﬁcient 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- ff that the amplitude is measured in absolute units (rs0 ). ferent heating rates, corresponding to χ = 1.5 (H = 0.002 rs0 v3 0 , blue curves) ff and χ = 3.9 (H = 0.004 rs0 v3 0 , red curves). Both runs are below the threshold ff 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 proﬁle. A running average of we ﬁrst calculate the running average a 50t of the shock Leg- the entropy gradient (temporal width 20tff0 ) appears as thick solid lines. endre coefﬁcients a over an interval 50tff0 , and then calculate ˆ the r.m.s. ﬂuctuation 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 coefﬁcients 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 ﬁnd that the amplitude of the = 1 and 2 tude with increasing heating rate. This conﬁrms 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 ﬂow 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 signiﬁcantly 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 ﬁnite 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 ﬂattens. 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 ﬂatten 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 1.2 -0.1 1 t = 20 tff0 0.5 0.8 rs0〈 ds/dr〉 θ,med 0.6 0 0.4 -0.5 0.2 0.1 1.0 t = 35 tff0 a2 / rs0 Mach number 0.8 0 0.6 0.4 0.5 rs0〈 ds/dr〉 θ,med 0.2 0 -0.5 1.0 t = 60 tff0 -1 0.8 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. 0.2 0 stants in the evolution of the run with ε = 0.15v2 0 , = 1 seed ff -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 ff the same time that convection has ﬁlled 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 signiﬁcant growth. The vertical lines show the vorticity of a convective ﬂow 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- ﬁgure 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 proﬁle in the forming neutron star. Within the framework 7. SUMMARY explored in this paper, we ﬁnd two extreme types of explo- We have investigated the effects of α-particle recombi- sion. In the ﬁrst, neutrino heating does most of the work, nation and neutrino heating on the hydrodynamics of core- with a signiﬁcant ﬁnal 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 ﬁnal 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. Signiﬁcantly 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. downﬂows, so that nuclear dissociation remains on balance 7. – Vortical motions with a Mach number ∼ 0.3-0.5 ﬁrst ap- an energy sink within the accretion ﬂow. 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 modiﬁed by heating. In contrast with the = 1, 2 modes of a a source of acoustic waves, which have a similar period to laminar accretion ﬂow, 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 difﬁcult 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 difﬁcult 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 ﬁnal push. Second, an tial scales. Volume-ﬁlling 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 signiﬁcant fraction of the heating oc- these issues is possible only with high-resolution 3D simula- curs in horizontal ﬂows 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 scientiﬁc 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 ﬁnd 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 ﬂow below the shock, to a volume-ﬁlling 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. APPENDIX A. ALPHA-PARTICLE ABUNDANCE IN NUCLEAR STATISTICAL EQUILIBRIUM We calculate the α-particle mass fraction Xα in nuclear statistical equilibrium by limiting the nuclear species to α-particles and free nucleons, and ﬁxing 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 ﬂow 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 h 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 3 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) 1 Xp + Xα = Ye . (A5) 2 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 ] -3 29 (a) (e) Prel / Pmat 1 28 10 η 27 0.1 26 (c) 1 ∂ Xα /∂ ln p, ∂ Xα /∂ ln ρ γ = (∂ ln p/∂ ln ρ)S (d) 1.6 S = 10 kB/nucleon (f) 3 (b) 1 kBT [MeV] eq 1.4 2 0 1 -1 1.2 eq 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 ﬂuid containing n, p, α, photons, and ﬁnite-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 eq a function of density, for ﬁxed 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 eq 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). eq 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 conﬁguration 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. B. TIME-INDEPENDENT FLOW EQUATIONS DESCRIBING INITIAL MODELS We write down the ordinary differential equations that are used to compute the initial ﬂow, and the density proﬁles 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 sufﬁce 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 and dp dρ + Aρ . Lα ≡ A p (B5) dr dr The coefﬁcients 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 ) r = . (B6) dr ρvr Eρ − pvr /ρ − Aρ + v2 ρvr E p − A p r The coefﬁcients in eqs. (B4) and (B5) work out to 1 Ep = (constant γ) (B7) (γ − 1)ρ p 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 speciﬁc internal energy is 1 3 p 3 3 kB T eint = 3prel + pmat =3 − 1 − Xα . (B9) ρ 2 ρ 2 4 mu One therefore ﬁnds 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 ﬂow 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 ﬁxed, 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) i r∗ where the sum is taken over the computational cells below the shock at our ﬁxed 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 coefﬁcient 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 . eq 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 from M˙ ρ1 (rs ) = 2 |v (r )| , (B13) 4πrs 1 s where vff (rs ) v1 (rs ) = − (B14) 1 + 2M−2 (rs )/(γ − 1) 1 is the upstream velocity at r = rs , while the upstream pressure satisﬁes ρ1 (rs )[v1 (rs )]2 p1 (rs ) = . (B15) γM2 (rs ) 1 ˙ 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 modiﬁcations described in Paper I. The rate of change of speciﬁc 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 ﬁxed 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 ﬂow. 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 veriﬁed 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.) REFERENCES 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

DOCUMENT INFO

Shared By:

Categories:

Tags:

Stats:

views: | 2 |

posted: | 10/20/2011 |

language: | English |

pages: | 22 |

OTHER DOCS BY hedongchenchen

How are you planning on using Docstoc?
BUSINESS
PERSONAL

By registering with docstoc.com you agree to our
privacy policy and
terms of service, and to receive content and offer notifications.

Docstoc is the premier online destination to start and grow small businesses. It hosts the best quality and widest selection of professional documents (over 20 million) and resources including expert videos, articles and productivity tools to make every small business better.

Search or Browse for any specific document or resource you need for your business. Or explore our curated resources for Starting a Business, Growing a Business or for Professional Development.

Feel free to Contact Us with any questions you might have.