A Parametric Study of the Acoustic Mechanism for Core-Collapse Supernovae

TL;DR

This study examines the conditions under which acoustic waves can revive stalled shock waves in core-collapse supernovae, revealing that a combination of neutrino heating and acoustic power is crucial for successful explosions.

Contribution

It introduces a parametric framework with a critical surface to evaluate the acoustic mechanism's effectiveness in supernova shock revival, extending previous models.

Findings

Large acoustic powers are needed in 1D to revive shocks due to neutrino cooling effects.

In 2D, the required acoustic powers align with previous estimates by Burrows et al.

Shock revival depends on the combined effect of neutrino heating and acoustic wave power.

Abstract

We investigate the criterion for the acoustic mechanism to work successfully in core-collapse supernovae. The acoustic mechanism is an alternative to the neutrino-heating mechanism. It was proposed by Burrows et al., who claimed that acoustic waves emitted by $g$-mode oscillations in proto-neutron stars (PNS) energize a stalled shock wave and eventually induce an explosion. Previous works mainly studied to which extent the $g$-modes are excited in the PNS. In this paper, on the other hand, we investigate how strong the acoustic wave needs to be if it were to revive a stalled shock wave. By adding the acoustic power as a new axis, we draw a critical surface, an extension of the critical curve commonly employed in the context of neutrino heating. We perform both 1D and 2D parametrized simulations, in which we inject acoustic waves from the inner boundary. In order to quantify the power of acoustic waves, we use the extended Myers theory to take neutrino reactions into proper account. We find for the 1D simulations that rather large acoustic powers are required to relaunch the shock wave, since the additional heating provided by the secondary shocks developed from acoustic waves is partially canceled by the neutrino cooling that is also enhanced. In 2D, the required acoustic powers are consistent with those of Burrows et al. Our results seem to imply, however, that it is the sum of neutrino heating and acoustic powers that matters for shock revival.

Tables

Table 1. Table 1: Success/Failure Score-sheet

$\dot{M}=1.0M_{☉}{s}^{-1}$
$L_{\nu }=4.5\times {10}^{52}\mathrm{erg}{s}^{-1}$		$L_{\nu }=4.0\times {10}^{52}\mathrm{erg}{s}^{-1}$
$\delta =0.40$	successful	$\delta =0.24$	successful
$\delta =0.39$	successful	$\delta =0.23$	successful
$\delta =0.38$	failed	$\delta =0.22$	failed
$\delta =0.37$	successful	$\delta =0.21$	successful
$\delta =0.36$	failed	$\delta =0.20$	failed
$\delta =0.35$	successful	$\delta =0.19$	failed
$\delta =0.34$	failed
$\delta =0.33$	failed

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{v})=0,$$

$$\frac{\partial \rho \mathbf{v}}{\partial t}+\nabla \cdot (\rho \mathbf{v}\mathbf{v}+P\mathbf{I})=-\rho \nabla \Phi ,$$

$$\frac{\partial \rho (e+\text{text@frac}12{\mathbf{v}}^2)}{\partial t}$$

$$\frac{\partial \rho Y_{\mathrm{e}}}{\partial t}+\nabla \cdot (\rho \mathbf{v}{Y}_{\mathrm{e}})=\rho \Gamma ,$$

$$\Phi =-\frac{G{M}_{\mathrm{PNS}}}{r},$$

$$\frac{\partial E_{\mathrm{dis}}}{\partial t}+\nabla \cdot {\mathbf{F}}_{\mathrm{dis}}=-D_{\mathrm{dis}},$$

$$E_{\mathrm{dis}}$$

$${\mathbf{F}}_{\mathrm{dis}}$$

$$D_{\mathrm{dis}}$$

$$+(\mathbf{m}$$

$$-(T$$

$$(m_r-m_{0r})(H-H_0)\sim 10^{36\text{–}38}\mathrm{erg}\mathrm{c}{\mathrm{m}}^{-2}{\mathrm{s}}^{-1}.$$

$$(m_r-m_{0r})T_0(s-s_0)$$

$$m_{0r}(T-T_0)(s-s_0)$$

$$(m_r-m_{0r})\frac{{\mu }_0}{m_{\mathrm{u}}}(Y_{\mathrm{e}}-Y_{e0})$$

$$m_{0r}\frac{\mu -{\mu }_0}{m_{\mathrm{u}}}(Y_{\mathrm{e}}-Y_{e0})$$

$$F_{\mathrm{dis},r}\sim 10^{36-38}\mathrm{erg}\mathrm{c}{\mathrm{m}}^{-2}{\mathrm{s}}^{-1}.$$

$${\dot{E}}_{\mathrm{aco}}\sim 10^{50\text{–}52}\mathrm{erg}{\mathrm{s}}^{-1}.$$

$$\frac{{\tau }_{\mathrm{adv}}}{{\tau }_{\mathrm{heat}}},$$

$$E_{\mathrm{turb}}=\frac{1}{2}{\int }_{\mathrm{gain}}\mathrm{d}V\rho \left(v_{\theta }^2+(v_r-⟨v_r⟩{)}^2\right),$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \mathbf{m}=0$$

$$\frac{\partial \mathbf{v}}{\partial t}+\mathbf{\zeta }+\nabla H-T\nabla s-\frac{\mu }{m_{\mathrm{u}}}\nabla Y_{\mathrm{e}}=-\nabla \Phi +\frac{1}{\rho }\mathbf{M}$$

$$\frac{\partial \rho s}{\partial t}+\nabla \cdot (\mathbf{m}s)=\frac{Q-\frac{\mu }{m_{\mathrm{u}}}\gamma -\mathbf{v}\cdot \mathbf{M}}{T}$$

$$\frac{\partial \rho Y_{\mathrm{e}}}{\partial t}+\nabla \cdot (\mathbf{m}{Y}_{\mathrm{e}})=\rho \Gamma$$

$$\mathrm{d}e$$

$$\mathrm{d}h$$

$$\frac{\partial }{\partial t}(\rho H-P)+\nabla \cdot (H\mathbf{m})+\mathbf{m}\cdot \nabla \Phi -Q=0.$$

$$C_i=0,{\mathbf{L}}_i={\mathbf{\lambda }}_i,S_i={\sigma }_i,\text{and}{G}_i={\gamma }_i,$$

$$\nabla \cdot ({\mathbf{m}}_0{H}_0)+{\mathbf{m}}_0\cdot \nabla \Phi -Q_0$$

$$\frac{\partial }{\partial t}(\rho H-P{)}_1+\nabla \cdot ({\mathbf{m}}_0{H}_1+{\mathbf{m}}_1{H}_0)+{\mathbf{m}}_1\cdot \nabla \Phi -Q_1$$

$$\frac{\partial }{\partial t}(\rho H-P{)}_2+\nabla \cdot ({\mathbf{m}}_0{H}_2+{\mathbf{m}}_1{H}_1+{\mathbf{m}}_2{H}_0)+{\mathbf{m}}_2\cdot \nabla \Phi -Q_2=0\text{(second-order)}.$$

$${\left(H-Ts-\frac{\mu }{m_{\mathrm{u}}}{Y}_{\mathrm{e}}\right)}_0{C}_0$$

$${\left(H-Ts-\frac{\mu }{m_{\mathrm{u}}}{Y}_{\mathrm{e}}\right)}_0{C}_1$$

$$+{\left(H-Ts-\frac{\mu }{m_{\mathrm{u}}}{Y}_{\mathrm{e}}\right)}_1{C}_0$$

$${\left(H-Ts-\frac{\mu }{m_{\mathrm{u}}}{Y}_{\mathrm{e}}\right)}_0{C}_2$$

Full text

methods: numerical – shock waves – supernovae: general

1 Introduction

Core-collapse supernovae (CCSNe) are deaths of massive stars. They are initiated by the gravitational collapse of an iron core at the end of the evolution and conclude with the explosion of outer envelopes with an energy of $\sim 10^{51}\,{\rm erg}$, which is accompanied with the formation of compact objects such as neutron stars and black holes. The explosion is assumed to be produced by the passage of the shock wave generated by the core bounce, which occurs when the central density exceeds the nuclear saturation density. The long-standing problem is that the shock wave stalls inside the core, consuming its energy to dissociate nuclei, and it cannot be clearly determined how this once-stalled shock is revived and eventually produces an explosion.

The current leading hypothesis is the neutrino-heating mechanism, in which the energy of $\sim 10^{53}\,{\rm erg}$ that is stored in the central portion of the core, or the proto-neutron star (PNS), is emitted in the form of neutrinos, and a fraction of the neutrinos is reabsorbed to heat and revive the shock.

According to recent studies (see Janka et al., 2016, for a review), multidimensional fluid instabilities such as neutrino-driven convection and standing accretion shock instability (SASI) are crucially important for the scenario: these instabilities induce turbulence and push the shock outward with thier turbulent pressure (Murphy et al., 2013; Couch & Ott, 2015); lateral motions in the turbulence increase the dwell time of fluid elements in the so-called gain region, where heating dominates cooling (Takiwaki et al., 2012); in multidimensions, both upward motions of hot and buoyant matter and down-flows of cold material are realized simultaneously, which enhances the efficiency of heating.

All these effects combined have recently produced successful shock revivals in many 2D and some 3D simulations, which failed earlier in 1D under spherical symmetry (Rampp & Janka, 2000; Liebendörfer et al., 2001; Thompson et al., 2003; Sumiyoshi et al., 2005). Some problems remain in these seemingly successful models, however. For one, results from different groups appear to be at odds with one another (Takiwaki et al., 2012; Bruenn et al., 2013; Dolence et al., 2015; Takiwaki et al., 2014; Tamborra et al., 2014; Müller, 2015). For another, the explsion energy obtained in these simulations is commonly lower than the canonical value of $\sim 10^{51}\,{\rm erg}$ (see, however, Bruenn et al., 2016). These problems may imply that some important physical process(es) is (are) still missing in the neutrino-heating mechanism.

The acoustic mechanism was proposed by Burrows et al. (2006, 2007a, 2007b) as an alternative to the neutrino-heating mechanim. In the context explained above, the “missing physics” is the oscillation of PNS. In this mechanism, turbulent accretion flows caused by SASI and/or convection beat the PNS anisotropically and repeatedly. These impulsive forces excite $g$-mode oscillations of PNS, which then emit acoustic waves into the accreting matter. The waves are steepened to form secondary shock waves as they propagate outward, and they deposit energy to shocked matter. The PNS in this way transduces the gravitational energy of the accretion flow into the energy of acoustic waves. This conversion of energy was claimed to last until an explosion is instigated, since the matter accretion, the energy source for acoustic waves, also continues until the explosion occurs. Another merit of the acoustic waves that the Burrows’ group claimed is that the waves do not go beyond the shock wave and hence deposit all their energies inside it.

One of the key issues in the acoustic mechanism is to what extent the oscillations of PNS are excited. In their papers, Burrows et al. (2006, 2007a, 2007b) found numerically that the energies of the excited $g$-modes are $10^{50}\text{--}10^{51}\,{\rm erg}$ with the $\ell=1$ mode being the most pronounced with a period of $\sim 3\,{\rm ms}$, where $\ell$ is the degree of the spherical harmonics. They also observed that these oscillations of PNS emit acoustic waves with powers of $\sim 10^{51}\,{\rm erg\,s^{-1}}$. Unfortunately, these results have not been reproduced in the simulations by other groups (Marek & Janka, 2009). Although in simulations of Marek & Janka (2009) the small central region is treated in 1D, they claimed and checked that the spherically symmetrized region was so small that their simulations could reproduce the core $g$-mode oscillation if it really occurred. It is true that such differences in numerical methods between Burrows’ group and others may not be ignored, but it is unclear why other groups fail to reproduce strong $g$-mode oscillations. Long-term computations required to obtain acoustic waves also hamper systematic studies by realistic simulations.

Yoshida et al. (2007) took a different, more phenomenological approach based on linear perturbation theory. They calculated forced oscillations of PNS by pressure fluctuations that are possibly imparted to the PNS surface by SASI. Employing the results of the numerical simulation of SASI by Ohnishi et al. (2006), they estimated that the energies of the excited $g$-modes are $\lesssim 10^{50}\,{\rm erg}$ and argued that these forced oscillations would inject energy at a rate of $\sim 10^{51}\,{\rm erg\,s^{-1}}$ as acoustic waves, which is comparable to what Burrows et al. (2006) obtained.

Then came a serious challenge from Weinberg & Quataert (2008), who analyzed the nonlinear three-mode couplings among $g$-modes, which transfer energy to two “daughter” modes with lower frequencies. The wave energies are assumed to be eventually dissipated and emitted as neutrinos. According to their calculations, the excitation of the most pronounced mode with $\ell=1$ is saturated at $\sim 10^{47-48}\,{\rm erg}$ for the steady energy feed at $\dot{E}=10^{50-51}\,{\rm erg\,s^{-1}}$ from the matter accretion, which, if we assume the acoustic damping rate of $10\,{\rm Hz}$, gives the acoustic power of $\sim 10^{48}\text{--}10^{49}\,{\rm erg\,s^{-1}}$. These values are $\sim 10\text{--}100$ times lower than the values given in Burrows et al. (2006). We note that the numerical simulations by the latter authors have most likely failed to capture these mode couplings, since one of the daughter modes of relevance for resonant couplings has such short wavelengths that the numerical resolutions were not sufficient. It should be mentioned that Yoshida et al. (2007) did not take the mode coupling into account, either. This is certainly a serious issue, but we recall that Weinberg & Quataert (2008) also made some assumptions in their analysis. For instance, although they assumed a steady energy injection from turbulent accretion flows to $g$-modes, this may not be a good approximation, since the forces excerted on the PNS suraface are more like a collection of impulsive hits, as envisaged by Yoshida et al. (2007). It also remains to be confirmed if the mode couplings neglected in their paper are indeed minor. Investigations by other independent groups are certainly desirable.

It is true that the best way is in principle either to improve the numerical resolution sufficiently or to conduct a fully nonlinear analysis of the mode couplings, but this is almost impossible for the moment. We hence take yet another way in this paper. Reversing the argument, we ask here what acoustic power is needed to revive the stalled shock wave. This is in accord with the spirit of the critical curve theory for the neutrino-heating mechanism (Burrows & Goshy, 1993), in which the critical neutrino luminosity for shock revival is considered as a function of the mass accretion rate. In this paper we add a new dimension in this theory and discuss the critical “surface.” Although the rotational velocity was introduced by Iwakami et al. (2014) to consider such a critical surface, the new dimension introduced in this paper is the intensity of acoustic waves. Then the critical surface will allow us to assess how intense acoustic waves need to be to obtain shock revival, which will in turn help us judge how promising the acoustic mechanism is in realistic settings. We note that in this framework of the critical surface we are in fact treating the acoustic and neutrino powers on an equal footing, which may be in contrast to the original claim by Burrows et al. that the acoustic waves are the dominant agent for shock revival. In this sence, our model may be referred to as a “hybrid” model.

In order to obtain the critical surface, we perform both 1D and 2D simulations under spherical and axial symmetries, respectively. Since the $g$-modes are non-spherically symmetric and hence the acoustic mechanism works only in multidimensions, the 1D simulations may not be realistic. Spherical acoustic waves are still conceivable, however, and they will be easier to analyze than non-spherical conterparts and hence are useful to capture the essential features in the propagation and energy deposit of acoustic waves.

This paper is organized as follows. In section 2 we describe the numerical methods. The results of simulations are shown in the subsequent sections. In section 3 we first summarize 1D models. Then 2D results are presented in section 4. In section 5 we give summary and some discussions.

2 Methods

In our models we prepare spherically symmetric steady accretion flows with a constant mass accretion rate $\dot{M}$ and neutrino luminosity $L_{\nu}$ and inject acoustic waves continuously from the inner boundary of the computation domain, which is located close to the PNS surface. For various combinations of $\dot{M}$ and $L_{\nu}$, we search the critical amplitudes of acoustic wave, which are the minimum amplitudes required for shock revival. We regard it as a successful shock revival if the shock reaches the radius of $500\,{\rm km}$ within $500\,{\rm ms}$ from the onset of the acoustic wave injection. We then draw the critical surface in the space spanned by $\dot{M}$, $L_{\nu}$, and the amplitude of the acoustic wave.

Basic equations in our calculations are inviscid hydrodynamics equations with neutrino emissions and absorptions as an energy sink/source and an equation for the electron fraction:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\rho$, $\boldsymbol{v}$, $P$, $e$, $Q$, $Y_{\rm e}$, $\Gamma$, $\Phi$, $G$, $M_{\rm PNS}$, and $r$ are the density, velocity, pressure, specific internal energy, net heating rate via neutrino emissions and absorptions, electron fraction, rate of change in $Y_{\rm e}$ by the neutrino reactions, gravitational potential by PNS, gravitational constant, mass of PNS, and distance from the center, respectively; $\boldsymbol{I}$ and $\boldsymbol{v}\boldsymbol{v}$ are the unit and dyadic tensors, respectively; the self-gravity of the accreting matter is ignored, and only the gravitational attraction by PNS is considered. In the following calculations, the PNS mass is fixed to $M_{\rm PNS}=1.4\,M_{\odot}$. We employ the so-called STOS equation of state (EOS) (Shen et al., 1998), which is based on the relativistic mean field theory and the Thomas-Fermi approximation. The light-bulb method in Ohnishi et al. (2006)111We drop $\pi$ from the denominator of equation (18) in their paper. is adopted to calculate $Q$ and $\Gamma$ in equations (3) and (4). We consider only the absorption and emission of neutrinos by free nucleons, $\nu_{\rm e}+{\rm n}\leftrightarrow{\rm e}^{-}+{\rm p}$ and $\bar{\nu}_{\rm e}+{\rm p}\leftrightarrow{\rm e}^{+}+{\rm n}$. The neutrino temperatures are set to $T_{\nu_{\rm e}}=4\,{\rm MeV}$ for $\nu_{\rm e}$ and $T_{\bar{\nu}_{\rm e}}=5\,{\rm MeV}$ for $\bar{\nu}_{\rm e}$.

We run 1D and 2D simulations on the spherical coordinates under spherical and axial symmetries, respectively. The inner boundary of the computational domain $r_{0}$ is fixed to the neutrinosphere $r_{\nu}$, which is defined in this paper to be the radius at the density of $10^{11}\,{\rm g\,cm^{-3}}$ in the initial condition, which is explained below. The radial mesh width $\Delta r_{i}$ at the $i$th grid point is set to $1\%$ of the radius $r_{i}$: $\Delta r_{i}=0.01r_{i}$. The number of radial grid points is $256$, but if the outermost radius $r_{256}$ is smaller than $500\,{\rm km}$, we increase the number to $320$ so that $r_{320}$ exceeds $500\,{\rm km}$. In 2D models, the entire meridian section is covered with the same radial grid points and $128$ $\theta$-grid points, the latter of which are deployed according to the Gaussian quadrature points and weights as in Sumiyoshi & Yamada (2012).

The hydrodynamical code employed in this paper is the same as that in Nagakura et al. (2014), except that only the point-mass gravity instead of the self-gravity is considered: the Harten-Lax-van Leer scheme (Harten et al., 1983) with the piecewise parabolic interpolation (Colella & Woodward, 1984) is used to evaluate the numerical flux; the time evolution is handled by the explicit, total-variation diminishing, third-order Runge-Kutta method.

The initial conditions are time-independent solutions of equations (1)–(5) for given combinations of constant mass accretion rate and neutrino luminosity. The numerical method to obtain these solutions is essentially the same as that in Yamasaki & Yamada (2006), except for two aspects: one is again the use of the point-mass gravity instead of the self-gravity, and the other is the definition of the neutrinosphere $r_{\nu}$, which in this paper is defined to be the radius at which the density is $10^{11}\,{\rm g\,cm^{-3}}$, wheres in Yamasaki & Yamada (2006) it was the radius, where the optical depth is $2/3$. Regardless, the neutrino luminosity is assumed to satisfy the following relation at the neutrinosphere: $L_{\nu}=\frac{7}{16}\sigma T_{\nu_{\rm e}}^{4}4\pi r_{\nu}^{2}$. In the upstream of the standing shock wave, on the other hand, we assume that matter is accreting spherically symmetric with the entropy per nucleon $s$, electron fraction $Y_{\rm e}$, and radial velocity $v_{r}$ being $s=3$ in units of the Boltzmann constant $k_{\rm B}$, $Y_{\rm e}=0.5$, and $v_{r}=\sqrt{2GM/r}$, respectively.

We inject acoustic waves from the inner boundary. We impose time-dependent boundary conditions there to generate outgoing sound waves. We assume a sinusoidal oscillation in the density, $\rho=\rho_{0}+\rho_{1}\sin(\omega t-kr)$, where $\rho_{0}$ is the density of the steady state at the inner boundary, whereas $\rho_{1}$, $\omega$, and $k$ are the amplitude, frequency, and wavenumber of the oscillation, respectively. We assume $\rho_{1}\propto\mathcal{P}_{0}(\mu)$, where $\mathcal{P}_{\ell}$ is the Legendre polynomial of $\ell$-th order and $\mu=\cos\theta$, in 1D simulations. Although $g$-mode oscillations, the main source of the acoustic waves, are intrinsically non-spherical, we consider 1D spherically symmetric acoustic waves with $\ell=0$ as well in this paper, since they elucidate the essential feature of the acoustic energy transport. In 2D simulations, on the other hand, we set $\rho_{1}\propto\mathcal{P}_{1}(\mu)$, although we can consider any non-zero $\ell$ in principle, since $\ell=1$ modes were the most prominent in Burrows et al. (2006). Throughout this paper, the normalized dimensionless amplitude of acoustic wave $\delta$ is defined as $\rho_{1}=\rho_{0}\mathcal{P}_{\ell}\delta$ for $\ell=0$ (1D) and $\ell=1$ (2D). The entropy per nucleon $s$ and electron fraction $Y_{\rm e}$ at the inner boundary are fixed to the steady-state values, reflecting the adiabatic character of acoustic waves. Other thermodynamic quantities such as temperature and pressure are determined by the EOS. The inner boundary condition for velocity is determined so that it should be consistent with the outgoing sound waves: $v_{r}=v_{0}+a\mathcal{P}_{\ell}\delta\sin(\omega t-kr)$, where $v_{0}$, $a=\sqrt{(\partial P/\partial\rho)_{s,Y_{\rm e}}}$, and $k=\omega/(v_{0}+a)$ are the velocity, sound speed, and wave number at the inner boundary in the steady state, respectively. The frequancy is set to $\omega=2\pi/(3\,{\rm ms})$, i.e., the oscillation period is $3\,{\rm ms}$, the value of the dominant $g$-mode oscillation in Burrows et al. (2006). Incidentally, the values of all quantities at the outer boundary are fixed to the values in the steady state.

For each combination of $\dot{M}$, $L_{\nu}$, and $\delta$, we run a simulation for $500\,{\rm ms}$. If the mean shock radius exceeds $500\,{\rm km}$ within this period, we interpret it as shock revival and consider that this model produces a successful explosion. This value of the radius is the same as the value adopted in Iwakami et al. (2014) and slightly more conservative compared with the value of $400\,{\rm km}$ employed in Nordhaus et al. (2010) and Hanke et al. (2012). The shock radius is defined in this paper as the radius at which the entropy per nucleon is $s=6\,k_{\rm B}$. This is roughly twice the entropy in the unshocked accretion flow, as explained above. We vary $\delta$ at regular intervals of $0.005$ (1D) or $0.01$ (2D), and search for the threshold of $\delta$ for shock revival for each combination of $\dot{M}$ and $L_{\nu}$. Connecting the points obtained this way, we draw the critical surface, a 2D analogue of the critical curve.

Before considering the results, we mention the dependence on the numerical resolution and the initial phase of the acoustic wave. In order to assess these effects, we conducted additional simulations. We first reduced the radial grid width by half in the 1D simulations with the mass accretion rate $\dot{M}=0.6\,M_{\odot}\,{\rm s^{-1}}$ and neutrino luminosity $L_{\nu}=4.0\times 10^{52}\,{\rm erg\,s^{-1}}$ and found that the critical surface is shifted only by $0.01\text{--}0.015$ in the positive direction of the $\delta$ axis. As for the choice of the initial phase of the oscillation, we ran four additional computations for the same 1D model with the original resolution, but with the phase being changed by $\pi/2$, and we confirmed that the critical surface is shifted only by $0.005$ in $\delta$. We thus concluded that the critical surface in 1D is determined fairly unambiguously. We also checked the numerical resolution in 2D simulations and confirmed that it is sufficient. In fact, doubling the number of grid points in the radial (angular) direction in the simulations with $\dot{M}=1.0\,M_{\odot}\,{\rm s^{-1}}$ and $L_{\nu}=4.5\times 10^{52}\,{\rm erg\,s^{-1}}$ lowered the critical surface only by $\sim 0.02$ ($\sim 0.01$) in $\delta$. Moreover, the change in initial phase by $\pi/2$ in the same model shifts the critical surface only by $\sim 0.01$ in $\delta$. This shows that the critical surface is probably also well determined in 2D.

3 Results in 1D

3.1 Critical Surface

The critical surface we obtained from the 1D simulations is shown in the upper panel of Figure 1. By definition, models with the parameters on or above this surface result in shock revival, whereas those beneath it fail. In the lower panel, three lines on the critical surface, each of which has an identical mass accretion rate, are projected onto the $\delta\text{--}L_{\nu}$ plain. Before discussing these results in detail, we first consider a representative model and examine its evolution.

In Figure 2 we show the temporal evolution of the radial velocity in the model with $\dot{M}=0.6\,M_{\odot}\,{\rm s^{-1}}$, $L_{\nu}=4.0\times 10^{52}\,{\rm erg\,s^{-1}}$, and $\delta=0.105$ as an example. This model is located on the critical surface. Discontinuous changes in color in the figure correspond to shock waves. We can clealy see radial oscillations of the primary-shock wave with a period of $\sim 70\,{\rm ms}$ in the lower panel. This is an oscillatory instability found by Fernández (2012), since the timescale of the advection from the shock to the point of maximum cooling is $\sim 24\,{\rm ms}$, indeed between $1/4$–$1/2$ of the period as expected in his paper. The upper panel, on the other hand, is a zoom-in to the designated area in the lower panel. It is evident in this panel that many shock waves are propagating outward on top of the oscillatory mode and periodically hit the primary-shock wave. These secondary shock waves are generated as a consequence of the steepening of acoustic waves injected from the inner boundary. This is understood from the fact that the interval between the consecutive shock waves is exactly equal to the period of the acoustic wave. Although these shocks are rather weak (the typical Mach number is $\sim 1.3$) in this model, each collision of a secondary shock with the primary shock causes a tremor in the latter, which is clearly visible in the upper panel.

The oscillatory motion of the primary-shock wave repeats itself with growing amplitudes. These overstable oscillations eventually lead to shock revival in this model, as demonstrated in the lower panel of Figure 2. Figure 3 shows other models with various acoustic amplitudes $\delta$. The inset is a zoom-in of the rectangular area in the main panel and displays tremblings caused by the secondary shock waves in these models. We can see from the figure that as the acoustic amplitude $\delta$ decreases, the time to shock revival increases, and eventually, no shock revival occurs, indicating that there is a critical amplitude for shock revival. This is common to other models with different combinations of $\dot{M}$ and $L_{\nu}$. The collection of these critical amplitudes gives the critical surface shown in Figure 1.

Although the acoustic amplitude is employed as the independent variable in the simulations, the acoustic power $\dot{E}_{\rm aco}$, which is the rate at which matter is heated by the acoustic wave and should be compared with the neutrino-heating rate, is more convenient for physical interpretations and comparisons with previous works (Burrows et al., 2006; Yoshida et al., 2007; Weinberg & Quataert, 2008). Since acoustic waves do not cross the shock and hence cannot escape from the postshock region, the acoustic power $\dot{E}_{\rm aco}$ should be equal to the acoustic luminosity $L_{\rm aco}$ at the inner boundary, which is the surface integral of the acoustic energy flux. It is well known that the energy density and flux of an acoustic wave are proportional to the amplitude squared as long as its amplitude is small enough to be in the linear regime. In our models, however, the amplitudes are not small in general and the acoustic waves may not be regarded as linear waves. Moreover, the fact that the acoustic waves are propagating on top of matter that is not at rest, but flows non-uniformly, complicates the evaluation of acoustic power even more. In order to handle these problems, we extend the Myers theory (Myers, 1986, 1991).

Myers derived the corollary of the energy conservation for disturbances in homentropic flows, i.e., acoustic waves (Myers, 1986), as well as in flows with inhomogeneous entropies (Myers, 1991). Here we extend the discussion of Myers (1991) in order to take the effects of neutrino reactions into account. As derived in appendix A, the resulting equation is expressed as follows:

[TABLE]

where each term is given as

[TABLE]

[TABLE]

[TABLE]

In the above equations, $m_{\rm u}$, $H=e+P/\rho+\text@frac{1}{2}{\boldsymbol{v}}^{2}$, ${\boldsymbol{m}}=\rho{\boldsymbol{v}}$, $T$, and $s$ are the atomic mass unit, specific stagnation enthalpy (or the Bernoulli function), mass flux, temperature, and specific entropy, respectively; $\mu$ is the chemical potential of electron-type neutrinos defined with the chemical potentials of electrons, protons, and neutrons, $\mu_{\rm e,p,n}$, as $\mu=\mu_{\rm e}+\mu_{\rm p}-\mu_{\rm n}$; $\boldsymbol{\zeta}$ is defined as $\boldsymbol{\zeta}=\boldsymbol{\omega}\boldsymbol{\times}\boldsymbol{v}$, where $\boldsymbol{\omega}=\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{v}$ is the vorticity. The quantities with and without subscript [math] stand for the unperturbed and perturbed variables, respectively. Although it is difficult to give an unambiguous interpretation to each term of Eq. (6), we regard $E_{\rm dis}$, $\boldsymbol{F}_{\rm dis}$, and $D_{\rm dis}$ as the density, flux, and dissipation, respectively, of the energy of acoustic waves with not necessarily small amplitudes. This may be justified by the facts that they obey the equation of a conservative form and that they are reduced to the well-known counterparts for linear waves if the amplitude is small and all neutrino contributions are turned off; they are hence natural extensions. See appendix A for more discussions.

One may think that the acoustic luminosity should be evaluated at the inner boundary, where the acoustic waves are generated artificially, but this may not be true, since the injected waves will be partially reflected back immediately after they enter the computational domain. This fact can be understood from Figure 4, which shows the acoustic luminosities defined at each radius as the surface integral of the radial component of Myers’ flux, $L_{\rm aco}(r)=4\pi r^{2}F_{{\rm dis},r}(r)$, for models sitting on the critical surface. We note that the luminosity, $L_{\rm aco}(r)$, at the radius $r$ is obtained by taking its average over $3\,{\rm ms}$, the period of the acoustic wave, from the instant when the acoustic wave just reaches the radius. One can see small transients on the first two grid points in all cases. We hence decided to use the acoustic luminosity obtained at the third grid point from the inner boundary, where the initial adjustment appears to have been already over, to estimate the truly injected acoustic power $\dot{E}_{\rm aco}$. Fernández (2012) gave a similar argument that the third grid point is the innermost point that is not significantly affected by the inner boundary. We note that we can recognize a common trend in this figure that the acoustic luminosities decrease with radius. This is particularly significant for models with small $L_{\nu}$’s. The reduction of the acoustic luminosity may be attributed to the dissipation term $D_{\rm dis}$ in equation (6).

It will be useful to give an order-of-magnitude estimate of the acoustic power here. Only the radial components are considered for vectors. The typical values of the temperature, specific entropy, electron fraction, chemical potential, specific stagnation enthalpy, and mass flux are $k_{\rm B}T_{0}\sim{\rm MeV}$, $s_{0}\sim 10\,k_{\rm B}$ per nucleon, $Y_{\rm e0}\sim 0.1$, $\mu_{0}\sim{\rm MeV}$, $H_{0}\sim 10^{19}\,{\rm erg\,g^{-1}}$, and $m_{0r}\sim-10^{18\text{--}19}\,{\rm g\,cm^{-2}}$, respectively, at $r\simeq 50\,{\rm km}$. As for the disturbances, the amplitudes are typically $\sim 10\%$ of the unperturbed counterparts, except for the mass flux, for which $m_{r}-m_{0r}$ is $\sim 1\text{--}10$ times larger than $m_{0r}$ and is positive. Then,

[TABLE]

Similarly,

[TABLE]

and

[TABLE]

Combining all these contributions, one obtains

[TABLE]

Recalling that $4\pi r^{2}\sim 10^{14}\,{\rm cm^{2}}$, one can estimate the acoustic power as

[TABLE]

Employing the acoustic power obtained in this way, we draw the critical surface in the parameter space spanned by $\dot{M}$, $L_{\nu}$, and $\dot{E}_{\rm aco}$, which is presented in the upper panel of Figure 5. In the lower panel, on the other hand, we also show three lines on the critical surface, each of which connects the results for the identical mass accretion rate, and which are are projected onto the $L_{\nu}\text{--}\dot{E}_{\rm aco}$ plane. It is apparent that the acoustic power required for shock revival increases as the neutrino luminosity $L_{\nu}$ decreases. This is a clear indication that the acoustic power indeed contributes to shock revival.

3.2 Energetics

In order to bolster this picture, we investigate the energetics in more detail. In Figure 6 we compare the net neutrino-heating rate, the acoustic power, and the sum of them (or the total heating rate) for some models on the critical surface. Here the net neutrino-heating rate is the volume integral of $Q$ in equation (3) over the gain region, which is the region where the heating by neutrino absorptions dominates the cooling by neutrino emissions. As expected intuitively, more acoustic power is needed as the neutrino luminosity decreases. It is found that the decrease in net neutrino-heating rate is almost compensated for by the increase in acoustic power, and the total heating rate does not change significantly for models in which the neutrino-heating rate dominates the acoustic power. One may be tempted to think that the explosion occurs if the total heating rate exceeds a certain threshold determined by the mass accretion rate, but this is not the case. In fact, for models in which the acoustic power is greater than the neutrino-heating rate, the total heating rate required for shock revival is no longer constant, but increases rather quickly as the neutrino luminosity decreases. This may imply that such large-power acoustic waves are inefficient in depositing their energy.

The neutrino cooling may be responsible for the lower efficiency of the acoustic heating at large amplitudes. Once formed, the secondary shock waves raise the matter temperature, and as a consequence, enhance the neutrino cooling, since it is roughly proportional to $T^{6}$. As explained in appendix A, the term $-(T-T_{0})(Q/T-Q_{0}/T_{0})$ in $D_{\rm dis}$ (see equation (LABEL:eq:daco)) describes the energy loss of disturbances owing to the neutrino cooling. This effect will be more efficient for higher-power acoustic waves, since they will produce stronger secondary shock waves and lead to higher temperatures.

In order to see this effect more quantitatively, we show the profiles of velocity, temperature, and entropy at different times in Figures 7, 8, and 9. In these figures the cooling and heating regions are colored in blue and red, respectively. They are divided by the so-called gain radius. The cooling region sits initially closer to the inner boundary, where the temperature is higher, as shown in panels (a) of Figures 7 and 8, which correspond to the time just after the onset of simulations. As the secondary shock propagates outward, the cooling layer is extended to the secondary shock at first, as seen in panel (b) of Figure 7. This is because the shock wave raises the temperature of the traversed matter, and as a result, enchances the neutrino cooling, thus reducing the efficiency of the acoustic heating. Although Burrows et al. claimed that the acoustic mechanism is more efficient than the neutrino-heating mechanism as all the energy is eventually consumed inside the primary-shock wave, this does not necessarily mean that they are used entirely for shock revival.

While cooling layers extended to just behind the secondary shock, heating layers surrounded by cooling layers sometimes appear, as shown in panels (b) of Figures 7 and 8 and in Figure 9. These inner heating layers coincide with troughs in the temperature. The rarefaction of flows that follows the compression by the secondary shock leads to lower temperatures and hence suppressions of the cooling, eventually producing the heating layer there. The effect of such a inner heating layer is small, however, and the secondary shock waves enhance the cooling as a whole. We also note that the rarefaction of flows is rather sensitive to the handling of the inner boundary condition and may possibly be an artifact of such treatments.

On the other hand, not all shock heating is spent on the neutrino cooling, either. Figure 8 shows in addition to the initial entropy distribution (panel (a)) a strong entropy production at the instance of the collision of the secondary shock with the primary shock (panel (b)). We note that this time is slightly later than the time in panel (b) of Figure 7. Since the primary shock is normally located far from the gain radius, the energy deposition associated with the collision does not lead to high enough temperatures for cooling to surpass heating, and the provided energy is not spent on neutrino emissions. If the primary shock is distant enough, the secondary shock ceases to convert the heating region into the cooling region well before it hits the primary shock, as demonstrated in Figure 9. This is again simply because the temperature does not become high enough, and this is the reason why most of the acoustic power is still available for shock revival, even though some of the deposited energy is spent on neutrino emissions.

In addition to the neutrino cooling we just discussed, the reflection of secondary shock waves may also be the cause of the inefficiency in the acoustic heating at high amplitudes. According to the theory of the Riemann problem, reverse shocks or rarefaction waves are formed in general and propagate inward when the secondary shock waves collide with the primary shock. In our models, only the rarefaction waves are formed. As a consequence of these reflections, not all the power is provided to the primary shock. The reflected waves will hit the PNS and may be recycled, however. If this is really the case, the reflection of the secondary shock waves may not reduce the efficiency of energy deposition so much as we see here, since such recycling may not be properly taken into account in the Dirichlet-type inner boundary employed in our simulations.

3.3 Diagnostics for Shock Revival

When does shock revival occur? In the context of the neutrino-heating mechanism, several diagnostics have been proposed so far to predict it (Thompson et al., 2005; Pejcha & Thompson, 2012; Murphy & Dolence, 2017), although it is known that none of them is perfect. The purpose of this section is not to seek something better, but to determine whether they are useful in the current context. We examine two often-used diagnostics here: the timescale ratio, and the antesonic factor.

One of the most frequently employed diagnostics is the ratio of advection timescale to heating timescale:

[TABLE]

where $\tau_{\rm adv}=\int_{r_{\rm gain}}^{r_{\rm shock}}{\rm d}r/|v_{r}|$ and $\tau_{\rm heat}=\int_{r_{\rm gain}}^{r_{\rm shock}}{\rm d}V\rho|\Phi|/Q$, with ${\rm d}V$ being the volume element. The radii $r_{\rm gain}$ and $r_{\rm shock}$ are the gain and shock radius, respectively. Thompson et al. (2005) was the first to claim that shock revival occurs when this ratio exceeds unity, which is intuitively understandable. Although this condition was originally meant for the neutrino-heating mechanism, it may be applicable to the hybrid of acoustic power and neutrino heating studied in this paper if one sees it as the acoustic-wave-assisted neutrino-heating mechanism. We show in Figure 10 the timescale ratio as a function of time for both successful and failed models. As can be seen, the difference between the successful and unsuccessful models is subtle: the ratio exceeds one sometimes even for the failed models, whereas some times of shock oscillations occur commonly in the successful models before shock revivals even after the ratio reaches unity. We hence conclude that this diagnostic is not capable of distinguishing the successful from the failed models.

One of the reason for this failure may be the obvious fact that only the neutrino-heating is considered in this diagnostic, although acoustic waves also contribute to the heating of matter in the present mechanism. Taking the acoustic heating into account in the discussions of heating timescales may not be so easy, however, since it is highly impulsive as the energy deposition occurs mainly when the secondary shock waves collide with the primary shock. Thus we do not pursue this issue further in this paper.

According to the antesonic condition (Pejcha & Thompson, 2012), an explosion should occur when the maximum value of antesonic factor, ${\rm max}(a^{2}/v_{\rm esc}^{2})$, in the downstream of the stalled shock wave exceeds a certain critical value, $\sim 0.2$, where $v_{\rm esc}$ is the escape velocity. Figure 11 shows this factor for some models just on and slightly below the critical surface. We plot the maximum value of $a^{2}/v_{\rm esc}^{2}$ attained by the time of shock revival in each model. We see from the figure that the maximum antesonic factor is larger in failed than in successful models more often than not, which is at odds with what the original theory posited. We note, however, that the antesonic condition is the condition for the steady-state solution not to exist, and as such, it may not be applicable to the current models in which the dynamical effects are also important. We hence conclude that the antesonic condition is not useful, either, in the mechanism considered here.

Since the diagnostics considered above were originally proposed in the context of the neutrino-heating mechanism, as mentioned repeatedly, it may not be surprising that they are not applicable to the acoustic-neutrino hybrid mechanism. Certainly, something better is necessary, but considering that they are known to be imperfect even in the neutrino heating, this is beyond the scope of this paper.

4 Results in 2D

It is true that the 1D models are convenient to understand the relevant physics, but we recall that the original acoustic mechanism works in multidimensional settings, since the $g$-mode oscillations are intrinsically non-spherical. In this section, we present the results of the 2D simulations and discuss how the dimensionality affects the critical surface.

Let us first look at the typical evolution of 2D acoustic explosion. In Figure 12 we show in color the entropy distributions in the meridian section at different times for the model with $\dot{M}=0.6\,M_{\odot}\,{\rm s^{-1}}$, $L_{\nu}=4.0\times 10^{52}\,{\rm erg\,s^{-1}}$, and $\delta=0.07$, which successfully leads to shock revival. It is observed from the figure that the initially spherical shock is deformed preferentially along the symmetry axis by large plumes that are produced by the dipolar acoustic waves injected from the inner boundary in the 2D model. The shock morphology changes in time, and eventually, shock revival occurs.

The time evolution of the primary shock is displayed for both successful and failed models in Figure 13. Both models show similar evolutions in the early phases: the shock radii expand for the first $\sim 100\,{\rm ms}$ and remain almost constant for the next $\sim 100\,{\rm ms}$. The late-phase evolutions are different, on the other hand: the shock stays at almost the same position until the end of the simulation in the failed model, while in the successful model it rapidly expands to reach $500\,{\rm km}$ by the time of $\sim 410\,{\rm ms}$. We note that the difference in $\delta$ is just $0.01$ between the successful and unsuccessful models.

Figure 14 shows the critical surface in the space spanned by $\dot{M}$, $L_{\nu}$, and $\delta$ for 2D models. We note that unlike in 1D, the critical surface is not defined in a clear-cut way in 2D. This is understood from Table 1, in which we present success or failure of shock revival for two series of models with different neutrino luminosities: $L_{\nu}=4.5\times 10^{52}\,{\rm erg\,s^{-1}}$ and $4.0\times 10^{52}\,{\rm erg\,s^{-1}}$. The mass accretion rate is fixed to $\dot{M}=1.0\,M_{\odot}{\rm s^{-1}}$. As $\delta$ increases in each series, a successful model appears at some point, which may be a possible critical point. For slightly higher values of $\delta$, however, we again find failure of shock revival. When we increase $\delta$ further, we eventually find only success. This is the ambiguity of the critical surface in 2D. The reason is that the shock revival takes place rather stochastically as a result of turbulence behind the primary shock, which is induced by convection and/or SASI. In such situations one may define the critical surface either as the surface below which all models fail to explode or as the surface above which all models explode. If we adopt the former, the critical surface passes through the points $(\dot{M},L_{\nu},\delta)=(1.0\,M_{\odot}{\rm s^{-1}},4.5\times 10^{52}\,{\rm erg\,s^{-1}},0.35)$ and $(1.0\,M_{\odot}{\rm s^{-1}},4.0\times 10^{52}\,{\rm erg\,s^{-1}},0.21)$, while it should pass through the points $(\dot{M},L_{\nu},\delta)=(1.0\,M_{\odot}{\rm s^{-1}},4.5\times 10^{52}\,{\rm erg\,s^{-1}},0.39)$ and $(1.0\,M_{\odot}{\rm s^{-1}},4.0\times 10^{52}\,{\rm erg\,s^{-1}},0.23)$ in the latter definition. In Figure 14 we adopt the former definition.

In the lower panel of Figure 14 we show in the $L_{\nu}\text{--}\delta$ plane some lines on the critical surface that have constant $\dot{M}$ for both 1D and 2D models. It is found that the critical amplitudes are not much different between 1D and 2D models, but are smaller in 2D than in 1D for large $L_{\nu}$ and vice versa for small $L_{\nu}$. This behavior may not be so important, however, since not the amplitude, but the acoustic power should be a more direct and hence a better measure for the shock revival.

Figure 15 shows the radial component of the Myers flux as a contour in the meridian section. The lateral flux is negligible compared to the radial flux and is not shown. The black thin semi-circle represents the initial shock radius in each panel. We find that almost everwhere inside the shock, the Myers flux is directed radially outward, and as expected, it has a dipolar angle-dependence, i.e., it is more intense close to the symmetry axis than near the equator. This suggests that the Myers flux gives an appropriate estimate of fluxes to acoustic waves with not necessarily small amplitudes. We also note that the negative acoustic fluxes in the vicinity of the inner boundary mean that acoustic waves are reflected inward there because of the steep density gradient in the background flow.

Since the Myers fluxes are positive in almost all directions, we can employ the acoustic luminosity $L_{\rm aco}$, the surface integral of the acoustic flux, to estimate the acoustic power. As in 1D, we also take the temporal average over the oscillation period, $3\,{\rm ms}$. The resulting acoustic luminosities are shown in Figure 16. We note that the acoustic luminosity at a radius $r$ is evaluated when the mean radius of the acoustic wave front exceeds $r$. Since it is almost constant in radius except on the first two radial grid points from the inner boundary, where some adjustments are taking place, we again adopt its value at the third grid point as the true injected acoustic power.

The critical surface for the acoustic power instead of the amplitude is presented in Figure 17. In contrast to the surface for the amplitude, the critical surface for the acoustic power in 2D is systematically lower than the 1D counterpart. This might to be though to be at odds with the previous findings that the critical surface for the acoustic amplitude is higher in 2D than in 1D at small $L_{\nu}$. The apparent contradiction is due to the different angular dependence of the acoustic waves between 1D and 2D. Since the angular integrations of the squared Legendre polynomials of our current interest are given as $\int_{-1}^{1}\mathcal{P}_{0}(\mu)^{2}{\rm d}\mu=2$ and $\int_{-1}^{1}\mathcal{P}_{1}(\mu)^{2}{\rm d}\mu=2/3$, the acoustic luminosities and hence powers as well are lower for the $\ell=1$ mode than for the $\ell=0$ mode if they have the same amplitude.

We study the energetics in more detail. As in the previous section, we compare in Figure 18 the sums of the neutrino-heating rate and acoustic power with the neutrino-heating rate alone for some models on the critical surface in both 1D and 2D. One can see again that shock revival in 2D requires less energy than in 1D. For the models marked with circles in the figure, the acoustic power needed for shock revival is much smaller than the neutrino-heating rate. In these models acoustic waves play a minor role as energy sources. Instead, they act as a driver of hydrodynamical instabilities, which enhance the neutrino heating. Since the shock revival occurs essentially by neutrino heating alone in these models, the effects of the turbulence may be estimated from the turbulent kinetic energy, dwell time, and neutrino-heating rate, which is demonstrated in Figure 19. The top panel of Figure 19 shows the turbulent kinetic energy $E_{\rm turb}$ defined as

[TABLE]

where $\langle v_{r}\rangle$ is the angle-averaged radial velocity, and the integral is performed over the gain region. It can be clearly seen that the injection of acoustic waves induces turbulent matter motions. The turbulence increases the dwell time in the gain region, and as a result, the gain mass, which is the mass in the gain region, also increases, as seen in the middle panel of Figure 19. This in turn raises the neutrino-heating rate integrated over the gain region, as is apparent in the lower panel of the same figure. Since such an enhancement is absent in the model without the injection of acoustic waves, we can conclude that acoustic waves are still playing an important role even with small amplitudes in enhancing the neutrino-heating rate via the fluid instability. Shock revival occurs essentially not via the acoustic mechanism, but via the neutrino-heating mechanism in these models.

For models with smaller $L_{\nu}$, the energy injection by acoustic waves plays a substantial role and it is the acoustic mechanism that gives rise to shock revival in these models. Here again the total powers required for shock revival are smaller in 2D than in 1D. The reason for the lower critical total powers in 2D is probably the enhancement of neutrino heating by the fluid instability again, although it is a minor player in this regime. In Figure 18 one also recognizes that the total critical power increases as the neutrino luminosity is decreased, which was also the case in 1D. Although we cannot conduct an analysis in 2D similar to the one given for 1D models in section 3, we infer that this is likely due to the reduced efficiency of the acoustic heating caused by the enhanced neutrino emissions by raised temperatures by the secondary shock waves as well as due to the partial reflection of acoustic waves.

5 summary and discussions

We performed 1D and 2D simulations of shock revival in the supernova core with the acoustic waves from PNS being taken into account phenomenologically. The 1D simulations were intended to capture the essential physics in the energy deposition by acoustic waves. For various combinations of mass accretion rate and neutrino luminosity, we obtained the critical amplitudes of an acoustic wave that devide successful from failed shock revival, and then drew the critical surface in the space of the neutrino luminosity $L_{\nu}$, mass accretion rate $\dot{M}$, and acoustic amplitude $\delta$. In the successful models, the primary stalled shock is repeatedly hit by secondary shocks, into which acoustic waves steepen. As a consequence, the primary shock develops overstable oscillations and eventually revives. It is a combination of neutrino heating and acoustic power that gives rise to shock revival, however. In this sense, the mechanism considered in this paper may be referred to as “hybrid.”

In order to discuss the energetics quantitatively, we extended the Myers energy flux for finite-amplitude acoustic waves to incorporate neutrino contributions. By redrawing the critical surface with the acoustic power estimated from the extended Myers flux, we found that the sum of the neutrino heating and the acoustic power is almost constant on the critical surface for given mass accretion rates, with the decrease in the former being nearly compensated for by the increase in the latter. It hence appears that the critical luminosity in the neutrino-heating mechanism can be replaced by this sum in the neutrino-acoustic hybrid mechanism. For low neutrino luminosities, however, more acoustic powers seem to be required than to merely compensate for the decrease in the neutrino heating. This is because large-amplitude acoustic waves tend to become strong shock waves, resulting in higher temperatures and thus spending the deposited energies on enhanced neutrino emissions. Incidentally, we applied two diagnostics for shock revival, which are commonly used in the context of the neutrino-heating mechanism. Neither of them was found to be useful in the present mechanism.

Next, we ran 2D simulations. We note that the acoustic mechanism is intrinsically multidimensional, since there are no spherical $g$-mode oscillations, which are the emitters of acoustic waves. Although the critical amplitudes of acoustic wave are derived from the simulations, they are not appropriate for the comparison with 1D results. The acoustic power is more suitable, and indeed, the critical surface in 2D is always lower in the acoustic power than the 1D counterpart. This is due to the multidimensional fluid instabilities, which are forcibly excited by acoustic waves and enhance the neutrino heating.

With the critical surface thus obtained, we revisited the numerical results of Burrows et al. (2006) from the viewpoint of energetics. Although the mass accretion rate and neutrino luminosity both vary with time in Burrows et al. (2006), we detected a representative combination of them as $\dot{M}\sim 0.1\,M_{\odot}\,{\rm s^{-1}}$ and $L_{\nu}\sim 2.0\times 10^{52}\,{\rm erg\,s^{-1}}$. The corresponding acoustic power is $\sim 4\times 10^{51}\,{\rm erg\,s^{-1}}$, which is by inspection much larger than the critical acoustic power obtained in this paper. Our model on the critical surface with $\dot{M}\sim 0.2\,M_{\odot}\,{\rm s^{-1}}$ and $L_{\nu}\sim 2.0\times 10^{52}\,{\rm erg\,s^{-1}}$ requires an acoustic power of $\sim 9\times 10^{50}\,{\rm erg\,s^{-1}}$ for shock revival. Note that the critical acoustic power decreases with the mass accretion rate. It is hence not surprising that Burrows et al. obtained explosions via the acoustic mechanism for such high acoustic powers in their simulations. Incidentally, the acoustic powers estimated by Yoshida et al. (2007) are close to the critical surface. On the other hand, the theory of Weinberg & Quataert (2008) predicts much smaller acoustic powers from the saturated $g$-mode oscillations, which are certainly insufficient for shock revival.

The critical amplitude of density perturbation should not exceed unity, since the density would become negative otherwise. This may give another interesting constraint. Linearly extrapolating the critical surface to $\delta=1$, we obtain the critical curve that runs through the points with $(\dot{M},L_{\nu})=(1.0\,M_{\odot}\,{\rm s^{-1}},\sim 2\times 10^{52}\,{\rm erg\,s^{-1}})$, $(0.6\,M_{\odot}\,{\rm s^{-1}},\sim 1\times 10^{52}\,{\rm erg\,s^{-1}})$, and $(0.2\,M_{\odot}\,{\rm s^{-1}},\sim 1\times 10^{52}\,{\rm erg\,s^{-1}})$. This may imply that models below these luminosities do not explode even if the acoustic waves power the stalled primary shock. We note that we considered only sinusoidal perturbations at the inner boundary in this paper. It may hence be true that the above estimate may not hold for other types of perturbations with different angular modes, oscillation periods, and so on, but we believe that it will not be changed by the order, since the density must be positive and hence there clearly exists a maximum fluctuation amplitude, regardless of the details of disturbances.

Although it is interesting to see that the extended Myers flux derived here serves well in estimating the energy flux of finite-amplitude acoustic waves and the critical surface obtained in the simplified settings provides useful conditions for shock revival, which seem to be consistent with realistic simulations, there are some caveats in the above assessment. First, our models neglected the turbulence in the postshock flows that should have existed before the injection of acoustic waves, since otherwise the $g$-mode oscillations could not have been excited in PNS in the first place. This problem should be also important in estimating the acoustic power, although it is beyond the scope of this paper. Second, our simulations are 2D at most. It is well known, however, that turbulence properties are qualitatively different between 2D and 3D. Since the inverse cascade develops in the 2D turbulence (Kraichnan, 1967), smaller turbulent eddies will be produced in 3D than in 2D (e.g., Couch, 2013; Takiwaki et al., 2014; Melson et al., 2015), which may result in reduced neutrino heating as well as weaker PNS oscillations, thus shifting the critical surface somewhat upward in 3D in terms of the acoustic power. Further investigations on these issues are certainly warranted. Although the neutrino heating is the most favored mechanism of CCSNe at present, we should not forget alternatives at any time.

The authors thank Masataka Ogane for fruitful disucussion. A.H. also thanks Bernhard Müller for useful suggestion and discussion. Numerical computations were carried out on Cray XC30 at Center for Computational Astrophysics, National Astronomical Observatory of Japan, and the supercomputer system A at High Energy Accelerator Research Organization (KEK, support by the Large Scale Simulation Program No. 14/15-17 (FY2014-2015), No. 15/16-08 (FY2015-2016), and No. 16/17-11 (FY2016-2017)). A.H. is supported by Advanced Leading Graduate Course for Photon Science (ALPS) in the University of Tokyo. H.N. was supported in part by JSPS Postdoctoral Fellowships for Research Abroad No. 27-348, and he was partially supported at Caltech through NSF award No. TCAN AST-1333520. This work was supported by the Grant-in-Aid for Scientific Research (B) (No. 16H03986) and Grant-in-Aid for Innovative Areas (No. 24103006) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan.

Appendix A Myers corollary to the theory of energy conservation

In this appendix, we extend the Myers corollary to the theory of energy conservation for finite-amplitude perturbations (Myers, 1986, 1991) and derive equations (6)–(LABEL:eq:daco). Although the content given below is almost a summary of Myers’ work (Myers, 1991), except for the incorporation of neutrino heating, we believe that it is worth reviewing.

Our basic equations are

[TABLE]

Equation (5) is used for the gravitational potential, where $\boldsymbol{m}=\rho\boldsymbol{v}$ is the mass flux, $\boldsymbol{\zeta}:={\boldsymbol{\omega}}\boldsymbol{\times}{\boldsymbol{v}}:=({\boldsymbol{\nabla}}\boldsymbol{\times}{\boldsymbol{v}})\boldsymbol{\times}{\boldsymbol{v}}$ and $\boldsymbol{\omega}$ is the vorticity; $C$, $\boldsymbol{L}$, $S$, and $G$ are the shorthand notations for the left-hand sides of the equations of continuity, Euler, entropy, and electron fractions, respectively; the corresponding right-hand side of the last three equations are denoted by $\boldsymbol{\lambda}$, $\sigma$, and $\gamma$, respectively. Equtions (A1)–(A4) are equivalent to equations (1)–(4) in the main body, except that we include the momentum transfer from the neutrinos, $\boldsymbol{M}$, in order to make the derivation below as general as possible. The relevant thermodynamic relations are

[TABLE]

where $T$, $s$, and $h=e+P/\rho$ are the temperature, entropy, and specific enthalpy, respectively; $\mu$ is defined as $\mu=\mu_{\rm e}+\mu_{\rm p}-\mu_{\rm n},$ with $\mu_{\rm e,p,n}$ being the chemical potentials of electron, proton, and neutron, respectively; $m_{\rm u}$ is the atomic mass unit. Using these relations, we can derive from our basic equations (A1)–(A4) the energy conservation law cast in the following form,

[TABLE]

Here the specific stagnation enthalpy (or the Bernoulli function) $h+\text@frac{1}{2}\boldsymbol{v}^{2}$ is denoted by $H$. We note that the identity $\boldsymbol{v}{\boldsymbol{\cdot}}\boldsymbol{\zeta}=\boldsymbol{v}{\boldsymbol{\cdot}}(\boldsymbol{\omega}\boldsymbol{\times}\boldsymbol{v})=0$ is used.

Consider a perturbative expansion of a quantity $q$ as follows: $q({\boldsymbol{r}},t)=q_{0}({\boldsymbol{r}})+\sum_{n=1}^{\infty}\delta^{n}q_{n}({\boldsymbol{r}},t)$. The subscript [math] denotes the unperturbed state with no disturbance, whereas the subscript $n$ represents the $n$th order perturbation. We note that the gravitational potential $\Phi$ is assumed to be determined by PNS in this paper and is hence not perturbed. Applying this expansion to each quantity in the equations given above and equating the terms of the same order, we obtain a sequence of equation systems that govern the perturbations at each order. We attach the subscript $i$ to the the shorthand notations introduced above,

[TABLE]

to represent the $i$th order perturbations to them. We also expand the energy-conservation equation (A7):

[TABLE]

and

[TABLE]

We rewrite equations further in the following form, expanding again thermodynamic quantities using the relations in equations (A5) and (A6), and using Maxwell’s relations obtained also from the same relations:

[TABLE]

[TABLE]

and

[TABLE]

where $E_{2}$, $\boldsymbol{F}_{2}$, and $D_{2}$ are given as

[TABLE]

[TABLE]

and

[TABLE]

We note that the zeroth- and first-order equations are trivially satisfied due to equations (A8). This is not true of the second-order equation, however, and the last line of equation (A14) remains, which may be interpreted as the energy conservation law for the first-order perturbation for the reasons given below. We note that $E_{2}$, $\boldsymbol{F}_{2}$ and $D_{2}$ contain only zeroth- and first-order quantities.

The first three terms on the right-hand side of equation (A15) are the well-known representation for the acoustic energy density in the homentropic flow. If the flow is not uniform in entropy, an extra contribution is expected from the $T{\rm d}s$ term in the thermodynamic relations. This is the origin of the fourth term in equation (A15). We note that we need to consider the product of the first-order perturbations in entropy and temperature so that it should not vanish after the average over the oscillation period, and indeed, $(\partial T/\partial s)_{P,Y_{\rm e}}s_{1}$ and $(\partial T/\partial Y_{\rm e})_{s,P}Y_{\rm e1}$ are the changes in temperature induced by the change in entropy and electron fraction, respectively. In a similar way, we sholud take into account the changes in both $Y_{\rm e}$ and $\mu$ for the contribution from the $\mu/m_{\rm u}{\rm d}Y_{\rm e}$ term. This is expressed as the fifth term in equation (A15). The same considerations can be applied to $\boldsymbol{F}_{2}$. The first term on the right-hand side in equation (A16) is again the well-known representation for the acoustic energy flux in the homentropic flow, in which inhomogeneities in the flow velocity are taken into account, whereas the second term originates from the changes in entropy, temperature, $Y_{\rm e}$, and $\mu$. The last term, $D_{2}$, is a residual that contains everything other than those included in $E_{2}$ and $\boldsymbol{F}_{2}$, and as such, it is admittedly the most difficult to interpret. Given the fact that $E_{2}$ and $\boldsymbol{F}_{2}$ can be regarded as the energy density and flux, however, one may interpret $D_{2}$ as a dissipation term. The neutrino cooling represented by the fourth term, $-T_{1}\sigma_{1}$, in $D_{2}$ clearly works that way: if $T_{1}$ is positive (negative), then the neutrino emission will be enhanced (suppressed), leading to negative (positive) values of $\sigma_{1}=\{(Q-\mu\gamma/m_{\rm u}-\boldsymbol{v}{\boldsymbol{\cdot}}\boldsymbol{M})/T\}_{1}$; this makes $-T_{1}\sigma_{1}$ always positive, implying that the neutrino cooling tends to reduce the perturbation energy. Other terms in $D_{2}$ may not be so easy to interpret, but we refer to $D_{2}$ as the dissipation term in this paper.

What we have done so far is an ordinary perturbative expansion of the basic equations up to the second order. The essential idea of Myers now plays a role. We first remark that the energy conservation equations up to the second order have the same structure. In his theory, Myers surmised that this is true to all orders and recast the exact energy conservation law into the following form:

[TABLE]

where $E_{\rm dis}$, $\boldsymbol{F}_{\rm dis}$, and $D_{\rm dis}$ are given as

[TABLE]

He then interpreted them as the density, flux, and dissipation of the energy for not necessarily small disturbances. We note that we modified the original expression to incorporate the neutrino heating in equation (A18). It is admittedly difficult to justify the interpretation unambiguously, but it may be somewhat comforting to point out (i) that in the limit of small perturbations, $E_{\rm dis}$, $\boldsymbol{F}_{\rm dis}$, and $D_{\rm dis}$ are reduced to $E_{2}$, $\boldsymbol{F}_{2}$, and $D_{2}$, respectively, and (ii) that the resulting equation for $E_{\rm dis}$, $\boldsymbol{F}_{\rm dis}$, and $D_{\rm dis}$ takes the conservative form. Ignoring the momentum transfer from neutrinos to matter, $\boldsymbol{M}$ and $\boldsymbol{M}_{0}$, which is well justified for the models considered in this paper, we finally obtain equations (6)–(LABEL:eq:daco). We note that the neutrino heating is accounted for by the first term in the second line of equation (A21).