IMPACTS OF COLLECTIVE NEUTRINO OSCILLATIONS ON CORE-COLLAPSE SUPERNOVA EXPLOSIONS

Yudai Suwa; Kei Kotake; Tomoya Takiwaki; Matthias Liebendörfer; Katsuhiko Sato

doi:10.1088/0004-637X/738/2/165

1. INTRODUCTION

Although the explosion mechanism of core-collapse supernovae is not yet completely understood, current multidimensional (multi-D) simulations based on refined numerical models show several promising scenarios. Among the candidates are the neutrino heating mechanism aided by convection and standing accretion shock instability (SASI; e.g., Marek & Janka 2009; Bruenn et al. 2009; Suwa et al. 2010), the acoustic mechanism (Burrows et al. 2007b), and the magnetohydrodynamic mechanism (e.g., Kotake et al. 2004, 2006; Obergaulinger et al. 2006; Burrows et al. 2007a; Takiwaki et al. 2009). Probably the best-studied one is the neutrino heating mechanism, whose basic concept was first proposed by Colgate & White (1966) and later reinforced by Bethe & Wilson (1985) to take the currently prevailing delayed form.

An important lesson from the multi-D simulations mentioned above is that hydrodynamic motions associated with convective overturn (Herant et al. 1994; Burrows et al. 1995; Janka & Mueller 1996; Fryer & Warren 2002, 2004) as well as the SASI (e.g., Blondin et al. 2003; Scheck et al. 2006; Ohnishi et al. 2006; Foglizzo et al. 2007; Murphy & Burrows 2008; Iwakami et al. 2008; Guilet et al. 2010; Fernández 2010) can help the onset of the neutrino-driven explosion, which otherwise generally fails in spherically symmetric one-dimensional (1D) simulations (Liebendörfer et al. 2001; Rampp & Janka 2002; Thompson et al. 2003; Sumiyoshi et al. 2005). This is mainly because the accretion timescale of matter in the gain region can be longer than in the 1D case, thereby enhancing the strength of neutrino–matter coupling there.

Neutrino-driven explosions have been obtained in the following state-of-the-art two-dimensional (2D) simulations. Using the MuDBaTH code, which includes one of the best available neutrino transfer approximations, Buras et al. (2006) first reported explosions for a non-rotating low-mass (11.2 M_☉) progenitor from Woosley et al. (2002), and then for a 15 M_☉ progenitor from Woosley & Weaver (1995) with a moderately rapid rotation imposed (Marek & Janka 2009). By implementing a multi-group flux-limited diffusion algorithm in the CHIMERA code (e.g., Bruenn et al. 2009), Yakunin et al. (2010) obtained explosions for non-rotating 12 M_☉ and 25 M_☉ progenitors from Woosley et al. (2002). More recently, Suwa et al. (2010) pointed out that a stronger explosion is obtained for a rapidly rotating 13 M_☉ progenitor from Nomoto & Hashimoto (1988) compared to the corresponding non-rotating model, when the isotropic diffusion source approximation (IDSA) for spectral neutrino transport (Liebendörfer et al. 2009) is implemented in the ZEUS code.

However, this success further begs new questions. First of all, the explosion energies obtained in these simulations are typically underpowered by one or two orders of magnitudes to explain the canonical supernova kinetic energy (∼10⁵¹ erg). Moreover, the softer nuclear equation of state (EOS) such as the Lattimer & Swesty (1991, hereafter LS) EOS with an incompressibility K = 180 MeV at nuclear densities is employed in those simulations. On top of evidence that favors a stiffer EOS based on nuclear experimental data (Shlomo et al. 2006), the soft EOS may not account for the recently observed massive neutron star of ∼2 M_☉ (Demorest et al. 2010; see the maximum mass for the LS180 EOS in O'Connor & Ott 2011). With a stiffer EOS, the explosion energy may be even lower than inferred from Marek & Janka (2009), who did not obtain the neutrino-driven explosion for their model with K = 263 MeV. What is then still missing? We may obtain the answer by using three-dimensional (3D) simulations (Nordhaus et al. 2010) or by taking into account new concepts, such as exotic physics in the core of the protoneutron star (Sagert et al. 2009), viscous heating by the magnetorotational instability (Thompson et al. 2005; Masada et al. 2011), or energy dissipation via Alfvén waves (Suzuki et al. 2008).

Joining in these efforts, we explore in this study the possible impacts of collective neutrino oscillations on energizing the neutrino-driven explosions. Collective neutrino oscillations, i.e., neutrinos of all energies that oscillate almost in phase, are attracting great attention, because they can induce dramatic observable effects such as a spectral split or swap (e.g., Raffelt & Smirnov 2007; Duan et al. 2008; Dasgupta et al. 2008; and references therein). These effects are predicted to emerge as distinct features in the energy spectra (see Duan et al. 2010; Dasgupta 2010; and references therein, for reviews of the rapidly growing research field). Among a number of important effects possibly created by self-interaction, we choose to consider the effect of spectral splits between electron- (ν_e) anti-electron neutrinos ( $\bar{\nu }_e$ ), and heavy-lepton neutrinos (ν_x, i.e., ν_μ, ν_τ, and their anti-particles) above a threshold energy (e.g., Fogli et al. 2007). Since ν_x have higher average energies than the other species in the postbounce phase, the neutrino flavor mixing would increase the effective energies of ν_e and $\bar{\nu }_e$ , and hence increase the neutrino heating rates in the gain region. A formalism to treat the neutrino oscillation using the Boltzmann neutrino transport is given in Yamada (2000) and Strack & Burrows (2005), but it is difficult to implement. To mimic the effects in this study, we perform the spectral swap by hand as a first step. By changing the average neutrino energy, $\langle\epsilon _{\nu _{x}}\rangle$ , as well as the position of the neutrino spheres ( $R_{\nu _x}$ ) in a parametric manner, we hope to constrain the parameter regions spanned by $\langle\epsilon _{\nu _{x}}\rangle$ and $R_{\nu _x}$ wherein the additional heating from collective neutrino oscillations could have impact on the explosion dynamics. Our strategy is as follows. We will first constrain the parameter regions to some extent by performing a number of 1D simulations. Here we also investigate the progenitor dependence using a suite of progenitor models (13, 15, 20, and 25 M_☉). After squeezing the condition in the 1D computations, we include the flavor conversions in 2D simulations to see their impact on the dynamics and also discuss how the critical condition for the collective effects in 1D can be subject to change in 2D.

The paper opens with descriptions of the initial models and the numerical methods, focusing on how to model the collective neutrino oscillations (Section 2). The main results are shown in Section 3. We summarize our results and discuss their implications in Section 4.

2. NUMERICAL METHODS

2.1. Hydrodynamics

The employed numerical methods are essentially the same as those in our previous paper (Suwa et al. 2010). For convenience, we briefly summarize them in the following. The basic evolution equations are written as

$\begin{equation} \frac{{d}\rho }{{d}t}+\rho \nabla \cdot \mathbf {v}=0, \end{equation} \tag{ 1 }$

$\begin{equation} \rho \frac{{d}\mathbf {v}}{{d}t}=-\nabla P -\rho \nabla \Phi, \end{equation} \tag{ 2 }$

$\begin{equation} \frac{{d}e^*}{{d}t}+ \nabla \cdot \left[\left(e^* + P\right) \mathbf {v}\right]= -\rho \mathbf {v} \cdot \nabla \Phi + Q_{E}, \end{equation} \tag{ 3 }$

$\begin{equation} \frac{{{d}Y_e}}{{d}t} = Q_N, \end{equation} \tag{ 4 }$

$\begin{equation} \bigtriangleup {\Phi } = 4\pi G \rho, \end{equation} \tag{ 5 }$

where $\rho, \mathbf {v}, P, \mathbf {v}, e^*, {\rm and}\, \Phi$ are density, fluid velocity, gas pressure including the radiation pressure of neutrinos, total energy density, and gravitational potential, respectively. The time derivatives are Lagrangian. As for the hydro solver, we employ the ZEUS-2D code (Stone & Norman 1992) which has been modified for core-collapse simulations (e.g., Suwa et al. 2007a, 2007b, 2009; Takiwaki et al. 2009). Q_E and Q_N (in Equations (3) and (4)) represent the change of energy and electron fraction (Y_e) due to interactions with neutrinos. To estimate these quantities, we implement spectral neutrino transport using the IDSA scheme (Liebendörfer et al. 2009). The IDSA scheme splits the neutrino distribution into two components, both of which are solved using separate numerical techniques. We apply the so-called ray-by-ray approach in which the neutrino transport is solved along a given radial direction assuming that the hydrodynamic medium for the direction is spherically symmetric. Although the current IDSA scheme does not yet include ν_x and the inelastic neutrino scattering with electrons, these simplifications save a significant amount of computational time compared to the canonical Boltzmann solvers (see Liebendörfer et al. 2009 for more details). Following the prescription in Müller et al. (2010), we improve the accuracy of the total energy conservation by using a conservation form in Equation (3), instead of solving the evolution of internal energy as originally designed in the ZEUS code. Numerical tests are presented in the Appendix.

The simulations are performed on a grid of 300 logarithmically spaced radial zones from the center up to 5000 km and 128 equidistant angular zones covering 0 ⩽ θ ⩽ π for 2D simulations. For the spectral transport, we use 20 logarithmically spaced energy bins ranging from 3 to 300 MeV.

2.2. Spectral Swapping

As mentioned in Section 1, we introduce a spectral interchange from heavy-lepton neutrinos (ν_μ, ν_τ, and their antineutrinos, collectively referred as ν_x hereafter) to electron-type neutrinos and antineutrinos, namely ν_x → ν_e and $\bar{\nu }_x \rightarrow \bar{\nu }_e$ . Instead of solving the transport equations for ν_x, we employ the so-called light-bulb approximation and focus on the optically thin region outside the neutrinosphere (e.g., Janka & Mueller 1996; Ohnishi et al. 2006).

According to Duan et al. (2010), the threshold energy, epsilon _th, is set to be 9 MeV, above which spectral swap takes place. Below the threshold, neutrino heating is estimated from the spectral transport via the IDSA scheme. Above the threshold, the heating rate is replaced by

$\begin{equation} Q_E\propto \int ^\infty _{\epsilon _{{\rm th}}} d\epsilon _\nu \epsilon ^3\left[j(\epsilon _\nu)+\chi (\epsilon _\nu)\right] f_{\nu }(r,\epsilon _\nu), \end{equation} \tag{ 6 }$

where j and χ are the neutrino emissivity and absorptivity, respectively, and f_ν(r, epsilon _ν) corresponds to the neutrino distribution function for ν_x with epsilon _ν being the energies of the electron neutrinos and antineutrinos. In the light-bulb approach, this is often approximated by the Fermi–Dirac distribution with a vanishing chemical potential (e.g., Ohnishi et al. 2006) as

$\begin{equation} f_{\nu }(r,\epsilon _\nu)=\frac{1}{e^{\epsilon _\nu /kT_{\nu _x}}+1}g(r), \end{equation} \tag{ 7 }$

where k and $T_{\nu _x}$ are the Boltzmann constant and the neutrino temperature, respectively. g(r) is the geometric factor, $g(r)=1-\left[1-(R_{\nu _x}/r)^2\right]^{1/2}$ , which is taken into account for the normalization, with $R_{\nu _x}$ the radius of the neutrinosphere. The neutrino luminosity of ν_x at infinity is given as

$\begin{equation} L_{\nu _x}=2.62\times 10^{52}\left(\frac{\langle\epsilon _{\nu _x}\rangle}{15 \,\mathrm{MeV}}\right)^4\left(\frac{R_{\nu _x}}{30 \,\mathrm{km}}\right)^2 \,\mathrm{erg \,s^{-1}}, \end{equation} \tag{ 8 }$

where $\langle\epsilon _{\nu _x}\rangle=\int ^\infty _0 d\epsilon _{\nu _x}\epsilon _{\nu _x}^3 f_\nu (\epsilon _{\nu _x})/\int ^\infty _0 d\epsilon _{\nu _x}\epsilon _{\nu _x}^2 f_\nu (\epsilon _{\nu _x})$ is the average energy of the emitted neutrinos. The position where spectral swapping sets in is fixed at 100 km (around the gain radius) and the onset time t_s is varied as a parameter: t_s = 100, 200, and 300 ms after bounce.

In fact, the threshold energy depends on the neutrino luminosities, spectra, and oscillation parameters (see, e.g., Duan et al. 2010 and references therein) with conserved net ν_e flux (i.e., the lepton number conservation). However, the conservation of the lepton number is too complicated to satisfy in the dynamical simulation because the neutrino spectrum and the luminosity evolve with time. In order to focus on the hydrodynamic features affected by the spectral modulation induced by the swapping, we simplify just a single threshold energy in this work.

To summarize, the parameters that we use to mimic spectral swapping are the following: (1) $R_{\nu _x}$ , the radius of the neutrinosphere of ν_x, (2) $\langle\epsilon _{\nu _x}\rangle$ , the average energy of ν_x, and (3) t_s, the time when spectral swapping sets in.

3. RESULT

3.1. One-dimensional Models

3.1.1. 1D without Spectral Swapping

In this subsection, we first outline the 1D collapse dynamics without spectral swapping. We take a 13 M_☉ progenitor (Nomoto & Hashimoto 1988) as a reference.

At around 112 ms after the onset of gravitational collapse, the bounce shock forms at a radius of ∼10 km with an enclosed mass of ∼0.7 M_☉.⁷ The central density at this time is ρ_c = 3.6 × 10¹⁴ g cm⁻³. The shock propagates outward but finally stalls at a radius of ∼100 km. Due to the decreasing accretion rate through the stalled shock, the shock can be still pushed outward. However, after some time, the shock radius begins to shrink. The ratio of the advection timescale, τ_adv, and the heating timescale, τ_heat, is an important indicator for the criteria of neutrino-driven explosion (Buras et al. 2006; Marek & Janka 2009; Suwa et al. 2010). In our 1D simulations, τ_adv/τ_heat is generally smaller than unity in the postbounce phase. This is the reason why our 1D simulations do not yield a delayed explosion. This is also the case for the other progenitors (15, 20, and 25 M_☉) investigated in this study. As for the accretion phase (later than ∼50 ms after the bounce), the typical neutrino luminosity at r = 5000 km is 3 × 10⁵² erg s⁻¹ for both ν_e and $\bar{\nu }_e$ , and the typical average energy is $\langle \epsilon _{\nu _e} \rangle\approx 9$ MeV and $\langle \epsilon _{\bar{\nu }_e}\rangle \approx 12$ MeV as shown in Figure 1. Figure 2 indicates the resultant neutrino luminosity spectrum 100 ms after the bounce.

**Figure 1.** Time evolution of the neutrino luminosity (top panel) and average energy (bottom panel) for ν_e (red solid line) and $\bar{\nu }_e$ (blue dashed line).
Download figure:
Standard image High-resolution image

**Figure 1.** Time evolution of the neutrino luminosity (top panel) and average energy (bottom panel) for ν_e (red solid line) and $\bar{\nu }_e$ (blue dashed line).
Download figure:
Standard image High-resolution image

**Figure 2.** Neutrino luminosity spectrum of ν_e (red solid line) and $\bar{\nu }_e$ (blue dashed line) without spectral swapping 100 ms after the bounce. For comparison, we show the injected luminosity spectrum of ν_x with $\langle\epsilon _{\nu _x}\rangle=15$ MeV, which will be swapped with the original spectrum of ν_e and $\bar{\nu }_e$ at _ν > 9 MeV for models including spectral swapping.
Download figure:
Standard image High-resolution image

3.1.2. 1D with Spectral Swapping

The investigated models with spectral swapping are summarized in Table 1. As already mentioned, the model parameters are the neutrinosphere radius ( $R_{\nu _x}$ ), the average energy of neutrinos ( $\langle \epsilon _{\nu _x}\rangle$ ), and the onset time of spectral swapping (t_s). The model names include these parameters: "NH13" represents the progenitor model, "R.." represents $R_{\nu _x}$ in units of km, "E.." represents $\langle \epsilon _{\nu _x}\rangle$ in MeV, "T.." represents t_s in ms, and "S" represents 1D (spherical symmetry).

Table 1. 1D Simulations

Model	Dimension	R_ν	$\langle\epsilon _{\nu _x}\rangle$	L_ν	t_s	Explosion	E^∞_diag	M^{t = t_s}₁₀	M^∞₁₀
		(km)	(MeV)	(10⁵² erg s⁻¹)	(ms)		(10⁵¹ erg)	(M_☉)	(M_☉)
NH13R10E15T100S	1D	10	15 MeV	0.29	100	No	...	1.18	...
NH13R10E17T100S	1D	10	17 MeV	0.48	100	No	...	1.18	...
NH13R10E18T100S	1D	10	18 MeV	0.60	100	No	...	1.18	...
NH13R10E19T100S	1D	10	19 MeV	0.75	100	Yes	1.00	1.18	1.14
NH13R10E20T100S	1D	10	20 MeV	0.92	100	Yes	1.49	1.18	1.12

NH13R20E13T100S	1D	20	13 MeV	0.66	100	No	...	1.18	...
NH13R20E13T150S	1D	20	13 MeV	0.66	150	No	...	1.21	...
NH13R20E13T200S	1D	20	13 MeV	0.66	200	No	...	1.25	...
NH13R20E14T100S	1D	20	14 MeV	0.88	100	No	...	1.18	...
NH13R20E14T150S	1D	20	14 MeV	0.88	150	No	...	1.21	...
NH13R20E14T200S	1D	20	14 MeV	0.88	200	No	...	1.25	...
NH13R20E15T100S	1D	20	15 MeV	1.16	100	Yes	0.97	1.18	1.15
NH13R20E15T150S	1D	20	15 MeV	1.16	150	Yes	0.54	1.21	<1.24
NH13R20E15T200S	1D	20	15 MeV	1.16	200	Yes	0.47	1.25	<1.26
NH13R20E21T100S	1D	20	21 MeV	4.47	100	Yes	5.56	1.18	1.07
NH13R20E22T100S	1D	20	22 MeV	5.39	100	Yes	6.50	1.18	1.07

NH13R28E13T100S	1D	28	13 MeV	1.29	100	No	...	1.18	...
NH13R29E13T100S	1D	29	13 MeV	1.38	100	No	...	1.18	...
NH13R30E11T100S	1D	30	11 MeV	0.76	100	No	...	1.18	...
NH13R30E11T150S	1D	30	11 MeV	0.76	150	No	...	1.21	...
NH13R30E11T200S	1D	30	11 MeV	0.76	200	No	...	1.25	...
NH13R30E12T100S	1D	30	12 MeV	1.07	100	No	...	1.18	...
NH13R30E12T150S	1D	30	12 MeV	1.07	150	No	...	1.21	...
NH13R30E12T200S	1D	30	12 MeV	1.07	200	No	...	1.25	...
NH13R30E13T100S	1D	30	13 MeV	1.48	100	Yes	0.85	1.18	<1.19
NH13R30E13T150S	1D	30	13 MeV	1.48	150	No	...	1.21	...
NH13R30E13T200S	1D	30	13 MeV	1.48	200	No	...	1.25	...
NH13R30E14T100S	1D	30	14 MeV	1.99	100	Yes	1.58	1.18	1.12
NH13R30E14T150S	1D	30	14 MeV	1.99	150	Yes	0.98	1.21	1.19
NH13R30E14T200S	1D	30	14 MeV	1.99	200	Yes	0.68	1.25	1.22
NH13R30E15T100S	1D	30	15 MeV	2.62	100	Yes	2.27	1.18	1.10
NH13R30E15T150S	1D	30	15 MeV	2.62	150	Yes	1.43	1.21	1.16
NH13R30E15T200S	1D	30	15 MeV	2.62	200	Yes	0.93	1.25	1.22
NH13R30E20T100S	1D	30	20 MeV	8.28	100	Yes	6.84	1.18	1.07

NH13R40E11T100S	1D	40	11 MeV	1.35	100	No	...	1.18	...
NH13R50E11T100S	1D	50	11 MeV	2.10	100	Yes	0.86	1.18	<1.18
NH13R60E11T100S	1D	60	11 MeV	3.03	100	Yes	1.48	1.18	1.12

Download table as: ASCII Typeset image

Figure 3 presents the time evolution of the mass shells for models NH13R30E12T100S and NH13R30E13T100S. The difference between these panels is the average energies of neutrinos, $\langle\epsilon _{\nu _x}\rangle=12$ MeV for the top panel and 13 MeV for the bottom panel. The thick solid lines represent the radial positions of the shock waves. Regardless of a small difference in $\langle\epsilon _{\nu _x}\rangle$ , model NH13R30E13T100S shows a shock expansion after the manual spectral swapping is switched on (see the thick line in the bottom panel of Figure 3), while the stalled shock does not revive for model NH13R30E12T100S (top panel). This suggests that there is a critical condition for successful explosion induced by spectral swapping. In the bottom panel, the regions enclosing the mass of M_r ∼ 1.2 M_☉ (thin black line) correspond to the so-called mass cut, which could be interpreted as the final mass of the remnant. The fact that a clear mass cut emerges in model NH13R30E13T100S indicates that a neutron star will be left behind in this model. Such a definite mass cut has been observed by Kitaura et al. (2006) who reported a successful neutrino-driven explosion (in 1D) for a lighter progenitor star, which is, however, difficult to realize for more massive stars in 2D (e.g., Figure 2 in Marek & Janka 2009 and Figure 1 in Suwa et al. 2010).

**Figure 3.** Time evolution of mass shells for NH13R30E12T100S (top) and NH13R30E13T100S (bottom). The black thin line corresponds to 1.2 M_☉ and the black thick line represents the shock wave position, respectively. The difference between these panels is the average energies of neutrinos, $\langle\epsilon _{\nu _x}\rangle=12$ MeV for the top panel and 13 MeV for the bottom panel.
Download figure:
Standard image High-resolution image

As a tool to measure the strength of an explosion, we define a diagnostic energy that refers to

$\begin{equation} E_\mathrm{diag}=\int _D dV \left(\frac{1}{2}\rho |v|^2+e-\rho \Phi \right), \end{equation} \tag{ 9 }$

where e is the internal energy, and D represents the domain in which the integrand is positive. Figure 4 shows the time evolution of E_diag for some selected models. The diagnostic energy increases with time for the green dotted line, but which decreases for the red line, noting that the difference between the pair of models is $\Delta \langle \epsilon _{\nu _x}\rangle = 1$ MeV. The blue dashed line (model NH13R30E15T100S) has $\langle\epsilon _{\nu _x}\rangle=15$ MeV and reaches larger E_diag than the green line (NH13R30E13T100S; $\langle\epsilon _{\nu _x}\rangle=13$ MeV). On the other hand, the later injection of spectral swapping leads to smaller E_diag, i.e., the brown dot-dashed line (t_s = 200 ms) shows smaller E_diag than the blue dashed line (t_s = 100 ms). For models that experience earlier spectral swapping with higher neutrino energy, the diagnostic energy becomes higher in an earlier stage, as expected.

**Figure 4.** Diagnostic energies as functions of time. Red solid, green dotted, blue dashed, and brown dot-dashed lines correspond to models NH13R30E12T100S, NH13R30E13T100S, NH13R30E15T100S, and NH13R30E15T200S, respectively. The ν_x average energies of three lines other than the brown dot-dashed line differ from each other, i.e., $\langle\epsilon _{\nu _x}\rangle=$ 12 MeV (red solid), 13 MeV (green dotted), and 15 MeV (blue dashed), respectively. As for the brown dot-dashed line, only the onset time of spectral swapping is different from the blue dashed line. The red solid line shows only an oscillation, while the other lines show increasing diagnostic energy.
Download figure:
Standard image High-resolution image

Looking at Figure 4 again, E_diag for the exploding models seems to saturate with time. These curves can be fitted by the following function,

$\begin{equation} E_\mathrm{diag}(t)=E_\mathrm{diag}^\infty (1-\mathrm{e}^{-at+b}), \end{equation} \tag{ 10 }$

where E^∞_diag is a converging value of E_diag, and a and b are the fitting parameters. As for NH13R30E13T100S, E^∞_diag = 8.5 × 10⁵⁰ erg. This fitting formula allows us to estimate the final diagnostic energy especially for the strongly exploding models whose diagnostic energy we cannot estimate in principle because the shock goes beyond the computational domains (r < 5000 km) before saturation.

Figure 5 shows the summary of 1D models for a given neutrino luminosity determined by $R_{\nu _x}$ and $\langle\epsilon _{\nu _x}\rangle$ (Equation (8)). The gray lines correspond to the neutrino luminosities determined by the pairs of $R_{\nu _x}$ and $\langle\epsilon _{\nu _x}\rangle$ which is (1–5) × 10⁵² erg s⁻¹ from bottom to top. Circles and crosses correspond to the exploding and non-exploding models, respectively. Not surprisingly, explosions are more easily obtained for higher neutrino luminosity.

As is well known, the combination of $\langle \epsilon _{\nu _x} \rangle$ and $L_{\nu _x}$ is an important quantity in diagnosing the success or failure of explosions, because the neutrino heating rate in the so-called gain region, Q⁺_ν, is proportional to $\langle\epsilon _{\nu _x}^2\rangle L_{\nu _x}$ (e.g., Equation (23) in Janka 2001).

Figure 6 shows E^∞_diag as a function of $\langle\epsilon _{\nu _x}\rangle^2L_{\nu _x}$ . Note in the plot that we set the horizontal axis not as $\langle\epsilon _{\nu _x}^2\rangle L_{\nu _x}$ but as $\langle\epsilon _{\nu _x}\rangle^2L_{\nu _x}$ so that we can deduce the following dependence more clearly and easily.⁸ In this figure, let us first focus on the red pluses, green crosses, and blue squares whose difference is characterized by t_s (2D results (filled circles) will be mentioned in a later section). The red (t_s = 100 ms), green (t_s = 150 ms), and blue (t_s = 200 ms) points have a clear correlation with $\langle\epsilon _{\nu _x}\rangle^2L_{\nu _x}$ . The orange and light-blue regions represent the non-exploding regions for the red and blue points, respectively. Both of them show that the minimum E^∞_diag decreases with t_s, indicating that the critical values of 〈 epsilon _ν〉²L_ν for explosion sharply depend on t_s. This is because the mass outside the shock wave gets smaller with time so that the minimum energy to blow up the star gets smaller, too. By the same reason, E_diag becomes larger as t_s becomes smaller given the same $\langle\epsilon _{\nu _x}\rangle^2L_{\nu _x}$ . To obtain a larger E^∞_diag, earlier spectral swapping is preferred.

Figure 7 shows the neutrino heating rate and the density distribution of NH13R30E13T100S for 10 ms and 250 ms after t_s (= 100 ms after the bounce). As the shock wave propagates outward, the density in the gain region drops sharply (e.g., 100–200 km, dashed blue line), leading to the suppression of the heating rate (dashed red line). This is the reason for the saturation in E_diag as shown in Figure 4.

The remnant mass is an important indicator to diagnose the consequences of the explosion in producing either a neutron star or a black hole. The last two lines in Table 1 show the integrated masses in the regions of ρ ⩾ 10¹⁰ g cm⁻³ at t = t_s and t = ∞. The latter one is estimated by fitting as

$\begin{equation} M_{10}(t)=M_{10}^\infty (1+\mathrm{e}^{-ct+d}),\vspace*{-5pt} \end{equation} \tag{ 11 }$

where c and d are the fitting parameters. For the exploding models, M^∞₁₀ becomes generally smaller than M^{t = t_s}₁₀ because of the mass ejection. Exceptions are weakly exploding models (NH13R20E15T150S, NH13R20E15T200S, NH13R30E13T100S, and NH13R50E11T100S), in which the mass accretion continues after t_s and eventually stops at late time (maximum masses are presented in Table 1). For the non-exploding models, the remnant mass simply increases with time. Regarding the 13 M_☉ progenitor investigated in this section, the remnant masses in models that produce strong explosion (E^∞_diag ≳ 10⁵¹ erg) are considerably smaller (1.1–1.2 M_☉) if compared to the typical mass of observed neutron stars ∼1.4 M_☉ (Lattimer & Prakash 2007). This may simply reflect the light iron core (∼1.26 M_☉) inherent to the progenitor or the existence of mass accretion induced by the matter fallback after the explosion. We now move on to investigate the progenitor dependence in the next section.

3.1.3. The Progenitor Dependence

In addition to the 13 M_☉ progenitor by Nomoto & Hashimoto (1988), we are going to investigate the progenitor dependence in 1D simulations. The computed models are NH15 (15 M_☉; Nomoto & Hashimoto 1988), s15s7b2 (15 M_☉; Woosley & Weaver 1995), s15.0 (15 M_☉), s20.0 (20 M_☉), and s25.0 (25 M_☉; Woosley et al. 2002), which are all listed in Table 2. The first set of characters for these models indicate the progenitors:

1.
NH: (Nomoto & Hashimoto 1988),
2.
WW: (Woosley & Weaver 1995), and
3.
WHW: (Woosley et al. 2002).

Table 2. Progenitor Dependence

Model	Dimension	R_ν	$\langle\epsilon _{\nu _x}\rangle$	L_ν	t_s	Explosion	E^∞_diag	M^{t = t_s}₁₀	M^∞₁₀
		(km)	(MeV)	(10⁵² erg s⁻¹)	(ms)		(10⁵¹ erg)	(M_☉)	(M_☉)
NH15R30E11T100S	1D	30	11 MeV	0.76	100	No	...	1.34	...
NH15R30E12T100S	1D	30	12 MeV	1.07	100	No	...	1.34	...
NH15R30E13T100S	1D	30	13 MeV	1.48	100	Yes	0.65	1.34	<1.38
NH15R30E14T100S	1D	30	14 MeV	1.99	100	Yes	2.17	1.34	1.25
NH15R30E15T100S	1D	30	15 MeV	2.62	100	Yes	3.73	1.34	1.21

WW15R30E11T100S	1D	30	11 MeV	0.76	100	No	...	1.40	...
WW15R30E12T100S	1D	30	12 MeV	1.07	100	No	...	1.40	...
WW15R30E13T100S	1D	30	13 MeV	1.48	100	No	...	1.40	...
WW15R30E14T100S	1D	30	14 MeV	1.99	100	Yes	1.94	1.40	1.31
WW15R30E15T100S	1D	30	15 MeV	2.62	100	Yes	3.41	1.40	1.25

WHW15R30E11T100S	1D	30	11 MeV	0.76	100	No	...	1.49	...
WHW15R30E12T100S	1D	30	12 MeV	1.07	100	No	...	1.49	...
WHW15R30E13T100S	1D	30	13 MeV	1.48	100	No	...	1.49	...
WHW15R30E14T100S	1D	30	14 MeV	1.99	100	No	...	1.49	...
WHW15R30E15T100S	1D	30	15 MeV	2.62	100	Yes	3.55	1.49	1.36

WHW20R30E11T100S	1D	30	11 MeV	0.76	100	No	...	1.45	...
WHW20R30E12T100S	1D	30	12 MeV	1.07	100	No	...	1.45	...
WHW20R30E13T100S	1D	30	13 MeV	1.48	100	Yes	0.99	1.45	...
WHW20R30E14T100S	1D	30	14 MeV	1.99	100	Yes	2.20	1.45	1.34
WHW20R30E15T100S	1D	30	15 MeV	2.62	100	Yes	3.61	1.45	1.29

WHW25R30E12T100S	1D	30	12 MeV	1.07	100	No	...	1.69	...
WHW25R30E13T100S	1D	30	13 MeV	1.48	100	No	...	1.69	...
WHW25R30E14T100S	1D	30	14 MeV	1.99	100	No	...	1.69	...
WHW25R30E15T100S	1D	30	15 MeV	2.62	100	Yes	0.73	1.69	<2.00
WHW25R30E16T100S	1D	30	16 MeV	3.39	100	Yes	5.92	1.69	1.49
WHW25R30E17T100S	1D	30	17 MeV	4.32	100	Yes	9.21	1.69	1.41

Download table as: ASCII Typeset image

Figure 8 depicts density profiles of these progenitors 100 ms after the bounce as a function of the enclosed mass (M_r). It can be seen that the density profiles for M_r ≲ 0.8 M_☉ are almost insensitive to the progenitor masses despite the difference in the pre-collapse phase (see, e.g., Figure 1 of Burrows et al. 2007b). On the other hand, the profiles of M_r ≳ 0.8 M_☉ differ between progenitors so that the critical heating rates and E^∞_diag are also expected to be different. In Figure 8, the envelope of WHW25 is shown to be the thickest, while the envelope of NH13 is the thinnest.

Figure 9 shows the critical heating rates as a function of the progenitor masses. In agreement with our intuition, the critical heating rates for models WHW25 and NH13 are in the high and low ends, respectively. However, the critical heating rate for model WHW20 is almost the same as the one for model NH13 although the envelope of model WHW20 is much thicker than that of model NH13 (see Figure 8). Our results show that the critical heating rate is indeed affected by the envelope mass; however, the relation is not one-to-one. It is also interesting to note that the critical heating rates for the 15 M_☉ progenitors of WW15, WHW15, and NH15 are different by a factor of ∼3, which may send a clear message that an accurate knowledge of supernova progenitors is also pivotal to pin down the supernova mechanism.

**Figure 9.** Critical heating rate, $\langle\epsilon _{\nu _x}\rangle^2L_{\nu _x}$ , as a function of the progenitor mass, M. The circles are progenitors from Nomoto & Hashimoto (1988), the square is from Woosley & Weaver (1995), and the crosses are from Woosley et al. (2002), respectively. The error bars represent the distance between the last failing and the first exploding model in our grid of models. The symbols are located at the centers of error bars. The error bar is small for model NH13 because we calculated a more refined grid of models for the 13 M_☉ progenitor (Table 1) than for the 15–25 M_☉ progenitors (Table 2).
Download figure:
Standard image High-resolution image

The integrated masses with ρ ⩾ 10¹⁰ g cm⁻³ for t = t_s and t = ∞ are listed in the last two lines in Table 2 and Figure 10. The tendencies are the same as those found with NH13. As for model WHW25, we obtain results with E^∞_diag > 10⁵¹ erg and M^∞₁₀ > 1.4 M_☉, simultaneously.

3.2. Two-dimensional Models

Here we discuss the effects of spectral swapping in 2D (axisymmetric) simulations. Since our 2D simulations, albeit utilizing the IDSA scheme, are still computationally expensive, it is not practical to perform a systematic survey in 2D as we have done in 1D simulations. Looking at Figure 9 again, we choose models WHW15 (Woosley et al. 2002) and NH13 (Nomoto & Hashimoto 1988), whose critical heating rates are in the high and low ends, respectively.

3.2.1. 2D without Spectral Swapping

The basic hydrodynamic picture is the same with 1D before the shock stall (e.g., until ≲ 10 ms after bounce). After that, convection as well as SASI sets in between the stalled shock and the gain radius, which leads to the neutrino-heated shock revival for model NH13 (e.g., Suwa et al. 2010). For model WHW15, the position of the stalled shock, following several oscillations begins to shrink at ≳ 400 ms after bounce.

Even after the shock revival, it should be emphasized that the shock propagation for model NH13 is the so-called passive one (Buras et al. 2006). This means that the amount of mass ejection is smaller than the accretion in the post-shock region of the expanding shock (see the motions of mass shells in the post-shock region of Figure 1 in Suwa et al. 2010). Some regions have a positive local energy (Equation (9)), but the volume integrated value is quite small, ≲ 10⁵⁰ erg at the maximum. In order to reverse the passive shock into an active one it is most important to energize the explosion in some way. Using these two progenitors that produce a very weak explosion (model NH13) and do not show even a shock revival (model WHW15), we hope to explore how the dynamics would change when spectral swapping is switched on.

3.2.2. 2D with Spectral Swapping

Table 3 shows a summary for our 2D models, in which the last character of each model (A) indicates "Axisymmetry." Models NH13A and WHW15A are 2D models without spectral swapping for NH13 and WHW15, respectively.

Table 3. 2D Simulations

Model	Dimension	R_ν	$\langle\epsilon _{\nu _x}\rangle$	L_ν	t_s	Explosion	E^∞_diag	M^{t = t_s}₁₀	M^∞₁₀
		(km)	(MeV)	(10⁵² erg s⁻¹)	(ms)		(10⁵¹ erg)	(M_☉)	(M_☉)
NH13A	2D	...	...	...	...	Yes	∼0.1 (oscillating)	...	...
NH13R30E11T100A	2D	30	11 MeV	0.76	100	No	...	1.18	...
NH13R30E12T100A	2D	30	12 MeV	1.07	100	Yes	0.45	1.18	<1.23
NH13R30E13T100A	2D	30	13 MeV	1.48	100	Yes	1.03	1.18	<1.18
NH13R30E15T100A	2D	30	15 MeV	2.62	100	Yes	2.33	1.18	1.10

WHW15A	2D	...	...	...	...	No	...	...	...
WHW15R30E13T100A	2D	30	13 MeV	1.48	100	No	...	...	...
WHW15R30E14T100A	2D	30	14 MeV	1.99	100	Yes	1.96	1.48	<1.52
WHW15R30E15T100A	2D	30	15 MeV	2.62	100	Yes	3.79	1.48	1.34

Download table as: ASCII Typeset image

As in 1D, the onset of spectral swapping is taken to be t_s = 100 ms after bounce. At this time, model NH13 shows the onset of gradual shock expansion with a small diagnostic energy of E_diag ∼ 3 × 10⁴⁹ erg; the shock radius is located at ∼300 km. As for model WHW15, there is no region with a positive local energy (e.g., Equation (9)) and the shock radius is ∼200 km. The density profile for this model is essentially the same as the one in the 1D counterpart (see Figure 8) but with small angular density modulations due to convection.

In Figure 6, red filled circles represent E^∞_diag for model NH13. It can be seen that the critical heating rate to obtain E^∞_diag ∼ 10⁵¹ erg is smaller for 2D than the corresponding 1D counterparts (compare the heating rates for $\langle\epsilon _{\nu _x}\rangle^2 L_{\nu _x} \sim 2.2\times 10^{54} \,\mathrm{MeV^2 \,erg \,s^{-1}}$ ). In fact, models with $\langle\epsilon _{\nu _x}\rangle^2 L_{\nu _x}\lesssim 2.2\times 10^{54} \,\mathrm{MeV^2 \,erg \,s^{-1}}$ fail to explode in 1D, but succeed in 2D (albeit with a relatively small E^∞_diag: less than 10⁵¹ erg). As opposed to 1D, it is rather difficult in 2D to determine a critical heating rate due to the stochastic nature of the explosion triggered by SASI and convection. In our limited set of 2D models, the critical heating rate is expected to be close to $\langle\epsilon _{\nu _x}\rangle^2 L_{\nu _x}\sim 1.5\times 10^{54} \,\mathrm{MeV^2 \,erg \,s^{-1}}$ , below which the shock does not revive (e.g., $\langle\epsilon _{\nu _x}\rangle^2 L_{\nu _x} \lesssim 1\times 10^{54} \,\mathrm{MeV^2 \,erg \,s^{-1}}$ is the lowest end in the horizontal axis in the figure).

As seen from Figure 6, E^∞_diag becomes visibly larger for 2D than 1D especially for a smaller $\langle\epsilon _{\nu _x}\rangle^2 L_{\nu _x}$ . As the heating rates become larger, the difference between 1D and 2D becomes smaller because the shock revival occurs almost in a spherically symmetric way (before SASI and convection develop nonlinearly). In Table 3, it is interesting to note that model NH13R30E11T100A fails to explode, while we observe shock revival for the corresponding model without spectral swapping (model NH13A). This is because the heating rate of model NH13R30E11T100A is smaller than that of NH13A due to the small $\langle\epsilon _{\nu _x}\rangle$ , which can make it more difficult to trigger ν_x explosions. On the other hand, if the energy gain due to the swap is high enough (i.e., for models with greater $\langle\epsilon _{\nu _x}\rangle$ than E12 in Table 3), the swap can facilitate explosions.

Figure 11 depicts the entropy distributions for models NH13A (top panel) and NH13R30E13T100A (bottom panel). It can be seen that model NH13A shows a unipolar-like explosion (see also Suwa et al. 2010), while model NH13R30E13A explodes in a rather spherical manner as mentioned above. Model NH13A experiences several oscillations aided by SASI and convection before explosion, while the stalled shock for model NH13R30E13T100A turns into expansion shortly after the onset of spectral swapping. In fact, the shapes of hot bubbles behind the expanding shock are shown to be barely changing with time (bottom panel), which indicates a quasi-homologous expansion of material behind the revived shock.

Figure 12 shows the time evolution of mass shells for models NH13A (thin-gray lines) and NH13R30E13T100A (thin-orange lines). Black and red thick lines represent the shock position at the north pole for each model. The mass shells for model NH13A continue to accrete to the PNS, since the shock passively expands outward as previously mentioned. Due to this continuing mass accretion, the remnant for this model would be a black hole instead of a neutron star. On the other hand, model NH13R30E13T100A shows a mass ejection with a definite outgoing momentum in the post-shock region so that the remnant could be a neutron star. Unfortunately, however, we cannot predict the final outcome due to the limited simulation time. A long-term simulation recently done in 1D (e.g., Fischer et al. 2010) should be indispensable also for our 2D case. This is, however, beyond the scope of this paper.

Here let us discuss the validity of the parameters for the spectral swap that we have assumed so far. For example, the criteria of explosion for model NH13R30E12T100A were $L_{\nu _x}\approx 1.07\times 10^{51}$ erg s⁻¹ and $\langle\epsilon _{\nu _x}\rangle\approx 12$ MeV. These values are even smaller than the typical values obtained in 1D Boltzmann simulations (e.g., Liebendörfer et al. 2004), which show $L_{\nu _x}\approx 2\times 10^{52}$ erg s⁻¹ and $\sqrt{\vphantom{A^A}\smash{\hbox{${\langle\epsilon _{\nu _x}^2\rangle}$}}}\approx 20$ MeV (i.e., $\langle\epsilon _{\nu _x}\rangle\approx 14$ MeV with a vanishing chemical potential) earlier in the postbounce phase. Therefore, spectral swapping, if it would work as we have assumed, may potentially assist explosions.

It should be noted that the critical heating rate in this study might be too small due to the approximation of the light-bulb scheme. In this scheme, we can include the geometrical effect of the finite size of the neutrinosphere as in Equation (7), but cannot include the back reaction by the matter, i.e., the absorption of neutrino. Some fraction of neutrinos, in fact, are absorbed in the gain region and the neutrino luminosity decreases with the radius. We omit this effect in this study so that the heating rate might be overestimated in the simulation with spectral swapping. Thus, a fully consistent simulation including spectral swapping is necessary for a more realistic critical heating rate, which is beyond the scope of this study.

Finally, we discuss the 15 M_☉ progenitor labeled WHW15. As mentioned, this progenitor fails to explode without spectral swapping even in 2D.⁹ Figure 13 shows the entropy distributions of WHW15A (left; non-exploding) and WHW15R30E15T100A (right; exploding) for 220 ms after the bounce (corresponding to 120 ms after t_s for model WHW15R30E15T100A). The model with $R_{\nu _x}=30$ km and $\langle\epsilon _{\nu _x}\rangle=14$ MeV does not explode in 1D but explodes in 2D (compare Tables 2 and 3). Again, the multidimensionality helps the onset of explosion. The critical heating rate in 2D is in the range of $2.5\le \langle\epsilon _{\nu _x}\rangle^2L_{\nu _x}/(10^{54}$ MeV² erg s⁻¹) ⩽ 3.9, while it is $3.9\le \langle\epsilon _{\nu _x}\rangle^2 L_{\nu _x}/(10^{54}$ MeV² erg s⁻¹) ⩽ 5.9 in 1D. Therefore, the critical heating rate in 2D can be by a factor of ∼2 smaller than in 1D. In 2D, the critical ν_x luminosity and average energy to obtain explosion are $L_{\nu _x}\sim 2\times 10^{52}$ erg s⁻¹ and $\langle\epsilon _{\nu _x}\rangle\sim 14$ MeV (corresponding to $\sqrt{\vphantom{A^A}\smash{\hbox{${\langle\epsilon _{\nu _x}^2\rangle}$}}}\sim 20$ MeV), which are close to the results obtained in a 1D Boltzmann simulation (Sumiyoshi et al. 2005) for a 15 M_☉ progenitor.¹⁰ The diagnostic energy as well as the estimated remnant masses are listed in the last three columns in Table 3. E^∞_diag (as well as M^∞_diag) for exploding models is shown to be larger than the model series of NH13. As a result, some of the 2D models for WHW15 produce strong explosions (E^∞_diag ∼ 10⁵¹ erg) while simultaneously leaving behind a remnant of 1.34–1.52 M_☉. We think that it is only a solution accidentally found by our parametric explosion models. However, again, the critical heating rates that require the neutrino-driven explosion to be assisted by spectral swapping are never far away from the ones obtained in the Boltzmann simulations. We hope that our exploratory results may give momentum to supernova theorists to elucidate the effects of collective neutrino oscillations in a more consistent manner.

**Figure 13.** Entropy distributions of WHW15A (left) and WHW15R30E15T100A (right) for 220 ms after the bounce.
Download figure:
Standard image High-resolution image

4. SUMMARY AND DISCUSSION

We performed a series of 1D and 2D hydrodynamic simulations of core-collapse supernovae with spectral neutrino transport via the IDSA scheme. To model spectral swapping, one of the possible outcomes of collective neutrino oscillations, we parameterized the onset time when the spectral swap begins, the radius where the spectral swap takes place, and the threshold energy above which the spectral interchange between heavy-lepton neutrinos and electron/anti-electron neutrinos occurs. By doing so, we systematically studied the shock evolution and matter ejection due to neutrino heating enhanced by spectral swapping. We also investigated the progenitor dependence using a suite of progenitor models (13, 15, 20, and 25 M_☉). With these computations, we found that there is a critical heating rate induced by spectral swapping to trigger explosions that differs between the progenitors. The critical heating rate is generally smaller for 2D than 1D due to the multidimensionality that enhances the neutrino heating efficiency (see also Janka & Mueller 1996). The remnant masses which range from 1.1–1.5 M_☉ depending on the progenitors, can be determined by the mass ejection driven by the neutrino heating. For our 2D model of the 15 M_☉ progenitor, we found a set of parameters that produces an explosion with a canonical supernova energy close to 10⁵¹ erg and at the same time leaves behind a remnant mass close to ∼1.4 M_☉. Our results suggest that collective neutrino oscillations have the potential to solve the supernova problem if they occur. These effects should be explored in a more self-consistent manner in hydrodynamic simulations.

Here it should be noted that the simulations in this paper are only a first step toward more realistic supernova modeling. For the neutrino transfer, we omitted the cooling of heavy-lepton neutrinos and the inelastic neutrino scattering by electrons. These omissions lead to an overestimation of the diagnostic energy and they also should relax the criteria for explosion. The ray-by-ray approximation may lead to an overestimation of the directional dependence of the neutrino anisotropies. A full-angle transport will give us a more correct answer (see Ott et al. 2008; Brandt et al. 2011). Moreover, due to the coordinate symmetry axis, the SASI develops preferentially along the axis; it could thus provide a more favorable condition for the explosion. As several exploratory simulations have been done recently (e.g., Iwakami et al. 2008; Scheidegger et al. 2008; Nordhaus et al. 2010), 3D supernova models are indeed necessary also to pin down the outcomes of spectral swapping.

Finally, we briefly discuss whether the oscillation parameters in this paper are really valid in view of recent work that focus on clarifying the still-veiled nature of collective neutrino oscillations. Following Duan et al. (2010), there are at least two conditions for the onset of collective neutrino oscillations in the case of inverted neutrino mass hierarchy.

The first criterion should be satisfied in the so-called bipolar regime of the collective oscillation. In this regime, the neutrino number density should exceed the critical value,

$\begin{eqnarray} n_{\bar{\nu }_e, \rm crit} &\simeq &\frac{1}{(\sqrt{1+\chi }-1)^2}\frac{\Delta m^2}{\sqrt{2}G_\mathrm{F} \langle\epsilon _{\bar{\nu }_e}\rangle}\nonumber \\ &\simeq & 1.4\times 10^{29} \,\mathrm{cm^{-3}}\left(\frac{0.2}{\chi }\right)^2\left(\frac{15 \,\mathrm{MeV}}{\langle\epsilon _{\bar{\nu }_e}\rangle}\right),\vspace*{-20pt} \end{eqnarray} \tag{ 12 }$

where χ is the fractional excess of neutrinos over antineutrinos, Δm² is the characteristic mass-squared splitting (a typical value of ∼2.4 × 10⁻³ eV² is employed here), and G_F is the Fermi coupling constant. By using our simulation results, we can estimate χ, which is often treated as a parameter (typically ∼0.01–0.25). The following estimation is given in Esteban-Pretel et al. (2007): $\chi \simeq F_{\nu _e}/F_{\bar{\nu }_e}-1$ in the case of vanishing $F_{\nu _x}$ , where $F_{\nu _i}$ is the number flux of ν_i. From Figure 14, it can be seen that χ ∼ 0.2–0.3 for 100–400 ms after bounce. Since the typical number density in the post-shock region (r ∼200–300 km) can be estimated as

$\begin{eqnarray} n_{\bar{\nu }_e}&=&\frac{L_{\bar{\nu }_e}}{4\pi r^2 c\langle\epsilon _{\bar{\nu }_e}\rangle} \sim 1.1\times 10^{31} \,\mathrm{cm}^{-3} \left(\frac{L_{\bar{\nu }_e}}{10^{52} \,\mathrm{erg \,s^{-1}}}\right)\nonumber\\[-2pt] && \times \left(\frac{100 \,\mathrm{km}}{r}\right)^2 \left(\frac{15 \,\mathrm{MeV}}{\langle\epsilon _{\bar{\nu }_e}\rangle}\right), \vspace*{-18pt} \end{eqnarray} \tag{ 13 }$

the first condition is satisfied.¹¹

**Figure 14.** Time evolution of $\chi =F_{\nu _e}/F_{\bar{\nu }_e}-1$ , where $F_{\nu _i}$ is the number flux of ν_i.
Download figure:
Standard image High-resolution image

**Figure 14.** Time evolution of $\chi =F_{\nu _e}/F_{\bar{\nu }_e}-1$ , where $F_{\nu _i}$ is the number flux of ν_i.
Download figure:
Standard image High-resolution image

The second criterion is related to the decoherence of collective oscillations by matter. In order to overwhelm the suppression by decoherence, the following condition should be satisfied:

$\begin{equation} n_{\bar{\nu }_e,{\rm crit}} \sim n_e, \end{equation} \tag{ 14 }$

where n_e is the number density of electrons where decoherence takes place. This is equivalent to

$\begin{equation} Y_{\bar{\nu }_e,{\rm crit}} \sim Y_e. \end{equation} \tag{ 15 }$

In our 1D simulation, $Y_{\bar{\nu }_e}\sim (0.1\hbox{--}0.2)\times Y_e$ for 100 km ≲ r ≲ r_sh, where r_sh is the shock radius.¹² Since this condition is barely satisfied, the collective oscillations in reality could modify the spectrum to some extent between heavy-lepton neutrinos and electron/anti-electron neutrinos; however, the full swapping assumed in this study may be exaggerated. Very recently,¹³Chakraborty et al. (2011a, 2011b) pointed out that the matter effect could fully suppress spectral swapping in the accretion phase using the 1D neutrino-radiation hydrodynamic simulation data of Fischer et al. (2010). However, the current understanding of the collective oscillation is not completed and calculations in this field employ several assumptions (e.g., single angle approximation; but see also Dasgupta et al. 2011 for more recent work). To draw a robust conclusion, one needs a more detailed study that includes the collective neutrino flavor oscillation in the hydrodynamic simulations in a more self-consistent manner, which we are going to undertake as a sequel of this study.

We thank K. Sumiyoshi, H. Suzuki, S. Yamada, and T. Yoshida for stimulating discussions. Numerical computations were in part carried on XT4 at CfCA of the National Astronomical Observatory of Japan. M.L. is supported by the Swiss National Science Foundation under grant Nos. PP00P2-124879 and 200020-122287. This study was supported in part by the Japan Society for Promotion of Science (JSPS) Research Fellowships (Y.S.), the Grants-in-Aid for the Scientific Research from the Ministry of Education, Science and Culture of Japan (Nos. 19104006, 19540309, and 20740150), and HPCI Strategic Program of Japanese MEXT.

APPENDIX: CODE VALIDITY

A1. Conservation of Energy

In this section, we demonstrate the conservation of physical quantities using the spherical collapse model (NH13). Figure 15 depicts the evolution of total binding energy of gravity (red line), total internal energy (green), total kinetic energy (blue), total trapped-neutrino energy (magenta), total energy leaked by neutrinos (cyan), and variation of overall energy (black dashed), respectively. The gravitational energy and total energy are negative and the absolute values are shown. The gravitational energy and internal energy dominate (with different sign) and reach ∼10⁵³ erg soon after bounce. Despite such an enormous energy change, the total energy varies only within ∼3 × 10⁴⁹ erg so that the violation of energy conservation remains < 0.03%. The energy of the trapped neutrinos decreases with the diffusion timescale, which leads to the PNS cooling. The kinetic energy rapidly drops because of the photodissociation of iron and the electron capture (ν_e emission) that is consistent with the shock stall. We have monitored these values in a 2D simulation and obtained a similar level of energy conservation.

**Figure 15.** Time evolution of gravitational energy (red), internal energy (green line), kinetic energy (blue line), trapped-neutrino energy (magenta line), energy released by neutrinos (cyan line), and summation of these energies (black dashed line). These quantities are determined by integration with respect to the volume included in our simulation except for energy released by neutrinos (magenta), which is $\int (L_{\nu _e}+L_{\bar{\nu }_e}) dt$ . Since gravitational energy and total energy are negative, the absolute values are shown. The violation of total energy (dashed line) remains < 3 × 10⁴⁹ erg, which is ∼ 0.03% of gravitational energy and internal energy after bounce (∼10⁵³ erg).
Download figure:
Standard image High-resolution image

A2. Comparison with AGILE

Here, we present the result of our numerical simulation with spherical symmetry and compare with the result of Adaptive Grid with Implicit Leap Extrapolation (AGILE)–IDSA code (Liebendörfer et al. 2009). AGILE is an implicit general relativistic hydrodynamics code that evolves the Einstein equations based on conservative finite differencing on an adaptive grid. We employ a 1D version of our ZEUS-2D code that has been developed to perform multidimensional supernova simulations.

We compare the evolution of a 13 M_☉ star from Nomoto & Hashimoto (1988) in Newtonian gravity from a pre-collapse model to 100 ms after bounce. We find good agreement between the results of the ZEUS-2D and AGILE during the early postbounce phase when the neutrino burst is launched and the accretion shock expands to its maximum radius. The hydrodynamic quantities are shown in Figures 16, 17, and 18.

**Figure 16.** Density (left top), entropy (right top), velocity (left bottom), and electron fraction (right bottom) as a function of enclosed mass for the result of ZEUS-2D (red lines) and AGILE (green lines). The comparison is shown at the time just after the bounce. A difference is seen in the entropy profile, which comes from the difference of the shock capturing scheme.
Download figure:
Standard image High-resolution image

**Figure 17.** Same as Figure 16 but for the time at 1 ms after bounce.
Download figure:
Standard image High-resolution image

**Figure 18.** Same as Figure 16 but for time 100 ms after bounce.
Download figure:
Standard image High-resolution image

IMPACTS OF COLLECTIVE NEUTRINO OSCILLATIONS ON CORE-COLLAPSE SUPERNOVA EXPLOSIONS

Article metrics

Permissions

Author e-mails

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION

2. NUMERICAL METHODS

2.1. Hydrodynamics

2.2. Spectral Swapping

3. RESULT

3.1. One-dimensional Models

3.1.1. 1D without Spectral Swapping

3.1.2. 1D with Spectral Swapping

3.1.3. The Progenitor Dependence

3.2. Two-dimensional Models

3.2.1. 2D without Spectral Swapping

3.2.2. 2D with Spectral Swapping

4. SUMMARY AND DISCUSSION

APPENDIX: CODE VALIDITY

A1. Conservation of Energy

A2. Comparison with AGILE

Footnotes

IMPACTS OF COLLECTIVE NEUTRINO OSCILLATIONS ON CORE-COLLAPSE SUPERNOVA EXPLOSIONS

Article metrics

Permissions

Share this article

Author e-mails

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION

2. NUMERICAL METHODS

2.1. Hydrodynamics

2.2. Spectral Swapping

3. RESULT

3.1. One-dimensional Models

3.1.1. 1D without Spectral Swapping

3.1.2. 1D with Spectral Swapping

3.1.3. The Progenitor Dependence

3.2. Two-dimensional Models

3.2.1. 2D without Spectral Swapping

3.2.2. 2D with Spectral Swapping

4. SUMMARY AND DISCUSSION

APPENDIX: CODE VALIDITY

A1. Conservation of Energy

A2. Comparison with AGILE

Footnotes