NEUTRINO-DRIVEN TURBULENT CONVECTION AND STANDING ACCRETION SHOCK INSTABILITY IN THREE-DIMENSIONAL CORE-COLLAPSE SUPERNOVAE

The top panel of Figure 2 shows the time evolution of the angle-averaged shock radius R_shock,avg in our three baseline-resolution simulations s27MRf_heat1.05, s27MRf_heat0.95, and s27MRf_heat0.8 with strong, medium, and weak neutrino heating, respectively. We show only the part of the evolution tracked in 3D. At early times (t−t_bounce ≲ 50–60 ms) the shock undergoes some transient oscillations as it relaxes on the 3D grid, which is reflected in R_shock,avg of all models.

Zoom In Zoom Out Reset image size

The average shock radius in model s27MRf_heat1.05 grows secularly from $\sim 105$ to 125 km in the first 90 ms of postbounce evolution. Subsequently, the shock expansion accelerates and R_shock,avg reaches ∼195 km by ∼130 ms after the bounce, which is when we stop following this model's evolution. The maximum shock radius at this time is ∼220 km. The expansion has become dynamical and is most likely transitioning to explosion. In contrast, models s27MRf_heat0.95 and s27MRf_heat0.8 do not show any sign of explosion within the simulated time. The average shock radius in these models decreases gradually until ∼160 ms after bounce, reaching ∼97 and ∼75 km, respectively. At this point, the silicon shell reaches the shock front, leading to a sudden decrease of the accretion rate (see the accretion rate shown in the top panel of Figure 2) and thus of the ram pressure experienced by the shock. This leads to a transient expansion of the shock by ∼10 km within ∼15 ms, after which it starts retreating again in both models and continues to do so until the end of our simulations. Due to the weaker heating in model s27MRf_heat0.8, R_shock,avg always remains smaller than in model s27MRf_heat0.95, but has the same qualitative evolution.

The shock radius evolution shown in Figure 2 can be directly linked to the strength of neutrino heating. We quantify the latter using a set of metrics shown in Figure 3: the integral net neutrino heating rate Q_net, the heating efficiency ($\eta ={Q}_{\mathrm{net}}{({L}_{{\nu }_{e}}+{L}_{{\bar{\nu }}_{e}})}^{-1}$, where ${L}_{{\nu }_{e}}$ and ${L}_{{\bar{\nu }}_{e}}$ are the electron neutrino and anti-electron neutrino luminosities incident from below the inner boundary of the gain region),¹¹ and the mass M_gain in the gain region. The oscillations in these quantities at early times are a combined artifact of the leakage/heating scheme, which is unreliable in highly dynamical situations, and of the shock settling on the 3D grid. Note that the outgoing luminosities are only mildly affected and the main effect comes from variations in the mean neutrino energies; see Ott et al. (2013). Similar features are present in the leakage simulations of Couch & O'Connor (2014). As expected, the larger the f_heat (see Equation (1)), the larger the integral net heating, heating efficiency, and the mass in the gain region. Note, however, how strong this relationship is: an increase of f_heat from 0.95 to 1.05 (∼10.5%) results in approximately twice as high Q_net, η, and M_gain around 50–100 ms after bounce. This is a consequence of the fact that more intense neutrino heating extends the region of net absorption to smaller radii. It also increases the thermal pressure and the vigor of turblence (and thus the effective turbulent ram pressure; Couch & Ott 2015) throughout this region. This, in turn, pushes the shock out, further increasing the volume of the gain region and leading to more net neutrino energy absorption. This nonlinear feedback shows just how extremely sensitive core-collapse supernovae near the critical line between explosion and no explosion are to the details of neutrino transport and neutrino–matter coupling.

Zoom In Zoom Out Reset image size

The general trends in neutrino heating with f_heat described in the above hold throughout the postbounce phase. However, as the shock radii in models s27MRf_heat0.95 and s27MRf_heat0.8 recede and their gain regions shrink, the values of their neutrino heating variables shown in Figure 3 approach each other. The sudden reduction of the ram pressure at the silicon interface, which has a significant effect on the shock radius (Figure 2), is barely noticeable in the neutrino heating.

We plot the average mass-weighted specific entropy in the gain region (s_gain) on the right ordinate of the bottom panel of Figure 3. In agreement with what was found in previous work (e.g., Hanke et al. 2012; Couch 2013; Dolence et al. 2013; Ott et al. 2013), s_gain is largely independent of the shock radius and the strength of neutrino heating in the postbounce pre-explosion phase simulated here. We attribute this to two competing effects that affect the averaged quantity s_gain: while strong neutrino heating (larger f_heat in our simulations) leads to locally higher specific entropy in the region of strongest heating, this is compensated for by the overall larger volume (and mass) of the gain region, which includes material of lower specific entropy that contributes to the average.

After considering the above range of indicative angle-averaged and/or volume-averaged quantities, it is now useful to study deviations from averaged dynamics. The center and bottom panel of Figure 2 depict the normalized rms angular deviation of the shock radius from its mean (R_{shock, avg.}), σ_shock, defined as

$$\begin{eqnarray}&&{\sigma }_{\mathrm{shock}}=\displaystyle \frac{1}{{R}_{\mathrm{shock},\mathrm{avg}.}}\sqrt{\displaystyle \frac{1}{4\pi }{\displaystyle \int }_{4\pi }{\left[{R}_{\mathrm{shock}}(\theta ,\phi )-{R}_{\mathrm{shock},\mathrm{avg}}\right]}^{2}d{\rm{\Omega }}},\end{eqnarray} \tag{ 2 }$$

and the ratio of maximum and mimimum shock radius R_shock,max./R_{shock, min.}, respectively. Both diagnostics yield qualitatively similar results, but the latter is more sensitive to small local variations. In the initial settling phase on the 3D grid, all three simulations shown in Figure 2 exhibit nearly identical shock deviations from sphericity. These are due to moderate-amplitude cubed (ℓ = 4–symmetric) shock oscillations as the models relax from spherical geometry to our 3D Cartesian grid. Subsequently, the shock deviations begin to differ between models. Model s27MRf_heat1.05 (strong neutrino heating) shows more or less steadily increasing asphericity as its shock gradually expands and develops large-scale deviations in the maximum shock radius R_shock,max. from R_shock,avg. driven by expanding localized high-entropy bubbles. However, the overall asymmetry, expressed by σ_shock in Figure 2, is relatively small, due to the strong neutrino heating that leads to a rather global shock expansion.

The shock asphericity in models s27MRf_heat0.95 and s27MRf_heat0.8 exhibits strong oscillations with clear (if temporally varying) periodicity—a tell-tale sign of active SASI. In model s27MRf_heat0.95 (moderate neutrino heating), the oscillations set in around ∼105 ms after the bounce while in model s27MRf_heat0.8 (weak neutrino heating), they are already present at ∼80 ms. In both models, the period of the oscillations changes when the silicon interface reaches the shock front around 160 ms after the bounce. We will analyze SASI in these models in more detail in Section 4.2.

Whenever our simulations experience strong neutrino heating, i.e., η ≳ 0.05 and Q_heat ≳10⁵² erg s⁻¹, we find neutrino-driven convection in the postshock region. This is quantified by the top panel of Figure 4, which shows the buoyant mass in the gain region, ${M}_{\mathrm{gain},{\upsilon }_{r}\gt 0}$, which we define as the mass of material with positive radial velocity. ${M}_{\mathrm{gain},{\upsilon }_{{\rm{r}}}\gt 0}$ correlates strongly with η and Q_net. Phases of strong neutrino heating (see Figure 3) correspond to strong neutrino-driven convection. Model s27MRf_heat1.05 undergoes convection throughout its postbounce evolution, while model s27MRf_heat0.95 exhibits strong neutrino-driven convection only until ∼110 ms after the bounce. Convective activity is clearly visible in the 2D (x–z plane) entropy color maps of models s27MRf_heat1.05 and s27MRf_heat0.95 at various postbounce times in Figure 5. We will further analyze neutrino-driven convection in our models in Section 4.3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

There are three defining characteristics of SASI: (1) low-(ℓ,m) oscillations of the shock front (e.g., Iwakami et al. 2008), (2) exponential growth of the (spherical harmonics) mode amplitudes in the linear phase (e.g., Blondin et al. 2003), and (3) saturation of the amplitudes once they reach the nonlinear phase (e.g., Guilet et al. 2010). In order to identify these features in our simulations, we decompose the shock front R_shock(θ, ϕ) into spherical harmonics:

$$\begin{eqnarray}&&{a}_{{\ell }m}=\displaystyle \frac{{(-1)}^{| m| }}{\sqrt{4\pi (2{\ell }+1)}}{\displaystyle \int }_{4\pi }{R}_{\mathrm{shock}}(\theta ,\phi ){Y}_{{\ell }}^{m}(\theta ,\phi )d{\rm{\Omega }},\end{eqnarray} \tag{ 3 }$$

where ${Y}_{{\ell }}^{m}$ are the standard real spherical harmonics (e.g., Boyd 2001). We employ the normalization convention used in Burrows et al. (2012), in which a₀₀ corresponds to the average shock radius, while a₁₁, a₁₀, and a_{1 − 1} correspond to the average x, z, and y Cartesian coordinates of the shock front, respectively.

Figure 6 depicts the normalized mode amplitudes ${a}_{{lm}}\cdot {a}_{00}^{-1}$ for ${\ell }=1$ (left panels) and ℓ = 2 (right panels) for the three previously introduced models with strong, medium, and weak neutrino heating. The mode amplitudes grow gradually in magnitude over time, which is reflected in the evolution of the angular deviation σ of the shock radius in Figure 2. The relative asphericity of the shock is increasing with time in all models.

Zoom In Zoom Out Reset image size

In model s27MRf_heat1.05, the ℓ = 1 mode amplitudes grow quickly at 50–80 ms after the bounce, exhibit ∼three periodic modulations with a period of ∼20 ms at nearly saturated magnitude, and then begin to increase to larger values. The ℓ = 2 modes start growing earlier, but show less clear periodicity. The evolution of the ℓ = 1 and ℓ = 2 modes suggest that some form of SASI is present in model s27MRf_heat1.05, but a look at the top row of specific entropy slices in Figure 5 reveals that violent neutrino-driven convection is active, fully developed, and driving the local deviation from spherical symmetry at late times in this model.

Models s27MRf_heat0.95 and s27MRf_heat0.8 exhibit strong SASI oscillations in their ℓ = 1 and ℓ = 2 modes that last for many cycles. The ℓ = 2 modes actually start growing first and the initial growth of all modes exhibits an exponential character until they reach saturation on a timescale of ∼50 ms. The oscillation period is ∼10 and ∼6 ms for ℓ = 1 and ℓ = 2, respectively.

From the x − z specific entropy slices shown in Figure 5, one notes that at ∼80 ms after the bounce there are signs of convection in model s27MRf_heat0.95, but no convective plumes are visible in model s27MRf_heat0.8 with the weakest neutrino heating. The entropy slices at 150 ms after the bounce show large shock deformations with ℓ = 2 symmetry and no clearly convective features in either model. Interestingly, in both models the ℓ = 2 modes get damped and overtaken by ℓ= 1 oscillations at ∼170 ms after the bounce (see Figure 6), when the silicon interface reaches the shock, leading to transient shock expansion. Accordingly, the late-time entropy slices of these models in Figure 5 exhibit predominantly ℓ = 1 asymmetry.

Although the accretion of the silicon interface damps the initially dominant ℓ = 2 modes significantly (see Figure 6), they again, but only episodically, reach large amplitudes at later postbounce times. This is uncharacteristic of linear growth of physical models and may possibly be due to nonlinear interactions with the then-dominant ℓ = 1 modes.

The ℓ = 1, m = {−1, 0, 1} modes shown in Figure 6 have different phases with respect to each other. This is suggestive of "spiral" SASI oscillations as identified, e.g., by Blondin & Mezzacappa (2007), Fernández (2010), and Iwakami et al. (2014). We analyze the vector

$$\begin{eqnarray}&&{{\boldsymbol{a}}}_{1}=\left({a}_{1-1},{a}_{10},{a}_{11}\right),\end{eqnarray} \tag{ 4 }$$

which gives the direction and magnitude of the ℓ = 1 shock deformation with respect to the center of the protoneutron star (Hanke et al. 2013). We visualize the time evolution of ${\boldsymbol{}}{a}_{1}/{a}_{00}$ with a line in 3D space in the top and bottom panels of Figure 7 for models s27MRf_heat0.95 and s27MRf_heat0.8, respectively. Each point on the graph is color coded according to postbounce time t − t_bounce. During the early postbounce evolution, $| {\boldsymbol{}}{a}_{1}| /{a}_{00}$ is small and does not exhibit any clear rotational patterns in either of the models. After the silicon interface has accreted through the shock, the ℓ = 1 modes reach large amplitudes. It is then (orange–red colors in Figure 7) that ${{\boldsymbol{a}}}_{1}/{a}_{00}$ clearly describes several complete spiral cycles in both models. This confirms the spiral nature of the late ℓ = 1 SASI, which is qualitatively very similar to what Hanke et al. (2013) and Couch & O'Connor (2014) found in their 3D simulations of the same progenitor.

Zoom In Zoom Out Reset image size

It is interesting to ask why we observe an early growth of the ℓ= 2 SASI mode in our simulations while ℓ = 1 is usually identified to be the most unstable SASI mode. We speculate, fueled by T. Foglizzo (2015, private communication), that one possible explanation may be related to the trend found by Foglizzo et al. (2007) that higher values of ℓ are favored when the shock radius is small. Just before the accretion of the silicon shell interface, the average shock radius R_shock,avg in models s27MRf_heat0.95 and s27MRf_heat0.8 is as small as 97 and 75 km, respectively. The reduction in ram pressure at the silicon interface lets the shock jump outward, possibly creating a situation more favorable for ℓ = 1 oscillations than before. This might be the reason for the sudden damping of the ℓ = 2 modes and the development of the ℓ = 1 oscillations.

In the simulations of the same 27 M_⊙ progenitor of Ott et al. (2013), Couch & O'Connor (2014), and Hanke et al. (2013), the ℓ = 1 modes reach large amplitudes before the accretion of the silicon interface. It dominates over ℓ = 2 at least in the early evolution in Couch & O'Connor (2014) and Ott et al. (2013; Hanke et al. 2013 do not provide ℓ = 2 amplitudes). In these simulations, the average shock radius is nearly always above 100 km in the early postbounce phase. It drops below this value early on in our present simulations with weak (s27MRf_heat0.8) and moderate (s27MRf_heat0.95) neutrino heating. Following the above argument, this may explain why only our simulations exhibit an initially predominantly ℓ = 2 SASI.

It is worth mentioning that, to the best of our knowledge, strong excitation of predominantly ℓ = 2 modes in the 3D case was observed only in the work of Takiwaki et al. (2012), who studied the 3D postbounce hydrodynamics in a 11.2 M_⊙ progenitor. However, in their simulation, this mode undergoes only 2–3 oscillations during the simulated time, whereas in our case, we observe ∼30 oscillation cycles before ℓ = 2-dominated dynamics ceases. The ℓ = 2 modes also reach large amplitudes in the 3D simulations of Iwakami et al. (2008), Ott et al. (2013), and Couch & O'Connor (2014), but their amplitudes generally do not exceed those of the ℓ = 1 modes.

Neutrino heating in the gain region establishes a negative radial entropy gradient (e.g., Herant et al. 1992) and thus can drive convection. In stable stars, convection occurs on a stationary background. Not so in the postshock region of a core-collapse supernova: material accreting through the stalled shock is advecting toward the protoneutron star with velocities up to a few percent of the speed of light. In order for convection to fully develop, convective plumes must not only be buoyant with respect to the rest frame of the background flow, but must be able to rise in the laboratory (coordinate) frame against the background advection stream.

Depending on the accretion rate (determined by the progenitor structure; e.g., O'Connor & Ott 2011), strength of neutrino heating (i.e., steepness of the entropy gradient), and initial size of the perturbations entering through the shock from which buoyant plumes can grow, one can identify three different regimes of convection: (1) dominance of advection, where plumes do not even become buoyant in the rest frame of the accretion flow; (2) plumes are buoyant in the rest frame of the accretion flow, but are still advected out of the gain region into the convectively stable cooling layer; (3) plumes are fully buoyant and rise against the accretion flow. As we shall see, our simulations cover all three of these regimes.

We analyze buoyant convection in our simulations with the Ledoux criterion (Ledoux 1947) and express its compositional dependence in terms of the lepton fraction Y_l:

$$\begin{eqnarray}&&{C}_{{\rm{L}}}=-{\left(\displaystyle \frac{\partial \rho }{\partial P}\right)}_{s,{Y}_{l}}\left[{\left(\displaystyle \frac{\partial P}{\partial {s}_{{\rm{T}}}}\right)}_{\rho ,{Y}_{l}}\left(\displaystyle \frac{{{ds}}_{{\rm{T}}}}{{dr}}\right)+{\left(\displaystyle \frac{\partial P}{\partial {Y}_{l}}\right)}_{\rho ,{s}_{{\rm{T}}}}\left(\displaystyle \frac{{{dY}}_{l}}{{dr}}\right)\right],\end{eqnarray} \tag{ 5 }$$

where ${s}_{{\rm{T}}}=s+{s}_{\nu }$ is the sum of entropies of the matter s and neutrino field ${s}_{\nu }$, while ${Y}_{l}={Y}_{e}+{Y}_{{\nu }_{e}}-{Y}_{\bar{{\nu }_{e}}}$ is the lepton fraction. Since our leakage/heating scheme does not track the neutrino distribution function, we set Y_l = Y_e and ${s}_{\nu }=0$ in Equation (5). This is a very good approximation in the gain region, where neutrinos are almost free streaming, but is less accurate in the protoneutron star where neutrinos are trapped at densities above $\sim {10}^{12}\;{\rm{g}}\;{\mathrm{cm}}^{-3}$. A fluid parcel is convectively stable if ${C}_{{\rm{L}}}\leqslant 0$ and unstable otherwise. In the latter case, the linear growth time of small perturbations to buoyant plumes is given, approximately, by the inverse of the Brunt–Väisälä (BV) frequency,

$$\begin{eqnarray}&&{\omega }_{\mathrm{BV}}=\mathrm{sgn}\left({C}_{{\rm{L}}}\right)\sqrt{\displaystyle \frac{\left|{C}_{{\rm{L}}}\right|g}{\rho }},\end{eqnarray} \tag{ 6 }$$

so ${\omega }_{\mathrm{BV}}\gt 0$ implies instability. Here, g is the local free-fall acceleration, which we approximate as $-{GM}(r){r}^{-2}$ in our postprocessing analysis, where M(r) is the mass enclosed within radius r. A similar approach was used in, e.g., Buras et al. (2006b), Takiwaki et al. (2012), and Ott et al. (2013).

In addition, we compute the Foglizzo χ parameter (Foglizzo et al. 2006),

$$\begin{eqnarray}&&\chi ={\displaystyle \int }_{{R}_{\mathrm{gain}}}^{{R}_{\mathrm{shock}}}\displaystyle \frac{{\omega }_{\mathrm{BV}}}{| {\upsilon }_{{\rm{r}}}| }{dr},\end{eqnarray} \tag{ 7 }$$

where ${\upsilon }_{{\rm{r}}}$ is the radial velocity in the gain region. χ can be interpreted as the ratio of the advection timescale to an average timescale of convective growth. Any small linear seed perturbation (coming, e.g., from turbulent convection in nuclear burning shells; e.g., Arnett & Meakin 2011; Couch & Ott 2013, 2015) accreting through the shock can at most grow by a factor of $\sim \mathrm{exp}(\chi )$ during its advection through the gain region. For such linear-scale perturbations, Foglizzo et al. (2006) found that $\chi \gtrsim 3$ is necessary for convection to develop in the gain region. The situation is different for large seed perturbations for which the time integral of buoyant acceleration is comparable to the advection velocity (Scheck et al. 2008). In this case, a seed perturbation may develop into a buoyant plume and stay in the gain region instead of being advected out. The results of Scheck et al. (2008) suggest that seed perturbations of $\sim 1\%$ may be sufficient to allow fully developed convection even when $\chi \lt 3$. Fernández et al. (2014) pointed out that χ is quite sensitive to the way it is calculated. We follow the recent works of Ott et al. (2013), Couch & O'Connor (2014), and Hanke et al. (2013), who all used instantaneous angle-averaged quantities to compute χ via Equation (7).

If convection develops (either in regime 2 or 3, which we introduced earlier in this section), its vigor can be measured using the anisotropic velocity ${\upsilon }_{\mathrm{aniso}}$ defined as (Takiwaki et al. 2012)

$$\begin{eqnarray}&&{\upsilon }_{\mathrm{aniso}}(r)=\sqrt{\displaystyle \frac{{\langle \rho \left[{\left({\upsilon }_{{\rm{r}}}-\langle {\upsilon }_{{\rm{r}}}{\rangle }_{4\pi }\right)}^{2}+{\upsilon }_{\theta }^{2}+{\upsilon }_{\varphi }^{2}\right]\rangle }_{4\pi }}{\langle \rho {\rangle }_{4\pi }}},\end{eqnarray} \tag{ 8 }$$

where ${\langle .\rangle }_{4\pi }$ denotes an angular average at a fixed radius r. ${\upsilon }_{\mathrm{aniso}}$ measures the magnitude of the velocity component that is not associated with a purely spherically symmetric radial background flow. ${\upsilon }_{\mathrm{aniso}}$ is high in regions of large angular variations in ${\upsilon }_{{\rm{r}}}$ and large nonradial velocities ${\upsilon }_{\theta }$ and ${\upsilon }_{\varphi }$.

Convective activity in our simulations can be diagnosed via Figure 4 (showing the amount of buoyant mass and Foglizzo χ), Figure 5 (showing color maps of 2D $x-z$ entropy slices at various postbounce times), and Figure 8 (showing the evolution of radial profiles of the angle-averaged BV frequency ${\omega }_{\mathrm{BV}}$ and ${\upsilon }_{\mathrm{aniso}}$).

Zoom In Zoom Out Reset image size

In all models, within milliseconds of the bounce, a convectively unstable region with a steep negative entropy gradient develops inside the radial shell ranging from $\sim 25$ to $\sim 40\;\mathrm{km}$ due to the propagation of the gradually weakening shock. In our simulations, this phase occurs already during the 1D evolution with GR1D (not shown here). This leads to the development of strong prompt convection within $\sim 20\;\mathrm{ms}$ after the start of the 3D simulations, as is evident from the ${\upsilon }_{\mathrm{aniso}}$ profiles shown in Figure 8. The χ parameter (Figure 4) is generally $\lt 3$ in all models, but prompt convection develops nevertheless from numerical perturbations, which are $\gtrsim 1\%$ at the time the profile is mapped from one dimension to three dimensions and settles on the 3D grid (see the discussion in Ott et al. 2013 about perturbations from the Cartesian computational grid). Prompt convection smoothes out the negative entropy gradient on a timescale of 5–10 ms, leading to a rapid weakening and then to complete disappearance of convection. The latter is most apparent from the dramatic decrease in buoyant mass shown in the top panel of Figure 4.

Deleptonization at the edge of the protoneutron star creates a negative lepton gradient within 30–40 km. It drives convection in the protoneutron star, setting in at 35–50 ms after the bounce (Figure 8). Protoneutron star convection (albeit modeled only schematically, given the limitations of our neutrino treatment; see Section 2) is similar in all models, since it is independent of neutrino heating in the gain region.

In model $s27\mathrm{MR}{f}_{\mathrm{heat}}1.05$, neutrino heating creates a negative entropy gradient in the region between $\sim 80\;\mathrm{km}$ and the shock, leading to a convectively unstable layer, as apparent from the upper left panel of Figure 8. This triggers and sustains convection in the postshock region starting at $t-{t}_{{\rm{b}}}\sim 50\;\mathrm{ms}$, at this early time aided by additional entropy perturbations coming from variations in the shock radius. The amount of buoyant mass (top panel of Figure 4) has a local maximum when convection first starts and exceeds this maximum only once the explosion begins to develop in model $s27\mathrm{MR}{f}_{\mathrm{heat}}1.05$. The Foglizzo χ parameter shown (bottom panel of Figure 4) suggests that much of the convection, while clearly visible in the entropy slices of this model shown in Figure 5, is not fully bouyant in the coordinate frame (regime 2). Only at $t\gtrsim 100\;\mathrm{ms}$ after the bounce does χ grow beyond the linear-theory threshold value of 3 and the amount of buoyant mass increases, indicating that convection is now fully buoyant and convective plumes begin to push out the stalled shock, driving both its expansion and asymmetry (regime 3). These general trends agree well with what was found by Burrows et al. (2012), Couch (2013), Dolence et al. (2013), Ott et al. (2013), and Couch & O'Connor (2014) for 3D simulations with strong neutrino heating that yielded explosions.

In model $s27\mathrm{MR}{f}_{\mathrm{heat}}0.95$ with moderate neutrino heating, the buoyant mass peaks when neutrino-driven convection first develops and then gradually declines with time. While we see clear signs of convection in the entropy snapshot at $\sim 80\;\mathrm{ms}$ after the bounce in Figure 5, convective plumes never become fully buoyant in the coordinate frame in this model, and regime 3 of fully developed buoyant convection is never reached. At $150\;\mathrm{ms}$ after the bounce, convection has all but disappeared and the buoyant mass has plummeted. At this point, SASI has taken over from neutrino-driven convection as the dominant hydrodynamical instability (see Section 4.2). It is the driving agent for the large anisotropic motions visible at late times in Figure 8.

Finally, in model $s27\mathrm{MR}{f}_{\mathrm{heat}}0.8$ with weak neutrino heating, convective instability is weak and only intermittent. As in the other models, the amount of buoyant mass peaks at $\sim 50\;\mathrm{ms}$ after the bounce, but convection weakens quickly and is almost gone at $80\;\mathrm{ms}$ after the bounce, as is obvious from the entropy snapshot of this model shown in Figure 5. SASI dominates the postbounce hydrodynamics in this model and is responsible for the strong anisotropic dynamics diagnosed via ${\upsilon }_{\mathrm{aniso}}$ in Figure 8 at later postbounce times.

Turbulence has recently become the center of attention in core-collapse supernova theory and simulation (Murphy & Meakin 2011; Murphy et al. 2013; Couch & Ott 2015). In the absence of very rapid core rotation and strong magnetic fields (the most likely scenario for the vast majority of massive stars; Heger et al. 2005; Ott et al. 2006), there is no physical source of viscosity in the postshock gain layer that could prevent neutrino-driven convection from developing into high Reynolds number turbulence (see Appendix A for a more detailed discussion of physical viscosity in the gain layer). Similarly, shear flows and entropy gradients due to periodic shock shape variations driven by SASI will seed turbulence behind the stalled shock.

A growing number of core-collapse supernova studies analyzing turbulence are showing that one of the key differences between 2D and 3D simulations is the well known (e.g., Kraichnan 1967) inverse and unphysical 2D turbulent cascade that drives kinetic energy toward large scales in 2D instead of toward small scales in three dimensions (e.g., Hanke et al. 2012; Couch 2013; Dolence et al. 2013; Takiwaki et al. 2014; Couch & O'Connor 2014; Couch & Ott 2015). Simulations suggest that kinetic energy at large scales is favorable for explosion, which may explain why 2D simulations appear to explode more easily than 3D simulations in many studies. Moreover, work by Murphy et al. (2013) and Couch & Ott (2015) demonstrated that the effective pressure generated by turbulent stress in the postshock region is an important contribution to the overall pressure behind the shock and likely pivotal in launching an explosion against the preshock ram pressure of accretion.

Turbulence in the postshock region of core-collapse supernovae is anisotropic in the radial direction and quasi-isotropic in nonradial motions (Murphy & Meakin 2011; Murphy et al. 2013; Handy et al. 2014; Couch & Ott 2015). It is mildly compressible (reaching pre-explosion Mach numbers of ∼0.3–0.5; Couch & Ott 2013) and only quasi-stationary. In the following, we focus on the kinetic energy spectra of turbulence in our simulations and compare neutrino-driven convection-dominated and SASI-dominated regimes of postbounce hydrodynamics.

We study the spectrum of turbulent motion in our simulations by decomposing the kinetic energy density of the nonradial motion into spherical harmonics on a spherical shell in the gain layer. Following previous work by Hanke et al. (2012), Couch (2013), Dolence et al. (2013), Couch & O'Connor (2014), and Handy et al. (2014), we define coefficients

$$\begin{eqnarray}&&{\epsilon }_{{\ell }m}(t)=\oint \sqrt{\rho (\theta ,\phi )}\;{\upsilon }_{t}{Y}_{{\ell }}^{m}(\theta ,\phi )d{\rm{\Omega }},\end{eqnarray} \tag{ 9 }$$

where ${\upsilon }_{t}=\sqrt{{\upsilon }_{\theta }^{2}+{\upsilon }_{\phi }^{2}}$ and where we average the $\sqrt{\rho }{\upsilon }_{t}$ part within the radial shell $r\in ({R}_{1},{R}_{2})$. In our analysis, we use ${R}_{1}=0.7{R}_{\mathrm{shock},\ \mathrm{min}}$, ${R}_{2}=0.8{R}_{\mathrm{shock},\mathrm{min}}$, where ${R}_{\mathrm{shock},\mathrm{min}}$ is the minimum shock radius at the time we carry out the spatial averaging (we also tested variations of R₁ and R₂, i.e., (0.7–0.9)${R}_{\mathrm{shock},\mathrm{min}}$ and $(0.6-0.8){R}_{\mathrm{shock},\mathrm{min}}$ and found no significant difference in the spectra). The total angular kinetic energy density at a given ℓ is then

$$\begin{eqnarray}&&E({\ell })=\displaystyle \sum _{m=-{\ell }}^{{\ell }}{\epsilon }_{{\ell }m}^{2}.\end{eqnarray} \tag{ 10 }$$

In order to calculate $E({\ell })$ at time t, we additionally average $E({\ell })$ over the time interval $(t-{\rm{\Delta }}t,t+{\rm{\Delta }}t)$, where we take ${\rm{\Delta }}t=5\;\mathrm{ms}$ in our analysis. We note that in the literature it is more common to express the turbulent energy spectrum in terms of the wave number k instead of ℓ. However, since we are decomposing the nonradial motion on a spherical shell, spherical harmonics are the natural choice. We expect $E({\ell })$ to be a power law $\propto {{\ell }}^{-\alpha }$, with α varying between different ranges in ℓ. Any power-law spectrum $E(k)\propto {k}^{-\alpha }$ corresponds to $E({\ell })\propto {{\ell }}^{-\alpha }$ in the limit of large ℓ (e.g., Chapter 21 of Peebles 1993) and, as pointed out by Hanke et al. (2012), the power-law indices of $E({\ell })$ and E(k) should correspond well to each other already at ${\ell }\gtrsim 4$.

Studies of 3D turbulent flows in various scenarios have shown that the spectrum of turbulent motion $E({\ell })$ consists of three different regions (e.g., Pope 2000). The energy of the turbulent flow is supplied in the energy-containing range at large spatial scales comparable to the size of the turbulent region by creating large-scale turbulent eddies with ℓ ~ a few. In the energy-containing range, $E({\ell })$ is typically nearly constant or increases mildly with ℓ. The inertial range is the range in ℓ in which energy cascades (i.e., is transferred) from large-scale eddies down to small scales and $E({\ell })$ decreases with ${{\ell }}^{-\alpha },\alpha \gt 1$. In the dissipation range, the dependence of $E({\ell })$ on ℓ is significantly steeper than in the inertial range, typically $E({\ell })\propto \mathrm{exp}(-{\ell })$ (e.g., Pope 2000). Our simulations do not contain any physical viscosity (which would, in any case, be extremely small in the postshock gain layer; see Appendix A) and dissipation is due to the numerical viscosity inherent to our hydrodynamics scheme.

In the Kolmogorov theory of isotropic, incompressible, stationary turbulence (e.g., Landau & Lifshitz 1959), $E({\ell })\propto {{\ell }}^{-5/3}$ in the inertial range. For the case of neutrino-driven convection in the gain layer, we expect a similar or even steeper scaling, since (1) turbulence is more or less isotropic in the nonradial directions considered here (Murphy et al. 2013), (2) turbulence has sufficient time to fully develop, since the pre-explosion stalled-shock phase lasts for many turnover cycles, and (3) a higher Mach-number (more compressible) flow generally leads to a more efficient turbulent cascade to small scales, and thus a steeper power law (e.g., Garnier et al. 2000).

The top panel of Figure 9 shows $E({\ell })$ at various postbounce times in model $s27\mathrm{MR}{f}_{\mathrm{heat}}1.05$, whose gain-layer hydrodynamics is dominated by neutrino-driven convection due to strong neutrino heating (see Section 4.3). While there are variations in $E({\ell })$ in the low-ℓenergy-containing range, at ${\ell }\gtrsim 10$ the spectra are quite steady after $t-{t}_{b}\sim 80\;\mathrm{ms}$, indicating that the flow is at least quasi-stationary at intermediate and small scales in this model. $E({\ell })$ should peak at ℓ, corresponding to the size of the convectively unstable gain region. At $90\;\mathrm{ms}$ after the bounce we infer from the top right panel of Figure 8(a) the radial extent of the turbulent region of $H\sim 70\;\mathrm{km}$ and a typical radius of $R\sim 90\;\mathrm{km}$ (the center of the convective region). The value of ℓ at which the spectrum $E({\ell })$ peaks should correspond to the number of eddies with diameter H that fit into the turbulent region, ${{\ell }}_{\mathrm{peak}}\sim (2\pi R)/H-1\approx 7$. This is close to what is realized by the spectrum at $90\;\mathrm{ms}$ after the bounce shown in Figure 9 for this model. At smaller scales (larger ℓ), the spectrum should first exhibit an extended inertial range region with $E({\ell })\propto {{\ell }}^{-5/3}$ before steepening in the dissipation range at very large ℓ. This, however, is not borne out by Figure 9. At intermediate ℓ of 10–40, the spectrum is much shallower than ${{\ell }}^{-5/3}$ and most consistent with ℓ⁻¹, and it steepens only at ${\ell }\gtrsim 40$ and quickly surpasses the ${{\ell }}^{-5/3}$ scaling. This kind of spectral behavior is qualitatively and quantitatively consistent with what was found for neutrino-driven turbulence in the simulations of Dolence et al. (2013), Couch & O'Connor (2014), and Couch & Ott (2015), who all used numerical methods and Cartesian grid setups very similar to ours.

Zoom In Zoom Out Reset image size

The bottom panel of Figure 9 shows $E({\ell })$ at various postbounce times in the SASI-dominated model $s27\mathrm{MR}{f}_{\mathrm{heat}}0.8$ (weak neutrino heating). Anisotropic motions in this model (and at postbounce times $\gtrsim 150\;\mathrm{ms}$ also in model $s27\mathrm{MR}{f}_{\mathrm{heat}}0.95$ whose $E({\ell })$ is not shown) are driven by entropy and vorticity perturbations caused by the SASI, which is much more intermittent than neutrino heating. This is reflected in the turbulent kinetic energy spectra that vary at all scales with postbounce time and do not reach the quasi-stationarity that we observe for neutrino-driven turbulent convection in the top panel of Figure 9. The variations in $E({\ell })$ in the SASI-dominated model can be directly correlated with the strength of SASI. For example, the overall magnitude of $E({\ell })$ grows from 90 to 150 ms after the bounce, which coincides with the increasing strength of SASI oscillations seen in Figure 6 for this model. At $\sim 180\;\mathrm{ms}$, $E({\ell })$ at large scales is decreased as a result of damped SASI oscillations shortly after the silicon interface advects through the shock (see Section 4.2). While there is much variation in the overall magnitude of $E({\ell })$, the scaling of $E({\ell })$ is significantly shallower than ${{\ell }}^{-5/3}$ and closer toℓ⁻¹ at the scales one would naively be tempted to identify with the inertial range. This is in agreement with the neutrino-driven turbulent convection case.

Several authors (e.g., Dolence et al. 2013; Couch & O'Connor 2014) have argued that the ℓ⁻¹ scaling observed in contemporary 3D simulations could be due to the physical nature of the postshock turbulent flow that deviates significantly from the assumptions of Kolmogorov turbulence. Our interpretation is different. An inertial range scaling with ${{\ell }}^{-\alpha }$ with $\alpha \leqslant 1$ is unphysical, since in the limit of infinite resolution, the integral turbulent energy is divergent. Neutrino-driven turbulence is essentially isotropic in the nonradial directions, is quasi-stationary, and is only mildly compressible. Local high-resolution studies of driven turbulence in this regime generally find an inertial range with $\alpha \simeq 5/3$ for the incompressible transverse flow component and $\alpha \gt 5/3$ for the compressible part (e.g., Schmidt et al. 2006). Those simulations and simulations of turbulence in other regimes (e.g., Porter et al. 1998; Sytine et al. 2000; Dobler et al. 2003; Kaneda et al. 2003; Haugen & Brandenburg 2004; Kritsuk et al. 2007; Federrath 2013), however, also find the appearance of a shallower region with $\alpha \sim 1$ near the end of the inertial range before the transition to the dissipation range. This corresponds to inefficient energy transport at these scales and is referred to as the bottleneck effect. This is understood to be a physical feature of turbulence that is related to a partial suppression of nonlinear interactions of turbulent eddies of different scale near the regime of strongest dissipation (She & Jackson 1993; Yakhot & Zakharov 1993; Falkovich 1994; Verma & Donzis 2007; Frisch et al. 2008).

Sytine et al. (2000) carried out a resolution study with local compressible (Mach 0.5) freely decaying turbulence simulations using the original PPM solver of Colella & Woodward (1984). Their Figure 11 shows that their local simulations with 1024³ and 512³ cells resolve an inertial range with $\alpha =5/3$. The bottleneck with $\alpha \lt 5/3$ appears at the end of this range. However, with decreasing resolution, the bottleneck shifts to progressively lower wavenumbers, consuming more and more of the resolved inertial range. Already at 256³, the inertial range is gone and energy injection and dissipation scales are joined directly with $1\lesssim \alpha \lt 5/3$. On the basis of their results and previous work by Porter et al. (1998), Sytine et al. (2000) argue that the numerical viscosity of their PPM scheme provides dissipation that directly affects the flow on spatial scales from 2 to ∼12 times the width of a computational cell. This should be the equivalent of the dissipation range. On somewhat larger scales, from ∼12 to ∼32 cell widths, the flow is still indirectly affected by the viscosity of PPM, creating the observed bottleneck effect. We point out that these studies focused on a specific regime of turbulence, freely decaying and isotropic, which is different from the one we observe in our simulations. In our simulations, turbulence is driven by buoyancy and is anisotropic. However, very recently, Radice et al. (2015) came to conclusions very similar to Sytine et al. (2000) also for driven anisotropic turbulence.

Our numerical hydrodynamics scheme is very similar to the PPM implementation used by Sytine et al. (2000), but likely more dissipative because we do not employ the original exact Riemann solver of Colella & Woodward (1984), but the more dissipative HLLE solver (see Section 2). The numerical viscosity of our scheme is thus larger than in the scheme of Sytine et al. (2000; see the comparison between HLLE and HLLC in Radice et al. 2015), and 32 cell widths is only a lower bound on the scale, which is affected by numerical viscosity in our simulations. In our fiducial medium resolution simulations for which we present $E({\ell })$ in Figure 9, the cell width is $\sim 1.4\;\mathrm{km}$ and the region that is turbulent has a radial extent of $\sim 70\;\mathrm{km}$ (see Figure 8). Hence, we have (in the best case) 70 km/1.4 km ≈ 50 cells covering the turbulent region (of which ∼32 are affected by numerical viscosity), which is much less than the 512 linear cell width needed by Sytine et al. (2000) to resolve an inertial range. We conclude that at the resolution employed here, we cannot reasonably expect to resolve the inertial range in the turbulent gain layer. All that we are seeing here, and that the simulations of Dolence et al. (2013), Couch & O'Connor (2014), and Couch & Ott (2015) show, is the contamination of the turbulent energy spectrum by numerical viscous effects all the way up to the energy-containing range. Turbulence is thus not resolved in these and in the present 3D simulations. This conclusion is further supported by the low numerical Reynolds number of $\mathrm{Re}\sim 70$ that we find in Appendix B for our simulations, suggesting that our simulations are somewhere between the perturbed laminar flow and turbulence. Couch & Ott (2015) estimated $\mathrm{Re}\sim 350$ via a simple comparison of the size of the convective region with linear grid spacing (e.g., Pope 2000). Using their approach, we find $\mathrm{Re}\sim 180$. Authors carrying out simulations on spherical grids have argued that they see α closer to $5/3$ and resolve the inertial range (Hanke et al. 2012; Handy et al. 2014). However, the angular and radial resolutions employed in these studies are significantly lower than the effective resolutions provided by our 3D Cartesian grids, and it is not clear how the turbulence could be resolved in their simulations if not in ours.

We explore the impact of numerical resolution in the regime of strong neutrino heating and neutrino-driven convection-dominated 3D hydrodynamics by running simulations of the s27 progenitor with a total of five different resolutions with the linear cell width in the postshock gain layer varying by almost a factor of four. Our baseline $s27\mathrm{MR}{f}_{\mathrm{heat}}1.05$ model has a resolution on the AMR level containing the postshock region and the shock with linear resolution ${{dx}}_{\mathrm{shock}}=1.416\;\mathrm{km}$. This correponds to an effective angular resolution at a radius of $100\;\mathrm{km}$ of $d(\theta ,\phi )=0\buildrel{\circ}\over{.} 81$. In models $s27\mathrm{ULR}{f}_{\mathrm{heat}}1.05$ ("ultra-low resolution"), $s27\mathrm{LR}{f}_{\mathrm{heat}}1.05$ ("low resolution"), $s27\mathrm{IR}{f}_{\mathrm{heat}}1.05$ ("intermediate resolution"), and $s27\mathrm{HR}{f}_{\mathrm{heat}}1.05$ ("high resolution"), this is $3.784$, $1.892$, $1.240$, and $1.064\;\mathrm{km}$, respectively. These correspond to effective angular resolutions at a radius of $100\;\mathrm{km}$ of $2\buildrel{\circ}\over{.} 15$, $1\buildrel{\circ}\over{.} 08$, $0\buildrel{\circ}\over{.} 81$, $0\buildrel{\circ}\over{.} 71$, and $0\buildrel{\circ}\over{.} 61$, for ULR, LR, IR, and HR, respectively (see also Table 1).

Figures 10 and 11 give a concise summary of the effects of resolution on the postbounce hydrodynamics and on the development of a neutrino-driven explosion. The overall trend is very clear: the lower the resolution, the larger the average shock radius, the higher the neutrino heating rate, the greater the heating efficiency, and the larger the mass in the gain layer.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While these overall trends are robust, there are some notable inconsistencies in the details. The mean specific entropy in the gain layer (bottom panel of Figure 11) appears almost completely independent of resolution. The asphericity of the shock front, measured by the normalized rms deviation ${\sigma }_{\mathrm{shock}}$ of the shock radius in the bottom panel of Figure 10 has no systematic resolution dependence in its magnitude and variations. The fiducial MR simulation is an outlier with the overall smallest ${\sigma }_{\mathrm{shock}}$. The shock radius, neutrino heating, heating efficiency, and mass in the gain region are very similar in the LR and MR models (differing in dx_shock by ∼30%), and at the end of its evolution the MR simulation actually has a slightly larger shock radius than its LR counterpart. On the other hand, the IR and HR simulations, which differ only by $\sim 15\%$ in resolution, are consistent with each other in all quantities except ${\sigma }_{\mathrm{shock}}$. The MR/IR simulation pair differs in resolution by $\sim 15\%$, yet their results are much farther apart than those of the LR/MR pair that differs by $\sim 30\%$ in resolution. These variations from the general trend are indicative of the possibility that many if not most (or all) of our simulations are not yet in the convergent regime. Perhaps much higher resolutions in the convectively unstable layer may be needed to accurately and in a converged manner capture the hydrodynamics of core-collapse supernovae dominated by neutrino-driven turbulent convection.

Hanke et al. (2012), Couch & O'Connor (2014), and Takiwaki et al. (2014), who carried out less extensive 3D parameter studies with similar or lower resolutions, found the same trends with resolution observed in our simulations. Handy et al. (2014), on the other hand, found improved conditions for explosion with increasing resolution. However, they studied angular grid spacings from 24° down to only 2°. Their highest resolution roughly corresponds to our ULR case. At such coarse resolutions, which suppress nonradial convective motions, it is not at all surprising that the conditions become more favorable for explosion as increasing resolution begins to allow nonradial motions. The Handy et al. (2014) simulations thus probe the resolution dependence of 3D postbounce hydrodynamics in a completely different regime than our simulations.

Figure 12 provides further evidence for why lower-resolution simulations are (artificially) favorable for neutrino-driven explosions. The lower the resolution, the larger the amount of buoyant mass (defined as the mass in the gain region with positive radial velocity) and the greater the amount of positive momentum in the gain region. The more mass is truly buoyant (and thus in regime 3 of neutrino-driven convection discussed in Section 4.3), the greater the neutrino heating rate and efficiency (see Figure 11). Note, however, that by comparing the total mass in the gain region given in the bottom panel of Figure 11 with the buoyant mass in the top panel of Figure 12, one finds that that the truly buoyant mass is at most $\sim 20\%$ of the mass in the gain region. We expect that this fraction will sensitively depend on progenitor structure and will be higher in progenitors with lower postbounce accretion rates than in the 27 ${M}_{\odot }$ progenitor that we study here.

Zoom In Zoom Out Reset image size

The bottom panel of Figure 12 shows the time evolution of the sum of the angle-averaged "turbulent" radial fluxes of enthalpy (F_C, also known as "convective flux") and kinetic energy (F_K) near the shock. We follow Hurlburt et al. (1986) and Handy et al. (2014) and define

$$\begin{eqnarray}{F}_{{\rm{C}}} & = & {\displaystyle \int }_{4\pi }\rho {\upsilon }_{{\rm{r}}}{\left(\epsilon +\displaystyle \frac{P}{\rho }\right)}^{\prime }{r}^{2}d{\rm{\Omega }},\\ {F}_{{\rm{K}}} & = & {\displaystyle \int }_{4\pi }\rho {\upsilon }_{{\rm{r}}}{\left(\displaystyle \frac{1}{2}{\upsilon }_{i}{\upsilon }_{i}\right)}^{\prime }{r}^{2}d{\rm{\Omega }},\end{eqnarray} \tag{ 11 }$$

where ρ is the density, ${\upsilon }_{{\rm{r}}}$ is the radial velocity, is the internal energy, P is the pressure, and ${\upsilon }_{i}$ is the ith component of velocity. All primed quantities represent variations about the angle-averaged mean so that, for instance, F_K measures the amount of the turbulent part of the specific kinetic energy (note that by construction $\bar{{\bar{\upsilon }}_{i}{\upsilon }_{i}^{\prime }}\equiv 0$, where $\bar{\cdot }$ denotes the angular average). We evaluate the angular integrals in Equation (11) at each time at a radius that corresponds to the instantaneous minimum shock radius. A number of studies (e.g., Burrows et al. 1995; Couch 2013; Dolence et al. 2013; Ott et al. 2013; Handy et al. 2014; Couch & Ott 2015) have argued that buoyant convective/turbulent bubbles are locally important in driving shock deformation and expansion. Murphy et al. (2013) and Couch & Ott (2015) have shown that the additional effective ram pressure due to turbulence is crucial for the relative ease of explosions in two dimensions and three dimensions compared with the 1D case. All these effects are related to the convective/turbulent flux of kinetic energy and enthalphy in the gain layer and near the shock (see Yamasaki & Yamada 2006). Figure 12 reveals that the sum ${F}_{{\rm{C}}}+{F}_{{\rm{K}}}$ near the shock decreases with increasing resolution, creating less favorable conditions for explosion.

The radial convective/turbulent fluxes are dominated by the flow at large and intermediate scales. In Figure 13, we plot angular turbulent kinetic energy spectra $E({\ell })$ (see Equation (10)) in the gain layer at $90\;\mathrm{ms}$ after the bounce for all resolutions. The left panel shows the plain $E({\ell })$ spectra, while the right panel shows "compensated" spectra that are rescaled by ${{\ell }}^{-5/3}$ as is customary in studies of Kolmogorov turbulence. A flat graph in the region where the inertial range is expected would indicate consistency with Kolmogorov turbulence. Given the spatial scale of the gain layer in our simulations, we would expect the energy-containing range to be around ${\ell }\sim 7$ (see Section 4.4), which should be followed by an inertial range with $E({\ell })\propto {{\ell }}^{-5/3}$ before dissipation sets in. None of our simulations, not even the HR case, exhibits any inertial range. Where the inertial range should be, $E({\ell })$ is most consistent with an ℓ⁻¹ scaling, which is indicative of a bottleneck due to viscous contamination because of insufficient numerical resolution (see Section 4.4).

Zoom In Zoom Out Reset image size

Figure 13 does not clearly show large differences of $E({\ell })$ in the energy-containing range with changing resolution. However, note that at low ℓ the spectra are not fully stationary (see Figure 9). One should also recall that we here project out the radial part and that the important radial kinetic energy and enthalpy fluxes decrease with increasing resolution, which indicates less total energy/power at large scales (Figure 12). The figure does, however, clearly demonstrate that transport of turbulent energy to small scales becomes increasingly efficient as the resolution is increased. The energy contained at large ℓ increases systematically with resolution and even appears to converge as the resolution gets close to the HR case. However, the resolution decrements between the various shown simulations are not constant and the three highest simulations differ only by $\sim 15\%$ in resolution, while MR and LR differ by $\sim 30\%$ and LR and ULR differ by a factor of two. Since no inertial range is realized, we do not consider any of our studied resolutions to be in the regime in which the flow is truly turbulent. The HR simulation, at $\sim 90\;\mathrm{ms}$ after the bounce, covers the entire postshock region with $\sim {240}^{3}$ computational cells, but only the outer $\sim 70\;\mathrm{km}$ are actually convectively unstable and are effectively covered by $70.0\;\mathrm{km}/1.064\;\mathrm{km}\approx 66$ linear cell widths. According to Sytine et al. (2000) this resolution may still be a factor of ≳7–8 too low for resolving the inertial range.

We investigate the resolution dependence in the weak neutrino heating, SASI-dominated case by comparing our baseline-resolution simulation $s27\mathrm{MR}{f}_{\mathrm{heat}}0.8$ with a simulation carried out with lower resolution, $s27\mathrm{LR}{f}_{\mathrm{heat}}0.8$, which uses the same resolution of the LR simulation as in the previous section. MR and LR resolutions differ by $\sim 30\%$ (see Table 1). Additional simulations with further decreased or increased resolution would be advisable but were not possible for the SASI-dominated case within the limitations of our computational resources.

The top panel of Figure 14 compares the evolution of the average shock radius in the MR and LR simulations. They are qualitatively and quantitatively nearly identical and significantly closer to each other than the LR and MR simulations in the strong neutrino heating case discussed in the previous Section 5.1. We also find (and show in Figure 3) that integral net neutrino heating, neutrino heating efficiency, and the mass in the gain region are very similar in the MR and LR models throughout the simulated postbounce time.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While the average shock radius evolves nearly identically in the MR and LR cases, we find that deviations from the average due to SASI oscillations are smaller in the LR case. This is apparent from the bottom panel of Figure 14, which shows the normalized rms deviation ${\sigma }_{\mathrm{shock}}$ of the shock radius from its angle-averaged value. The oscillations in ${\sigma }_{\mathrm{shock}}$ are due to SASI, and their amplitudes are much smaller in the LR simulation. Figure 15 depicts the evolution of the normalized ${\ell }=2,m=2$ and ${\ell }=1,m=-1$ amplitudes (Equation (3)) of the shock front as representative examples of the ${\ell }=\{1,2\}$ mode families in the LR and MR simulations. The evolution of these modes is qualitatively similar in both LR and MR simulations, but the LR simulation shows systematically lower mode amplitudes in both ${\ell }=1$ and ${\ell }=2$ until $\sim 160\;\mathrm{ms}$ after the bounce. At that time, the silicon interface advects through the shock, leading to its transient expansion, and to a profound change in the SASI mode structure (see Section 4.2). In the LR simulation, the ${\ell }=2$ mode amplitudes decay less than in the MR case, but the ${\ell }=1$ modes do not grow as strongly as in the MR case. This deviation between MR and LR SASI dynamics has an effect on the average shock radius, whose MR and LR evolutions depart from each other toward the end of the LR simulation at $\sim 200\;\mathrm{ms}$ after the bounce.

Our results show that the weak neutrino heating, SASI-dominated regime of 3D postbounce hydrodynamics is sensitive to resolution and this sensitivity is strongest in the development and nonlinear dynamics of SASI. Sato et al. (2009) have shown that for SASI to reach convergence, the numerical resolution must be sufficiently high to capture the full advective-acoustic cycle of entropy/vorticity perturbations that advect through the postshock region, are reflected at the protoneutron star, and propagate back up to the shock. The LR simulation (${{dx}}_{\mathrm{shock}}=1.892\;\mathrm{km}$) has evidently too low resolution, but since we only have two resolutions at hand, we cannot with confidence say that the MR simulation (${{dx}}_{\mathrm{shock}}=1.416\;\mathrm{km}$) is in the convergent regime for SASI.