Three-Dimensional Boltzmann-Hydro Code for core-collapse in massive stars III. A New method for momentum feedback from neutrino to matter

TL;DR

This paper introduces a novel hybrid neutrino-matter coupling method in multi-dimensional supernova simulations, significantly improving momentum conservation accuracy and enabling direct modeling of proto-neutron star kicks.

Contribution

A new hybrid approach for neutrino-matter coupling that enhances momentum conservation in core-collapse supernova simulations with full Boltzmann neutrino transport.

Findings

Improved momentum conservation accuracy in simulations.

First direct modeling of proto-neutron star kicks.

Effective handling of diffusion limit with discrete ordinate method.

Abstract

We present a new method for neutrino-matter coupling in multi-dimensional radiation-hydrodynamic simulations of core-collapse supernova (CCSN) with the full Boltzmann neutrino transport. The development is motivated by the fact that the accurate conservation of momentum is required for reliable numerical modelings of CCSN dynamics including a recoil of proto-neutron star (PNS). The new method is built on a hybrid approach, in which we use the energy-momentum tensor of neutrinos to compute the momentum feedback from neutrino to matter in the optically thick region while we employ the collision term in the optically thin region. In this method we utilize a general relativistic description of radiation-hydrodynamics with angular moments, which allows us to evaluate the momentum feedback from neutrino to matter without inconsistency with our Boltzmann solver. We demonstrate that the new method substantially improves the accuracy of linear momentum conservation in our CCSN simulations under reasonable angular resolutions in momentum space, alleviating the difficulty in giving the diffusion limit precisely with the discrete ordinate (Sn) method. It is the first-ever demonstration that the PNS kick can be handled directly and properly in multi-dimensional radiation-hydrodynamic simulations with the full Boltzmann neutrino transport.

Equations

$$T_{(\mathrm{rad})\beta ;\alpha }^{\alpha }=G_{\beta },$$

$$T_{(\mathrm{rad})\beta }^{\alpha }\equiv \sum _{i=1}\int f_{(\mathrm{i})}{p}^{\alpha }{p}_{\beta }d{V}_p,$$

$$G_{\beta }\equiv \sum _{i=1}\int p_{\beta }\nu S_{\mathrm{rad}(\mathrm{i})}d{V}_p,.$$

$$\displaystyle(\mbox{\boldmath$e$}_{z})^{\beta}T^{\alpha}_{({\rm rad})\beta;\alpha}=(\mbox{\boldmath$e$}_{z})^{\beta}G_{\beta}.$$

$$\displaystyle\frac{1}{r^{2}{\rm sin}{\theta}}\partial_{\alpha}\biggl{(}r^{2}{\rm sin}{\theta}(\mbox{\boldmath$e$}_{z}\cdot\mbox{\boldmath$T$})^{\alpha}\biggl{)}=\mbox{\boldmath$e$}_{z}\cdot\mbox{\boldmath$G$},$$

$$\displaystyle\frac{1}{r^{2}{\rm sin}{\theta}}\partial_{j}\biggl{(}r^{2}{\rm sin}{\theta}(\mbox{\boldmath$e$}_{z}\cdot\mbox{\boldmath$T$})^{j}\biggl{)}=\mbox{\boldmath$e$}_{z}\cdot\mbox{\boldmath$G$},$$

$$\displaystyle 2\pi R^{2}\int_{0}^{\pi}{\rm sin}{\theta}(\mbox{\boldmath$e$}_{z}\cdot\mbox{\boldmath$T$})^{r}(R,\theta)d\theta=$$

$$\displaystyle 2\pi\int_{0}^{R}\int_{0}^{\pi}r^{2}{\rm sin}{\theta}(\mbox{\boldmath$e$}_{z}\cdot\mbox{\boldmath$G$})(r,\theta)d\theta dr.$$

$$H_i=\frac{1}{3\kappa }{\nabla }_{i}J+\mathcal{O}(\frac{1}{{\kappa }^2}),$$

$$K_{ij}=\frac{1}{3}{\gamma }_{ij}J+\mathcal{O}(\frac{1}{{\kappa }^2}),$$

$${\nabla }_z{K}_{zz}=\kappa H_z+\mathcal{O}(\frac{1}{\kappa }).$$

$$T_{(\mathrm{tot})\beta ;\alpha }^{\alpha }=0,$$

$$\displaystyle\mbox{\boldmath$T$}_{\rm(tot)}=\mbox{\boldmath$T$}_{\rm(rad)}+\mbox{\boldmath$T$}_{\rm(mat)}.$$

$$T_{(\mathrm{mat})\beta ;\alpha }^{\alpha }=-G_{\beta }.$$

$$T_{(\mathrm{mat})\beta ;\alpha }^{\alpha }=-T_{(\mathrm{rad})\beta ;\alpha }^{\alpha }.$$

$$T_{(\mathrm{mat})\beta ;\alpha }^{\alpha }=-\lambda T_{(\mathrm{rad})\beta ;\alpha }^{\alpha }-(1-\lambda )G_{\beta }.$$

$$\lambda ({\rho }_{\mathrm{ave}})=\max \left(\min (1,\frac{{\rho }_{\mathrm{ave}}-{\rho }_L}{{\rho }_H-{\rho }_L}),0\right),$$

$${\gamma }_i^{\beta }{T}_{(\mathrm{mat})\beta ;\alpha }^{\alpha }=-\lambda {\gamma }_i^{\beta }{T}_{(\mathrm{rad})\beta ;\alpha }^{\alpha }-(1-\lambda ){\gamma }_i^{\beta }{G}_{\beta }.$$

$${\gamma }_i^{\beta }{T}_{(\mathrm{rad})\beta ;\alpha }^{\alpha }$$

$$E=T_{(\mathrm{rad})}^{\alpha \beta }{n}_{\alpha }{n}_{\beta },$$

$$F_i=-T_{(\mathrm{rad})}^{\alpha \beta }{n}_{\alpha }{\gamma }_{\beta i},$$

$$P_i^j=T_{(\mathrm{rad})}^{\alpha \beta }{\gamma }_{\alpha i}{\gamma }_{\beta }^j,$$

$${\partial }_t{F}_i=\frac{F_i^{(n+1)}-F_i^{(n)}}{\Delta t},$$

$${\gamma }_i^{\beta }{T}_{(\mathrm{rad})\beta ;\alpha }^{\alpha }$$

$$n^{\beta }{T}_{(\mathrm{mat})\beta ;\alpha }^{\alpha }=-n^{\beta }{G}_{\beta },$$

Full text

supernovae: general—neutrinos—hydrodynamics

1 Introduction

First-principles numerical modeling of core-collapse supernova (CCSN) is an indispensable approach to reveal two unresolved mechanisms: explosion and neutron star (NS) kick. These two issues are intimately related with each other and should be addressed simultaneously. One of the numerical challenges in addressing the issue is a requirement of high accuracy under multi-scale, multi-dimensional (multi-D) and multi-physics circumstances, since small numerical errors could crucially affect the final outcome. In this paper we present a novel method to improve the accuracy of linear momentum conservation in Boltzmann-radiation-hydrodynamics simulations of CCSNe.

We first give a brief overview of our CCSN code and some relevant differences from others. The most important feature in our code is that multi-D Boltzmann equations for neutrino transport are solved by discretizing the full phase space (3 in space and 3 in momentum space). The basic framework of our code was developed in Sumiyoshi & Yamada (2012), which treated the non-relativistic Boltzmann equation for a given matter profile (see also Sumiyoshi et al. (2015)). Nagakura et al. (2014) extended it to a special relativistic treatment, in which relativistic corrections were taken into account to all orders of $v/c$ in a non-conventional way. In the same paper we also coupled the Boltzmann solver with a Newtonian self-gravity hydrodynamics, which employs essentially the same numerical method as in Nagakura et al. (2011) except for special relativistic terms. Nagakura et al. (2017) further developed a moving mesh technique to track proper motions of PNS. The technique is built on a general relativistic (GR) description of radiation-hydrodynamics, in which we use a conservative form of GR Boltzmann equation in Shibata et al. (2014). It is probably a unique approach to implement the moving mesh covariantly. More recently we entered the phase of scientific runs of axisymmetric CCSN simulations on $\sim 10$ PFLOPS K-computer (Nagakura et al., 2018b; Harada et al., 2018).

On the other hand, other state-of-the-art numerical simulations have employed some approximations for neutrino transport (see e.g., (Lentz et al., 2015; Takiwaki et al., 2016; Kuroda et al., 2016; Roberts et al., 2016; Bruenn et al., 2016; Pan et al., 2016; Müller et al., 2017; Just et al., 2018; O’Connor & Couch, 2018; Pan et al., 2018; Glas et al., 2018; Kuroda et al., 2018; Takiwaki & Kotake, 2018; Vartanyan et al., 2019) ). Most of them take advantage of (angular) moment approaches one way or another. Integrating out the angular degrees of freedom in the Boltzmann equation, one obtains an infinite hierarchy of equations for all ranks of angular moments. For practical reasons, it is truncated at a certain level: most of the schemes normally treat moments up to the 0th or 1st order. Importantly, the 2nd angular moment is related with the energy-momentum tensor of neutrinos, and its divergence describes the conservation of energy and momentum in neutrino transport (see e.g., Shibata et al. (2011)). It is hence straightforward for these schemes to satisfy the conservation of energy and momentum simultaneously in simulations111Strictly speaking, the statement depends on individual numerical schemes..

Contrary to the moment method, Boltzmann solvers do not guarantee the conservation of energy and momentum in general. The conservative form of GR Boltzmann equation implemented in our code satisfies only the conservation of number and the conservation of the energy and momentum is violated in general once the equation is finite-differenced. In fact, the difficulty of handling multiple conservation laws (number, energy and momentum) in the Boltzmann neutrino transport has been discussed over the past few decades (see e.g., Mezzacappa & Bruenn (1993); Liebendörfer et al. (2004); Cardall et al. (2013)). Unfortunately, however, no satisfactory solution has been demonstrated in realistic CCSN simulations so far.

We found in one of our recent simulations that violation of the linear momentum conservation in the Boltzmann neutrino transport could manifest itself in hydrodynamics, i.e., artificial acceleration of PNS to unphysically large velocities. In this paper we present a novel approach to reduce the error in the momentum conservation and prevent the artificial acceleration of PNS. We will not employ the moment method but cling to solving the Boltzmann equation and change the treatment of the feedback from neutrino to matter. We will also show that the previous treatment, in which the collision term is directly integrated for the feedback, generates a systematic error in the momentum exchange between neutrino and matter in the optically thick region (see in Sec. 3). We investigate the cause of the error, which, it turns out, is intimately related with the problem of the Boltzmann solver in the diffusion limit. The new method alleviates the shortcoming of the Boltzmann solver and allows us to overcome the problem with computationally feasible angular resolutions in momentum space.

This paper is organized as follows. In Sec. 2 we briefly summarize the issue of the unphysical PNS acceleration that we encountered in one of our axisymmetric CCSN simulations. Then, we investigate the cause of the PNS acceleration by systematically carrying out Boltzmann simulations for a frozen fluid background in Sec. 3. In Sec. 4 we describe the connection between the artificial PNS acceleration and the problem of the Boltzmann solver in the diffusion limit, and then we present the new method in detail in Sec. 5. We examine its validity using a series of axisymmetric CCSN simulations in Sec. 6. We summarize our conclusions in Sec. 7. Throughout this paper, Greek and Latin subscripts denote space-time and space components, respectively. We use the metric signature of $-+++$. Unless otherwise stated, we work in units with $c=G=1$, where $c$ and $G$ are the light speed and gravitational constant, respectively.

2 Summary of the problem

In one of our latest axisymmetric CCSN simulations, we found that a PNS was accelerated to velocities more than $1000{\rm km/s}$ by the end of the simulation ($300$ms after bounce). It should be noted, however, that we did not find such an unphysically large PNS kick in other simulations and this happens only in a special situation. In Sec. 3, we investigate the cause of this anomalous PNS acceleration and also clarify the reason why other models do not have the same issue. As we shall see, the large PNS acceleration turns out to be a numerical artifact due to an error in the momentum feedback from neutrino to matter in the optically thick region. Before going into details of the analysis, we summarize the essential results of the simulation, particularly focusing on the PNS kick.

In the simulation, most parts of the numerical setup are the same as those used in Nagakura et al. (2018b). We employ one of the most realistic nuclear equations-of-state (Togashi & Takano, 2013; Togashi et al., 2017). Neutrino-matter interactions are based on those given in Bruenn (1985) except for the recent improvements (Nagakura et al., 2018a). Note that we also incorporated nonisoenergetic scatterings on electrons and positrons, bremsstrahlung by nucleon collisions and electron-positron pair processes together with their inverse reactions (see e.g., Sumiyoshi et al. (2005); Sumiyoshi & Yamada (2012) and references therein). We employ a 11.2 $M_{\sun}$ progenitor in Woosley et al. (2002). 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 and deploy $384(r)\times 128(\theta)$ grid points. Neutrino energy space is discretized non-uniformly with $20$ energy grid points. The lowest-energy cell covers $0\leq\varepsilon\leq 1{\rm MeV}$ and the rest of the mesh covers $1\leq\varepsilon\leq 300{\rm MeV}$ logarithmically. In our code, the polar and azimuthal angles $(\tilde{\theta},\tilde{\phi})$ in neutrino momentum space are locally measured from the radial direction. The angular space is covered with $10(\tilde{\theta})\times 6(\tilde{\phi})$ grid points over the entire solid angle. The grid structure is described in the appendix of Sumiyoshi & Yamada (2012). In our Boltzmann solver, we distinguish three neutrino species: electron-type neutrinos $\nu_{\rm e}$, electron-type anti-neutrinos $\bar{\nu}_{\rm e}$ and all the others collectively denoted by $\nu_{x}$.

In the simulation, a rapid shock expansion occurs in the northern hemisphere at $\sim 180$ms after bounce and then the north-south asymmetry persists. Figure 1 displays the trajectory of shock radius, in which the shaded region corresponds to the range of shock radius. Then, the PNS receives a linear momentum in the opposite direction to the stronger shock expansion. We display the PNS velocity and position in Fig. 2 as a function of time. Note that the velocity of PNS exceeds $1000{\rm km/s}$ and still keeps rising at the end of our simulation ($t=300{\rm ms}$). This unphysically large acceleration is highly likely an numerical artifact.

Before we move on to the in-depth analysis of the cause of the large PNS acceleration, we would like to emphasize that the moving mesh technique developed in Nagakura et al. (2017) has nothing to do with the current issue and it works quite well in the simulation. Note that many other state-of-the-art simulations treat PNS motions quite pragmatically. For instance, some simulations sphericalize either gravitational field or matter distributions or both (Müller et al., 2012; Takiwaki et al., 2014; Müller & Janka, 2015; Müller, 2015; Lentz et al., 2015; Bruenn et al., 2016; Radice et al., 2017; Vartanyan et al., 2018; Müller et al., 2018; Vartanyan et al., 2019), others excise the PNS region from their computational domain (Hanke et al., 2013; Wongwathanarat et al., 2013; Gessner & Janka, 2018).

Note that Scheck et al. (2006) implemented their own moving mesh technique to treat the proper motions of PNS. However, it was applied only to hydrodynamics but not to neutrino transport; they also excised the interior of PNS from the computational domain, in sharp contrast to our computations, in which neutrino transport and hydrodynamics are treated fully consistently on the moving mesh and the whole PNS is included in the computational domain. Nordhaus et al. (2010) handled the NS kick directly in their radiation hydrodynamic simulations with the multi-group flux limited diffusion (MGFLD) approximation without excising the interior of PNS. However, their treatment of neutrino transport is not self-consistent: they added neutrino luminosities by hand on top of the solution of MGFLD neutrino transport in order to promote shock revival. Brandt et al. (2011) carried out axisymmetric CCSN simulations with their own neutrino transport solver without excising the interior of PNS. They solved Boltzmann neutrino transport in the optically thin regime, meanwhile the MGFLD transport was used in the optically thick side. In addition they did not include velocity dependent terms in their transport solver. These are again sharp contrast to our computations. Note also that strong asymmetric shock expansions and NS proper motions were not observed in their simulations.

In these simulations, kick velocities of PNS have been estimated by post-process calculations, i.e., possible feedbacks from the PNS kick to neutrino transport and hydrodynamics are entirely neglected. It is true that these treatments allow CCSN simulations to avoid various practical problems arising from the PNS motion, but it is necessary to justify the prescriptions. Since our method directly handles the PNS motion in a self-consistent manner, it has a potential to serve as a reference model to validate other treatments once the current issue is addressed properly.

3 Cause of the unphysically large acceleration of PNS

In this section, we investigate the cause of the unphysically large acceleration of PNS in detail. To identify the source of the problem clearly, we run additional simulations, in which we freeze the time evolution of matter and compute only steady neutrino distributions on top of it. As shown below, these simulations allow us to disentangle neutrino transport from the complexity of CCSN dynamics and easily quantify the error of linear momentum conservation in our Boltzmann solver. The detailed analysis finally leads us to the conclusion that the large PNS acceleration is a indeed numerical artifact caused by an error in the momentum exchange between neutrino and matter. We also find that the error is negligibly small if the matter profile inside of PNS is almost spherically symmetric. This is the main reason why other CCSN simulations did not yield such a numerical artifact. In this section, we describe these analyses in detail.

By taking the same strategy in Sumiyoshi et al. (2015); Richers et al. (2017), we compute steady neutrino distribution functions by evolving only the neutrino distributions on top of the fixed matter background in the flat spacetime and neglecting the velocity dependence. The matter profile is taken from a snapshot of our axisymmetric CCSN simulation that yielded the large PNS kick. We select the time $t=200$ms, since a runaway of PNS was initiated around this time (see Fig. 2). In these simulations, we also put another simplification: we do not include electron(positron)-neutrino scatterings to save the computational time. Since the opacity of these non-isoenergetic scatterings is subdominant compared with other reactions, this simplification does not affect any conclusions in this analysis.

At first, we explain a reason why the analysis on steady states is useful. We take advantage of the following fact that the advection and the neutrino-matter interactions are balanced with each other for all conserved quantities. As already explained in Sec. 1, however, grid-based Boltzmann solvers like ours fail to guarantee all the balances simultaneously in general. In fact, since we obtain the neutrino distribution function by solving the conservative form of Boltzmann equation, the balance is satisfied for neutrino numbers but it is not the case for energy and momentum. As we will show below, we can easily quantify the violation of momentum conservation in the steady-state analysis. It also provides us with crucial information to judge whether such errors in the Boltzmann solver account for the large PNS accleration observed in the CCSN simulation.

We first prepare some useful equations to check the violation of linear momentum conservation in the steady-state. We start with the equation for the conservation law of the energy-momentum of neutrinos:

[TABLE]

where the left- and right hand sides are defined as

[TABLE]

In these expressions, $f$ and $p^{\mu}$ denote the distribution function and four momentum of neutrinos, respectively. $S_{\rm{rad}}$ originates from the collision term for neutrino-matter interactions in the Boltzmann equations (see also Eq. (1) in Nagakura et al. (2017)); $dV_{p}(=\nu{\rm sin}\bar{\theta}d\nu d\bar{\theta}d\bar{\phi})$ denotes the invariant volume in the neutrino momentum space. Note that we adopt the spherical coordinates in the neutrino momentum space and $\bar{\theta}$ and $\bar{\phi}$ stand for the polar and azimuthal angles, respectively, while $\nu$ denotes the energy of neutrino. $\nu$, $\bar{\theta}$ and $\bar{\phi}$ are measured in the fluid-rest frame222We neglect the fluid-velocity dependence in all simulations discussed in this section, which implies that there is no difference between the fluid-rest and laboratory frames. Note that our new method can be applied even when the velocity-dependent terms are present. We hence specify the frame in defining the momentum variables.. The subscript ”$\rm{i}$” indicates the neutrino species.

Because of axisymmetry, we only focus on the linear momentum parallel to the symmetry axis (z-axis). We make an inner product of the unit 4-vector aligned with the z coordinate, \mbox{\boldmathe}_{z}, with Eq. (1):

[TABLE]

Since \mbox{\boldmathe}_{z} is a killing vector in the flat spacetime, we can rewrite the equation as:

[TABLE]

which describes the linear momentum conservation of neutrinos with respect to the z-direction. In this expression, we employ the spherical coordinates in real space since it is adopted in our Boltzmann solver. By further imposing the steady-state condition in Eq. (5), it can be written as:

[TABLE]

where $j$ stands for the spatial components: $j=\{r,\theta,\phi\}$.

We evaluate each term in Eq. (6) by using the steady profile of $f$ computed with the Boltzmann solver. We then take the sum of all terms on the left hand side of the equation, which is, hereafter, denoted by ”flux-term” while the right hand side of the equation is referred to as ”collision-term”. Their spatial distributions are displayed in Fig. 3. In this figure, we selected two radii inside the PNS ($r=5$ and $10$km), and show the profile as a function of the polar angle. As displayed in this figure, the two terms are almost identical at both radii, which indicates that the balance is well achieved: indeed, the difference is less than $1\%$ with our standard angular resolution in momentum space. As we will show below, however, the precision is not sufficient to make reliable estimations of the PNS acceleration.

We proceed further to assess directly the error in the linear momentum of PNS. We take the volume integral of Eq. (6) over the spatial range $r\leq R$ to obtain the following relation:

[TABLE]

The left hand side of this equation represents the (radially) cumulative profile of the ”flux-term” while the right hand side is the ”collision-term” counterpart. The axisymmetry is also used in obtaining this expression. Importantly, the dimension in the equation is ”[dyn]”, which means that this equation allows us to directly quantify the error in the linear momentum of neutrinos.

In Fig. 4, the flux- (left hand side of Eq. (7)) and collision- (right hand side of Eq. (7)) terms are displayed as solid and dashed lines, respectively. Red color denotes the standard resolution. Note that we run a higher angular resolution simulation, in which we use $14(\bar{\theta})\times 8(\bar{\phi})$ angular grid points: the result is shown in blue color. The right panel focuses on the radial range $20\leq R\leq 40$km for the same quantities. Note that the radial profiles of the flux term for the two angular resolutions are almost identical: the blue solid line masks the red one. This result tells us that the flux term can be computed with a sufficiently high precision even in the standard angular resolution. On the other hand, some angular-resolution dependence can be clearly seen in the collision term. Importantly, the difference between the flux and collision terms is smaller for the higher angular resolution than for the standard one, which is the evidence that the error of linear momentum is smaller in the higher angular resolution in momentum space as expected.

The angular-resolution study reveals that the collision term contains the dominant error in the linear momentum conservation. Since the neutrino transport and hydrodynamics are coupled through the collision term in our CCSN simulations, the errors may be the cause of the PNS acceleration. Here we quantify to what extent the errors are transferred from neutrino to matter using the result of the steady-state simulations. As shown in the right panel of Fig. 4, the difference between the flux and collision terms at the edge of PNS surface ($\sim 40$km) is an order of $\sim 10^{42}$dyn, which corresponds to the PNS acceleration of $\sim 10^{9}{\rm cm/s}^{2}$. This value is consistent with the PNS motion we found in our CCSN simulation (see Fig. 2). We hence conclude that the strong PNS acceleration is a numerical artifact caused by the errors in the momentum exchange between neutrino and matter.

We performed the same analysis for the axisymmetric steady-state simulation with the spherically symmetric matter profile and confirmed that both the flux and collision terms in Eq. (7) were zero to the machine precision. This implies that the artificial PNS acceleration appears once the matter profile strongly deviates from spherical symmetry, which occurs almost at the same time when the shock on the north side expands rapidly in the simulation that gives rise to the unphysically large acceleration of PNS . It is also important to note that the error may have increased with time due to the fact that the PNS motion itself enhances the asymmetry in the matter distribution. In other CCSN simulations, on the other hand, we found less asymmetric matter profiles in the post-shock flow including the PNS. This is the main reason why we did not encounter the numerical artifact of large PNS acceleration in these models.

The analysis in this section shows that the inaccuracy in the momentum conservation is primarily due to the lack of angular resolution in momentum space. This study also indicates that just increasing the angular resolution in momentum space will suppress the artificial acceleration of PNS (see Fig. 4). However, it is not the best solution at the moment, since the CCSN simulations with the full Boltzmann neutrino transport are numerically expensive even under the current standard resolution. In addition, we need to know how fine an angular resolution is required prior to the beginning of simulations, since it may depend on models.

4 Neutrino transport in the diffusion regime

The steady-state analysis in the previous section reveals that the large PNS acceleration is due to the violation of linear momentum conservation in the optically thick region ($r\lesssim 10$km). We also find that the flux-term defined as the left hand side of Eq. (7) is more accurate than the collision-term for the same angular resolution. We take advantage of this property in the new method (see Sec. 5 for details of the method). Before entering into details of the new method, we describe how the different resolution dependences of the flux- and collision-terms are obtained in our Boltzmann solver.

It should be noted that it is not easy for the angular-grid-based Boltzmann solvers to accurately handle the neutrino transport in the diffusion regime: more specifically, it is difficult to evaluate the neutrino flux (i.e., the 1st moment) precisely (see e.g., Mezzacappa & Bruenn (1993); Liebendörfer et al. (2004)). In our previous paper (Sumiyoshi & Yamada, 2012), we conducted detailed tests of our Boltzmann solver in the diffusion regime, in which we studied 1D, 2D and 3D diffusions of a Gaussian packet in a uniform matter background and found in fact that our Boltzmann solver tends to overestimate the neutrino flux due to the numerical diffusion arising from the finite resolution. As shown below, the different angular-resolution dependences of the flux- and collision-terms are intimately related with this problem of the Boltzmann neutrino transport in the diffusion regime.

In this regime, neutrinos interact with matter very frequently and, as a consequence, their angular distribution is almost isotropic. This fact allows us to express the 1st ($H_{i}$) and 2nd ($K_{ij}$) moments in terms of the 0th moment ($J$) as follows (see, e.g., Shibata et al. (2011)):

[TABLE]

where $\kappa$ and $\gamma_{ij}$ denote the total opacity and the 3-metric for flat space, respectively. From these equations, we can also obtain the following relation:

[TABLE]

It is important to realize that this is the same as the diffusion limit of Eq. (6)333See also Berninger et al. (2012) to derive the diffusion equation from the Boltzmann equation.. By virtue of this correspondence, we can understand the angular-resolution dependences of the flux- and collision-terms as follows.

Eq. (10) indicates that the flux- and collision-terms in Eq. (7) are nothing but the 2nd and 1st moments, respectively, in the diffusion limit. This implies that the collision-term in Eq. (7) is more difficult to evaluate precisely with the angular-grid-based Boltzmann solvers in the diffusion limit, since we derive it mainly from the 1st moment. To the contrary, the 0th and 2nd moments are more accurate, since they are dominated by the isotropic component of the angular distribution in this regime. This is the main reason why the larger angular-resolution dependence is observed for the collision-term than for the flux-term. It should be now clear that the use of the 2nd moment instead of the 1st one for the evaluation of the momentum exchange between neutrino and matter in the diffusion regime is the key to the solution of the large-kick problem.

5 New method

The analyses in Sections 3 and 4 reveal the problem in the original treatment of the momentum feedback from neutrino to matter. The direct integration of the collision term is essentially equivalent to utilizing the 1st moment of neutrino angular distribution in the diffusion regime. The angular-grid-based Boltzmann solvers like ours have difficulties in reproducing the 1st moment accurately, compared with the 0th and 2nd moments. It is important to note that the above problem have been already recognized in some previous papers: for instance, Pan et al. (2016) does not use the collision term but compute the neutrino-pressure gradient instead to evaluate the momentum feedback from neutrino to matter in their simulations. It should be noted that we can not utilize the same formulation to address the current issue, since their prescription is based on the diffusion approximation, which is not compatible with the Boltzmann solver. Moreover, our Boltzmann solver is formulated in the fully general relativistic framework and it is a highly non-trivial issue how to evaluate the momentum feedback without inconsistency. In this section, we present a new method, which satisfies all requirements stated above.

The essence of our new method is as follows. We modify only the treatment of the momentum feedback from neutrino to matter without changing the numerical algorithms to solve the Boltzmann equation, the hydrodynamic equations and the Poisson equation for Newtonian gravity (see Fig. 6 in Nagakura et al. (2014)). In the original treatment, we first compute the momentum exchange from individual neutrino-matter interactions, and then their sum is given to the source term of the Euler equation (see Sec. 4 in Nagakura et al. (2017)444It should be noted that there are several typos in Nagakura et al. (2017). On the right hand side of Eq.(31), $G_{i}$ should be replaced with $\alpha\sqrt{\gamma}G_{i}$. In Eq.(43), $G^{r},G^{\theta}$ and $G^{\phi}$ should be replaced with $G_{r},G_{\theta}$ and $G_{\phi}$, respectively. Finally, the sign on the right hand side of Eq.(44) is negative.). In the new method, on the other hand, we evaluate the momentum exchange by taking an appropriate average of the values obtained in two different ways, one of which is the original treatment and the other is the employment of the energy-momentum tensor of neutrinos (see below for more details).

In the new method, the latter treatment is adopted only in the optically thick region with the original treatment still being used in other regions. There are mainly two reasons why we take this approach. The first one is that, in the optically thin region, the accuracy of the 0th and 2nd moments is not much different from that of the 1st moment. In fact, Richers et al. (2017) showed that these moments computed by our Boltzmann solver with the standard angular resolution deviate from those provided by the Monte Carlo transport (which is capable of evaluating moments with higher precision than our method) by $\sim 10\%$ in the optically thin region. The second reason is that the original treatment is preferable for hydrodynamics in the optically thin region. This is simply because the original treatment guarantees that the feedback to hydrodynamics becomes small in the optically thin region just as expected while the new treatment may generate some artificial momentum exchange from numerical errors. Although the resolution dependence and the effect of momentum feedback on CCSN dynamics should be investigated in detail even under the original treatment to make reliable CCSN models, it is not a main subject in this study and will be conducted elsewhere.

Now we describe the new method in detail. We emphasize that the formulation is consistent with the fully general relativistic treatment of neutrino transport in our code and ensures the accurate momentum feedback from neutrino to matter even in the diffusion regime. We start with the total energy-momentum conservation:

[TABLE]

where the total energy-momentum tensor is decomposed as

[TABLE]

In the above equation, \mbox{\boldmathT}_{\rm(mat)} denotes the energy-momentum tensor of matter. From Eqs. (1), (11) and  (12), we can obtain the following relation:

[TABLE]

On the other hand, we can simply rewrite Eq. (11) as

[TABLE]

Mathematically speaking, Eqs. (13) and (14) are equivalent to each other. It is not true in practice, though, since Eq. (1) is not exactly satisfied in our Boltzmann solver (see in Sec. 1). Hence, in the new method, we mix the two equations by introducing a parameter $\lambda$ $(0\leq\lambda\leq 1)$,

[TABLE]

It is chosen so that $\lambda\rightarrow 1$ in the optically thick limit whereas $\lambda\rightarrow 0$ in the opposite limit. Since the optical depth is correlated with the matter density, we determine $\lambda$ as a function of $\rho_{\rm ave}$, which denotes the angle average of matter density. According to the steady-state analyses in Sec. 3, the main region violating the balance between the flux and collision terms is the interior of PNS with $\rho_{\rm ave}\gtrsim 10^{13}{\rm g/cm}^{3}$. We hence take the following function:

[TABLE]

where we set $\rho_{H}=10^{13}{\rm g/cm}^{3}$ and $\rho_{L}=7\times 10^{12}{\rm g/cm}^{3}$, respectively. Although there may be better choices and it is possible to change $\rho_{H}$ and $\rho_{L}$ depending on the angular resolution, we do not explore them further here. The important thing is that the new method with the above choice significantly improves the accuracy of the momentum conservation and addresses well the issues of the unphysical acceleration of PNS (see Sec. 6).

Using Eq. (15), we rewrite the Euler equation in the new method as:

[TABLE]

where $\gamma^{\beta}_{\hskip 2.84526pt\alpha}$ denotes the projection tensor restricted to the spacelike 3-dimensional hypersurface. Below, we describe a procedure to compute $\gamma^{\beta}_{\hskip 2.84526pti}T^{\alpha}_{({\rm rad})\beta;\alpha}$ (the first term of right hand side of the equation), while we refer the reader to Nagakura et al. (2017) for the procedure computing the 2nd term. It is convenient to use a moment formalism of Shibata et al. (2011), since it is fully general relativistic and is compatible with our Boltzmann solver:

[TABLE]

where $\alpha$, $\beta^{i}$ and $\gamma$ denote the lapse function, the shift vector and the determinant of 3-metric ($\gamma_{ij}$), respectively. Other quantities, $E,F_{i}$ and $P_{i}^{\hskip 1.42262ptj}$ are defined as

[TABLE]

where $n^{\alpha}$ denotes the unit 4-vector normal to the spacelike 3-dimensional hypersurface. We evaluate each term on the right hand side of Eq. (18) in our code as follows.

Suppose we know the distribution function $f$ at $t=t^{(n)}$ and $t^{(n+1)}$, where the subscripts represent the time steps. $F_{i}$ can be obtained at both time steps by using Eq. (20). We then evaluate the time derivative of $F_{i}$ as

[TABLE]

where $F_{i}^{(n)}$ and $F_{i}^{(n+1)}$ denote the values of $F_{i}$ at $t=t^{(n)}$ and $t^{(n+1)}$, respectively; $\Delta t$ is the time difference between the two time steps. Note that $\partial_{t}\gamma$ is zero in our current CCSN simulations, since we assume that the spacetime is flat. In a dynamical spacetime, on the other hand, $\partial_{t}\gamma$ will be given by solving the Einstein equation possibly under the 3+1 formalism, which we will couple to our CCSN code in the future. The lapse function $\alpha$ is also set to 1 in the current CCSN code. Note that the shift vector $\beta^{i}$ is not zero and time-dependent quantity in our code since we utilize it to deal with the motion of PNS (see Nagakura et al. (2017)). Only the value of $\beta^{i}$ at $t=t^{(n)}$ is used to evaluate $\beta$-related terms on the right hand side of Eq. (18). The values of $E,F_{i}$ and $P_{i}^{\hskip 1.42262ptj}$ in the same equation are evaluated at $t=t^{(n+1)}$, on the other hand. All spatial derivatives are evaluated as the central difference.

We expect that the above formulation will improve not only the exchange of linear momentum between neutrino and matter but also the total momentum conservation in the radiation-hydrodynamic simulations. This is because our new method is essentially reduced to solving the total momentum conservation equation in the optically thick region, where $\lambda$ is set to 1 in Eq. (17). In addition, the dominant component on the right hand side of Eq. (18) is the neutrino pressure gradient, i.e.,

[TABLE]

in the optically thick region, as long as the motion of PNS is slow, $V_{\rm NS}/c\ll 1$, where $V_{\rm NS}$ denotes the speed of the proper motion of PNS with respect to the laboratory frame. Recall that the neutrino pressure is evaluated much more accurately than the flux in the diffusion regime with our code. Note also that Eq. (23) is essentially the same as the treatment of momentum feedback in Pan et al. (2016). This implies that our new formulation is an extension from the previous treatment that reproduces the diffusion limit but incorporates the all other effects such as deviations from the limit as well as proper motions of PNS, which could not be handled by the previous prescription, without losing the consistency with the general relativistic Boltzmann neutrino transport.

Note also that our method is fundamentally different from the two-moment method, since we evaluate $P_{i}^{\hskip 1.42262ptj}$ based on the solution of Boltzmann equation555Note that the closure relation in the optically thick region is well established in the two-moment method. Thus, the difference of $P_{i}^{\hskip 1.42262ptj}$ between our method and the two-moment method may be small.. On the other hand, our method does not perfectly guarantee the conservation of momentum since we keep using a original Boltzmann solver in the optically thin region. Note that if we used the prescription for the optically thick region also in the optically thin region, it would enforce the multiple conservation laws simultaneously. However, as mentioned already, that would sacrifice the accuracy and may generate different problems from the lack of precision in the moments in the optically thin region.

Finally we make a comment on the energy feedback from neutrino to matter. Since the energy equation has nothing to do with the artificial PNS acceleration, we leave untouched the original treatment in the energy equation:

[TABLE]

i.e., we compute the feedback by directly integrating the neutrino-matter interactions.

6 Examinations

In this section we examine the validity of our new method by running some axisymmetric CCSN simulations. As shown below, we find that the new method successfully suppresses the artificial acceleration of PNS. We also show that the new method improves the accuracy of the total linear momentum conservation in our CCSN code.

As the first examination, we run four short-term ($8$ms) simulations. We name them ”Old-ST”, ”Old-HR”, ”New-ST” and ”New-HR”, respectively. The term of ”Old” denotes the simulations under the previous treatment, while ”New” is given to the models employing the new method. ”ST” and ”HR” represent the standard and high angular resolutions, respectively, in the former of which $10(\tilde{\theta})\times 6(\tilde{\phi})$ and in the latter $14(\tilde{\theta})\times 8(\tilde{\phi})$ grid points are deployed, respectively. We prepare the initial condition by taking the matter and neutrino profiles at $t=200$ms in the original simulation. In the ”HR” models, we set the initial angular distribution of neutrinos by linearly interpolating the original data. In these simulations, we restore all input physics, velocity-dependent terms and the moving-mesh technique.

Fig. 5 displays the time evolution of the velocity of PNS ($V_{\rm PNS}$) for all models. Different colors distinguish the previous and new treatments while the line types represent the angular resolutions. As clearly seen in this figure, our new method (blue lines) suppresses the PNS acceleration. We also find that the difference of $V_{\rm PNS}$ between New-ST and New-HR is much smaller than that in ”Old” models. This indicates that the accuracy of momentum conservation is improved by our new method. Note that the result in Old-HR is close to those in ”New” models. This fact implies that even the previous prescription will obtain the convergence to the true result by increasing angular resolutions. Another important finding is that New-ST may decelerate the PNS a bit excessively, since the result of New-HR comes between those of Old-HR and New-ST. This may indicate that the new treatment with the standard angular resolution would underestimate the velocity of PNS. This fact should be kept in mind in our future CCSN simulations.

Finally we carry out a long-term simulation (up to $t=300$ms) with the new method. The simulation is started at $t=150$ms from the result of the original simulation. We choose this time because it is slightly before the start of the runaway motion of PNS in the original computation. Fig. 6 displays the time evolution of the position and $V_{\rm PNS}$ of PNS. We also show the same quantities in the original simulation as red lines. As can be seen in this figure, our new treatment works quite well and no numerical problems are observed in this long-term simulation.

7 Conclusions

We develop a new method for the neutrino-matter coupling to improve the accuracy of the linear momentum conservation in multi-D CCSN simulations with the angular-grid-based Boltzmann solver for neutrino transport. This effort is motivated by the requirement of extremely high accuracy of CCSN simulations for reliable modelling of CCSN explosions and PNS kicks under restricted computational resources. In fact, we observed an unphysically large acceleration of PNS in one of our axisymmetric CCSN simulations and suspected that it was due to the violation of linear momentum conservation in our Boltzmann solver.

To see the problem clearly, we perform Boltzmann simulations on top of the fixed matter profile. This study reveals that the large acceleration of PNS is due to errors in the momentum feedback from neutrino to matter in the optically thick region. Note that our Boltzmann solver works basically well (see Fig.3): indeed, it satisfies the linear momentum conservation within $\lesssim 1\%$ for the standard angular resolution in momentum space. However, this accuracy is not sufficient to handle the PNS motion when the matter profile becomes strongly asymmetric. We find that the error of total linear momentum leads to the PNS acceleration of $\sim 10^{9}{\rm cm/s}^{2}$, which accounts indeed for the large acceleration of PNS observed in one of our CCSN simulations (see Fig. 4). We also find that increasing the angular resolution in momentum space alone reduces the PNS acceleration, which indicates that the inaccuracy in the momentum conservation is due to the lack of angular resolution in momentum space for neutrinos. This means that we can avoid the problem simply by increasing the angular resolution. However, it is not the best solution, since it is inhibitingly costly in terms of the numerical resource. Thus, the development of the new method for the improvement of the neutrino-matter coupling is one of necessary tasks in our CCSN project.

The basic idea in the new method is that our Boltzmann solver accurately evaluates the momentum feedback from neutrino to matter in the optically thick region with the currently available angular resolution in momentum space if the energy-momentum tensor of neutrinos is employed. This is due to the fact that the energy momentum tensor is dominated by the isotropic component of the angular distribution of neutrino. We also identify the source of errors in the previous treatment, in which the collision term is directly integrated. The problem is that it essentially utilizes the 1st moment of the angular distribution in momentum space to evaluate the momentum feedback from neutrino to matter. It should be noted that the angular-grid-based Boltzmann solvers tend to generate larger relative errors in the 1st moment than in the 0th and 2nd moments, since the former depends solely on the deviation from isotropy, which is very small in the diffusion regime; on the other hand, the latter two are mainly determined by the dominant isotropic component and can be obtained accurately in the same regime. This is the main reason why the original treatment has higher sensitivity to the angular resolution in momentum space.

Although the problem of the Boltzmann neutrino transport in the diffusion regime has been already recognized and possible solutions were proposed for some numerical schemes (see e.g., Pan et al. (2016)), they took full advantage of the diffusion approximation, which is not compatible with our Boltzmann solver. Moreover, we need to ensure that the solution should be consistent with the general relativistic formulation of our code. This is not a trivial issue and needs the new method, which satisfies all these requirements.

In the new method, the momentum feedback from neutrino to matter is computed based on the temporal and spatial derivatives of energy-momentum tensor (see Eq. (18)). Although the new treatment can be applied in all regimes in principle, we adopt it only in the optically thick region. In the optically thin region, we continue to use the original treatment, in which we compute the momentum exchange between neutrino and matter by integrating the collision terms for individual neutrino-matter interactions and then give their sum to the Euler equation. In fact, the original method is more suitable in the optically thin region, since the momentum exchange decreases just as expected as the neutrino-matter interaction becomes inefficient, while it is not necessarily guaranteed in the new treatment.

We find that the new method works quite well. It successfully suppresses the unphysically large PNS acceleration even with the standard angular resolution in momentum space (see Fig. (6)). The angular-resolution study also shows that the difference between the standard- and high-resolution simulations is smaller in the new method than in the original treatment (see Fig. (5)).

Importantly, we still have a coherent PNS motion in the same simulation but with the new treatment (see Fig. 6). We believe that this recoil of PNS is real, i.e., this is the first successful direct treatment of PNS kick in multi-D radiation-hydrodynamic simulations with the full Boltzmann neutrino transport. The impact on CCSN dynamics will be discussed in a separate paper. Last but not least, the new method will be useful to improve the accuracy of the momentum exchange between radiation and matter for other radiation-hydrodynamic codes that adopt multi-angle treatments such as those by Livne et al. (2004); Ott et al. (2008) for neutrinos and by Jiang et al. (2014a, b); Roth & Kasen (2015); Hopkins & Grudić (2018) for photons. The versatility is one of the important merits in the new method.

We acknowledge Adam Burrows for fruitful discussion. The numerical computations were performed on the supercomputers at K, at AICS, FX10 at Information Technology Center of Nagoya University. Large-scale storage of numerical data is supported by JLDG constructed over SINET4 of NII. H.N. was supported by Princeton University through DOE SciDAC4 Grant DE-SC0018297 (subaward 00009650). This work was also supported by Grant-in-Aid for the Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan (15K05093, 25870099, 26104006, 16H03986, 17H06357, 17H06365), HPCI Strategic Program of Japanese MEXT and K computer at the RIKEN (Project ID: hpci 160071, 160211, 170230, 170031, 170304, hp180179, hp180111, hp180239).