Simulations of core-collapse supernovae in spatial axisymmetry with full Boltzmann neutrino transport

TL;DR

This paper presents the first axisymmetric supernova simulations with full Boltzmann neutrino transport, revealing detailed neutrino distributions and the impact of nuclear EOS on shock evolution.

Contribution

It introduces a full 5D neutrino transport simulation in axisymmetry, incorporating relativistic effects and comparing different nuclear equations-of-state.

Findings

Shock wave reaches ~700km in LS EOS simulation

More vigorous turbulence in LS EOS leads to higher neutrino heating

Non-axisymmetric neutrino angular distributions affect flux evolution

Abstract

We present the first results of our spatially axisymmetric core-collapse supernova simulations with full Boltzmann neutrino transport, which amount to a time-dependent 5-dimensional (2 in space and 3 in momentum space) problem in fact. Special relativistic effects are fully taken into account with a two-energy-grid technique. We performed two simulations for a progenitor of 11.2M, employing different nuclear equations-of-state (EOS's): Lattimer and Swesty's EOS with the incompressibility of K = 220MeV (LS EOS) and Furusawa's EOS based on the relativistic mean field theory with the TM1 parameter set (FS EOS). In the LS EOS the shock wave reaches ~700km at 300ms after bounce and is still expanding whereas in the FS EOS it stalled at ~200km and has started to recede by the same time. This seems to be due to more vigorous turbulent motions in the former during the entire post-bounce phase, which leads to higher neutrino-heating efficiency in the neutrino-driven convection. We also look into the neutrino distributions in momentum space, which is the advantage of the Boltzmann transport over other approximate methods. We find non-axisymmetric angular distributions with respect to the local radial direction, which also generate off-diagonal components of the Eddington tensor. We find that the r {\theta}-component reaches ~10% of the dominant rr-component and, more importantly, it dictates the evolution of lateral neutrino fluxes, dominating over the {\theta}{\theta}-component, in the semi-transparent region. These data will be useful to further test and possibly improve the prescriptions used in the approximate methods.

Equations

$$M^{\mu \nu }(\epsilon )\equiv \frac{1}{\epsilon }\int f(\epsilon ,{\Omega }_m)p^{\mu }{p}^{\nu }d{\Omega }_m,$$

$$k^{ij}(\epsilon )\equiv \frac{P^{ij}(\epsilon )}{E(\epsilon )},$$

$$P^{ij}(\epsilon )\equiv {\gamma }_{\mu }^i{\gamma }_{\nu }^j{M}^{\mu \nu }(\epsilon ),$$

$$E(\epsilon )\equiv n_{\mu }{n}_{\nu }{M}^{\mu \nu }(\epsilon ),$$

$$P_{M1}^{ij}(\epsilon )=\frac{3\zeta (\epsilon )-1}{2}{P}_{\mathrm{thin}}^{ij}(\epsilon )+\frac{3(1-\zeta (\epsilon ))}{2}{P}_{\mathrm{thick}}^{ij}(\epsilon ),$$

$$\zeta (\epsilon )=\frac{3+4\bar{F}(\epsilon {)}^2}{5+2\sqrt{4-3\bar{F}(\epsilon {)}^2}}.$$

$$\displaystyle\bar{F}(\varepsilon)=\Biggl{(}\frac{h_{\mu\nu}H^{\mu}(\varepsilon)H^{\nu}(\varepsilon)}{J(\varepsilon)^{2}}\Biggr{)}^{1/2},$$

$$J(\epsilon )=u_{\mu }{u}_{\nu }{M}^{\mu \nu }(\epsilon ),$$

$$H^{\mu }(\epsilon )=-h_{\alpha }^{\mu }{u}_{\beta }{M}^{\alpha \beta }(\epsilon ),$$

$$P_{\mathrm{thick}}^{ij}(\epsilon )=J(\epsilon )\frac{{\gamma }^{ij}+4V^i{V}^j}{3}+H^i(\epsilon )V^j+V^i{H}^j(\epsilon ),$$

$$P_{\mathrm{thin}}^{ij}(\epsilon )=E(\epsilon )\frac{F^i(\epsilon )F^j(\epsilon )}{F(\epsilon {)}^2},$$

$$F^i(\epsilon )=-{\gamma }_{\mu }^i{n}_{\nu }{M}^{\mu \nu }(\epsilon ).$$

$${\partial }_t(F^r)\sim -{\partial }_r(E{k}^{rr})-\frac{1}{r}{\partial }_{\theta }(E{k}^{r\theta }),$$

$${\partial }_t(F^{\theta })\sim -{\partial }_r(E{k}^{r\theta })-\frac{1}{r}{\partial }_{\theta }(E{k}^{\theta \theta }),$$

Full text

supernovae: general—neutrinos—hydrodynamics

1 Introduction

The theoretical study of the explosion mechanism of core-collapse supernovae (CCSNe) has heavily relied on numerical simulations. This is mainly because nearby CCSNe are rare (van den Bergh & Tammann, 1991; Cappellaro et al., 1993; Tammann et al., 1994; Reed, 2005; Diehl et al., 2006; Maoz & Badenes, 2010; Li et al., 2011) and, in fact, SN1987A is the only one close enough to extract some useful information on what happened deep inside the massive star from, among other things, the detection of neutrinos (Bionta et al., 1987; Hirata et al., 1987). Since the CCSNe are intrinsically multi-scale, multi-physics and multi-dimensional (multi-D) phenomena, their mechanism can be addressed only with detailed numerical simulations.

Unfortunately, even the most advanced multi-D simulations of CCSNe employed approximations one way or another in their numerical treatment of neutrino transport (Marek & Janka, 2009; Müller et al., 2012; Bruenn et al., 2013; Takiwaki et al., 2014; Bruenn et al., 2016; Dolence et al., 2015; Lentz et al., 2015; Melson et al., 2015; Kuroda et al., 2016; Skinner et al., 2016; O’Connor & Couch, 2015; Pan et al., 2016; Just et al., 2015; Summa et al., 2016; Roberts et al., 2016; Andresen et al., 2017; Burrows et al., 2016). Most of them somehow integrated out the angular degrees of freedom in momentum space or neglected non-radial fluxes in neutrino transport. Ott et al. (2008) is the only exception, in which they conducted time-dependent 5-dimensional simulations in spatial axisymmetry. However, they ignored relativistic corrections completely, dropping all fluid-velocity-dependent terms, which are crucial for qualitatively correct descriptions of the angular distribution of neutrinos in momentum space (see e.g., Buras et al. (2006); Lentz et al. (2012)).

The best way to calibrate all these approximate methods should be to compare them with simulations that solve full Boltzmann equations, retaining the angular degree of freedom, for neutrino transport. Under axisymmetry in space, this is possible now indeed and we have achieved such simulations with the K computer in Japan, one of the currently available best supercomputers with $\sim 10$PFLOPS. The validation of our Boltzmann solver has been conducted in a series of papers: the standard tests in static matter distributions meant for radiation transport codes were done in Sumiyoshi & Yamada (2012); Nagakura et al. (2014) coupled the Boltzmann solver with a hydrodynamics code of their own construction and tested in dynamical settings the capability of the integrated code in treating special relativistic effects; Nagakura et al. (2017), on the other hand, tested a new module implemented for the tracking of the motion of a proto neutron star (PNS) with a moving grid; very recently Richers et al. (2017) made a detailed comparison with another Boltzmann solver based on the Monte Carlo method using snapshots from our 2D and 1D simulations and calculating steady-state neutrino distribution functions in the static fluid backgrounds. Having established reliability of our code with these test computations, we now proceed to present the first series of multi-D simulations of CCSNe with the full Boltzmann neutrino transport. In this paper we also pay attention to the neutrino angular distributions in momentum space, which are what the Boltzmann solver is meant for in the first place. Throughout this paper, Greek and Latin subscripts denote spacetime and spatial components, respectively. We use the metric signature $-+++$. Unless otherwise stated, we work in units with $c=G=1$, where $c$ is the speed of light and $G$ is the gravitational constant.

2 Methods and Models

We solve numerically the equations of neutrino-radiation hydrodynamics. We apply the so-called discrete-ordinate (DO) method to the Boltzmann equations for neutrino transport, taking fully into account special relativistic effects by virtue of a two-energy-grid technique (Nagakura et al., 2014). It has already incorporated general relativistic capabilities as well, a part of which is utilized to track the proper motion of PNS (Nagakura et al., 2017). The hydrodynamics and self-gravity are still Newtonian: the so-called central scheme of second-order accuracy in both space and time is employed for the former and the Poisson equation is solved for the latter.

It should be noted that our treatment of neutrino transport is essentially different from other approximate methods such as the M1 scheme that are commonly employed in the currently most elaborate supernova simulations and are based on the truncated moment formalism one way or another. It is combined with the ray-by-ray approximation in some applications (see e.g., Müller et al. (2012)). In the moment formalism the Boltzmann equation is angle-integrated in momentum space to obtain an infinite number of equations for angular moments, which are then truncated at some order somehow (see Sec. 4 for more details). Such approximations reduce the computational cost drastically. On the other hand, they inevitably introduce the so-called closure relation among lo-order moments, which are the artificial prescriptions to make the truncated equations self-contained. Although the validity of those prescriptions has been assessed for spherically symmetric cases in the literature (Richers et al., 2017), it remains to be demonstrated in multi-dimensional and, more importantly, dynamical settings. In sharp contrast, our approach does not employ any such artificial prescription in the neutrino transport except for the finite-discretization of the Boltzmann equation, which is all but mandatory for this sort of simulations.

We adopt spherical coordinates $(r,\theta)$ covering $0\leq r\leq 5000{\rm km}$ and $0^{\circ}\leq\theta\leq 180^{\circ}$ in the meridian section. We deploy $384(r)\times 128(\theta)$ grid points. Momentum space is also discretized non-uniformally with $20$ energy mesh points covering $0\leq\varepsilon\leq 300{\rm MeV}$ and $10(\bar{\theta})\times 6(\bar{\phi})$ angular grid points over the entire solid angle. The polar and azimuthal angles $(\bar{\theta},\bar{\phi})$ are locally measured from the radial direction. Three neutrino species are distinguished: electron-type neutrinos $\nu_{\rm e}$, electron-type anti-neutrinos $\bar{\nu}_{\rm e}$ and all the others collectively denoted by $\nu_{x}$.

We pick up a non-rotating progenitor model of 11.2 $M_{\odot}$ from Woosley et al. (2002). We employ two nuclear EOS’s: Lattimer & Swesty’s EOS with the incompressibility of $K=220{\rm MeV}$ (Lattimer & Swesty, 1991) and Furusawa’s EOS derived from H. Shen’s relativistic mean-field EOS with the TM1 parameter set (Furusawa et al., 2011, 2013); the former is softer than the latter (see Sumiyoshi et al. (2004)). In the following, they are referred to as the ”LS” and ”FS” EOS’s, respectively111The maximum gravitational masses at zero temperature and non-rotating neutron stars are $2.02M_{\odot}$ for LS EOS and $2.21M_{\odot}$ for FS EOS, respectively.. The choice of EOS’s is simply based on the fact that most of previous simulations employed one of these EOS’s. We are currently running similar simulations, but with another EOS: Togashi’s nuclear EOS based on the variational method with realistic nuclear potentials (Togashi & Takano, 2013) extended by Furusawa et al. (2017) to sub-nuclear densities; it takes into account the full ensemble of heavy nuclei in nuclear-statistical equilibrium (NSE). The results will be reported elsewhere (Nagakura et al., 2018). Neutrino-matter interactions are based on those given by Bruenn (1985) but we have implemented the up-to-date electron capture rates for heavy nuclei (Juodagalvis et al., 2010; Langanke & Martínez-Pinedo, 2000; Langanke et al., 2003); they are calculated based on the abundance of heavy nuclei obtained in the FS EOS; the same rates are employed in the LS EOS model just for simplicity; note also that the LS EOS employs a single-nucleus approximation and the detailed information on the population of various nuclei is unavailable. In the current simulations we incorporated the non-isoenergetic scatterings on electrons and positrons as well as the bremsstrahlung in nucleon collisions. We refer readers to Nagakura et al. (2014, 2017); Sumiyoshi & Yamada (2012) for more details of our code.

We start the simulations in spherical symmetry and switch them to axisymmetric computations at $\sim 1{\rm ms}$ after core bounce when a negative entropy gradient starts to develop behind the shock wave. We seed by hand at this point of time perturbations of $0.1\%$ in the radial velocities at $30\leq r\leq 50{\rm km}$, where convection is expected to occur (see Fig. 2). Note that we do not explicitly consider possible turbulent motions that have already existed in the progenitors before collapse. We then expect in non-rotating models that non-radial motions develop initially in the convectively unstable region, and then spread in the rest of the post-shock flow. Each model is run up to $t=300{\rm ms}$ after bounce.

3 Dynamics

As displayed in Fig. 1(a), the shock wave produced at core bounce expands rather gradually with time for the LS EOS and its maximum radius reaches $\sim 700{\rm km}$ at $t=300{\rm ms}$. For the FS EOS, on the other hand, the shock wave stalls at $r\sim 200{\rm km}$ at $t\sim 100{\rm ms}$ and then starts to recede at $t\sim 250{\rm ms}$ and shrinks back to $r\sim 100{\rm km}$ by $t\sim 300{\rm ms}$. Although the time evolutions of the average shock radii of the two models are quite similar to each other until $\sim 60{\rm ms}$ after bounce and their deviations become remarkable thereafter, some differences have already appeared in the post-shock flows by this time in fact.

In the top two panels of Fig. 2, we compare the angle-averaged amplitudes of lateral velocity for the two models. The more reddish the color is, the stronger the lateral motions are. It is apparent that they become appreciable initially at $t\sim 10{\rm ms}$ in a region at $r\sim 30{\rm km}$ almost simultaneously, which marks the onset of the prompt convection. Although the turbulent region extends upwards in both models, the amplitudes of the lateral velocity are larger for the LS EOS than for the FS EOS, indicating that the prompt convection is more vigorous in the former. This trend persists until much later times, though, as is also evident from the figure.

The difference in the strength of the prompt convection may be understood from the difference in the Brunt-Väisälä frequencies, which are compared in the lower panels of Fig. 2. Reddish colors again imply more rapid (exponential) growths of the convection. As expected, the unstable region emerges at $r\sim 30-50{\rm km}$ immediately after the switch to the 2D computations for both models. Although this strongly unstable region persists until $t\sim 15{\rm ms}$ at about the same location for both models, the maximum frequency is larger for the LS EOS. This difference can be traced back to the difference in photodissociations of heavy nuclei by shock heating. In fact, they are stronger in the LS EOS and, as a result, the shock is weakened more severely, producing steeper negative entropy gradients in this case. The initial fluctuations produced this way are carried upwards by acoustic waves, which are also stronger in the LS EOS. As a consequence, the turbulent motions are more vigorous for the LS EOS than for the FS EOS as already mentioned, the fact that has an important implication for later evolutions of the shock waves.

It is interesting that the neutrino luminosities ($L$) and mean energies ($E_{m}$, defined as the ratio of energy density to number density) are almost identical between the two cases (Fig. 1(b)). It should be noted, however, that the neutrino-heating efficiency is different, being higher for the LS EOS (see solid lines in Fig. 1(c)). This is mainly because the total baryon mass in the gain region, where heating dominates over cooling and the net heating occurs, is consistently larger for the LS EOS than for the FS EOS (dashed lines in the same panel). This in turn seems to be a consequence of the turbulent motions that are more vigorous for the LS EOS as we mentioned in the previous paragraphs.

Figure 3 compares the entropy and velocity distributions between the two models at $t=200{\rm ms}$. Their post-shock morphologies are quite similar to each other and only the scales are different. In fact, the convection is dominant over SASI in most of the post-bounce phase for both models (see the $\chi$ parameter (Foglizzo et al., 2006; Iwakami et al., 2014) in Fig. 1(d)). In the same panel, we also show the ratio of the advection timescale ($T_{\rm adv}=M_{\rm g}/\dot{M}$ with $M_{\rm g}$ and $\dot{M}$ denoting the mass in the gain region and the mass accretion rate, respectively) to the heating timescale ($T_{\rm heat}=|E_{\rm tot}|/\dot{Q}_{\nu}$ with $E_{\rm tot}$ and $\dot{Q}_{\nu}$ being the total energy and the heating rate in the gain region, respectively) as solid lines. One can see that it is consistently larger for the LS EOS than for the FS EOS, meaning that the former has more favorable conditions for shock revival than the latter.

The decline of this ratio near the end of the simulation for the LS EOS in spite of a continuous growth of the maximum shock radius is an artifact originated from our choice of the minimum shock radius in the evaluation of the ratio. As displayed in Fig. 1(a), the minimum shock radius is still decreasing with time at the end of the simulation. Then the volume of gain region is underestimated and, as a result, $T_{\rm heat}$ is overestimated. The fact that the ratio occasionally exceeds unity but still yields no shock revival for the FS EOS indicates that the criterion is not a rigorous condition, which is understood also from the uncertainty in its definition just mentioned. We do not intend to discuss the applicability of the diagnostics any further in this paper but we still think it is useful in judging, albeit roughly, which model is closer to shock revival.

4 $\nu$-Distributions in Momentum Space

Next we turn our attention to novel features of the neutrino distributions in momentum space. We find in our calculations significant non-axisymmetry with respect to the radial direction in the neutrino angular distributions. It is produced by lateral inhomogeneities in matter, which are in turn generated by hydrodynamical instabilities. The asymmetry hence appears inevitably in multi-D simulations.

Figure 4 shows as an example the angular distributions of $\nu_{\rm e}$ with an energy of $\varepsilon=11.1{\rm MeV}$ at three different radial positions. Each surface displays the neutrino distribution function for different propagation directions normalized by the maximum value in the fluid-rest frame. Colors of the surfaces denote the locations on an arbitrarily chosen radial ray. The angular distribution is almost isotropic at $r=23{\rm km}$ (red surface) while they become forward peaked (green and blue surfaces) as the radius increases, a fact that is well known. What is really new here is that they are non-axisymmetric with respect to the radial direction, which is more apparent in Fig. 5, in which the isotropic contributions are subtracted from the original distributions and the resultant ones are re-normalized by their maximum values. Note that the feature is robust, occurring irrespective of neutrino energies or species.

It should be mentioned, however, that the non-axisymmetric angular distributions obtained in the current simulations still have a symmetry with respect to the azimuthal angle ($\bar{\phi}$) in momentum space. This is due to the fact that these are non-rotating models and there is a mirror symmetry with respect to the plane spanned by \mbox{\boldmath\bar{e}}_{r} and \mbox{\boldmath\bar{e}}_{\theta} in momentum space in the absence of rotation. Once rotation is taken into account, the symmetry is lost even in (spatial) axisymmetry. This is the reason why we do not assume this symmetry in our code. In 3D simulations, no symmetry remains in the angular distribution in momentum space. Its characterization is an interesting subject of spatially 3D supernova simulations with multi-angle neutrino transport, which are currently being undertaken and will be reported elsewhere later.

The multi-angle treatment of neutrino transport in our simulations enables us to evaluate the so-called Eddington tensor ($k^{ij}$), which characterizes these non-axisymmetric angular distributions more quantitatively. The Eddington tensor is obtained from the neutrino distribution function ($f$) as follows: we first define the second angular moment $M^{\mu\nu}$ as

[TABLE]

where $p^{\mu}$ is the four-momentum of neutrino and $\varepsilon$ and $\Omega_{m}$ are the corresponding energy and solid angle measured in the fluid-rest frame; then the Eddington tensor $k^{ij}$ is given as

[TABLE]

where $P^{ij}$ and $E$ are defined from $M^{\mu\nu}$ as

[TABLE]

with $n_{\mu}$ and $\gamma_{\hskip 2.84526pt\mu}^{i}$ (= $\delta_{\hskip 2.84526pt\mu}^{i}+n^{i}n_{\mu})$ being the unit vector orthogonal to a hypersurface of constant coordinate time and the projection tensor onto this hypersurface, respectively.

We pay particular attention here to one of the off-diagonal components of the Eddington tensor, $k^{r\theta}$, which are zero in spherical symmetry in space, i.e., they are a measure of genuinely multi-dimensional transfer. The left panel in Fig. 6(a) shows $k^{r\theta}$ for $\nu_{\rm e}$ with the mean energy at each point. As expected, it is almost zero inside the PNS, where matter is opaque enough to make the neutrino distribution isotropic. It becomes non-zero outside the PNS, however, and increases with radius in accord with the appearance of the non-axisymmetric structures in the neutrino angular distribution (see Fig. 4). In fact, the $k^{r\theta}$ corresponds to the mode with $\ell=2,m=1$ in the spherical harmonics expansion of the distribution function.

The right panel in Fig. 6(a) compares $k^{r\theta}$ obtained from our simulation with that which is evaluated according to the M1 prescription: the Eddington tensor in the M1 prescription ($k^{ij}_{\rm M1}$) is obtained by replacing $P^{ij}$ in Eq. (3) with

[TABLE]

where $\zeta$ is referred to as the variable Eddington factor, which we set as

[TABLE]

In this expression, $\bar{F}$ denotes the so-called flux factor, which is the energy-flux normalized with the energy density in the fluid-rest frame. The flux factor that we use in this paper is measured in the fluid-rest frame (see Shibata et al. (2011) for another option);

[TABLE]

where $J$ and $H^{\mu}$ can be expressed in terms of $M^{\mu\nu}$ as

[TABLE]

with $u^{\mu}$ and $h_{\hskip 2.84526pt\nu}^{\mu}(=\delta_{\hskip 2.84526pt\nu}^{\mu}+u^{\mu}u_{\nu})$ being the fluid four velocity and the projection tensor onto the fluid-rest frame, respectively. The optically thick and thin limits of $P^{ij}$ are denoted by $P^{ij}_{\rm thick}$ and $P^{ij}_{\rm thin}$ (Just et al., 2015; Shibata et al., 2011; O’Connor & Couch, 2015; Kuroda et al., 2016), which are written as

[TABLE]

where $V^{i}$ denotes the three dimensional vector of fluid velocity. $F^{i}$ can be expressed in terms of $M^{\mu\nu}$ as

[TABLE]

As clearly seen in this panel, the values of $k^{r\theta}$ are substantially different between the two cases. We find that such discrepancies in $k^{r\theta}$ are rather generic, being insensitive to the choice of the prescription for the Eddington factor (see Just et al. (2015) for various options). They are also systematic in the sense that the increase in the number of grid points in the M1 prescription does not reduce the difference. This is in contrast to our approach, in which the accuracy is simply improved with the resolution.

Moreover, we find in $k^{r\theta}$ an intriguing correlation/anti-correlation between $\nu_{\rm e}$ and $\bar{\nu}_{\rm e}$. The two panels of Fig. 6(b) compare $k^{r\theta}$ for $\nu_{\rm e}$ and $\bar{\nu}_{\rm e}$ with the same energy of $\varepsilon=8.5{\rm MeV}$. As can be seen in these panels, they are anti-correlated with each other in the vicinity of PNS ($\lesssim 50{\rm km}$) whereas they are positively correlated at larger radii ($>80{\rm km}$). The anti-correlation is particularly remarkable for low-energy neutrinos with $\lesssim 10{\rm MeV}$. We find that the sign of $k^{r\theta}$ roughly coincides with that of the lateral neutrino flux, which is shown in Fig. 6(c). In fact, it is apparent that the lateral flux is oriented in the opposite directions for $\nu_{\rm e}$ and $\bar{\nu}_{\rm e}$. This is in turn due to the Fermi-degeneracy of $\nu_{\rm e}$ at $r\lesssim 30{\rm km}$, which produces opposite trends in the number densities of $\nu_{\rm e}$ and $\bar{\nu}_{\rm e}$. Since neutrinos flow from high to low $\nu$ number density regions in the diffusion regime, the fluxes of $\nu_{\rm e}$ and $\bar{\nu}_{\rm e}$ should be naturally anti-correlated as a result of the opposite trend in the number densities of $\nu_{\rm e}$ and $\bar{\nu}_{\rm e}$. We do not know for the moment how this anti-correlation in the fluxes is transferred to that in $k^{r\theta}$. It will be necessary to analyze more in detail the equations of motion for higher moments including $k^{r\theta}$.

Importantly, the anti-correlation is then carried to larger radii by the radial flux and remains non-vanishing even at $r\sim 50{\rm km}$, where $\nu_{\rm e}$ is no longer degenerate. On the other hand, at even larger radii, where matter is optically thin to neutrinos, $k^{r\theta}$ is correlated with the local lateral velocity of matter due to relativistic aberration. Note that this positive correlation at large distances is less remarkable than the anti-correlation in the vicinity of PNS (see the equatorial region in Fig. 6(b)), since the angular distribution is no longer determined locally and the correlation is somewhat smeared out.

As will be discussed in Sec. 6, the appropriate treatment of $k^{r\theta}$ is related with the accurate calculation of the neutrino flux, in particular its lateral component (see Eqs. (11) and (12)). It is true that these correlation/anti-correlation look rather minor but they may play an important role through the lateral fluxes of neutrinos. In fact they clearly indicate the intricacy of neutrino transport in non-spherically dynamical settings. It will be interesting to see how well the M1 scheme can reproduce these features and to conceive possible improvements of its prescription.

5 Angular Resolution in momentum space

This study is the first ever attempt to perform spatially 2D supernova simulations with multi-angle and multi-energy neutrino transport, taking into account all special relativistic effects completely. It is a legitimate concern, however, that the current simulations may not have a sufficient numerical resolution especially in momentum space (Richers et al., 2017). In this section we hence discuss this resolution issue, focusing on the angular resolution in momentum space.

For that purpose we perform a new high-resolution simulation for the early post-bounce phase, whereas for the discussion of the late post-bounce phase we employ the results of our previous analyses (Richers et al., 2017) of time-independent solutions of the Boltzmann equations for neutrinos in given matter distributions; close comparisons were made with the data obtained with Monte Carlo Simulations (Richers et al., 2015). Note that although the use of the time-independent solutions for the fixed matter distributions enabled us to conduct rigorous comparisons, its applicability may be limited to the late post-bounce phase, where the time scale of variations in the background is indeed long. For the earlier phase, however, we need to consider time-dependent solutions. We hence run a higher-resolution simulation, in which the time evolutions of both neutrino and matter distributions are computed for only $15$ms from the bounce with the LS EOS. We compare the results so obtained with the original ones to see to what extent the angular resolution could affect the outcome. Note, however, that the comparisons are not so clear-cut as in the previous paper, since the matter dynamics in this phase is chaotic and small perturbations induced by the change in the angular resolution modify not only the neutrino distributions but also the matter configurations in the background substantially.

Richers et al. (2017) demonstrated that our Boltzmann solver tends to underestimate the forward peak in the angular distributions of neutrinos in momentum space at large radii if the number of the angular mesh points is not large enough. This is actually just as expected and was indeed pointed out by Yamada et al. (1999) in their 1D study. As a matter of fact, neutrinos are moving almost radially at large distances from the neutrino sphere no matter what happens to them at small radii and if the angular spread becomes smaller than the smallest width of the angular bin employed in the Boltzmann solver, it is no longer resolved.

Such properties of our Boltzmann solver should have some implications for the success or failure of explosion in our simulations, since the underestimation of the forward peak in the angular distribution in momentum space leads in turn to the overestimation of the local number density of neutrinos and, as a result, the overestimation of neutrino heating in the gain region. On the other hand, Richers et al. (2017) also found that the finite energy resolution tends to underestimate the neutrino heating. We then surmise from these results that the volume-integrated net energy deposition in the gain region is probably underestimated in the current simulations by a few percent.

For the study of the resolution dependence in the early post-bounce phase, we conduct a high-resolution simulation for a short period as mentioned earlier. This time the matter distribution is not fixed but calculated just as in the ordinary run. We deploy $14(\bar{\theta})\times 10(\bar{\phi})$ angular grid points over the entire solid angle while space and energy grids are unchanged from the normal run. In Fig. 7, we compare the radial profiles of some angle-averaged quantities at $15$ms after bounce between the models with the normal and high angular resolutions. As can be seen in this figure, the prompt shock wave is a bit faster and reaches a larger radius in the high-resolution model than in the normal-resolution model (upper left panel); in association with this, the deleptonization behind the shock is slightly stronger in the former around $20\leq r\leq 40{\rm km}$ (upper right). These are all attributed to the fact that the high-resolution simulation experiences a stronger prompt convection. This is indeed corroborated both in the fluid velocity and their lateral component in the convectively unstable region: they are a little larger in the high-resolution simulation consistently. As mentioned earlier, however, matter motions in this region are stochastic due to the chaotic nature of convection. The results would be different substantially if, for example, the initial time is changed even slightly. It is also difficult to isolate the influence of the angular resolution on the neutrino transport alone. More detailed resolution studies in dynamical settings will be reported elsewhere. With these caveats in mind, we will further compare some quantities of relevance in neutrino transport.

Figure 8 displays the radial profiles along two radial rays with $\theta=\pi/4$ (left column) and $\theta=3\pi/4$ (right column) of some relevant quantities in the $\nu_{e}$ distribution at $15$ms after bounce. The neutrino energy is set to the average value at each point. In the top panels, the flux factors ($\bar{F}$) defined in Eq. (7) are shown. One immediately recognizes that it is systematically smaller for the high-resolution case in the post-shock region. This is not directly related with the angular resolution, though. Instead it is simply because the shock radius is larger in the high-resolution run and, as a result, the flux factor increases more slowly from the optically thick limit ($\bar{F}=0$) to the thin limit ($\bar{F}=1$). On the other hand, the flux factor is always smaller for the normal case than for the high-resolution case at large radii. This is a direct resolution effect, i.e., the low-resolution simulation fails to reproduce the forward peak in the angular distribution at large radii.

The rr components of the Eddington tensor, $k^{rr}$, are shown in the middle panels of Fig. 8. It is observed that they also increase a bit more slowly initially in the high-resolution run. This is again a mere consequence of the larger shock radius in that case. In these panels, we also display as additional dotted lines the same components of the Eddington tensor that are obtained with the M1 prescription. Except in the inner optically thick region, they are always slightly greater than those obtained with the Boltzmann code for both resolutions. Considering the result in Richers et al. (2017) that low-resolution computations with the Boltzmann solver tend to underestimate $k^{rr}$, one may think that the results of the M1 prescription is closer to the true values. It should be noted, however, that the differences found here in $k^{rr}$ between the Boltzmann and M1 results are larger than those obtained in Richers et al. (2017) (see Fig.17 in their paper). This may imply that the M1 prescription has its own problem in reproducing $k^{rr}$ for highly-time dependent and highly-inhomogeneous matter distributions considered here. This issue will be further studied in our forthcoming paper. It is incidentally pointed out that the M1 prescription needs the flux factor to obtain the Eddington tensor (see Eqs. (5) and (6)). In the present comparison it is provided by the Boltzmann solver although it should be calculated on its own in the actual simulations with the M1 approximation. It is hence desirable to make comparisons, employing the results of such M1 simulations, which is another subject worth further investigations.

The bottom panels in Fig. 8 are again the Eddington tensors but for the $r\theta$ component $k^{r\theta}$ this time. It should be noted first that $k^{r\theta}$ is very sensitive to the matter motion in the background. As a result, their profiles are quite different between the normal and high-resolution simulations and it is rather difficult to discuss the convergence in the current dynamical setting. Nevertheless, it is evident that the Boltzmann and M1 results are substantially different from each other even qualitatively in the semi-transparent region although they agree in both the optically thin and thick limits irrespective of resolutions. This is indeed consistent with the findings by Richers et al. (2017), who also came to the same conclusion that the difference in $k^{r\theta}$ between the Boltzmann transport with multi-angles and the M1 prescription in the semi-transparent regime is intrinsic and never reduced by increasing resolution. As will be demonstrated in Sec. 6, inaccurate $k^{r\theta}$ may give a $\sim 10\%$ level of errors in the neutrino luminosity and, more importantly, will lead to qualitatively wrong lateral fluxes of neutrinos in the semi-transparent region.

In Fig. 9, we compare the angular distributions in momentum space obtained with the two resolution. Note that, the isotropic contributions are subtracted as previously in these pictures so that the anisotropies could be better recognized. In panel (a), the purple surface is identical to the one presented in Fig. 5, while the black surface is the high-resolution counterpart. In Fig. 9(b), we change the viewing angle to facilitate readers’ understanding of the non-axisymmetric features. As mentioned above, since the matter distributions in the background are different between the two cases, the neutrino angular distributions differ qualitatively. It is important, however, that the degree of asymmetry is even more prominent in the high-resolution simulation. This is again consistent with the finding in Richers et al. (2017) that $k^{r\theta}$ tends to be underestimated in low-angular resolution simulations (see the right panel of Fig.15 in their paper).

6 Possible implications of off-diagonal components on supernova dynamics

The existence of the non-axisymmetric features in the angular distributions of neutrinos and the appearance of the non-vanishing off-diagonal components of the Eddington tensor as a result are the main novel findings in this paper. The legitimate question then is how significant they are for supernova dynamics. In order to fully address this issue, it is required to run additional simulations with some approximate neutrino transport scheme such as the ray-by-ray and/or M1 methods, which either completely ignore or employ a makeshift prescription for these non-axisymmetric features, for the same progenitor, resolution, EOS and input physics and make a detailed comparison, which is certainly beyond the scope of this paper. Instead, in this section, we compare different components of the Eddington tensor quantitatively and discuss how the off-diagonal components might become important.

Note first that the equations for both the zeroth and first moments of the angular distribution include in principle all components of Eddington tensor (see, e.g., Eqs. (3.37) and (3.38) in Shibata et al. (2011)). It should be also pointed out that reaction rates of some neutrino-matter interactions such as non-isoenergetic scatterings and pair processes depend on higher-order moments including the Eddington tensor. The neglect of them may have some implications for CCSNe dynamics. Although this is an interesting issue and is in fact on our to-do-list, in the following, we will limit our discussion to the advection part of the neutrino transport.

The principal part of the equations for the first angular moment or the flux can be approximately written as (see also Eq. (3.38) in Shibata et al. (2011))

[TABLE]

where we ignore collision terms and assume that the spacetime is flat and the background matter is axisymmetric and non-rotating. The off-diagonal component of Eddington tensor $k^{r\theta}$ appears in the second and first terms on the right hand side of Eqs. (11) and (12), respectively. Note that it does not show up in the principal part of the zeroth-order equation for the energy density.

In Fig. 10, we display radial profiles of the absolute values of the ratios of $\partial_{\theta}(Ek^{r\theta})/r$ to $\partial_{r}(Ek^{rr})$ (upper panels) and $\partial_{r}(Ek^{r\theta})$ to $\partial_{\theta}(Ek^{\theta\theta})/r$ (lower panels) on two radial rays with $\theta=\pi/4$ and $3\pi/4$ at $15$ms after bounce. In this analysis, we consider electron-type neutrinos alone, and their energy is set to the mean energy at each point. The results of both the normal- and high-resolution simulations are presented for comparison.

As seen in the upper panels, the radial flux is in general dictated mainly by $k^{rr}$ with $k^{r\theta}$ being at most $10\%$. This is certainly not a large value but still may not be ignored, since, as Burrows et al. (2016) claims, an accumulation of seemingly minor effects may turn out to be crucially important. On the other hand, $k^{r\theta}$ plays more important roles in the equation for the lateral component of neutrino flux as demonstrated in the bottom panels. In fact, the ratio of radial gradient of $Ek^{r\theta}$ to the lateral gradient of $Ek^{r\theta}/r$ exceeds unity in some post-shock regions. This is also the case for the result of the high-resolution simulation although the radial profiles themselves are quite different from those in the normal-resolution run, which is a consequence of the fact that matter distributions in the background become different between the two cases.

As discussed in sections 4 and 5, the M1 prescription is not very successful in reproducing $k^{r\theta}$ in the semi-transparent region particularly in non-spherical settings. Although there is no artificially preferred direction in the M1 transport unlike in the ray-by-ray approximation, the lateral neutrino flux may be still inaccurate. It is misleading to argue that the Eddington tensor is reproduced again very well in the transparent regime with its off-diagonal component becoming negligible compared with the dominant $k^{rr}$. This is because the errors in the semi-transparent region will not be confined there and spread to transparent region in time. The errors in the flux will lead to those in the Eddington tensor through the closure relation, which will again contribute to errors in the flux. This may eventually affect CCSNe dynamics. The quantitative assessment of this effect requires detailed comparisons in collaboration with other groups and is much beyond the scope of this first report of our new simulations.

It is finally mentioned that the above analysis is based on the result of the early post-bounce phase, in which the semi-transparent region is highly dynamical owing to the prompt convection, and the $k^{r\theta}$ effect may be much smaller in the later phase. The errors in early times have some influences on the evolution in later times in principle, though. It should be also added that convections in the proto-neutron star and other hydrodynamical instabilities such as SASI and convections in the heating region occur more often than not even in the late phase. It is repeated that the quantitative assessments are certainly in order and will be studied in subsequent papers.

7 Comparison with previous works

In this section, we attempt to make a comparison of our results with other CCSNe simulations. The same progenitor model has been employed by many authors so far (Müller et al., 2012; Takiwaki et al., 2014; Summa et al., 2016). It is mentioned first that our results are qualitatively in line with them in that softer EOS’s are advantageous for shock revival. It should be pointed out, however, that there are some studies, in which softer EOS including the LS EOS have smaller shock radii initially than the stiffer ones (see, e.g., Fischer et al. (2014)), in apparent contradiction with our results.

According to Fischer et al. (2014), the difference in the shock trajectory originates mainly not from the stiffness of EOS but from the treatment of electron captures on heavy nuclei: representative heavy nuclei tend to be smaller in the softer LS EOS than in the stiffer STOS EOS, which is essentially the same as our FS EOS except for the single-nucleus approximation in the former, resulting in the greater deleptonization in the LS EOS during the collapse phase; this in turn leads to the smaller inner core and hence the weaker prompt shock wave for the LS EOS. It should be recalled, however, that the electron capture rates employed in our simulation with the LS EOS are the same as those for the simulation with the FS EOS. As a result, the effects just mentioned are not taken into account in our current simulations and the shock trajectories reflect the difference in the stiffness of EOS’s alone.

The treatment of nuclear weak interactions consistent with the EOS employed is important to compute CCSNe dynamics accurately. We stress that the current approximate treatment is meant just for simplicity in models with EOS’s that employ the single-nucleus approximation. We believe that multi-nucleus EOS’s are indispensable for the quantitative study of the nuclear weak interactions mentioned above. Such a study indeed under way (Nagakura et al., 2018) with the multi-nucleus extension by Furusawa et al. (2017) of Togashi’s EOS (Togashi & Takano, 2013), which is based on the variational method for realistic nulcear potentials.

It is also important to point out that the shock expansion in our model looks less energetic than those in other simulations with the same progenitor model (see e.g., Takiwaki et al. (2014)). It is difficult to pin down the cause of the discrepancy, since there are many differences in input physics as well as numerical methods for hydrodynamics and neutrino transport, but the ray-by-ray approximation employed for neutrino transport in their simulations may be one of the main causes of the difference. In fact, Skinner et al. (2016) pointed out that the ray-by-ray approximation tends to artificially facilitate explosion in 2D, enhancing sloshing motions in axisymmetry. A similar concern was also expressed by Sumiyoshi et al. (2015), in which they showed that the asymmetry in the neutrino heating tends to be overestimated in the ray-by-ray approximation. More detailed comparisons in collaborations with other groups are required to substantiate the claim, though.

8 Summary and Discussion

We have presented the first report of spatially axisymmetric CCSNe simulations with the full Boltzmann neutrino transport. We have found both similarities and differences between the two models with two different nuclear EOS’s. On the one hand, the neutrino luminosities and mean energies as well as the post-shock morphologies except the scale are very similar between the two. This seems to be a consequence of the cancellation of the stronger bounce that would be expected in the softer LS EOS by the greater electron captures that produced the smaller inner core in the LS EOS model. On the other hand, the neutrino-heating efficiency and the mass in the gain region are consistently higher for the LS EOS. This seems to be due to more vigorous turbulent motions in the post-shock flow for the LS EOS than for the FS EOS, the fact which results in the greater expansion of the shock wave: it has reached $\sim 700{\rm km}$ by $300{\rm ms}$ after bounce and its maximum radius is still growing.

By virtue of the multi-angle treatment in our simulations, we have found interesting features in the neutrino distribution in momentum space, such as the lack of axisymmetry with respective to the local radial direction and the non-vanishing off-diagonal component of the Eddington tensor. With an aid of our previous analyses in Richers et al. (2017) and an additional high-resolution simulation for the early post-bounce phase, we have estimated that the current simulations may have underestimated the neutrino-heating rate by a few percent owing to rather low angular and energy resolutions in momentum space. The possible effects of the off-diagonal component of the Eddington tensor, $k^{r\theta}$, on neutrino transport have been also discussed quantitatively: it plays a non-negligible role for the time evolutions of neutrino fluxes; it may give a $\sim 10\%$ level of contribution to the neutrino luminosity and, more importantly, can be a dominant factor for the time evolution of lateral flux in the semi-transparent region.

We have found an interesting correlation/anti-correlation in $k^{r\theta}$ between $\nu_{e}$ and $\bar{\nu}_{e}$ depending on the radius. It is related with the lateral fluxes of these neutrinos. It will be interesting to see how well the M1 approximation fares in reproducing these features and hence the lateral fluxes. The close comparison between our Boltzmann solver and other approximate methods possibly in collaboration with other groups will be indispensable to assess critically and quantitatively the significance of the findings in this paper for the CCSNe dynamics. It will also enable us to calibrate and possibly improve the prescriptions, which should be given by hand in approximate transport schemes. This is indeed important practically, since our method is very costly in terms of required numerical resources.

We have made an attempt to compare our results with those obtained by other groups for the same progenitor model. We have found that the general trend that softer EOS’s are favorable for shock revival is also true of our simulations. On the other hand, the continuous shock expansion observed for the softer LS EOS looks less energetic than that found by others. Although this seems to be consistent with the finding by Skinner et al. (2016) that the ray-by-ray approximation in spatial axisymmetry may artificially enhance shock revival, more detailed comparisons are certainly necessary to draw some conclusions.

There are also certainly many other issues remaining to be addressed. The top priority is to make detailed comparisons with other approximate methods to assess the importance of multi-angle treatments for supernova dynamics by possibly collaborating with other groups. We will also proceed to explore other progenitors with different masses. The EOS dependence should be further clarified. Rotation is another concern, since the angular distribution in momentum space is then qualitatively changed: e.g., the principal axis will not be aligned with the radial direction in general and another off-diagonal component, $k^{r\phi}$, will no longer be vanishing. We are currently implementing general relativity in our code to investigate its influences, which are expected to be non-negligible. The angular distributions for different species of neutrinos we obtained in this study are valuable in their own right for e.g. the analysis of collective oscillations of neutrino flavors (Duan et al., 2010; Mirizzi, 2013; Capozzi et al., 2017; Izaguirre et al., 2017), which feed on the differences in the angular distributions among different neutrino species. They are currently being investigated and the results will be reported elsewhere.

H.N. acknowledges to C. D. Ott, S. Richers, L. Roberts, D. Radice, M. Shibata, Y. Sekiguchi, K. Kiuchi and T. Takiwaki for valuable comments and discussions. The numerical computations were performed on K computer, at AICS, FX10 at Information Technology Center of Tokyo University, SR16000 at YITP of Kyoto University, and 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, the XC30 and the general common use computer system at the Center for the Computational Astrophysics, CfCA, the National Astronomical Observatory of Japan. Large-scale storage of numerical data is supported by JLDG constructed over SINET4 of NII. H.N and S.F were supported in part by JSPS Postdoctoral Fellowships for Research Abroad No. 27-348 and 28-472 and H.N was partially supported at Caltech through NSF award No. TCAN AST-1333520. This work was supported by Grant-in-Aid for the Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan (15K05093, 24103006, 24105008, 24740165, 24244036, 25870099, 26104006, 16H03986, 17H06357, 17H06365) and HPCI Strategic Program of Japanese MEXT and K-computer at the RIKEN and Post-K project (Project ID: hp 140211, 150225, 160071, 160211, 170230, 170031, 170304).