{'fraction_non_alphanumeric': 0.03277005116379621, 'fraction_numerical': 0.02642862290120343, 'mean_word_length': 4.568161300030093, 'pattern_counts': {'":': 1, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the Impact of Three Dimensions in Simulations of Neutrino-Driven Core-Collapse Supernova Explosions

We present 1D, 2D, and 3D hydrodynamical simulations of core-collapse supernovae including a parameterized neutrino heating and cooling scheme in order to investigate the critical core neutrino luminosity (L_crit) required for explosion. In contrast to some previous works, we find that 3D simulations explode later than 2D simulations, and that L_crit at fixed mass accretion rate is somewhat higher in 3D than in 2D. We find, however, that in 2D L_crit increases as the numerical resolution of the simulation increases. In contrast to some previous works, we argue that the average entropy of the gain region is in fact not a good indicator of explosion but is rather a reflection of the greater mass in the gain region in 2D. We compare our simulations to semi-analytic explosion criteria and examine the nature of the convective motions in 2D and 3D. We discuss the balance between neutrino-driven-buoyancy and drag forces. In particular, we show that the drag force will be proportional to a buoyant plume's surface area while the buoyant force is proportional to a plume's volume and, therefore, plumes with greater volume-to-surface area ratios will rise more quickly. We show that buoyant plumes in 2D are inherently larger, with greater volume-to-surface area ratios, than plumes in 3D. In the scenario that the supernova shock expansion is dominated by neutrino-driven buoyancy, this balance between buoyancy and drag forces may explain why 3D simulations explode later than 2D simulations and why L_crit increases with resolution. Finally, we provide a comparison of our results with other calculations in the literature.


INTRODUCTION
Multidimensional phenomena play a critical role in the corecollapse supernova (CCSN) mechanism. Instabilities such as proto-neutron star convection (Epstein 1979;Burrows & Fryxell 1993), neutrino-driven convection (Herant et al. 1994;Burrows et al. 1995;Janka & Mueller 1996;Murphy et al. 2013), and the standing accretion shock instability (SASI, Blondin et al. 2003;Blondin & Mezzacappa 2006) have abetted explosions in 2D (Marek & Janka 2009;Müller et al. 2012b,a) for progenitors that refuse to explode in spherical symmetry (Rampp & Janka 2000;Liebendörfer et al. 2001;Thompson et al. 2003;Suwa et al. 2010). The additional degrees of freedom afforded by multiple dimensions can also increase the dwell time of matter in the post-shock region where the accreting matter experiences a net gain of neutrino energy resulting in an increased efficiency of neutrino heating. Taken together, these multidimensional effects lower the critical neutrino luminosity threshold for which explosions are obtained when comparing 2D to 1D (Murphy & Burrows 2008). In spherical symmetry, the radial stability of the supernova shock has been studied in detail (Burrows & Goshy 1993;Yamasaki & Yamada 2007;Pejcha & Thompson 2012;Fernández 2012). The current state of the research into the CCSN mechanism has been reviewed by Janka (2012), Burrows (2013), .
We lack a first principles understanding of how the aforementioned multidimensional effects result in a lowered critical threshold for explosion in 2D, however, it is attractive to extrapolate the trend to 3D and postulate that the threshold for explosion will be reduced even further, perhaps permitting energetic explosions for the vast majority of progenitors. Nordhaus et al. (2010) found precisely this using a simplified prescription for 1 Hubble Fellow neutrino heating and cooling and deleptonization. Nordhaus et al. found that the critical luminosity in 3D was reduced by 15% -25% as compared to 2D simulations for a 15 M progenitor. Using a similar parametric approach, Hanke et al. (2012) attempted to reproduce the results of Nordhaus et al. for both the 15 M progenitor and a 11.2 M progenitor. Contrary to the findings of Nordhaus et al., Hanke et al. find that there is little difference between the critical luminosities in 2D and 3D, while recovering the result that the critical luminosity in 2D is significantly lower than that in 1D (Murphy & Burrows 2008). While differences exist between the approaches of Nordhaus et al. and Hanke et al., the exact cause of their disparate results is unclear. Burrows et al. (2012) report that the results of Nordhaus et al. were beset by inaccuracies in the gravity solver in CASTRO that have since been corrected, and very recently Dolence et al. (2013) present new 3D CASTRO simulations showing faster shock expansion in 3D than in 2D. The important question of whether the threshold for explosion is lower in three dimensions is, thus, still an open one.
We are on the precipice of achieving 3D simulations of corecollapse supernovae with full spectral neutrino transport and adequate resolution, though a number of 3D CCSN simulations employing various approximations have already been accomplished. The 3D simulations to-date have been facilitated by approximations, typically to the neutrino transport, or by low resolution. Here we briefly mention a non-exhaustive list of 3D results relevant to CCSNe extant in the literature. The 3D calculations of Mueller & Janka (1997) and Khokhlov et al. (1999) did not include the effects of neutrinos and excised the PNS from the domain. A number of early 3D simulations utilized the smoothed-particle hydrodynamics approximation (Fryer & Warren 2002, 2004Fryer & Young 2007), which has certain disadvantages over grid-based codes (see, e.g., Plewa 2001;Agertz et al. 2007;McNally et al. 2012;Sijacki et al. 2012). Studies of the SASI have made various approxima-tions to achieve 3D simulations (Blondin & Mezzacappa 2007;Iwakami et al. 2008Iwakami et al. , 2009Fernández 2010). Simulations neglecting the effects of neutrinos and employing a simplified equation of state (EOS) have been used to study the amplification of magnetic fields in 3D (Endeve et al. 2010(Endeve et al. , 2012. Hammer et al. (2010), using a neutrino "lightbulb" scheme (Scheck et al. 2006), followed the evolution of the supernova explosion all the way through the envelope of the progenitor in 3D and examined the asymmetric development of instabilities. Some studies have focussed on highly-magnetized progenitors with approximate (or no) treatments of neutrino effects (Kuroda & Umeda 2010;Winteler et al. 2012). The hydrodynamical kick imparted to the PNS has been studied in 3D by Wongwathanarat et al. (2010Wongwathanarat et al. ( , 2013. The rotational stability of the PNS has been explored in 3D by Ott et al. (2005). Many 3D calculations have focussed on the emergent gravitational wave signal from CCSNe (Ott et al. 2007(Ott et al. , 2011(Ott et al. , 2012Scheidegger et al. 2008Scheidegger et al. , 2010bKotake et al. 2009;Müller et al. 2012c). Ott et al. (2013) examined the development of the SASI in 3D general relativistic simulations with neutrino leakage. Some studies have focused also on the neutrino signal from 3D CC-SNe (Lund et al. 2012;Ott et al. 2012). Takiwaki et al. (2012) present low-resolution 3D simulations with spectral neutrino transport using the isotropic diffusion source approximation (IDSA, Liebendörfer et al. 2009). And recently fully-3D Boltzmann transport for neutrinos has been developed by Sumiyoshi & Yamada (2012).
In this paper, we describe our multidimensional study of neutrino-driven CCSN explosions using a parameterization of the neutrino effects similar to that of Nordhaus et al. (2010) and Hanke et al. (2012). We find that the delay time until explosion for a given neutrino luminosity is greater in 3D than in 2D, i.e., L crit is greater in 3D than in 2D. In Section 2 we describe our computational approach. In Section 3, we present our results. In Section 4 we discuss the dependence of the critical neutrino luminosity on dimensionality and resolution and suggest that our results can be understood by considering the balance between buoyant and drag forces acting on neutrino-driven bubbles. We also examine the difference in the character of the shock motion between 2D and 3D. We demonstrate the resolution dependence and convergence of our results in Section 5. In Section 6 we discuss the implications of our results and we conclude in Section 7.
2. NUMERICAL METHOD Our numerical simulation approach is similar to that described by Couch (2013). We solve the Eulerian equations of hydrodynamics, where ρ is the mass density, v v v the velocity vector, P the pressure, Φ the gravitational potential, E the total specific energy, H is the specific neutrino heating, and C is the specific neutrino cooling. We use the directionally-unsplit hydrodynamics solver provided by the FLASH simulation framework (Dubey et al. 2009, , Lee et al., in prep.) to solve equations (1) -(3). We use third-order piecewise-parabolic spatial reconstruction (PPM, Colella & Woodward 1984) and a hybrid Riemann for the 3D model with L νe,52 = 1.7 at 850 ms post-bounce. The shock front is also volume-rendered in gray-scale. The blue sphere is a isodensity contour marking the edge of the proto-neutron star. The explosion develops in a non-symmetric manner with the PNS recoiling in the direction opposite the dominant direction of the explosion. The high-entropy buoyant plumes display a great deal of small-scale structure in 3D.
solver that uses the HLLE solver inside of shocks and the HLLC solver in smooth flow. We use a "hybrid" slope limiter that applies the monotonized central (mc) limiter to linear wave families and the more diffusive minmod limiter to nonlinear, self-steepening wave families. We use a monopole approximation to calculate the self-gravity of the flow. To facilitate comparison with Nordhaus et al. (2010) and Hanke et al. (2012), we use the Shen et al. (1998)  Our approach for neutrino heating and cooling is that described by Murphy & Burrows (2008), the same approach used by Nordhaus et al. (2010) and Hanke et al. (2012). The neutrino heating and cooling are given by, and where L νe is the electron neutrino luminosity (note it is assumed that Lν e = L νe ), T ν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, τ νe is the electron neutrino optical depth, and T is the matter temperature. In all of our simulations, T νe is set to 4 MeV. We approximate the neutrino optical depth by a piecewise fit based only on the density (see Couch 2013). This avoids the need to calculate radial integrals for the optical depth and is justified because the factor e −τν e is included in equations (4)  We follow the approach proposed by Liebendörfer (2005) for following the evolution of the electron fraction, Y e . In this approach, calibrated with 1D Boltzmann transport simulations, Y e is dependent only on density. This is strictly only applicable during the pre-bounce collapse phase, however, we continue to use the density-dependent electron fraction approach postbounce, as done by Nordhaus et al. (2010). We have found that neglecting any changes in Y e post-bounce results in substantially earlier explosions for a given neutrino luminosity. We do not include the entropy changes due to deleptonization given in Liebendörfer (2005).
We use 1D spherical, 2D cylindrical, and 3D Cartesian geometries with adaptive mesh refinement (AMR) as implemented in FLASH via PARAMESH (v.4-dev, MacNeice et al. 2000). For this study we use a fiducial resolution at the maximum refinement level of 0.7 km in each direction. We limit the maximum refinement level with radius such that a pseudologarithmic radial grid spacing is obtained. Our refinement limiter takes the form where ∆x i is the grid spacing in the i-direction at refinement level , r is the spherical radius, and η is a parameter that sets the effective angular spacing. If equation (6) is not satisfied by a given AMR block, further refinement of that block is prohibited. For our fiducial resolution we set the finest grid spacing to 0.7 km and η = 1.25%, resulting in an effective "angular" resolution of 0. • 54. In 1D, the simulated domain spans 0 km to 5000 km, in 2D the domain is 0 km to 5000 km in cylindrical radius, R, and -5000 km to 5000 km in z, and in 3D the domain is -5000 km to 5000 km in each Cartesian dimension. At the outer spatial limits of the domain, we set boundary conditions that apply power-law profiles to density and velocity that approximate the stellar envelope outside the domain. Such boundary conditions are critically important to the results of the present study as simple "outflow" boundary conditions overestimate the mass accretion rate at late times, altering the explosion time for near-critical luminosities. This is because "outflow" boundary conditions enforce a zero-gradient condition for the flow variables which mimics a flat density (etc.) profile outside the simulation domain artificially enhancing the mass flux into the domain from the boundary. We use the 15 M progenitor of Woosley & Weaver (1995) in all of our simulations.

RESULTS: EXPLOSION TIMES
We have run a series of 1D, 2D, and 3D simulations in which we varied the driving neutrino luminosity. We start in the precollapsed progenitor phase and follow the evolution through collapse, bounce, shock stagnation and eventual revival. Figure  1 shows a volume rendering of entropy and the shock surface in a 3D simulation at 850 ms post-bounce. In Table 1 we give the explosion delay times for our series of simulations and Figure 2 shows the average shock radii as a function of time post-bounce for a number of our simulations. Figure 3 shows the critical luminosity curves as functions of both postbounce explosion time and mass accretion rate at explosion. · · · · · · · · · · · · 2.1 2.2 · · · · · · · · · · · · 2.  Figure 2. Average shock radii as a function of time relative to bounce for three neutrino luminosities in 2D and 3D. Also shown for comparison is the shock radius from the 1D simulation with L νe,52 = 2.3. Universally, the shock expands more rapidly in 2D than in 3D. Increasing the resolution in 2D delays explosion, as shown by the cyan curve.
We consider a model to have exploded once the average shock radius exceeds 400 km and does not subsequently fall back below this value (as in Nordhaus et al. 2010;Hanke et al. 2012), though other metrics, such as reaching a critical value of the ratio of advection time to heating time in the gain region (e.g., Fernández 2012) or satisfying the 'ante-sonic' condition (Pejcha & Thompson 2012) may be used (for a comparison of the difference between these metrics, see Dolence et al. 2013). We find that the critical luminosity curve is lowered in multidimensional simulations as compared with sphericallysymmetric simulations, consistent with all previous similar studies (Murphy & Burrows 2008;Nordhaus et al. 2010;Hanke et al. 2012;Couch 2013). When comparing 2D to 3D, however, we find interesting and heretofore unprecedented behavior: at 1D, 0.5 km 1D, 0.7 km 2D, 0.5 km 2D, 0.7 km 3D, 0.7 km Figure 3. Critical neutrino luminosity curves in both mass accretion rate (left) and post-bounce time (right) for our set of 1D, 2D, and 3D simulations with the 15 M progenitor. We reproduce the common result that the critical curves are lower in 2D than in 1D, however, we find that for our fiducial resolution of 0.7 km the 3D critical curves are higher than the 2D curves. Thus, our results indicate that obtaining explosion is more difficult in 3D. Critical curves for increased resolution 1D and 2D simulations are also shown. Increasing resolution in 1D results in almost no difference in the curves, whereas in 2D the higher-resolution simulations explode later and the curves are much closer to the 3D counterparts.
our fiducial resolution the 2D simulations consistently explode earlier than 3D simulations at the same neutrino luminosity. Figure 2 shows that for a given neutrino luminosity the average shock radius expands more quickly in 2D than 3D. In 2D, the explosion time for a given luminosity is sensitive to the grid resolution used. Increasing the finest grid resolution to 0.5 km and reducing the radial refinement limiter, η, to 0.94% results in a 2D criticality curve much more similar to the 3D curve at the fiducial resolution of 0.7 km (and η = 1.25%). Increasing the resolution in 1D simulations results in almost no change in the explosion times as a function of neutrino luminosity. This very importantly indicates that the cause of the resolution dependence is connected to an intrinsically multidimensional process. At present, we lack the necessary computational resources to carry out a resolution study in 3D, though our results certainly indicate the necessity of such a study. Hanke et al. (2012) do carry out a resolution study including 3D and also find that the explosion times are very dependent on grid spacing. For simulations with 400 unequally-spaced radial zones, Hanke et al. find that increasing the angular resolution of their spherical grid results in earlier explosions in 2D but later explosions in 3D. Considering simulations with at least 600 radial zones, the dependence of the explosion times on resolution is not always consistent amongst their results; increasing the angular resolution delays explosion in 2D for some of their models. It is important to note that their fiducial resolution (3. • ) is substantial coarser than ours (∼ 0. • 54) and the finest resolution they use in 3D (1. • 5) is still coarser than our resolution.

INTERPRETATION AND ANALYSIS 4.1. Shock Expansion Driven by Buoyancy
So then, what is the explanation of our results? That is, why is it that our 2D simulations explode earlier than our 3D? and why does increasing the resolution in 2D result in later explosions? As we will argue in the following sections, our results indicate that neutrino-driven buoyant convection is the dominant instability that encourages shock expansion for this progenitor, particularly in 3D. Similar arguments have  Table 1). The left panel shows the 3D data and the right shows the 2D data revolved about the symmetry axis. As discussed in the text, 2D simulations show buoyant plumes that have much smaller surface area-to-volume ratios than for 3D. been made recently by Burrows et al. (2012) and Murphy et al. (2013). In this picture, accreting gas in the gain region absorbs neutrino energy eventually becoming buoyant and rises toward the shock where this buoyant energy is used to push the shock further out. The speed at which a plume rises will be determined by the competition between the plume's buoyancy, determined by the amount of neutrino energy it has absorbed, and the drag force from cold gas traveling downward through the gain region (see, e.g., Thompson 2000;Dolence et al. 2013). A plume's buoyancy will be proportional to its volume since the neutrino energy absorption rate will scale as the solid angle of the plume times its neutrino optical depth. The drag force pushing back against the buoyant plume will scale as the surface area of the plume. Therefore, smaller plumes with greater surface area-to-volume ratios will rise more slowly than larger plumes.
To see this, consider a spherical buoyant bubble of radius r b at a distance from the coordinate origin R b . The instantaneous buoyant force on this bubble is provided by neutrino radiation and so must be equal to the neutrino radiation force: whereσ ν is an effective neutrino crosssection that includes any geometric constants, ρ b is the density of the bubble, and c is the speed of light. The instantaneous drag force on the bubble is , whereC d is a drag coefficient that contains any geometric constants, and v is the bubble's velocity relative to the background accretion flow. The ratio of the buoyant force to the drag force on the bubble is then Thus, this ratio increases with bubble size, r b , and larger bubbles or plumes will rise faster than smaller. In 2D, the typical plume size is larger than in 3D for multiple reasons. First, the axial symmetry intrinsic to 2D cylindrical coordinates means that off-axis plumes are really rings. In 3D, such plume rings are unstable and will break up into many smaller plumes. In Figure 4 we demonstrate this difference in plume scale with constant entropy contours in 3D (left) and 2D revolved around the symmetry axis (right). In 2D the forced symmetry results in very large "3D" plumes whereas in 3D no such large scale plumes can exist. These contours are plotted at the respective explosion times in correspondence with the right panel of Figure 5. Second, the symmetry axis in 2D encourages the growth of large plumes along it, due either to the action of low-order modes of the SASI or that of buoyant convection. In Section 4.2 we show that the amplitudes of loworder spherical harmonic modes of the shock deformation are reduced in 3D as compared to the 2D case. Third, the "inverse" turbulent energy cascade in 2D (Kraichnan 1967) will pump energy to larger scales whereas in 3D the "forward" energy cascade will send energy to smaller scales (see Section 4.4).
In Figure 5 we compare entropy slices between 3D and 2D simulations with L νe,52 = 1.7. The left half of Figure 5 is at a time of 100 ms post-bounce and the right half is at the time of explosion, 821 ms for 3D and 388 ms for 2D. At 100 ms, the 3D simulation shows a shock that is still nearly spherical and developing convection behind it. The largest convective plumes are just reaching the shock perturbing it's spherical structure and stochastically pushing it out in radius. The 2D simulation is similar, but the developing convection in the gain region is visibly more coherent and vortex-like. Again, in 2D these vortical convective cells truly represent three-dimensional rings. Also in 2D, the influence of the symmetry axis is already apparent as the shock is becoming elongated along it, particularly along the southern axis. At the time of explosion, the 2D and 3D structures have diverged significantly. The 2D explosion shows the characteristic "bipolar" shock structure and is dominated by about three large buoyant plumes, one each along the axes and one, somewhat smaller buoyant plume north of the equator. The 3D structure is quite different: many smaller buoyant plumes exist and are more evenly distributed in solid angle and the shock, while certainly not perfectly spherical, does not show large-amplitude low-order deformation as in the 2D case. The salient point is that Figures 4 and 5 clearly demonstrate that the surface area-to-volume ratio of the buoyant plumes is much higher in 3D than in 2D. Thus, due to the greater amount of drag relative to buoyant force the plumes will rise more slowly in 3D than in 2D and will, therefore, encourage a slower growth of the average shock radius and later explosion times.

Character of the Shock Motion
Multidimensional phenomena such as neutrino-driven convection and the SASI naturally result in aspherical shock motion in CCSNe. Linear analysis indicates that the most unstable modes of the SASI will be low-order, ≈ 1 (Foglizzo et al. 2007;Guilet & Foglizzo 2012). Three-dimensional simulations (Iwakami et al. 2008) and 2D R − φ axisymmetric simulations (Blondin & Shaw 2007) tailored to study the SASI show that the SASI develops an m = 1 "spiral" mode that can be considered the superposition of two or more out-of-phase = 1 modes (Fernández 2010). Neutrino-driven convection will excite much higher-modes of the shock motion and we have in the previous sections argued that such buoyant convection is the dominant instability that encourages expansion of the shock in the present simulations. In this section we justify this assertion by analyzing the shock motion via spherical harmonic decomposition. Our approach is comparable to that of Burrows et al. (2012). A very thorough study of the character of the shock motion in 3D is also presented by Ott et al. (2013). Their 3D GR simulations with neutrino leakage are the most physically-detailed for which the shock spherical harmonics have been investigated. They find that high-order neutrinodriven convection dominates the low-order development of the SASI, which is evident but weak in their calculations.
We can decompose some scalar quantity, X, into spherical harmonic components with coefficients where the spherical harmonics are and The spherical harmonics are computed from the associated Legendre polynomials, P m , and we use the physics convention for spherical coordinates: θ is the polar angle and φ is the azimuthal angle. In order to study the shock behavior, we decompose the shock surface R s (θ, φ) according to equation (8) and define a spherical harmonic "power": where we follow Burrows et al. (2012) in normalizing the coefficients by a factor (−1) |m| / 4π(2 + 1) such that a 00 = R s , a 11 = x s , a 1−1 = y s , and a 10 = z s . In 2D, all m = 0 components cancel so that P ( ) = a 2 . We use an accurate shock finding algorithm in FLASH to track and flag zones that are within the shock (Balsara & Spicer 1999). We show in Figure 6 the first three spherical harmonic powers of the shock surface, normalized by the average shock radius, for both 2D and 3D simulations at various neutrino Figure 5. Entropy pseudo color plots for 3D and 2D simulations with L νe,52 = 1.7 at 100 ms post-bounce (left) and time of explosion (right). For 3D, three orthogonal slice planes are shown. By 100 ms the shock structure in 3D is still very spherical and high-entropy buoyant plumes are just starting to reach the shock. In 2D at 100 ms, elongation along the symmetry axis is already evident, particularly in the southern hemisphere, and the convective structures are larger and more coherent. At explosion time, 388 ms and 821 ms in 2D and 3D, respectively, the character of the shock and the buoyant convection behind it is completely divergent between 2D and 3D. The 2D explosion is characteristically dipolar and dominated by a few, large buoyant plumes (or arguably = 1 SASI). In 3D, the shock does not show a dominant low-order shape and the convective plumes show much more small scale structure. There are also a greater number of low-entropy down flows in 3D. These fundamental differences between 2D and 3D result in convective plumes that have much greater surface area-to-volume ratios in 3D. This results in a greater drag-to-buoyant force ratio which, for buoyancy-dominated shock expansion, leads to slower shock expansion in 3D. luminosities as functions of time. We find that the amplitudes of the low-order harmonics are significantly reduced in 3D as compared to 2D, particularly prior to the onset of explosion (marked by vertical dotted lines in Figure 6). In 3D, the amplitudes are especially small prior to accretion of the Fe/Si interface at around 125 ms. When this occurs, the higher entropy in the Si shell encourages rapid shock expansion (Fig.  2) that excites substantial shock deformation. This transient high-amplitude spike fades away, however, after the entirety of the Fe/Si interface is accreted through the shock. We find a general trend that the later the explosion time, the higher the P (1, 2, 3) at explosion. For example, the 3D model with L νe,52 = 2.1 explodes with very small P (1, 2, 3) and only after explosion sets in do the amplitudes grow to larger values.
Of note, and not evident in Figure 6, is that between 50 and 100 ms we find a low-amplitude spiral mode of the shock deformation. This spiral mode is damped once neutrino-heated plumes reach and impinge upon the shock at around 100 ms (see Fig. 5). Also, we have computed the spherical harmonic powers of the shock deformation for the higher-resolution 2D simulations and find no significant difference from the fiducial resolution 2D simulations.

Explosion Indicators
A number of quantities have been identified as possible indicators or predictors of explosion in core-collapse supernovae. Nordhaus et al. (2010) argued that a higher average entropy in the gain region correlates with likelihood of explosion and found a clear dimensional hierarchy in this quantity: 3D > 2D > 1D. Hanke et al. (2012), however, suggest that there exists only a small, possibly insignificant, dependence of the gain region average entropy on dimensionality. In Figure 7 we plot the mass-averaged entropy in the gain region for our 2D and 3D simulations. We define the gain region as any part of the simulation domain that has a net positive neutrino heating as determined by the difference between equations 4 and 5. We find that there is a clear hierarchy between 2D and 3D with 3D having significantly higher average entropies than 2D. The 1D average entropy, however, does not fit the dimensional hierarchy. Our results point out that average entropy in the gain region is not a good indicator of explosion likelihood as we find that 2D explodes before 3D yet still has smaller average entropies. The average entropy in the gain region is rather a measure of the time-integrated specific energy absorbed. Thus, the lower average entropies in 2D reflect that the gain region . Normalized power in the first three spherical harmonic modes of the shock deformation in 2D (red) and 3D (blue) for three neutrino luminosities. The "power" is the squared spherical harmonic coefficient summed over all m's, so in 2D where the m's cancel P ( ) = a 2 0 . The amplitudes for 3D are generally smaller than for 2D. The data necessary to compute the shock spherical harmonics for the 3D model with L νe,52 = 1.7 prior to t = 0.3 s was, mistakenly, not saved. in 2D contains a greater amount of mass than in 3D, as shown in Figure 8. This higher gain region mass is also reflected in a higher integrated heating rate (Fig. 9). A higher net heating rate should, intuitively, result in a greater rate of shock expansion. Pejcha & Thompson (2012) suggest that instability to runaway shock expansion sets in when the maximum value of the sound speed squared to the escape velocity squared reaches a critical value, the so-called 'antesonic' condition. For an isothermal equation of state, the antesonic condition is max(c 2 s /v 2 esc ) = 0.19, where c s is the adiabatic sound speed and v esc is the escape speed. For a polytropic equation of L ν,52 = 1.7 L ν,52 = 1.9 L ν,52 = 2.1 3D 2D 1D, 2.3 Figure 8. Mass in the gain region as a function of time. Gain region mass is higher in 2D than in 3D, resulting in a higher net heating rate (Fig. 9) and a lower mass-averaged entropy (Fig. 7). state, the antesonic condition becomes max(c 2 s /v 2 esc ) = 0.19Γ, where Γ is the polytropic index. For the completely general EOS we use, we modify the polytropic antesonic condition to be max(c 2 s /v 2 esc γ c ) = 0.19, where γ c is now the varying adiabatic index given by the EOS.
In Figure 10 we plot the left hand side of equation (12)  L ν,52 = 1.7 L ν,52 = 1.9 L ν,52 = 2.1 3D 2D 1D, 2.3 Figure 9. Net neutrino heating rate in the gain region. Because of the greater gain region mass in 2D, the net heating is also higher. The net heating rises dramatically at late times (thick blue line) because the average temperature in the gain region drops, reducing the neutrino cooling rate (eq. 5). L ν,52 = 1.7 L ν,52 = 1.9 L ν,52 = 2.1 3D 2D, 0.7 km 2D, 0.5 km 1D, 2.3 Figure 10. Maximum value of the ratio of the squared adiabatic sound speed to the squared escape velocity times adiabatic index, the 'antesonic' condition. We find that the antesonic condition is a good predictor of proximity to explosion, although our data indicate that the critical value of ∼ 0.2 suggested by Pejcha & Thompson (2012) is somewhat low. We also plot the data for the high-resolution 2D simulations. The antesonic condition is a better predictor of the delayed explosion time for the high-resolution simulations than any of the other integral quantities plotted in Figs. 7 -9. We also show the critical ratio for the 1D model with L νe,52 = 2.3 (black dotted line). The 1D data look very similar to the 3D data for L νe,52 = 1.7 and peak at explosion time.
shock radius exceeds 400 km. Also, the antesonic value for our 2D simulations exceeds the critical value before the 3D simulations, echoing our result that 2D explodes before 3D.
In Figure 10 we also plot the antesonic value for 2D simulations with higher resolution. Interestingly, of the quantities shown in Figures 7 -10 the antesonic condition shows the most noticeable differences between the 0.7 km and 0.5 km resolution 2D simulations. The difference is most notable for the simulations with L νe,52 = 1.7 where we see that based on the antesonic condition we would expect the higher resolution simulation to explode later as, in fact, it does.

Non-Radial Motion and Turbulence
Non-radial motion, particularly on large scales, has been suggested as a primary factor resulting in easier explosions in L ν,52 = 1.7 L ν,52 = 1.9 L ν,52 = 2.1 Figure 11. Non-radial components of the kinetic energy in the gain region. The transverse kinetic energy is higher in 2D than in 3D, but this is again a reflection of the greater mass in the gain region for 2D. multidimensional simulations as compared with sphericallysymmetric calculations . Such motion increases the matter dwell times in the gain region and, thus, the net neutrino heating rate. A measure of non-radial motion is the gain region kinetic energy in transverse, θ/φ-direction, motion, which we plot in Figure 11. The kinetic energy of transverse motion in the gain region is seen to corollate with increased neutrino luminosity. The transverse kinetic energy is also typically higher in 2D than in 3D. Turbulence has been suggested to play an important role in the supernova shock expansion in multidimensional simulations (Murphy & Meakin 2011;Murphy et al. 2013). The nature of turbulence between 2D and 3D, however, is fundamentally different. The best-known example of this fundamental difference is the so-called 'inverse' energy cascade in 2D: energy is transported to large scales in 2D whereas in 3D energy is transported to smaller scales (Kraichnan 1967). The characteristic power-law slope of the energy cascade is -5/3 in either the spherical harmonic mode, , or wavenumber, k. Also in 2D, enstrophy, a quantity proportional to the squared vorticity, is transported to smaller scales in a so-called forward cascade with a characteristic power-law index of -3.
In order to study the character of the turbulence in our simulations we have computed the θ-direction kinetic energy spectrum in spherical harmonic basis where the energy in spherical harmonic mode is and the coefficients a m are given by equation (8) with X(θ, φ) = √ ρ(r, θ, φ)v θ (r, θ, φ). Transverse kinetic energy spectra for our multidimensional simulations are shown in Figure 12 where we restrict the evaluation to a ∼10 km shell centered at r = 150 km. We also average the spectra over 10 ms centered on the respective labeled times. The differing character of turbulence between 2D and 3D is evident. Beginning around = 10, which we identify as the driving scale of the turbulence, the 3D spectrum is consistent with a −5/3 power-law and the 2D spectra roughly follow a −3 power-law, consistent with our expectations based on turbulence theory (Kolmogorov 1941;Kraichnan 1967). Also evident is that there is much more energy at small scales (large 's) in 3D than in  Figure 12. Energy spectra of the θ-direction kinetic energy, a proxy for the turbulent kinetic energy, in spherical harmonic basis calculated in a ∼10 km wide shell centered on 150 km. We plot spectra for both 2D and 3D, and the spectra are averaged over 10 ms centered on the times indicated. For 3D we find an inertial range from about = 10 − 70 over which the spectrum is roughly consistent with a −5/3 power-law. Dissipation sets in around ∼ 70 where the spectrum falls off more steeply than −5/3 . In all cases, the 2D spectra show more energy at = 1 than the 3D spectrum, reflective of the inverse energy cascade in 2D. At around the driving scale of ∼ 10, the 2D spectra follow a −3 power-law consistent with the forward enstrophy cascade. The dissipation scale for the higher resolution 2D simulations is at noticeably larger .
2D, while 2D has more energy at large scales ( = 1 − 3). In a real sense, the different natures of these energy spectra reflect the different characteristic buoyant plume sizes seen in our simulations: plumes are smaller and more numerous in 3D as compared with 2D. This is then directly related to our argument that the smaller plumes in 3D will experience a greater drag-to-buoyant force, slowing their average ascent relative to the 2D case, resulting in a less-rapid average shock expansion in 3D. Figure 12 also shows energy spectra for 2D at different times and at different resolutions. Note that for all the spectra plotted the average shock radius is approximately the same at the time shown (250-275 km, see Fig. 2). Considering the evolution of the 2D energy spectra with time, the inverse energy cascade is evident in the increase in energy by an order of magnitude at = 1 while the spectrum in the inertial range, > 10, remains relatively unchanged in the later time. The lower-resolution (0.7 km) 2D spectrum shows two important distinctions. First, dissipation sets in at smaller , consistent with the assumption that grid spacing determines the dissipation scale. Second, at 200 ms there is substantially more energy on large scales ( = 1, 2) than for the 0.5 km resolution case.
It is of curious, perhaps surprising, significance that the energy spectra seem to match the expected power-law slopes so well. The power-law scaling in the inertial ranges of -5/3 for 3D and -3 for 2D are based on the assumption of fullydeveloped isotropic turbulence (Kolmogorov 1941;Kraichnan 1967). The turbulence in our simulations is not isotropic but convective. There is no model of buoyant convective turbulence that would suggest any other power-law scalings, but that they seem to be the same is an interesting coincidence.

RESOLUTION DEPENDENCE
As discussed in the previous sections, we find that increasing the resolution in 2D simulations results in later explosion times. This is an expected result of equation (7) that smaller buoyant plumes rise more slowly than large buoyant plumes due to the greater amount of drag force relative to buoyant force. The different convective plume size is evident by the differences in the energy spectra in Figure 12, but also by visual inspection of the 2D data. In Figure 13 we show a side-by-side comparison of the 2D simulations with different resolution at two different post-bounce times. What is evident in this comparison is that, even at 200 ms, the higher-resolution simulation shows a larger number of low-entropy down flows and deep penetration of these down flows. The larger number of down flows breaks up the rising buoyant plumes while their deeper penetration further increases the surface are of the plumes and, therefore, the drag force felt by the plume. This is especially evident at explosion time (right half of Fig. 13).
This is simple to understand by considering the primary instabilities that will contribute to break-up of the buoyant plumes. The plumes are inherently Rayleigh-Taylor unstable 2 and will also develop parasitic Kelvin-Helmholtz instabilities as they rise through the post-shock accretion flow. In the inviscid limit, the growth rate for both of these instabilities is inversely proportional to the grid scale (Chandrasekhar 1961;Youngs 1984), thus the higher-resolution simulation experiences faster growth of these instabilities that impede the rising plumes, slowing overall shock expansion. Hanke et al. (2012) report that increasing the resolution in 3D results in delayed explosions while increasing resolution in 2D results in earlier explosions. Their 2D resolution study results are contradictory to our simulation results. Hanke et al.'s conclusion that increasing resolution in 2D results in earlier explosions is based on their simulations using only 400 radial zones. When considering their results with 600 or 800 radial zones, however, the dependence on resolution is less clear; for some sets of 2D models with more than 400 radial zones increasing the angular resolution results in later explosions, as we have found. Hanke et al. suggest that this behavior with increased radial resolution is due to a artificial density peak in the cooling region that grows with increasing radial resolution, enhancing total cooling. We do not see such an artificial density peak in our simulations. It is also worth noting that our fiducial "angular" resolution of 0. • 54 is substantially finer than the fiducial resolution considered by Hanke et al. (3. • ) and comparable to the finest resolution they use in 2D (0. • 5). Nordhaus et al. (2010) and Dolence et al. (2013) use a fiducial resolution greater than ours (0.5 km) but do not conduct resolution studies so it is unknown how their results would vary with resolution.
In order to further investigate the resolution dependence of our results, and whether they are "converged" with respect to grid resolution, we have carried out additional simulations with varied resolution. In 2D we further increased the effective angular resolution from 0. • 4 to 0. • 27 while keeping the maximum resolution fixed at 0.5 km. This moved the transition from ∆x i = 0.5 km to ∆x i = 1.0 km from around a radius of 100 km to 200 km. In 3D, due to the great expense of higher resolution simulations, we reduced the angular resolution from 0. • 54 to 1. • 07 while keeping the maximum resolution fixed at 0.7 km. This moved the resolution transition from 0.5 km to 1.0 km from a radius of 100 km to around 50 km. The shock radius curves for calculations of various resolutions are shown in Figure 14. All these simulations use L νe,52 = 1.7. We see The labels indicate maximum resolution as well as effective angular resolution for each simulation shown. For 2D convergence with respect to shock radius history is obtained at a resolution of 0. • 40 while in 3D such convergence is already achieved at resolutions of 1. • 07.
that in 2D resolutions of 0. • 4 and 0. • 27 result in very similar shock expansion histories and explosion times, indicating that the 2D results are "converged" at an angular resolution of 0. • 4. In 3D, the 1. • 07 resolution simulation also results in a similar shock history and explosion time as the 0. • 54 simulation. This is encouraging and interesting because it implies that the rate of convergence is faster in 3D than in 2D. Thus our fiducial 3D resolution is sufficient to obtain convergence with respect to explosion time. We note that convergence is a difficult notion in astrophysical simulations that involve instabilities and turbulence such as ours. Our simulations are not converged in the true sense that changing the resolution would yield indistinguishable results. For this to be the case, we would need to reach numerical Reynolds numbers comparable to the physical Reynolds numbers of the flow. We are no where near achieving this criterion in our calculations. We have, however, achieved convergence with resolution of the result we seek: the explosion times.

DISCUSSION
The results of previous studies similar to ours have been contradictory. The Princeton group has found that 3D simulations explode sooner than 2D (Nordhaus et al. 2010;Dolence et al. 2013) while the Garching group finds little difference in the explosion times between 3D and 2D . The initial aspiration of our study was to cast a tie-breaking vote in this discussion but instead it has raised new questions by yielding a third result: our 3D simulations explode later than our 2D simulations. The source of the disparity in the results from the different groups is still unclear. There are significant differences in numerical approach between the Princton and Garching groups. The Princeton group uses the new code CAS-TRO ) which implements a directionallyunsplit piecewise parabolic method with the Riemann solver of Bell et al. (1989). CASTRO uses patch-based adaptive mesh refinement (c.f., Berger & Oliger 1984;Berger & Colella 1989) as provided by the Boxlib framework. The Princeton simulations are run in 1D spherical, 2D cylindrical, and 3D Cartesian coordinates. The Garching group uses the venerable hydrodynamics code PROMETHEUS (Fryxell et al. 1991) which implements directionally-split PPM with a 'two shock' Riemann solver (Colella & Glaz 1985) in smooth flow and the HLLE Riemann solver in shocks. PROMETHEUS does not use AMR but instead relies on 1D, 2D, and 3D spherical coordinates with non-equidistant radial spacing. Though perhaps the most important differences between the approaches of the Princeton and Garching groups are those dealing with approximations to the neutrino physics. Nordhaus et al. (2010) and Dolence et al. (2013) follow closely the approach of Murphy & Burrows (2008): deleptonization is approximated using the density-dependent parameterization of Liebendörfer (2005), both pre-and post-bounce, and post-bounce neutrino heating and cooling are considered locally based on rates derived by Janka (2001). Heating and cooling are shut off at high density by adding a e −τν term to the rates. Hanke et al. (2012) follow a very similar approach except for the following. They follow collapse and bounce to 15 ms post-bounce with full neutrino transport in 1D and do not employ the parameterization of Liebendörfer (2005) at all. The subsequent multidimensional evolution uses the same heating and cooling rates as in Nordhaus et al. (2010) but the neutrino optical depth is computed by integrating the appropriate opacity in radius whereas the Princeton group uses a density-dependent parameterization for the optical depth (J. Murphy, private communication, see also Couch 2013). In their multidimensional post-bounce calculations, the Garching group follows the evolution of the inner core in 1D to avoid issues with grid convergence.
Our approach is very similar to those of Nordhaus et al. and 2D cylindrical,and 3D Cartesian with AMR. We use the Liebendorfer parameterized deleptonization pre-and postbounce and approximate the neutrino optical depth with a density-dependent piecewise fit, although we scale optical depth differently than the Princeton group resulting in a different normalization of the critical curves (see also the discussion in Hanke et al.). Notable differences in our scheme and that of the Princeton group are we use an oct-tree, block-structured AMR package that does not sub-cycle in time, i.e. all refinement levels advance with the same time step size. CASTRO uses adaptive time refinement allowing coarser resolution levels to advance with larger time steps. And we use a "hybrid" Riemann solver that uses the HLLC solver in smooth flows and the HLLE solver in shocks. Another possibly significant difference amongst all three groups is in the EOS implementation. While each group is using the Shen et al. high-density baryonic EOS for the referenced simulations, the construction of the tables actually used is no doubt somewhat different. There could also be important differences in the details of the implementations of the monopole gravity solver amongst the different codes. What is clearly mandated by the disparity in the results is a rigorous code-to-code comparison.
Our results, particularly in 3D, support the conjecture that, for this progenitor and treatment of neutrino physics, the SASI is subdominant to neutrino-driven convection in advancing the average shock radius (Burrows et al. 2012;Murphy et al. 2013). We are, however, very hesitant to extrapolate this conclusion to all possible progenitors. For instance, Müller et al. (2012a) show clear evidence for a strong SASI in the 2D explosion of a 27 M progenitor using conformally-flat general relativistic dynamics and full neutrino transport. Although recently Ott et al. (2013) have simulated the same 27 M progenitor with full 3D GR and a multispecies neutrino leakage scheme and found that neutrino-driven convection becomes the dominant instability in exploding models.
In the scenario that neutrino-driven convection dominates, the shock expansion can be roughly described by buoyant plumes plumes rising against the drag force exerted by the post-shock accretion flow . The buoyant force exerted by the neutrino radiation on a plume will be proportional to the plumes subtended solid angle times the optical depth of the plume multiplied by the neutrino luminosity, i.e., the plume's volume. The drag force exerted by the accretion down flows will be proportional to the plume's surface area. A natural result of this model is that a single plume with volume V will rise more quickly than two plumes with volumes V /2 due to the greater amount of total surface area, and hence drag force, for the two plumes [see equation (7)]. As discussed at length in Section 4, 3D simulations develop smaller, more numerous buoyant plumes than 2D simulations. Thus the balance between buoyancy and drag explains why the 3D models explode later than the 2D models. This picture also accounts for the 2D resolution dependence as higher-resolution will result in plumes that have greater surface area-to-volume ratios. 7. CONCLUSIONS We have conducted a 1D, 2D, and 3D parameter study of neutrino-driven core-collapse supernovae designed to explore the difference that 3D makes on the explosion characteristics, particularly time of explosion relative to 1D and 2D. We find, as have a number of previous studies, that the so-called 'critical curve' in the neutrino luminosity-mass accretion rate at explosion time plane is significantly lowered in multiple dimensions relative to the spherical-symmetric case. A novel result of our study is we find that 3D explosions occur later than 2D explosions at the same neutrino luminosity, i.e., the 3D critical curve is higher than the 2D critical curve. We find that the 2D results are resolution-dependent: increasing the resolution in 2D delays explosion pushing the high-resolution 2D critical curve very near to the 3D critical curve (Fig. 3).
We suggest that our results can be explained by the competition between buoyancy and drag. The shock expansion in our simulations is dominated by the action of neutrino-driven buoyant convection (see also Burrows et al. 2012;Murphy et al. 2013). In this case, the shock expansion can be fit by the motion of a buoyant plume rising through the shocked accretion flow . The buoyant force acting on the plume, provided by the absorbed neutrino radiation, is proportional to a plume's volume whereas the drag force resulting from the downward-flowing accretion is proportional to a plume's area. Thus, a plume's ascension velocity increases with increasing volume-to-surface area ratio. Rapid shock expansion, therefore, is best abetted by plumes that subtend large solid angles. We have shown that 3D simulations naturally result in many more, smaller-solid angle plumes than comparable 2D simulations. We posit that this, then, is why our 3D simulations explode later than our 2D simulations and why higher-resolution 2D simulations also result in later explosions.
We examined several differences between our 2D and 3D simulations. In Section 4.2 we explored the character of the shock motion by calculating the first few spherical harmonic powers of the shock deformation in 2D and 3D. We find that the amplitudes of the low-order, = 1, 2, modes of the shock motion that are often associated with the SASI are much reduced in 3D relative to 2D. This result is in qualitative agreement with Burrows et al. (2012).
When considering possible indicators of explosion, we find that the "antesonic" condition of Pejcha & Thompson (2012) correlates well with proximity to explosion. We find that the mass-averaged entropy in the gain region is typically higher in 3D than in 2D, but our results indicate that is is not a good indicator of proximity to explosion, as it was suggested to be by Nordhaus et al. (2010). Instead this higher average entropy, which is proportional to the time-integrated neutrino heating per mass, only reflects that there is less mass in the gain region in 3D as compared to 2D. The greater gain region mass in 2D also results in a greater net heating rate, which also aids explosion relative to 3D.
The character of the non-radial motion and turbulence between our 2D and 3D simulations is distinct. In 3D, the forward energy cascade transports energy to smaller and smaller scales until dissipation do to finite grid size sets-in. In 2D, the forward enstrophy cascade results in a much less efficient transport of energy to small scales while the inverse energy cascade actually transports energy from the driving scales to larger scales. This behavior is conducive to explosion as it encourages both the growth of the low-order SASI and the development of large buoyant plumes.

COUCH
The cause of the differences in results among the groups attempting similar multidimensional parameter studies (e.g., Nordhaus et al. 2010;Hanke et al. 2012;Dolence et al. 2013, this work) remain to be explained. The simplified physics employed by such studies would make a detailed code-to-code comparison effort far more straight-forward than a comparison that included neutrino transport. Such a code-to-code comparison will hopefully be seen as a high priority and be accomplished in the near future.
The approximations we employ are admittedly crude. They have tremendous merit, however, in that they facilitate 3D parameter studies of CCSNe. While certain conclusions about CCSNe based on studies such as the present should be made with caution, the 3D results appearing in the literature todate make one undeniable point: 3D core-collapse supernovae are fundamentally and dramatically different than 2D corecollapse supernovae. The absence of forced-symmetry in 3D makes an enormous impact on the character of the shock motion and the development of neutrino-driven convection and turbulence. Our results indicate, however, that 3D alone may not be the key to successful, robust explosions in more realistic simulations. Still, the enormous difference between 2D and 3D CCSN simulations emphasizes and underlines the need for fully 3D simulations by this community.
The author thanks Evan O'Connor, Chris Daley, Carlo Graziani, Cal Jordan, Jeremiah Murphy, and Christian Ott for helpful and insightful conversations. The author is especially grateful to Todd Thompson for very valuable comments on the manuscript. 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