3D-simulations of SASI in Supernovae

# Three-Dimensional Simulations of Standing Accretion Shock Instability in Core-Collapse Supernovae

Wakana Iwakami1 , Kei Kotake1 1 , Naofumi Ohnishi1 1 , Shoichi Yamada1 1 and Keisuke Sawada1 Department of Aerospace Engineering, Tohoku University, 6-6-01 Aramaki-Aza-Aoba, Aoba-ku, Sendai, 980-8579, Japan Division of Theoretical Astronomy, National Astronomical Observatory of Japan, 2-21-1, Osawa, Mitaka, Tokyo, 181-8588, Japan Max-Planck-Institut für Astrophysik, Karl-Schwarzshild -Str. 1, D-85741, Garching, Germany Center for Research Strategy and Support, Tohoku University, 6-6-01 Aramaki-Aza-Aoba, Aoba-ku, Sendai, 980-8579, Japan Science & Engineering, Waseda University, 3-4-1 Okubo, Shinjuku, Tokyo, 169-8555, Japan Advanced Research Institute for Science and Engineering, Waseda University, 3-4-1 Okubo, Shinjuku, Tokyo, 169-8555, Japan
1affiliationmark:
###### Abstract

We have studied non-axisymmetric standing accretion shock instability, or SASI, by 3D hydrodynamical simulations. This is an extention of our previous study on axisymmetric SASI. We have prepared a spherically symmetric and steady accretion flow through a standing shock wave onto a proto-neutron star, taking into account a realistic equation of state and neutrino heating and cooling. This unperturbed model is supposed to represent approximately the typical post-bounce phase of core-collapse supernovae. We then have added a small perturbation (1%) to the radial velocity and computed the ensuing evolutions. Not only axisymmetric but non-axisymmetric perturbations have been also imposed. We have applied mode analysis to the non-spherical deformation of the shock surface, using the spherical harmonics. We have found that (1) the growth rates of SASI are degenerate with respect to the azimuthal index of the spherical harmonics , just as expected for a spherically symmetric background, (2) nonlinear mode couplings produce only modes for the axisymmetric perturbations, whereas modes are also generated in the non-axisymmetric cases according to the selection rule for the quadratic couplings, (3) the nonlinear saturation level of each mode is lower in general for 3D than for 2D because a larger number of modes are contributing to turbulence in 3D, (4) low modes are dominant in the nonlinear phase, (5) the equi-partition is nearly established among different modes in the nonlinear phase, (6) the spectra with respect to obey power laws with a slope slightly steeper for 3D, and (7) although these features are common to the models with and without a shock revival at the end of simulation, the dominance of low modes is more remarkable in the models with a shock revival.

supernovae: collapse — neutrinos — hydrodynamics — instability
slugcomment:

## 1 Introduction

Many efforts have been made for the multi-dimensional modeling of core-collapse supernovae (see Woosley & Janka (2005); Kotake et al. (2006) for reviews), urged by accumulated observational evidences that core-collapse supernovae are globally aspherical commonly (Wang et al., 1996, 2001, 2002). Various mechanisms to produce the asymmetry have been discussed so far: convection (e.g., Herant et al. (1994); Burrows et al. (1995); Janka & Mueller (1996)), magnetic field and rapid rotation (see, e.g., Kotake et al. (2006) for collective references), standing (/stationary/spherical) accretion shock instability, or SASI, (Blondin et al., 2003; Scheck et al., 2004; Blondin & Mezzacappa, 2006; Ohnishi et al., 2006, 2007; Foglizzo et al., 2006), and g-mode oscillations of protoneutron stars (Burrows et al., 2006). Note, however, that most of them have been investigated with two-dimensional (2D) simulations.

Recently a 3D study on SASI was reported by Blondin & Mezzacappa (2007). In 2D, the shock deformation by SASI is described with the Legendre polynomials , or the spherical harmonics with . Various numerical simulations have demonstrated unequivocally that the mode is dominant and a bipolar sloshing of the standing shock wave occurs with pulsational strong expansions and contractions along the symmetry axis (Blondin et al., 2003; Scheck et al., 2004; Blondin & Mezzacappa, 2006; Ohnishi et al., 2006, 2007). In 3D, on the other hand, Blondin & Mezzacappa (2007) observed the dominance of a non-axisymmetric mode with , , which produces a single-armed spiral in the later nonlinear phase. They claimed that this ”spiral SASI” generates a rotational flow in the accretion (see also Blondin & Shaw (2007) for 2D computations in the equatorial plane) and that it may be an origin of pulsar spins.

There are many questions on 3D SASI remaining to be answered yet, however. For example, we are interested in how the growth of SASI differs between 3D and 2D among other things. In particular, the change in the saturation properties should be made clear. This will be the focus of this paper. Another issue will be the generation of rotation in the accretion flow by SASI (Blondin & Mezzacappa, 2007). Its efficiency and possible correlation with the net linear momentum should be studied more in detail and will be the subject of our forthcoming paper (Iwakami et al., 2007). Incidentally, it is noted that the neutrino heating and cooling were entirely ignored and the flow was assumed to be isentropic in Blondin & Mezzacappa (2007), but the implementation of these physics is helpful for considering the implications not only for the shock revival but also for the nucleosynthesis (Kifonidis et al., 2006).

In this paper, we have performed 3D hydrodynamic simulations, employing a realistic equation of state (Shen et al., 1998) and taking into account the heating and cooling of matter via neutrino emissions and absorptions on nucleons as done in our 2D studies (Ohnishi et al., 2006, 2007). The inclusion of the neutrino heating allows us to discuss how the critical luminosity for SASI-triggered explosion could be changed in 3D from those in 2D. To answer the questions we raised above, we have varied the initial perturbations as well as the neutrino luminosity and compared the growth of SASI between 2D and 3D in detail, conducting a mode analysis not only for the linear phase but also for the nonlinear saturation phase.

The plan of this paper is as follows: In Section 2, we describe the models and numerical methods. We show the main numerical results in Section 3. We conclude the paper in Section 4.

## 2 Numerical Models

The numerical methods employed in the present paper are based on the code ZEUS-MP/2 (Hayes et al., 2006), which is a computational fluid dynamics code for the simulation of astrophysical phenomena, parallelized by the MPI (message-passing) library. The ZEUS-MP/2 code employs the Eulerian hydrodynamics algorithms based on the finite-difference method with a staggered mesh. In this study, we have modified the original code substantially according to the prescriptions in our preceding 2D simulations (Ohnishi et al., 2006, 2007).

We consider spherical coordinates with the origin at the center of the proto-neutron star. The basic evolution equations describing accretion flows of matter attracted by a proto-neutron star and irradiated by neutrinos emitted from the proto-neutron star are written as follows,

 dρdt+ρ∇⋅\boldmathv=0, (1)
 ρd\boldmathvdt=−∇P−ρ∇Φ−∇⋅Q, (2)
 ρddt(eρ)=−P∇⋅\boldmathv+QE−Q:∇v, (3)
 dYedt=QN, (4)
 Φ=−GMinr, (5)

where , and are the density, velocity, internal energy, pressure, electron fraction, and gravitational potential, respectively. is the gravitational constant. The self-gravity of matter in the accretion flow is ignored. is the artificial viscous tensor. and represent the heating/cooling and electron source/sink via neutrino absorptions and emissions by free nucleons. The Lagrangian derivative is denoted by . The tabulated realistic equation of state based on the relativistic mean field theory (Shen et al., 1998) is implemented according to the prescription in Kotake et al. (2003). The mass accretion rate and the mass of the central object are fixed to be  s and , respectively. The neutrino heating is estimated under the assumptions that neutrinos are emitted isotropically from the central object and that the neutrino flux is not affected by local absorptions and emissions (See Ohnishi et al. (2006)). We consider only the interactions of electron-type neutrinos and anti-neutrinos. Their temperatures are also assumed to be constant and set to be  MeV and  MeV, which are the typical values in the post-bounce phase. The neutrino luminosity is varied in the range of  ergs s.

The spherical polar coordinates are adopted. In the radial direction, the computational mesh is nonuniform, while the grid points are equally spaced in other directions. We use 300 radial mesh points to cover , where is the radius of the inner boundary located roughly at the neutrino sphere and is the radius of the outer boundary, at which the flow is supersonic. 30 polar and 60 azimuthal mesh points are used to cover the whole solid angle. In order to see if this angular resolution is sufficient, we have computed a model with the mesh points and compared it to the counterpart with the mesh points. As shown later in Appendix B, the results agree reasonably well with each other both in the linear and nonlinear phases. Although the computational cost does not allow us to do the convergence test further, we think that the resolution of this study is good enough.

We use an artificial viscosity of tensor-type described in Appendix A instead of the von Neumann and Richtmyer-type that was originally employed in ZEUS-MP/2. For 3D simulations with the spherical polar mesh, we find it far better to employ the former than the latter to prevent the so-called carbuncle instability (Quirk, 1994), which we observe around the shock front near the symmetry axis, . With the original artifitial viscosity, an appropriate dissipation is not obtained in the azimuthal direction for the shear flow resulting from the converging accretion particularly when a fine mesh is used (Stone & Norman, 1992). We have applied this artificial viscosity also to axisymmetric 2D simulations and reproduced the previous results (Ohnishi et al., 2006).

Figure 1 shows the radial distributions of various variables for the unperturbed flows. The spherically symmetric steady accretion flow through a standing shock wave is prepared in the same manner as in Ohnishi et al. (2006). Behind the shock wave, the electron fraction is less than 0.5 owing to the electron capture, and a region of negative entropy gradient with positive net heating rates is formed for the neutrino luminosities,  ergs s. The values of these variables on the ghost mesh points at the outer boundary are fixed to be constant in time while on the ghost mesh points at the inner boundary they are set to be identical to those on the adjacent active mesh points, except for the radial velocity, which is fixed to the initial value both on the inner and outer boundaries.

In order to induce the non-spherical instability, we have added a radial velocity perturbation, , to the steady spherically symmetric flow according to the following equation,

 vr(r,θ,ϕ)=v1Dr(r)+δvr(θ,ϕ), (6)

where is the unperturbed radial velocity. In this study, we consider three types of perturbations: (1) an axisymmetric () single-mode perturbation,

 δvr(θ,ϕ)∝√34πcosθ⋅v1Dr(r), (7)

(2) a non-axisymmetric perturbation with ,

 δvr(r,θ,ϕ)∝[√34πcosθ−√38πsinθcosϕ]⋅v1Dr(r), (8)

and (3) a random multimode perturbation,

 δvr(θ,ϕ)∝rand⋅v1Dr(r)(0≤rand<1), (9)

where ”rand” is a pseudo-random number. The perturbation amplitude is set to be less than 1% of the unperturbed velocity. It is noted that there is no distinction between and modes as long as the initial perturbation is added only to the radial velocity, as is the case in this paper. To put it another way, the modes contribute equally. Hence they are referred to as the mode in the following. On the other hand, differences do show up when the initial perturbation is added also to the non-radial velocity components, which case is important in discussing the origin of pulsar spins proposed by Blondin & Mezzacappa (2007) and will be the subject of our forthcoming paper (Iwakami et al., 2007). All the models presented in this paper are summarized in table 1.

In the next section, we perform the mode analysis as follows. The deformation of the shock surface can be expanded as a linear combination of the spherical harmonics components :

 RS(θ,ϕ)=∞∑l=0l∑m=−lcmlYml(θ,ϕ), (10)

where is expressed by the associated Legendre polynomial and a constant given as

 Yml=KmlPml(cosθ)eimϕ, (11)
 Kml=√2l+14π(l−m)!(l+m)!. (12)

The expansion coefficients can be obtained as follows,

 cml=∫2π0dϕ∫π0dθsinθRS(θ,ϕ)Ym∗l(θ,ϕ), (13)

where the superscript * denotes complex conjugation.

## 3 Results and Discussions

### 3.1 Axisymmetric Single-Mode Perturbation

Before presenting the results of 3D simulations, we first demonstrate the validity of our newly developed 3D code. We compare the axisymmetric flows obtained by 2D and 3D simulations for the axisymmetric () single-mode perturbation. By 2D simulations we mean that axisymmetry is assumed and computations are done in the meridian section with all -derivatives being dropped, whereas in 3D simulations we retain all these derivatives and 3D mesh is employed. This validation is important because numerical errors may induce azimuthal motions even for the axisymmetric initial conditions. Hence we have to confirm that azimuthal errors do not appear in the non-axisymmetric simulation.

Figure 2 shows the time evolutions of the average shock radius . The 2D and 3D results are displayed in Figure 2 (a) and (b), respectively. The average shock radius is obtained from the expansion coefficient in Eq. (10) multiplied by . The solid, broken, and dotted lines correspond to the neutrino luminosities of and ergs s, respectively. It should be first mentioned that one cannot expect an perfect agreement between two computations of the exponential growth of the instability followed by the turbulent motions through mode couplings as considered here. It is still obvious that the results of the 2D and 3D simulations agree on essential features. In particular, both in the 2D and 3D results, we see that the is settled to an almost constant value for and ergs s, whereas it continues to increase for ergs s.

Figure 3 presents the normalized amplitudes as a function of time for Model I ( ergs s). The 2D and 3D results are displayed in Figure 3 (a) and (b), respectively. The red, blue, black, and gray solid lines correspond to the modes of and , respectively. As can be seen, the time evolution can be divided into two phases. One is the linear growth phase, in which the amplitude of the initially imposed mode with grows exponentially. It lasts for 100 ms. It is also found that higher modes are generated by mode couplings and grow exponentially during this phase. Then starts the second phase, which is characterized by the saturation of amplitudes of the order of 0.1. In this phase, the accretion flow becomes turbulent. It is interesting to note that in this nonlinear saturation phase the amplitude of the mode with is almost of the same order as the initially imposed mode with and is dominant over other modes, which fact was observed in Ohnishi et al. (2006).

Most important of all, however, is the fact that all the modes are not produced, implying that the perturbed flow retains axisymmetry, which is a necessary condition for the numerical study of non-axisymmetric instability. Although the results of the 2D and 3D simulations are not identical again, the essential features such as the linear growth rate of the initial perturbation () and the production of other modes via nonlinear mode couplings as well as the saturation levels in the nonlinear phase are all in good agreement for the two cases. The quantitative differences between the 2D and 3D results mainly originate from the difference in the time steps, which depend on the minimum grid width. Note that in addition to the 300 radial and 30 polar mesh points used for the 2D computations, 60 azimuthal mesh points are employed in the 3D simulatioins and, as a result, smaller time steps are usually taken for the 3D computations.

### 3.2 Non-Axisymmetric Perturbation with l=1

Now we discuss the results of genuinely 3D simulations, in which the non-axisymmetric perturbation with () is added to the axisymmetric () perturbation.

Figure 4 shows the time evolutions of some of the normalized amplitudes for Model II with the neutrino luminosity ergs s. The red, yellow, blue, light blue, green, black, brown, violet, and pink solid lines denote the amplitudes of the modes with and , respectively. In the linear phase, the initially imposed modes with and grow exponentially just as in the axisymmetric single-mode perturbation. It should be noted that the growth rate of mode is identical to that of mode. This is just as expected for the spherically symmetric background, for which the modes with the same but different are degenerate in the linear eigen frequency.

It is observed in the figure that the modes with are first produced by the nonlinear mode-couplings of . Then the modes with can be produced by the couplings of the first-order modes with and the second-order modes with . Even higher order modes are produced subsequently toward the non-linear saturation but are not shown here. Although the branching ratios should be investigated more in detail before giving any conclusion, the above sequence of the mode generations strongly suggests that the nonlinear coupling is mainly of quadratic nature.

In the nonlinear phase that begins at ms, these mode amplitudes are saturated and settled into a quasi-steady state. As in the axisymmetric case, the modes both with and are dominant over other modes in this stage for the non-axisymmetric case. One thing to be noted here, however, is the fact that the saturated amplitudes for Model II are lower than those for Model I in general (see Figure 3 (b)). This is, we think, because the turbulent energy is shared by a larger number of modes in the non-axisymmetric case than in the axistmmetric case. To demonstrate this more clearly, we have added the random perturbation to the radial velocity in the axisymmetric Model I at ms, by which time the axisymmetric nonlinear turbulence has fully grown. We refer to this model as Model III. The entire time evolutions of the normalized amplitudes for Model III are shown in Figure 5. It is clear that the axisymmetric mode amplitudes are reduced after the additional perturbation is given and the non-axisymmetric modes grow to the saturation level at ms. The power spectra of the turbulent motions will be discussed more in detail later.

### 3.3 Random Multi-Mode Perturbations

#### 3.3.1 Dynamical Features and Critical Luminosity

Having understood the elementary processes of the linear growth and nonlinear mode couplings in the previous section, we will now discuss the 3D SASI induced by the random multimode perturbations, which are supposed to be closer to the reality. In this subsection, we will show the typical dynamics, paying attention to the time evolution of the shock wave and, hence, to the critical luminosity, at which the stalled shock is revived.

Figure 6 shows the time variations of the average shock radius for Models IV, V and VI with the neutrino luminosities, and ergs s, respectively. It is found that for and ergs s, the shock is settled to a quasi-steady state after the linear growth whereas it continues to expand for ergs s, which is also true for the axisymmetric counterpart displayed in the right panel of the figure for comparison. This implies that the critical neutrino luminosity is not much affected by the change from 2D to 3D. In the following we refer to the models with and ergs s as the non-explosion models and the model with ergs s as the explosion model and look into their dynamical features in turn. The mode analyses will be done in sections 3.3.2 and 3.3.3.

Figure 7 shows the sideviews of the iso-entropy surfaces together with the velocity vectors in the meridian section at four different times for the non-explosion model with ergs s (Model IV). The hemispheres () of eight iso-entropy surfaces are superimposed on one another and the outermost surface almost corresponds to the shock front. The initial perturbations grow exponentially in the linear phase ((a) ms). At the end of the linear phase, a blob of high entropy is formed and grows subsequently ((b) ms). High entropy blobs are continuously generated and the nonlinear phase begins ((c) ms). Circulating flows occur inside the blobs while high-velocity accretions onto a proto-neutron star surround the blobs. These blobs expand and shrink repeatedly, being distorted, splitted and merged with each other inside the shock wave ((d) ms). Reflecting these complex motions, the shock surface oscillates in all directions but remains alomost spherical for the non-axisymmetric model. This is in sharp contrast to the axisymmetric case, in which large-amplitude oscillations occur mainly in the direction of the symmetry axis (Blondin et al. (2003); Ohnishi et al. (2006)).

Figure 8 displays the snapshots of the iso-entropy surfaces and the velocity vectors for the explosion model with ergs s (Model VI). As in the non-explosion model, the sequence of events starts with the linear growth of the initial perturbation inside the shock wave ((a) ms). In the explosion model, however, many high entropy blobs emerge simultaneously much earlier on ((b) ms) than in the non-explosion model. Then, these blobs repeat expansions and contractions, being distorted, splitted and merged together just as in the non-explosion counterpart ((c) ms). After that, some of the blobs get bigger, absorbing other blobs ((d) ms) and, as a result, the shock radius increases almost continuously up to the end of the computation ( ms) as already demonstrated in Figure 6 (a). At this point of time, the shock surface looks ellipsoidal rather than spherical. However, the major axis is not necessarily aligned with the symmetry axis and the flow is not symmetric with respect to this major axis, either.

#### 3.3.2 Mode Spectra

Next we look into the spectral intensity more in detail. As a standard case, we take Model IV with ergs s. Figure 9 shows the time evolutions of the normalized amplitudes with and compares them with the axisymmetric counterparts in Model VII. Note that the modes also exist in the non-axisymmetric model, which are not shown in the figure but will be discussed in the following paragraphs. As can be clearly seen, the amplitude of each mode grows exponentially until 100 ms, which is the linear phase. Note in particular that the growth rate and oscillation frequency for the mode are the same as those obtained for the model with the single mode perturbation in the previous section. After 100 ms, the mode amplitudes are saturated and the evolution enters the nonlinear phase with a clear dominance of the modes with . It is also evident in the figure that the saturation level is lower in general for the non-axisymmetric case than for the axisymmetric one, which confirms the claim in the previous section based on the results for the single-mode perturbation.

In Figure 10, we display the snapshots at four different times of the spectral distributions both for the non-axisymmetric (the left panels) and axisymmetric (the right panels) cases. The uppermost panels correspond to the linear phase and the intensity is distributed rather uniformly over all modes. As the time passes, however, the amplitudes of low modes grow much more rapidly than those with higher ’s (the second and third rows in the figure). One can see a similarity in the time evolutions between the non-axisymmetric and axisymmetric models. After the nonlinear phase starts (the lowermost panels), the growths of all modes are saturated and the spectral distributions are settled to be quasi-steady. It is again obvious in these panels that the low modes are dominant in the nonlinear phase and the saturation level is lower in the non-axisymmetric case.

Now we turn our attention to the quasi-steady turbulence in the nonlinear phase. Shown as a function of in Figure 11 are the power spectra averaged over the interval from  ms to 400 ms. Both the non-axisymmetric and axisymmetric cases are displayed by the symbols, and , respectively. In the left panel, modes with different ’s are plotted separately whereas they are summed up in the right panel. It is found that the time-averaged power spectra are not much different among the modes with the same but with a different . This implies that the equi-partition of the turbulence energy is nearly established among the modes with the same . This will have an important ramification for the origin of pulsar spin and will be discussed in our forthcoming paper (Iwakami et al., 2007).

The right panel of the figure demonstrates that the time-averaged power spectrum summed over for the non-axisymmetric case is not much different from that for the axisymmetric counterpart. This means that the total turbulence energy does not differ very much between the axisymmetric and non-axisymmetric cases. As a result, the power of each mode in the non-axisymmetric case is smaller roughly by a factor of than that of the same mode in the axsymmetric case. This is the reason why we have observed that the saturation level of the nonlinear SASI is smaller in the non-axisymmetric case than in the axisymmetric case and the shock surface oscillates in all directions with smaller amplitudes in the non-axisymmetric flow whereas it sloshes in the direction of the symmetry axis with larger amplitudes in the axisymmetric case. The difference in the fluctuations of the average shock radius seen in Figure 6 (a) and (b) is also explained in the same manner.

Another interesting feature observed in the figure is the fact that the time-averaged power spectra obey a power-law at for both the non-axisymmetric and axisymmetric models. Two straight lines in Figure 11 (a) are the fits to the data for , each obtained for the axisimmetric and non-axisymmetric models. The powers are found to be and for the non-axisymmetric and axisymmetric cases, respectively. Although the difference is almost unity, which originates from the difference in multiplicity of the modes with the same , the slope excluding this effect is still a little bit steeper in the non-axisymmetric case as seen in Figure 11 (b). At the moment we do not know if this difference in slope is real or not and the power itself is also remaining to be explained theoretically.

#### 3.3.3 The Explosion Models

So far we have been discussing the non-explosion models, in which the SASI is saturated in the non-linear phase and is settled to a quasi-steady state. In considering the possible consequences of SASI after the supernova explosion occurs, however, it is also interesting to investigate the explosion models, in which a shock revival occurs as a result of SASI.

Figure 12 shows the time evolutions of the normalized amplitudes for the explosion model (Model VI) and compares it with the axisymmetric counterpart (Model VIII). The feature of the explosion models common to the axisymmetric and non-axisymmetric cases is as follows: The nonlinear stage is devided into two phases. The earlier phase is quite similar to the nonlinear phase that we have seen so far for the non-explosion models. The later phase, on the other hand, has a distinctive feature that the mode becomes much more remarkable than other modes and the oscillation period tends to be longer as the shock radius gets larger and the advection-acoustic cycle, supposedly an underlying mechanism of the instability, takes longer. The later prominence of mode is intriguing and is important to discuss the kick velocity of pulsar quantitatively. The theoretical account, however, is yet to be given. The differences between the non-axisymmetric and axisymmetric models, such as the saturation level, are just taken over to the explosion models.

The power spectra averaged over the time intervals of  ms and  ms are presented in Figure 13 (a) and (b), respectively. As in the non-explosion models, the equi-partition among different modes is again established for the explosion models. This is true even in the later nonlinear phase as seen in the right panel of the figure. As a result, the saturation level is smaller for the non-axisymmetric case than for the axisymmetric case as mentioned just above. The power spectra in the earlier nonlinear phase look much similar to those found in the non-explosion models, with the power law being satisfied for . In the later nonlinear phase, on the other hand, the power law is extended to much lower ’s. This is related to the late-time prominence of the mode as mentioned above. In fact, the amplitudes in the lower portion of the power spectrum are enhanced in . The mechanism of this amplification is remaining to be understood but might be related with the volume-filling thermal convection advocated for the late stage of the convective instability in the supernova core by Kifonidis et al. (2006). As a result of this evolution of the spectrum, the shock tends to be ellipsoidal as the shock expands.

### 3.4 Neutrino Heating

The SASI is supposed to be an important ingredient not only for the pulsar’s proper motions and spins but also for the explosion mechanism itself. In this section, we look into the neutrino heating in the non-axisymmetric SASI.

We show in Figure 14 the color contours of the net heating rate in the meridian section for the non-explosion model with ergs s (non-axisymmetric Model IV and axisymmetric Model VII) as well as for the non-axisymmetric explosion model with ergs s (Model VI ). In the early phase, the cooling region with negative net heating rates exists around the proto neutron star while the heating region extends over the cooling region up to the stalled shock wave in all the cases. As the time passes and the SASI grows, a pocket of regions with positive but relatively low net heating rates appear. These regions correspond to the high entropy blobs (high entropy rings for the axisymmetric case), where the circulating flow exists observed in Figure 7 (b), (c), and (d). Since the neutrino emission in these regions is more efficient than the surroundings, the net heating rate is rather low.

Although the critical neutrino luminosties are not much different between the non-axisymmetric and axisymmetric cases, the spatial distributions of neutrino heating are different. In the axisymmetric case the shock wave oscillates up and down whereas it moves in all directions in the non-axisymmetric case. The oscillation amplitudes are larger in the former than in the latter in general as repeated mentioned. Reflecting this difference in the shock motions, the neutrino heating is enhanced in mainly in the polar regions in the axisymmetric case while in the non-axisymmetric case the heating rate is affected by SASI chiefly around the high entropy blobs. These are clearly seen in Figure 14 (see also Figure 7).

As described in section 3.3.1, the generation of the high entropy blobs is started around the end of the linear phase and continues during the nonlinear phase. Although the blobs repeat expansions and contractions, being sometimes merged or splitted, the turbulent motions together with the neutrino heating are finally settled to a quasi-steady state in the non-explosion model. For the explosion model, on the other hand, the number and volume of high entropy blobs are increased much more quickly and, as a result, the heating region prevails, pushing the shock wave outwards and narrowing the cooling region near the proto neutron star.

Figure 15 shows the time evolutions of the net heating rates integrated over the region inside the shock wave. For the non-explosion model (Model IV, solid line), the volume-integrated heating rate grows until ms but is saturated thereafter, whereas it increases continuously for the explosion model (Model VI, dashed-dotted line). These behaviors of the volume-integrated heating rate are in accord with the time evolutions of the average shock radius as shown in Figure 6 (a). It is clear from the figure that the heating rates are mainly affected by SASI during the nonlinear phase. For comparison, the spherically symmetric counterparts, Models IX with ergs s and X with ergs s are also presented as dashed and dotted lines, respectively. Note that the last model does not lead to the shock revival even with this high neutrino luminosity. It is found that the SASI enhances the neutrino heating both for the explosion and non-explosion models.

## 4 Conclusion

In this paper we have studied the non-axisymmetric SASI by 3D hydrodynamical simulations, taking into account the realistic EOS and neutrino-heating and cooling. We have added various non-axisymmetric perturbations to spherically symmetric steady flows that accrete through a standing shock wave onto a proto-neutron star, being irradiated by neutrinos emitted from the proto-neutron star. The mode analysis has been done for the deformation of the shock surface by the spherical harmonics expansion. After confirming that our 3D code is able to reproduce for the axisymmetric perturbations the previous results on 2D SASI that we obtained in Ohnishi et al. (2006), we have done genuinely 3D simulations and found the followings.

First, the model of the initial perturbation with and have demonstrated that the non-axisymmetric SASI proceeds much in the same manner as the axisymmetric SASI: the linear phase, in which the initial perturbation grows exponentially, lasts for about ms and is followed by the nonlinear phase, in which various modes are produced by nonlinear mode couplings and their amplitudes are saturated. It has been found that the critical neutrino luminosity, above which the shock revival occurs, is not much different between 2D and 3D. For the neutrino luminosities lower than the critical value, the SASI is settled to a quasi-steady turbulence. We have found that the saturation level of each mode in the non-axisymmetric SASI is lower in general than that of the axisymmetric counterpart. This is mainly due to the fact that the number of the modes contributing to the turbulence is larger for the non-axisymmetric case. The sequence of the mode generation, on the other hand, strongly suggests that the nonlinear mode coupling is chiefly of quadratic nature.

Second, the simulations with the random multi-mode perturbations being imposed initially have shown that the dynamics in the linear phase is essentially a superimposition of those of single-modes. Toward the end of the linear phase, high entropy blobs are generated continuously and grow, starting the nonlinear phase. We have observed that these blobs repeat expansions and contrations, being merged and splitted from time to time. In the non-explosion models, these nonlinear dynamics lead to the saturation of mode amplitudes and the quasi-steady turbulence. For the explosion models, on the other hand, the production of blobs proceeds much more quickly, followed by the oligarchic evolution, with a relatively small number of great blobs absorbing smaller ones, and as a result the shock radius increases almost monotonically. The spectral analysis has clearly demonstrated that low modes are predominant in the nonlinear phase just as in the axisymmetric case. We have also shown that the equi-partition is nearly established among different modes in the quasi-steady turbulence and that the spectrum summed over in the non-axisymmetric case is quite similar to the axisymmetric counterpart. This implies that the larger number of modes in the non-axisymmetric case is the main reason why the amplitude of each mode is smaller in 3D than in 2D. The power spectrum is approximated by the power law for . Although the slope is slightly steeper for the non-axisymmetric models, whether the difference is significant or not is unknown at present.

We have seen in the explosion models, on the other hand, that the oscillation period of each mode becomes longer in the late nonlinear phase as the shock radius gets larger and the advection-acoustic cycle becomes longer. What is more interesting is the fact that in this late phase the dominance of low modes becomes even more remarkable. Although this may be related with the volume-filling thermal convection, the mechanism is yet to be revealed. This spectral evolution leads to the global deformation of the shock surface to an ellipsoidal shape, whose major axis is not necessarily aligned with the symmetry axis, however.

Finally, we have presented the neutrino heating in the 3D SASI. It has been seen that the volume-integrated heating rate is affected mainly in the nonlinear phase. The comparison with the spherically symmetric counterparts has confirmed that the SASI enhances the neutrino heating also in the non-axisymmetric case. Although the critical neutrino luminosity in the non-axisymmetric SASI is not much changed from that for the axisymmetric case, the spatial distribution of neutrino heating is different in the non-linear phase. Relatively narrow regions surrounding high entropy blobs are efficiently heated for the non-axisymmetric case, while wider regions near the symmetry axis are heated strongly in accord with the sloshing motion of shock wave along the symmetry axis for the axisymmetric case. For the non-explosion models, the high entropy blobs produced by neutrino heatings occupy the inside of shock wave, repeating expansion and contraction and being splitted and merged intermittently, but the flow is finally settled to a quasi-static state. For the explosion models, on the other hand, the high entropy blobs are generated much more quickly and extend the heating region, pushing the shock outwards and narrowing the cooling region near the proto-neutron star.

In the present study we have not observed a persistent segregation of angular momentum in the accretion flow such as found by Blondin & Mezzacappa (2007); Blondin & Shaw (2007)) although instantaneous spiral flows are frequently seen. As discussed above, the equi-partition is nearly established among different modes in our models. It should be emphasized here again, however, that we have added the initial perturbations only to the radial velocity in this study and, as a result, the modes with are equally contributing. We will defer the analysis of the models with the initial non-axisymmetric perturbations being added also to the azimuthal component of velocity to our forthcoming paper (Iwakami et al., 2007), in which we will also discuss a possible correlation between the kick velocity and spin of neutron stars if they are produced by the 3D SASI indeed. A lot of issues on SASI are still remaining to be studied. In particular, the influences of rotation and magnetic field are among the top priorities and will be addressed elsewhere in the near future.

One of the authors (WI) expresses her sincere gratitude to assoc. prof. Kazuyuki Ueno, assistant prof. Michiko Furudate, and the member of Sawada and Ueno laboratory for continuing encouragements and advices. She would like to thank also Yoshitaka Nakashima of Tohoku University for his advices on the visuallization by AVS. KK expresses his thanks to Katsuhiko Sato for his continuous encouragements. Numerical cumputations were performed on the Altix3700Bx2 at the Institute of Fluid Science, Tohoku University as well as on VPP5000 and the general common use computer system at the center for Computational Astrophysics, CfCA, the National Astronomical Observatory of Japan. This study was partially supported by Program for Improvement of Research Environment for Young Researchers from Special Coordination Funds for Promoting Science and Technology (SCF), the Grants-in-Aid for the Scientific Research (No.S19104006, No.S14102004, No.14079202, No.14740166) and Grant-in-Aid for the 21st century COE program “Holistic Research and Education Center for Physics of Self-organizing Systems” of Waseda University by the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan.

## Appendix A Tensor Artificial Viscosity

We present the tensor-type artificial viscosity, , that we use in this paper. It is a direct extention to 3D of the one employed by Stone & Norman (1992) for 2D simulations. Following the notations in Stone & Norman (1992), is written as follows:

 Q={l2ρ∇⋅v[∇v−13(∇⋅v)e]if   ∇⋅v<00otherwise, (A1)

where denotes a constant with a dimension of length, is the symmetrized velocity-gradient tensor, is the unit tensor. Dropping the off-diagonal components, we utilize only the diagonal components of in this study, which are written in the finite difference form as

 (A2)

where stands for the density at the site specified by and is written as

 li,j,k=max(C2 dx1ai,C2 dx2aj,C2 dx3ak). (A3)

Here is a dimensionless constant controlling the number of grid points, over which a shock is spread, and , , and are the grid widths of the “a-mesh” at the -, -, and -th grid points in the , , and directions, respectively. The “a-mesh” and “b-mesh” are distinguished in ZEUS-MP/2 and are defined on the cell-edge and cell-center, respectively (see Stone & Norman (1992) for more details).

The artificial viscous stress and artificial dissipation in the momentum equation (2) and energy equation (3) are given as follows,

 (∇⋅Q)(1)=1g22g31∂∂x1(g22g31Q11),(∇⋅Q)(2)=1g2g232∂∂x2(g232Q22)+Q11g2g32∂g32∂x2,(∇⋅Q)(3)=1g2g32∂Q33∂x3, (A4)
 Q:∇v=l2ρ∇⋅v13[(∇v(11)−∇v(22))2+(∇v(11)−∇v(33))2+(∇v(33)−∇v(22))2], (A5)

where , , and are the metric components for the spherical coordinates, (see again Stone & Norman (1992)). These equations are discretized as

 (∇⋅Q)1,i,j,k=g2b2ig31biQ11i,j,k−g2b2i−1g31bi−1Q11i−1,j,kg2a2ig31aidx1bi,(∇⋅Q)2,i,j,k=g32b2jQ22i,j,k−g32b2j−1Q22i,j−1,kg2big32a2idx2bj+Q11i,j,k+Q11i,j−1,k2g2big32aj(∂g32aj∂x2),(∇⋅Q)3,i,j,k=Q33i,j,k−Q33i,j,k−1g2bi,g32bjdx3bk, (A6)
 (Q:∇v)i,j,k=l2i,j,kdi,j,k(∇⋅v)i,j,k13[[(∇v(11))i,j,k−(∇v(22))i,j,k]2+[(∇v(11))i,j,k−(∇v(33))i,j,k]2+[(∇v(33))i,j,k−(∇v(22))i,j,k]2], (A7)

where , , and are , , and defined on the“a-mesh”, whereas , , and are those defined on the“b-mesh”. , , and represent the width of the “b-mesh” at the -, -, and -th grid points in the , , and directions, respectively. Finally the velocity-gradient tensor , and are given by the following equations:

 (∇v(11))i,j,k=v1i+1,j,k−v1i,j,kdx1ai,(∇v(22))i,j,k=v2i,j+1,k−v2i,j,kg2bidx2aj+v1i,j,k+v1i+1,j,k2g2bi(∂g2bi∂x1),(∇v(33))i,j,k=v3i,j,k+1−v3i,j,kg2big32bjdx3ak+v2i,j,k+v2i,j+1,k2g31big32bj(∂g32bj∂x2)+v1i,j,k+v1i+1,j,k2g31bi(∂g31bi∂x1). (A8)

## Appendix B Numerical Convergence Tests

In order to see if the numerical resolution employed in the main body is sufficient, we increase the number of angular grid points and compare the results. Figure 16 (a) shows the time evolutions of the normalized amplitudes in the linear phase. In this comparison, we impose the perturbation initially. We refer to the model with mess points as MESH0 and to that with grid points as MESH1in the figure. It is found that the linear growth rates agree with each other reasonably well although the coarser mesh slightly overestimates the growth time. Figure 16 (b) presents the power spectra that are time-averaged over the nonlinar phase. The random perturbation is imposed in this case. It is again clear that the results for MESH0 is in good agreement with those for MESH1.

It should be mentioned that for MESH0 it takes 32 parallel processors about 1.5 days to compute the evolution up to ms, while MESH1 needs 22 days even for 128 parallel processors. This is partly because of the difference in the Courant numbers, which are set to be 0.5 for MESH0 but to be 0.1 for MESH1 to better use the tensor-type artificial viscosity. Although this severe limit of CPU time does not allow us to do more thorough convergence tests, we think, based on the results of the tests shown above, that our results given in this paper are credible.

## References

• Blondin et al. (2003) Blondin, J. M., Mezzacappa, A., & DeMarino, C. 2003, ApJ, 584, 971
• Blondin & Mezzacappa (2006) Blondin, J. M., & Mezzacappa, A. 2006, ApJ, 642, 401
• Blondin & Mezzacappa (2007) Blondin, J. M., & Mezzacappa, A. 2007, Nature, 445, 58
• Blondin & Shaw (2007) Blondin, J. M., & Shaw, S. 2007, ApJ, 656, 366
• Burrows et al. (1995) Burrows, A., Hayes, J., & Fryxell, B. A. 1995, ApJ, 450, 830
• Burrows et al. (2006) Burrows, A., Livne, E., Dessart, L., Ott, C. D., & Murphy, J. 2006, ApJ, 640, 878
• Foglizzo et al. (2006) Foglizzo, T., Scheck, L., & Janka, H.-T. 2006, ApJ, 652, 1436
• Hayes et al. (2006) Hayes, J. C., Norman, M. L., Fiedler, R. A., Bordner, J. O., Li, P. S., Clark, S. E., ud-Doula, A., & Mac Low, M.-M. 2006, ApJS, 165, 188 S. E., & Spruit, H. C. 2005, ApJ, 626, 350
• Herant et al. (1994) Herant, M., Benz, W., Hix, W. R., Fryer, C. L., & Colgate, S. A. 1994, ApJ, 435, 339
• Janka & Mueller (1996) Janka, H.-T., & Mueller, E. 1996, A&A, 306, 167
• Iwakami et al. (2007) Iwakami, W., Kotake, K., Ohnishi, N. and Yamada, S. 2007, in preparation
• Kifonidis et al. (2006) Kifonidis, K., Plewa, T., Scheck, L., Janka, H.-T., Müller, E. 2006, A&A, 453, 661
• Kotake et al. (2003) Kotake, K., Yamada, S., & Sato, K. 2003, ApJ, 595, 304
• Kotake et al. (2006) Kotake, K., Sato, K., & Takahashi, K. 2006, Reports of Progress in Physics, 69, 971
• Ohnishi et al. (2006) Ohnishi, N., Kotake, K., & Yamada, S. 2006, ApJ, 641, 1018
• Ohnishi et al. (2007) Ohnishi, N., Kotake, K., & Yamada, S. 2007, accepted ApJTohline, J. E., & Burrows, A. 2005, ApJ, 625, L119
• Quirk (1994) Quirk, J. 1994, Internal Journal for Numerical Methods in Fluids, 18, 555
• Scheck et al. (2004) Scheck, L., Plewa, T., Janka, H.-T., Kifonidis, K., Müller, E. 2004, Physical Review Letters, 92, 011103
• Shen et al. (1998) Shen, H., Toki, H., Oyamatsu, K., Sumiyoshi, K. 1998, Nuclear Physics, A637, 43, 109, 301
• Stone & Norman (1992) Stone, J. M. & Norman, M. L. 1992, ApJS, 80, 753
• Woosley & Janka (2005) Woosley, S., & Janka, T. 2005, Nature Physics, 1, 147
• Wang et al. (1996) Wang, L., Wheeler, J. C., Li, Z., & Clocchiatti, A. 1996, ApJ, 467, 435
• Wang et al. (2001) Wang, L., Howell, D. A., Höflich, P., & Wheeler, J. C. 2001, ApJ, 550, 1030
• Wang et al. (2002) Wang, L., et al. 2002, ApJ, 579,
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters