THE DEPENDENCE OF THE NEUTRINO MECHANISM OF CORE-COLLAPSE SUPERNOVAE ON THE EQUATION OF STATE

The exact process that halts the collapse of the core of a massive star at the end of its life and drives a successful supernova explosion, releasing a multitude of neutrinos and ejecta with around 10⁵¹ erg of kinetic energy, is not yet fully-understood. The delayed neutrino-heating mechanism of core-collapse supernovae (CCSNe; Colgate & White 1966; Bethe & Wilson 1985), the leading candidate for the explosion mechanism, fails to consistently produce energetic explosions for the wide range of possible progenitor masses in simulations that include high-fidelity treatments of neutrino transport. Neutrinos are an attractive explanation for the explosion mechanism since the bulk of the gravitational binding energy of the progenitor core, some 10⁵³ erg, will be radiated away as neutrinos as the proto-neutron star (PNS) forms and cools; only about 10⁵¹ erg of this energy needs to be transferred to the collapsing stellar material in order to explain typical CCSNe energies. Only a fraction, however, of the core's original gravitational binding energy, perhaps a few times 10⁵² erg, will be released as neutrino radiation in the first second following the collapse and bounce of the core, approximately the timescale on which a successful explosion must occur. Neutrino radiation-hydrodynamic simulations of core collapse show that this is challenging to achieve, particularly in 1D simulations. It is now clear that multidimensional effects, such as PNS convection (Epstein 1979; Burrows & Fryxell 1993; Dessart et al. 2006; Burrows et al. 2007a), neutrino-driven convection (Janka & Mueller 1996), the standing accretion shock instability (SASI, Blondin et al. 2003), and turbulence (Murphy & Meakin 2011), play a critical role in the success of the neutrino mechanism in driving explosions. Together, these multidimensional effects can push some progenitors over the critical threshold, resulting in somewhat marginal explosions in 2D in certain simulations (Marek & Janka 2009; Bruenn et al. 2009; Suwa et al. 2010; Mueller et al. 2012). This result is dependent on details of the numerical scheme and treatment of neutrino transport; the 2D simulations of Burrows et al. (2006, 2007b) do not find neutrino-driven explosions for any progenitors.² Based on a series of general relativistic 1D simulations including Boltzmann transport for neutrinos, Lentz et al. (2012) suggest that the discrepancies in the results from the various groups simulating CCSNe may be due in large part to differences in how three major aspects of the physics are treated: general relativistic versus Newtonian dynamics; the neutrino interactions and opacities considered, particularly the inclusion of inelastic neutrino scattering; and observer frame corrections in the neutrino transport.

All of this indicates, unsurprisingly, that simulations of the neutrino mechanism are sensitive to microphysics, and details of how the microphysics is treated. It is still uncertain, however, if the variation in the handling of microphysics is large enough to explain the dearth of energetic explosions; that is, if a more accurate handling of the details of the neutrino physics will result in successful explosions with around 10⁵¹ erg of kinetic energy for a wide range of progenitor masses. Nordhaus et al. (2010) have suggested a larger impact on the success of the neutrino mechanism results from fully 3D simulations. Using a simplified neutrino heating and cooling scheme rather than expensive neutrino transport, Nordhaus et al. (2010) find that the driving neutrino luminosity necessary to obtain an explosion is about 20% lower in 3D than in 2D. This difference may be larger than what we could expect to gain from higher-fidelity treatments of neutrino transport or the inclusion of general relativity in simulations. Hanke et al. (2012) sought to reproduce this result. Employing a neutrino heating/cooling approach similar to Nordhaus et al. (2010), Hanke et al. (2012) find that the critical neutrino luminosity in 3D is not significantly lower than that of 2D. This raises the questions of just how beneficial 3D simulations are to the success of the neutrino mechanism and what is the cause for the differences in the results of Nordhaus et al. (2010) and Hanke et al. (2012)? As 3D simulations become more feasible and available, we may be able to answer the question of the importance of 3D. Recently, Takiwaki et al. (2012) reported an explosion for an 11.2 M_☉ progenitor in a 3D simulation performed using the isotropic diffusion source approximation for neutrinos (IDSA; Liebendörfer et al. 2009). Compared to 2D, Takiwaki et al. found that the convection below the gain region was more vigorous in 3D enhancing the neutrino luminosity beyond the 2D case. The somewhat modest resolution they were able to afford, however, means that the issue of the importance of 3D simulations cannot be settled based on of their results alone.

Besides the neutrino physics, an important microphysical dependence of the CCSN mechanism is the equation of state (EOS). The EOS for matter at densities and temperatures relevant to CCSNe is still an active area of research (e.g., Hempel & Schaffner-Bielich 2010; Shen et al. 2010, 2011a, 2011b, 2011c; Furusawa et al. 2011). Our understanding of nuclear interactions at these densities and temperatures is incomplete and a number of models exist. Testing the various EOS models with experiment and observation is challenging, but a few constraints do exist. Hempel et al. (2012) review some of the constraints on the EOS parameters from laboratory experiments. In Table 1, we summarize the constraints on the nuclear EOS parameters and the respective numbers for five EOS models used in CCSN simulations: Lattimer & Swesty with incompressibility parameters of 180 MeV (LS180), 220 MeV (LS220), and 375 MeV (LS375), Shen et al. (STOS), and the recent EOS of Hempel & Schaffner-Bielich (2010) using the TMA model for nuclear interactions (HS (TMA)). Most values for the limits are taken from experiment (for references to the experimental data see Hempel et al. 2012), except for the skewness coefficient, K', which is only a theoretical estimate. Based on the measurements listed in the table, only LS220 falls within or near the limits for each of the EOS parameters.

Observations of neutron stars and pulsars can also place constraints on the EOS for hot, dense matter. Using a parameterized EOS, Steiner et al. (2010) analyze a small sample of neutron stars and determine a most likely mass–radius relationship for cold neutron stars. Their results favor a radius of 11–12 km for a 1.4 M_☉ neutron star, and a maximum neutron star mass in the range 1.9–2.2 M_☉. In more recent work, Steiner et al. (2012) place broader limits of the radius of a 1.4 M_☉ NS at 10.4–12.9 km, while the analysis of Lattimer & Lim (2012) suggests very tight limits on NS radii of 11–12 km. The maximum mass of a cold neutron star is an important constraint for the EOS models. The observations of PSR J1614-2230 by Demorest et al. (2010) place a very tight limit on the mass of the neutron star in this system of 1.97 ± 0.04 M_☉. Given an EOS model, solutions to the Tolman–Oppenheimer–Volkov equation yield the gravitational mass–radius relationship for cold neutron stars. Figure 1 shows the mass–radius relationships for the EOS of Lattimer & Swesty (1991) and Shen et al. (1998), along with the mass limits for PSR J1614-2230. Based on the measurement of Demorest et al. (2010), LS180 is ruled out as it does not produce a large enough maximally-massed neutron star. The Shen et al. EOS produces a maximum neutron star mass well above the limit placed by PSR J1614-2230 (2.22 M_☉), but STOS predicts a radius for a 1.4 M_☉ NS of around 15 km, well outside the limits estimated by Steiner et al. (2010). The theoretical predictions for the radius of a 1.4 M_☉ NS of Hebeler et al. (2010) are somewhat broader, 9.7–13.9 km, but STOS still lies outside of this range.³ Only LS220 satisfies the constraint on the maximum NS mass from PSR J1614-2230 and produces a 1.4 M_☉ NS with a radius (∼12.7 km) close to the upper estimate of Steiner et al. (2010) and well within the limits of Hebeler et al. (2010). Figure 1 also shows the mass–radius relationship for LS375. This EOS easily satisfies the requirement on the maximum neutron star mass, but falls above the radius limits of Steiner et al. (2010) and has an incompressibility (375 MeV) well above that found experimentally. We note that the symmetry energy coefficient, J, is not very well-constrained experimentally and the limit given in Table 1 should be taken with caution (see discussion in Lattimer & Lim 2012).

Zoom In Zoom Out Reset image size

In this work, we study the dependence of the neutrino mechanism on the EOS employed. Using 1D and 2D hydrodynamic simulations with simplified neutrino transport, we determine the explosion times for a 15 M_☉ progenitor as a function of neutrino luminosity for three EOS: Lattimer & Swesty (1991) with K = 180 MeV and K = 220 MeV, and Shen et al. (1998). In 1D only, we also run simulations using LS375. Our treatment of the neutrino physics is similar to those used by Murphy & Burrows (2008), Nordhaus et al. (2010), and Hanke et al. (2012): we assume local neutrino heating and cooling, based on the rates derived by Janka (2001), and a constant neutrino luminosity. We find that the explosion times for a given neutrino luminosity are significantly dependent on the EOS used, with the general trend of easier explosions for softer EOS. Thus, the Lattimer & Swesty models result in easier explosions than the Shen et al. EOS, with LS180 leading to the earliest explosions at a given neutrino luminosity. Early 1D simulations using a simple, parameterized EOS showed that stronger shocks result for softer EOS (Baron et al. 1985a, 1985b). Our results show that this trend is recovered for more realistic nuclear EOS and in 2D simulations.

The EOS has been shown to influence the results of CCSN simulations in a number of contexts. Thompson et al. (2003) explored the sensitivity of their 1D neutrino radiation hydrodynamic simulations to the value of K for the Lattimer & Swesty EOS. They found that the overall evolution during the first 200 ms post-bounce is not very dependent on K: the resulting temperatures were slightly higher and the emergent neutrino luminosities were at most 9% larger for LS180 as compared to LS375. Sumiyoshi et al. (2005) compared LS180 and STOS in 1D general relativistic simulations. While neither model explodes in their simulations, they find differences in composition resulting from the different symmetry energies (see J in Table 1) and in temperature evolution, with LS180 giving slightly larger temperatures and, therefore, slightly higher neutrino luminosities. The higher neutrino luminosity and heating rate found for LS180 is somewhat counteracted by the more compact proto-neutron, a result of both lower incompressibility and smaller symmetry energy slope, L. Sumiyoshi et al. (2005) find that the differences between the two models grow with time, particularly after 300 ms post-bounce. Marek et al. (2009) explored the dependence on the gravitational wave signal produced by 2D simulations on the EOS, comparing LS180 and the stiffer EOS of Hillebrandt & Wolff (1985). They find that the matter-generated gravitational wave signal is not very dependent on the EOS, but that the differences in the anisotropic neutrino emission gravitational wave signal, which can dominate the matter gravitational wave signal, are significantly different between LS180 and Hillebrandt & Wolff EOS. In failed explosions, the formation of black holes has been found to be highly dependent on the EOS (Sumiyoshi et al. 2007; Fischer et al. 2009; O'Connor & Ott 2011; Hempel et al. 2012). Generally, a stiffer EOS results in a greater time to black hole formation following bounce.

This paper is organized as follows. In Section 2, we describe the details of our numerical approach. In Section 3, we describe the results of our simulations and the dependence of the explosions on the EOS. In Section 4, we conclude and discuss how our work relates to recent work on neutrino-driven supernovae.

We conduct our study using the FLASH⁴ simulation framework (Dubey et al. 2009) to solve the Eulerian equations of hydrodynamics,

$$\begin{eqnarray} \frac{\partial \rho }{\partial t} + \pmb {\nabla } \cdot (\rho \pmb { v}) &= 0, \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \frac{\partial \rho {\pmb {v}}}{\partial t} + \pmb {\nabla } \cdot (\rho \pmb { v v}) + \pmb {\nabla } P &= -\rho \pmb {\nabla } \Phi, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \frac{\partial \rho E}{\partial t} + \pmb {\nabla } \cdot [(\rho E + P){\pmb {v}}] &= \rho {\pmb {v}} \cdot \pmb {\nabla } \Phi + \rho (\mathcal {H-C}), \end{eqnarray} \tag{ 3 }$$

where ρ is the mass density, v the velocity vector, P the pressure, Φ the gravitational potential, E the total specific energy, and $(\mathcal {H-C})$ is the difference between neutrino heating and cooling. We use the directionally-unsplit hydrodynamics solver in FLASH to solve Equations (1)–(3) (D. Lee et al., in preparation; see also Lee & Deane 2009). The unsplit solver in FLASH has many options for the spatial reconstruction scheme, Riemann solver, slope limiter, etc. For this work, we use third-order piecewise-parabolic spatial reconstruction (PPM; Colella & Woodward 1984), a "hybrid" slope limiter, and the HLLC Riemann solver. Inside of shocks, we use the HLLE Riemann solver in order to avoid the odd–even decoupling instability in 2D (Quirk 1994). This approach allows us to avoid the smearing of contact discontinuities that the HLLE solver is prone to while benefiting from the HLLE solver's immunity to odd–even decoupling at shocks. The "hybrid" slope limiter (Balsara 2004) applies the monotonized central (mc) limiter to linear wave families (such as density) and the robust, slightly more diffusive minmod limiter to non-linear, self-steepening wave families. The gravitational potential, Φ, is calculated using the new multipole Poisson solver available in FLASH4 (see Section 8.10.2.2 of the FLASH4 User's Guide), however, we set ℓ_max = 0, yielding monopole gravity.

Following the approaches of Murphy & Burrows (2008), Nordhaus et al. (2010), and Hanke et al. (2012), we use the neutrino heating and cooling rates derived by Janka (2001),

$$\begin{eqnarray} \mathcal {H} &=& 1.544\times 10^{20} \left(\frac{L_{\nu _e}}{10^{52}\ {\rm erg\ s^{-1}}} \right) \left(\frac{T_{\nu _e}}{4\ {\rm MeV}} \right)^2 \nonumber \\ &&\times \left(\frac{100\ \rm km}{r} \right)^2 (Y_p + Y_n) e^{-\tau _{\nu _e}}\ \left[ \rm {\frac{erg}{g \ s}} \right], \end{eqnarray} \tag{ 4 }$$

and

$$\begin{equation} \mathcal {C} = 1.399\times 10^{20} \left(\frac{T}{2\ \rm MeV} \right)^6 (Y_p + Y_n) e^{-\tau _{\nu _e}}\ \left[ \rm { \frac{erg}{g\ s}} \right], \end{equation} \tag{ 5 }$$

where $L_{\nu _e}$ is the electron neutrino luminosity (note it is assumed that $L_{\bar{\nu }_e}=L_{\nu _e}$), $T_{\nu _e}$ is the electron neutrino temperature, r is the spherical radius, (Y_p + Y_n) is the sum of the neutron and proton number fractions, $\tau _{\nu _e}$ is the electron neutrino optical depth, and T is the matter temperature. In all of our calculations, $T_{\nu _e}$ is set to 4 MeV. The heating term, Equation (4), is applied only after core bounce has occurred in our simulations. The factor $e^{-\tau _{\nu _e}}$ is included to provide an exponential cutoff of the heating and cooling at high density where the material is optically-thick to neutrinos. As such, we do not calculate the optical depth self-consistently but instead rely on a piece-wise fit of the optical depth as a function of density. In Figure 2, we compare our piece-wise fitting approach for $\tau _{\nu _e}$ with the integrated electron neutrino optical depth,

$$\begin{equation} \tau _{\nu }(r) = {-}\int ^{r}_\infty \kappa _{\rm eff} dr, \end{equation} \tag{ 6 }$$

computed using the effective opacity given by Janka (2001), $\kappa _{\rm eff} = 1.5\times 10^{-7} \rho _{10} {({1}/{2})}(Y_p+Y_n)(k T_{\nu _e}/4\ {\rm MeV})^2$ cm⁻¹. As shown, the piece-wise fit closely approximates the actual integrated optical depth. We also show in Figure 2 an exponential cutoff based directly on the density. We see that our piecewise approach better approximates the integrated optical depth than this simple density-based approach. We stress that the factor $e^{-\tau _{\nu _e}}$ in Equations (4) and (5) is used simply as a cutoff for the heating and cooling functions at high-density and the exact values of the cutoff are less important than consistent application and use in all of the models considered in a study such as the present one.

Zoom In Zoom Out Reset image size

We calculate the evolution of the electron fraction, Y_e, according to the parameterization of Liebendörfer (2005) in which Y_e is dependent only on density. This approach is appropriate only for the pre-bounce collapse phase, however, we employ the parameterization post-bounce as well in order to maintain consistency with Nordhaus et al. (2010). With no post-bounce evolution of the electron fraction, the cooling above the neutrinosphere is grossly underestimated, leading to significantly earlier and easier explosions, as mentioned by Hanke et al. (2012). We do not include the entropy increases due to neutrino thermalization as prescribed by Liebendörfer (2005).

We run our simulations in 1D spherical and 2D cylindrical geometries using adaptive mesh refinement (AMR) as implemented in FLASH via PARAMESH (v.4-dev; MacNeice et al. 2000). For this work, we have extended FLASH's unsplit solver to spherical and cylindrical coordinates using an approach similar to those described in Skinner & Ostriker (2010) and Appendix A of Mignone et al. (2005). Our fiducial set of simulations has a minimum grid-spacing of 0.5 km in each direction, obtained using 8 AMR levels. The grid in our 1D spherical simulations covers radii from 0 km to 5000 km. In 2D cylindrical runs, we simulate a grid that spans 5000 km in R and 10,000 km in z. We do not remove the PNS from our grid and follow the evolution from progenitor phase through collapse and bounce and into the shock stagnation/revival phase. For boundary conditions at the outer edge of the domain, we assume radial power-law profiles that approximate the stellar envelope. We find that simple "outflow" boundary conditions, which enforce a zero-gradient condition, are not appropriate for the core-collapse context because the amount of stellar material entering the outer edge of the domain will be overestimated, artificially enhancing the mass accretion rate at late times. By experimentation, we have found that "outflow" boundary conditions can, in this way, suppress explosions for neutrino luminosities near critical. For all of our simulations, we use the 15 M_☉ progenitor of Woosley & Weaver (1995).

We have incorporated the finite temperature EOS routines of O'Connor & Ott (2010) into the FLASH framework.⁵ We use three different EOS models in our simulations, the models of Lattimer & Swesty (1991) with incompressibility, K, of 180 MeV and 220 MeV, and that of Shen et al. (1998). Some parameters for these EOS are listed in Table 1.

3.1. Dependence on EOS

We ran a series of core collapse simulations in 1D and 2D in which we varied the driving neutrino luminosity. For models that explode, we measure the post-bounce time of the explosion and the mass accretion rate at the time of explosion. We consider a model to have exploded if the average shock radius exceeds 400 km and does not subsequently fall back below 400 km. We measure the mass accretion rate at a radius of 500 km in both 1D and 2D simulations, though the measured mass accretion rates are identical since the 2D solution remains spherically-symmetric outside of the shock. Figure 3 shows our measured mass accretion rates as a function of time post-bounce for the three EOS we consider. Our mass accretion rate is very similar to that of Hanke et al. (2012) for their 15 M_☉ progenitor model. The differences in mass accretion history between STOS and LS can be attributed to models using STOS collapsing and bouncing a little faster than models using LS; bounce occurs 50 ms earlier for STOS as compared to LS. Table 2 summarizes the simulations we ran and the resulting explosion times and mass accretion rates at the time of explosion. In Figure 4, we plot the driving neutrino luminosities as a function of explosion time and mass accretion rate at the time of explosion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Time Post-bounce and Accretion Rates at Time of Explosion for the 15 M_☉ Progenitor with 0.5 km Resolution

	LS180		LS220		STOS
$L_{\nu _e}$a	texpb	$\dot{M}_{\hbox{exp}}$c	texp	$\dot{M}_{\hbox{exp}}$	texp	$\dot{M}_{\hbox{exp}}$
(1052 erg s−1)	(ms)	(M☉ s−1)	(ms)	(M☉ s−1)	(ms)	(M☉ s−1)
1D
2.0	...	...	...	...	...	...
2.1	...	...	...	...
2.2	689	0.192	816	0.169	...	...
2.3	374	0.254	571	0.211	943	0.153
2.4	213	0.289	402	0.251	829	0.170
2.5	193	0.312	196	0.311	380	0.262
2.7	175	0.317	178	0.320	216	0.310
2.9					200	0.314
2D
1.3	813	0.182	...	...	...	...
1.5	277	0.279	368	0.262	732	0.192
1.7	203	0.307	246	0.295	374	0.263
1.8			194	0.313
1.9	178	0.315			249	0.298
2.0	179	0.318	171	0.322	217	0.313
2.1			171	0.322	201	0.317
2.2			167	0.324	193	0.322
2.4					183	0.322
2.6					168	0.328

Notes. ^aElectron-neutrino luminosity. ^bTime after bounce of onset of explosion. A "..." symbol indicates that the model does not explode during the simulated period of evolution. ^cMass accretion rate at onset of explosion.

Download table as:  ASCIITypeset image

Our results show that the Lattimer & Swesty EOS explode more easily than that of Shen et al., with LS180 resulting in the earliest explosions for a given neutrino luminosity (lowest curves in Figure 4). The results then follow a basic trend where the stiffer the EOS is, the harder it is to drive an explosion (LS180 being the softest EOS and STOS being the stiffest EOS we consider). This trend holds in both 1D and 2D simulations. For high neutrino luminosities resulting in earlier explosions, the curves in Figure 4 for LS180 and LS220 converge, showing very little difference in time of explosion at these luminosities. Figure 5 shows the (average) shock radii in 1D (2D) for the three EOS and several neutrino luminosities. The 1D simulations show a clear delineation of models into the three categories of the radial shock instability (Fernández 2012): oscillatory/stable,⁶ oscillatory/unstable, non-oscillatory/unstable. In 2D, ℓ > 0 modes of the SASI dominate the purely radial (ℓ = 0) shock oscillations seen in 1D, leading to a more chaotic variation in the average shock radius in time. In Figure 6, we compare the radial entropy profiles in 1D and 2D for LS220 and STOS at 50 ms and 200 ms post-bounce for $L_{\nu _e, 52} = L_{\nu _e}/(10^{52}\ {\rm erg\ s^{-1}}) = 2.0$. For this luminosity, both 2D simulations result in rather prompt explosions, whereas the 1D simulations fail to explode entirely, with the shock stalling around 200 km. The entropy scatter plots for the 2D simulations shown in Figure 6 also demonstrate the entropy re-distribution due to convection and turbulence.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 7–9 show the evolution of the entropy in 2D for $L_{\nu _e} = 1.3\times 10^{52}\ {\rm erg\ s^{-1}}$. For this luminosity, only LS180 results in an explosion within 1 s post-bounce. Figure 7 corresponds to 100 ms post-bounce. At this stage, convection has set in at the post-shock region as well as the PNS cooling region (faint blue area in the figure). Some extension along the axis is evident, particularly for LS180, but a strong ℓ = 1 SASI is not yet apparent. By 300 ms (Figure 8), however, the SASI is very obvious and is reaching large amplitude for LS180, though less so for LS220 and STOS. At 600 ms (Figure 9), large-amplitude ℓ = 1 SASI motions are clear for LS180 and LS220, but the SASI amplitude remains modest for STOS. While LS180 and LS220 show shock expansion, particularly along the axis, STOS shows a fairly stationary shock, in an average sense, at around 200 km.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our results show a clear trend with the incompressibility, K, of the EOS: higher values of K result in later explosions. In order to test the importance of K in determining the time to explosion, we have also run a series of 1D simulations using the Lattimer & Swesty EOS with K = 375 MeV. The results are shown in Figure 4, where we plot the neutrino luminosity versus mass accretion rate at time of explosion for our 1D simulations. As can be seen, the trend with K does not survive: the simulations with LS375 explode earlier (at higher $\dot{M}$) than STOS. At low-luminosities, near critical, LS375 explodes later than LS220 and LS180. At higher luminosities, LS375 roughly converges with LS220, and all three versions of the LS EOS converge at $L_{\nu _e,52} = 2.5$. This indicates that the time to explosion for a given EOS is not dependent on K alone, but must also depend on the higher-order EOS parameters listed in Table 1.

Why, given an identical neutrino luminosity, do differences in the baryonic EOS result in different evolutions of the supernova shock? That is, what mechanism couples the behavior of the PNS, where the differences between the EOS are important, to the motion of the shock? The difference in time-to-explosion we find for the various EOS is likely due to EOS-dependent behavior of the SASI. Whether the underlying SASI mechanism is solely acoustic (Blondin & Mezzacappa 2006) or advective-acoustic (Foglizzo et al. 2007), the magnitude of the acoustic energy flux from the PNS is a determining factor in its growth (Guilet & Foglizzo 2012).⁷ A soft EOS, such as LS180 and LS220, results in a more compact, higher-density PNS as compared to a stiff EOS like STOS. This is true for the central core of the PNS as well as the "deceleration" region at the edge of the PNS where advected perturbations are turned into sound waves. The EOS dependence we find can likely be explained by a difference in the acoustic flux generated at the edge of the PNS: the denser PNSs produced by LS180 and LS220 generate a larger outbound acoustic flux than that of STOS and, thus, lead to larger amplitude excitations of the SASI. Dense PNSs are, in effect, better sonic transducers. This will hold true in both 2D and 1D, where the SASI is purely radial. Differences in the growth of the SASI also account for the similarity in time-to-explosion for LS180 and LS220 at high $L_{\nu _e}$: these prompt explosions do not depend so much on the SASI.

To illustrate the EOS dependence of the PNS structure, we show in Figure 10 the density and temperature profiles from 1D simulations with L_{ν, 52} = 2.3 at the time of explosion. Among the LS EOS, the central PNS density increases with the incompressibility parameter while STOS yields the smallest central PNS density. Note that the central densities we find will be smaller than those from simulations that include more realistic neutrino transport since the prescription we use prevents the PNS from cooling at all. The most important region of the profiles generating more of the acoustic flux associated with the growth of the SASI is in the deceleration region outside the edge of the PNS around 75 km. This radius roughly corresponds to the neutrinosphere, at a density of around 10¹¹ g cm⁻³. At this radius, there is a clear ordering in density that corresponds to the inverse of the ordering in explosion time we find: STOS < LS375 < LS220 < LS180. A simple measure of the acoustic flux moving against the accretion flow at some radius is given by (Foglizzo & Tagger 2000),

$$\begin{equation} F^- = \dot{M} c_S^2 \frac{(1-\mathcal {M})^2}{\mathcal {M}} \left(\frac{\delta P^-}{\gamma _c P} \right)^2, \end{equation} \tag{ 7 }$$

where $\dot{M}$ is the mass accretion rate, c_S is the adiabatic sound speed, $\mathcal {M}$ is the mach number, P is the pressure, δP⁻ is the incoming pressure perturbation, and γ_c is the adiabatic index. In Figure 11, we plot this outward-propagating acoustic flux for 1D simulations at a radius of 75 km assuming a constant pressure perturbation amplitude δP⁻/P ∼ 1%. We find a clear trend in the magnitude of this measure of the acoustic flux and explosion delay times amongst the different EOS: larger typical values of F⁻₇₅ correlate with earlier explosion times. Even for non-exploding models (top panel of Figure 11), the "softer" EOS have greater acoustic fluxes. The peaks in the acoustic flux correspond with minima of the shock radius. The acoustic flux curves for LS375 and STOS are very similar, though the peaks for LS375 reach larger amplitudes. This suggests that the difference between the explosion times of LS375 and STOS (see Figure 4) is not due entirely to a greater acoustic flux for LS375. Here, the difference in mass accretion history between the EOS plays a role. The slightly longer period of iron core accretion (see Figure 3), due to different bounce times between LS and STOS, delays explosion for STOS relative to LS. We note also that since our prescription for the neutrino physics prevents the PNS from cooling, it does not contract as quickly as it would in a more realistic simulation. This could affect the emergent acoustic flux.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The higher mass accretion rate of STOS at early times (Figure 3) can also account for differences in time-to-explosion for high-$L_{\nu _e}$ simulations. This higher accretion rate prior to 100 ms post-bounce results in a greater mass of iron being dissociated at the shock and, therefore, more energy loss. In the high-$L_{\nu _e}$ models, the shock begins running away to large radii at the moment the mass accretion rate drops precipitously (see Figure 5), corresponding to the end of the accretion of the iron core. Since this is delayed slightly for STOS, the shock begins its runaway expansion later as compared with LS180 and LS220. This difference will be less important at later times and cannot explain the differences in time-to-explosion between LS180 and LS220 (and LS375 in 1D) which we find for lower-$L_{\nu _e}$ models.

The EOS-dependence of the acoustic flux from the PNS atmosphere, and the concomitant feedback on the growth of the SASI, impacts the turbulent energy on large scales in our 2D simulations. In Figure 12, we plot the non-radial part of the kinetic energy as a function of spherical harmonic order, ℓ. We follow Hanke et al. (2012) and define the ℓ-dependent energy as

$$\begin{equation} E(\ell) = \sum _{m=-\ell }^\ell \left|\oint \sqrt{\rho } v_\theta Y_\ell ^m(\theta ,\phi) d\Omega , \right|^2, \end{equation} \tag{ 8 }$$

where the Y^m_ℓ are the spherical harmonics and v_θ is the velocity in the θ-direction. We evaluate the transverse kinetic energy in a ∼10 km shell centered at 150 km and average over 10 ms centered at 200 ms post-bounce to reduce noise in the spectra. Similar to Hanke et al., we see two different power-law scalings of the energy. Above the driving scale of about ℓ ≈ 20, the spectra follow a power-law of ∼ℓ⁻³, consistent with the forward enstrophy cascade. Below ℓ ≈ 20, the spectra follow more closely a scaling of ℓ^−5/3, reflecting the inverse energy cascade (Kraichnan 1967). At large ℓ (i.e., small scales), there is very little difference among the spectra for different EOS, it is only at small ℓ's that significant and important differences appear. The energy on these scales follows the same trend as the acoustic flux: LS180 > LS220 > STOS. Kinetic energy on large scales is conducive to explosion (see Hanke et al. 2012), and we would thus expect models with greater energy at small ℓ's to explode earlier, which is precisely the trend we see. This is also consistent with the EOS-dependence of the acoustic flux we show in Figure 11 as the SASI grows fastest for small ℓ.

Zoom In Zoom Out Reset image size

3.2. Criterion for Transition to Runaway Shock Expansion

Recently, Pejcha & Thompson (2012) and Fernández (2012) have studied the conditions that lead to explosions in simple, parameterized models of stalled shocks in the CCSN context. Using principles of stellar pulsation theory, Fernández (2012) develops an instability criterion that closely resembles the critical condition found in previous simulations: a greater advection time through the gain region than through the cooling region. Pejcha & Thompson (2012) find a criterion based on the ratio of the sound speed to the escape velocity, the so-called "ante-sonic condition." The ante-sonic condition for a polytropic EOS is

$$\begin{equation} {\rm max} \left(\frac{c_S^2}{v_{\rm esc}^2} \right) = 0.18 \Gamma , \end{equation} \tag{ 9 }$$

where v_esc is the escape velocity, and Γ is the polytropic index. In Figure 13, we plot the quantity (c²_S/v_esc²)/(0.18γ_c) as a function of space and time for the 1D simulation using the STOS EOS with $L_{\nu _e,52} = 2.3$, along with the shock radius as a function of time for this model. For the general EOS, we substitute the spatially-varying adiabatic index as returned by the EOS, γ_c, for the polytropic index in the ante-sonic condition. In Figure 13, instability of the type described by Pejcha & Thompson (2012) occurs for values greater than one. The highest values for the ante-sonic relation occur when the shock reaches radial maxima, and thus we plot the times when the shock has just started to collapse back during its radial oscillations. As shown, the ante-sonic condition is only satisfied on the last, explosive expansion of the shock. This is true for all of our 1D models except the $L_{\nu _e,52}=2.4$ model using STOS; in this model, we see the shock expand well beyond 400 km and then fall back inward one last time before runaway shock expansion occurs (see Figure 5). During the shock's first expansion beyond 400 km in this model, the ante-sonic condition is met but the shock nevertheless contracts again.

Zoom In Zoom Out Reset image size

In 2D, we find that the ante-sonic condition is generally met in post-shock regions that have expanded to large radii. This can happen long before the average shock radius exceeds 400 km when we consider the model to have successfully "exploded." For radial rays along which the ante-sonic condition has been met at some point, the shock radius tends to remain above 400 km.

3.3. Dependence on Resolution

Hanke et al. (2012) report that their simulations, which employ the same parametric neutrino heating/cooling as this study, are significantly resolution-dependent. In their 2D simulations, finer resolution led to earlier explosions for the same driving neutrino luminosity. We have examined the dependence of both our 1D and 2D results on resolution and show the results in Figure 14. While the development of small-scale instabilities and turbulence will depend heavily on the resolution used, for the present study, the important question is whether or not changing the resolution significantly alters the time at which a model explodes for a given neutrino luminosity. To determine this, we ran several models using STOS with an increased minimum grid spacing at the maximum level of refinement of 0.7 km, versus 0.5 km for our fiducial models. We find that the results for explosion times as a function of neutrino luminosity are essentially identical in 1D and not significantly different in 2D. As shown in Figure 14, the 1D curves lie on top of one another and the 2D curves are very similar. For 2D, the small differences are that at lower driving luminosities (and, hence, later explosion times), the low-resolution models explode earlier, while at higher luminosities the lower resolution models explode later. Thus, the high-resolution and low-resolution curves cross each other in Figure 14.

Zoom In Zoom Out Reset image size

It is difficult to make a direct comparison between our results and those of Hanke et al. (2012) because the nature of the computational grids employed is quite different. Hanke et al. (2012) use 2D spherical-polar coordinates with non-equidistant spacing in radius. Here, we use 2D cylindrical R−z coordinates with adaptive mesh refinement.⁸ Hanke et al. use a fiducial resolution of 3° in angle and 400 zones in radius. At a radius of 100 km, near the base of the gain region, this yields a linear grid spacing in the angular direction of 1.67 km, more than three times our fiducial resolution. For a spherical grid, this grid spacing will increase with radius, while in our simulations we allow refinement to the highest level anywhere in the computational domain, allowing us to achieve 0.5 km resolution for the entire post-shock region at all times. The highest resolution simulations of Hanke et al. use 05 angular resolution, resulting in 0.28 km linear resolution at 100 km and 1.1 km linear resolution at 400 km. This resolution is likely comparable to our fiducial resolution, the linear grid spacing being smaller at the base of the gain region and larger at radii where shock oscillations reach the greatest extent.

We should not expect hydrodynamic simulations that involve instabilities and turbulence, such as these, to converge in spatial resolution, but the fact that we do not see significant differences in the explosion times at different resolutions indicates that our results are converged with respect to the result we seek.

We have conducted a parameter study exploring the dependence of the neutrino mechanism of core-collapse supernovae on the EOS. We varied the driving neutrino luminosity and EOS in several 1D and 2D simulations of stellar core collapse and measured the resulting explosion times and mass accretion rates at the time of explosion. Table 2 and Figure 3 summarize our results. We find that for every EOS, 2D explosion are obtained more easily than in 1D, and that models using the Lattimer & Swesty EOS explode more easily than models using the EOS of Shen et al., with LS180 resulting in slightly earlier explosions than LS220. Thus, our results show a general trend that softer EOS lead to easier explosions. The "stiffness" of the EOS cannot be characterized solely by the nuclear incompressibility, K, as exemplified by our 1D simulations using LS375. This model has a substantially higher incompressibility than STOS (375 MeV versus 281 MeV). We find that STOS, despite having a lower incompressibility, explodes later than LS375 for all neutrino luminosities. This highlights that the higher-order EOS parameters, such as symmetry energy coefficient and slope, play an important role in determining the stiffness of an EOS. Our results for explosion delay times follow a trend in resulting neutron star radii (see Figure 1). EOS that produce larger typical neutron star radii explode later in our simulations. The dependence of the time-to-explosion on the EOS is likely due to the dependence of the acoustic energy flux from the PNS on the stiffness of the EOS. The softer EOS of Lattimer & Swesty result in more compact PNS (Figure 10) that are more effective at translating advected perturbations into outgoing sound waves (Figure 11). This then results in more acoustic excitation of the standing accretion shock instability and an overall faster rate of shock expansion. This should be true no matter if the underlying SASI mechanism is acoustic (Blondin & Mezzacappa 2006) or advective-acoustic (Foglizzo et al. 2007) since both mechanisms depend on the acoustic flux from the PNS. The enhanced action of the SASI is reflected by a greater amount of transverse kinetic energy at large scales for softer EOS (Figure 12). For high-$L_{\nu _e}$ models that do not depend strongly on the SASI, we attribute the difference between STOS and LS to the higher mass accretion rate just after bounce for STOS.

In this study, we use a parameterized neutrino "light-bulb" with a fixed neutrino luminosity. There is no back-reaction from the EOS on the neutrino luminosity, allowing us to examine the dependence of the resulting explosion times on the EOS alone. Previous studies have shown that the neutrino luminosity can be significantly dependent on the EOS (Thompson et al. 2003; Sumiyoshi et al. 2005; Janka 2012), and this is a very important consequence of the choice of EOS, but our study shows that even for a fixed neutrino luminosity, the likelihood of a given model to explode also depends on the EOS.

The EOS for matter at temperatures and densities relevant to CCSNe is still an active area of research. As summarized in the Introduction and Table 1, available experiments and observations that constrain the EOS parameters favor LS220. LS180 is ruled out as it cannot produce a ∼2 M_☉ neutron star as required by the observations of PSR J1614-2230 (Demorest et al. 2010). STOS yields neutron stars of typical mass (around 1.4 M_☉) with too large a radius (∼15 km) to account for the results of Steiner et al. (2010), who estimate typical NS radii of 11–12 km. Several of the EOS parameters listed in Table 1 for STOS also fall outside of the experimental and theoretical estimates. So, LS220 is the EOS that satisfies, or nearly satisfies, the largest number of constraints. Our results show, additionally, that LS220 is intermediate between STOS and LS180 in favoring explosion for a given neutrino luminosity. Until we are more certain about the EOS at temperatures and densities relevant to CCSN simulations, LS220 seems a good choice. Some recent works have chosen to use LS220 on these same grounds (Ott et al. 2012; Dessart et al. 2012).

The author is grateful to Craig Wheeler, Christian Ott, Miloš Milosavljević, and Dongwook Lee for helpful discussions and to the anonymous referee for comments that improved this work. The author thanks Evan O'Connor for providing help with implementing the EOS routines used in this work, for providing EOS data, and for very insightful discussions. Support for this work was provided by NASA through Hubble Fellowship grant No. 51286.01 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. The software used in this work was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago. This research used computational resources at ALCF at ANL, which is supported by the Office of Science of the US Department of Energy under contract No. DE-AC02-06CH11357. The author acknowledges the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing high-performance computing, visualization, and data storage resources that have contributed to the research results reported within this paper.