A COMPARISON OF TWO- AND THREE-DIMENSIONAL NEUTRINO-HYDRODYNAMICS SIMULATIONS OF CORE-COLLAPSE SUPERNOVAE

Ever since the first numerical simulation of core-collapse supernovae (CCSNe; Colgate & White 1966), the neutrino-driven mechanism has been the leading candidate of the explosion mechanism for more than four decades. In the long history, a very important lesson we have learned from Rampp & Janka (2000), Liebendörfer et al. (2001), Thompson et al. (2003), and Sumiyoshi et al. (2005) near the Millennium, is that the spherically symmetric one-dimensional (1D) form of this mechanism fails to explode canonical massive stars. Supported by accumulating supernova observations of the blast morphology (e.g., Wang & Wheeler 2008 and references therein), a number of multi-dimensional (multi-D) hydrodynamic simulations have been reported so far, which gives us confidence that hydrodynamic motions associated with neutrino-driven convection (e.g., Herant et al. 1994; Burrows et al. 1995; Janka & Müller 1996; Fryer et al. 2002 and see collective references in Burrows et al. 2012; Murphy et al. 2013; Couch 2013) and standing-accretion-shock instability (SASI; e.g., Blondin et al. 2003; Scheck et al. 2004, 2006; Ohnishi et al. 2006, 2007; Ott et al. 2008; Murphy & Burrows 2008; Foglizzo et al. 2006, 2007, 2012; Iwakami et al. 2008, 2009, 2013; Endeve et al. 2010, 2012; Fernández & Thompson 2009a, 2009b; Fernández 2010; Hanke et al. 2012) can help the onset of neutrino-driven explosions.

In fact, a growing number of neutrino-driven explosions have been recently obtained in the state-of-the-art two-dimensional (2D) simulations, in which spectral neutrino transport is solved with different levels of sophistication (e.g., Buras et al. 2006; Marek & Janka 2009; Müller et al. 2012, 2013; Bruenn et al. 2013; Suwa et al. 2010, 2011, 2013; Janka 2012 for a review). This success is, however, accompanying new questions. Among them,⁴ three-dimensional (3D) effects on the neutrino-driven mechanism are attracting paramount attention (e.g., see Burrows 2013 and Kotake et al. 2012 for a review). Unfortunately, however, experimental 3D models that employed a light-bulb scheme (e.g., Murphy & Burrows 2008) have provided divergent results so far. The basic result of Nordhaus et al. (2010), who were the first to point out that 3D leads to easier explosions than 2D, has been supported by the follow-up studies (Burrows et al. 2012; Dolence et al. 2013), but not by Hanke et al. (2012) and Couch (2013). On top of the urgent task of making a detailed comparison between these idealized models, self-consistent 3D simulations should be done in order to have the final word on the 3D effects.

At present, there are only a few 3D CCSN simulations including spectral neutrino transport (Hanke et al. 2013; Takiwaki et al. 2012). Very recently, Hanke et al. (2013) succeeded in performing 3D simulations with detailed neutrino transport for a 27 M_☉ star. In addition to the first discovery regarding the violent SASI activity in self-consistent 3D models, their results illuminate the importance of going beyond the prevalent light-bulb scheme; only by doing so, the nonlinear couplings, such as between core-contraction of the proto-neutron star (PNS), the accretion neutrino luminosity, and the multi-D hydrodynamic feedback of neutrino-driven convection and the SASI, can be self-consistently determined. On the other hand, there is a very high computational cost allowed Hanke et al. (2013) to focus on a single (self-consistent) 3D model, and it has not yet been clarified whether 3D helps or harms the onset of neutrino-driven explosions compared to 2D.

To address this question, we investigate how the explosion dynamics will differ from 3D to 2D by systematically changing numerical resolutions and initial seed perturbations in multi-D radiation-hydrodynamic simulations. For the multi-group neutrino transport, the isotropic diffusion source approximation (IDSA) scheme (Liebendörfer et al. 2009) is implemented in a ray-by-ray manner (e.g., see Takiwaki et al. 2012 for more details), which we have updated to achieve a very high computational efficiency. As in Takiwaki et al. (2012), we focus here on the evolution of an 11.2 M_☉ star of Woosley et al. (2002). We choose this lighter progenitor because the shock revival occurs relatively earlier after bounce (e.g., Buras et al. 2006) compared to more massive progenitor models as employed in Hanke et al. (2013). The updated transport scheme, together with the employed earlier-to-explode progenitor, allows us to conduct a systematic numerical study, for the first time, in both 2D and 3D (to the best of our knowledge) in the context of self-consistent neutrino-driven supernova models.

Here, we briefly summarize several major updates of the code that we have implemented after our previous work (Takiwaki et al. 2012) in which the spectral neutrino transport scheme IDSA (Liebendörfer et al. 2009) was implemented in the ZEUS-MP code (Hayes et al. 2006).

In the original IDSA scheme, a steady-state approximation (∂f^s/(∂t) = 0) is assumed. Here, f^s represents the streaming part of the neutrino distribution function (e.g., Liebendörfer et al. 2009). Then, one should deal with a Poisson-type equation to find the solution of f^s (e.g., Equation (10) in Liebendörfer et al. 2009). This is computationally expensive, because a collective data-communication is required on the message passing interface (MPI) routines for all of the processors (along the given radial direction in the ray-by-ray approximation) to solve the Cauchy problem.

To get around the problem, one needs to solve the evolution of f^s as,

$$\begin{eqnarray} \frac{\partial \mathcal {E}^{\rm s}}{c\partial {t}} &+& \frac{1}{r^2}\frac{\partial }{\partial r}r^2 \mathcal {F}^{\rm s} = \mathcal {S}[{j}, {\chi }, {\Sigma }], \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \mathcal {E}^{\rm s} &\equiv & \frac{1}{2}\int d{\mu }\,{f}^{{\rm s}}, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \mathcal {F}^{\rm s} &\equiv & \frac{1}{2}\int \mu d{\mu }\,{f}^{{\rm s}}, \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} \mathcal {S} &\equiv & -(\hat{j}+\hat{\chi })\mathcal {E}^{\rm s}+\Sigma, \end{eqnarray} \tag{ 4 }$$

where $\mathcal {E}^{\rm s}$ and $\mathcal {F}^{\rm s}$ correspond to the radiation energy and flux of the streaming particle, and $\mathcal {S}$ represents the source term that is a functional of the neutrino emissivity (j), absorptivity (χ), and the isotropic diffusion term (Σ) all defined in the laboratory frame, respectively. From local hydrodynamic quantities (density, Y_e, entropy), the source term of Equation (1) ($\mathcal {S}[{j}, {\chi }, {\Sigma }]$) can be determined. For closure, we use a prescribed relation between the radiation energy and flux as $\mathcal {F}^{\rm s}/\mathcal {E}^{\rm s}= ({1}/{2}) (1 +\sqrt{\vphantom{A^A}\smash{{{1- [{R_{\nu }}/\max \left(r,R_\nu \right)]^2}}}})$ with R_ν being the radius of an energy-dependent scattering sphere (see Equation (11) in Liebendörfer et al. 2009). Since the cell-centered value of the flux, $\mathcal {F}^{\rm s}$, is obtained by the prescribed relation, the cell-interface value is estimated by the first-order upwind scheme assuming that the flux is out-going along the radial direction. With the numerical flux, the transport equation of $\mathcal {E}^{\rm s}$ (Equation (1)) now expressed in a hyperbolic form is numerically solved. This modification does not produce any significant changes in the numerical results (see T. Takiwaki et al. 2014, in preparation for more details), however, the computational cost becomes more than 10 times smaller than that in the previous treatment. The velocity dependent terms (O(v/c)) are only included (up to the leading order) in the trapped part of the distribution function (Equation (15) in Liebendörfer et al. 2009). Concerning heavy-lepton neutrinos ($\nu _x = \nu _{\mu },\nu _{\tau },\bar{\nu }_{\mu }, \bar{\nu }_{\tau }$), we employ a leakage scheme to include the ν_x cooling via pair, photo, and plasma processes (e.g., Rosswog & Liebendörfer 2003; Itoh et al. 1989). 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. To improve the accuracy of total energy conservation, we follow the prescription proposed by Müller et al. (2010). For the calculations presented here, self-gravity is computed by a Newtonian monopole approximation. We use the equation of state (EOS) by Lattimer & Swesty (1991) with a compressibility modulus of K = 180 MeV (LS180).

Our fiducial 3D models are computed on a spherical polar grid with a resolution of n_r × n_θ × n_ϕ = 320 × 64 × 128, in which unequal spatial radial zones cover from the center to an outer boundary of 5000 km. The radial grid is chosen such that the resolution Δr is better than 2 km in the PNS interior and typically better than 5 km outside the PNS. For the spectral transport, we use 20 logarithmically spaced energy bins ranging from 3 to 300 MeV, and we take a ray-by-ray approximation (e.g., Buras et al. 2006; Bruenn et al. 2013), in which a ray is cast for every angular zone. In all the multi-D runs, the innermost 5 km is computed in spherical symmetry to avoid excessive time-step limitations. Seed perturbations for aspherical instabilities are imposed by hand at 10 ms after bounce by introducing random perturbations of 1% in velocity behind the stalled shock.

To test the sensitivity of the supernova dynamics to numerical resolutions, we compute 3D model series with lower angular resolutions, namely, half or quarter of the (equidistant) mesh numbers in the azimuthal direction (n_r × n_θ × n_ϕ = 320 × 64 × 64 and n_r × n_θ × n_ϕ = 320 × 64 × 32). In 2D simulations, we vary the mesh numbers in the lateral direction as n_r × n_θ = 320 × 64, n_r × n_θ = 320 × 128 and n_r × n_θ = 320 × 256, respectively (see Table 1). In the table, model 3D-H-1 differs from model 3D-H-2 (and 3D-H-3 etc.) only in the random seed perturbations (with the perturbation amplitudes being the same in all cases). Note that the lowest-resolution 3D model in this work corresponds to the best-resolution model in Takiwaki et al. (2012). By using the fastest K computer in Japan, it typically took 1.3 months (equivalently ∼4 million core-hour computing-time) for each of our 3D fiducial models.

Table 1. Model Summary

Model	nr × nθ × nϕ	Angular	t400	t400, av	σ[t400]	tend	Rmax	Ediag	(l, m)max	|cl, m|max
		Resolution (°)	(ms)	(ms)	(ms)	(ms)	(km)	(1050 erg)
3D-H-1	320 × 64 × 128	28 × 28	223	236	24	284	550	0.15	1,0	0.02
3D-H-2	320 × 64 × 128	28 × 28	216			369	850	0.25	2,0	0.02
3D-H-3	320 × 64 × 128	28 × 28	269			269	400	0.15	2,0	0.02
3D-M-1	320 × 64 × 64	28 × 56	192	194	1.2	269	600	0.25	1,0	0.01
3D-M-2	320 × 64 × 64	28 × 56	194			319	700	0.30	1,0	0.02
3D-M-3	320 × 64 × 64	28 × 56	195			279	700	0.27	1,0	0.01
3D-L-1	320 × 64 × 32	28 × 113	193	188	4	314	1000	0.60	2,1	0.009
3D-L-2	320 × 64 × 32	28 × 113	188			273	800	0.45	2,1	0.009
3D-L-3	320 × 64 × 32	28 × 113	183			273	700	0.35	2,1	0.01
2D-H-1	320 × 256	07	138	151	7	300	1100	0.65	1,0	0.05
2D-H-2	320 × 256	07	154			329	1200	0.45	2,0	0.08
2D-H-3	320 × 256	07	159			311	1200	0.76	2,0	0.1
2D-H-4	320 × 256	07	156			368	1500	0.6	2,0	0.05
2D-H-5	320 × 256	07	150			345	1200	0.7	2,0	0.1
2D-M-1	320 × 128	14	140	142	10	369	1700	0.5	1,0	0.1
2D-M-2	320 × 128	14	125			319	1800	0.8	2,0	0.06
2D-M-3	320 × 128	14	151			400	1700	0.85	1,0	0.1
2D-M-4	320 × 128	14	144			469	2300	1.0	1,0	0.1
2D-M-5	320 × 128	14	152			400	1800	1.0	2,0	0.05
2D-L-1	320 × 64	28	137	133	6	387	2100	1.1	2,0	0.08
2D-L-2	320 × 64	28	137			395	2200	1.26	2,0	0.1
2D-L-3	320 × 64	28	126			483	2800	1.3	1,0	0.09
2D-L-4	320 × 64	28	140			559	2400	1.3	1,0	0.09
2D-L-5	320 × 64	28	125			569	2500	1.3	2,0	0.05

Notes. Note that "H," "M," and "L" appended to our models stand for high, moderate, and low azimuthal angular resolutions, respectively. The third number (−i) of each model, which runs from 1 to 3 for 3D models and from 1 to 5 for 2D models, represents the difference only in the random seed perturbations (with the perturbation amplitudes being the same in all cases). t₄₀₀ represents the time when the average shock radius touches a radius of 400 km, and t_{400, av} and σ[t₄₀₀] denote the model average and the dispersion of t₄₀₀. t_end denotes the time of the end of the simulation. R_max and E_diag are the averaged shock radius and the diagnostic explosion energy at t_end. (ℓ, m)_max denotes the set of spherical harmonics mode when the normalized amplitude takes the maximum value (|c_{ℓ, m}|_max) during the linear SASI phase.

Download table as:  ASCIITypeset image

As summarized in Table 1, all the computed models including nine 3D models and fifteen 2D models exhibit shock revival, leading to the possibility of explosion. Before going into detail how the explosion dynamics and stochasticity are different in 2D and 3D, we briefly outline the hydrodynamic features using model 3D-H-1 as an example.

The top panel of Figure 1 shows the blast morphology of model 3D-H-1 at t_pb = 230 ms (postbounce) when the revived shock is reaching an angle-averaged radius of 400 km (e.g., red dashed line in the right panel of Figure 2). As seen from the side wall panels, a bipolar explosion is obtained for this model. The bottom left panel (red regions) shows that the ratio of the residency timescale to the neutrino-heating timescale (e.g., Equations (6) and (7) in Takiwaki et al. 2012) exceeds unity behind the shock, which presents evidence that the shock revival is driven by the neutrino-heating mechanism. The bottom right panel of Figure 1 depicts spatial distribution of the net neutrino-heating rate at t_pb = 150 ms. Small scale inhomogeneities (colored as red or yellow) are seen, which predominantly comes from neutrino-driven convection and anisotropies of the accretion flow, but the shape of the gain region is very close to being spherical before the onset of an explosion. This suggests that the bipolar geometry of the shock is produced not by the global anisotropy of the neutrino heating in the vicinity of the neutrino sphere, but by multi-D effects such as by neutrino-driven convection and the SASI in the gain regions after the explosion (gradually) sets in.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Reflecting the stochastic nature of the multi-D neutrino-driven explosions, the blast morphology changes from model to model. The left and middle panels of Figure 2 show that a stronger explosion is obtained toward the north direction (model 3D-H-2) and the south pole (model 3D-H-3), which is only different from model 3D-H-1 (e.g., Table 1) in terms of the imposed initial random perturbations. Note that, due to the use of the spherical coordinates, we cannot omit the possibility that the polar axis still gives a special direction in our 3D simulations. But more importantly, our results show that the final fate of the 3D and 2D models whether an explosion is obtained or not, is little affected by the stochasticity of the explosion geometry.

In fact, the right panel of Figure 2 shows the evolution of the average shock radius for our 1D (blue line), 2D (green lines), and 3D (red lines) models. Before the onset of shock revival (before 100 ms after bounce), the evolution of the shocks are all similar to that of the 1D model (blue line). After that, our results show that the shock expansion is systematically more energetic in 2D (green lines) than in 3D (red lines). This feature is qualitatively consistent with Takiwaki et al. (2012) and also with Hanke et al. (2013) who recently reported 2D versus 3D comparison based on a single 3D model but employing a more detailed neutrino transport than ours.

Due to our lack of necessary computational resources, our 3D models should be terminated typically before t_end ≲ 300 ms postbounce (e.g., in Table 1) but we expect them to produce explosions subsequently, seeing a continuous shock expansion out to a radius of 500 km in 3D (and 700 km in 2D). Given the same numerical resolution (e.g., model series 3D-L in Table 1), the average shock radii in this study is smaller than those in Takiwaki et al. (2012), in which the cooling by heavy-lepton neutrinos was not taken into account. Due to the inclusion of the ν_x cooling, the (angle-averaged) $\bar{\nu }_e$ luminosity decreases more quickly after bounce (compare Figure 3 and Figure 14 of Takiwaki et al. 2012), which leads to the less energetic shock expansion in this study. It should be mentioned that, by comparing our ν_x luminosity estimated by the leakage scheme with that obtained by the work of Buras et al. (2006) with detailed neutrino transport, the peak luminosity is more than 20% smaller in our case. Such underestimation of cooling by heavy-lepton neutrinos should lead to artificially larger maximum shock extent (R_max ∼ 260 km, blue line in the right panel of Figure 2) compared to R_max ∼ 170 km in Buras et al. (2006). We have to emphasize that the use of the leakage scheme, together with the omission of inelastic neutrino scattering on electrons and general relativity (GR) effects in the present scheme, is likely to facilitate artificially easier explosions. Regarding our 2D models, the relatively earlier shock revival (∼100 ms postbounce) coincides with the decline of the mass accretion rate onto the central PNS. This could be the reason that the timescale is similar to that in Müller et al. (2012) who reported 2D (GR) models for the same progenitor model with detailed neutrino transport.

As seen from Figure 3, the angle-averaged neutrino luminosity (〈L_ν〉) and the mean neutrino energy ($\langle \epsilon _{\nu } \rangle = \int E^3 \mathcal {F}^{\rm s}dE /\int E^2 \mathcal {F}^{\rm s}dE$, where E is neutrino energy) are barely affected by the imposed initial perturbations (presumably at a few-percent levels in amplitude). This again supports our finding that the explosion stochasticity is very influential in determining the blast morphology but not the working of the neutrino-heating mechanism.

Zoom In Zoom Out Reset image size

From the bottom panel of Figure 3, it can be seen that the overall trend in the neutrino luminosities and the mean energies is similar between our 3D and 2D models. The neutrino luminosities in the 2D model (green lines) show a short-time variability (with periods of milliseconds to ≳ 10 ms) after around 100 ms postbounce. Such fast variations in the postbounce luminosity evolution have been already found in previous 2D studies (e.g., Ott et al. 2008; Marek et al. 2009). This is caused by the modulation of the mass accretion rate due to convective plumes and downflows hitting onto the PNS surface (see also Lund et al. 2012 and Tamborra et al. 2013 about the detectability of these neutrino signals). It is interesting to note that such a fast variability is less pronounced in our 3D model (red lines in the bottom panel). This is qualitatively consistent with Lund et al. (2012) who analyzed the neutrino signals from 2D and 3D models, in which an approximate neutrino transport was solved (Wongwathanarat et al. 2010) as in Scheck et al. (2006).

Figure 4 shows the evolution of the average PNS radius for the 1D (blue line), 2D (green line), and 3D models (red line) that are defined by a fiducial density of 10¹¹ g cm⁻³. The PNS contraction is similar among the 1D, 2D, and 3D models. Although the PNS contraction potentially affects the evolution of the shock (Hanke et al. 2013; Suwa et al. 2013), in our cases, it is unchanged by the difference of the dimension and is not the main agent to explain the divergence of the shock evolution in 1D, 2D, and 3D. The PNS contraction is slightly stronger in the later postbounce phase in 1D (≳ 150 ms postbounce, compare blue with green and red lines) compared to 2D and 3D because no shock revival was obtained in the 1D model, and heavier PNS and slightly deeper gravitational potential are obtained compared to that of the multi-D models. In the figure, three more lines (solid, dashed, dotted gray lines) are plotted, in which we estimate the evolution of the PNS radius based on the fitting formula (Equation (1) of Scheck et al. 2006) by changing a final radius of PNS R_f for a given set of an exponential timescale of t_ib = 1 s and an initial radius of PNS R_i = 85 km. As can be seen, the dashed gray line (R_f = 12 km) can most closely reproduce our results, which is just between the slow and fast contraction investigated in the work by Hanke et al. (2013).

Zoom In Zoom Out Reset image size

The top panel of Figure 5 shows angle-averaged entropy profile at 100 ms postbounce, after the difference of the subsequent shock evolution between our 1D, 2D, and 3D models becomes remarkable (e.g., right panel in Figure 2). As has been studied in detail since the 1990s (e.g., Herant et al. 1994; Burrows et al. 1995; Janka & Müller 1996), buoyancy-driven convection supported by turbulence (e.g., Murphy et al. 2013) transports heat radially outward, leading to a more extended entropy profile in the 2D (green line in the panel) and 3D (red line) model compared to the 1D model (blue line; see also Hanke et al. 2012).

Zoom In Zoom Out Reset image size

The bottom panel of Figure 5 compares the turbulent energy spectra of the anisotropic velocity (Takiwaki et al. 2012) as a function of wavenumber (k) between our 2D (green line) and 3D (red line) model. The green and red lines cross at around k_cross = 0.02 km⁻¹ (corresponding to ∼50 km in spatial scale), above which the amplitude for the 3D model (red line) dominates that of the 2D model (green), and below which the amplitude for the 2D model dominates that of the 3D model. This is qualitatively consistent with the previous results using the light-bulb method (e.g., Hanke et al. 2012; Burrows et al. 2012; Dolence et al. 2013; Couch 2013; Fernández et al. 2013). By comparing the 2D and 3D curve with spectral slopes labeled by the corresponding exponent (−5/3: dotted blue line, −3: dotted black line), the power-low dependence (∝k^−5/3) approximately holds for k_int ≲ 0.02 km⁻¹ in the 2D model (green line) presumably as a result of the inverse energy cascade (Kraichnan 1967). Above k_cross, the slope of the 3D model (red line) becomes closer to k^−5/3 (blue line), while the spectrum slope of the 2D model drops much more steeply with the wave number as k⁻³ (gray line). These features are again in accord with the previous studies mentioned above (e.g., Hanke et al. 2012; Burrows et al. 2012; Dolence et al. 2013; Couch 2013; Fernández et al. 2013).

As shown in the middle panel of Figure 2, model 3D-H-3 produces a one-sided explosion toward the south pole during the simulation time. Reflecting the unipolarity, the average shock radius is smaller than for the other 3D models (compare dotted line (in red) with solid and dashed red line (in red) in the right panel of Figure 2). To quantify the vigor of the shock expansion, t₄₀₀ is a useful quantity that was defined in Hanke et al. (2012) as the moment of time when the shock reaches an average radius of 400 km. In fact, as seen from Table 1, t₄₀₀ of model 3D-H-3 (t₄₀₀ ≈ 270 ms) is delayed about 50 ms compared to those of models 3D-H-1 and 3D-H-2 (t₄₀₀ ≈ 220 ms).

The model average of t₄₀₀ (e.g., t_{400, av} in Table 1) clearly shows that models with higher numerical resolutions lead to a slower onset of the shock revival in both our 3D and 2D models. This feature is qualitatively consistent with the 2D self-consistent models by Marek & Janka (2009) and with the 3D idealized models by Hanke et al. (2013) and Couch (2013).

In Table 1, σ₄₀₀ represents the model dispersion of t₄₀₀, which varies much more stochastically in 3D models with different numerical resolution (from 1.2 to 24 ms) than those in 2D models (from 6 to 10 ms). For model 3D-H-3, the shock revival is the most delayed (e.g., Table 1) and the shock expansion is weakest among the computed models (dotted red line in Figure 1). Nagakura et al. (2013) proposed that shock revival is very sensitive to the imposed seed perturbations near the stalled shock. Following the hypothesis, we speculate that the influence of seed perturbations is seen most remarkably in our weakest explosion model.

In Table 1, E_diag denotes the diagnostic energy defined as the total energy (internal plus kinetic plus gravitational) integrated over all matter where the sum of the corresponding specific energies is positive. We include recombination energy in internal energy (Bruenn et al. 2013). Reflecting the earlier shock revival, the diagnostic energy is systematically bigger in 2D than in 3D models. When we terminated the simulation, these diagnostic energies were typically on the order of ∼10⁴⁹ erg and ∼10⁵⁰ erg for our 3D and 2D models, respectively. It should be noted that this quantity is estimated at the end of simulation, t_end. In order to compare the diagnostic energy with the observed kinetic explosion energy (∼10⁵¹ erg), a much longer-term simulation including improved microphysics, GR, and nuclear burning would be needed.

Recently, it is enthusiastically discussed whether the neutrino-driven convection or the SASI plays a more crucial role in facilitating neutrino-driven explosions (Foglizzo et al. 2006; Dolence et al. 2013; Murphy et al. 2013; Burrows et al. 2012; Hanke et al. 2012, 2013). The left panel of Figure 6 shows the time evolution of the Foglizzo parameter χ (Foglizzo et al. 2006). As seen, χ continuously exceeds the critical number of 3 (Foglizzo et al. 2006) rather shortly after bounce (∼40 ms), marking the transition to the nonlinear phase. The earlier onset of the shock revival and the absence of clear features of the SASI in the linear phase (see discussion below) could indicate that neutrino-driven convection dominates over the SASI when the explosion sets in (e.g., at t₄₀₀ in Table 1).

Zoom In Zoom Out Reset image size

The right panel of Figure 6 shows the evolution of the normalized coefficient of spherical harmonics of the shock surface for model 3D-H-1. No clear feature of the linear growth of the SASI (≲ 40 ms postbounce) was obtained for the 11.2 M_☉ star explored in this work. In both our 3D and 2D models, the qualitative behaviors of the harmonic modes seen in Figure 6 are less sensitive to the employed numerical resolution and seed random perturbations (see the Appendix). The dominant channels are of low modes (e.g., (ℓ, m)_max in Table 1), the amplitudes of which (e.g., |c_{ℓ, m}|_max in the table) are systematically larger in 2D than in 3D. This is qualitatively consistent with previous 3D simulations employing different numerical setups (e.g., Nordhaus et al. 2010; Hanke et al. 2012, 2013; Burrows et al. 2012).

After the linear phase comes to an end at around 100 ms postbounce, the trajectory of the revived shock shows a wider diversity (e.g., Figure 1) depending on the employed numerical resolution and seed perturbations.

Studying an 11.2 M_☉ and changing numerical resolutions and seed perturbations systematically in the multi-D simulations employing the updated IDSA scheme, we studied how the postbounce dynamics are different in 2D and 3D models. All of the computed models exhibit the neutrino-driven revival of the stalled bounce shock, leading to the possibility of an explosion. Though the blast morphology changes from model to model reflecting the stochastic nature of multi-D neutrino-driven explosions, it was found that the final fate of these multi-D models whether an explosion is obtained or not, is little affected by the explosion stochasticity at least in the current investigated progenitor model. In line with some previous studies, higher numerical resolutions lead to slower onset of the shock revival in both 3D and 2D. Our results systematically showed that the revived shock expands more energetically in 2D than in 3D.

The caveats of our 3D models include the ray-by-ray approximation, the use of the softer EOS, and the omission of detailed neutrino reactions and GR (e.g., Kuroda et al. 2012; Ott et al. 2013). Keeping our efforts to improve them, it is important to study the dependence of progenitors (e.g., Buras et al. 2006; Bruenn et al. 2013) and EOS (e.g., Marek & Janka 2009; Suwa et al. 2013) on the neutrino-driven mechanism in 3D computations. A number of exciting issues still remain to be investigated, such as gravitational-wave signatures (e.g., Kotake 2013), neutrino emission and its detectability (e.g., Lund et al. 2012), and the possibility of 3D SASI flows generating pulsar kicks and spins (e.g., Wongwathanarat et al. 2013). Shifting from individuals to populations of 3D models, a rush of 3D explorations with increasing sophistication are now going to shed light on these fascinating riddles (hopefully in the near future) with increasing supercomputing resources on our side.

We are thankful to K. Sato and S. Yamada for continuing encouragements. Numerical computations were carried on in part on the K computer at the RIKEN Advanced Institute for Computational Science (Proposal number hp120285), XC30/XT4 at NAOJ, and on SR16000 at YITP at Kyoto University. This study was supported in part by Grants-in-Aid for Scientific Research (Nos. 23540323, 23340069, 24244036, and 25103511) and by the HPCI Strategic Program of the Japanese MEXT.

In this Appendix, we briefly report on the dependence on the amplitudes of the initial seed perturbations. We have added a radial velocity perturbation, δv_r(r, θ, ϕ), to the profile obtained by 1D simulation, $v_r^{1D}(r)$, according to the equation $\delta v_r = p_{\rm amp}\,{\rm rand} \times v_r^{1D}$ where rand is a pseudorandom number that takes a value from −1 to 1 and p_amp represents the absolute amplitude of the seed perturbations that we take as 1% in the models discussed in the main section. Let us firstly note that the amplitude of seed perturbations, a_{ℓ, m}, assumed in this study scales as p_amp/N^1/2, where a_{ℓ, m} denotes the amplitude of the (ℓ, m) component of the seed perturbations at a given radius (i.e., $a_{\ell,m} = \int d \phi \int \sin \theta d \theta (v_r^{1D}+\delta v_r)Y_{\ell,m}^{*}$, where $Y_{\ell,m}^{*}$ is conjugate of spherical harmonic function) and N = n_θ × n_ϕ represents the total angular mesh number where n_θ, n_ϕ is the mesh number in the θ and ϕ direction, respectively (e.g., in our 3D model with highest resolution, n_θ = 64, n_ϕ = 128). Therefore, the seed amplitudes assumed in this work depend on numerical resolution (i.e., N).

If we keep the seed amplitudes (a_{ℓ, m}) constant for models with different resolutions, will our results drastically change? Since we cannot rerun 3D models due to the limited computational resources, we compute a number of 2D models to answer this question. First of all, let us discuss 2D models with different seed amplitudes. The top panel of Figure 7 shows the evolution of c_{1, 0} (the normalized harmonic amplitudes of shock position) for 2D models with p_amp = 1% (red lines) or $0.176 \% (= 1\%/\sqrt{32})$ (green lines) for five different realizations of initial random perturbation. The reduced amplitude is determined by the ratio of total grid number of 2D-H models and 3D-H models (i.e., $1/\sqrt{32} = \sqrt{\vphantom{A^A}\smash{{{ n_\theta |_{2D-H}/(n_\theta n_\phi)|_{3D-H}}}}}$). One might guess that c_{ℓ, m}(t) could evolve as a_{ℓ, m}exp (t/t₀) with t₀ representing the duration of the linear growth rate of hydrodynamic instabilities including the SASI or neutrino-driven convection. If this were the case, the linear growth amplitude (before 35 ms postbounce in the top panel) should be higher by about $\sqrt{32}$ times for red lines compared to green lines. But as is shown, this is not the case.⁵ After the early rising phase (about 40 ms postbounce in the panel), the saturation amplitudes are shown to be insensitive to the initial perturbation amplitudes (for the initial strength employed above). Remembering that neutrino-driven convection is likely to dominate over the SASI in the nonlinear regime for the 11.2 M_☉ progenitor employed in this work, the delay of t₄₀₀ (in the nonlinear phase) for models with higher numerical resolutions cannot be simply ascribed to the difference of the seed amplitudes.

Zoom In Zoom Out Reset image size

The bottom panels of Figure 7 show the evolution of c_{1, 0} for different resolutions (H, M, and L) between our 2D and 3D models. For a given numerical resolution, t₄₀₀ of the chosen models in these plots is close to a median of t₄₀₀ for each model series. Note that in these models, the initial seed amplitudes are dependent on both resolution and dimension. As can be seen, in both 2D and 3D, the amplitudes in the linear phase (<40 ms postbounce) are comparable for models with different resolutions. In the nonlinear regime, monotonic dependence of the nonlinear evolution on the initial seed amplitudes (between H, M, and L models) cannot be found. These results show that the findings in this work are less sensitive to the assumed initial seed perturbations.