The Boltzmann-radiation-hydrodynamics Simulations of Core-collapse Supernovae with Different Equations of State: The Role of Nuclear Composition and the Behavior of Neutrinos

Akira Harada; Hiroki Nagakura; Wakana Iwakami; Hirotada Okawa; Shun Furusawa; Kohsuke Sumiyoshi; Hideo Matsufuru; Shoichi Yamada

doi:10.3847/1538-4357/abb5a9

1. Introduction

Core-collapse supernovae (CCSNe) are considered to be the explosive death of massive stars. The energy source of the explosion is the gravitational energy released when the stellar iron core collapses to form a neutron star (NS; Baade & Zwicky 1934). It amounts to ∼10⁵³ erg. CCSNe produce some heavy elements in the universe. In addition, the recent discovery and understanding of binary NS mergers imply that some r-process elements are produced there (Abbott et al. 2017a, 2017b; Tanaka et al. 2017). Therefore, it is important to understand CCSNe to explain elemental evolution in the universe.

The standard scenario for the explosion is the neutrino heating mechanism (see Janka 2012 for a review). The stellar iron core eventually collapses. When the central density goes beyond the nuclear density, the matter suddenly becomes stiff, and the bounce shock is formed. This bounce shock propagates outward, consuming its energy by photodissociating the accreting heavy nuclei such as iron, and it finally stops propagating. SN modelers have been investigating how to revive this stalled shock. The leading hypothetical mechanism for shock revival is the neutrino heating mechanism (Wilson 1985). After the core bounce, a central hydrostatic object called a protoneutron star (PNS) is formed. The PNS is still hot and contains a lot of protons. The energy and the lepton number of the PNS are carried away by the neutrinos to evolve into the NS. In the neutrino heating mechanism, these emitted neutrinos are absorbed by matter just behind the shock and heat up the matter. Then, the shock restarts its propagation.

The progress of CCSN simulations started with 1D spherical symmetry. From the first stellar core-collapse simulation by Colgate & Johnson (1960), CCSN simulations have been continuously improved. The works of Liebendörfer et al. (2001b) and Sumiyoshi et al. (2005) employed very sophisticated simulation codes, which solve the Boltzmann equation for neutrino transport with general relativity (Mezzacappa & Bruenn 1993; Yamada 1997; Yamada et al. 1999; Liebendörfer et al. 2001a). These sophisticated simulations did not show a successful shock revival. Nowadays, it is concluded that spherically symmetric CCSNe do not explode except for a special progenitor (Kitaura et al. 2006).

From the middle of the 1990s, 2D axisymmetric simulations have been performed (Herant et al. 1994; Burrows et al. 1995; Buras et al. 2003, 2006; Marek & Janka 2009; Müller et al. 2012; Bruenn et al. 2013; Dolence et al. 2015; Bruenn et al. 2016; Summa et al. 2016; O'Connor & Couch 2018). On this stage, shock revivals are found in simulations. These simulations revealed that multidimensional effects such as turbulence help the neutrino heating and shock revival.

The Boltzmann neutrino transport, which is employed in, e.g., Yamada et al. (1999), Burrows et al. (2000), and Liebendörfer et al. (2001a) under spherical symmetry, requires a lot of numerical resources; thus, almost all multidimensional simulations so far use approximate neutrino transport. For example, the flux-limited diffusion scheme (Burrows et al. 2006, 2007), two-moment scheme (M1 closure: Kuroda et al. 2012, Just et al. 2015; variable Eddington factor: Rampp & Janka 2002), and isotropic diffusion source approximation (Liebendörfer et al. 2009) are utilized. In addition, the ray-by-ray(-plus) approach (Rampp & Janka 2002; Buras et al. 2006) is sometimes used. As the employed approximation schemes are different among SN modelers, the outcomes such as the explodability and explosion energy of the simulations are also different. Some collaboration works have reported on a code comparison project, and they show basic agreement in 1D simulations or in the early stages of the postbounce dynamics in multiple dimensions (Yamada et al. 1999; Liebendörfer et al. 2005; Müller et al. 2010; Skinner et al. 2016; Cabezón et al. 2018; Just et al. 2018; O'Connor et al. 2018; Glas et al. 2019; Pan et al. 2019). They guarantee that the basic ingredients are correctly included in the codes, but the origin of some differences is yet to be answered.

Even though the shock revives in 2D simulations, there are still some problems. The simulated explosion energy is about 10⁵⁰ erg, at least at several hundreds of milliseconds after the core bounce. This is much smaller than the observed value of 10⁵¹ erg. Current development of computational resources allows SN modelers to perform long-term simulations exceeding 1 s after the core bounce, and some of these simulations indicate that the continued evolution of the explosion energy may reach the observed explosion energy (Bruenn et al. 2013, 2016). On the other hand, such long-term and relatively slow heating may be problematic in explaining another observable, the amount of ${}^{56}\mathrm{Ni}$ . Suwa et al. (2019) suggested from the analytic model calibrated by numerical simulations that the observed amount of ${}^{56}\mathrm{Ni}$ is reproduced if the heating rate exceeds ${ \mathcal O }({10}^{51})\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ , about 10 times larger than what is expected in simulations. Therefore, the problem related to the explosion energy, or the heating rate, still remains.

Recent SN modelers are going to 3D simulations, thanks to increasing computational power (Fryer & Warren 2002; Takiwaki et al. 2012; Hanke et al. 2013; Tamborra et al. 2013, 2014; Takiwaki et al. 2014; Lentz et al. 2015; Melson et al. 2015; Müller 2015; Kuroda et al. 2016; Ott et al. 2018; Vartanyan et al. 2019). This is because nature is, of course, 3D, and hence 3D simulations are appropriate for comparing to observations such as the anisotropic distribution of ejected elements, the aspherical morphology of SN remnants, and so on (Wongwathanarat et al. 2017; Ono et al. 2020). However, although some simulations show shock revival, problems similar to those in 2D simulations still remain.

Using SN simulations, the roles of microphysics such as nuclear EOSs are investigated. Various nuclear EOS models have been developed for SN simulations (e.g., Lattimer & Swesty 1991; Shen et al. 1998; Furusawa et al. 2011, 2013, 2017a, 2017b; Togashi et al. 2017; Sumiyoshi et al. 2019). The effects of these EOSs on the SN dynamics have been discussed (e.g., Marek et al. 2009; Couch 2013; Suwa et al. 2013; Fischer et al. 2014; Schneider et al. 2019); those on multimessenger signals, i.e., neutrinos and gravitational waves, have been examined (e.g., Pan et al. 2018; Nakazato & Suzuki 2020). These results imply that soft EOSs such as the Lattimer–Swesty EOS (Lattimer & Swesty 1991) tend to show robust shock revival.

We should note that the failure to revive shocks in 1D simulations is concluded from the general relativistic Boltzmann-radiation-hydrodynamics simulations. By omitting the approximations in the governing equations, we can robustly assess the effects of the uncertainty in the input physics: the nuclear EOSs and microphysical reactions. Because there still remain some problems in 2D simulations, it is important to use first-principles simulations, which do not use approximations for the governing equations, to understand the origin of the problem.

We utilized improved computational resources to solve the Boltzmann-radiation hydrodynamics under axisymmetry, instead of performing 3D simulations. The development of the Boltzmann-radiation-hydrodynamics code is reported in Sumiyoshi & Yamada (2012) and Nagakura et al. (2014, 2017) and briefly explained in Section 2. We employed the S_N method for neutrino transport, and the comparison with Monte Carlo transport is presented in Richers et al. (2017). Furthermore, Chan & Mueller (2020) also developed a multidimensional Boltzmann-neutrino-transport solver, indicating the necessity of the exact treatment of neutrino transport.

We report the results of thorough analyses of the simulations using the Boltzmann-radiation-hydrodynamics code with different nuclear EOSs in this paper. We have reported the results of the simulations with the code in Nagakura et al. (2018) and Harada et al. (2019) so far. The initial work, Nagakura et al. (2018), demonstrated the novel features of our Boltzmann-radiation-hydrodynamics code with simulations with different nuclear EOSs: the differences in the dynamical features were discussed briefly; the momentum space distributions of neutrinos were deeply analyzed, but only limited time snapshots and spatial points were considered. Therefore, thorough analyses of the dynamical features and neutrino distributions at various times and spatial regions are necessary; we analyzed the effects of the different EOSs and the behavior of the neutrinos deeply and widely.

The simulation time is extended to confirm the fate of the shock. In the previous paper, the simulation time was not enough to see whether the shock is revived or not. Thus, we ran additional simulation continuing from the previous paper and show the results here. Thanks to the extended simulation, the shock revival is more noticeable than in the previous paper.

This paper is structured as follows. In Section 2, a summary of our code and the models are presented. Next, in Section 3, we discuss the comparison of the postbounce hydrodynamics and some observable neutrino quantities from the viewpoint of the EOSs. Then, in Section 4, we give detailed analyses of the neutrino quantities: the neutrino fluxes and the Eddington tensors. Finally, we conclude our paper in Section 5. The units with $c=G={\hslash }=1$ are considered in this paper unless otherwise stated. The signature of the spacetime metric is $(-+++)$ . Greek and Latin indices run over 0–3 for spacetime and 1–3 for space, respectively.

2. Numerical Modeling

The Boltzmann-radiation-hydrodynamics code used in this paper solves the directly discretized Boltzmann equation (the so-called S_N method). The details of the code are described in Sumiyoshi & Yamada (2012) and Nagakura et al. (2014, 2017). The Boltzmann equation for neutrinos with a general metric is cast in the conservative form (Shibata et al. 2014):

$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{1}{\sqrt{-g}}\displaystyle \frac{\partial }{\partial {x}^{\alpha }}\left|\Space{0ex}{3.25ex}{0ex},| ,{}_{{q}_{i}}\right.\left[\left({e}_{(0)}^{\alpha }+\displaystyle \displaystyle \sum _{i=1}^{3}{{\ell }}_{(i)}{e}_{i}^{\alpha }\right)\sqrt{-g}f\right]\\ -\displaystyle \frac{1}{{\epsilon }^{2}}\displaystyle \frac{\partial }{\partial \epsilon }({\epsilon }^{3}f{\omega }_{(0)})+\displaystyle \frac{1}{\sin {\theta }_{\nu }}\displaystyle \frac{\partial }{\partial {\theta }_{\nu }}(\sin {\theta }_{\nu }f{\omega }_{({\theta }_{\nu })})\\ +\displaystyle \frac{1}{{\sin }^{2}{\theta }_{\nu }}\displaystyle \frac{\partial }{\partial {\phi }_{\nu }}(f{\omega }_{({\phi }_{\nu })})={S}_{\mathrm{rad}},\end{array}\end{eqnarray} \tag{ 1 }$

$\begin{eqnarray}&&{{\ell }}_{(i)}:= (\cos {\theta }_{\nu },\sin {\theta }_{\nu }\cos {\phi }_{\nu },\sin {\theta }_{\nu }\sin {\phi }_{\nu }),\end{eqnarray} \tag{ 2 }$

$\begin{eqnarray}&&{\omega }_{(0)}:= {\epsilon }^{-2}{p}^{\alpha }{p}_{\beta }{{\rm{\nabla }}}_{\alpha }{e}_{(0)}^{\beta },\end{eqnarray} \tag{ 3 }$

$\begin{eqnarray}&&{\omega }_{({\theta }_{\nu })}:= \displaystyle \sum _{i=1}^{3}{\omega }_{i}\displaystyle \frac{\partial {{\ell }}_{(i)}}{\partial {\theta }_{\nu }},\end{eqnarray} \tag{ 4 }$

$\begin{eqnarray}&&{\omega }_{({\phi }_{\nu })}:= \displaystyle \sum _{i=2}^{3}{\omega }_{i}\displaystyle \frac{\partial {{\ell }}_{(i)}}{\partial {\phi }_{\nu }},\end{eqnarray} \tag{ 5 }$

$\begin{eqnarray}&&{\omega }_{i}:= {\epsilon }^{-2}{p}^{\alpha }{p}_{\beta }{{\rm{\nabla }}}_{\alpha }{e}_{(i)}^{\beta },\end{eqnarray} \tag{ 6 }$

where x^α, g, ${e}_{(\mu )}^{\alpha }$ , f, p^α, $\epsilon := -{p}_{\mu }{e}_{(0)}^{\mu }$ , ${\theta }_{\nu }$ , ϕ_ν, and S_rad are the spatial coordinates, the determinant of the metric ${g}_{\mu \nu }$ , the tetrad bases, the distribution function, the four-momentum of the neutrino, the energy of the neutrino, the zenith and azimuthal angle of the neutrino flight direction with respect to the tetrad basis, and the collision term, respectively. The zeroth tetrad basis ${e}_{(0)}^{\mu }$ is chosen to be the normal vector to the spatial hypersurface ${n}^{\mu }=({\alpha }^{-1},-{\alpha }^{-1}{\beta }^{i})$ , where α and βⁱ are the lapse function and the shift vector. According to the 3 + 1 decomposition, the spacetime metric is decomposed as

$\begin{eqnarray}&&{{ds}}^{2}={g}_{\mu \nu }{{dx}}^{\mu }{{dx}}^{\nu }\end{eqnarray} \tag{ 7 }$

$\begin{eqnarray}&&=-{\alpha }^{2}{{dt}}^{2}+{\gamma }_{{ij}}({{dx}}^{i}+{\beta }^{i}{dt})({{dx}}^{j}+{\beta }^{j}{dt}),\end{eqnarray} \tag{ 8 }$

where ${\gamma }_{\mu \nu }={g}_{\mu \nu }+{n}_{\mu }{n}_{\nu }$ is the spatial metric. Because we choose polar coordinates $(r,\theta ,\phi )$ , the other tetrad bases are chosen to be

$\begin{eqnarray}&&{e}_{(1)}^{\alpha }={\gamma }_{{rr}}^{-1/2}{\partial }_{r},\end{eqnarray} \tag{ 9 }$

$\begin{eqnarray}&&{e}_{(2)}^{\alpha }=-\displaystyle \frac{{\gamma }_{r\theta }}{\sqrt{{\gamma }_{{rr}}({\gamma }_{{rr}}{\gamma }_{\theta \theta }-{\gamma }_{r\theta }^{2})}}{\partial }_{r}+\sqrt{\displaystyle \frac{{\gamma }_{{rr}}}{{\gamma }_{{rr}}{\gamma }_{\theta \theta }-{\gamma }_{r\theta }^{2}}}{\partial }_{\theta },\end{eqnarray} \tag{ 10 }$

$\begin{eqnarray}&&{e}_{(3)}^{\alpha }=\displaystyle \frac{{\gamma }^{r\phi }}{\sqrt{{\gamma }^{\phi \phi }}}{\partial }_{r}+\displaystyle \frac{{\gamma }^{\theta \phi }}{\sqrt{{\gamma }^{\phi \phi }}}{\partial }_{\theta }+\sqrt{{\gamma }^{\phi \phi }}{\partial }_{\phi },\end{eqnarray} \tag{ 11 }$

where ${\partial }_{i}$ are the coordinate bases of the vector. Although these expressions can be applied to the general, curved spacetime, we only focus on the flat spacetime with drifting coordinates. The coordinates move with the PNS for their centers to match. This is realized by the flat spacetime with the nonzero shift vector. This frame is nothing but an acceleration frame. The method to determine the shift vector is discussed in Nagakura et al. (2017). Because the coordinate drift is not significant in the current models as discussed in Section 3.1, we just call it the laboratory frame. Although this equation describes the time evolution of the six-dimensional distribution function, we impose axisymmetry due to limited computational resources, and as a consequence, the distribution function is a function in the five-dimensional phase space (two in configuration space and three in momentum space).

The neutrino reactions are based on the standard set (Bruenn 1985): the (anti)electron capture on the nucleon and nuclei, the scattering off the nucleon and nuclei, and the pair neutrino production from the electron–positron pair. However, there are several modifications: the electron capture on heavy nuclei is updated according to Juodagalvis et al. (2010), Langanke & Martínez-Pinedo (2000), and Langanke et al. (2003); the inelastic scattering off electrons is incorporated; and nucleon–nucleon bremsstrahlung is implemented. Contrary to the recent work (Nagakura et al. 2019b), the electron capture on light nuclei is not incorporated. Instead, it is artificially replaced by the weak interactions with the free nucleons, which are the constituents of light nuclei. This treatment is similar to the DV112 model in Nagakura et al. (2019a). We set the same rates for the reactions involving the heavy-lepton-type neutrinos ${\nu }_{\mu }$ , $\bar{{\nu }_{\mu }}$ , ${\nu }_{\tau }$ , and $\bar{{\nu }_{\tau }}$ (but see Bollig et al. 2017). They are collectively denoted by ${\nu }_{x}$ . Therefore, we solve the Boltzmann equations for ${\nu }_{{\rm{e}}}$ , $\bar{{\nu }_{{\rm{e}}}}$ , and ${\nu }_{x}$ .

As for the hydrodynamics, we solve the Newtonian hydrodynamics equations in the acceleration frame:

$\begin{eqnarray}&&{\partial }_{t}(\sqrt{-g}\rho )+{\partial }_{i}(\sqrt{-g}\rho {v}^{i})=0,\end{eqnarray} \tag{ 12 }$

$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}(\sqrt{-g}\rho {v}_{r})+{\partial }_{i}\left(\sqrt{-g}(\rho {v}_{r}{v}^{i}+p{\delta }_{r}^{i})\right)\\ =\,\sqrt{-g}\rho \left(-{\partial }_{r}\psi +r{\left({v}^{\theta }\right)}^{2}+r{\sin }^{2}\theta {\left({v}^{\phi }\right)}^{2}+\displaystyle \frac{2p}{r\rho }\right)\\ -\sqrt{-g}{G}_{r}-\sqrt{-g}\rho \dot{{\beta }_{r}},\end{array}\end{eqnarray} \tag{ 13 }$

$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}(\sqrt{-g}\rho {v}_{\theta })+{\partial }_{i}\left(\sqrt{-g}(\rho {v}_{\theta }{v}^{i}+p{\delta }_{\theta }^{i})\right)\\ =\,\sqrt{-g}\rho \left(-{r}^{2}{\partial }_{\theta }\psi +\sin \theta \cos \theta {\left({v}^{\phi }\right)}^{2}+\displaystyle \frac{p\cos \theta }{\rho \sin \theta }\right)\\ -\sqrt{-g}{G}_{\theta }-\sqrt{-g}\rho \dot{{\beta }_{\theta }},\end{array}\end{eqnarray} \tag{ 14 }$

$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}(\sqrt{-g}\rho {v}_{\phi })+{\partial }_{i}\left(\sqrt{-g}(\rho {v}_{\phi }{v}^{i}+p{\delta }_{\phi }^{i})\right)\\ =\,-\sqrt{-g}\rho {\partial }_{\phi }\psi -\sqrt{-g}{G}_{\phi }-\sqrt{-g}\rho \dot{{\beta }_{\phi }},\end{array}\end{eqnarray} \tag{ 15 }$

$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}\left(\sqrt{-g}(e+\displaystyle \frac{1}{2}\rho {v}^{2})\right)+{\partial }_{i}\left(\sqrt{-g}(e+p+\displaystyle \frac{1}{2}\rho {v}^{2}){v}^{i}\right)\\ =\,-\sqrt{-g}\rho {v}^{i}{\partial }_{i}\psi -\sqrt{-g}{G}^{t}-\sqrt{-g}\rho {v}^{i}\dot{{\beta }_{i}},\end{array}\end{eqnarray} \tag{ 16 }$

and

$\begin{eqnarray}&&{\partial }_{t}(\sqrt{-g}\rho {Y}_{{\rm{e}}})+{\partial }_{i}(\sqrt{-g}\rho {Y}_{{\rm{e}}}{v}^{i})=-\sqrt{-g}({{\rm{\Gamma }}}_{{\nu }_{{\rm{e}}}}-{{\rm{\Gamma }}}_{\bar{{\nu }_{{\rm{e}}}}}).\end{eqnarray} \tag{ 17 }$

The symbols ρ, vⁱ, p, e, ${Y}_{{\rm{e}}}$ , and ψ represent the density, the velocity, the pressure, the internal energy, the electron fraction, and the Newtonian gravitational potential, which obeys

$\begin{eqnarray}&&{\rm{\Delta }}\psi =4\pi \rho ,\end{eqnarray} \tag{ 18 }$

respectively. In these Newtonian hydrodynamics equations, we consider the flat spacetime metric, i.e., $\sqrt{-g}={r}^{2}\sin \theta$ . The exchange of energy and momentum between matter and neutrinos is

$\begin{eqnarray}&&{G}^{\mu }=\int {p}^{\mu }\epsilon {S}_{\mathrm{rad}}{{dV}}_{p},\end{eqnarray} \tag{ 19 }$

where ${{dV}}_{p}$ is the invariant volume element in momentum space. The reaction rate of the electron-type neutrinos and antineutrinos are

$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\nu }_{{\rm{e}}}/\bar{{\nu }_{{\rm{e}}}}}={m}_{{\rm{u}}}\int \epsilon {S}_{\mathrm{rad},{\nu }_{{\rm{e}}}/\bar{{\nu }_{{\rm{e}}}}}{{dV}}_{p},\end{eqnarray} \tag{ 20 }$

where ${m}_{{\rm{u}}}$ and ${S}_{\mathrm{rad},{\nu }_{{\rm{e}}}/\bar{{\nu }_{{\rm{e}}}}}$ are the atomic mass unit and the collision term for the indicated neutrino species.

The EOSs adopted in this paper are the Furusawa–Shen (FS; Furusawa et al. 2011, 2013) EOS and the Lattimer–Swesty (LS; Lattimer & Swesty 1991) EOS. The FS EOS is based on the Shen EOS (Shen et al. 1998), and nuclear statistical equilibrium (NSE) is considered for the ensemble of nuclei in order to calculate the thermodynamical and statistical properties of nonuniform matter. The Shen EOS, the basis of the FS EOS, models the strong interaction by the relativistic mean field theory with the TM1 parameter set. On the other hand, the LS EOS is based on the liquid drop model of the nuclei and the Skyrm-type interaction. As for the composition, Lattimer & Swesty (1991) assume that the heavy nuclei are represented by a single nuclear species (single nuclear approximation, SNA), and only the alpha particle is considered as the light nuclei. Among the EOSs offered by Lattimer & Swesty (1991), we choose the EOS with the incompressibility parameter $K=220\,\mathrm{MeV}$ . Because our simulation code uses tabulated EOSs, we convert the subroutine EOS originally provided by Lattimer & Swesty (1991) to the tabulated EOS.⁹ For a given set of density, temperature, and electron fraction, the electron capture rates on the heavy nuclei are calculated based on the heavy nucleus abundance of the FS EOS. The same rate as a function of the thermodynamic quantities is used for both models with the FS and LS EOSs. Although the more sophisticated Furusawa–Togashi (FT) EOS is available (Togashi et al. 2017; Nagakura et al. 2019a), we do not employ it in this paper. Because the weak interaction rates are different between the EOSs in this paper and the FT EOS, we could not focus on the properties of nuclear matter if we compared the simulations with the LS/FS and FT EOSs.

The formulation described above is numerically evolved with the following schemes. The numerical flux of the Boltzmann equation is evaluated by the combinations of the upwind and central difference scheme according to the local mean free path (MFP). Specifically, we carefully discretize the Boltzmann equation to guarantee the steady-state infinite homogenous solution, i.e., a constant distribution function with respect to the spacetime and momentum. The equation is evolved semi-implicitly, thus we need to solve large linear coupled equations. The Bi-CGSTAB method (Saad 2003) is utilized for the matrix inversion with the point Jacobi preconditioner. For the hydrodynamics equations, the Harten–Lax–van Leer (HLL) scheme (Harten et al. 1983) with piecewise-parabolic interpolation (Colella & Woodward 1984) determines the numerical flux. The second-order Runge–Kutta method is adopted for the time integration. The Poisson equation for the gravitational potential ψ is solved by direct multiplication of the inverse matrix of the discretized Laplacian operator. The inverse matrix is calculated using the MICCG method (Nagakura et al. 2011).

We run the two-dimensional axisymmetric simulations from the beginning of the prompt convection. The progenitor model is the $11.2\,{M}_{\odot }$ model in Woosley et al. (2002). Until the negative entropy gradient is formed, one-dimensional simulations are performed from the onset of the collapse. When the negative entropy gradient is formed at ${t}_{\mathrm{pb}}\sim 0.6\,\mathrm{ms}$ , where ${t}_{\mathrm{pb}}$ is the postbounce time, the hydrodynamical quantities and neutrino distributions are mapped to the two-dimensional simulations with seed perturbations of $0.1 \%$ in the radial velocity inside the region $30\leqslant r\leqslant 50\,\mathrm{km}$ . Because the shock radius is $\sim 30\,\mathrm{km}$ at that time, this is the preshock velocity perturbation. The computational domain is $0\leqslant r\leqslant 5000\,\mathrm{km}$ and $0\leqslant \theta \leqslant \pi$ . The radial and zenith coordinates are divided into 384 and 128 grids, respectively. Although our code can solve Equation (15) simultaneously (Harada et al. 2019; Iwakami et al. 2020), we do not solve it for the nonrotating and axially symmetric simulation in this paper. The neutrino energies range from 0 to $300\,\mathrm{MeV}$ and are divided into 20 grids. The ${\theta }_{\nu }$ and ${\phi }_{\nu }$ cover a full solid angle and are divided into 10 and 6 bins, respectively.

The LS and FS models are followed up to $400\,\mathrm{ms}$ and $300\,\mathrm{ms}$ after the core bounce, respectively. The LS model was followed until $300\,\mathrm{ms}$ in Nagakura et al. (2018). The time was not enough to judge whether the shock is revived or not, and hence we extend the time for the LS model. We do not extend the time for the FS model because the FS model clearly fails to explode.

3. Dynamical Features

In this section, we discuss the dynamical features of our two simulations. First, we give an overview of the time evolutions of the shock waves, neutrino luminosities and energies, entropies, and the PNS motion in Section 3.1. Next, we give a detailed analysis of the deformation of the shock and the convection in Section 3.2. In Section 3.3, we investigate the origin of the difference in the shock evolution between the LS and FS models using the timescale ratio. Then, we discuss the influence of the EOSs on the structure of the PNS in Section 3.4.

3.1. Overview of The Dynamics

First of all, we present several quantities to summarize the supernova dynamics in Figure 1: the shock evolutions, the neutrino luminosities, the neutrino mean energies, and the rms energies. Throughout this paper, the shock radius is defined by the outermost radius where the radial velocity gradient becomes negative. The maximum shock radius reaches $\sim 1000\,\mathrm{km}$ for the LS model. The FS model shows a contracting mean shock radius despite the neutrino luminosities and mean energies for the FS model being slightly higher than those for the LS model. Because the maximum and mean shock radii show continuous expansion, we regard the LS model as a successful shock revival model.

**Figure 1.** Time evolution of several key quantities. For each panel, red and blue represent the LS and the FS models, respectively. The upper-left panel shows the shock (solid lines and filled regions) and the gain (dashed lines) radii. The solid lines are the mean shock radii, and the upper and lower edges of the filled regions indicate the maximum and the minimum shock radii. The upper-right, lower-left, and lower-right panels represent the neutrino luminosities, the mean energies, and the rms energies for ${\nu }_{{\rm{e}}}$ (solid line), $\bar{{\nu }_{{\rm{e}}}}$ (dashed line), and ${\nu }_{x}$ (dashed–dotted line), respectively. Both luminosities and energies are measured at $r=500\,\mathrm{km}$ . The upper-right panel is divided into two parts, in order to show the neutronization burst for ${\nu }_{{\rm{e}}}$ and the later light curve at once. The upper and lower halves have different vertical scales each other.
Download figure:
Standard image High-resolution image

To illustrate the dynamical features of our simulations, we show the snapshots of the entropy and the speed $v:= \sqrt{{({v}^{r})}^{2}+{({{rv}}^{\theta })}^{2}}$ at different times for the LS and FS models in Figures 2 and 3, respectively. Note that ${v}^{\theta }$ is the angular velocity, which is usually measured in units of $\mathrm{rad}\,{{\rm{s}}}^{-1}$ . At the onset of the prompt convection ( ${t}_{\mathrm{pb}}=10.5\,\mathrm{ms}$ ), a slightly stronger convective motion is seen in the LS model than in the FS model. As a consequence, the LS model shows a more aspherical, violent shock deformation at ${t}_{\mathrm{pb}}=17.5\,\mathrm{ms}$ , while the sizes of the shock themselves are similar for both models. Then, the shock radii of both models expand gradually as seen in the snapshots at ${t}_{\mathrm{pb}}=60\,\mathrm{ms}$ and ${t}_{\mathrm{pb}}=100\,\mathrm{ms}$ . Again, the shape of the shock in the LS model is more aspherical than that in the FS model. At ${t}_{\mathrm{pb}}=200\,\mathrm{ms}$ and ${t}_{\mathrm{pb}}=300\,\mathrm{ms}$ , the LS model indicates continuous shock expansion. Although the FS model shows vigorous shock deformation at ${t}_{\mathrm{pb}}=200\,\mathrm{ms}$ , the shock contracts after that, as seen in the snapshot at ${t}_{\mathrm{pb}}=300\,\mathrm{ms}$ .

**Figure 2.** Snapshots of the entropy and the velocity of the LS model. For each panel, the left half shows the entropy distribution and the right half displays the speed distribution. The time after the bounce is indicated at the bottom left of each figure.
Download figure:
Standard image High-resolution image

**Figure 3.** Same as Figure 2, except that the snapshots of the FS model are shown.
Download figure:
Standard image High-resolution image

In Figure 4, we show the time evolution of the entropy along the north ( $\theta =0$ ) and the south ( $\theta =\pi$ ) poles for both models up to the end of the simulations. Again, we can see that the LS model shows gradual but continual shock expansion, especially along the north pole, while the FS model does not show shock expansion.

**Figure 4.** Time evolution of the entropy profiles along the north and south poles. The left panel shows the LS model, while the right panel represents the FS model. The horizontal scale is different between the two panels because the simulations are stopped at different times: $400\,\mathrm{ms}$ for the LS model and $300\,\mathrm{ms}$ for the FS model.
Download figure:
Standard image High-resolution image

Figure 5 shows the explosion energies of the LS and FS models. The total energy density is defined by

$\begin{eqnarray}&&{e}_{\mathrm{tot}}:= \displaystyle \frac{1}{2}\rho {v}^{2}+{e}_{\mathrm{th}}+\rho \psi \end{eqnarray} \tag{ 21 }$

in this paper. Here,

$\begin{eqnarray}&&{e}_{\mathrm{th}}:= \rho \displaystyle \sum _{i\in \mathrm{nuclear}\,\mathrm{species}}\displaystyle \frac{3}{2}\displaystyle \frac{{k}_{{\rm{B}}}T}{{m}_{{\rm{u}}}}\displaystyle \frac{{X}_{i}}{{A}_{i}}+{{aT}}^{4}+\left({e}_{{{\rm{e}}}^{-}{{\rm{e}}}^{+}}-\rho {Y}_{{\rm{e}}}\displaystyle \frac{{m}_{{\rm{e}}}{c}^{2}}{{m}_{{\rm{u}}}}\right)\end{eqnarray} \tag{ 22 }$

is the thermal energy defined in Bruenn et al. (2016); X_i, A_i, a, ${e}_{{{\rm{e}}}^{-}{{\rm{e}}}^{+}}$ , and ${m}_{{\rm{e}}}$ are the mass fraction and the mass number of nuclear species i, the radiation constant, the internal energy density for electrons and positrons including rest-mass energy, and the mass of the electron, respectively. For the moment, we explicitly write the light speed c to emphasize that the rest-mass energy of the electron is subtracted. In Figure 5, we consider two kinds of explosion energies: the total positive energy ${E}_{+}$ and the diagnostic explosion energy ${E}_{\mathrm{diag}}$ . The total positive energy is the total energy integrated over the region where the total energy is positive:

$\begin{eqnarray}&&{E}_{+}:= {\int }_{{e}_{\mathrm{tot}}\gt 0}{e}_{\mathrm{tot}}{dV}.\end{eqnarray} \tag{ 23 }$

The diagnostic explosion energy is similar to the total positive energy, but the integrand contains the nuclear binding energy:

$\begin{eqnarray}&&{E}_{\mathrm{diag}}:= {\int }_{{e}_{\mathrm{tot}}\gt 0}({e}_{\mathrm{tot}}+[{e}_{\mathrm{bind}}({}^{56}\mathrm{Fe})-{e}_{\mathrm{bind}}]){dV}.\end{eqnarray} \tag{ 24 }$

The binding energy of iron ${e}_{\mathrm{bind}}({}^{56}\mathrm{Fe})=8.8\,\mathrm{MeV}\times \rho /{m}_{{\rm{u}}}$ is the nuclear binding energy if all nucleons were used to form ${}^{56}\mathrm{Fe}$ . The binding energy ${e}_{\mathrm{bind}}$ is the actual nuclear binding energy evaluated by the internal energy density subtracted by the thermal energy ${e}_{\mathrm{th}}$ .

**Figure 5.** Explosion energies of the LS (red) and FS (blue) models. The solid and dashed lines represent the total positive energy ${E}_{+}$ and the diagnostic explosion energy ${E}_{\mathrm{diag}}$ , respectively. Although the legend for the FS model is indicated, the explosion energy for the FS model is zero.
Download figure:
Standard image High-resolution image

**Figure 5.** Explosion energies of the LS (red) and FS (blue) models. The solid and dashed lines represent the total positive energy ${E}_{+}$ and the diagnostic explosion energy ${E}_{\mathrm{diag}}$ , respectively. Although the legend for the FS model is indicated, the explosion energy for the FS model is zero.
Download figure:
Standard image High-resolution image

Although the diagnostic explosion energy for the LS model is much smaller than ${10}^{51}\,\mathrm{erg}$ at the end of the simulation, it seems to grow at the rate of $\sim 2\times {10}^{50}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ . On the other hand, the explosion energies of the FS model is zero as this model fails to explode.

The central PNSs move due to the recoil of the asymmetric expansion of the matter, but the kick velocities are small. Figure 6 indicates the trajectories of the PNS centers along the symmetry axis for both EOS models. Tracking the trajectory of the PNS center is one of the unique features of our code because we utilize the acceleration frame. As shown in the figure, the offsets are a few kilometers and the velocities are $\sim 10\,\mathrm{km}\,{{\rm{s}}}^{-1}$ . They are not very violent; thus, the differences in the acceleration and the laboratory frames for both models are small. Therefore, we ignore the difference between these frames in this paper. The larger kick is possibly seen after the time we stopped the simulation, or with simulations with other settings such as the progenitor, the rotation, the EOSs, and so on. An example is discussed in Nagakura et al. (2019b).¹⁰

3.2. Deformation of the Shock with Convection

The snapshot at ${t}_{\mathrm{pb}}=17.5\,\mathrm{ms}$ for the LS model (the top-right panel of Figure 2) indicates that the shock in the LS model deforms significantly. To compare the deformation of the shock, we decompose the shock shape into Legendre polynomials. Figure 7 indicates the time evolution of the expansion coefficients normalized by the mean shock radii. The expansion coefficients are defined by

$\begin{eqnarray}&&{a}_{{\ell }}:= \displaystyle \frac{2{\ell }+1}{2}{\int }_{0}^{\pi }{r}_{\mathrm{sh}}(\theta ){P}_{{\ell }}(\cos \theta )\sin \theta d\theta ,\end{eqnarray} \tag{ 25 }$

where ${r}_{\mathrm{sh}}(\theta )$ and ${P}_{{\ell }}(\cos \theta )$ are the shock radius as a function of θ and Legendre polynomial of degree ℓ, respectively. The mean shock radius is nothing but a₀. Here, we show ${\ell }=1$ , 2, and 3 modes in the figure. The mode amplitudes are much larger for the LS model than the FS model.

**Figure 7.** Legendre expansion coefficients of the shock radius. The left and right panels indicate the coefficients for the LS and FS models, respectively. The different colors correspond to different modes: red, blue, and green for the ${\ell }=1$ , 2, and 3 modes, respectively. Each coefficient is normalized by the mean shock radius, i.e., the ${\ell }=0$ mode coefficient a₀.
Download figure:
Standard image High-resolution image

Although oscillations of the shock shape in the large-scale modes at the early postbounce phase are seen in Figure 7, this is not evidence of the strong standing accretion shock instability (SASI) motion. The SASI is the instability of the standing shock by definition. Because the top-right panel of Figure 2 is in the prompt convection phase with an initially expanding shock, this is not the SASI. This deformation of the shock arises from the interaction with the accretion flow outside the shock and the convection inside the shock. When the rising convective parcels reach the shock surface, they are refracted along the shock. The newly accreting flow is involved with this refracting flow, and then it overrides the neighboring shock front. This motion produces violent deformation in the top-right panel of Figure 2.

The deformation of the shock at the standing shock phase, around ${t}_{\mathrm{pb}}\sim 100\,\mathrm{ms}$ , originates from the neutrino-driven convection. Figure 8 shows Foglizzo's χ parameter (Foglizzo et al. 2006); neutrino-driven convection develops when $\chi \gtrsim 3$ , while the SASI develops when $\chi \lesssim 3$ . Here, χ parameter is defined by

$\begin{eqnarray}&&\chi ={\int }_{\langle {r}_{\mathrm{gain}}\rangle }^{\langle {r}_{\mathrm{sh}}\rangle }\displaystyle \frac{{\omega }_{\mathrm{BV},1{\rm{D}}}}{\langle {v}^{r}\rangle }{dr},\end{eqnarray} \tag{ 26 }$

where

$\begin{eqnarray}&&{\omega }_{\mathrm{BV},1{\rm{D}}}=\sqrt{\left|\left(\displaystyle \frac{1}{\langle \gamma \rangle \langle p\rangle }\displaystyle \frac{\partial \langle p\rangle }{\partial r}-\displaystyle \frac{1}{\langle \rho \rangle }\displaystyle \frac{\partial \langle \rho \rangle }{\partial r}\right)\displaystyle \frac{\partial \psi }{\partial r}\right|}\end{eqnarray} \tag{ 27 }$

is the Brunt–Väisälä frequency, and γ is the adiabatic index. The angle brackets represent the angular average of a physical quantity: $\langle \bullet \rangle := \int \bullet d{\rm{\Omega }}/4\pi$ . The original definition in Foglizzo et al. (2006) is meant to be applied under spherical symmetry, and hence Equation (26) is evaluated from the angular-averaged profiles similarly to Fernández et al. (2014) and Iwakami et al. (2014). Note that the χ-parameter in Figure 1 of Nagakura et al. (2018) is defined by the integral from the average gain radius to the maximum shock radius, and hence it is slightly different from the χ parameter in this paper. In Figure 8, we consider the time after ${t}_{\mathrm{pb}}=50\,\mathrm{ms}$ because the gain region appears only after that time. Figure 8 clearly shows that $\chi \gt 3$ for both LS and FS models before ${t}_{\mathrm{pb}}\sim 200\,\mathrm{ms}$ , and hence, neutrino-driven convection develops first in both models. After ${t}_{\mathrm{pb}}\sim 200\,\mathrm{ms}$ , the χ parameter in the FS model becomes smaller than 3. Hence, the SASI may develop in the FS model, but the highly developed turbulence makes it difficult to distinguish the class of fluid instability, the convection, or SASI, at that time.

**Figure 8.** Time evolution of Foglizzo's χ parameter. The red and blue lines represent the LS and FS models, respectively. The horizontal dotted line indicates $\chi =3$ , where the class of the fluid instability bifurcates.
Download figure:
Standard image High-resolution image

**Figure 8.** Time evolution of Foglizzo's χ parameter. The red and blue lines represent the LS and FS models, respectively. The horizontal dotted line indicates $\chi =3$ , where the class of the fluid instability bifurcates.
Download figure:
Standard image High-resolution image

Compared to Müller et al. (2012) and Takiwaki et al. (2012, 2014), the ${\ell }=1$ mode amplitudes of the shock in Figure 7 are $\sim 10\mbox{--}{10}^{2}$ times larger. We impose the initial perturbations on our simulations in a different way from their simulations. The radial velocity perturbation with an amplitude of $0.1 \%$ is imposed in the preshock region at ${t}_{\mathrm{pb}}\sim 0.6\,\mathrm{ms}$ in our simulations, while the velocity perturbation with an amplitude of 1% is employed in the postshock region at ${t}_{\mathrm{pb}}=10\,\mathrm{ms}$ in Takiwaki et al. (2014), for example. According to Abdikamalov et al. (2016), the energy of the preshock perturbation (squared velocity) is amplified by a factor of ∼2 when passing the shock front, and it plays a minor role here. Owing to the prompt convection, the perturbation initially given at ${ \mathcal O }(0.1)\mathrm{ms}$ grows much larger than 10× at ${t}_{\mathrm{pb}}=10\,\mathrm{ms}$ as estimated in the following. The convective growth rate, or Brunt–Väisälä frequency, is approximately (see Figure 2 in Nagakura et al. 2018)

$\begin{eqnarray}\begin{array}{rcl}{\omega }_{\mathrm{BV}} & \sim & \sqrt{\displaystyle \frac{{{GM}}_{\mathrm{PNS}}}{{r}^{2}}\displaystyle \frac{\partial s/{k}_{{\rm{B}}}}{\partial r}}\\ & = & 2.4\,{\mathrm{ms}}^{-1}\\ & & \times {\left(\displaystyle \frac{{M}_{\mathrm{PNS}}}{0.7{M}_{\odot }}\right)}^{1/2}{\left(\displaystyle \frac{r}{40\mathrm{km}}\right)}^{-1}{\left(\displaystyle \frac{\partial (s/{k}_{{\rm{B}}})/\partial r}{{10}^{-6}{\mathrm{cm}}^{-1}}\right)}^{1/2},\end{array}\end{eqnarray} \tag{ 28 }$

where the gravitational acceleration is approximated by

$\begin{eqnarray}&&\displaystyle \frac{\partial \psi }{\partial r}\sim \displaystyle \frac{{{GM}}_{\mathrm{PNS}}}{{r}^{2}},\end{eqnarray} \tag{ 29 }$

and ${M}_{\mathrm{PNS}}$ is the PNS mass. As long as the perturbation is linear, the initial perturbation at ${ \mathcal O }(0.1)\,\mathrm{ms}$ grows $\exp (10\,\mathrm{ms}\,\times {\omega }_{\mathrm{BV}})\sim 3\,\times \,{10}^{10}$ times larger at $10\,\mathrm{ms}$ . Such a large amplification breaks the linear approximation, and the actual amplification should be smaller than this estimation. It is still certain that the initial $0.1 \%$ velocity perturbation at 0.6 ms like ours grows to exceed 1% velocity perturbation in Takiwaki et al. (2014) at $10\,\mathrm{ms}$ . Because the initial perturbation is effectively large in this work, the prompt convection and hence the shock deformation are much larger than other works. Although the central iron core is convectively stable, a possible slight perturbation imparted by the convective silicon layer can be amplified during core collapse (Hanawa & Matsumoto 2000). Hence, the early initial perturbation employed in this paper is possible.

3.3. Timescale Ratio and the Effects of the Nuclear Composition of the EOSs

Next, we evaluate the timescale ratio ${\tau }_{\mathrm{adv}}/{\tau }_{\mathrm{heat}}$ to see how close our two simulations are to success. If this ratio exceeds unity, successful shock revival occurs. Here, ${\tau }_{\mathrm{adv}}:= {M}_{\mathrm{gain}}/{\dot{M}}_{\mathrm{sh}}$ and ${\tau }_{\mathrm{heat}}:= | {E}_{\mathrm{gain}}| /{Q}_{\mathrm{gain}}$ are the advection timescale through the gain layer and the heating timescale, respectively, where ${M}_{\mathrm{gain}}$ , ${\dot{M}}_{\mathrm{sh}}$ , ${E}_{\mathrm{gain}}$ , and ${Q}_{\mathrm{gain}}$ are the gain mass, the mass accretion rate at the shock, the total energy in the gain layer, and the net heating rate in the gain layer, respectively. The mass accretion rate is defined as ${\dot{M}}_{\mathrm{sh}}:= -{\int }_{\mathrm{shock}}4\pi \rho {v}^{r}{r}^{2}d{\rm{\Omega }}$ . In this paper, we define ${E}_{\mathrm{gain}}$ as

$\begin{eqnarray}&&{E}_{\mathrm{gain}}:= {\int }_{\mathrm{gain}}{e}_{\mathrm{tot}}{dV},\end{eqnarray} \tag{ 30 }$

where ${e}_{\mathrm{tot}}$ is defined in Equation (21).

The definition of the timescale ratio in this paper is improved from that in the previous paper (Nagakura et al. 2018). They are different in two respects: the definition of the total energy and the domain of the integral are different. In the previous paper, the total energy includes the nuclear binding energy, whose offset we are free to choose. The total energy defined by Equation (21) with the thermal energy defined by Equation (30) is free from such ambiguity. The integral in the previous paper is the radial integral of the angular-averaged quantities from the mean gain radius to the minimum shock radius. On the other hand, the angular dependence of the gain radii and shock radii is taken into account in this paper.

The timescale ratio in the LS model exceeds unity for most of the time after $\sim 70\,\mathrm{ms}$ as indicated in Figure 9. On the other hand, the FS model shows a ratio smaller than unity for almost the entirety of the simulation period. This implies that the LS model is favorable for a successful explosion. It is worth noting that the timescale ratio of the FS model exceeds unity in Nagakura et al. (2018) due to the inappropriate treatment of thermal energy. Thanks to Equation (30), more reliable evaluations of the timescale ratio are obtained.

**Figure 9.** Evolution of the timescale ratio ${\tau }_{\mathrm{adv}}/{\tau }_{\mathrm{heat}}$ and related quantities. For all panels, the red and blue lines correspond to the LS and FS models, respectively. The top-left panel indicates the timescale ratio ${\tau }_{\mathrm{adv}}/{\tau }_{\mathrm{heat}}$ with the horizontal dotted line indicating unity. The top-right panel shows the advection timescale (solid lines) and the heating timescale (dashed lines) themselves. The middle-left, middle-right, bottom-left, and bottom-right panels display the gain mass, the mass accretion rate, the total energy in the gain region, and the net heating rate in the gain region, respectively.
Download figure:
Standard image High-resolution image

A closer look reveals to us that the heating timescales are similar for both models, while the advection timescales are much larger for the LS model. For the advection timescales, the mass accretion rate is not that different between the two models, but the gain mass is much larger for the LS model. For the heating timescales, the total energy in the gain region ${E}_{\mathrm{gain}}$ is slightly lower for the LS model, but this is compensated by the slightly higher net heating rate ${Q}_{\mathrm{gain}}$ .

The reason for the difference in the gain mass is the strength of the turbulent motion. The right-side parts of each panel of Figures 2 and 3 show the speed of matter $v=\sqrt{{({v}^{r})}^{2}+{({{rv}}^{\theta })}^{2}}$ . From the snapshots at ${t}_{\mathrm{pb}}=10.5\,\mathrm{ms}$ , the speed of matter at the prompt convection in the LS model is more vigorous than that in the FS model. This stronger prompt convection in the LS model gives stronger seed perturbations of the later neutrino-driven convection. This develops stronger turbulence. This mechanism was also suggested in the previous paper (Nagakura et al. 2018), but we present a more detailed description here.

The origin of the stronger prompt convection is the difference in the radial gradient of the entropy profiles. The convectively unstable regions at ${t}_{\mathrm{pb}}=10.5\,\mathrm{ms}$ in Figures 2 and 3 are up to a few tens of kilometers. Figure 10 shows the entropy profiles at the time when the shock waves have swept the convectively unstable region. Due to the weakening of the propagating shock, the radial gradients of the entropy are negative in $35\,\mathrm{km}\lesssim r\lesssim 50\,\mathrm{km}$ for both models. The gradient is steeper for the LS model than for the FS model. Note that the value of the entropy gradient discussed in Section 3.2 is estimated from this figure.

**Figure 10.** Radial entropy profiles at ${t}_{\mathrm{pb}}\simeq 1.7\,\mathrm{ms}$ , when the shock waves have swept the convectively unstable region. The red and blue lines are for the LS and FS models, respectively.
Download figure:
Standard image High-resolution image

**Figure 10.** Radial entropy profiles at ${t}_{\mathrm{pb}}\simeq 1.7\,\mathrm{ms}$ , when the shock waves have swept the convectively unstable region. The red and blue lines are for the LS and FS models, respectively.
Download figure:
Standard image High-resolution image

Because an important player in the shock weakening is the photodissociation of heavy nuclei, we compare the nuclear composition between the LS and FS models in Figure 11. The time when the shock sweeps the region $35\,\mathrm{km}\lesssim r\lesssim 50\,\mathrm{km}$ is merely $0.9\,\mathrm{ms}\lesssim {t}_{\mathrm{pb}}\lesssim 1.7\,\mathrm{ms}$ . Figure 11 shows the average mass number $\langle A{\rangle }_{\mathrm{sh}}$ and the average mass fraction $\langle {X}_{A}{\rangle }_{\mathrm{sh}}$ of the matter slightly outside the shock during this period for both the LS and FS models. Also, the mass fraction of alpha particles $\langle {X}_{\alpha }{\rangle }_{\mathrm{sh}}$ is shown. The angular-averaged quantity of the accreting matter $\langle \bullet {\rangle }_{\mathrm{sh}}$ is defined as ${\int }_{{r}_{\mathrm{sh}}}\bullet d{\rm{\Omega }}/{\int }_{{r}_{\mathrm{sh}}}d{\rm{\Omega }}$ , and the domain of the integral is the surface slightly outside the shock. Here, $\langle A{\rangle }_{\mathrm{sh}}$ for the LS EOS is the mass number of a representative heavy nucleus, while that for the FS EOS is the average mass number of the nuclei in the NSE. For both EOSs, $\langle {X}_{A}{\rangle }_{\mathrm{sh}}$ is the mass fraction of the heavy nuclei.

**Figure 11.** Nuclear composition of accreting matter just above the shock for the LS and FS models during the period relevant to the growth rate of the prompt convection. The top panel shows the mass number $\langle A{\rangle }_{\mathrm{sh}}$ . The middle panel displays the mass fraction $\langle {X}_{A,\alpha }{\rangle }_{\mathrm{sh}}$ of the representative heavy nucleus (solid lines) and the alpha particle (dashed lines). The horizontal dotted line indicates unity. The bottom panel indicates the estimated nuclear binding energy per nucleon of the accreting matter. For all panels, the red and blue lines correspond to the LS and FS models, respectively.
Download figure:
Standard image High-resolution image

**Figure 11.** Nuclear composition of accreting matter just above the shock for the LS and FS models during the period relevant to the growth rate of the prompt convection. The top panel shows the mass number $\langle A{\rangle }_{\mathrm{sh}}$ . The middle panel displays the mass fraction $\langle {X}_{A,\alpha }{\rangle }_{\mathrm{sh}}$ of the representative heavy nucleus (solid lines) and the alpha particle (dashed lines). The horizontal dotted line indicates unity. The bottom panel indicates the estimated nuclear binding energy per nucleon of the accreting matter. For all panels, the red and blue lines correspond to the LS and FS models, respectively.
Download figure:
Standard image High-resolution image

The difference in the nuclear composition seems to play a key role in creating the difference in the entropy gradient. The heavy nuclei and alpha particles are finally dissolved into nucleons and consume the shock energy. The heavy nuclei and alpha particles in the LS model are more and less abundant than in the FS model, respectively. Hence, which model consumes more energy is not apparent. In order to estimate the loss of shock energy, Figure 11 also shows the total nuclear binding energy per nucleon of the accreting matter. Here, the binding energy of the heavy nuclei per nucleon is approximately set to $8.8\,\mathrm{MeV}$ , the value for ${}^{56}\mathrm{Fe}$ , because the binding energy of the nuclei around ${}^{56}\mathrm{Fe}$ is insensitive to the mass number. The binding energy of the alpha particle per nucleon is $7.1\,\mathrm{MeV}$ . Therefore, the total binding energy per nucleon estimated here is $\langle {X}_{A}{\rangle }_{\mathrm{sh}}\times 8.8\,\mathrm{MeV}+\langle {X}_{\alpha }{\rangle }_{\mathrm{sh}}\times 7.1\,\mathrm{MeV}$ . The figure indicates that the total binding energy is slightly higher for the LS model. It implies that more shock energy is lost due to photodissociation during that period, and the shock is weakened more rapidly in the LS model. This seems to produce the steeper negative radial gradient of the entropy to drive the prompt convection: the Brunt–Väisälä frequency, i.e., the convective growth rate, gets larger in the LS model. This appears to be the origin of the strong prompt convection in the LS model.

Although the mass accretion rate at the shock ${\dot{M}}_{\mathrm{sh}}$ shown in Figure 9 is smooth, that inside the shock is modulated by the shock motion. Precisely, the mass accretion rate at the shock is measured slightly outside the shock, where the flow is smooth. Figure 12 shows the mass accretion rate $\dot{m}$ in the LS model as a function of time and radius. Here, the mass accretion rate is measured at a constant radius: $\dot{m}:= -4\pi {r}^{2}\langle \rho \rangle \langle {v}^{r}\rangle$ , where $\langle \bullet \rangle := {\int }_{r}\bullet d{\rm{\Omega }}/{\int }_{r}d{\rm{\Omega }}$ means the angular average. The shock significantly recedes at ${t}_{\mathrm{pb}}\sim 250\,\mathrm{ms}$ . This results in significant mass accretion inside the shock at that time in Figure 12. Owing to the released gravitational energy of this strong accretion, the ${\nu }_{{\rm{e}}}$ luminosity increases momentarily as seen in Figure 1.

**Figure 12.** Mass accretion rate $\dot{m}$ in the LS model as a function of time and radius. This is evaluated from the angle-averaged density and radial velocity.
Download figure:
Standard image High-resolution image

3.4. Structure of the PNS and the Influence of the EOSs

Figure 13 displays the evolutions of the neutrinosphere radii and temperatures. The neutrinosphere radius is defined as the radius where the angle-averaged total (absorption+scattering) optical depth for the mean-energy neutrino is 2/3, and the neutrinosphere temperature is the angle-averaged temperature at the neutrinosphere radius. Here, we attempt to illustrate the similarity of the neutrino luminosities and differences in the mean and rms energies of neutrinos between the LS and FS models in Figure 1. The neutrinosphere temperature is a useful indicator of the neutrino mean energies. The neutrinosphere radius is slightly smaller for the FS model, but the neutrinosphere temperature is higher for the FS model. Because they compensate each other, the neutrino luminosities for both models are similar. The difference in the temperature explains the difference in the mean and rms energies between the two models.

Although the LS EOS is known to be a soft EOS, the outer part of the PNS in the LS model is less compact than the FS model. One may expect a smaller PNS radius, a smaller neutrinosphere radius, and hence a higher neutrinosphere temperature with the LS EOS compared with those with the FS EOS, which is a stiffer EOS. Figure 13, however, shows the opposite result. To understand this, we show the time evolutions of the central density in Figure 14. The structure of the PNS, namely the radii, the enclosed masses, and the compactness at the angular-averaged densities of 10¹¹, 10¹², 10¹³, and ${10}^{14}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ , is also shown. Here, the compactness is defined as

$\begin{eqnarray}&&{\xi }_{\rho }=\displaystyle \frac{{M}_{\rho }/{M}_{\odot }}{{r}_{\rho }/10\,\mathrm{km}},\end{eqnarray} \tag{ 31 }$

where M_ρ and r_ρ are the enclosed mass and radius at the angle-averaged density ρ. The central densities for the LS model are larger than those for the FS model, as expected. The radii and the enclosed masses at $\rho ={10}^{14}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ are similar and larger for the LS model, respectively. Thus, the LS model shows higher compactness at this density. This is again consistent with the stiffness of the EOS. On the other hand, the radii and the enclosed masses at $\rho ={10}^{13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ are similar and smaller for the LS model, respectively. The resultant compactness for the LS model is lower than that for the FS model. At both $\rho ={10}^{12}$ and ${10}^{11}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ , the radii, the enclosed mass, and the compactness for the LS model are larger, smaller, and lower, respectively. The compactness for the LS model is higher only near the nuclear densities.

**Figure 14.** Evolution of the central density and the PNS compactness for the LS (red) and FS (blue) models. The upper-left panel shows the central density evolution. The upper-right and lower-left panels indicate the radii and enclosed masses where the densities are 10¹¹ (solid), 10¹² (dashed), 10¹³ (dashed–dotted), and 10¹⁴ (dotted) ${\rm{g}}\,{\mathrm{cm}}^{-3}$ . The lower-right panel indicates the compactness of each density defined by Equation (31). Each quantity is calculated from the 1D averaged radial profiles.
Download figure:
Standard image High-resolution image

The higher compactness in the FS model originates from the lower stiffness of the EOS at the densities around ${10}^{13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}\lesssim \rho \lt 5\times {10}^{13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ . We show the pressure versus the density in the top panel of Figure 15. The pressures for both EOSs are almost identical at densities lower than $\sim {10}^{13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ and higher for the FS EOS at densities higher than $\sim {10}^{14}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ . At around $\sim 4\times {10}^{13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ , the pressure is lower for the FS EOS. We show an effective adiabatic index $\gamma := \partial \mathrm{ln}P/\partial \mathrm{ln}\rho$ in the middle panel, although the "stiffness" of the EOS is determined by the adiabatic index ${(\partial \mathrm{ln}P/\partial \mathrm{ln}\rho )}_{s}$ . Here, the entropy is not necessarily constant for the differentiation, and hence what is shown in the middle panel is nothing but the slope of the top panel. We use this effective γ as an indicator of the stiffness. For densities higher than $\sim 5\times {10}^{13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ , the effective γ is higher for the FS EOS, indicating that the FS EOS is stiff. However, for density between $\sim {10}^{13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ and $\sim 5\times {10}^{13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ , the effective γ is lower for the FS EOS. This softness leads to the more compact structure of the PNS at the region where the density is lower than ${10}^{13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ for the FS model as shown in Figure 14.

**Figure 15.** Angular-averaged thermal properties and compositions vs. the density at ${t}_{\mathrm{pb}}\simeq 200\,\mathrm{ms}$ for the LS (dashed) and FS (solid) models. The top panel shows the relation between the density and the pressure, and the middle panel indicates the effective adiabatic index γ (see the text). The bottom panel displays the corresponding composition: green for the neutron, magenta for the proton, cyan for the deuteron, yellow for the triton, red for the light nuclei, and black for the heavy nuclei. The light nuclei here are the nuclei with $Z\lt 6$ other than deuteron, triton, helion, and the alpha particle; the heavy nuclei are nuclei with $Z\geqslant 6$ . For the FS EOS, the mass fractions of the helion and alpha particle are less than 10⁻².
Download figure:
Standard image High-resolution image

**Figure 15.** Angular-averaged thermal properties and compositions vs. the density at ${t}_{\mathrm{pb}}\simeq 200\,\mathrm{ms}$ for the LS (dashed) and FS (solid) models. The top panel shows the relation between the density and the pressure, and the middle panel indicates the effective adiabatic index γ (see the text). The bottom panel displays the corresponding composition: green for the neutron, magenta for the proton, cyan for the deuteron, yellow for the triton, red for the light nuclei, and black for the heavy nuclei. The light nuclei here are the nuclei with $Z\lt 6$ other than deuteron, triton, helion, and the alpha particle; the heavy nuclei are nuclei with $Z\geqslant 6$ . For the FS EOS, the mass fractions of the helion and alpha particle are less than 10⁻².
Download figure:
Standard image High-resolution image

The mass accretion rate might influence the PNS compactness, but it is minor. If the mass accretion rate is low, the thermal energy provided to the PNS is low, and hence the PNS becomes compact (Nagakura et al. 2019a). Although the mass accretion rate of the FS model is very similar to or slightly higher than that of the LS model until $\sim 200\,\mathrm{ms}$ after the bounce as shown in Figure 9, the PNS in the FS model is more compact. This implies that the effects of mass accretion rate is minor compared to that of the nuclear composition.

The low stiffness at the densities discussed above for the FS EOS originates from the difference in the composition. In the bottom panel of Figure 15, we show the mass fractions of nuclei that are larger than 1%. In the high-density region $\rho \gt {10}^{14}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ , only neutrons and protons are the constituents of matter. In the subnuclear density region, light and heavy nuclei start to appear in the FS model, while no light and heavy nuclei appear for the LS model. This is simply because the NSE is considered in the FS EOS and light and heavy nuclei can appear, while only a single representative heavy nucleus and alpha particle are considered in the LS EOS. When searching for the free energy minimum to construct the EOS, there are more degrees of freedom for the FS EOS than those for the LS EOS due to the degrees of freedom in the composition. Therefore, the achieved free energy is smaller for the FS EOS, and hence the pressure and the effective γ are smaller.

A caveat is that the updated version of the FS EOS (Furusawa et al. 2017a) might be stiffer than the FS EOS employed in this paper, though it is probably still softer than the LS EOS. The temperature in the region where heavy nuclei appear in the FS model is $\sim 20\,\mathrm{MeV}$ . The heavy nuclei are dissolved into free nucleons with temperatures above $18\,\mathrm{MeV}$ in the updated FS EOS (Furusawa et al. 2017a). This effect reduces the degree of freedom of the heavy nuclei, and hence, the updated FS EOS might be stiffer than the FS EOS in this paper. Although this melting of the heavy nuclei is not included in our current simulation, the degree of freedom of the composition of the light nuclei is probably still enough to make the updated FS EOS softer than the LS EOS at subnuclear densities. It is not only the stiffness but also the opacity that might be affected by the nuclear composition in the EOS. One might think that the difference in the neutrino mean and rms energies between the LS and FS models shown in Figure 1 is attributable to the different ${Y}_{{\rm{p}}}$ because the reaction rate of the electron capture on protons is reduced with small ${Y}_{{\rm{p}}}$ . The opacity of the electron capture reaction in our simulations, however, does not fully take the composition shown in Figure 15 into account. Although the reduction of the proton mass fraction is due to the presence of the light nuclei, the weak interaction rate is evaluated as if the light nuclei were dissolved into free nucleons as described in Section 2. Hence, the proton mass fraction that was used to evaluate the electron capture rate in the FS model is not as low as ${Y}_{{\rm{p}}}$ in Figure 15 actually.

Even if we use ${Y}_{{\rm{p}}}$ in Figure 15 to evaluate the electron capture rate of free protons, the effects on the mean and rms energies are still limited. The neutrinosphere for the mean-energy neutrinos is located at $\sim {10}^{11}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ . At this density, the nuclei are almost dissolved, and the ${Y}_{{\rm{p}}}$ in the FS model is close to that in the LS model in Figure 15. The mean and rms energies are mainly determined by the neutrinosphere as discussed so far, and hence, the mean-energy neutrinos are hardly affected by the presence of light nuclei in this case. Besides, the mean and rms energies are also affected by the low-density region. In the low-density region, the nuclei are completely dissolved, and hence ${Y}_{{\rm{p}}}$ is similar for the LS and FS models. Because the optical depth is determined by the integral of the opacity outside the given radius, the location of the neutrinosphere is almost determined by the opacity in the region that has nothing to do with the ${Y}_{{\rm{p}}}$ reduction in Figure 15. We also note that the mean and rms energies are modulated during the propagation from the neutrinosphere to the radius where the neutrino quantities are measured. This is because the neutrinos with energies higher than the mean still interact with matter outside the mean-energy neutrinosphere, although the mean-energy neutrinos are decoupled from matter. This modulation is also determined by the opacity in the low-density region.

In general, the existence of light and heavy nuclei can itself influence the opacity (e.g., Sumiyoshi & Röpke 2008; Fischer et al. 2016; Nagakura et al. 2019a). The cross sections for the weak interactions of light and heavy nuclei are different from those of the nucleons (Fischer et al. 2016). If we consider the opacity of the weak interactions to be consistent with the composition shown in Figure 15, the mean and rms energies of neutrinos might be modified as discussed in Nagakura et al. (2019a). This effect is, however, beyond the scope of this paper.

4. Neutrino Distributions

In this section, we discuss the physical quantities related to the neutrino angular distributions. The neutrino flux, the first angular moment, is discussed in Section 4.1. The Eddington tensor, the second angular moment, is discussed in Section 4.2. Because we solve the Boltzmann equations for neutrinos, we can provide the information of the angular distributions.

4.1. Neutrino Flux

Figure 16 shows the number density and direction of the number flux of ${\nu }_{{\rm{e}}}$ in both the laboratory and fluid rest frames at different scales at ${t}_{\mathrm{pb}}=10.5\,\mathrm{ms}$ . For the laboratory frame, the direction of fluid velocity is also shown. First, let us focus on the zoom-in snapshot (the left panel).

At the central $r\lesssim 20\,\mathrm{km}$ region, the neutrino fluxes have various directions in the laboratory frame, while they are almost radially directed in the fluid rest frame. Because the central region is opaque, the neutrinos move in tandem with matter, and hence their directions in the laboratory frame are not uniform. On the other hand, the radial flux in the fluid rest frame comes from the diffusion of the neutrinos. The lateral motion seen in the laboratory frame is purely a consequence of neutrinos and matter comoving. Note that due to the diffusion flux, the directions of the fluid velocity and the neutrino flux are slightly different. Nagakura et al. (2019c) discussed that a precise evaluation of the flux in the employed code is difficult. The fixed-length arrows in Figure 16 do not show the magnitude of the diffusion flux, and it is small at the very central region actually. Therefore, the evaluation of such small flux is perhaps not precise.

At the region where $20\,\mathrm{km}\lesssim r\lesssim 50\,\mathrm{km}$ , neutrinos also move outward, but the lateral component alternately changes between clockwise (reddish arrows) and counterclockwise (bluish arrows) in the fluid rest frame. This is because of the convection and diffusion. The time considered here is the very early stage of the prompt convection. The region $20\,\mathrm{km}\lesssim r\lesssim 50\,\mathrm{km}$ is convectively unstable. The radial gradient of the ${\nu }_{{\rm{e}}}$ number density is negative on average there. Due to the prompt convection, some fluid parcels rise and some sink in this region. The rising fluid parcels hence have a large ${\nu }_{{\rm{e}}}$ number, whereas the sinking parcels have a small ${\nu }_{{\rm{e}}}$ number. Thus, the finger-like pattern in the number density shown in the figure develops. Due to the difference in number density, diffusion flux emerges and neutrinos flow from rising parcels to sinking parcels in the fluid rest frame as clearly seen in the figure. In the laboratory frame, the neutrino flux is also affected by the matter velocity.

Although the flux in the inner region is determined by the diffusion and matter velocity, the neutrinos simply flow almost radially in the outer region. This is shown in the right panel of Figure 16. The zoom-in region shown in the left panel is indicated by the orange rectangle in the right panel. Because the outer region of the orange rectangle is relatively optically thin, the neutrinos freely stream, and the flux is radially directed.

Figure 17 shows the ${\nu }_{{\rm{e}}}$ number densities and fluxes at different times. Basically, the behavior of the neutrino flux is similar to the snapshot at ${t}_{\mathrm{pb}}=10.5\,\mathrm{ms}$ : the direction is determined by the diffusion and matter velocity at the inner region and by the free-streaming at the outer region. However, at later times, the prompt convection diminishes and the region just outside the PNS is convectively stable. As a consequence, the neutrino number density is almost isotropic, and the diffusion flux is directed almost radially in the fluid rest frame. The flux in the laboratory frame is also almost radially directed but dragged by the matter velocity to some extent especially inside the PNS. The matter velocity inside the PNS originates from the PNS convection driven by the ${Y}_{{\rm{e}}}$ gradient.

**Figure 17.** Same as Figure 16, except for the time of snapshots: ${t}_{\mathrm{pb}}=17.5\,\mathrm{ms}$ (upper left), $60\,\mathrm{ms}$ (upper right), 100 ms (lower left), and 200 ms (lower right).
Download figure:
Standard image High-resolution image

The diffusion flux of $\bar{{\nu }_{{\rm{e}}}}$ is different from that of ${\nu }_{{\rm{e}}}$ due to the positive radial gradient of the $\bar{{\nu }_{{\rm{e}}}}$ number density. Figure 18 indicates the number density and direction of the number flux of $\bar{{\nu }_{{\rm{e}}}}$ . The flux in the upper-left panel shows a different pattern from the left panel of Figure 16. Although it is similar in that the lateral component changes alternately in the region where $20\,\mathrm{km}\lesssim r\lesssim 50\,\mathrm{km}$ , the clockwise (reddish) and counterclockwise (bluish) pattern is opposite to the pattern seen in the ${\nu }_{{\rm{e}}}$ flux. The flux in the fluid rest frame is directed inward around the center. These differences originate from the positive radial gradient of number density there. Due to the positive radial gradient, the sinking fluid parcels have a large $\bar{{\nu }_{{\rm{e}}}}$ number, while the rising parcels have a small $\bar{{\nu }_{{\rm{e}}}}$ number as shown in Figure 18. The lateral component of the $\bar{{\nu }_{{\rm{e}}}}$ diffusion flux hence shows the opposite sign to the ${\nu }_{{\rm{e}}}$ diffusion flux.

The positive radial gradient of $\bar{{\nu }_{{\rm{e}}}}$ comes from the degeneracy of ${\nu }_{{\rm{e}}}$ . The electron-type neutrinos are in β equilibrium and degenerate at the center, and hence, their antiparticles are suppressed. Around the center, the matter density decreases, and the chemical potential of ${\nu }_{{\rm{e}}}$ also decreases. The temperature increases around the center owing to the neutrino diffusion (Burrows et al. 1981; Bethe 1990). The decreasing chemical potential and increasing temperature make electron-type neutrinos nondegenerate, and hence electron-type antineutrinos start to appear. The degeneracy of ${\nu }_{{\rm{e}}}$ decreases with the radius, then the number density of $\bar{{\nu }_{{\rm{e}}}}$ increases with the radius.

Except for the different pattern in the $\bar{{\nu }_{{\rm{e}}}}$ number density and the diffusion flux, the behavior of the $\bar{{\nu }_{{\rm{e}}}}$ number flux is similar to that of the ${\nu }_{{\rm{e}}}$ number flux. In the inner region, the flux is determined by the diffusion and the matter velocity. In the outer region, $\bar{{\nu }_{{\rm{e}}}}$ streams freely. These behaviors are also seen in the snapshots at different times.

Figure 19 displays the number density and direction of the number flux of ${\nu }_{x}$ . Again, how the number flux of ${\nu }_{x}$ is determined is similar to that of ${\nu }_{{\rm{e}}}$ and $\bar{{\nu }_{{\rm{e}}}}$ : the flux is determined by the diffusion and motion of matter in the optically thick region and the free-streaming in the optically thin region. The distribution of the flux of ${\nu }_{x}$ itself is similar to that of $\bar{{\nu }_{{\rm{e}}}}$ . Because the chemical potential is zero, the distribution of ${\nu }_{x}$ in the optically thick region is determined by the temperature. Therefore, the peak in the ${\nu }_{x}$ number density lies not in the center but around the center, and the radial gradient of the number density is positive around the center. Thus, the deformation by convection results in the diffusion flux from the sinking fluid parcels to the rising parcels. However, compared with $\bar{{\nu }_{{\rm{e}}}}$ , the peak is located slightly closer to the center due to the zero chemical potential, and hence, the resultant pattern in the distribution and diffusion flux is also slightly different.

Finally, let us compare the neutrino fluxes of the LS and FS models briefly. Basically, the overall behavior is almost the same for both models, though we show no figures: the direction of the flux is determined by the diffusion and matter velocity at the inner region and by the free-streaming at the outer region. Due to the difference in EOSs, the distributions of the background fluid are different as discussed in Section 3, and hence the resultant flux distributions are also different. Here we only focus on the early stage, when the turbulence is not fully developed. This is because the comparison at the late stage is not meaningful: the chaotic nature of the completely developed turbulence makes the flux patterns too different to be compared.

Figure 20 shows the radial profile of the angular rms fluctuation of ${Y}_{{\rm{e}}}$ in order to show how a large area is mixed by convection. Here, the angular rms fluctuation is defined as $\sqrt{\langle {({Y}_{{\rm{e}}}-\langle {Y}_{{\rm{e}}}\rangle )}^{2}\rangle }$ with $\langle \bullet \rangle := {\int }_{r}\bullet d{\rm{\Omega }}/{\int }_{r}d{\rm{\Omega }}$ , where r is given as a constant radius. The central region and the unshocked accretion flow are spherically symmetric, and hence the rms fluctuation is zero. On the other hand, the region $20\,\mathrm{km}\lesssim r\lesssim 50\,\mathrm{km}$ is mixed by the convection, and hence the rms fluctuation becomes large. The size of the region where the rms fluctuation is large in the FS model is smaller than that in the LS model because the prompt convection is stronger in the LS model.

As a result, the lateral flux by diffusion emerges from a larger region for the LS model as shown in Figure 21. We choose $r=20\,\mathrm{km}$ and $r=40\,\mathrm{km}$ for the radii where the angular profiles of the lateral flux are shown. This is because these radii are appropriate to see the difference in the strength of the prompt convection: the convective parcels in the FS model are almost confined between these radii, while those in the LS model overshoot these radii. As expected, the diffusion flux in the FS model is much smaller than that in the LS model in those regions. In the region where convection develops for both models, e.g., at the radius $r=30\,\mathrm{km}$ , the orders of magnitude of the lateral diffusion fluxes are similar, though we do not show in the figure.

4.2. Eddington Tensor

Let us now look into the Eddington tensor, the second moment of the distribution function. Here, we obey the definition of the Eddington tensor presented in Shibata et al. (2011),¹¹ considering a 3 + 1 decomposition of the spacetime. First, the second moment of the distribution function is defined as

$\begin{eqnarray}\begin{array}{rcl}{M}^{\alpha \beta }(\epsilon ) & := & \displaystyle \int f\delta \left(\displaystyle \frac{{\epsilon }^{3}}{3}-\displaystyle \frac{{\epsilon }^{{\prime} 3}}{3}\right){p}^{{\prime} \alpha }{p}^{{\prime} \beta }{{dV}}_{{p}^{{\prime} }}\\ & = & \displaystyle \frac{1}{\epsilon }\displaystyle \int {{fp}}^{{\prime} \alpha }{p}^{{\prime} \beta }d{{\rm{\Omega }}}_{p}^{{\prime} }\left|{}_{{\epsilon }^{{\prime} }=\epsilon }\right.,\end{array}\end{eqnarray} \tag{ 32 }$

where ${\epsilon }^{{\prime} }:= -{p}^{{\prime} }\cdot {e}_{(0)}$ is the neutrino energy measured in the fluid rest frame (Thorne 1981); ${{dV}}_{{p}^{{\prime} }}$ is the invariant volume element of the momentum space. The integral in the second line is evaluated with the condition that the energy measured in the fluid rest frame is the constant . Next, we define the spatial–spatial and temporal–temporal projections of the second moment, which are nothing but the stress tensor and the energy density, respectively, as

$\begin{eqnarray}&&{P}^{{ij}}(\epsilon ):= {\gamma }^{i}{}_{\alpha }{\gamma }^{j}{}_{\beta }{M}^{\alpha \beta }(\epsilon )\end{eqnarray} \tag{ 33 }$

and

$\begin{eqnarray}&&E(\epsilon ):= {n}_{\alpha }{n}_{\beta }{M}^{\alpha \beta }(\epsilon ),\end{eqnarray} \tag{ 34 }$

where ${\gamma }_{\alpha }^{i}$ and n_α are the spatial projection tensor and the normal vector to the spatial hypersurface, respectively. Finally, the Eddington tensor ${k}^{{ij}}(\epsilon )$ is defined by

$\begin{eqnarray}&&{k}^{{ij}}(\epsilon ):= \displaystyle \frac{{P}^{{ij}}(\epsilon )}{E(\epsilon )}.\end{eqnarray} \tag{ 35 }$

For the later convenience, we call this ${k}^{{ij}}(\epsilon )$ the "Boltzmann–Eddington tensor."

We compare the Boltzmann–Eddington tensor and the Eddington tensor calculated from the M1-closure approximation. The M1-closure scheme gives the Eddington tensor from the energy flux and energy density of neutrinos to close the moment equations of the radiative transfer up to the first order. Hence, the Eddington tensor in the M1-closure scheme plays a key role in evolving the energy flux of neutrinos. Because a direct comparison of the results of the radiation-hydrodynamics simulation with the Boltzmann neutrino transport and those with the M1-closure scheme is not in the scope of this paper, we compare the Eddington tensors. Hereafter, we call the Eddington tensor calculated with the M1-closure prescription the "M1–Eddington tensor" ${k}_{{\rm{M}}1}^{{ij}}(\epsilon )$ .

Here, we follow the M1-closure scheme suggested in Shibata et al. (2011). In the M1-closure scheme, two limiting cases are interpolated to obtain an approximate value of the stress tensor ${P}_{{\rm{M}}1}^{{ij}}(\epsilon )$ :

$\begin{eqnarray}&&{P}_{{\rm{M}}1}^{{ij}}(\epsilon ):= \displaystyle \frac{3\zeta (\epsilon )-1}{2}{P}_{\mathrm{thin}}^{{ij}}(\epsilon )+\displaystyle \frac{3(1-\zeta (\epsilon ))}{2}{P}_{\mathrm{thick}}^{{ij}}(\epsilon ),\end{eqnarray} \tag{ 36 }$

where ${P}_{\mathrm{thin}}^{{ij}}$ and ${P}_{\mathrm{thick}}^{{ij}}$ are the optically thin and thick limits of the stress tensor, respectively; $\zeta (\epsilon )$ is the Eddington factor, which is defined as the eigenvalue of the Eddington tensor along the flux direction. The M1-closure scheme gives the Eddington factor by an analytic function of the flux factor $\bar{F}(\epsilon ):=$ $\sqrt{{h}_{\alpha \beta }{H}^{\alpha }(\epsilon ){H}^{\beta }(\epsilon )/J{\left(\epsilon \right)}^{2}}$ :

$\begin{eqnarray}&&\zeta (\epsilon )=\displaystyle \frac{3+4\bar{F}{\left(\epsilon \right)}^{2}}{5+2\sqrt{4-3\bar{F}(\epsilon )}}.\end{eqnarray} \tag{ 37 }$

Here, $J(\epsilon )\,:=\,{u}_{\alpha }{u}_{\beta }{M}^{\alpha \beta }(\epsilon )$ and ${H}^{i}(\epsilon ):= -{h}_{\alpha }^{i}{u}_{\beta }{M}^{\alpha \beta }(\epsilon )$ are the energy density and flux in the fluid rest frame, respectively; ${u}^{\alpha }$ is the four-velocity of the fluid; ${h}_{\alpha \beta }:= {g}_{\alpha \beta }+{u}_{\alpha }{u}_{\beta }$ is the projection tensor onto the fluid rest frame. With this stress tensor ${P}_{{\rm{M}}1}^{{ij}}(\epsilon )$ , we define the M1–Eddington tensor as ${k}_{{\rm{M}}1}^{{ij}}(\epsilon ):= {P}_{{\rm{M}}1}^{{ij}}(\epsilon )/E(\epsilon )$ .

The optically thin and thick limits of the stress tensor are given as follows:

$\begin{eqnarray}&&{P}_{\mathrm{thin}}^{{ij}}(\epsilon ):= E(\epsilon )\displaystyle \frac{{F}^{i}(\epsilon ){F}^{j}(\epsilon )}{F{\left(\epsilon \right)}^{2}}\end{eqnarray} \tag{ 38 }$

at the optically thin limit and

$\begin{eqnarray}&&{P}_{\mathrm{thick}}^{{ij}}(\epsilon ):= J(\epsilon )\displaystyle \frac{{\gamma }^{{ij}}+4{V}^{i}{V}^{j}}{3}+{H}^{i}(\epsilon ){V}^{j}+{V}^{i}{H}^{j}(\epsilon )\end{eqnarray} \tag{ 39 }$

at the optically thick limit. Here, some projected quantities are also utilized: ${F}^{i}(\epsilon ):= -{\gamma }_{\alpha }^{i}{n}_{\beta }{M}^{\alpha \beta }(\epsilon )$ is the spatial-temporal projection, or the energy flux, in the laboratory frame; ${V}^{i}:= {u}^{i}/{u}^{0}$ is the three-velocity of the fluid.

In the following, we evaluate the M1–Eddington tensor from the second moment of the distribution function ${M}^{\alpha \beta }(\epsilon )$ : we first obtain the energy densities and fluxes $E(\epsilon )$ , $J(\epsilon )$ , ${F}^{i}(\epsilon )$ , and ${H}^{i}(\epsilon )$ by projecting ${M}^{\alpha \beta }$ then calculate the stress tensor ${P}_{{\rm{M}}1}^{{ij}}(\epsilon )$ from Equations (36)–(39). In the standard M1-closure scheme, however, we would not have the variables in both the laboratory and fluid rest frames. Usually, we only have the energy density and flux in the laboratory frame ( $E(\epsilon )$ and ${F}^{i}(\epsilon )$ ) instead of the entire second moment. Then, we guess the energy density and flux in the fluid rest frame ( ${J}_{\mathrm{guess}}(\epsilon )$ and ${H}_{\mathrm{guess}}^{i}(\epsilon )$ ), construct the stress tensor ${P}_{{\rm{M}}1,\mathrm{guess}}^{{ij}}(\epsilon )$ via Equation (36), construct the second moment ${M}_{\mathrm{guess}}^{\alpha \beta }(\epsilon )$ (see equation (3.29) in Shibata et al. 2011) from $E(\epsilon )$ , ${F}^{i}(\epsilon )$ , and ${P}_{{\rm{M}}1,\mathrm{guess}}^{{ij}}(\epsilon )$ , and obtain new guesses of ${J}_{\mathrm{guess}}(\epsilon )\,={u}_{\alpha }{u}_{\beta }{M}_{\mathrm{guess}}^{\alpha \beta }(\epsilon )$ and ${H}_{\mathrm{guess}}^{i}(\epsilon )=-{h}_{\alpha }^{i}{u}_{\beta }{M}_{\mathrm{guess}}^{\alpha \beta }(\epsilon )$ in the fluid rest frame. By iterating this correction numerically until the guess converges, we finally obtain the consistent energy density $J(\epsilon )$ , flux ${H}^{i}(\epsilon )$ , and stress tensor ${P}_{{\rm{M}}1}^{{ij}}(\epsilon )$ . Although we compare this iteratively obtained M1–Eddington tensor and our M1–Eddington tensor obtained from the second moment, the difference is only $\lesssim {10}^{-4}$ at ${t}_{\mathrm{pb}}=100\,\mathrm{ms}$ . Therefore, we do not use the iteratively obtained stress tensor for the M1–Eddington tensor.

It is worth noting that Equation (39) is derived from the perturbative expansion of the angular moments of the distribution function with respect to the local MFP. The second- or higher-order terms of the local MFP are ignored in the M1 prescription. Instead of including the higher-order terms in the semitransparent cases, the M1 scheme takes their effect into account by interpolating with the optically thin limit. Besides, some terms are ignored in Equation (39) compared to the original expression (Equation (6.19) in Shibata et al. 2011). This is just because such prescription makes the Boltzmann– and M1–Eddington tensors match better.

The M1–Eddington tensor shows overestimated filamental patterns compared with the Boltzmann–Eddington tensor. The Boltzmann– and M1–Eddington tensors and their difference for ${\nu }_{{\rm{e}}}$ are displayed in Figure 22. In the circles in the middle-left panel of Figure 22 ( ${\nu }_{{\rm{e}}}$ for the LS model at ${t}_{\mathrm{pb}}=100\,\mathrm{ms}$ ), red or white filaments are seen on the bluish background. Although the filaments are almost white or pale red ( ${k}^{r\theta }\gtrsim 0$ ) in the panel for the Boltzmann–Eddington tensor, the filaments in the panel for the M1–Eddington tensor are red ( ${k}_{{\rm{M}}1}^{r\theta }\gt 0$ ). The circles in the middle-right panel ( ${\nu }_{{\rm{e}}}$ for FS model at ${t}_{\mathrm{pb}}=100\,\mathrm{ms}$ ) show white or blue filaments on the reddish background. These filaments in the Boltzmann–Eddington tensor are white ( ${k}^{r\theta }\sim 0$ ), while those in the M1–Eddington tensor are blue ( ${k}_{{\rm{M}}1}^{r\theta }\lt 0$ ). Similar overestimated filamental patterns are also found in other regions or other panels.

**Figure 22.** Comparison of the rθ components of the Boltzmann–- and M1–Eddington tensors for ${\nu }_{{\rm{e}}}$ . The left column is for the LS model, while the right column is for the FS model. The top, middle, and bottom rows are taken from the snapshots at t_pb = 17.5 ms, 100 ms, and 300 ms, respectively. For each panel, the left and middle portions are the Boltzmann– and the M1–Eddington tensors, respectively; the right portion is the difference between them. They are all measured at the mean energies. The black circles in the middle panels indicate the filamental patterns discussed in the text. Note that the values are multiplied by a factor of 10 for the display.
Download figure:
Standard image High-resolution image

This overestimation comes from the limitation of the approximation. The central idea of the M1-closure scheme is the interpolation between ${P}_{\mathrm{thick}}^{{ij}}$ and ${P}_{\mathrm{thin}}^{{ij}}$ to approximate P^ij: the inequality $| {P}_{\mathrm{thick}}^{r\theta }| \leqslant | {P}^{r\theta }| \leqslant | {P}_{\mathrm{thin}}^{r\theta }|$ is assumed. In reality, $| {P}_{\mathrm{thin}}^{r\theta }| \leqslant | {P}_{\mathrm{thick}}^{r\theta }|$ holds in the overestimated filaments, and hence, the assumption of the M1-closure relation is violated.

The overestimation of $| {P}_{\mathrm{thick}}^{r\theta }|$ originates from the ignored higher-order terms of the local MFP. The origin of the filamental shapes seen in Figure 22 is the lateral matter velocity: the distribution of the neutrinos are slightly distorted by the lateral matter velocity in the semitransparent region via the scattering or the emission. In the M1-closure scheme, the matter motion is encoded through Vⁱ in Equation (39). Equation (39) includes up to first-order terms of the local MFP, and higher-order terms are ignored. In the filamental patterns in Figure 22, the ${H}^{r}{V}^{\theta }$ term becomes too large because the matter moves laterally in the semitransparent regions: H^r is no longer the first order of the MFP and ${V}^{\theta }$ is nonnegligible. This term should be canceled by the higher-order terms, but they are ignored in this formulation. Hence, $| {P}_{\mathrm{thick}}^{r\theta }|$ becomes too large for the assumption of the M1-closure scheme to hold. Similar discussions are also presented in Harada et al. (2019).

Richers et al. (2017) showed that the limited angular resolution of the S_N solver makes the off-diagonal component of the Eddington tensor underestimated. Thus, higher-resolution simulations probably show less difference between the Boltzmann– and M1–Eddington tensors. The possible violation of the assumption of the M1-closure relation, the inequality $| {P}_{\mathrm{thick}}^{r\theta }| \lt | {P}_{\mathrm{thin}}^{r\theta }|$ , is still problematic, however.

One of the peculiar features in Figure 22 is the very (negative) large value of ${k}_{{\rm{M}}1}^{r\theta }$ colored in dark blue in the top-left panel (LS model at ${t}_{\mathrm{pb}}=17.5\,\mathrm{ms}$ ). This panel corresponds to the first phase of the prompt convection, and the shock expanding from the equator rides over the shock in the upper region. A significant lateral motion to the negative θ direction arises there, and hence, such matter motion produces a large lateral flux to the negative θ direction. This lateral flux results in the large ${F}^{r}{F}^{\theta }$ and hence the large ${P}_{\mathrm{thin}}^{r\theta }$ . Besides, the large ${V}^{\theta }$ makes ${P}_{\mathrm{thick}}^{r\theta }$ large. Both of them result in the large ${k}_{{\rm{M}}1}^{r\theta }$ . Note that the large ${F}^{\theta }$ does not necessarily result in the large ${P}^{r\theta }$ . This is because the ${P}^{r\theta }$ is related to the quadrupole moment of the distribution function, while the Fⁱ is the dipole moment of the distribution.

The Eddington tensors of $\bar{{\nu }_{{\rm{e}}}}$ displayed in Figure 23 show similar background and opposite filamental patterns to those of ${\nu }_{{\rm{e}}}$ . The basic pattern of the signature, which is the background pattern behind the filamental patterns, is opposite between ${\nu }_{{\rm{e}}}$ (Figure 22) and $\bar{{\nu }_{{\rm{e}}}}$ (Figure 23). For example, the signs in the equatorial and northern/southern regions of the upper-left panel of Figure 23 ( $\bar{{\nu }_{{\rm{e}}}}$ for LS model at ${t}_{\mathrm{pb}}=17.5\,\mathrm{ms}$ ) are positive and negative, respectively, contrary to those of Figure 22; the sign in the upper region of the middle-left panel of Figure 23 ( $\bar{{\nu }_{{\rm{e}}}}$ for LS model at ${t}_{\mathrm{pb}}=100\,\mathrm{ms}$ ) is negative while that of Figure 22 for ${\nu }_{{\rm{e}}}$ is positive. On the other hand, the filamental patterns for ${\nu }_{{\rm{e}}}$ and $\bar{{\nu }_{{\rm{e}}}}$ are similar. The filaments that are enclosed by the black circles in the middle panels of Figure 22 ( ${\nu }_{{\rm{e}}}$ at ${t}_{\mathrm{pb}}=100\,\mathrm{ms}$ ) are also found in the middle panels of Figure 23 ( $\bar{{\nu }_{{\rm{e}}}}$ at ${t}_{\mathrm{pb}}=100\,\mathrm{ms}$ ).

The $r\theta$ component of the Eddington tensor is determined by the emission from the optically thick region, whose signature is different between ${\nu }_{{\rm{e}}}$ and $\bar{{\nu }_{{\rm{e}}}}$ , and the motion of matter, whose signature is the same between the two neutrino species. The different signs of the background pattern discussed above originate from the degeneracy of ${\nu }_{{\rm{e}}}$ discussed in Section 4.1. In the optically thick region, the number density of $\bar{{\nu }_{{\rm{e}}}}$ is large in the region where the ${\nu }_{{\rm{e}}}$ number density is small, and the lateral component of the diffusion flux of $\bar{{\nu }_{{\rm{e}}}}$ is opposite to that of ${\nu }_{{\rm{e}}}$ . This is imprinted not only in the flux but also the $r\theta$ component of the Eddington tensor. For the semitransparent to optically thin region, though ${\nu }_{{\rm{e}}}$ is no longer degenerate, the pattern generated in the optically thick region is transported and shows the same signature as the inner regions. The background pattern is hence determined by the emission from the optically thick region. On the other hand, the filamental patterns are induced by the motion of matter. Since the motion of matter is common for both neutrino species, the filamental patterns for ${\nu }_{{\rm{e}}}$ and $\bar{{\nu }_{{\rm{e}}}}$ are similar.

The Eddington tensors for ${\nu }_{x}$ are shown in Figure 24. The filamental pattern is similar to ${\nu }_{{\rm{e}}}$ and $\bar{{\nu }_{{\rm{e}}}}$ . This is because the motion of matter again determines the signature there. However, the pattern in the other regions does not have a clear (anti)correlation with ${\nu }_{{\rm{e}}}$ and $\bar{{\nu }_{{\rm{e}}}}$ . Because the number density is solely determined by the temperature and has nothing to do with the chemical potential of ${\nu }_{{\rm{e}}}$ , the resultant pattern is again not so related to that of ${\nu }_{{\rm{e}}}$ and $\bar{{\nu }_{{\rm{e}}}}$ .

**Figure 24.** Same as Figures 22 and 23, except that the displayed neutrino species is ${\nu }_{x}$ .
Download figure:
Standard image High-resolution image

5. Summary and Conclusions

In this paper, we have discussed the postbounce dynamics and the neutrino properties of the models with different EOSs. They are summarized as follows:

1.
The model with the LS EOS shows shock revival, while that with the FS EOS does not. The neutrino luminosities, energies, heating rates, the mass accretion rates, and total energies in the gain region of the two models are similar, while the gain mass and strength of the turbulence are different. This originated from the difference in the nuclear composition of the accreting matter. The estimated binding energy is larger for the LS model, and hence more energy is consumed to photodissociate them. It results in a steeper entropy gradient, which drives stronger prompt convection. It enhances the later neutrino-driven convection.
2.
The structure of the PNS implies the importance of the composition. The central region of the PNS in the LS model is more compact than that in the FS model. However, the outer region in the LS model is slightly less compact than that in the FS model. This is understood from the effective stiffness. The nuclear composition under the NSE has more degrees of freedom than that under the SNA. Because the FS EOS considers the NSE, the free energy in the FS model is lower than that in the LS model. Therefore, the pressure gets lower in the FS model, and hence, the FS EOS at subnuclear densities is softer than the LS EOS.
3.
The neutrino flux is determined by the diffusion of the neutrinos, whose direction is along the gradient of the neutrino number density, and the matter velocity via the Lorentz transformation. The ${\nu }_{{\rm{e}}}$ number density decreases with the radius, and hence the diffusion flux is directed from the rising fluid parcels to the sinking parcels when the prompt convection develops. On the other hand, the $\bar{{\nu }_{{\rm{e}}}}$ number density increases with the radius in the region where prompt convection occurs. This is because the ${\nu }_{{\rm{e}}}$ number density and hence the degeneracy of ${\nu }_{{\rm{e}}}$ decreases with the radius. Therefore, the diffusion flux is directed from sinking to rising fluid parcels. For the heavy-lepton-type neutrinos, their number density is determined by the zero-chemical-potential thermal distribution: it increases with the temperature. Therefore, the diffusion flux is similar to that of electron-type antineutrinos, although is slightly different because the peak in the number density lies in the convectively unstable region.
4.
We compared the Boltzmann– and M1–Eddington tensors and found that the effect of the matter velocity on the Eddington tensor is overestimated in the M1–Eddington tensor. This is mainly because the higher-order terms of the local MFP ignored in Equation (39) is too large for the M1 approximation to hold. The Eddington tensor itself is affected by both the diffusion of neutrinos and the matter velocity. The contribution of the neutrino diffusion is opposite between ${\nu }_{{\rm{e}}}$ and $\bar{{\nu }_{{\rm{e}}}}$ due to the ${\nu }_{{\rm{e}}}$ degeneracy. The contribution of the diffusion of ${\nu }_{x}$ has nothing to do with ${\nu }_{{\rm{e}}}$ and $\bar{{\nu }_{{\rm{e}}}}$ because the distribution of ${\nu }_{x}$ is solely determined by the temperature, not the neutrino chemical potential. On the other hand, the contribution of the matter velocity is common between them. This difference explains the distribution of the $r\theta$ component of the Eddington tensors of ${\nu }_{{\rm{e}}}$ , $\bar{{\nu }_{{\rm{e}}}}$ , and ${\nu }_{x}$ .

We have conducted thorough analyses of the simulations with a nonrotating $11.2\,{M}_{\odot }$ progenitor with the LS and FS EOSs in this paper, but systematic studies about the input physics would be also required. We are further running the Boltzmann-radiation-hydrodynamics code with a variety of input physics: different progenitors, different initial rotations, and different EOSs. Collecting deep analyses of these simulations, we can obtain valuable clues for understanding the explosion mechanism of CCSNe.

In addition, we are continuously developing the code. For 2D simulations, improvement in the feedback from neutrino–matter interaction and its influence on the PNS kick are discussed in Nagakura et al. (2019b, 2019c). These simulations employ one of the currently most realistic EOS models (Togashi et al. 2017). We have developed a 3D version of the Boltzmann-radiation-hydrodynamics code, and it can follow the first few tens of milliseconds (Iwakami et al. 2020). A general relativistic version of the code is also under development. The spherical polar coordinate employed in our code has some difficulty at the center and along the pole when running the numerical relativity, but Baumgarte et al. (2013) has suggested a method to avoid it. We are developing the general relativistic version of our code using their technique. With the coming next-generation supercomputers, we hope to run the general relativistic Boltzmann-radiation-hydrodynamics code without any spatial symmetry. The simulations with the code would provide deep understanding of the effect of each physical process on the dynamics of SN explosions.

We thank K. Nakazato, Y. Suwa, and T. Takiwaki for fruitful discussions. This research used high-performance computing resources of the K-computer and the FX10 of the HPCI system provided by the R-CCS and the University of Tokyo through the HPCI System Research Project (Project ID: hp 140211, 150225, 160071, 160211, 170230, 170031, 170304, 180111, 180239), SR16000 and XC40 at YITP of Kyoto University, SR16000 and Blue Gene/Q at KEK under the support of its Large Scale Simulation Program (14/15-17, 15/16-08, 16/17-11), Research Center for Nuclear Physics (RCNP) at Osaka University, and the XC30 and the general common-use computer system at the Center for Computational Astrophysics, CfCA, the National Astronomical Observatory of Japan. The large-scale storage of numerical data is supported by JLDG constructed over SINET of NII. This work was supported in part by a Grant-in-Aid for Scientific Research on Innovative areas "Gravitational wave physics and astronomy: Genesis" (17H06357, 17H06365) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan, and in part a Grant-in-Aid for Scientific Research (C; 15K05093, 19K03837, B; 20H01905) and Young Scientists (Start-up, JP19K23435) from the Japan Society for the Promotion of Science (JSPS). This work was also partly supported by research programs at K-computer of the RIKEN R-CCS, HPCI Strategic Program of Japanese MEXT, Priority Issue on Post-K-computer (Elucidation of the Fundamental Laws and Evolution of the Universe), Joint Institute for Computational Fundamental Sciences (JICFus), and by MEXT as "Program for Promoting Researches on the Supercomputer Fugaku" (Toward a unified view of the universe: from large scale structures to planets). A.H. was supported in part by the Advanced Leading Graduate Course for Photon Science (ALPS) in the University of Tokyo and a Grant-in-Aid for JSPS Research Fellows (JP17J04422). K.S. acknowledges the Particle, Nuclear and Astro Physics Simulation Program (2019-002, 2020-004) at KEK.

Software: matplotlib (Hunter 2007).

The Boltzmann-radiation-hydrodynamics Simulations of Core-collapse Supernovae with Different Equations of State: The Role of Nuclear Composition and the Behavior of Neutrinos

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Numerical Modeling

3. Dynamical Features

3.1. Overview of The Dynamics

3.2. Deformation of the Shock with Convection

3.3. Timescale Ratio and the Effects of the Nuclear Composition of the EOSs

3.4. Structure of the PNS and the Influence of the EOSs

4. Neutrino Distributions

4.1. Neutrino Flux

4.2. Eddington Tensor

5. Summary and Conclusions

Footnotes

The Boltzmann-radiation-hydrodynamics Simulations of Core-collapse Supernovae with Different Equations of State: The Role of Nuclear Composition and the Behavior of Neutrinos

Article metrics

Permissions

Share this article

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Numerical Modeling

3. Dynamical Features

3.1. Overview of The Dynamics

3.2. Deformation of the Shock with Convection

3.3. Timescale Ratio and the Effects of the Nuclear Composition of the EOSs

3.4. Structure of the PNS and the Influence of the EOSs

4. Neutrino Distributions

4.1. Neutrino Flux

4.2. Eddington Tensor

5. Summary and Conclusions

Footnotes