A NEW MULTI-ENERGY NEUTRINO RADIATION-HYDRODYNAMICS CODE IN FULL GENERAL RELATIVITY AND ITS APPLICATION TO THE GRAVITATIONAL COLLAPSE OF MASSIVE STARS

There are major differences compared to our previous code (Kuroda et al. 2012). On the one hand, we evolved an energy-integrated ("gray") neutrino radiation field with an approximate description of neutrino–matter interaction based on the neutrino leakage scheme, and on the other hand we now solve the spectral neutrino transport where the source terms are treated self-consistently following a standard procedure of the M1 closure scheme (Shibata et al. 2011). For convenience, we briefly summarize the basic equations of our newly developed code in the following (see Shibata et al. 2011 and Cardall et al. 2013a for the complete derivation).

The total stress–energy tensor ${T}_{(\mathrm{total})}^{\alpha \beta }$ is expressed as

$$\begin{eqnarray}&&{T}_{(\mathrm{total})}^{\alpha \beta }={T}_{(\mathrm{fluid})}^{\alpha \beta }+\int d\varepsilon \displaystyle \sum _{\nu \in {\nu }_{e},{\bar{\nu }}_{e},{\nu }_{x}}{T}_{(\nu ,\varepsilon )}^{\alpha \beta },\end{eqnarray} \tag{ 1 }$$

where ${T}_{(\mathrm{fluid})}^{\alpha \beta }$ and ${T}_{(\nu ,\varepsilon )}^{\alpha \beta }$ are the stress–energy tensor of fluid and the energy-dependent neutrino radiation field, respectively. Note in the above equation, summation is taken for all species of neutrinos (${\nu }_{e},{\bar{\nu }}_{e},{\nu }_{x}$) with ${\nu }_{x}$ representing heavy lepton neutrinos (i.e., ${\nu }_{\mu },{\nu }_{\tau }$ and their anti-particles), and ε represents neutrino energy measured in the comoving frame with the fluid. For simplicity, the neutrino flavor index ν is omitted below.

With introducing radiation energy (${E}_{(\varepsilon )}$), radiation flux $({F}_{(\varepsilon )}^{\mu })$, and radiation pressure $({P}_{(\varepsilon )}^{\mu \nu })$, measured by an Eulerian observer or (${J}_{(\varepsilon )}$, ${H}_{(\varepsilon )}^{\mu }$ and ${L}_{(\varepsilon )}^{\mu \nu }$) measured in a comoving frame, ${T}_{(\varepsilon )}^{\mu \nu }$ can be written in covariant form as

$$\begin{eqnarray}&&{T}_{(\varepsilon )}^{\mu \nu }={E}_{(\varepsilon )}{n}^{\mu }{n}^{\nu }+{F}_{(\varepsilon )}^{\mu }{n}^{\nu }+{F}_{(\varepsilon )}^{\nu }{n}^{\mu }+{P}_{(\varepsilon )}^{\mu \nu },\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&={J}_{(\varepsilon )}{u}^{\mu }{u}^{\nu }+{H}_{(\varepsilon )}^{\mu }{u}^{\nu }+{H}_{(\varepsilon )}^{\nu }{u}^{\mu }+{L}_{(\varepsilon )}^{\mu \nu }.\end{eqnarray} \tag{ 3 }$$

In the above equations, ${n}^{\mu }=(1/\alpha ,-{\beta }^{k}/\alpha )$ is a unit vector orthogonal to the space-like hypersurface and ${u}^{\mu }$ is the four velocity of fluid. In the truncated moment formalism (Thorne 1981; Shibata et al. 2011), one evolves radiation energy (${E}_{(\varepsilon )}$) and radiation flux (${F}_{(\varepsilon )}^{\alpha }$) in a conservative form and ${P}_{(\varepsilon )}^{\mu \nu }$ is determined by an analytic closure relation (e.g., Equation (6)). The evolution equations for ${E}_{(\varepsilon )}$ and ${F}_{(\varepsilon )}^{\alpha }$ are given by

$$\begin{eqnarray}&&{\partial }_{t}\sqrt{\gamma }{E}_{(\varepsilon )}+{\partial }_{i}\sqrt{\gamma }(\alpha {F}_{(\varepsilon )}^{i}-{\beta }^{i}{E}_{(\varepsilon )})+\sqrt{\gamma }\alpha {\partial }_{\varepsilon }(\varepsilon {\tilde{M}}_{(\varepsilon )}^{\mu }{n}_{\mu })\\ &&\quad =\sqrt{\gamma }(\alpha {P}_{(\varepsilon )}^{{ij}}{K}_{{ij}}-{F}_{(\varepsilon )}^{i}{\partial }_{i}\alpha -\alpha {S}_{(\varepsilon )}^{\mu }{n}_{\mu }),\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&{\partial }_{t}\sqrt{\gamma }{{F}_{(\varepsilon )}}_{i}+{\partial }_{j}\sqrt{\gamma }(\alpha {{P}_{(\varepsilon )}}_{i}^{j}-{\beta }^{j}{{F}_{(\varepsilon )}}_{i})-\sqrt{\gamma }\alpha {\partial }_{\varepsilon }(\varepsilon {\tilde{M}}_{(\varepsilon )}^{\mu }{\gamma }_{i\mu })\\ &&\quad =\sqrt{\gamma }[-{E}_{(\varepsilon )}{\partial }_{i}\alpha +{{F}_{(\varepsilon )}}_{j}{\partial }_{i}{\beta }^{j}+(\alpha /2){P}_{(\varepsilon )}^{{jk}}{\partial }_{i}{\gamma }_{{jk}}+\alpha {S}_{(\varepsilon )}^{\mu }{\gamma }_{i\mu }],\end{eqnarray} \tag{ 5 }$$

respectively. Here γ is the determinant of the three metric $\gamma \equiv {\rm{det}}({\gamma }_{{ij}})$ and ${S}_{(\varepsilon )}^{\mu }$ is the source term for neutrino matter interactions (see Appendix A for the currently implemented processes). ${\tilde{M}}_{(\varepsilon )}^{\mu }$ is defined by ${\tilde{M}}_{(\varepsilon )}^{\mu }\equiv {M}_{(\varepsilon )}^{\mu \alpha \beta }{{\rm{\nabla }}}_{\beta }{u}_{\alpha }$, where ${M}_{(\varepsilon )}^{\mu \alpha \beta }$ denotes the third rank moment of the neutrino distribution function (see Shibata et al. 2011 for the explicit expression).

By adopting the M1 closure scheme, the radiation pressure can be expressed as

$$\begin{eqnarray}&&{P}_{(\varepsilon )}^{{ij}}=\displaystyle \frac{3{\chi }_{(\varepsilon )}-1}{2}{P}_{\mathrm{thin}(\varepsilon )}^{{ij}}+\displaystyle \frac{3(1-{\chi }_{(\varepsilon )})}{2}{P}_{\mathrm{thick}(\varepsilon )}^{{ij}},\end{eqnarray} \tag{ 6 }$$

where ${\chi }_{(\varepsilon )}$ represents the variable Eddington factor, and ${P}_{\mathrm{thin}(\varepsilon )}^{{ij}}$ and ${P}_{\mathrm{thick}(\varepsilon )}^{{ij}}$ correspond to the radiation pressure in the optically thin and thick limit, respectively. They are written in terms of ${J}_{(\varepsilon )}$ and ${H}_{(\varepsilon )}^{\mu }$ (Shibata et al. 2011). Following Minerbo (1978), Cernohorsky & Bludman (1994), and Obergaulinger & Janka (2011), we take the variable Eddington factor ${\chi }_{(\varepsilon )}$ as

$$\begin{eqnarray}&&{\chi }_{(\varepsilon )}=\displaystyle \frac{5+6{\bar{F}}_{(\varepsilon )}^{2}-2{\bar{F}}_{(\varepsilon )}^{3}+6{\bar{F}}_{(\varepsilon )}^{4}}{15},\end{eqnarray} \tag{ 7 }$$

where

$$\begin{eqnarray}&&{\bar{F}}_{(\varepsilon )}^{2}\equiv \displaystyle \frac{{h}_{\mu \nu }{H}_{(\varepsilon )}^{\mu }{H}_{(\varepsilon )}^{\nu }}{{J}_{(\varepsilon )}^{2}}.\end{eqnarray} \tag{ 8 }$$

In Equation (8), ${h}_{\mu \nu }\equiv {g}_{\mu \nu }+{u}_{\mu }{u}_{\nu }$ is the projection operator. As we will discuss later, by the definition of ${\bar{F}}_{(\varepsilon )}$ in Equation (8), one can appropriately reproduce several important neutrino behaviors, for example, neutrino trapping in the rapidly collapsing opaque core. We iteratively solve the simultaneous Equations (7)–(8) to find the converged solution of ${\chi }_{(\varepsilon )}$.

The hydrodynamic equations are written in a conservative form as

$$\begin{eqnarray}&&{\partial }_{t}{\rho }_{*}+{\partial }_{i}({\rho }_{*}{v}^{i})=0,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\partial }_{t}\sqrt{\gamma }{S}_{i}+{\partial }_{j}\sqrt{\gamma }({S}_{i}{v}^{j}+\alpha P{\delta }_{i}^{j})\\ &&\quad =-\sqrt{\gamma }[{S}_{0}{\partial }_{i}\alpha -{S}_{k}{\partial }_{i}{\beta }^{k}-2\alpha {S}_{k}^{k}{\partial }_{i}\phi \\ &&\qquad +\left.\alpha {e}^{-4\phi }({S}_{{jk}}-P{\gamma }_{{jk}}){\partial }_{i}{\tilde{\gamma }}^{{jk}}/2+\alpha \displaystyle \int d\varepsilon {S}_{(\varepsilon )}^{\mu }{\gamma }_{i\mu }\right],\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\partial }_{t}\sqrt{\gamma }\tau +{\partial }_{i}\sqrt{\gamma }(\tau {v}^{i}+P({v}^{i}+{\beta }^{i}))\\ &&\quad =\sqrt{\gamma }[\alpha {{KS}}_{k}^{k}/3+\alpha {e}^{-4\phi }({S}_{{ij}}-P{\gamma }_{{ij}})\tilde{{A}^{{ij}}}\\ &&\qquad -\left.{S}_{i}{D}^{i}\alpha +\alpha \displaystyle \int d\varepsilon {S}_{(\varepsilon )}^{\mu }{n}_{\mu }\right],\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{\partial }_{t}({\rho }_{*}{Y}_{e})+{\partial }_{i}({\rho }_{*}{Y}_{e}{v}^{i})=\sqrt{\gamma }\alpha {m}_{{\rm{u}}}\int \displaystyle \frac{d\varepsilon }{\varepsilon }({S}_{({\nu }_{e},\varepsilon )}^{\mu }-{S}_{({\bar{\nu }}_{e},\varepsilon )}^{\mu }){u}_{\mu },\end{eqnarray} \tag{ 12 }$$

where ${\rho }_{*}=\rho \sqrt{\gamma }W$, ${S}_{i}=\rho {{hWu}}_{i}$, ${S}_{{ij}}=\rho {{hu}}_{i}{u}_{j}+P{\gamma }_{{ij}}$, ${S}_{k}^{k}={\gamma }^{{ij}}{S}_{{ij}}$, ${S}_{0}=\rho {{hW}}^{2}-P$, and $\phi ={\rm{log}}(\gamma )/12$. ρ is the rest mass density, W is the Lorentz factor, $h=1+e+P/\rho$ is the specific enthalpy, ${v}^{i}={u}^{i}/{u}^{t}$, $\tau ={S}_{0}-\rho W$, ${Y}_{e}\equiv {n}_{e}/{n}_{b}$ is the electron fraction (n_X is the number density of X), e and P are the specific internal energy and pressure of matter, respectively, and ${m}_{{\rm{u}}}$ is the atomic mass unit. $P(\rho ,s,{Y}_{e})$ and $e(\rho ,s,{Y}_{e})$ are given by an equation of state (EOS) with s denoting the entropy per baryon. We employ an EOS by Lattimer & Douglas Swesty (1991, hereafter LS220; see Section 5.1 for more details). In the right-hand side (rhs) of Equation (11), Dⁱ represents the covariant derivative with respect to the three metric ${\gamma }_{{ij}}$.

As explained in Section 2.1, the formalism of our code that treats the radiation-hydrodynamics equations in a conservative form is suitable to satisfy the energy conservation of the total system (neutrinos and matters; see also Kuroda & Umeda 2010 and Kuroda et al. 2012 for more details). Let us first show how the energy conservation law is obtained in our code.

To focus only on the energy exchange between the matter and neutrino radiation field, we omit the gravitational source term in the following discussion. Then, the equations of energy conservation of matter and neutrinos (e.g., Equations (4) and (11)) become

$$\begin{eqnarray}&&{\partial }_{t}\sqrt{\gamma }\tau +{\partial }_{i}\sqrt{\gamma }(\tau {v}^{i}+P({v}^{i}+{\beta }^{i}))=\int d\varepsilon \sqrt{\gamma }\alpha {S}_{(\varepsilon )}^{\mu }{n}_{\mu },\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\partial }_{t}\sqrt{\gamma }{E}_{(\varepsilon )}+{\partial }_{i}\sqrt{\gamma }(\alpha {F}_{(\varepsilon )}^{i}-{\beta }^{i}{E}_{(\varepsilon )})+\sqrt{\gamma }\alpha {\partial }_{\varepsilon }(\varepsilon {\tilde{M}}_{(\varepsilon )}^{\mu }{n}_{\mu })\\ &&\quad =-\sqrt{\gamma }\alpha {S}_{(\varepsilon )}^{\mu }{n}_{\mu }.\end{eqnarray} \tag{ 14 }$$

From the above two equations, one can readily see that the total energy (sum of matter and neutrinos) contained in the computational domain, ${E}_{\nu m}\equiv \int {{dx}}^{3}\sqrt{\gamma }(\tau +\int d\varepsilon {E}_{(\varepsilon )})$, is conserved in our basic equations as long as there is no net energy flux through the numerical and momentum space boundaries (i.e., $\int d\varepsilon [{\partial }_{\varepsilon }(\varepsilon {\tilde{M}}_{(\varepsilon )}^{\mu }{n}_{\mu })]=0$).

The lepton number conservation needs to be satisfied with good accuracy because it determines the PNS mass and the postbounce supernova dynamics. We here explain how we treat it in our code. As for the electron and neutrino number conservation, the basic equations are given by

$$\begin{eqnarray}&&{\partial }_{t}\left(\displaystyle \frac{{\rho }_{*}{Y}_{e}}{{m}_{{\rm{u}}}}\right)+{\partial }_{i}\left(\displaystyle \frac{{\rho }_{*}{Y}_{e}{v}^{i}}{{m}_{{\rm{u}}}}\right)\\ &&\quad =\sqrt{\gamma }\alpha \displaystyle \int \displaystyle \frac{d\varepsilon }{\varepsilon }({S}_{({\nu }_{e},\varepsilon )}^{\mu }-{S}_{({\bar{\nu }}_{e},\varepsilon )}^{\mu }){u}_{\mu },\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\partial }_{t}({q}_{(\varepsilon )}^{0}\sqrt{\gamma }\alpha )+{\partial }_{i}({q}_{(\varepsilon )}^{i}\sqrt{\gamma }\alpha )-\sqrt{\gamma }\alpha {\partial }_{\varepsilon }(\varepsilon {q}_{(\varepsilon )}^{\alpha \beta }){{\rm{\nabla }}}_{\beta }{u}_{\alpha }\\ &&\quad =-\displaystyle \frac{\sqrt{\gamma }\alpha }{\varepsilon }{S}_{(\varepsilon )}^{\mu }{u}_{\mu },\end{eqnarray} \tag{ 16 }$$

where

$$\begin{eqnarray}&&{q}_{(\varepsilon )}^{\alpha }\equiv -{\varepsilon }^{-1}{T}_{(\varepsilon )}^{\alpha \beta }{u}_{\beta }=({{ \mathcal J }}_{(\varepsilon )}{u}^{\alpha }+{{ \mathcal H }}_{(\varepsilon )}^{\alpha }),\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{q}_{(\varepsilon )}^{\alpha \beta }\equiv {\varepsilon }^{-1}{T}_{(\varepsilon )}^{\gamma \beta }{h}_{\gamma }^{\alpha }=({{ \mathcal H }}_{(\varepsilon )}^{\alpha }{u}^{\beta }+{{ \mathcal L }}_{(\varepsilon )}^{\alpha \beta }),\end{eqnarray} \tag{ 18 }$$

with $({ \mathcal J },{{ \mathcal H }}^{\alpha },{{ \mathcal L }}^{\alpha \beta })\equiv {\varepsilon }^{-1}(J,{H}^{\alpha },{L}^{\alpha \beta })$. The conservation equation for neutrinos (16) corresponds to Equation (3.22) (divided by ε) in Shibata et al. (2011). Since the neutrino number density measured by an Eulerian observer is expressed as

$$\begin{eqnarray}&&{n}_{\nu ,{\rm{Euler}}}=\int d\varepsilon {q}_{(\nu ,\varepsilon )}^{0}\sqrt{\gamma }\alpha ,\end{eqnarray} \tag{ 19 }$$

the equation of the total lepton number conservation becomes

$$\begin{eqnarray}&&{\partial }_{t}({\rho }_{*}{Y}_{l})+{\partial }_{i}\left({\rho }_{*}{Y}_{e}{v}^{i}+{m}_{{\rm{u}}}\sqrt{\gamma }\alpha \int d\varepsilon [{q}_{({\nu }_{e},\varepsilon )}^{i}-{q}_{({\bar{\nu }}_{e},\varepsilon )}^{i}]\right)=0.\end{eqnarray} \tag{ 20 }$$

Here the total lepton fraction Y_l is defined by

$$\begin{eqnarray}{Y}_{l} & = & {Y}_{e}+{Y}_{{\nu }_{e}}-{Y}_{{\bar{\nu }}_{e}}\\ & = & {Y}_{e}+\displaystyle \frac{{m}_{{\rm{u}}}}{\rho {u}^{t}}\displaystyle \int d\varepsilon ({q}_{({\nu }_{e},\varepsilon )}^{0}-{q}_{({\bar{\nu }}_{e},\varepsilon )}^{0}).\end{eqnarray} \tag{ 21 }$$

From Equation (20), one can readily see that the total lepton number is conserved irrespective of the included neutrino matter interaction processes in case there is no net flux through the numerical and energy space boundaries. The distribution of Y_l into the each component (e.g., ${Y}_{e},{Y}_{{\nu }_{e}},{Y}_{{\bar{\nu }}_{e}}$) is determined by the details of the implemented microphysics, which should be checked carefully and will be reported in Section 5.

2.2.1. Neutrino Number Transport in the Diffusion Limit

In the collapsing iron core, it is well known that the central core becomes opaque to neutrinos due mainly to scattering off heavy nuclei when the central density exceeds ∼10^11–12 g cm⁻³ (Sato 1975) and neutrinos are trapped with matter afterward. In the diffusion limit at large neutrino opacities, the trapped neutrinos move with the matter velocity v_i for an Eulerian observer. Thus, their advection equation can be described with the same form of Y_e as

$$\begin{eqnarray}&&{\partial }_{t}{\rho }_{*}{Y}_{\nu }+{\partial }_{i}{\rho }_{*}{Y}_{\nu }{v}^{i}=(\mathrm{Source}\;{\rm{terms}}).\end{eqnarray} \tag{ 22 }$$

Because the source terms in the above equation amount to equal the negative value of the source terms in the electron number conservation equation, the core lepton number is conserved with good accuracy until it gradually decreases by diffusion in the PNS cooling phase. Since the central core mass depends on the core lepton number, the CCSN simulation should be able to capture this important phenomena appropriately.

In our formalism, however, this is not a trivial problem because we solve the energy–momentum conservation in Equations (4)–(5) and not the neutrino number conservation in Equation (16). In this section, we check whether our basic equations can reproduce the neutrino diffusion limit adequately, i.e., they reach asymptotically to Equation (22).

For simplicity, we assume in the following that the spacetime is flat and the typical velocity of the matter field (v) is much smaller than the speed of light (slow motion limit; neglecting terms higher than the second order with respect to $(v/c)$).

Let us first check whether Equation (16) can satisfy the local neutrino number conservation in the trapped region. In this limit, the neutrino number density at each energy bin q⁰ and its flux qⁱ approach

$$\begin{eqnarray}&&{q}^{0}={ \mathcal J }{u}^{t}+{{ \mathcal H }}^{0}\sim { \mathcal J }+\displaystyle \frac{{{ \mathcal H }}_{i}{v}^{i}}{W}\sim { \mathcal J },\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{q}^{i}={ \mathcal J }{u}^{i}+{{ \mathcal H }}^{i}\sim { \mathcal J }{v}^{i}+{{ \mathcal H }}^{i}\sim {q}^{0}({v}^{i}+{{ \mathcal H }}^{i}/{ \mathcal J }).\end{eqnarray} \tag{ 24 }$$

From these relations, it is obvious that the neutrino number density at each energy bin q⁰ is transferred with the matter velocity vⁱ plus the diffusion velocity ${{ \mathcal H }}^{i}/{ \mathcal J }$ and the equation of the total lepton number (Equation (20)) in the slow motion limit becomes

$$\begin{eqnarray}&&{\partial }_{t}({\rho }_{*}{Y}_{l})+{\partial }_{i}\left({\rho }_{*}{Y}_{e}{v}^{i}+{m}_{{\rm{u}}}\sqrt{\gamma }\alpha \displaystyle \int d\varepsilon [{q}_{({\nu }_{e},\varepsilon )}^{0}-{q}_{({\bar{\nu }}_{e},\varepsilon )}^{0}]{v}^{i}\right)\\ &&\quad ={\partial }_{t}({\rho }_{*}{Y}_{l})+{\partial }_{i}({\rho }_{*}{Y}_{l}{v}^{i})=0,\end{eqnarray} \tag{ 25 }$$

demonstrating that Equation (20) satisfies the local lepton number conservation in the trapped region.

Next, we take the diffusion limit of Equation (4). In this limit, ${H}^{i}/J$ should approach 0 (i.e., the radiation flux (Hⁱ) in the comoving frame vanishes). From this, the following relation can be derived:

$$\begin{eqnarray}&&{H}^{i}\sim -{{Ev}}^{i}+{F}^{i}-{P}^{{ij}}{v}_{j}\sim 0\quad \Longrightarrow \quad {F}^{i}\sim \displaystyle \frac{4}{3}{{Ev}}^{i},\end{eqnarray} \tag{ 26 }$$

where we take a simple closure relation ${P}^{{ij}}=\frac{E}{3}{\delta }^{{ij}}$ in the diffusion limit. Inserting this into the left-hand side of Equation (4), one can get

$$\begin{eqnarray}&&{\partial }_{t}E+{\partial }_{i}{F}^{i}+{\partial }_{\varepsilon }(-\varepsilon {P}^{{ij}}{\partial }_{i}{v}_{j})\\ &&\quad \sim {\partial }_{t}E+{\partial }_{i}\left(\displaystyle \frac{4}{3}{{Ev}}^{i}\right)-\displaystyle \frac{E}{3}{\delta }^{{ij}}{\partial }_{i}{v}_{j}-\varepsilon ({\partial }_{\varepsilon }{P}^{{ij}}){\partial }_{i}{v}_{j}\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\sim \;{\partial }_{t}E+{\partial }_{i}\left(\displaystyle \frac{4}{3}{{Ev}}^{i}\right)-{\partial }_{i}\left(\displaystyle \frac{1}{3}{{Ev}}^{i}\right)-\varepsilon ({\partial }_{\varepsilon }{P}^{{ij}}){\partial }_{i}{v}_{j}\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&\sim \;{\partial }_{t}E+{\partial }_{i}({{Ev}}^{i})-\varepsilon ({\partial }_{\varepsilon }{P}^{{ij}}){\partial }_{i}{v}_{j}.\end{eqnarray} \tag{ 29 }$$

Moving from the rhs of Equation (27) to that of Equation (28), we assumed that E has almost no spatial gradient well below the neutrino spheres (at high opacities) in the prebounce core. This is satisfied quite well as shown in the literature (Bruenn 1985), in which a nearly flat ${Y}_{{\nu }_{e}}$ profile is shown in their standard models within the central core with mass ∼0.5–0.6${M}_{\odot }$ after the central density exceeds ∼5 × 10¹³ g cm⁻³. Note that in Equation (28) the third term $-{\partial }_{i}({{Ev}}^{i}/3)$, which is originally a part of the advection terms in the energy space ${\partial }_{\varepsilon }(-\varepsilon {P}^{{ij}}{\partial }_{i}{v}_{j})$, balances with part of the spatial advection term ${\partial }_{i}(4{{Ev}}^{i}/3)$ and they lead to the second term in Equation (29). From this, one can clearly see how the apparent advection speed of E at each energy bin ε approaches the matter velocity vⁱ in the diffusion limit. Then, assuming that the neutrino number density at each energy bin can be approximately expressed as

$$\begin{eqnarray}&&{q}^{0}\sim { \mathcal J }\sim { \mathcal E }-2{{ \mathcal F }}^{i}{v}_{i}\sim { \mathcal E },\end{eqnarray} \tag{ 30 }$$

in the slow motion limit, when we divide Equation (29) by ε and integrate it in energy space, it can be summarized as

$$\begin{eqnarray}&&\int d\varepsilon [{\partial }_{t}{q}^{0}+{\partial }_{i}({q}^{0}{v}^{i})-({\partial }_{\varepsilon }{P}^{{ij}}){\partial }_{i}{v}_{j}]={\partial }_{t}{n}_{\nu }+{\partial }_{i}{n}_{\nu }{v}^{i}.\end{eqnarray} \tag{ 31 }$$

From this, one can clearly see that the neutrino number density ${n}_{\nu }(=\rho {Y}_{\nu })$ is transported with the same matter velocity vⁱ in the diffusion limit.

We employed a standard high-resolution, shock-capturing scheme and utilize the HLL (Harten–Lax–van Leer) scheme (Harten et al. 1983) to evaluate the numerical fluxes in space (Kuroda & Umeda 2010). As for the fastest and slowest characteristic wave speeds of the radiation field system (Equations (4) and (5)), we again use the same definition as in Kuroda et al. (2012; see also Shibata et al. 2011) and connect ${\lambda }_{\mathrm{rad},\mathrm{thin}}$ and ${\lambda }_{\mathrm{rad},\mathrm{thick}}$ smoothly via the variable Eddington factor χ as

$$\begin{eqnarray}&&{\lambda }_{\mathrm{rad}}=\displaystyle \frac{3\chi -1}{2}{\lambda }_{\mathrm{rad},\mathrm{thin}}+\displaystyle \frac{3(1-\chi )}{2}{\lambda }_{\mathrm{rad},\mathrm{thick}},\end{eqnarray} \tag{ 36 }$$

where ${\lambda }_{\mathrm{rad},\mathrm{thin}}$ and ${\lambda }_{\mathrm{rad},\mathrm{thick}}$ are the wave speed in the optically thin and thick limits, respectively.

To enforce the numerical flux of the radiation field in the opaque region as it asymptotically approaches the diffusion limit, we evaluate the energy flux (${F}_{\mathrm{hll}}^{0}$) and the momentum flux (${F}_{\mathrm{hll}}^{i}$) as

$$\begin{eqnarray}&&{F}_{\mathrm{hll}}^{0}=\displaystyle \frac{{\tilde{\lambda }}_{+}{F}_{L}^{0}-{\tilde{\lambda }}_{-}{F}_{R}^{0}+\epsilon {\tilde{\lambda }}_{-}{\tilde{\lambda }}_{+}({Q}_{R}^{0}-{Q}_{L}^{0})}{{\tilde{\lambda }}_{+}-{\tilde{\lambda }}_{-}},\end{eqnarray} \tag{ 37 }$$

and

$$\begin{eqnarray}{F}_{\mathrm{hll}}^{i} & = & \displaystyle \frac{{\epsilon }^{2}({\tilde{\lambda }}_{+}{F}_{L}^{i}-{\tilde{\lambda }}_{-}{F}_{R}^{i})+\epsilon {\tilde{\lambda }}_{-}{\tilde{\lambda }}_{+}({Q}_{R}^{i}-{Q}_{L}^{i})}{{\tilde{\lambda }}_{+}-{\tilde{\lambda }}_{-}}\\ & & +(1-{\epsilon }^{2})\displaystyle \frac{{F}_{L}^{i}+{F}_{R}^{i}}{2},\end{eqnarray} \tag{ 38 }$$

respectively (Audit et al. 2002; O'Connor & Ott 2013). Here, ${Q}_{L/R}^{\alpha }$ and ${F}_{L/R}^{\alpha }$ are the conservative variables and their corresponding fluxes, respectively, with L/R denoting the left/right states for the Riemann problem. All the radiation-hydrodynamical variables are defined at the cell center. For those cell centered (primitive) variables, we use a monotonized central reconstruction, which has second order accuracy in space, and obtain left/right states at the cell surface (see Kuroda & Umeda 2010 for a more detailed explanation). After the reconstruction, all the characteristic wave speeds of the matter and radiation fields are evaluated.

is a modification parameter to fit the numerical flux to the diffusion limit, which we take as

$$\begin{eqnarray}&&\epsilon ={\rm{min}}\left(1,\displaystyle \frac{1}{\kappa {\rm{\Delta }}x}\right),\end{eqnarray} \tag{ 39 }$$

where κ is the total opacity and ${\rm{\Delta }}x$ is the grid width (Audit et al. 2002; O'Connor & Ott 2013).

Regarding the advection term in energy space in all conservation equations (Equations (4), (5), and (16)), we define the advection fluxes at the interface of the energy bin as the same as in Müller et al. (2010) so that all energy-integrated advection terms will vanish. This can be achieved since both terms, ${\tilde{M}}^{\alpha }$, appearing in energy and momentum conservation equations, and $\varepsilon {q}^{\alpha \beta }$ in number conservation equation, are expressed in terms of linear combinations of radiation momenta, J, ${H}^{\alpha }$, ${L}^{\alpha \beta }...$, and we therefore can define their cell surfaced values with an appropriate weighting function to suppress the violation of all conservations simultaneously.

For all orders of the radiation momenta $X\in \{J,{H}^{\alpha },{L}^{\alpha \beta }...\}$, the following conditions need to be satisfied for number conservation,

$$\begin{eqnarray}&&\int d\varepsilon {\partial }_{\varepsilon }X=0\ \ \ \Longrightarrow \ \ \displaystyle \sum _{i}{\rm{\Delta }}{\varepsilon }_{i}\displaystyle \frac{{X}_{i+1/2}-{X}_{i-1/2}}{{\rm{\Delta }}{\varepsilon }_{i}}=0,\end{eqnarray} \tag{ 40 }$$

and for energy and momentum conservations,

$$\begin{eqnarray}\displaystyle \int d\varepsilon {\partial }_{\varepsilon }(\varepsilon X) & = & \displaystyle \int d\varepsilon (\varepsilon {\partial }_{\varepsilon }X+X)=0\\ & \Longrightarrow & \displaystyle \sum _{i}{\rm{\Delta }}{\varepsilon }_{i}\left[{\varepsilon }_{i}\displaystyle \frac{{X}_{i+1/2}-{X}_{i-1/2}}{{\rm{\Delta }}{\varepsilon }_{i}}+{X}_{i}\right]=0,\end{eqnarray} \tag{ 41 }$$

respectively. The rhs of the above equations represent the finite difference expressions with i and $i+1/2$ denoting the ith energy bin and the interface between the i- and $(i+1)\mathrm{th}$ energy bins, respectively. ${\rm{\Delta }}{\varepsilon }_{i}$ is the energy grid width ${\rm{\Delta }}{\varepsilon }_{i}={\varepsilon }_{i+1/2}-{\varepsilon }_{i-1/2}$. It is straightforwad to show that Equation (40) can be automatically satisfied for any cell surfaced quantities as long as they vanish at the outer boundary in the energy space. By introducing a definition ${X}_{i+1/2}\equiv {X}_{i}^{{\rm{L}}}+{X}_{i+1}^{{\rm{R}}}$, the rhs of Equation (41) can be expressed as

$$\begin{eqnarray}&&\displaystyle \sum _{i}[-({\varepsilon }_{i+1}-{\varepsilon }_{i}){X}_{i}^{{\rm{L}}}-({\varepsilon }_{i}-{\varepsilon }_{i-1}){X}_{i}^{{\rm{R}}}+{\rm{\Delta }}{\varepsilon }_{i}{X}_{i}]=0.\end{eqnarray} \tag{ 42 }$$

As in Müller et al. (2010), we can get the solution of Equation (42) as

$$\begin{eqnarray}&&{X}_{i}^{{\rm{L}}}\equiv \displaystyle \frac{{\rm{\Delta }}{\varepsilon }_{i}}{{\varepsilon }_{i+1}-{\varepsilon }_{i}}{X}_{i}{\xi }_{i},\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&{X}_{i}^{{\rm{R}}}\equiv \displaystyle \frac{{\rm{\Delta }}{\varepsilon }_{i}}{{\varepsilon }_{i}-{\varepsilon }_{i-1}}{X}_{i}(1-{\xi }_{i}),\end{eqnarray} \tag{ 44 }$$

where ${\xi }_{i}$ is a weighting function and is expressed as

$$\begin{eqnarray}&&{\xi }_{i}\equiv \displaystyle \frac{{f}_{i+1/2}^{\sigma }}{{f}_{i-1/2}^{\sigma }+{f}_{i+1/2}^{\sigma }}.\end{eqnarray} \tag{ 45 }$$

In this study, we used a "Harmonic" interpolation ($\sigma =1$) for the energy density $({f}_{i+1/2}^{\sigma })$ as

$$\begin{eqnarray}&&{f}_{i+1/2}^{\sigma }\equiv {\left[{\left(\displaystyle \frac{{E}_{({\varepsilon }_{i})}}{{\varepsilon }_{i}^{3}}\right)}^{1-{r}_{i+1/2}}{\left(\displaystyle \frac{{E}_{({\varepsilon }_{i+1})}}{{\varepsilon }_{i+1}^{3}}\right)}^{{r}_{i+1/2}}\right]}^{\sigma },\end{eqnarray} \tag{ 46 }$$

with ${r}_{i+1/2}=({\varepsilon }_{i+1/2}-{\varepsilon }_{i})/({\varepsilon }_{i+1}-{\varepsilon }_{i})$ (see Müller et al. 2010 for more details). Like these, just only by evaluating an appropriate radiation momenta ${X}_{i+1/2}$ at the energy bin surface, we can simultaneously vanish Equations (40)–(41) numerically and maintain energy, momentum, and number conservations.

After the explicit update, we solve the following simultaneous equation (e.g., Equation (35));

$$\begin{eqnarray}{\boldsymbol{f}}({{\boldsymbol{P}}}^{n+1}) & \equiv & \displaystyle \frac{{\boldsymbol{Q}}({{\boldsymbol{P}}}^{n+1})-{{\boldsymbol{Q}}}^{*}}{{\rm{\Delta }}t}+{{\boldsymbol{S}}}_{{\rm{adv,e}}}({{\boldsymbol{P}}}^{n+1})\\ & & +{{\boldsymbol{S}}}_{\nu {\rm{m}}}({{\boldsymbol{P}}}^{n+1})=0,\end{eqnarray} \tag{ 47 }$$

where ${\boldsymbol{Q}}$, ${{\boldsymbol{S}}}_{{\rm{adv,e}}}$, and ${{\boldsymbol{S}}}_{\nu {\rm{m}}}$ are expressed in terms of the primitive variables ${\boldsymbol{P}}$

$$\begin{eqnarray}&&{\boldsymbol{P}}=\left[\begin{array}{c}\rho \\ {u}_{i}\\ s\\ {Y}_{e}\\ {E}_{(\nu ,\varepsilon )}\\ {{F}_{(\nu ,\varepsilon )}}_{i}\end{array}\right],\end{eqnarray} \tag{ 48 }$$

at the $(n+1)\mathrm{th}$ time step. Obviously, the baryon number density does not change in this step, i.e., ${\rho }^{*}={\rho }^{n+1}$. To get the solutions of the above simultaneous equation, we employ the Newton–Raphson iteration scheme with the inversion of the following matrix until a sufficient convergence is achieved:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{f}}({{\boldsymbol{P}}}^{I})}{\partial {\boldsymbol{P}}}\delta {{\boldsymbol{P}}}^{I}=-{\boldsymbol{f}}({{\boldsymbol{P}}}^{I}),\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&{{\boldsymbol{P}}}^{I+1}={{\boldsymbol{P}}}^{I}+\delta {{\boldsymbol{P}}}^{I},\end{eqnarray} \tag{ 50 }$$

for $I=0,1,2.....,$ with the initial condition ${{\boldsymbol{P}}}^{0}={{\boldsymbol{P}}}^{*}$. As for the convergence criteria for the Newton–Raphson iteration, we monitor

$$\begin{eqnarray}&&\displaystyle \frac{| \delta {{\boldsymbol{P}}}^{I}| }{| {{\boldsymbol{P}}}^{I}| }\lt \mathrm{tol},\end{eqnarray} \tag{ 51 }$$

where tol represents a tolerance and we typically set $\mathrm{tol}={10}^{-8}$.

We note that even if we evaluate the advection terms in energy space explicitly in time as ${{{\boldsymbol{S}}}_{{\rm{adv,e}}}}^{n}$, Equations (47)–(51) do not change except the term ${{\boldsymbol{S}}}_{{\rm{adv,e}}}({{\boldsymbol{P}}}^{n+1})$ in Equation (47) is now moved to the explicit part.

Our method for time update from the nth to $(n+1)\mathrm{th}$ time step is summarized  as follows and schematically drawn in Figure 1.

• 1.  
  Geometrical Update. We first evolve all the BSSN and the gauge variables ${\boldsymbol{G}}$ = {${\tilde{\gamma }}_{{ij}}$, ${\tilde{A}}_{{ij}}$, ϕ, K, ${\tilde{{\rm{\Gamma }}}}^{i}$, α, ${\beta }^{i}$} from the nth to the $(n+1)\mathrm{th}$ time step.
• 2.  
  Explicit Update. In the second step, all the advection in space $({{\boldsymbol{S}}}_{\mathrm{adv},{\rm{s}}})$ and gravitational source $({{\boldsymbol{S}}}_{\mathrm{grv}})$ terms are added to obtain the fractional time step values ${{\boldsymbol{Q}}}^{*}$ and their consistent primitive variables ${\boldsymbol{P}}({{\boldsymbol{G}}}^{n+1},{{\boldsymbol{Q}}}^{*})=(\rho ,{u}_{i},s,{Y}_{e},{E}_{(\nu ,\varepsilon )},{{F}_{(\nu ,\varepsilon )}}_{i})$. Advection of the total lepton number (Equation (20)) is also performed simultaneously to evaluate ${Y}_{l}^{*}$.
• 3.  
  Implicit Update. Finally, quantities in the fractional time step (${{\boldsymbol{Q}}}^{*}$) are implicitly updated to those in the $(n+1)\mathrm{th}$ time step (${{\boldsymbol{Q}}}^{n+1}$) by the Newton–Raphson iteration until a sufficient convergence is obtained with a constraint that the local lepton fraction does not change (i.e., ${Y}_{l}^{n+1}={Y}_{l}^{*}$ since ${\rho }^{n+1}={\rho }^{*}$) in the diffusion region.

Zoom In Zoom Out Reset image size

Here we shortly summarize some technical details to achieve a sufficient lepton number conservation in practical core-collapse simulations. As we have already mentioned in Section 2.2, neutrino number conservation is formally satisfied by solving the energy–momentum conservation Equations (4)–(5). Especially, in the neutrino trapping regime, solving the energy conservation (Equation (4)) is practically identical to solving the neutrino number conservation (Equation (16)) which was proven in Section 2.2.1. However, in the finite difference method, it does not guarantee a perfect match between Equations (4)–(5) and (16) due to a discretization error. To minimize the difference, we add the following constraint to Equation (51)

$$\begin{eqnarray}&&\displaystyle \frac{| {Y}_{l}^{I}-{Y}_{l}^{*}| }{{Y}_{l}^{I}}\lt \mathrm{tol},\end{eqnarray} \tag{ 52 }$$

as another criterion to exit the Newton–Raphson iteration. Because of this additional criterion, the lepton number conservation is also satisfied when the Newton–Raphson iteration converges. Note that the local baryon numbers do not change in Equation (47), i.e., ${Y}_{l}^{*}={Y}_{l}^{n+1}$ should be met if the lepton number is conserved.

We also adopt another prescription in which ${Y}_{l}^{n+1}$ is used to achieve a better convergence for the Newton–Raphson method. In some cases, Y_e^I happens to go beyond the range of the employed EOS table during the iteration, especially when the Jacobian matrix is not well evaluated. In those cases, we use ${Y}_{l}^{n+1}$ to reset Y_e^I as below:

$$\begin{eqnarray}&&{Y}_{e}^{I}={Y}_{l}^{n+1}-{Y}_{{\nu }_{e}}^{I}+{Y}_{{\bar{\nu }}_{e}}^{I}.\end{eqnarray} \tag{ 53 }$$

With this reset, we found that Y_e^I hardly goes beyond the range of the EOS table and also does not take an unreasonable value for the supernova core. In addition, the number of the Newton–Raphson iteration can sometimes be reduced. This is because the reset value Y_e^I automatically satisfies the lepton number conservation ${Y}_{l}^{I}={Y}_{e}^{I}+{Y}_{{\nu }_{e}}^{I}-{Y}_{{\bar{\nu }}_{e}}^{I}={Y}_{l}^{n+1}$, which leads to faster convergence.

Without these two prescriptions, we observed that the central lepton fraction at bounce deviates maximally ∼0.05 from the value immediately after neutrino trapping sets in (${\rho }_{c}\gtrsim {10}^{12}\;{\rm{g}}\;{{\rm{cm}}}^{-3}$). Note that these ad hoc constraints are introduced just to find a more efficient path toward convergence in the implicit time update. The criterion for convergence is solely given by Equation (51). Since we do not include Equation (20) into (47), there is no inconsistency between the number of simultaneous equations ${\boldsymbol{f}}(P)$ and the unknown variables ${{\boldsymbol{P}}}^{n+1}$.

We first perform a diffusion wave test, by which we check the validity of our flux implementation in the diffusion limit. Following Pons et al. (2000), a Dirac δ-function-type radiation source is initially located at r = 0. Then we follow the diffusion of the source into the optically thick medium with zero absorptivity (${\kappa }_{a}=0$) and high scattering opacity (${\kappa }_{s}={10}^{2}$, 10⁵). The source term in Equations (4)–(5) thus becomes

$$\begin{eqnarray}&&{S}^{\mu }=-{\kappa }_{s}{H}^{\mu }.\end{eqnarray} \tag{ 54 }$$

The 3D computational domain is covered by 100 equidistant Cartesian zones in each direction (${\boldsymbol{x}}\in [-0.5,0.5]$) for a model with ${\kappa }_{s}={10}^{2}$. This corresponds to a Peclet number ${Pe}={\kappa }_{s}{\rm{\Delta }}x=1$. While in the other model with ${\kappa }_{s}={10}^{5}$ we raise the nested level to 2 to achieve sufficient resolution at the center and to check that the nested structure does not interfere with the radiation propagation in the diffusion limit. In this model, 32³ numerical cells cover the base domain ${\boldsymbol{x}}\in [-0.5,0.5]$ with a nested level of 2. The minimum grid width and Peclet number at the center are thus ${\rm{\Delta }}x=1/128\sim 0.78$ and ${Pe}={10}^{5}/128\sim 780$, respectively. Since Pe is much higher than unity, in Equation (39) approaches zero, which enables us to check whether the radiation transport in the diffusion limit is solved appropriately.

The analytical solution of the zeroth and first order radiation momenta profiles at time T (Pons et al. 2000) is given as

$$\begin{eqnarray}&&E(r,T)={\left(\displaystyle \frac{{\kappa }_{s}}{T}\right)}^{3/2}{\rm{exp}}\left(\displaystyle \frac{-3{\kappa }_{s}{r}^{2}}{4T}\right),\end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray}&&{F}_{r}(r,T)=\displaystyle \frac{r}{2T}E,\end{eqnarray} \tag{ 56 }$$

respectively. From Figure 2, it can be seen that our code can reproduce the analytical results quite well.

Zoom In Zoom Out Reset image size

Next, we move on to a shadow casting test to check the ability of our code and whether the anisotropic radiation field can be appropriately evolved in the free streaming limit. The initial setup is essentially the same as in Kanno et al. (2013). We set a perfect absorbing region at $x\in [2,3]$ and $\{y,z\}\in [-2,2]$ inside the numerical domain $x\in [0,12]$ and $\{y,z\}\in [-4,4]$. We impose a constant radiation $(E,{F}_{x})=(1,{f}_{\mathrm{max}})$ from the left boundary. Numerical resolution is ${\rm{\Delta }}x(={\rm{\Delta }}y={\rm{\Delta }}z)=12/192$ and we set ${f}_{\mathrm{max}}=\;0.999$. Within the perfect absorbing region, the absorbing opacity is set as ${\kappa }_{a}{\rm{\Delta }}x\sim {10}^{10}$, otherwise ${\kappa }_{a}=0$. The source term in Equations (4)–(5) thus becomes

$$\begin{eqnarray}&&{S}^{\mu }=-{\kappa }_{a}((J-{J}^{\mathrm{eq}}){u}^{\mu }+{H}^{\mu }),\end{eqnarray} \tag{ 57 }$$

with ${J}^{\mathrm{eq}}=0$ and ${u}^{\mu }=(1,0,0,0)$. In Figure 3, we show three different time snapshots (T = 3.5, 7.5, and 15.5) of the radiation energy density in a logarithmic scale. The absorbing region is expressed by white dashed line and the radiation shock front (x = T) is denoted by a vertical green line. From this figure, we see that the absorbing condition works appropriately in our scheme and the radiation front propagates with the light speed v = 1. Furthermore, the region beyond the absorbing box is barely contaminated by the radiation. These features are in good agreement with Kanno et al. (2013).

Zoom In Zoom Out Reset image size

The next test is a spherical expansion of the radiation field from a point-like source into an optically thin medium. The aim of this test is to check whether our code using the Cartesian coordinates can maintain the sphericity of the field during the expansion. By following Just et al. (2015), we define the radiation source ${ \mathcal S }$ with radius ${r}_{{ \mathcal S }}=1.5$ which centers at the origin of the 2D Cartesian domain with ${\boldsymbol{x}}\in [-7.5,7.5]$. We also define the purely absorbing region ${ \mathcal A }$ centered at ${\boldsymbol{x}}=(3.5,0)$ with radius ${r}_{{ \mathcal A }}=1$. Inside those two regions, the absorptivity ${\kappa }_{a}({\boldsymbol{x}})$ and equilibrium energy density ${J}^{\mathrm{eq}}({\boldsymbol{x}})$ are set as

$$\begin{eqnarray}{\kappa }_{a}({\boldsymbol{x}})=\left\{\begin{array}{ll}10{\rm{exp}}\{-{(4\sqrt{{x}^{2}+{y}^{2}}/{r}_{{ \mathcal S }})}^{2}\}, & {\boldsymbol{x}}\in { \mathcal S }\\ 10, & {\boldsymbol{x}}\in { \mathcal A },\end{array}\right.\end{eqnarray} \tag{ 58 }$$

and

$$\begin{eqnarray}{J}^{\mathrm{eq}}({\boldsymbol{x}})=\left\{\begin{array}{ll}{10}^{-1}, & {\boldsymbol{x}}\in { \mathcal S }\\ 0, & {\boldsymbol{x}}\in { \mathcal A },\end{array}\right.\end{eqnarray} \tag{ 59 }$$

respectively. In other regions, we set ${\kappa }_{a}({\boldsymbol{x}})=0$ and ${J}^{\mathrm{eq}}({\boldsymbol{x}})=0$. We neglect the scattering opacity ${\kappa }_{{\rm{s}}}({\boldsymbol{x}})$. The source term in Equations (4)–(5) is thus the same as Equation (57). Regarding the initial setup for the zeroth and first order momenta of radiation field, we assume an arbitrary dilute radiation field as

$$\begin{eqnarray}&&E={10}^{-9}\end{eqnarray} \tag{ 60 }$$

$$\begin{eqnarray}&&| F| /E={10}^{-10}.\end{eqnarray} \tag{ 61 }$$

The 2D numerical domain is covered by 196 equally distant Cartesian zones in each direction with two nested levels, i.e., ${\rm{\Delta }}x$ varies from 15/784 to 15/196.

In this test, we use the following formula for the variable Eddington factor χ (Levermore 1984):

$$\begin{eqnarray}&&\chi =\displaystyle \frac{3+4{\bar{F}}^{2}}{5+2\sqrt{4-3{\bar{F}}^{2}}}.\end{eqnarray} \tag{ 62 }$$

Only in this propagation test, we take two possible options for evaluating the characteristic wave speeds of the radiation field in the optically thin limit ${\lambda }_{\pm }$. Assuming a flat spacetime and using the following closure relation in the optically thin limit

$$\begin{eqnarray}&&{P}^{\mu \nu }=E\displaystyle \frac{{F}^{\mu }{F}^{\nu }}{{F}_{k}{F}^{k}},\end{eqnarray} \tag{ 63 }$$

the first one is (Shibata et al. 2011)

$$\begin{eqnarray}&&{\lambda }_{\pm }^{{\rm{S}},i}={\rm{max/min}}\left[\displaystyle \frac{{{EF}}^{i}}{{F}_{k}{F}^{k}},\pm \displaystyle \frac{{F}^{i}}{\sqrt{{F}_{k}{F}^{k}}}\right],\end{eqnarray} \tag{ 64 }$$

and the second one is (Skinner & Ostriker 2013)

$$\begin{eqnarray}&&{\lambda }_{\pm }^{{\rm{SO}},i}=\left\{{\mu }^{i}\bar{F}\pm \sqrt{\displaystyle \frac{2}{3}(4-3{\bar{F}}^{2}-\sqrt{4-3{\bar{F}}^{2}})+2{{\mu }^{i}}^{2}(2-{\bar{F}}^{2}-\sqrt{4-3{\bar{F}}^{2}})}\right\}/\sqrt{4-3{\bar{F}}^{2}}.\end{eqnarray} \tag{ 65 }$$

Here, ${\mu }^{i}\equiv \mathrm{cos}\theta$ is determined by the angle θ which defines the orientation of the energy flux ${\boldsymbol{F}}$ relative to the interface normal xⁱ. Note that in the rest of the paper we use only Equation (64) for the free streaming part.

Also of note, we limit the flux factor $f\equiv \sqrt{{\gamma }^{{ij}}{F}_{i}{F}_{j}}/E$ to be less than the maximum allowed value ${f}_{\mathrm{max}}$ by modifying the first-order moment as

$$\begin{eqnarray}&&{F}_{i}\to {\rm{min}}({f}_{\mathrm{max}}/f,1){F}_{i},\end{eqnarray} \tag{ 66 }$$

in every time step. Note that setting ${f}_{\mathrm{max}}\lt 1$ means that we add contribution from the isotropic radiation pressure ${P}_{\mathrm{thick}}^{{ij}}$ according to Equation (6). The reason of this modification will be discussed later in this section.

In Figure 4, we show results for three different models taking different values for ${f}_{\mathrm{max}}$ and using different evaluation formulae for ${\lambda }_{\pm }^{{\rm{S}}/\mathrm{SO}}$. The upper two panels (panels (a) and (b)) and the bottom one (panel (c)) are models using Equations (64) and (65), respectively. In panels (a) and (c), we take ${f}_{\mathrm{max}}=1$, while ${f}_{\mathrm{max}}=0.93$ is taken in panel (b). The color represents the energy density E in a logarithmic scale at four different time slices ($T=1,3,5,7$). The inner two squares represent the boundary of the nested structure and the dotted circle shows the purely absorbing region ${ \mathcal A }$. From models (a) and (c), which use ${f}_{\mathrm{max}}=1$, we find non-isotropy which can be seen as corrugated patterns behind the radiation front, obviously with the exception of the absorbing circle ${ \mathcal A }$ and its shadowing region. Furthermore, in model (a) at T = 7, another remarkable violation of isotropy is seen at 45° away from the coordinate axises. Note that from snapshot at T = 7 in panel (a), this violation seems to be associated with the nested structure. We, however, confirmed that this is not the artifact of that by performing a run with a zero nested level. By comparing models (a) and (c), the violation seen at 45° is eliminated in model (c). Model (c) takes into account the propagation angle (${\mu }^{i}$) relative to the interface normal more explicitly than model (a) when we evaluate ${\lambda }_{\pm }^{\mathrm{SO}}$. We therefore consider a more accurate evaluation in the characteristic wave speeds ${\lambda }_{\pm }$ for the numerical flux (Equations (37)–(38)) is key to follow the radiation in an optically thin limit properly. While in panel (b), in which we use ${f}_{\mathrm{max}}=\;0.93$, we do not see any spherical symmetry breaking.

Zoom In Zoom Out Reset image size

We also apply this test to a supernova core profile. We fix the hydrodynamical background, which has a typical core profile at ${T}_{\mathrm{pb}}\sim 100$ ms, and follow only the propagation of neutrino radiation. The neutrino transport part is identical to our practical calculation reported later in Section 5 and all relevant neutrino–matter interactions, gravitational redshift, and Doppler shift terms are taken into account. There is thus a transition from an optically thick to a thin regime. We calculate two models with different values for ${f}_{\mathrm{max}}$. ${\lambda }_{\pm }^{{\rm{S}},{\rm{i}}}$ (Equation (64)) is used in both models to evaluate the characteristic wave speeds. In Figure 5, we show the radial profiles of (electron-type) neutrino luminosity for two cases: one is calculated with ${f}_{\mathrm{max}}=\;0.999$ (black diamonds) and the other is with ${f}_{\mathrm{max}}=\;0.93$ (red filled triangles). As can be clearly seen, beyond the shock position ($R\gtrsim 100$ km) where it is optically thin regime, a large scatter is seen at R ≳ 300–400 km for a model with ${f}_{\mathrm{max}}=\;0.999$. The spatial location of the highly deviated radiation is again concentrated $\sim 45^\circ$ away from the coordinate axes, which is the same as in the previous simple propagation test (panel (a) in Figure 4). Meanwhile, in model with ${f}_{\mathrm{max}}=\;0.93$, there is little deviation and the surface-integrated local luminosity stays nearly constant in the transparent region, as it should be.

Zoom In Zoom Out Reset image size

From these two tests, we are currently enforced to set a ceiling value for the flux factor $f\leqslant {f}_{\mathrm{max}}=\;0.93$ in the practical core-collapse simulation to follow spherical-like propagation stably. Since ${f}_{\mathrm{max}}=1$ is the physically correct value, one may suspect that our modification can be an obstacle for comparison of the emergent radiation profile with previous studies. From previous fully relativistic Boltzmann transport simulation (Liebendörfer et al. 2001), however, it was shown that the flux factor higher than, e.g., ∼0.93 appears only beyond R ≳ 500–1000 km. Therefore, even if we evaluate the emergent neutrino profile at the same radius R = 500 km as previous studies, it is altered only slightly by a few percent and the modification cannot be a significant obstacle in this study. Obviously, the SN hydrodynamics itself is not affected by our artificial treatment since it is applied only to the optically thin region.

To check the energy-coupling terms for gravitational redshift and Doppler shift, $\sqrt{\gamma }\alpha {\partial }_{\varepsilon }(\varepsilon {\tilde{M}}_{(\varepsilon )}^{\mu }{n}_{\mu })$ and $\sqrt{\gamma }\alpha {\partial }_{\varepsilon }(\varepsilon {\tilde{M}}_{(\varepsilon )}^{\mu }{\gamma }_{i\mu })$ in Equations (4)–(5), we repeat the same tests performed in Müller et al. (2010) and O'Connor (2014). We consider the propagation of radiation from a sphere with radius R = 10 km through a curved spacetime and sharp velocity profile in spherical symmetry. Regarding the sharp velocity profile, we take the following:

$$\begin{eqnarray}{u}_{r}=\left\{\begin{array}{cl}0, & r\leqslant 70\;{\rm{km}}\\ -0.2c\left(\displaystyle \frac{r-70\;{\rm{km}}}{10\;{\rm{km}}}\right), & 70\;{\rm{km}}\leqslant r\leqslant 80\;{\rm{km}}\\ -0.2c{\left(\displaystyle \frac{80{\rm{km}}}{r}\right)}^{2}, & r\geqslant 80\;{\rm{km}}.\end{array}\right.\end{eqnarray} \tag{ 67 }$$

Here, u_r is the radial component of u_i. As for the curved spacetime, we take the same energy density profile for matter which is used in the second test problem in Section 4.3 and solve the Hamiltonian constraint to obtain the conformal factor ψ and lapse α. When we solve the Hamiltonian constraint, we assume a zero velocity profile (${v}_{i}=0={\beta }_{i}$) and the conformally flat approximation (${\gamma }_{{ij}}={\psi }^{4}{\tilde{\gamma }}_{{ij}}={\psi }^{4}{\delta }_{{ij}}$). As for the initial neutrino profile within the radiating sphere, we take the zero chemical potential ${\mu }_{\nu }=0$ and a temperature of 5 MeV, i.e., $\langle \varepsilon \rangle \sim 15.7\;{\rm{MeV}}$, in the free streaming limit $E=\sqrt{{\gamma }_{{ij}}{F}^{i}{F}^{j}}$. Here $\langle \varepsilon \rangle$ is the mean energy of neutrinos as measured in the comoving frame. Outside of the central radiating sphere, we choose an arbitral dilute radiation field with $\langle \varepsilon \rangle \sim 15.7\;{\rm{MeV}}$. We do not evolve neutrinos inside the radiating sphere and follow the propagation only at $R\geqslant 10$ km. The 3D computational domain is a cubic box with 3000 km width (i.e., the outer boundary is at the radius of 1500 km from the origin) and nested boxes with 7 refinement levels are embedded. Each box contains 64³ cells and the minimum grid size near the origin is ${\rm{\Delta }}x=366$ m. As for the energy grid of the neutrino radiation field, we use logarithmically spaced 20 energy bins $({N}_{\varepsilon }=20)$ which center from $\varepsilon =1$ to 300 MeV. We calculate two models with and without the sharp velocity profile. Both models are calculated in the curved spacetime.

For the given velocity and spacetime profiles, the free streaming neutrinos propagate by obeying the following relations (Müller et al. 2010; O'Connor 2014):

$$\begin{eqnarray}&&W\alpha (1+{v}_{r})\langle \varepsilon \rangle ={\rm{const,}}\end{eqnarray} \tag{ 68 }$$

$$\begin{eqnarray}&&\sqrt{\gamma }\alpha {r}^{2}{\gamma }^{{rr}}\int d\varepsilon {F}_{r}\equiv {L}_{\mathrm{eul}}={\rm{const,}}\end{eqnarray} \tag{ 69 }$$

where W is the Lorentz factor, and ${v}_{r}={u}_{r}/{u}^{t}$. $\sqrt{\gamma }={\psi }^{6}$ and ${\gamma }^{{rr}}={\psi }^{-4}$ are the determinant and rr-component of the three metrics, respectively, in the conformally flat approximation. ${L}_{\mathrm{eul}}$ is the luminosity measured by an Eulerian observer. Equation (69) is derived from the stationary solution of Equation (4) with neglecting the source terms.

In the top panel of Figure 6, we show the background profiles for ψ, α, and $1-{u}_{r}$. $\langle \varepsilon \rangle$ and ${L}_{\mathrm{eul}}$ are plotted in the middle and bottom panels, respectively. In the middle panel, our results (filled triangles) and analytical ones (solid lines) are plotted. The analytical expression for $\langle \varepsilon \rangle$ is derived from Equation (68), where the constant in the rhs is evaluated at R = 10 km. Numerical results are taken after they settle down in nearly stationary.

Zoom In Zoom Out Reset image size

As seen in the middle panel, our results reproduce the analytical ones quite well except at the shock front. We consider that the less agreement at the shock front originates from the coarse spatial resolution. The finite difference expression for the term ${{\rm{\nabla }}}_{\alpha }{u}_{\beta }$ in the energy-coupling terms cannot capture the sharp velocity profile with enough accuracy. From the bottom panel, we can find the Eulerian luminosity stays nearly constant with a few percent deviation beyond $R\gtrsim 20$ km. Note that the zig-zag pattern seen in ${L}_{\mathrm{eul}}$ is an artifact of the mesh refinement.

We employ a 15M_⊙ progenitor (model "s15s7b2" in Woosley & Weaver 1995) and follow core-collapse, bounce, and initial postbounce phase up to ∼50 ms. We use the EOS of Lattimer & Douglas Swesty (1991; LS EOS) with an incompressibility parameter of K = 220 MeV. Even though our choice of K = 220 MeV is different from the one (K = 180 MeV) used in Liebendörfer et al. (2005), Sumiyoshi et al. (2005), Müller et al. (2010), Müller & Janka (2014), and Thompson et al. (2003) showed that differences in the neutrino profiles among models with the different K of LS EOS are a few percent around core bounce and less than ∼10% for the first ∼200 ms after bounce. We thus consider that the different choice of K barely disturbs the aim of our comparison study.

As for neutrino opacities, the standard weak interaction set in Bruenn (1985) and Rampp & Janka (2002) plus nucleon–nucleon bremsstrahlung (Hannestad & Raffelt 1998) is taken into account (see Table 1 and Appendix A for more details). For simplicity, we neglect higher harmonic angular dependence of the reaction angle when we calculate the source terms for neutrino electron scattering, thermal pair production, and the annihilation of neutrinos, and nucleon–nucleon bremsstrahlung (i.e., ${B}_{(\varepsilon ),\mathrm{nes}/\mathrm{tp}}^{0,\mu }$ and ${C}_{(\varepsilon ),\mathrm{nes}/\mathrm{tp}}^{1,\mu }$ are set to be 0; see Appendices A.3 and A.4).

Table 1.  The Opacity Set Included in this Study and their References

Process	Reference	Summarized In
$n{\nu }_{e}\leftrightarrow {e}^{-}p$	Bruenn (1985), Rampp & Janka (2002)	Appendix A.1
$p{\bar{\nu }}_{e}\leftrightarrow {e}^{+}n$	Bruenn (1985), Rampp & Janka (2002)	Appendix A.1
${\nu }_{e}A\leftrightarrow {e}^{-}A^{\prime}$	Bruenn (1985), Rampp & Janka (2002)	Appendix A.1
$\nu p\leftrightarrow \nu p$	Bruenn (1985), Rampp & Janka (2002)	Appendix A.2
$\nu n\leftrightarrow \nu n$	Bruenn (1985), Rampp & Janka (2002)	Appendix A.2
$\nu A\leftrightarrow \nu A$	Bruenn (1985), Rampp & Janka (2002)	Appendix A.2
$\nu {e}^{\pm }\leftrightarrow \nu {e}^{\pm }$	Bruenn (1985)	Appendix A.3
${e}^{-}{e}^{+}\leftrightarrow \nu \bar{\nu }$	Bruenn (1985)	Appendix A.4
${NN}\leftrightarrow \nu \bar{\nu }{NN}$	Hannestad & Raffelt (1998)	Appendix A.5

Note. Note that ν, in neutral current reactions, represents all species of neutrinos (${\nu }_{e},{\bar{\nu }}_{e},{\nu }_{x}$) with ${\nu }_{x}$ representing heavy lepton neutrinos (i.e., ${\nu }_{\mu },{\nu }_{\tau }$ and their anti-particles).

Download table as:  ASCIITypeset image

In this study, we investigate two models, 1DG15 and 3DG15, both of which are computed in the 3D Cartesian coordinates. While model 3DG15 is calculated without any artificial constraints, model 1DG15 is considered to mimic the 1D model by artificially suppressing the non-radial component of the flow velocity u_i as

$$\begin{eqnarray}&&{u}_{i}=\displaystyle \frac{{x}^{j}{u}_{j}}{| x{| }^{2}}{x}^{i}.\end{eqnarray} \tag{ 70 }$$

Although this artificial elimination could potentially lead to a shift of the kinetic energy into the thermal one, our previous study (Kuroda et al. 2012) showed that the violation of the momentum constraint is negligible during the early postbounce phase (${T}_{\mathrm{pb}}\lesssim 100$ ms where ${T}_{\mathrm{pb}}$ denotes the postbounce time).

The 3D computational domain is a cubic box with 8000 km width (i.e., the outer boundary is at the radius of 4000 km from the origin) and nested boxes with 5 refinement levels, at the beginning of calculation, to 8 refinement levels, when the central rest mass density reaches 5 × 10¹³ g cm⁻³, are embedded without any spatial symmetry. Each box contains 64³ cells and the minimum grid size near the origin at bounce is ${\rm{\Delta }}x=488$ m. In the vicinity of the stalled shock front $R\sim 100$ km, our resolution achieves ${\rm{\Delta }}x\sim 3.6$ or 7.2 km, i.e., the effective angular resolution becomes 3.6 km/100 km $\sim 2^\circ$ or $\sim 4^\circ$, which is considered to be a rather too coarse resolution to follow a nearly spherical structure in the Cartesian grids. Compared to recent 3D-GR studies (Ott et al. 2012; Kuroda et al. 2014), the resolution is still approximately two times coarser. The time step size is ${\rm{\Delta }}t={10}^{-7}\;{\rm{s}}$ which corresponds to ∼0.06 in the CFL number. The main reason for taking such a small time step is to get a rapid convergence during the iteration in our implicit time update part (see Section 3.3), which is not restricted by the CFL condition (of the explicit update in a pure hydrodynamics case). As for the energy grid of the neutrino radiation field, we use logarithmically spaced 20 energy bins $({N}_{\varepsilon }=20)$ which center from $\varepsilon =1$ to 300 MeV.

In this section we start to make a detailed comparison first before bounce (Section 5.2.1) and then after bounce (Section 5.2.2) between our code (models 1DG15 and 3DG15), AGILE-BOLTZTRAN (model ABG15), VERTEX-PROMETHEUS (model VXG15), and (partly) from VERTEX-COCONUT (model BMG15). In the following, we call the results from the latter three codes the reference results.

5.2.1. Before Bounce

In the upper panel of Figure 7, we plot the central (matter) entropy s, electron fraction Y_e, and the total lepton fraction ${Y}_{l}={Y}_{e}+{Y}_{\nu }$ as a function of the central (rest mass) density ${\rho }_{{\rm{c}}}$ for model 1DG15 (black), ABG15 (blue), and VXG15 (green), respectively. The lower panel of Figure 1 shows the evolution of ${\rho }_{{\rm{c}}}$ and the deviation of the ADM mass ${\rm{\Delta }}{M}_{\mathrm{ADM}}$ from its initial value. Here ${M}_{\mathrm{ADM}}$ is given by

$$\begin{eqnarray}{M}_{\mathrm{ADM}} & = & \displaystyle \int {{dx}}^{3}\left[\left({S}_{0}+\displaystyle \int d\varepsilon {E}_{(\nu ,\varepsilon )}\right){e}^{-\phi }\right.\\ & & +\left.\displaystyle \frac{{e}^{5\phi }}{16\pi }\left({\tilde{A}}_{{ij}}{\tilde{A}}^{{ij}}-\displaystyle \frac{2}{3}{K}^{2}-{\tilde{\gamma }}^{{ij}}{\tilde{R}}_{{ij}}{e}^{-4\phi }\right)\right]\\ & & +\displaystyle \int d\sigma \left[{e}^{-\phi }\left({S}^{i}+\displaystyle \int d\varepsilon {{F}_{(\nu ,\varepsilon )}}^{i}\right)\cdot \hat{{\boldsymbol{n}}}\right],\end{eqnarray} \tag{ 71 }$$

where the second line denotes energy loss due to momentum and neutrino energy flux through the numerical boundary and $\hat{{\boldsymbol{n}}}$ represents a unit normal vector to the surface element $d\sigma$. In the above surface integration, we neglect energy loss due to gravitational wave emission since it is negligibly small ($\sim {10}^{-11}{M}_{\odot }{c}^{2}$; e.g., Scheidegger et al. 2010) compared to the violation of the ADM mass in the CCSN environment (see, e.g., Kotake 2013 for a review).

Zoom In Zoom Out Reset image size

The top panel of Figure 1 shows that the neutrino trapping starts ${\rho }_{{\rm{c}}}\sim 2\times {10}^{12}$ g cm⁻³ for model 1DG15 and the lepton fraction remains constant with Y_l = 0.323 afterward. The evolution of our central Y_e and Y_l (black lines) is quantitatively in good agreement with VXG15 (green linel see also Buras et al. 2006b; Müller et al. 2010) and ABG15 (blue line). As already explained in Section 2.2, we solve the advection equation (Equation (20)) of the total lepton number (density) $\rho {Y}_{l}$ as a constraint to ensure the lepton number conservation. Because of the treatment, the evolution of the lepton/electron fraction is in excellent agreement with the reference results. As already pointed out by O'Connor (2014), this also relies on the accurate implementation of inelastic neutrino–electron scattering, energy–bin coupling, and the appropriate closure relation.

After the core deleptonization ceases (i.e., Y_l stays nearly constant with increasing central density), the inner core evolves almost adiabatically and the entropy remains nearly constant as it should be. A small breaking of the adiabaticity (decrease by 3.8% in the central entropy) is seen in model 1DG15 before bounce (see also O'Connor 2014). However, the (time-averaged) value $s\sim 1.3\;{k}_{{\rm{B}}}\;{{\rm{baryon}}}^{-1}$ is in roughly good agreement with the reference results (see also Liebendörfer et al. 2005; Sumiyoshi et al. 2005; Buras et al. 2006b; Müller et al. 2010) and this would not have a big impact on the subsequent core evolution due to the short simulation time in this study.

As for the total energy conservation (bottom panel of Figure 3), it is maximally violated with the amount of ${\rm{\Delta }}{M}_{\mathrm{ADM}}\sim 4\times {10}^{-4}\;{M}_{\odot }$ for model 1DG15 (i.e., $\sim 8\times {10}^{50}$ erg). The violation at bounce is slightly worse than that (${{\rm{\Delta }}}_{\mathrm{ADM}}\sim 5\times {10}^{50}$ erg) of the VERTEX-CoCoNut code (Müller et al. 2010), which should be improved and kept much smaller in more precise CCSN modeling. As one would anticipate, the violation is bigger for model 1DG15 (dashed black line) with the artificial elimination of the non-radial velocity than for model 3DG15 without (dashed red line). In our previous study (with the gray M1 scheme; Kuroda et al. 2012), the violation of the ADM mass is typically one order of magnitude smaller than that for the corresponding 3D model with (approximately) twice higher resolution. Because the computational time for the implicit update in the current code is very expensive (using our best available resources), we are now forced to employ a quite coarse resolution. It is important to clarify how the violation of the total energy conservation would be improved with increasing numerical resolution. We leave this for future work.

Figure 8 shows a spectral shape of the neutrino distribution function $f(\nu ,\varepsilon )$ (filled triangles),

$$\begin{eqnarray}&&f(\nu ,\varepsilon )=\displaystyle \frac{{J}_{(\nu ,\varepsilon )}}{4\pi {\varepsilon }^{3}},\end{eqnarray} \tag{ 72 }$$

which is estimated at the innermost grid point of model 1DG15 when the central density ${\rho }_{c}$ reaches ${10}^{\mathrm{11,12},\mathrm{13,14}}$ g cm⁻³, and at bounce, respectively. Solid curves represent the Fermi–Dirac distribution at equilibrium. From the left panel, it can be seen that β-equilibrium is achieved for electron-type neutrino (${\nu }_{e}$) at ${10}^{12}\;{\rm{g}}\;{{\rm{cm}}}^{-3}\lesssim {\rho }_{c}\lesssim {10}^{13}\;{\rm{g}}\;{{\rm{cm}}}^{-3}$, which is consistent with the neutrino trapping density as shown in Figure 7. In Bruenn (1985), the β-equilibrium for ${\nu }_{e}$ was obtained after ${\rho }_{c}=2.46\times {10}^{12}$ g cm⁻³ in their model that corresponds to ours ("standard" model). Note that at bounce, electron-type neutrinos, in the low-energy range $\lesssim 5\;{\rm{MeV}}$, slightly violate the Fermi blocking, i.e., the distribution function $f({\nu }_{e})$ exceeds one. The excess, however, is ≲10⁻⁵ and we consider that it is a negligible amount. As shown from the middle (${\bar{\nu }}_{e}$) and right panels (${\nu }_{X}$) of Figure 4, other neutrino species are thermalized only after ${\rho }_{{\rm{c}}}$ exceeds 10¹⁴ g cm⁻³. These features are quantitatively consistent with Bruenn & Mezzacappa (1997) and Rampp & Janka (2002). It may be surprising that regardless of the use of different EOS and different hydrodynamics codes, the trapping density of ${\nu }_{e}$ (${\rho }_{\mathrm{trap}}=2\sim 3\times {10}^{12}$ g cm⁻³) in modern simulations (e.g., top panel of Figure 3) is in good agreement with the pioneering work in the 1980s (Bruenn 1985). It should be also noted that for the more accurate determination of the core deleptonization the improved electron capture rates (Langanke & Martínez-Pinedo 2000; Juodagalvis et al. 2010) need to be implemented as in Langanke et al. (2003) and Lentz et al. (2012a, 2012b).

Zoom In Zoom Out Reset image size

So far, we have shown that our M1 scheme can capture several important phenomena regarding deleptonization, such as the neutrino trapping and the conservation of the lepton fraction in the diffusion region. As already denoted in Section 2.1, this is not trivial in the finite difference method especially when one transports conservative radiation moments $({E}_{(\varepsilon )},{{F}_{(\varepsilon )}}_{i})$ instead of the corresponding comoving variables of $({{ \mathcal J }}_{(\varepsilon )},{{{ \mathcal H }}_{(\varepsilon )}}_{i})$. A key is finding an appropriate Eddington factor ${\chi }_{(\varepsilon )}$ by which the neutrino energy flux approaches ${{F}_{(\varepsilon )}}^{i}\to 4/3{E}_{(\varepsilon )}{v}^{i}$ in the diffusion limit (see Equation (26)). Since this relation can be achieved only when ${P}_{(\varepsilon )}^{{ij}}={E}_{(\varepsilon )}{\gamma }^{{ij}}/3$ holds, ${\chi }_{(\varepsilon )}$ and $\bar{F}$ should approach 1/3 and 0 (e.g., our closure relation (Equation (6))) in the limit, respectively.

To show that both $\bar{F}=\sqrt{{h}_{\mu \nu }{H}_{(\varepsilon )}^{\mu }{H}_{(\varepsilon )}^{\nu }/{J}_{(\varepsilon )}^{2}}$ and ${{F}_{(\varepsilon )}}^{i}$ actually approach 0 and $4/3{E}_{(\varepsilon )}{v}^{i}$ in the opaque region, we plot in Figure 9 the radial profiles of $\langle {F}_{r}\rangle /\langle E\rangle$ (solid lines) and $\langle {H}_{r}\rangle /\langle J\rangle$ (dashed–dotted lines) at different ${\rho }_{c}$. Here $\langle X\rangle \equiv \int d\varepsilon X$ represents the energy integration of X. At ${\rho }_{{\rm{c}}}={10}^{12}$ g cm⁻³, both solid and dashed–dotted green lines almost coincide. However, as the infalling matter velocity comes closer to being relativistic (${\rho }_{{\rm{c}}}\gtrsim {10}^{13}$ g cm⁻³), both lines start to split especially within $R\lesssim 70$ km. In the central region ($R\lesssim 70$ km), $\langle {H}_{r}\rangle /\langle J\rangle$ approaches 0 toward the center, whereas $\langle {F}_{r}\rangle /\langle E\rangle$ becomes negatively large with the peak being around $R\sim 30$ km and then converges to zero to the center. We also plot the radial velocity of matter (V_r) measured in the Eulerian frame and multiplied by 4/3 (blue diamonds) at ${\rho }_{{\rm{c}}}={10}^{14}$ g cm⁻³.

Zoom In Zoom Out Reset image size

Because of our appropriate evaluation for the Eddington factor, the flux factor measured in the Eulerian frame (approximated here by $\langle {F}_{r}\rangle /\langle E\rangle$) nicely matches with $4{V}_{r}/3$ in the optically thick region. This neutrino advection is essentially important for the radiation energy (${E}_{(\varepsilon )}$) to move with the same velocity vⁱ with matter in the opaque region (see Equations (26)–(29)). Here we shortly comment on the definition of the flux factor $\bar{F}$. In our previous study (Kuroda et al. 2012), we employed the definition $\bar{F}\equiv \sqrt{{\gamma }^{{ij}}{F}_{i}{F}_{j}}/E$, which is one of the candidates for $\bar{F}$ (Shibata et al. 2011). As can be clearly seen from the split between the solid and dashed–dotted lines in Figure 9, we show that our previous choice of the flux factor is not adequate because the optically thick medium moves (albeit mildly) relativistically in the collapsing core.

5.2.2. After Bounce

In Figure 10, we show the radial profiles of various quantities (the rest mass density ρ, radial velocity ${v}_{{\rm{r}}}$, matter temperature T, entropy s, averaged neutrino energy E_ν, and electron/neutrino fraction ${Y}_{e}/{Y}_{l}$) at the selected time slices shortly after bounce. Here, we define the average neutrino energy E_ν as

$$\begin{eqnarray}&&{E}_{\nu }\equiv \displaystyle \frac{\int d\varepsilon \;{\varepsilon }^{3}\;f(\nu ,\varepsilon )}{\int d\varepsilon \;{\varepsilon }^{2}f(\nu ,\varepsilon )}.\end{eqnarray} \tag{ 73 }$$

Zoom In Zoom Out Reset image size

When the central density exceeds nuclear saturation densities (top left panel of Figure 10), the core bounces because of the strong repulsive force of nuclear matter. Then the bounce shock is formed (green line, left middle panel of Figure 10) which propagates outward with dissociating infalling heavy nuclei into free protons and neutrons. The production of enormous free protons/neutrons significantly enhances the electron capture process, ${e}^{-}p\to {\nu }_{e}n$, behind the stalling shock. Immediately after bounce, those high-energy neutrinos (top right panel) are still trapped inside the optically thick medium behind the shock. The medium, however, quickly becomes transparent to neutrinos and those neutrinos are liberated suddenly as a burst. This neutronization burst enhances further electron capture behind the shock due to the continuous deviation from the β-equilibrium, leading to a characteristic trough in the Y_e profile seen at $R\sim 50$ km behind the shock (right bottom panel of Figure 10). Due to the energy loss by the photodissociation of the iron nuclei and the rapid neutrino leakage, the bounce shock stalls at $R\sim 70$ km within ${T}_{\mathrm{pb}}\sim 3$ ms (left panels of Figure 10). Such dynamical features are commonly seen in the reference models (Liebendörfer et al. 2005; Müller et al. 2010). This indicates that our M1 scheme can capture the basic behaviors of the neutrino propagation from the optically thick to thin medium (otherwise it would result in either the absence or the different position of the Y_e trough).

How the neutronization burst is produced is more clearly depicted in Figure 11 where we plot the radial profiles of the neutrino (${\nu }_{e}$) energy flux (solid lines) and the radial velocity of matter (dashed–dotted lines). At ${T}_{\mathrm{pb}}=0.7$ ms, enormous neutrinos are still trapped and confined behind the shock, which is shown as a sharp peak in the energy flux around $10\lesssim R\lesssim 20$ km. At ${T}_{\mathrm{pb}}=1.7$ ms, these neutrinos overtake the shock front because of the lowering opacity outward. Then the pulse of the neutrino burst propagates freely to the optically thin region (time label larger than 4) and eventually emerges out of the computational domain (labels 9 and 10). These profiles are consistent with the results of Bruenn & Haxton (1991), Thompson et al. (2003), and Rampp & Janka (2002).

Zoom In Zoom Out Reset image size

Now we move on to make the code comparison in the postbounce phase. Figure 12 compares the evolution of the shock radii for five models; two from our code (1DG15 (black line) and 3DG15 (red line)), one from AGILE-BOLTZTRAN (ABG15, blue line), and two from VERTEX with PROMETHEUS (VXG15, green line) and with CoCoNuT (BMG15, light blue line). Note that the final simulation time (${T}_{\mathrm{sim}}$) is 50 ms for model 1DG15 (black line), whereas ${T}_{\mathrm{sim}}$ is 32 ms for model 3DG15 (red line) simply limited by our available computational resources.

Zoom In Zoom Out Reset image size

The deviation seen in model 3DG15 (red line in Figure 12) from the rest of the 1D models is remarkable especially after ${T}_{\mathrm{pb}}\sim 5$ ms. This is because the bounce shock expands more energetically in 3D pushed primarily by prompt convection behind the shock. Using the same progenitor, Müller et al. (2012b) showed that the average shock radius becomes larger in 2D (their model G15) than in 1D (their model G15–1D) at ${T}_{\mathrm{pb}}\sim 70$ ms because of the hot-bubble convection starts which is seeded during the deceleration of the prompt shock. The Cartesian coordinates have intrinsic quadrupole perturbation and it significantly affects the growth of the prompt convection. The postbounce time, when the significant multidimensionality appears, thus differs between our model and the ones using the spherical polar coordinates. In addition, our coarse numerical resolution might also lead to the earlier appearance of the initial convection because of even larger seed perturbations.

The larger shock radius in model 3DG15 than that in 1DG15 is also consistent with our previous result with the leakage scheme (Kuroda et al. 2012; see also Hanke et al. 2012; Couch 2013for extensive discussions about the dimensional dependence on the postbounce dynamics). Now let us focus on our (pseudo-) 1D model (black line in Figure 12). The shock radius of our code is in good agreement with the reference results exceptionally before ${T}_{\mathrm{pb}}\lesssim 20$ ms, whereas the shock radius tends to be smaller until the end of the simulation time. We consider that the difference could primarily come from the use of the Cartesian coordinates with low numerical resolution and not from the neutrino–matter interaction terms. This is because our calculation in 1D spherical coordinates with using the same neutrino–matter interaction terms shows a good agreement in the shock evolution (see Appendix B). We give a more detailed discussion elsewhere below.

In Figures 13 and 14, we show various quantities from our code (labeled by "$\langle 1\mathrm{DG}15\rangle$" and "$\langle 3\mathrm{DG}15\rangle$" only in Figure 14), AGILE-BOLTZTRAN ("ABG15"), and VERTEX-PROMETHEUS ("VXG15") at two different postbounce times.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At 3 ms after bounce (Figure 13), the (angle-averaged) radial position of the stalled shock (bottom left panel) is $R\sim 70$ km for model 1DG15 (thick solid line). As seen, the velocity profile matches more closely the profile from AGILE-BOLTZTRAN (ABG15) than from VERTEX-PROMETHEUS (VXG15). This is also the case for the profiles of the density (top left panel), entropy (top right panel), and Y_e (bottom right panel). It should be noted that the more recent results from the VERTEX-PROMETHEUS code with an improved GR potential (Marek et al. 2006) agree very well with the AGILE-BOLTZTRAN code, hence with our code. Therefore, we mainly compare to model ABG15 in the following.

Looking at Figure 13 more closely, one can see that the profiles of our entropy (top right panel) and Y_e (bottom right panel) differ appreciably from model ABG15 especially in the region behind the stalled shock ($R\lesssim 70$ km) and above the unshocked inner core ($R\gtrsim 10$ km). Let us remark that the early postbounce evolution starts from the shock formation, followed by the emergence of the neutronization burst, until the shock stall is numerically most challenging. The code differs from the shock-capturing scheme as well as the treatment of GR, and the accuracy of the neutrino transport schemes could potentially impact the radiation-hydrodynamics evolution at the transient phase (i.e., 3–5 ms after bounce). Another remarkable difference is seen in the electron fraction. Among three results, only our result shows negative bump in Y_e profile at $10\lesssim R\lesssim 20$ km. However, this bump disappears in our 1D spherically symmetric test problem in which we solve the advection terms in energy space ${{\boldsymbol{S}}}_{{\rm{adv,e}}}$ explicitly in time (see Appendix B for more details).

Regarding the shock capturing, AGILE-BOLTZTRAN uses an artificial viscosity type with the second order accuracy in space, whereas our code employs the approximate Riemann solver (HLLE scheme like the VERTEX code) with second-order accuracy in space for both radiation-hydrodynamical and geometrical variables. Thus, the hydrodynamics part of our code is slightly more accurate than AGILE-BOLTZTRAN. On the other hand, the use of the approximate closure relation apparently falls behind the Boltzmann code especially in the semi-transparent region. Above all, the use of the Cartesian coordinates, which is very common in full GR simulations,¹¹ makes the comparison to the genuine "1D" results of the reference models (based on the 1D Lagrangian code (AGILE) and the multi-D code using the polar coordinates (VERTEX)) even more challenging.

At ${T}_{\mathrm{pb}}=50$ ms (Figure 14), the differences between our pseudo-1D model (1DG15, thick line) and the reference results still remain to be seen rather remarkably in the postshock region ($R\gtrsim 100$ km); however, this is not surprising given the different shock evolution (Figure 12). Here we consider that the numerical resolution in the postshock region sensitively affects the shock evolution. In the current resolution, the typical grid size of our nested box is ∼7.8 km at $120\lesssim R\lesssim 240$ km ($\sim 4^\circ$ resolution). As shown in Figure 12 (red line), the deviation of the shock radius from the reference models becomes remarkable at ${T}_{\mathrm{pb}}\gtrsim 20$ ms, which roughly coincides with the time when the shock reaches the coarser level of the nested grid. There the shock front is resolved only by a few grid cells. We consider that at least a factor of two or more higher resolution is required to reproduce 1D results, i.e., to recover the sphericity of the system in the Cartesian coordinates. However, this is not an easy task as we will discuss the code performance in Section 6.

From Figure 14, one can also see that the profile of the density, velocity, entropy, and Y_e at the central region ($R\lesssim 60$ km) agrees with the reference results roughly within an accuracy of ∼10%. Given the insufficient numerical resolution, this good agreement in the semi-transparent region would support the validity of the prescribed closure scheme, which is in line with the recent results by O'Connor (2014) and Just et al. (2015). Figure 14 also shows that the profile of model 3DG15 (thin solid line) differs from that of model 1DG15 (thick solid line).

The entropy profile for model 3DG15 (thin line) behind the shock ($R\lesssim 130$ km) becomes flatter compared to the 1D models, which is due to convective mixing behind the shock. As a result, the profiles of Y_e and entropy in model 3DG15 become slightly closer to the reference results, though the similarity is just a coincidence and is not meaningful.

Figure 15 shows the comparison of the neutrino luminosity L_ν which is the surface integral of the energy-integrated comoving energy flux ${H}_{(\varepsilon )}^{\mu }$ through the surface of cubic box with 1000 km width (i.e., approximately 500 km from the origin). The peak luminosity, ${L}_{{\nu }_{e},{\rm{peak}}}=4.1\times {10}^{53}$ erg s⁻¹, well agrees with $3.85\times {10}^{53}$ (ABG15), $3.8\times {10}^{53}$ (VXG15), and $4.3\times {10}^{53}$ (BMG15), respectively. After the neutronization phase when the 1D core enters to a quasi-static phase (${T}_{\mathrm{pb}}\gtrsim 20$ ms), the anti-electron and heavy lepton neutrino luminosities show quite consistent behaviors with the reference models and the differences between our results and, e.g., ABG15 are as small as ∼10%.

Zoom In Zoom Out Reset image size

Regarding the ${\nu }_{e}$ luminosity, we see a systematically lower value than the other three reference models. The difference reaches ∼30% (or ∼10⁵² erg s⁻¹) compared to ABG15 at ${T}_{\mathrm{pb}}=50$ ms. As a consequence, ${\nu }_{e}$ and ${\bar{\nu }}_{e}$ luminosities become comparable at ${T}_{\mathrm{pb}}\sim 45$ ms, which is significantly earlier than the reference values ${T}_{\mathrm{pb}}\sim 70$ ms in ABG15, BMG15, and also in a model of O'Connor (2014). This means that the lepton number loss from PNS is less than those previous studies since it can be measured roughly by the time integration of ${L}_{{\nu }_{e}}/{\varepsilon }_{{\nu }_{e}}-{L}_{{\bar{\nu }}_{e}}/{\varepsilon }_{{\bar{\nu }}_{e}}$. However, our 1D spherical test, albeit in the Newtonian limit, shown in Appendix B, does not exhibit such inconsistency and shows quite reasonable neutrino profiles with a previous study. We therefore consider that the reason for less agreement, especially seen in ${L}_{{\nu }_{e}}$, mainly comes from the spatial advection term in the Cartesian coordinates and not from the local neutrino matter interaction terms.

We also compare the neutrino spectral difference in Figure 16. In the figure, we show time evolutions of the rms energies of emergent neutrinos ${E}_{\nu ,{\rm{rms}}}$ measured at the surface of cubic box with 1000 km width. For ${E}_{\nu ,{\rm{rms}}}$, we use the same definition as in Liebendörfer et al. (2005) and it is defined as below:

$$\begin{eqnarray}&&{E}_{\nu ,{\rm{rms}}}\equiv \sqrt{\displaystyle \frac{\int d\varepsilon f(\nu ,\varepsilon ){\varepsilon }^{4}}{\int d\varepsilon f(\nu ,\varepsilon ){\varepsilon }^{2}}}.\end{eqnarray} \tag{ 74 }$$

Again in Figure 16, we plot averaged value $\langle {E}_{\nu ,{\rm{rms}}}\rangle$ over the surface of the cubic box. From the figure, we find good agreement with the reference codes and differences are less than ∼1 MeV for ${\nu }_{e}$ and ${\bar{\nu }}_{e}$ and $\lesssim 2\;{\rm{MeV}}$ for ${\nu }_{X}$ after the neutronization phase ceases ${T}_{\mathrm{pb}}\gtrsim 20$ ms. In our model 1DG15, we see a spurious second peak in the ${E}_{{\nu }_{X},{\rm{rms}}}$ profile around ${T}_{\mathrm{pb}}\sim 13$ ms. However, the second peak disappears in our model 3DG15 and we thus consider that it is most likely due to our artificial treatment for the matter velocity and is insignificant.

Zoom In Zoom Out Reset image size

We consider that neutrino absorption and emission: by free neutrons ${\nu }_{e}n\leftrightarrow {e}^{-}p$, by free protons ${\bar{\nu }}_{e}p\leftrightarrow {e}^{+}n$, and by heavy nuclei ${\nu }_{e}A\leftrightarrow {e}^{-}A^{\prime}$. The four vector source term ${S}^{\mathrm{nae},\mu }$ is given by

$$\begin{eqnarray}{S}_{(\nu ,\varepsilon )}^{\mathrm{nae},\mu } & = & {\kappa }_{(\nu ,\varepsilon )}^{\mathrm{nae}}\left[\left(\displaystyle \frac{4\pi {\varepsilon }^{3}}{{\rm{exp}}\{\beta (\varepsilon -{\mu }_{\nu })\}+1}-{J}_{(\nu ,\varepsilon )}\right){u}^{\mu }\right.\\ & & -{H}_{(\nu ,\varepsilon )}^{\mu }],\end{eqnarray} \tag{ 85 }$$

where ${\kappa }^{\mathrm{nae}}$ is the opacity, $\beta =1/{k}_{{\rm{B}}}T$ with ${k}_{{\rm{B}}}$ the Boltzmann's constant, and ${\mu }_{{\nu }_{e}}=-{\mu }_{{\bar{\nu }}_{e}}={\mu }_{e}-{\mu }_{p}+{\mu }_{n}$ is the chemical potential of neutrinos in thermal equilibrium with matter. Here, ${\mu }_{e}$, ${\mu }_{n}$, and ${\mu }_{p}$ are the chemical potentials of electrons, neutrons, and protons, including the rest mass energy of each particle, respectively. For each reaction (${\nu }_{e}n\leftrightarrow {e}^{-}p$, ${\bar{\nu }}_{e}p\leftrightarrow {e}^{+}n$ and ${\nu }_{e}A\leftrightarrow {e}^{-}A^{\prime}$), we first evaluate absorptivity $1/\lambda [{s}^{-1}]$ and then emissivity $j[{s}^{-1}]$.

Absorptivity for reaction ${\nu }_{e}n\leftrightarrow {e}^{-}p$ is expressed as (Bruenn 1985; Rampp & Janka 2002)

$$\begin{eqnarray}\displaystyle \frac{1}{{\lambda }_{\varepsilon }} & = & \displaystyle \frac{c}{{({\hslash }c)}^{4}}\displaystyle \frac{{G}_{F}^{2}}{\pi }({g}_{V}^{2}+3{g}_{A}^{2}){\eta }_{{np}}[1-{F}_{{e}^{-}}(\varepsilon +Q)]\\ & & \times (\varepsilon +Q)\sqrt{{(\varepsilon +Q)}^{2}-{m}_{e}^{2}{c}^{4}},\end{eqnarray} \tag{ 86 }$$

where we employed g_V = 1 and g_A = 1.23 as form factors resulting from the virtual strong interaction processes. ${F}_{x}(\varepsilon )\equiv {[1+{\rm{exp}}(\beta (\varepsilon -{\mu }_{x}))]}^{-1}$ is the Fermi distribution function of fermion x with energy ε. ${G}_{F}(=8.957\times {10}^{-44}$ MeV cm³) is the Fermi constant.

For the reaction ${\bar{\nu }}_{e}p\leftrightarrow {e}^{+}n$,

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\lambda }_{\varepsilon }}\;\\ &&=\;\left\{\begin{array}{l}\begin{array}{l}\displaystyle \frac{c}{{({\hslash }c)}^{4}}\displaystyle \frac{{G}_{F}^{2}}{\pi }({g}_{V}^{2}+3{g}_{A}^{2}){\eta }_{{pn}}[1-{F}_{{e}^{+}}(\varepsilon -Q)]\\ \quad (\varepsilon -Q)\sqrt{{(\varepsilon -Q)}^{2}-{m}_{e}^{2}{c}^{4}},\end{array}\\ \qquad {\rm{for}}\ \varepsilon \geqslant {m}_{e}{c}^{2}+Q\\ 0,\\ \quad {\rm{otherwise,}}\end{array}\right.\end{eqnarray} \tag{ 87 }$$

where $Q={m}_{n}{c}^{2}-{m}_{p}{c}^{2}$ is the rest mass energy difference of the neutron and proton.

Finally, for the reaction ${\nu }_{e}A\leftrightarrow {e}^{-}A^{\prime}$,

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\lambda }_{\varepsilon }}\;\\ &&=\;\left\{\begin{array}{l}\begin{array}{l}{e}^{-\beta {\rm{\Delta }}}\displaystyle \frac{c}{{({\hslash }c)}^{4}}\displaystyle \frac{{G}_{F}^{2}}{\pi }{g}_{A}^{2}\displaystyle \frac{2}{7}{N}_{p}{N}_{h}{n}_{A}[1-{F}_{{e}^{-}}(\varepsilon +Q^{\prime} )]\\ \quad (\varepsilon +Q^{\prime} )\sqrt{{(\varepsilon +Q^{\prime} )}^{2}-{m}_{e}^{2}{c}^{4}},\end{array}\\ \qquad {\rm{for}}\ \varepsilon \geqslant {m}_{e}{c}^{2}-Q^{\prime} \\ 0,\\ \quad {\rm{otherwise,}}\end{array}\right.\end{eqnarray} \tag{ 88 }$$

where $Q^{\prime} \equiv {\mu }_{n}-{\mu }_{p}+{\rm{\Delta }}$ is the mass difference between the initial and the final states through the reaction. By following Bruenn (1985) and Rampp & Janka (2002), we employed ${\rm{\Delta }}=3\;{\rm{MeV}}$ and

$$\begin{eqnarray}{N}_{p}=\left\{\begin{array}{ll}0, & {\rm{for}}\ Z\lt 20\\ Z-20, & {\rm{for}}\ 20\lt Z\lt 28\\ 8, & {\rm{for}}\ Z\gt 28,\end{array}\right.\end{eqnarray} \tag{ 89 }$$

$$\begin{eqnarray}{N}_{h}=\left\{\begin{array}{ll}6, & {\rm{for}}\ N\lt 34\\ 40-N, & {\rm{for}}\ 34\lt N\lt 40\\ 0, & {\rm{for}}\ N\gt 40,\end{array}\right.\end{eqnarray} \tag{ 90 }$$

for the number of protons N_p and holes N_h in the dominant GT resonance for the electron capture.

As for the values ${\eta }_{{pn}}$ and ${\eta }_{{np}}$, we adopted those proposed in Bruenn (1985) and Rampp & Janka (2002),

$$\begin{eqnarray}&&{\eta }_{{pn}}=\displaystyle \frac{{n}_{n}-{n}_{p}}{{\rm{exp}}(\beta ({\mu }_{n}-{\mu }_{p}-Q))-1},\end{eqnarray} \tag{ 91 }$$

$$\begin{eqnarray}&&{\eta }_{{np}}=\displaystyle \frac{{n}_{p}-{n}_{n}}{{\rm{exp}}(-\beta ({\mu }_{n}-{\mu }_{p}-Q))-1},\end{eqnarray} \tag{ 92 }$$

while we set

$$\begin{eqnarray}&&{\eta }_{{pn}}={n}_{p},\end{eqnarray} \tag{ 93 }$$

$$\begin{eqnarray}&&{\eta }_{{np}}={n}_{n},\end{eqnarray} \tag{ 94 }$$

in the non-degeneracy regime where ${\mu }_{n}-{\mu }_{p}-Q\lt 0.01\;{\rm{MeV}}$ is met.

Once we evaluate the absorptivity, the emissivity j and the opacity κ are obtained as (Bruenn 1985)

$$\begin{eqnarray}&&{j}_{\varepsilon }={e}^{\beta (\varepsilon -{\mu }_{\nu })}\displaystyle \frac{1}{{\lambda }_{\varepsilon }},\end{eqnarray} \tag{ 95 }$$

and

$$\begin{eqnarray}&&{\kappa }_{(\varepsilon )}^{\mathrm{nae}}={j}_{\varepsilon }+\displaystyle \frac{1}{{\lambda }_{\varepsilon }}.\end{eqnarray} \tag{ 96 }$$

The four vector source term for isoenergetic scattering of neutrinos off free nucleons and heavy nuclei is written as (Shibata et al. 2011)

$$\begin{eqnarray}&&{S}_{(\varepsilon )}^{\mathrm{iso},\mu }=-{\chi }_{(\varepsilon )}^{\mathrm{iso}}{H}_{(\varepsilon )}^{\mu }.\end{eqnarray} \tag{ 97 }$$

To evaluate ${\chi }_{(\varepsilon )}^{\mathrm{iso}}$ from the isoenergetic scattering kernel ${R}^{\mathrm{iso}}(\varepsilon ,\omega )$, we expand it into a Legendre series in terms of ω (here ω is the cosine of the scattering angle) up to the first order as

$$\begin{eqnarray}&&{R}^{\mathrm{iso}}(\varepsilon ,\omega )\approx \displaystyle \frac{1}{2}{{\rm{\Phi }}}_{(\varepsilon )}^{\mathrm{iso},0}+\displaystyle \frac{3}{2}\omega {{\rm{\Phi }}}_{(\varepsilon )}^{\mathrm{iso},1},\end{eqnarray} \tag{ 98 }$$

where ${{\rm{\Phi }}}_{(\varepsilon )}^{\mathrm{iso},0}$ and ${{\rm{\Phi }}}_{(\varepsilon )}^{\mathrm{iso},1}$ are the zeroth and first order of scattering kernels, respectively. After angular integration of the angular-dependent source term with respect to ω, ${\chi }_{(\varepsilon )}^{\mathrm{iso}}$ is expressed as below (Shibata et al. 2011):

$$\begin{eqnarray}&&{\chi }_{(\varepsilon )}^{\mathrm{iso}}={\varepsilon }^{2}({{\rm{\Phi }}}_{(\varepsilon )}^{\mathrm{iso},0}-{{\rm{\Phi }}}_{(\varepsilon )}^{\mathrm{iso},1}).\end{eqnarray} \tag{ 99 }$$

For the scattering process on the free nucleon (n/p), the zeroth- and first-order kernels become

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{(\varepsilon )}^{\mathrm{iso},0}=\displaystyle \frac{2\pi }{{({hc})}^{3}}\displaystyle \frac{{G}_{F}^{2}}{2\pi h}{\eta }_{{NN}}({({h}_{V}^{N})}^{2}+3{({h}_{A}^{N})}^{2}),\end{eqnarray} \tag{ 100 }$$

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{(\varepsilon )}^{\mathrm{iso},1}=\displaystyle \frac{2\pi }{{({hc})}^{3}}\displaystyle \frac{{G}_{F}^{2}}{2\pi h}{\eta }_{{NN}}({({h}_{V}^{N})}^{2}-{({h}_{A}^{N})}^{2}).\end{eqnarray} \tag{ 101 }$$

Obviously, ${{\rm{\Phi }}}_{(\varepsilon )}^{\mathrm{iso},0/1}$ has a dimension MeV⁻² s⁻¹ and ${\chi }_{(\varepsilon )}^{\mathrm{iso}}$ thus has a dimension s⁻¹. In the above, N takes n (neutron) or p (proton) and h_V^N and h_A^N are defined as (Bruenn 1985)

$$\begin{eqnarray}&&{h}_{V}^{n}=-\displaystyle \frac{1}{2},\end{eqnarray} \tag{ 102 }$$

$$\begin{eqnarray}&&{h}_{V}^{p}=\displaystyle \frac{1}{2}-2{{\rm{sin}}}^{2}{\theta }_{W},\end{eqnarray} \tag{ 103 }$$

$$\begin{eqnarray}&&{h}_{A}^{n}=-\displaystyle \frac{1}{2}{g}_{A},\end{eqnarray} \tag{ 104 }$$

$$\begin{eqnarray}&&{h}_{A}^{p}=\displaystyle \frac{1}{2}{g}_{A},\end{eqnarray} \tag{ 105 }$$

where ${\theta }_{W}$ is the Weinberg angle and we adopt the value ${{\rm{sin}}}^{2}{\theta }_{W}=0.2325$. By following Mezzacappa & Bruenn (1993c) and Rampp & Janka (2002), we evaluate ${\eta }_{{NN}}$ as

$$\begin{eqnarray}&&{\eta }_{{NN}}={n}_{N}\displaystyle \frac{{\xi }_{N}}{\sqrt{1+{\xi }_{N}^{2}}},\end{eqnarray} \tag{ 106 }$$

$$\begin{eqnarray}&&{\xi }_{N}=\displaystyle \frac{3}{2\beta {E}_{N}^{F}},\end{eqnarray} \tag{ 107 }$$

$$\begin{eqnarray}&&{E}_{N}^{F}=\displaystyle \frac{{({\hslash }c)}^{2}}{2{m}_{N}{c}^{2}}{(3{\pi }^{2}{n}_{N})}^{2/3}.\end{eqnarray} \tag{ 108 }$$

Next, for coherent scattering on heavy nuclei, the zeroth- and first-order kernels become (Bruenn 1985)

$$\begin{eqnarray}{{\rm{\Phi }}}_{(\varepsilon )}^{\mathrm{iso},0} & = & \displaystyle \frac{2\pi }{{({hc})}^{3}}\displaystyle \frac{{G}_{F}^{2}}{4\pi h}{n}_{A}{\left({{AC}}_{V}^{0}-\left(\displaystyle \frac{1}{2}A-Z\right){C}_{V}^{1}\right)}^{2}\\ & & \times \left(\displaystyle \frac{2y-1+{e}^{-2y}}{{y}^{2}}\right),\end{eqnarray} \tag{ 109 }$$

$$\begin{eqnarray}{{\rm{\Phi }}}_{(\varepsilon )}^{\mathrm{iso},1} & = & \displaystyle \frac{2\pi }{{({hc})}^{3}}\displaystyle \frac{{G}_{F}^{2}}{4\pi h}{n}_{A}{\left({{AC}}_{V}^{0}-\left(\displaystyle \frac{1}{2}A-Z\right){C}_{V}^{1}\right)}^{2}\\ & & \times \left(\displaystyle \frac{2-3y+2{y}^{2}-(2+y){e}^{-2y}}{{y}^{3}}\right),\end{eqnarray} \tag{ 110 }$$

with ${C}_{V/A}^{0}=({h}_{V/A}^{p}+{h}_{V/A}^{n})/2$, ${C}_{V/A}^{1}=({h}_{V/A}^{p}-{h}_{V/A}^{n})$, $y=4b{\varepsilon }^{2}$, and $b=4.8\times {10}^{-6}{A}^{2/3}$. Here n_A, A, and Z denote the number density of heavy nuclei, the atomic number, and the charge, respectively.

Inelastic scattering of neutrinos off electrons plays important role to help neutrinos escape more freely from the center as a consequence of down-scattering and it thus enhances deleptonization of central core. The source term ${S}_{(\varepsilon )}^{\mathrm{nes},\mu }$ is expressed in terms of the collision integral ${B}_{(\varepsilon ,{\rm{\Omega }})}^{\mathrm{nes}}$ as (Shibata et al. 2011)

$$\begin{eqnarray}&&{S}_{(\varepsilon )}^{\mathrm{nes},\mu }={\varepsilon }^{3}\int d{\rm{\Omega }}{B}_{(\varepsilon ,{\rm{\Omega }})}^{\mathrm{nes}}({u}^{\mu }+{l}^{\mu }),\end{eqnarray} \tag{ 111 }$$

where ${l}^{\mu }$ is a unit normal four vector orthogonal to ${u}^{\mu }$. The collision integral ${B}_{(\varepsilon ,{\rm{\Omega }})}^{\mathrm{nes}}$ along the propagation direction Ω is expressed as

$$\begin{eqnarray}{B}_{(\varepsilon ,{\rm{\Omega }})}^{\mathrm{nes}} & = & \displaystyle \int \varepsilon {^{\prime} }^{2}d\varepsilon ^{\prime} d{\rm{\Omega }}^{\prime} [f(\varepsilon ^{\prime} ,{\rm{\Omega }}^{\prime} )\{1-f(\varepsilon ,{\rm{\Omega }})\}{R}^{\mathrm{nes},\mathrm{in}}(\varepsilon ,\varepsilon ^{\prime} ,\omega )\\ & & -\{1-f(\varepsilon ^{\prime} ,{\rm{\Omega }}^{\prime} )\}f(\varepsilon ,{\rm{\Omega }}){R}^{\mathrm{nes},\mathrm{out}}(\varepsilon ,\varepsilon ^{\prime} ,\omega )],\end{eqnarray} \tag{ 112 }$$

where ${R}^{\mathrm{nes},\mathrm{in}/\mathrm{out}}(\varepsilon ,\varepsilon ^{\prime} ,\omega )$ is the angular-dependent inward/outward scattering kernel and ω is the cosine of scattering angle, i.e., angle between Ω and ${\rm{\Omega }}^{\prime}$. With expanding the angular-dependent inward/outward scattering kernel ${R}^{\mathrm{nes},\mathrm{in}/\mathrm{out}}(\varepsilon ,\varepsilon ^{\prime} ,\omega )$ into a Legendre series with respect to ω and taking up to the first order as

$$\begin{eqnarray}&&{R}^{\mathrm{nes},\mathrm{in}/\mathrm{out}}(\varepsilon ,\varepsilon ^{\prime} ,\omega )\approx \displaystyle \frac{1}{2}{{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{\mathrm{nes},\mathrm{in}/\mathrm{out},0}+\displaystyle \frac{3}{2}\omega {{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{\mathrm{nes},\mathrm{in}/\mathrm{out},1},\end{eqnarray} \tag{ 113 }$$

and with decomposing the neutrino distribution function after scattering into isotropic and non-isotropic parts as

$$\begin{eqnarray}&&f(\varepsilon ^{\prime} ,{\rm{\Omega }}^{\prime} )\approx {f}^{0}(\varepsilon ^{\prime} )+{f}^{1,\mu }(\varepsilon ^{\prime} ){l}_{\mu }^{\prime },\end{eqnarray} \tag{ 114 }$$

the final expression of the scattering integral is described as 

$$\begin{eqnarray}&&{B}_{(\varepsilon ,{\rm{\Omega }})}^{\mathrm{nes}}=f(\varepsilon ,{\rm{\Omega }})({A}_{(\varepsilon ),\mathrm{nes}}^{0}+{B}_{(\varepsilon ),\mathrm{nes}}^{0,\mu }{l}_{\mu })+{C}_{(\varepsilon ),\mathrm{nes}}^{0}+{C}_{(\varepsilon ),\mathrm{nes}}^{1,\mu }{l}_{\mu }.\end{eqnarray} \tag{ 115 }$$

In the above, we used the same notations used in Bruenn (1985; see their Equations (108)–(113)) for ${A}_{(\varepsilon ),\mathrm{nes}}^{0}$, ${B}_{(\varepsilon ),\mathrm{nes}}^{0,\mu }$, ${C}_{(\varepsilon ),\mathrm{nes}}^{0}$, and ${C}_{(\varepsilon ),\mathrm{nes}}^{1,\mu }$ with the exception of ${B}_{(\varepsilon ),\mathrm{nes}}^{0,\mu }$ and ${C}_{(\varepsilon ),\mathrm{nes}}^{1,\mu }$, which have three spatial components (the zeroth component is determined from orthogonality ${H}^{\mu }{u}_{\mu }=0$) and of ${B}_{(\varepsilon ),\mathrm{nes}}^{0,\mu }$, which is a factor of 3 larger than that of Equation (109) in Bruenn (1985). They are explicitly written as

$$\begin{eqnarray}{A}_{(\varepsilon ),\mathrm{nes}}^{0} & = & -\displaystyle \frac{2\pi }{{({hc})}^{3}}\displaystyle \int \varepsilon {^{\prime} }^{2}d\varepsilon ^{\prime} [{f}^{0}(\varepsilon ^{\prime} ){{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{\mathrm{nes},\mathrm{in},0}\\ & & +(1-{f}^{0}(\varepsilon ^{\prime} )){{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{\mathrm{nes},\mathrm{out},0}],\end{eqnarray} \tag{ 116 }$$

$$\begin{eqnarray}&&{B}_{(\varepsilon ),\mathrm{nes}}^{0,\mu }=-\displaystyle \frac{2\pi }{{({hc})}^{3}}\int \varepsilon {^{\prime} }^{2}d\varepsilon ^{\prime} {f}^{1,\mu }(\varepsilon ^{\prime} )({{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{\mathrm{nes},\mathrm{in},1}-{{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{\mathrm{nes},\mathrm{out},1}),\end{eqnarray} \tag{ 117 }$$

$$\begin{eqnarray}&&{C}_{(\varepsilon ),\mathrm{nes}}^{0}=\displaystyle \frac{2\pi }{{({hc})}^{3}}\int \varepsilon {^{\prime} }^{2}d\varepsilon ^{\prime} {f}^{0}(\varepsilon ^{\prime} ){{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{\mathrm{nes},\mathrm{in},0},\end{eqnarray} \tag{ 118 }$$

$$\begin{eqnarray}&&{C}_{(\varepsilon ),\mathrm{nes}}^{1,\mu }=\displaystyle \frac{2\pi }{{({hc})}^{3}}\int \varepsilon {^{\prime} }^{2}d\varepsilon ^{\prime} {f}^{1,\mu }(\varepsilon ^{\prime} ){{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{\mathrm{nes},\mathrm{in},1},\end{eqnarray} \tag{ 119 }$$

where we employ

$$\begin{eqnarray}&&{f}^{0}(\varepsilon )=\displaystyle \frac{{J}_{(\varepsilon )}}{4\pi {\varepsilon }^{3}},\end{eqnarray} \tag{ 120 }$$

$$\begin{eqnarray}&&{f}^{1,\mu }(\varepsilon )=\displaystyle \frac{3{H}_{(\varepsilon )}^{\mu }}{4\pi {\varepsilon }^{3}},\end{eqnarray} \tag{ 121 }$$

for the isotropic and non-isotropic parts of the distribution function (see also Shibata et al. 2011). For an explicit evaluation for ${{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{{\rm{nes,in/out}},0/1}$, we refer the reader to Yueh & Buchler (1976) and Bruenn (1985).

Since coefficients (116)–(119) do not depend on the propagation direction Ω, the final integral Equation (111) with imposing Equation (115) becomes

$$\begin{eqnarray}{S}_{(\varepsilon )}^{\mathrm{nes},\mu } & = & ({J}_{(\varepsilon )}{u}^{\mu }+{H}_{(\varepsilon )}^{\mu }){A}_{(\varepsilon ),\mathrm{nes}}^{0}\\ & & +({H}_{\alpha (\varepsilon )}{u}^{\mu }+{{L}_{\alpha }^{\mu }}_{(\varepsilon )}){B}_{(\varepsilon ),\mathrm{nes}}^{0,\alpha }\\ & & +4\pi {\varepsilon }^{3}{u}^{\mu }{C}_{(\varepsilon ),\mathrm{nes}}^{0}\\ & & +\displaystyle \frac{4\pi {\varepsilon }^{3}}{3}{{h}_{\alpha }}^{\mu }{C}_{(\varepsilon ),\mathrm{nes}}^{1,\alpha },\end{eqnarray} \tag{ 122 }$$

where ${h}_{\mu \nu }\equiv {g}_{\mu \nu }+{u}_{\mu }{u}_{\nu }$ is the projection operator. Since ${{\rm{\Phi }}}^{\mathrm{nes}}$ has a unit [cm³ s⁻¹], ${A}_{(\varepsilon ),\mathrm{nes}}^{0}$, ${B}_{(\varepsilon ),\mathrm{nes}}^{0,\mu }$, ${C}_{(\varepsilon ),\mathrm{nes}}^{0}$, and ${C}_{(\varepsilon ),\mathrm{nes}}^{1,\mu }$ thus have a unit [s⁻¹] which is required to let the dimension of ${S}_{(\varepsilon )}^{\mathrm{nes},\mu }$ to be [MeV³ s⁻¹].

As for the thermal pair production/annihilation process of neutrinos, we take the same approach as the neutrino electron scattering process. The collision integral ${B}_{(\varepsilon ,{\rm{\Omega }})}^{\mathrm{tp}}$ is expressed as

$$\begin{eqnarray}{B}_{(\varepsilon ,{\rm{\Omega }})}^{\mathrm{nes}} & = & \displaystyle \int \varepsilon {^{\prime} }^{2}d\varepsilon ^{\prime} d{\rm{\Omega }}^{\prime} [\{1-\bar{f}(\varepsilon ^{\prime} ,{\rm{\Omega }}^{\prime} )\}\\ & & \times \{1-f(\varepsilon ,{\rm{\Omega }})\}{R}^{\mathrm{tp},\mathrm{pro}}(\varepsilon ,\varepsilon ^{\prime} ,\omega )\\ & & -\bar{f}(\varepsilon ^{\prime} ,{\rm{\Omega }}^{\prime} )f(\varepsilon ,{\rm{\Omega }}){R}^{\mathrm{tp},\mathrm{ann}}(\varepsilon ,\varepsilon ^{\prime} ,\omega )],\end{eqnarray} \tag{ 123 }$$

where $\bar{f}$ denotes the anti-neutrino distribution function and ${R}^{\mathrm{tp},\mathrm{pro}/\mathrm{ann}}(\varepsilon ,\varepsilon ^{\prime} ,\omega )$ is the angular-dependent production/annihilation kernel. We again expand the kernel ${R}^{\mathrm{tp},\mathrm{pro}/\mathrm{ann}}(\varepsilon ,\varepsilon ^{\prime} ,\omega )$ up to the first order in ω as

$$\begin{eqnarray}&&{R}^{\mathrm{tp},\mathrm{pro}/\mathrm{ann}}(\varepsilon ,\varepsilon ^{\prime} ,\omega )\approx \displaystyle \frac{1}{2}{{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{\mathrm{tp},\mathrm{pro}/\mathrm{ann},0}+\displaystyle \frac{3}{2}\omega {{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{\mathrm{tp},\mathrm{pro}/\mathrm{ann},1},\end{eqnarray} \tag{ 124 }$$

and the final expression of the scattering integral is thus described as 

$$\begin{eqnarray}&&{B}_{(\varepsilon ,{\rm{\Omega }})}^{\mathrm{tp}}=f(\varepsilon ,{\rm{\Omega }})({A}_{(\varepsilon ),\mathrm{tp}}^{0}+{B}_{(\varepsilon ),\mathrm{tp}}^{0,\mu }{l}_{\mu })+{C}_{(\varepsilon ),\mathrm{tp}}^{0}+{C}_{(\varepsilon ),\mathrm{tp}}^{1,\mu }{l}_{\mu },\end{eqnarray} \tag{ 125 }$$

which has exactly the same expression as that of the neutrino electron scattering (Equation (115)). Again, coefficients have the same notations as in Bruenn (1985; see their Equations (117)–(121)) and described as

$$\begin{eqnarray}{A}_{(\varepsilon ),\mathrm{tp}}^{0} & = & -\displaystyle \frac{2\pi }{{({hc})}^{3}}\displaystyle \int \varepsilon {^{\prime} }^{2}d\varepsilon ^{\prime} [\{1-{\bar{f}}^{0}(\varepsilon ^{\prime} )\}{{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{\mathrm{tp},\mathrm{pro},0}\\ & & +{\bar{f}}^{0}(\varepsilon ^{\prime} ){{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{\mathrm{tp},\mathrm{ann},0}],\end{eqnarray} \tag{ 126 }$$

$$\begin{eqnarray}&&{B}_{(\varepsilon ),\mathrm{tp}}^{0,\mu }=\displaystyle \frac{2\pi }{{({hc})}^{3}}\int \varepsilon {^{\prime} }^{2}d\varepsilon ^{\prime} {\bar{f}}^{1,\mu }(\varepsilon ^{\prime} )({{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{\mathrm{tp},\mathrm{pro},1}-{{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{\mathrm{tp},\mathrm{ann},1}),\end{eqnarray} \tag{ 127 }$$

$$\begin{eqnarray}&&{C}_{(\varepsilon ),\mathrm{tp}}^{0}=\displaystyle \frac{2\pi }{{({hc})}^{3}}\int \varepsilon {^{\prime} }^{2}d\varepsilon ^{\prime} \{1-{\bar{f}}^{0}(\varepsilon ^{\prime} )\}{{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{\mathrm{tp},\mathrm{pro},0},\end{eqnarray} \tag{ 128 }$$

$$\begin{eqnarray}&&{C}_{(\varepsilon ),\mathrm{tp}}^{1,\mu }=-\displaystyle \frac{2\pi }{{({hc})}^{3}}\int \varepsilon {^{\prime} }^{2}d\varepsilon ^{\prime} {\bar{f}}^{1,\mu }(\varepsilon ^{\prime} ){{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{\mathrm{tp},\mathrm{pro},1}.\end{eqnarray} \tag{ 129 }$$

For an explicit evaluation for ${{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{{\rm{tp,pro/ann}},0/1}$, we refer the reader to Bruenn (1985). The final expression of the source term ${S}_{(\varepsilon )}^{\mathrm{tp},\mu }$ is simply obtained by replacing coefficients in Equation (122) with those of thermal process, i.e., with Equations (126)–(129).

The collision integral for Nucleon-nucleon Bremsstrahlung ${B}_{(\varepsilon ,{\rm{\Omega }})}^{\mathrm{br}}$ has the same expression as that of the thermal pair production/annihilation of neutrinos and therefore the same notations used in Appendix A.4 can be directly applicable. ${B}_{(\varepsilon ,{\rm{\Omega }})}^{\mathrm{br}}$ is written as

$$\begin{eqnarray}{B}_{(\varepsilon ,{\rm{\Omega }})}^{\mathrm{br}} & = & \displaystyle \int \varepsilon {^{\prime} }^{2}d\varepsilon ^{\prime} d{\rm{\Omega }}^{\prime} [\{1-\bar{f}(\varepsilon ^{\prime} ,{\rm{\Omega }}^{\prime} )\}\\ & & \times \{1-f(\varepsilon ,{\rm{\Omega }})\}{R}^{\mathrm{br},\mathrm{pro}}(\varepsilon ,\varepsilon ^{\prime} ,\omega )\\ & & -\bar{f}(\varepsilon ^{\prime} ,{\rm{\Omega }}^{\prime} )f(\varepsilon ,{\rm{\Omega }}){R}^{\mathrm{br},\mathrm{ann}}(\varepsilon ,\varepsilon ^{\prime} ,\omega )],\end{eqnarray} \tag{ 130 }$$

where ${R}^{\mathrm{br},\mathrm{pro}/\mathrm{ann}}(\varepsilon ,\varepsilon ^{\prime} ,\omega )$ is the production/annihilation kernel. As for the Bremsstrahlung process, we again expand the kernel ${R}^{\mathrm{br},\mathrm{pro}/\mathrm{ann}}(\varepsilon ,\varepsilon ^{\prime} ,\omega )$ into a Legendre series, but take only the zeroth order term in ω for simplicity. Then the kernel becomes

$$\begin{eqnarray}&&{R}^{\mathrm{br},\mathrm{pro}/\mathrm{ann}}(\varepsilon ,\varepsilon ^{\prime} ,\omega )\approx \displaystyle \frac{1}{2}{{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{\mathrm{br},\mathrm{pro}/\mathrm{ann},0}.\end{eqnarray} \tag{ 131 }$$

Coefficients used in the final expression for the collision integral and the source term are simply evaluated by replacing ${{\rm{\Phi }}}^{\mathrm{tp}}$ in Equations (126)–(129) with ${{\rm{\Phi }}}^{\mathrm{br}}$.

Isotropic production kernel ${{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{\mathrm{br},\mathrm{pro},0}$ can be evaluated by following manner.

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{(\varepsilon ,\varepsilon ^{\prime} )}^{\mathrm{br},\mathrm{pro},0}\equiv {e}^{-\beta (\varepsilon +\varepsilon ^{\prime} )}\displaystyle \sum _{i\in {\rm{nn,pp,np}}}{\rm{min}}({\phi }_{\mathrm{br}}^{{\rm{D}},i}(\varepsilon ,\varepsilon ^{\prime} ),{\phi }_{\mathrm{br}}^{\mathrm{ND},i}(\varepsilon ,\varepsilon ^{\prime} )).\end{eqnarray} \tag{ 132 }$$

In the above equation, indexes D and ND denote the degenerate and non-degenerate limit of free nucleons, respectively, and (nn, pp, np) represent bremsstrahlung due to neutron–neutron, proton–proton, and neutron–proton pair, respectively. ${\phi }_{\mathrm{br}}^{{\rm{D}},i}(\varepsilon ,\varepsilon ^{\prime} )$ and ${\phi }_{\mathrm{br}}^{\mathrm{ND},i}(\varepsilon ,\varepsilon ^{\prime} )$ are expressed as

$$\begin{eqnarray}{\phi }_{\mathrm{br}}^{{\rm{D}},i}(\varepsilon ,\varepsilon ^{\prime} ) & = & {({\hslash }c)}^{6}\displaystyle \frac{{G}^{2}}{4}{\left(\displaystyle \frac{2\pi f}{{m}_{\pi }}\right)}^{4}\\ & & \times \displaystyle \frac{{m}_{i}^{4}}{12{\pi }^{9}}\displaystyle \frac{1}{{\beta }^{2}(\varepsilon +\varepsilon ^{\prime} )}\displaystyle \frac{(4{\pi }^{2}+{\beta }^{2}{(\varepsilon +\varepsilon ^{\prime} )}^{2})}{1-{e}^{-\beta (\varepsilon +\varepsilon ^{\prime} )}}\displaystyle \frac{{\alpha }_{i}}{{S}_{i}}\\ & & \times \left[\displaystyle \frac{2{\hslash }c{(3{\pi }^{2}{X}_{i}{n}_{b})}^{\frac{1}{3}}}{{({\hslash }c)}^{6}}\right],\end{eqnarray} \tag{ 133 }$$

$$\begin{eqnarray}{\phi }_{\mathrm{br}}^{\mathrm{ND},i}(\varepsilon ,\varepsilon ^{\prime} ) & = & {({\hslash }c)}^{6}\displaystyle \frac{{G}^{2}}{4}{\left(\displaystyle \frac{2\pi f}{{m}_{\pi }}\right)}^{4}\\ & & \times \displaystyle \frac{2{m}_{i}^{3/2}}{{\pi }^{11/2}}\displaystyle \frac{1}{{\beta }^{1/2}{(\varepsilon +\varepsilon ^{\prime} )}^{2}}\sqrt{1+\beta (\varepsilon +\varepsilon ^{\prime} )}\\ & & \times \displaystyle \frac{{\alpha }_{i}}{{S}_{i}}{X}_{i}^{2}{n}_{b}^{2}.\end{eqnarray} \tag{ 134 }$$

In the above equations, ${\phi }_{\mathrm{br}}^{{\rm{D}}/\mathrm{ND},i}(\varepsilon ,\varepsilon ^{\prime} )$ have a dimension in [cm³ s⁻¹], f = 1 and ${G}^{2}=c{({\hslash }c)}^{2}{G}_{F}^{2}\quad =\quad 1.55\quad \times \quad {10}^{-33}$ [MeV⁻² cm³ s⁻¹]. We employed values for ${\alpha }_{i}$, X_i, m_i, and S_i as

$$\begin{eqnarray}{\alpha }_{i} & = & \left(\begin{array}{c}3{g}_{A}^{2}\\ 3{g}_{A}^{2}\\ \displaystyle \frac{7}{2}{g}_{A}^{2}\end{array}\right),\ {X}_{i}=\left(\begin{array}{c}{X}_{n}\\ {X}_{p}\\ {\rm{min}}({X}_{n},{X}_{p})\end{array}\right),\ {m}_{i}=\left(\begin{array}{c}{m}_{n}\\ {m}_{p}\\ \sqrt{{m}_{n}{m}_{p}}\end{array}\right),\\ {S}_{i} & = & \left(\begin{array}{c}4\\ 4\\ 1\end{array}\right),\ \mathrm{for}\ i=\left(\begin{array}{c}{nn}\\ {pp}\\ {np}\end{array}\right).\end{eqnarray} \tag{ 135 }$$

After deriving the production kernel, we can obtain the annihilation one by using a following relation

$$\begin{eqnarray}&&{R}^{\mathrm{br},\mathrm{ann}}(\varepsilon ,\varepsilon ^{\prime} ,\omega )={e}^{\beta (\varepsilon +\varepsilon ^{\prime} )}{R}^{\mathrm{br},\mathrm{pro}}(\varepsilon ,\varepsilon ^{\prime} ,\omega ).\end{eqnarray} \tag{ 136 }$$

By combining all the source terms mentioned above, the final expression becomes

$$\begin{eqnarray}{S}_{(\varepsilon )}^{\mu } & = & {\kappa }_{(\varepsilon )}^{\mathrm{nae}}\{({J}_{(\varepsilon )}^{\mathrm{eq}}-{J}_{(\varepsilon )}){u}^{\mu }-{H}_{(\varepsilon )}^{\mu }\}-{\chi }_{(\varepsilon )}^{\mathrm{iso}}{H}_{(\varepsilon )}^{\mu }\\ & & +\displaystyle \sum _{X\in {\rm{nes,tp,br}}}\{({J}_{(\varepsilon )}{u}^{\mu }+{H}_{(\varepsilon )}^{\mu }){A}_{(\varepsilon ),{\rm{X}}}^{0}+4\pi {\varepsilon }^{3}{u}^{\mu }{C}_{(\varepsilon ),{\rm{X}}}^{0}\\ & & +\left.({H}_{\alpha (\varepsilon )}{u}^{\mu }+{{L}_{\alpha }^{\mu }}_{(\varepsilon )}){B}_{(\varepsilon ),{\rm{X}}}^{0,\alpha }+\displaystyle \frac{4\pi {\varepsilon }^{3}}{3}{{h}_{\alpha }}^{\mu }{C}_{(\varepsilon ),{\rm{X}}}^{1,\alpha }\right\},\end{eqnarray} \tag{ 137 }$$

where ${J}_{(\varepsilon )}^{\mathrm{eq}}$ is written by

$$\begin{eqnarray}&&{J}_{(\varepsilon )}^{\mathrm{eq}}=\displaystyle \frac{4\pi {\varepsilon }^{3}}{{\rm{exp}}\{\beta (\varepsilon -{\mu }_{\nu })\}+1}.\end{eqnarray} \tag{ 138 }$$

In the end, we plot inverse mean free paths of each process for reference. Assumed hydrodynamical profiles are $\rho ={10}^{10}$ g cm⁻³, T = 0.638 MeV, and Y_e = 0.43 (Figure 17), which corresponds to the central profile of "s15s7b2" in Woosley & Weaver (1995), $\rho ={10}^{11}$ g cm⁻³, T = 1.382 MeV, and Y_e = 0.4 (Figure 18), which is the same condition used in Figure 36 in Bruenn (1985), and $\rho ={10}^{13}$ g cm⁻³, T = 3 MeV, and Y_e = 0.25 (Figure 19), which is a typical profile after neutrino trapping during collapse phase. As in Bruenn (1985), we assumed the final state occupancy of the neutrinos is zero. More explicitly, we plot $1/{\lambda }_{\mathrm{mfp}}$ [cm⁻¹] expressed by

$$\begin{eqnarray}&&{\lambda }_{\mathrm{mfp}}^{-1}\\ =\;\displaystyle \frac{1}{c}\left\{\begin{array}{cc}{\kappa }_{(\nu ,\varepsilon )}^{\mathrm{nae}}, & (\mathrm{neutrino}\ \mathrm{absorption}\ \mathrm{and}\ \mathrm{emission})\\ {\chi }_{(\varepsilon )}^{\mathrm{iso}}, & (\mathrm{isoenergy}\ \mathrm{scattering})\\ {A}_{(\varepsilon ),\mathrm{nes}}^{0}, & (\mathrm{neutrino}\ \mathrm{electron}\ \mathrm{scattering})\\ {A}_{(\varepsilon ),\mathrm{tp}}^{0}, & (\mathrm{thermal}\ \mathrm{pair}\ \mathrm{production}\ \mathrm{and}\ \mathrm{annihilation})\\ {A}_{(\varepsilon ),\mathrm{br}}^{0}, & (\mathrm{nucleon}\mbox{--}\mathrm{nucleon}\ \mathrm{bremsstrahlung}).\end{array}\right.\end{eqnarray} \tag{ 139 }$$

Since the isoenergy scattering rate has the same value for all neutrino species, we only plotted for electron-type neutrinos. Note that the inverse mean free path for bremsstrahlung is multiplied by 10¹⁵ (Figure 17), 10¹⁰ (Figure 18), and 10⁵ (Figure 19), and for the thermal process it is multiplied by 10¹⁰ (Figures 17–19).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size