PARAMETRIC STUDY OF FLOW PATTERNS BEHIND THE STANDING ACCRETION SHOCK WAVE FOR CORE-COLLAPSE SUPERNOVAE

The numerical method is the same as that used in our previous paper (Iwakami et al. 2008). The ZEUS-MP/2 code (Hayes et al. 2006) is modified to solve the accretion flow irradiated by neutrinos emitted from the proto-neutron star (PNS); the tabulated EOS, by Shen et al. (1998), is implemented in the code; the light bulb approximation (Ohnishi et al. 2006, for details) is adopted for the neutrino heating and cooling. It is assumed that neutrinos are emitted isotropically from PNS with prescribed fluxes, and the matter outside PNS is optically thin. Only electron-type neutrinos and anti-neutrinos are considered in this study. Their temperatures are also assumed to be constant in time and set to the typical values in the post-bounce phase, $T_{\nu _{\mathrm{e}}} = 4$ MeV and $T_{\bar{\nu }_{\mathrm{e}}} = 5$ MeV. The mass of PNS is fixed to M_PNS = 1.4 M_☉. As model parameters, the neutrino luminosity and mass accretion rate are varied in the range of L_ν = 2.0–6.0 × 10⁵² erg s⁻¹ and $\dot{M} = 0.2\hbox{--}1.0\, M_{\odot }$ s⁻¹, respectively. All models presented in this paper are summarized in Table 1. More than three different realizations are computed for each model.

Table 1. Summary of Parameters for All Models

Model	$\dot{M}$a	Lνb	rinc	routd	texpe
	(M☉ s−1)	(1052 erg s−1)	(km)	(km)	(ms)
A	0.2	2.0	29	581	⋅⋅⋅
B	0.2	2.5	33	655	⋅⋅⋅
C	0.2	3.0	36	712	270–352
D	0.6	4.0	41	822	⋅⋅⋅
E	0.6	4.5	44	872	⋅⋅⋅
F	0.6	5.0	46	920	201–208
G	1.0	5.0	46	919	⋅⋅⋅
H	1.0	5.5	49	965	⋅⋅⋅
I	1.0	6.0	51	1007	358–411

Notes. ^aThe mass accretion rate. ^bThe neutrino luminosity. ^cThe radial position of the inner boundary. ^dThe radial position of the outer boundary. ^eThe time for the shock wave to reach r_out.

Download table as:  ASCIITypeset image

A staggered grid is used in the ZEUS-MP/2. The spherical polar coordinates (r, θ, ϕ) are adopted. The computational domain covers the entire solid angle. The radial inner boundary is located at r_in ∼ 28–50 km, and the outer boundary is located at r_out ∼ 500–1000 km. The computational region is divided into 300 × 30 × 60 grid cells in r × θ × ϕ directions, as explained below. The radial size, Δr_i, and position, r_i, of the i-th cell are given as

$$\begin{equation} \Delta r_{i} = r_c \Delta r_{i-1},\ \ \ r_i = \Delta r_i/(r_c -1),\ \ \ (i=1, 2, \cdot \cdot \cdot, N), \end{equation} \tag{ 1 }$$

where N denotes the total number of grid cells, and r_c is the common ratio of the geometric sequence fixed to be 1.01; that is, the radial resolution is 1% of the radius. The inner boundary of the computational domain is located at ρ = 10¹¹ g cm⁻³ in the spherically symmetric, steady, unperturbed states. It roughly coincides with the neutrino sphere and hence depends on L_ν and $\dot{M}$. The locations of outer boundary r_out are also different among models, since N is fixed to 300. The values of r_in and r_out in each model are listed in Table 1.

The boundary conditions for other quantities are set as follows. The velocity, density, internal energy, and electron fraction are fixed to their initial values at the outer boundary. At the inner boundary, on the other hand, the density, internal energy, and electron fraction are assumed to have a vanishing gradient in the direction normal to the inner boundary, and the radial velocity is fixed to the initial value, whereas θ and ϕ components are set so that the conservation of angular momentum should be conserved.

The initial conditions, or spherically symmetric, steady, shocked accretion flows, are obtained in the same manner as Yamasaki & Yamada (2006), except that we do not take self-gravity into account and determine the location of the inner boundary based on the density, not on the optical depth, for simplicity. The profiles of radial velocity and entropy in the initial unperturbed states are displayed in Figure 1. As observed in the upper panels, the standing shock wave abruptly decelerates the inward flow and the shocked matter slows down further as it approaches the PNS. As shown in the lower panels, on the other hand, all of the models considered in this paper have a region with a negative entropy gradient, which corresponds to the gain region, as is well known. The combinations of L_ν and $\dot{M}$ adopted in this paper hence give the models that are unstable to convection in the classical sense. Discussions on this issue based on the so-called χ parameter (Foglizzo et al. 2006) will be given later. It will also be useful to mention that the parameter combinations correspond to those that give the oscillatory, unstable modes in the linear analysis by Yamasaki & Yamada (2007), which the authors interpreted as SASI. Furthermore, we confirm that our models are stable to radial perturbations, also consistent with the results of Yamasaki & Yamada (2007). In order to induce non-spherical instabilities, we add cell-wise, random perturbations within 1% to the radial velocity at the beginning of computation.

Zoom In Zoom Out Reset image size

We employ mode analysis to judge which flow pattern is dominant. The deformation of shock surface can be expanded as a linear combination of the spherical harmonic components,

$$\begin{equation} R_\mathrm{sh}(\theta, \phi, t) = \sum ^{\infty }_{l=0} \sum ^{l}_{m=-l} c^m_{l} (t) \, Y^m_{l}(\theta, \phi), \end{equation} \tag{ 2 }$$

where $Y^m_{l}$ is expressed by the associated Legendre polynomial, $P^m_{l}$, as

$$\begin{equation} Y^m_{l} = K^m_{l} P^m_{l}(\cos \theta) \, e^{im\phi }, \ \ \ K^m_{l} = \sqrt{\frac{2l+1}{4\pi }\frac{(l-m)!}{(l+m)!}}. \end{equation} \tag{ 3 }$$

The expansion coefficients can be obtained as follows,

$$\begin{equation} c^m_{l} (t) =\int ^{2\pi }_0 \! \! \! \! d\phi \! \int ^{\pi }_0 \! \! d\theta \, \sin \theta \, R_\mathrm{sh}(\theta, \phi, t) \, Y^{m*}_{l} (\theta, \phi), \end{equation} \tag{ 4 }$$

where the superscript * denotes complex conjugation.

We define the following quantities:

$$\begin{equation} A_l(t) = \sqrt{\Sigma ^l_{m=-l} |c^m_l(t)/c^0_0(t)|^2}, \end{equation} \tag{ 5 }$$

$$\begin{equation} A_{1, 2}(t) = \sqrt{\Sigma ^2_{l=1}\Sigma ^l_{m=-l} |c^m_l(t)/c^0_0(t)|^2}, \end{equation} \tag{ 6 }$$

$$\begin{equation} A_{4, 5} (t) = \sqrt{\Sigma ^5_{l=4}\Sigma ^l_{m=-l} |c^m_l(t)/c^0_0(t)|^2}, \end{equation} \tag{ 7 }$$

and also utilize their time averages:

$$\begin{equation} \bar{A}_l = \frac{1}{T}\int ^{t_e}_{t_s} A_l(t) dt, \end{equation} \tag{ 8 }$$

$$\begin{equation} \bar{A}_{1,2} = \frac{1}{T}\int ^{t_e}_{t_s} A_{1,2}(t) dt, \end{equation} \tag{ 9 }$$

$$\begin{equation} \bar{A}_{4,5} = \frac{1}{T}\int ^{t_e}_{t_s} A_{4,5}(t) dt, \end{equation} \tag{ 10 }$$

where T = t_e − t_s is the integral time, t_s is the starting time, and t_e is the ending time for integration.

The classical condition of convective instability is the existence of the negative entropy gradient. In the presence of matter flow, however, this condition is not sufficient. In fact, matter passes through a convectively unstable region in a finite time, in general. If this advection time is shorter than the growth time of convection, then the instability is suppressed. Foglizzo et al. (2006) propose a new criterion that takes the advection of matter into account. They introduce the χ parameter, which is the ratio of the advection time to the local timescale of buoyancy and is defined as

$$\begin{equation} \chi \equiv \int _{r_{\rm gain}}^{r_{\rm sh}} \left|\frac{N}{u_r}\right| dr, \end{equation} \tag{ 11 }$$

where the integration is done over the gain region with r_gain being the gain radius, i.e., the boundary between the cooling and heating region, and r_sh being the shock radius. u_r is the radial velocity. The Brunt-Väisälä frequency N is defined as

$$\begin{equation} N^2\equiv \left[ \frac{1}{p}\left(\frac{\partial p}{\partial S}\right)_{\rho,Y_e}\frac{dS }{dr} +\frac{1}{p}\left(\frac{\partial p}{\partial Y_e}\right)_{\rho, S}\frac{dY_e}{dr} \right] \frac{g}{\Gamma _1}, \end{equation} \tag{ 12 }$$

$$\begin{equation} \Gamma _1 = \left(\frac{\partial \ln {p}}{\partial \ln {\rho }}\right)_{S, Y_e}, \end{equation} \tag{ 13 }$$

where p, ρ, S, Y_e, and G are the pressure, density, entropy, electron fraction, and gravitational constant, respectively. g is the gravitational acceleration and is approximately given in the gain region as g = (GM_PNS/r²), in which M_PNS is the PNS mass. Substituting Equation (13) into Equation (12), we obtain the following expression,

$$\begin{equation} N^2 =\left|\frac{1}{\Gamma _1 p}\frac{dp}{dr}-\frac{1}{\rho }\frac{d\rho }{dr}\right|g, \end{equation} \tag{ 14 }$$

which is actually used for the calculation of χ in the following.

The above formula, obtained originally in the linear analysis (Foglizzo et al. 2006), is meant for application to 1D steady states. In the literature, however, it has been also applied to somehow defined mean flows in the full-fledged turbulent induced by convection (e.g., Fernández et al. 2013; Couch & O'Connor 2013). Following these practices, for non-exploding models in this paper, we define the mean flow by time-averaging the quantities in the flow from the onset of the quasi-steady phase up to the end of computation:

$$\begin{equation} \bar{q}(r,\theta,\phi)=\frac{1}{T}\int ^{t_e}_{t_s} q(r,\theta,\phi,t) dt, \end{equation} \tag{ 15 }$$

where q(r, θ, ϕ, t) means an arbitrary quantity (i.e., ρ, u_r, and so on). In applying the above criterion, we further take an angle-average over the entire solid angle at each radius,

$$\begin{equation} \bar{q}_{\rm 1D}(r) = \frac{1}{4\pi }\int _0^{4\pi } \bar{q}(r,\theta,\phi) d\Omega. \end{equation} \tag{ 16 }$$

The parameter $\bar{\chi }_{\rm 1D}$ is then obtained as

$$\begin{equation} \bar{\chi }_{\rm 1D} = \int _{\bar{r}_{\rm gain 1D}}^{\bar{r}_{\rm sh 1D}} \left|\frac{\bar{N}_{\rm 1D}(r)}{\bar{u}_{r {\rm 1D}}(r)}\right| dr, \end{equation} \tag{ 17 }$$

$$\begin{equation} \bar{N}^2_{\rm 1D} =\left|\frac{1}{\bar{\Gamma }_{1\rm 1D} \bar{p}_{\rm 1D}}\frac{d\bar{p}_{\rm 1D}}{dr}-\frac{1}{\bar{\rho }_{\rm 1D}}\frac{d\bar{\rho }_{\rm 1D}}{dr}\right|g, \end{equation} \tag{ 18 }$$

where the quantities with a bar are calculated for the angle-averaged mean flow. The thermodynamical variables, $\bar{p}_{\rm 1D}$ and $\bar{\Gamma }_{1\rm 1D}$, are given by the EOS table as a function of $\bar{\rho }_{\rm 1D}$, $\bar{e}_{\rm 1D}$, and $\bar{Y}_{e\rm 1D}$.

In this section, we describe the basic characteristics of each model presented in Table 2.

Table 2. Flow Patterns and Some Parameters in the Semi-Nonlinear and Nonlinear Phases

Model	Pattern	χa	rgain b	rshc	Pattern	$\bar{\chi }$ d	$\bar{r}_{\rm gain}$ e	$\bar{r}_{\rm sh}$ f	$\bar{A}_{1,2}/\bar{A}_{4,5}$ g
	(Semi-Nonlinear)		(km)	(km)	(nonlinear)		(km)	(km)
A	SL/SP	1.2	52	68	SL/SP	2.4/2.2	49/49	107/104	9.6/7.3
B	SP	2.7	56	89	SPB	3.8	55	127	5.1
C	BB	4.9	63	112	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅
D	SP	1.6	71	92	SP	2.1	69	139	6.7
E	SP	2.5	75	109	BB	3.6	74	142	3.0
F	BB	4.1	82	132	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅
G	SP	1.4	78	97	SP	1.6	77	143	7.4
H	SP	2.0	82	110	SPP	2.5	81	182	6.0
I	SP	2.9	87	128	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅

Notes. The flow patterns are classified into sloshing motion (SL), spiral motion (SP), multiple buoyant bubbles (BB), spiral motion with buoyant-bubble formation (SPB), and spiral motion with pulsating change of rotational velocities (SPP). For the explosion models, only flow patterns in the semi-nonlinear phase are written in this table. All values listed above are averaged among different realizations. ^aχ parameter in the initial flow. ^bThe gain radius in the initial flow. ^cThe shock radius in the initial flow. ^dχ parameter in the time- and angle-averaged flow during the nonlinear phase. ^eThe gain radius in the time- and angle-averaged flow during the nonlinear phase. ^fThe shock radius in the time- and angle-averaged flow during the nonlinear phase. ^gThe ratio of the time-averaged mode amplitudes of A_{1, 2} to A_{4, 5}.

Download table as:  ASCIITypeset image

The time evolutions of the shock radius are shown in Figure 2, obtained from the expansion coefficient $c_0^0$ in Equation (4). The solid and dashed lines correspond to the non-explosion models, and the dotted lines correspond to the explosion models. In the non-explosion models, the average shock radius is larger for higher neutrino luminosities. In the explosion models, we observe the monotonic increase of the shock radius up to the outer boundary.

Zoom In Zoom Out Reset image size

In both cases, the shock radius is initially almost constant for a certain period. In this phase, the perturbation grows exponentially from the initial small value with no substantial feedback to the background flow. In this paper, we refer to this period as the linear phase. In the following period, which we call the semi-nonlinear phase, there is a discernible increase of the average shock radius that continues to the first peak in the non-explosion models. As shown later in Figure 5, the ending time of the semi-nonlinear phase corresponds to the first turning point from exponential growth to its saturation or relaxation in the time evolution of mode amplitudes. Nonlinear mode couplings cannot be ignored any more in this phase, and a dominance of a particular mode becomes evident. The semi-nonlinear phase is followed by the nonlinear phase, in which the hydrodynamical instabilities (convection and/or SASI) are saturated, and quasi-steady states are established with some intermittent activities in the non-explosion models. In the following, we focus on the flow patterns realized in the semi-nonlinear and nonlinear phases.

Here, we introduce three well-known post-shock flow patterns: sloshing motion, spiral motion, and formation of multiple buoyant bubbles with high entropies. Some snapshots of entropy distribution in a cross-section are presented as contours in Figure 3. The inner- and outermost contour lines correspond to the PNS surface and the shock front, respectively. The sloshing and spiral motions are accompanied by global deformations of shock front, whereas the formation of buoyant bubbles generates shock deformations on a smaller scale. The sloshing motion oscillates the shock front along a certain axis. Figure 3(a) shows an example of the sloshing motion, which occurs along the x-axis in this case. The spiral motion is intrinsically non-axisymmetric with rotations of matter around a certain axis, although the accretion is spherically symmetric outside of the shock wave. Figure 3(b) presents such an example with the rotation axis almost perpendicular to the x − z-plane. As shown in Figure 3(c), on the other hand, the formation of high-entropy bubbles in the gain region tends to deform the shock wave rather locally.

Zoom In Zoom Out Reset image size

We find from our simulations that the post-shock flows in the semi-nonlinear and nonlinear phases are classified into one of the above three patterns: sloshing motion (abbreviated as SL in the following), spiral motion (SP), and buoyant-bubble formation (BB), or into the following intermediate patterns: spiral motion with bubble formation (SPB) and spiral motion with pulsating changes of rotational velocity (SPP). Interestingly, the flow pattern is not always identical between the semi-nonlinear and nonlinear phases. Hence, we summarize the results for the two phases separately in Table 2. In order to facilitate the grasp of the trend in the pattern realization, we also show the same results in the $L_\nu -\dot{M}$ diagram in Figure 4. In this figure, the flow patterns given on the left and right sides of arrows correspond to the flow patterns realized in the semi-nonlinear and nonlinear phases, respectively. As discussed later, the pattern realized in the post-shock flow is rather robust; there are some exceptions, in which pattern change is observed when the initial perturbation is changed. In fact, two patterns on the same site in Figure 4 implies that one of them is realized depending on the initial perturbation. Note that for the explosion models the flow patterns are determined only in the semi-nonlinear phase, since the quasi-steady states are not established after the semi-nonlinear phase.

Zoom In Zoom Out Reset image size

In Figure 4, the dotted curve roughly indicates the critical luminosity estimated from the simulations in this study. Although it agrees well with the results obtained by Yamasaki & Yamada (2006) for almost the same settings, it is generically higher than those found in more realistic simulations (Nordhaus et al. 2010; Hanke et al. 2012; Couch 2013). Considering the simplifications adopted in these papers, we think that this difference is not surprising.

We find that SL and SP tend to appear in the higher $\dot{M}$ and lower L_ν, whereas BB prevails in the lower $\dot{M}$ and higher L_ν. This is consistent with Burrows (2012) and Hanke et al. (2013). Burrows (2012) performed parameterized simulations for a 15 M_☉ progenitor model and showed that there is no evidence of vigorous sloshing motions in 3D. Note that $\dot{M}$ is rather low in their models. Hanke et al. (2013) demonstrated, by their own 3D simulations for 25 and 27 M_☉ progenitors, that higher mass accretion rates lead to violent spiral motions, which are quenched when the accretion rate drops drastically with the infall of the Si/SiO interface.

Looking at Figure 4 in more detail, we recognize the following. First, the indeterminacy of flow pattern is observed only for Model A, in which both L_ν and $\dot{M}$ are small. This may be a numerical artifact, since we find only SP in two higher-resolution simulations, as will be shown later in Section 3.4. Considering the robustness in other cases, however, we think it is more likely that the indeterminacy is real and is yet another form of the intermediate case. Second, we find different flow patterns between the semi-nonlinear and nonlinear phases, only near the critical neutrino luminosity. Although SP always prevails in the semi-nonlinear phase, in the nonlinear phase, BB is mixed with SP to form SPB in Model B and is fully dominant in Model E, whereas another intermediate pattern SPP is realized in Model H. The characteristics of each flow pattern and the diversity in the nonlinear phase near the critical neutrino luminosity are discussed in the next section in more detail.

In this section, we discuss the dynamical features, based on which we classify the flow patterns. We first focus on the semi-nonlinear phase and then turn to the nonlinear phase.

3.2.1. Semi-Nonlinear Phase

In this paper, the semi-nonlinear phase is defined to be the period following the linear growth phase, in which there is a discernible increase of the average shock radius up to the first turning point from exponential growth to its saturation or relaxation in the time evolution of mode amplitudes. In this phase, nonlinear mode-couplings are no longer ignored, and a dominant mode emerges as a consequence. We find that one of the three basic flow patterns: SL, SP, BB, is dominant, depending on the mass accretion rate and the neutrino luminosity. We pay attention to the features in the temporal evolution of each flow pattern in the following.

Figure 5 shows the time evolutions of the mode amplitudes introduced in Section 2.2 as Equation (5). The inset in each figure is a zoom-up of the semi-nonlinear phase. We clearly see the dominance of the l = 1 mode (SL or SP) with a periodic oscillation superimposed on an exponential growth in Figures 5(a), (b), (d), (e), and (g)–(i). In Figures 5(c) and (f), on the other hand, several modes grow monotonically with similar growth rates, which is reminiscent of the non-oscillatory unstable modes that are identified as convective modes in their linear analysis by Yamasaki & Yamada (2007). We hence identify the latter two cases as BB, which is confirmed by the entropy distributions in the post-shock flows.

Zoom In Zoom Out Reset image size

In the literature, the χ parameter proposed by Foglizzo et al. (2006) is frequently used as a measure of convective instability. Since this criterion is based on the linear analysis, it is rather ambiguous for what flow the χ parameter should be calculated, except for the linear phase. For the analysis of the semi-nonlinear phase, we employ the unperturbed states. The values of χ, calculated by Equations (11) and (14), are listed in Table 2. It is found that BB is realized for χ ≳ 4 and SP or SL is observed otherwise, which seems to be consistent with the original criterion, χ > 3, although the critical value may be a little larger in our models. It is also noted that the critical χ values in our models are much larger than that found in Yamasaki & Yamada (2007), the reason for which is not clear for the moment.

SL and SP may be discriminated by the temporal evolution of the orientation of the instantaneous rotation axis for the post-shock flow. Note that, in the linear phase, only the pattern is rotating and matter has a negligible angular momentum in the spiral SASI. In the semi-nonlinear phase (and non-linear phase as well), on the other hand, matter also rotates in the same direction as the pattern (Blondin & Mezzacappa 2007). The sloshing motion does not generate angular momentum in principle. We hence expect that the instantaneous rotation axis is stable in SP with the orientation being fixed longer than the instantaneous rotation period, whereas it changes orientation stochastically in SL. Such a dichotomy is actually observed in Figure 6, where we show the orientation of the instantaneous rotation axis, which is denoted by θ and ϕ in the polar coordinate system and is calculated from the instantaneous angular momentum integrated over the post-shock flow. In panels 6(b), (d), (e), (g), (h), and (i), we see that the orientation of the instantaneous rotation axis remains almost fixed for a while and changes its direction from time to time rather suddenly, which is the feature we expect for SP. In panel 6(a), on the other hand, we observe a random walk of the rotation axis, which we suppose to characterize SL. We hence identify the flow patterns for B0, D0, E0, G0, H0, and I0 as SP and that for A0 as SL in the semi-nonlinear phase.

Zoom In Zoom Out Reset image size

3.2.2. Nonlinear Phase

After the instability grows in the linear and semi-nonlinear phases, it enters the fully nonlinear phase and is saturated. It is a non-trivial task to objectively classify various flow patterns in this phase but the ratio of the time-averaged mode amplitude of l = 1, 2 to that of l = 4, 5 may be useful, since strong sloshing and spiral motions are expected to produce large saturation amplitudes of low l modes, whereas multiple high-entropy buoyant-bubble formations will generate high l modes. Moreover, we find mixed patterns (SPP and SPB) in the nonlinear phase, which have intermediate ratios as expected. Details will follow shortly.

First, we discuss SL, sloshing pattern. High saturation amplitudes of low l modes (Figure 5(a)), extremely low angular momenta (Figure 7(a)), and unstable rotation axes (Figure 6(a)) indicate that the SL is dominant for Model A0. In fact, the ratio of $\bar{A}_{1,2}$ to $\bar{A}_{4,5}$ is pretty large, ∼10, as shown in Figure 5, in which we present the temporal evolutions of various modes. Figure 7 shows the time evolutions of the magnitude of the angular momentum integrated over the post-shock flow. Actually, the direction of sloshing changes from time to time in 3D as observed around t = 300 ms in Model A0. It is important that SL is maintained after such rather violent transitions.

Zoom In Zoom Out Reset image size

Next, we focus on SP, spiral pattern. High saturation amplitudes of low l modes (Figures 5(d) and (g)) are common with SL but high angular momenta (Figures 7(d) and (g)) as well as stable rotation axis (Figures 6(d) and (g)) characterize SP. The rotation axis shifts gradually for Model D0, and it is almost fixed during ∼1 s for Model G0. The magnitude of angular momentum in the post-shock flow is proportional to the mass accretion rate in Models D and G, for which the shock radii are almost identical. The results are consistent with the analytical study by Guilet & Fernández (2013). Model H0 is a peculiar case. The magnitude of angular momentum changes periodically, i.e., rapid and slow rotations occur alternatively with a period of ∼150 ms (Figure 7(h)). The orientation of the rotation axis varies stochastically at the local minimum of the angular momentum (Figures 6(h) and 7(h)). We hence refer to this flow pattern as SPP, the spiral motions with pulsating change of rotational velocities. Reflecting this situation, the amplitudes of all modes, including l = 0, oscillate in phase (Figures 2(c) and 5(h)). Such large time-variations of l = 2 mode may give additional features to the gravitational wave signal. Interestingly, Fernández & Thompson (2009) reported a similar phenomenon in their axisymmetric 2D simulations of post-shock flows, which lack spiral motions. In fact, they found intermittent sloshing motions, depending on the dissociation energy. This might give us a clue to the understanding of SPP.

We now turn to BB, buoyant-bubble formations. Model E0 is prototypical. After the amplitude of the dominant spiral mode reaches ∼10% of the shock radius, the flow pattern changes from SP to BB and the amplitudes of various modes are reduced to ∼4% (Figure 5(e)). In this nonlinear phase, multiple buoyant-bubbles are formed repeatedly and collide with the shock wave randomly. As a result, large-scale deformations of the shock wave, as we see them in SL or SP, are suppressed in BB. This is most clearly indicated by the ratio of $\bar{A}_{1,2}$ to $\bar{A}_{4,5}$, which is much smaller than those for SL or SP. Moreover, buoyant bubbles disrupt the regular circulations around PNS. As a result, the angular momentum is not so high (Figure 7(e)), and the orientation of rotation axis changes rapidly in a short period (Figure 6(e)). This unstable behavior of the rotation axis is another characteristic that we adopt in judging the flow pattern as BB. Our previous 3D simulations did not generate strong spiral motions as found by Blondin & Mezzacappa (2007), and the saturation amplitudes, a few percent, were much smaller than theirs. The reason for this apparent discrepancy is now clear: we worked in the BB regime.

In applying the criterion for convection with the χ parameter to the nonlinear phase of non-explosion models, we employ the time- and angle-averaged flows that are evaluated according to Equations (15) and (16) at t_s = 0.4 s and t_e = 1.0s. In Figure 8, we show the entropy $\bar{S}$, the net cooling rate $\bar{Q}_\nu$, and the ratio of the Brunt–Väisälä frequency to the radial velocity, $|\bar{N}/\bar{u}_r|$, as a function of radius for Models A0, B0, D0, E0, G0, and I0 in Figure 8. The average gain radius $\bar{r}_{\rm gain}$ is located at $\bar{Q}_\nu =0$, and the maximum shock radius $\bar{r}_{\rm sh}$ is defined to be the radius with $\bar{S}= 3.1$. One observes that a larger portion of the gain region has rather flat entropy distributions in BB and SPB than in other flow patterns. This region corresponds to convective activities. One also finds that the ratio of the Brunt-Väisälä frequency to the radial velocity becomes large in this region. As a consequence, $\bar{\chi }$ obtained by the integration of this ratio over the gain region is greater for BB than for SP/SL. The numerical values are listed in Table 2. It is recognized that $\bar{\chi }$ for SPB is a bit larger than that for BB, which is at odds with the expectation. This is likely due to insufficient grid resolutions. As a matter of fact, $\bar{\chi }$ for BB tends to increase with the spatial resolution, as shown in Table 3, and will be discussed later in Section 3.4 in more in detail.

Zoom In Zoom Out Reset image size

Table 3. Flow Patterns and Physical Parameters in Higher Resolution

Model	Pattern	Pattern	$\bar{\chi }$	$\bar{A}_{1,2}/\bar{A}_{4,5}$
	(Semi-Nonlinear)	(Nonlinear)
A0-3	SL/SP	SL/SP	2.3–2.4/2.2	9.5–9.6/7.3
A4	SP	SP	2.4	7.2
A5	SP	SP	2.3	7.4
B0-3	SP	SPB	3.7–3.8	4.9–5.3
B4	SP	SPB	3.7	4.4
B5	SP	SPB	3.7	3.6
D0-3	SP	SP	1.9–2.3	6.6–6.8
D4	SP	SP	2.0	6.8
D5	SP	SP	1.9	6.8
E0-3	SP	BB	3.5–3.6	2.8–3.2
E4	SP	BB	3.8	2.1
E5	SP	BB	4.1	1.8
G0-3	SP	SP	1.5–1.8	7.4–7.5
G4	SP	SP	1.8	7.4
G5	SP	SP	1.5	7.5
H0-3	SP	SPP	2.4–2.6	5.9–6.1
H4	SP	SPP	2.3	6.3
H5	SP	SPP	2.4	5.3

Notes. Models A0–3, A4, and A5 correspond to the results of 300 × 30 × 60, 300 × 50 × 100, and 300 × 60 × 120 mesh points, respectively, and other models are also named in a similar way. In the normal resolution models, the ranges of the χ parameter taken among different realizations are listed in this table. The meanings of symbols are referred in Table 2.

Download table as:  ASCIITypeset image

Next, we consider the intermediate pattern between SP and BB, that is, SPB. Model B0 is the prototype. The ratio of $\bar{A}_{1, 2}$ to $\bar{A}_{4, 5}$ for this model lies between the values for SP and BB as expected (Figure 5(b)). When the angular momentum is relatively large (Figure 7(b)), the rotation axis gradually changes its orientation as in SP for most of time (Figure 6(b)). When the angular momentum becomes small, as observed at t ∼ 0.6 s and t ∼ 0.8–1.0 s (Figure 7(b)), however, the rotation axis is turned as in BB (Figure 6(b)). In these periods, spiral features disappear and convective ones emerge instead. We hence classify this flow pattern in Model B0 as SPB.

Finally, we discuss the explosion models: C0, F0, and I0. In these models it takes the shock wave only a few hundred milliseconds to reach the outer boundary of the computational domain, which is located at r ∼ 500–1000 km. We may extract some more information from the late evolutions of these models. The low l modes tend to have greater amplitudes and longer oscillation periods in the explosion models than in the non-explosion models (Figures 5(c), (f) and (i)). This feature was reported in other works (e.g., Burrows 2012). The angular momentum generated in the spiral modes tend to be larger as the shock radius increases in the explosion models. It is also positively correlated with the mass accretion rate (Figures 7(c), (f), and (i)). These are consistent with the results in Guilet & Fernández (2013). The amplification of angular momentum in the nonlinear phase was also found in Blondin & Mezzacappa (2007). The unstable behavior of the rotation axis is observed in all of the explosion models and is slower in Model I0 than in Models C0 and F0 (Figures 6(c), (f), and (i)).

In this section, we discuss the reproducibility of the flow patterns, imposing different random perturbations on the initial flows. Since the post-shock flow is highly stochastic, it is non-trivial whether the flow patterns we have discussed are robust or not. Three more realizations are presented for each model here. Taking Model A as an example, we name these additional Models as A1, A2, and A3, while the original model is referred to as A0. Figures 9–17 show the results of these additional realizations for Models A−I. In the figures, the time evolutions of the mode amplitudes, the orientation, and magnitudes of angular momentum are displayed in the upper, middle, and lower panels, respectively. The different realizations reproduce the same flow patterns except for a few cases, which will be discussed in the following.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 9 indicates that there are two flow patterns, SL and SP, which might be possible in Model A. In fact, Model A1 has essentially the same features as Model A0: the oscillatory exponential growths in the linear and semi-nonlinear phases and the higher saturation amplitudes of l = 1, 2 modes (Figures 5(a) and 9(a)), the unstable behavior of the rotational axis (Figures 6(a) and 9(d)), and the lower magnitude of angular momentum (Figures 7(a) and 9(g)). Model A1 is hence classified as SL. We find, on the other hand, that Model A3 has clear SP features: the gradual change of the rotational axis (Figure 9(f)) and the higher magnitude of angular momentum (Figure 9(i)). Model A2 occupies a place in between: the SL in the semi-nonlinear and nonlinear phases until around 300 ms turn into SP after that (Figures 9(b), (e), and (h)). It is thus evident that small differences in the initial perturbations lead to the bifurcation of the flow pattern in the semi-nonlinear and nonlinear phases. It seems that SP is a little bit more stable than SL. In fact, we add more seven realizations for this model and find that SL appears in four out of the total 11 realizations. It is also noted that the direction of the sloshing motion tends to be aligned with the z-axis, which may be a numerical artifact of the spherical coordinates.

Figure 16 shows different realizations of Model H, in which we consistently observe that SP in the semi-nonlinear phase is changed to SPP in the nonlinear phase. The stability of the rotational axis is different among H0, H1, H2, and H3, though. The rotational axis is almost fixed in H1 (Figure 16(d)), whereas its orientation is mildly varied for H2 (Figure 16(e)), and drastically changed for H0 and H3 (Figures 6(h) and 16(f)). The time evolutions of the mode amplitude and angular momentum indicate that large fluctuations of their values are observed as a common feature of these models, while local minimum values of angular momentum tend to be larger for Models H1 and H2 (Figures 16(g) and (h)) than for Models H0 and H3 (Figures 7(b) and 16(i)), which correspond to the less-frequent changes of the rotational axis for Models H1 and H2. As mentioned in previously, the orientation of the rotational axis changes when the rotation almost stops. The local minimum values of angular momentum are rather stochastic.

For the rest of the models, Figures 10–15, and 17 show that the qualitative features of the flow patterns are essentially unchanged among different realizations. As an example, we take a look at Model G in Figure 15. SP always emerges with an almost fixed rotation axis both in the semi-nonlinear and nonlinear phases (Figures 6(g) and 15(d)–(f)). The angular momentum has similar magnitudes (Figures 7(g) and 15(g)–(i)). The orientation of the rotational axis is random, on the other hand, which is a consequence of the intrinsically stochastic nature of the post-shock flow. Similar things can be said irrespective of their particular flow patterns for other models (B, C, D, E, F, and I): the flow pattern itself is identical among different realizations. Some properties, such as the direction of SL/SP and the instantaneous shock geometry of BB, have stochasticity. In the explosion models, in particular, the time of the shock wave arrival at the outer boundary (Table 1) and the magnitude of angular momentum behind the shock wave at that time (Figures 17(g)–(i)) are quite different from realization to realization.

For the study of turbulent flows, high enough numerical resolutions are critically important. In this section, we discuss the dependence of the results on the space and time resolutions. So far we have used a mesh having 300 × 30 × 60 cells on the spherical coordinate system and adopted the Courant number of 0.5. Because of the limitation of numerical resources available to us we add two more models with higher resolutions: one with 300 × 50 × 100 mesh points and the other model with 300 × 60 × 120 mesh points; the Courant number is set to 0.2 in the former model and 0.1 for the latter. We refer to these two cases by attaching 4 and 5, respectively, to each model name.

The flow patterns, average χ parameter $\bar{\chi }$, and the ratio of average mode amplitudes $\bar{A}_{1,2}/\bar{A}_{3,4}$ are listed in Table 3. We found that all of the flow patterns obtained with the higher resolution are the same as the ones with the normal resolution both in the semi-nonlinear and nonlinear phases. The values of $\bar{\chi }$ and $\bar{A}_{1,2}/\bar{A}_{3,4}$ for SL, SP, and SPP are also consistent. In Figures 18, 19, and 20, we show the time evolutions of the mode amplitudes, orientations of the rotation axis, and magnitudes of angular momentum, respectively. Their qualitative features are unchanged with the high resolutions except for Model A, in which the flow pattern appearing in the higher resolution computations, Models A4 and A5, is not SL but always SP as in A3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

However, there are some differences between the normal and high resolutions. One is the value of $\bar{A}_{1,2}/\bar{A}_{4,5}$ in the buoyant-bubble formations. In fact, as the number of angular mesh points increases, the absolute values of the ratio $\bar{A}_{1,2}/\bar{A}_{4,5}$ decrease only in SPB and BB, as shown in Table 3. The angle-averaged χ parameters $\bar{\chi }$ also get larger in BB, which may suggest that the parasitic instability might be better computed in the higher resolution models (Guilet et al. 2010). In fact, the capture of small structures generated by Kelvin–Helmholtz and/or Rayleigh–Taylor instabilities will require high spatial resolutions. This is consistent with the just-mentioned fact that the ratio $\bar{A}_{1,2}/\bar{A}_{4,5}$ tends to decrease as the angular mesh points are increased. Note however, that this will pose no problem in the identification of the flow pattern. In fact, the linear analysis for $\dot{M}=1.0\,M_\odot$ demonstrated that modes with l ∼ 6 are the most unstable (Yamasaki & Yamada 2007), which are still enough to be resolved even with the ordinary resolution of 300 × 30 × 60 mesh. Although this is based on the linear analysis, the dominant unstable modes will not be much different for the mean flows.

Another difference we find in the high-resolution simulations is the onset time of the semi-nonlinear phase. The semi-nonlinear phase is the region that the specific modes are selected and exponentially grow, while all modes induced by random perturbation grow in linear phase. In some models, it is found that the onset of the semi-nonlinear phase is delayed in the high-resolution simulations. The nonlinear phase that follows the semi-nonlinear phase also commences later, as is clearly shown in Figures 18(d), (g), and (i). The end of the linear phase is most easily recognized as SP in Figures 19(a), (b), (d), (e), (g)–(i). In fact, the time, at which θ and ϕ start to be constant, marks the beginning of the semi-nonlinear phase. It is also found that the delay of the semi-nonlinear phase is longer for the models with high mass accretion rates. Note, however, that this may be simply due to the fact that, when the resolution is increased, the shares of large-scale modes in the initial cell-by-cell random perturbations become smaller. Nevertheless, the dominant flow pattern is insensitive to these differences. The growth rate of the dominant mode in the semi-nonlinear phase is unaffected either, as seen for SP in Models A, B, D, E, H, and I, although it is slightly lower in Model G5 than in G0–G4.

It is admitted that the tests presented here are not sufficient to prove the numerical convergence. It is not easy to increase the resolution further, however, due to the numerical cost. We should emphasize that, in this paper, we are interested only in relatively large and coherent features, such as global shock oscillations and formations of large buoyant bubbles. It is a future task to understand how these coherent structures are formed from a mixture of various modes through nonlinear interactions.