Simulations of the Early Postbounce Phase of Core-collapse Supernovae in Three-dimensional Space with Full Boltzmann Neutrino Transport

The explosion mechanism of core-collapse supernovae has been studied for a long time (see Janka 2012; Burrows 2013; Janka et al. 2016; Müller 2016, for reviews). Although it is not completely understood yet, accurate simulations in one spatial dimension (1D) under spherical symmetry have shown that the delayed explosion scenario by neutrino heating fails for most of the progenitor models (e.g., Liebendörfer et al. 2001; Sumiyoshi et al. 2005), whereas multidimensional simulations with some approximate neutrino transport have brought successful explosions (e.g., Lentz et al. 2015; Bruenn et al. 2016; Pan et al. 2016; Takiwaki et al. 2016; Müller et al. 2017; Just et al. 2018; Vartanyan et al. 2019, for recent papers). However, it has not been clearly demonstrated whether the approximations for neutrino transport are really justified in the transition region between the optically thin and thick limits, particularly in the multidimensional settings. In fact, the turbulence energized by the neutrino-driven convection and the standing accretion shock instability (SASI) form complex structures in density, electron fraction, and temperature in the transition region. It is hence very interesting to see how the neutrinos are transported through such a chaotic environment.

Motivated by these arguments, we have developed the Boltzmann radiation-hydrodynamics code in recent years: the basic framework of our code is given in Sumiyoshi & Yamada (2012); the nonrelativistic Boltzmann equation with the standard weak interactions is solved on some fixed matter distributions derived from supernova simulations (Sumiyoshi et al. 2015); the Boltzmann neutrino transport solver is extended so that it could treat the special relativistic effects to all orders of v/c, where v and c are the speeds of matter and light, respectively; the Boltzmann solver is coupled to a hydrodynamics solver with Newtonian self-gravity (Nagakura et al. 2014); a moving-mesh technique to follow the proper motion of proto–neutron stars (PNSs) is implemented in the general relativistic Boltzmann solver in the 3+1 formalism (Nagakura et al. 2017); very recently a new method to improve the accuracy of the momentum feedback from neutrino to matter has been proposed by Nagakura et al. (2019c); and weak interactions for light nuclei are added to the code in Nagakura et al. (2019a).

These codes have been employed in spatially two-dimensional (2D) simulations under axisymmetry to study the equation-of-state (EOS) dependence (Nagakura et al. 2018) of and rotational effects (Harada et al. 2019) on dynamics, as well as possible early PNS kicks (Nagakura et al. 2019d) and fast neutrino flavor conversions (Delfan Azari et al. 2019; Nagakura et al. 2019b). In particular, Nagakura et al. (2018) and Harada et al. (2019) revealed nonaxisymmetric features in the neutrino distribution function in momentum space and some differences in the Eddington tensors between the Boltzmann simulations and the M1 closure approximation. In this paper, we present for the first time the results of radiation hydrodynamic simulations in three spatial dimensions (3D) without any symmetry with our Boltzmann radiation-hydrodynamics code. We pay particular attention to the neutrino angular distributions in the early postbounce phase when the prompt convection grows. We compare the Eddington tensor obtained in the Boltzmann simulations with those evaluated in the M1 closure approximation, employing a new analysis method.

This paper is organized as follows. We explain the basic equations, weak interactions, and numerical setup used in this paper in Section 2. We first describe the 3D features in matter flows and radiation fields, comparing them with the 1D and 2D counterparts in Section 3. We then conduct the new analysis with the eigenvalues and eigenvectors of the Eddington tensor in Section 4. We conclude this paper in Section 5. Throughout the paper, the metric signature is − + + +, and we use the units with c = G = h = 1 unless otherwise stated, where G and h are the gravitational constant and Planck's constant, respectively.

We use the Boltzmann radiation-hydrodynamics code for core-collapse simulations (Sumiyoshi & Yamada 2012; Nagakura et al. 2014). The coordinate system in phase space (r, θ, ϕ, , θ_ν, ϕ_ν) is shown in Figure 1, where r, θ, ϕ, , θ_ν, and ϕ_ν denote radius, polar and azimuthal angles in space, neutrino energy, and polar and azimuthal angles in momentum space, respectively. The special relativistic Boltzmann equation in the laboratory frame,

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial f}{\partial t}+\displaystyle \frac{{\mu }_{\nu }}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}({r}^{2}f)+\displaystyle \frac{\sqrt{1-{\mu }_{\nu }^{2}}\cos {\phi }_{\nu }}{r\sin \theta }\displaystyle \frac{\partial }{\partial \theta }(\sin \theta f)\\ \quad \ +\ \displaystyle \frac{\sqrt{1-{\mu }_{\nu }^{2}}\sin {\phi }_{\nu }}{r\sin \theta }\displaystyle \frac{\partial f}{\partial \phi }+\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial {\mu }_{\nu }}[(1-{\mu }_{\nu }^{2})f]\\ \quad \ -\ \displaystyle \frac{\sqrt{1-{\mu }_{\nu }^{2}}}{r}\displaystyle \frac{\cos \theta }{\sin \theta }\displaystyle \frac{\partial }{\partial {\phi }_{\nu }}(\sin {\phi }_{\nu }f)={\left(\displaystyle \frac{\delta f}{\delta t}\right)}_{\mathrm{col}},\end{array}\end{eqnarray} \tag{ 1 }$$

is solved by the discrete ordinate S_N method for three neutrino species: ν_e, ${\bar{\nu }}_{e}$, and ν_x; f and ${\mu }_{\nu }(=\cos {\theta }_{\nu })$ are the neutrino distribution function in phase space and cosine of the polar angle in momentum space, respectively. Using the finite-volume method, we also solve the Newtonian compressible hydrodynamics equations along with the time evolution equation of electron number density:

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{Q}}+{\partial }_{j}{{\boldsymbol{U}}}^{j}={{\boldsymbol{W}}}_{h}+{{\boldsymbol{W}}}_{i},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\boldsymbol{Q}}=\sqrt{g}\left(\begin{array}{c}\rho \\ \rho {v}_{r}\\ \rho {v}_{\theta }\\ \rho {v}_{\phi }\\ e+\tfrac{1}{2}\rho {v}^{2}\\ \rho {Y}_{e}\end{array}\right),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\boldsymbol{U}}}^{j}=\sqrt{g}\left(\begin{array}{c}\rho {v}^{j}\\ \rho {v}_{r}{v}^{j}+P{\delta }_{r}^{j}\\ \rho {v}_{\theta }{v}^{j}+P{\delta }_{\theta }^{j}\\ \rho {v}_{\phi }{v}^{j}+P{\delta }_{\phi }^{j}\\ \left(e+P+\tfrac{1}{2}\rho {v}^{2}\right){v}^{j}\\ \rho {Y}_{e}{v}^{j}\end{array}\right),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{{\boldsymbol{W}}}_{h}=\sqrt{g}\left(\begin{array}{c}0\\ -\rho {\partial }_{r}{\rm{\Psi }}+\rho r{\left({v}^{\theta }\right)}^{2}+\rho r{\sin }^{2}\theta {\left({v}^{\phi }\right)}^{2}+\tfrac{2P}{r}\\ -\rho {r}^{2}{\partial }_{\theta }{\rm{\Psi }}+\rho \sin \theta \cos \theta {\left({v}^{\phi }\right)}^{2}+\tfrac{P\cos \theta }{\sin \theta }\\ -\rho {\partial }_{\phi }{\rm{\Psi }}\\ -\rho {v}^{j}{\partial }_{j}{\rm{\Psi }}\\ 0\end{array}\right),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{\boldsymbol{W}}}_{i}=\sqrt{g}\left(\begin{array}{c}0\\ -{G}_{r}\\ -{G}_{\theta }\\ -{G}_{\phi }\\ {G}_{t}\\ -{\rm{\Gamma }}\end{array}\right),\end{eqnarray} \tag{ 6 }$$

where ρ, v^j, e, Y_e, P, Ψ, ${\delta }_{i}^{j}$, and $g\ (={r}^{2}\sin \theta )$ are the mass density, fluid velocity, internal energy density, electron fraction, pressure, Newtonian gravitational potential, Kronecker's delta, and volume factor in the spherical coordinates, respectively. Here the index j runs over the values 1, 2, and 3, which correspond to the three components of spatial coordinates. The numerical flux is determined by the HLL scheme with the piecewise-parabolic interpolation, and the time integration is done with the third-order Runge–Kutta method. Unlike in the papers by Nagakura et al. (2018) and Harada et al. (2019), we ignore the term of the coordinate acceleration  W_a given in Nagakura et al. (2017) and assume the monopole gravity: the gradient of Ψ is obtained as follows:

$$\begin{eqnarray}&&{\partial }_{j}{\rm{\Psi }}=\left(\displaystyle \frac{M(r)}{{r}^{2}},0,0\right),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&M(r)={\int }_{0}^{r}{dr}\int d{\rm{\Omega }}\ {r}^{2}\rho (r,\theta ,\phi ),\end{eqnarray} \tag{ 8 }$$

where $d{\rm{\Omega }}\ (=\sin \theta d\theta d\phi )$ is the differential solid angle in the spherical coordinates. Equation (6) presents the feedback from neutrino to hydrodynamics,

$$\begin{eqnarray}&&{G}^{\alpha }\equiv \displaystyle \sum _{s}{G}_{s}^{\alpha },\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{G}_{s}^{\alpha }\equiv \int {p}_{s}^{\alpha }{\left(\displaystyle \frac{\delta f}{\delta \lambda }\right)}_{\mathrm{col}(s)}{{dV}}_{p},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\rm{\Gamma }}\equiv {{\rm{\Gamma }}}_{{\nu }_{e}}-{{\rm{\Gamma }}}_{{\bar{\nu }}_{e}},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{s}\equiv \int {\left(\displaystyle \frac{\delta f}{\delta \lambda }\right)}_{\mathrm{col}(s)}{{dV}}_{p},\end{eqnarray} \tag{ 12 }$$

where ${p}_{s}^{\alpha }$, λ, and dV_p denote the four-momentum of neutrino of species s, affine parameter, and invariant volume in momentum space, respectively. The index α runs over the four values 0, 1, 2, and 3, corresponding to the time and space components. The collision term (δf/δλ)_col is related to (δf/δt)_col in the laboratory frame and (δf/δτ)_col in the fluid rest frame as follows:

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\delta f}{\delta \lambda }\right)}_{\mathrm{col}}=\epsilon {\left(\displaystyle \frac{\delta f}{\delta t}\right)}_{\mathrm{col}}={\epsilon }_{\mathrm{FR}}{\left(\displaystyle \frac{\delta f}{\delta \tau }\right)}_{\mathrm{col}},\end{eqnarray} \tag{ 13 }$$

where τ and _FR are the proper time of fluid element and the neutrino energy measured in the fluid rest frame, respectively.

Zoom In Zoom Out Reset image size

We consider the following neutrino reactions:

the absorption/emissions,

$$\begin{eqnarray}&&(\mathrm{ecp}):\ {e}^{-}+p\,\longleftrightarrow \,{\nu }_{e}+n,\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&(\mathrm{aecp}):\ {e}^{+}+n\,\longleftrightarrow \,{\bar{\nu }}_{e}+p,\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&(\mathrm{eca}):\ {e}^{-}+A\,\longleftrightarrow \,{\nu }_{e}+A^{\prime} ,\end{eqnarray} \tag{ 16 }$$

the scatterings,

$$\begin{eqnarray}&&(\mathrm{nsc}):\ {\nu }_{s}+N\,\longleftrightarrow \,{\nu }_{s}+N,\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&(\csc ):\ {\nu }_{s}+A\,\longleftrightarrow \,{\nu }_{s}+A,\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&(\mathrm{esc}):\ {\nu }_{s}+e\,\longleftrightarrow \,{\nu }_{s}+e,\end{eqnarray} \tag{ 19 }$$

and the pair processes,

$$\begin{eqnarray}&&(\mathrm{pap}):\ {e}^{-}+{e}^{+}\,\longleftrightarrow \,{\nu }_{s}+{\bar{\nu }}_{s},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&(\mathrm{nbr}):\ N+N\,\longleftrightarrow \,N+N+{\nu }_{s}+{\bar{\nu }}_{s},\end{eqnarray} \tag{ 21 }$$

where the subscript s again denotes the neutrino species. The detailed expressions of the collision terms are given in Sumiyoshi & Yamada (2012). Most of the reaction rates are taken from Bruenn (1985). The weak interactions of electron neutrino (ecl) with light nuclei,

$$\begin{eqnarray}&&(\mathrm{elpp}):\ {\nu }_{e}+{}^{2}{\rm{H}}\,\longleftrightarrow \,{e}^{-}+p+p,\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&(\mathrm{el}2{\rm{h}}):\ {\nu }_{e}+n+n\,\longleftrightarrow \,{e}^{-}+{}^{2}{\rm{H}},\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&(\mathrm{el}3{\rm{h}}):\ {\nu }_{e}+{}^{3}{\rm{H}}\,\longleftrightarrow \,{e}^{-}+{}^{3}\mathrm{He},\end{eqnarray} \tag{ 24 }$$

and those of antielectron neutrino (aecl):

$$\begin{eqnarray}&&(\mathrm{ponn}):\ {\bar{\nu }}_{e}+{}^{2}{\rm{H}}\,\longleftrightarrow \,{e}^{+}+n+n,\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&(\mathrm{po}2{\rm{h}}):\ {\bar{\nu }}_{e}+p+p\,\longleftrightarrow \,{e}^{+}+{}^{2}{\rm{H}},\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&(\mathrm{po}3{\rm{h}}):\ {\bar{\nu }}_{e}+{}^{3}\mathrm{He}\,\longleftrightarrow \,{e}^{+}+{}^{3}{\rm{H}},\end{eqnarray} \tag{ 27 }$$

are also taken into account. We use the multicompositional EOS at subnuclear densities (Furusawa et al. 2017) based on the supranuclear density EOS calculated with the variational method (Togashi & Takano 2013) and the new tables for electron capture by heavy and light nuclei and positron capture by light nuclei (Nagakura et al. 2019c).

In this study, we employ the 11.2M_⊙ progenitor model in Woosley et al. (2002). The spatial domain of 0 ≤ r ≤ 5000 km is divided into 384 radial grid cells in 1D simulations. The number is reduced to 256, and the spatial domain covers only 0 ≤ r ≤ 200 km in multidimensional simulations. The entire solid angle in space is divided into 48 grid cells in θ for 2D and 48 × 96 grid cells in θ × ϕ for 3D simulations. The neutrino energy in the range of 0 ≤  ≤ 300 MeV is divided into 16 cells. The entire solid angle in momentum space is divided into six cells in θ_ν for 2D and 6 × 6 cells in θ_ν × ϕ_ν for 3D simulations. More detailed discussions on the numerical resolution are provided in Appendix A.

The procedure of our 3D simulations is as follows:

• 1.  
  the corresponding 1D simulation is done from the onset of the core-collapse to a shock stagnation;
• 2.  
  the time-dependent boundary data are extracted from the 1D results for the use in multidimensional simulations;
• 3.  
  the multidimensional simulation is started with an introduction of 1% random velocity perturbations on the spherically symmetric flow when a negative entropy gradient emerges for the first time;
• 4.  
  the values of conserved variables at the outer boundary are obtained at every time step by linearly interpolating the boundary data obtained in step 2.

Furthermore, two kinds of coarse graining are implemented to relax the Courant–Friedrichs–Lewy condition for numerical stability around the singular points in spherical coordinates: near the coordinate origin (0 < r < 8 km), the physical quantities are averaged over the solid angle; in the vicinity of the polar axis, the conserved variables are averaged over every eight, four, and two cells in the ϕ direction for the first, second, and third θ cells from the polar axis, respectively (Figure 2). Note that the moving-mesh technique is not activated in this study.

Zoom In Zoom Out Reset image size

We present some noteworthy features in matter flows and neutrino radiation fields in 1D, 2D, and 3D simulations. Figures 3 and 4 show, respectively, the isosurfaces of entropy with the fluid velocity vectors and those of neutrino number density with the average velocity vectors for ν_e, ${\bar{\nu }}_{e}$, and ν_x (first, second, and third rows) at t = 10 ms. The average velocity of neutrino ${v}_{\nu }^{i}$ is defined as

$$\begin{eqnarray}&&{v}_{\nu }^{i}=\displaystyle \frac{{{ \mathcal F }}^{i}}{{ \mathcal N }}.\end{eqnarray} \tag{ 28 }$$

The neutrino number density ${ \mathcal N }$ and number flux ${{ \mathcal F }}^{i}$ are given as

$$\begin{eqnarray}&&{ \mathcal N }=\int {\epsilon }^{2}d\epsilon \int d{{\rm{\Omega }}}_{\nu }\ f(\epsilon ,{{\rm{\Omega }}}_{\nu }),\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{{ \mathcal F }}^{i}=\int {\epsilon }^{2}d\epsilon \int d{{\rm{\Omega }}}_{\nu }\ f(\epsilon ,{{\rm{\Omega }}}_{\nu })\ {n}^{i},\end{eqnarray} \tag{ 30 }$$

where nⁱ and $d{{\rm{\Omega }}}_{\nu }(=\sin {\theta }_{\nu }d{\theta }_{\nu }d{\phi }_{\nu })$ are the unit vector in space and the differential solid angle in the spherical coordinates, respectively. Here the integrations in Equations (29) and (30) are done in the laboratory frame. The isosurfaces are cut away between ϕ = 225° and 360°. The vectors are superimposed only on the meridian plane in the right half part. For better visibility, the vectors of the average velocity for ${\bar{\nu }}_{e}$ are not shown in r < 20 km. The shock wave is located at r = 70 km, corresponding to the surfaces of orange spheres. Unlike the results in 1D (Figure 3(a)), the prompt convection grows inside the shock wave in the multidimensional computations. Three vortex rings grow in 2D (Figure 3(b)), whereas multiple round convective vorticies are formed in 3D (Figure 3(c)). These convective flows affect the transport of ν_e, ${\bar{\nu }}_{e}$, and ν_x. Multidimensional structures are also developed in the number densities of neutrinos in both 2D (Figures 4(b), (e), and (h)) and 3D (Figures 4(c), (f), and (i)). In the central convective region, neutrinos move in various directions with matter. In the outer region, they propagate outward in both 2D (Figures 4(b), (e), and (h)) and 3D (Figures 4(c), (f), and (i)) anisotropically.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5 shows the isosurfaces of the density, the electron fraction, the entropy with the fluid velocity vectors, and the number density with the average velocity vectors for ν_e, ${\bar{\nu }}_{e}$, and ν_x on the equatorial plane for 3D. The isosurfaces are cut away above the equatorial plane. The vectors are superimposed on the cut planes. The integrations in Equations (29) and (30) are done in both the laboratory (LB) and fluid rest (FR) frames. The random perturbation added initially by hand grows exponentially in the region of the negative entropy gradient, and complex structures in density, electron fraction, entropy, and velocity develop as a result on the equatorial plane in the 3D simulation (Figures 5(a), (d), and (g)). The spatial distributions of number density and average velocity for ν_e, ${\bar{\nu }}_{e}$, and ν_x are also affected by them (Figures 5(b), (e), and (h)). Although the radiation fields evolve differently among three species of neutrinos owing to their interactions with matter in different ways, they have a common feature. The average neutrino velocity measured in the laboratory frame agrees roughly with the matter velocity in the central convective region (Figures 5(b), (e), (g), and (h)), whereas the average neutrino velocity in the fluid rest frame is quite small in the same region (Figures 5(c), (f), and (i)). This means that neutrinos are tightly coupled to matter in this optically thick region. The two-energy grid approach enables us to capture this effect to the full order of v/c (Nagakura et al. 2014). On the other hand, neutrinos begin to propagate outward freely in the outer region. Our code can capture the transition between the optically thick and thin limits.

Zoom In Zoom Out Reset image size

Figure 6 shows the isosurfaces of neutrino distribution function f at r = 30 km, which is in the diffusion regime. The directions of the three orthogonal axes P_X, P_Y, and P_Z in momentum space are chosen to agree with those of the θ-, ϕ-, and r-axes in space, respectively. The angle from the P_Z axis is defined as θ_ν, and the angle from the P_X axis on the P_X–P_Y plane is denoted by ϕ_ν (see Figure 1). The distance from the origin corresponds to the neutrino energy , and the length of the axes is 20 MeV in these plots. The neutrino distribution function is axisymmetric around the P_Z axis in 1D (Figures 6(a) and (d)), symmetric with respect to the P_Z–P_X plane in 2D (Figures 6(b) and (e)), whereas it does not have any symmetry in 3D (Figures 6(c) and (f)).

Zoom In Zoom Out Reset image size

The Eddington tensor is defined as the pressure tensor divided by the energy density. In the Boltzmann neutrino radiation transport simulations, we can directly calculate it from the distribution function (see Appendix B). The moment method is one of the common approximations for neutrino transport. In the M1 closure approximation, in particular, the evolution equations of the energy density E and flux Fⁱ are solved with an artificial closure relation for the pressure tensor P^ij (e.g., Thorne 1981; Shibata et al. 2011). In the following paragraph, we compare the Eddington tensor k^ij in Equation (B7) obtained with the Boltzmann simulation and the Eddington tensor ${k}_{M1}^{{ij}}$ estimated with the M1 approximation as in Equation (B13). Note that both k^ij and ${k}_{M1}^{{ij}}$ are measured in the laboratory frame throughout this paper.

Figure 7 shows the typical radial distributions of diagonal (top panels) and off-diagonal (bottom panels) components of the Eddington tensors k^ij and ${k}_{M1}^{{ij}}$ at ${\epsilon }_{\mathrm{FR}}=\langle {\epsilon }_{\mathrm{FR}}\rangle$ for ν_e at t = 10 ms in 1D (left panels), 2D (middle panels), and 3D (right panels) simulations. Here the average energy is defined as

$$\begin{eqnarray}&&\langle \epsilon \rangle =\displaystyle \frac{{ \mathcal E }}{{ \mathcal N }},\end{eqnarray} \tag{ 31 }$$

where ${ \mathcal N }$ is the number density in Equation (29) and ${ \mathcal E }$ is the energy density given as

$$\begin{eqnarray}&&{ \mathcal E }=\int {\epsilon }^{2}\,d\epsilon \int d{{\rm{\Omega }}}_{\nu }\ \epsilon f(\epsilon ,{{\rm{\Omega }}}_{\nu }).\end{eqnarray} \tag{ 32 }$$

The integration above is done in the fluid rest frame to calculate _FR. The Eddington tensors obtained in the Boltzmann simulations and in the M1 closure approximation are drawn with the solid and dashed lines, respectively. At t = 10 ms, the radial distributions of diagonal components of k^ij at ${\epsilon }_{\mathrm{FR}}=\langle {\epsilon }_{\mathrm{FR}}\rangle$ are almost the same among the 1D, 2D, and 3D cases (Figures 7(a), (b), and (c), respectively). They are almost 1/3 in the optically thick region extending up to 60 km. They start to deviate from 1/3 around the shock wave, with k^rr decreasing initially around 60 km and then increasing monotonically with radius. The other diagonal components have the opposite trend. On the other hand, the off-diagonal components of k^ij are quite different among the 1D, 2D, and 3D cases (Figure 7(d), (e), and (f), respectively). All components in 1D and k^rϕ and k^θϕ in 2D are zero identically owing to the symmetry of the neutrino distribution function in momentum space (Figures 6(d) and (e)). Only k^rθ in 2D and all off-diagonal components in 3D have nonvanishing values.

Zoom In Zoom Out Reset image size

Now we compare the Boltzmann Eddington tensor k^ij with the M1 counterpart ${k}_{M1}^{{ij}}$. The diagonal components kⁱⁱ differ from ${k}_{M1}^{{ii}}$ around the shock wave at t = 10 ms (Figures 7(a), (b), and (c)). In fact, the latter change monotonically with radius. It is also found that the absolute value of ${k}_{M1}^{\theta \phi }$ tends to be smaller than those of other off-diagonal components ${k}_{M1}^{r\theta }$ and ${k}_{M1}^{r\phi }$ in both the transition and optically thin regions (Figure 7(f)). This is understandable if one recalls that the radial flux dominates over other components in these regions. We find a large difference between k^θϕ and ${k}_{M1}^{\theta \phi }$ in the same regions (Figure 7(f)). This suggests that the interpolation between the optically thick and thin limits employed in the M1 prescription is not very good for this component in the 3D simulation.

In this section, we apply a new analysis in the comparison of the Eddington tensors obtained from the Boltzmann simulations and the M1 closure approximation. The Eddington tensor k^ij in Equation (B7) can be written as

$$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{k}} & = & {k}^{{rr}}({{\boldsymbol{e}}}_{r}\otimes {{\boldsymbol{e}}}_{r})+{k}^{r\theta }({{\boldsymbol{e}}}_{r}\otimes {{\boldsymbol{e}}}_{\theta })+{k}^{r\phi }({{\boldsymbol{e}}}_{r}\otimes {{\boldsymbol{e}}}_{\phi })\\ & & +\ {k}^{\theta r}({{\boldsymbol{e}}}_{\theta }\otimes {{\boldsymbol{e}}}_{r})+{k}^{\theta \theta }({{\boldsymbol{e}}}_{\theta }\otimes {{\boldsymbol{e}}}_{\theta })+{k}^{\theta \phi }({{\boldsymbol{e}}}_{\theta }\otimes {{\boldsymbol{e}}}_{\phi })\\ & & +\ {k}^{\phi r}({{\boldsymbol{e}}}_{\phi }\otimes {{\boldsymbol{e}}}_{r})+{k}^{\phi \theta }({{\boldsymbol{e}}}_{\phi }\otimes {{\boldsymbol{e}}}_{\theta })+{k}^{\phi \phi }({{\boldsymbol{e}}}_{\phi }\otimes {{\boldsymbol{e}}}_{\phi }),\end{array}\end{eqnarray} \tag{ 33 }$$

where  e_r,  e_θ, and  e_ϕ denote the basis vectors of the spherical coordinate system. Since k^ij is a real symmetric tensor, it has real eigenvalues and eigenvectors, which are easily obtained with the Jacobi method described in Appendix C. The diagonalized tensor  D of  k in Equation (C6) is expressed as

$$\begin{eqnarray}&&{\boldsymbol{D}}={\lambda }^{1}({{\boldsymbol{e}}}_{1}\otimes {{\boldsymbol{e}}}_{1})+{\lambda }^{2}({{\boldsymbol{e}}}_{2}\otimes {{\boldsymbol{e}}}_{2})+{\lambda }^{3}({{\boldsymbol{e}}}_{3}\otimes {{\boldsymbol{e}}}_{3}),\end{eqnarray} \tag{ 34 }$$

where λ^j and  e_j denote the jth eigenvalue and the corresponding normalized eigenvector of  k, respectively. Then, the rotation matrix  V in Equation (C7) is defined as

$$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{V}} & = & {V}^{r1}({{\boldsymbol{e}}}_{r}\otimes {{\boldsymbol{e}}}_{1})+{V}^{r2}({{\boldsymbol{e}}}_{r}\otimes {{\boldsymbol{e}}}_{2})+{V}^{r3}({{\boldsymbol{e}}}_{r}\otimes {{\boldsymbol{e}}}_{3})\\ & & +\ {V}^{\theta 1}({{\boldsymbol{e}}}_{\theta }\otimes {{\boldsymbol{e}}}_{1})+{V}^{\theta 2}({{\boldsymbol{e}}}_{\theta }\otimes {{\boldsymbol{e}}}_{2})+{V}^{\theta 3}({{\boldsymbol{e}}}_{\theta }\otimes {{\boldsymbol{e}}}_{3})\\ & & +\ {V}^{\phi 1}({{\boldsymbol{e}}}_{\phi }\otimes {{\boldsymbol{e}}}_{1})+{V}^{\phi 2}({{\boldsymbol{e}}}_{\phi }\otimes {{\boldsymbol{e}}}_{2})+{V}^{\phi 3}({{\boldsymbol{e}}}_{\phi }\otimes {{\boldsymbol{e}}}_{3}),\end{array}\end{eqnarray} \tag{ 35 }$$

where (V^rj, V^θj, V^ϕj) are the r-, θ-, and ϕ-components of the jth eigenvector. Using these eigenvalues and eigenvectors, we visualize the Eddington tensor as an ellipsoid in momentum space,

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\boldsymbol{p}}\cdot {{\boldsymbol{e}}}_{1}}{{\lambda }^{1}}\right)}^{2}+{\left(\displaystyle \frac{{\boldsymbol{p}}\cdot {{\boldsymbol{e}}}_{2}}{{\lambda }^{2}}\right)}^{2}+{\left(\displaystyle \frac{{\boldsymbol{p}}\cdot {{\boldsymbol{e}}}_{3}}{{\lambda }^{3}}\right)}^{2}=1,\end{eqnarray} \tag{ 36 }$$

where  p is the momentum vector. Note that the P_X, P_Y, and P_Z axes in momentum space are parallel to the basis vectors  e_θ,  e_ϕ, and  e_r, respectively (Figure 1). The surface of the triaxial ellipsoid is also represented parametrically as

$$\begin{eqnarray}&&\left(\begin{array}{c}{p}^{x}\\ {p}^{y}\\ {p}^{z}\end{array}\right)=\left(\begin{array}{ccc}{V}^{\theta 1} & {V}^{\theta 2} & {V}^{\theta 3}\\ {V}^{\phi 1} & {V}^{\phi 2} & {V}^{\phi 3}\\ {V}^{r1} & {V}^{r2} & {V}^{r3}\end{array}\right)\left(\begin{array}{l}{\lambda }^{1}\,\cos \,u\cos \,v\\ {\lambda }^{2}\,\cos \,u\sin \,v\\ {\lambda }^{3}\,\sin \,u\end{array}\right),\end{eqnarray} \tag{ 37 }$$

where (p^x, p^y, p^z) are the coordinates of a point on the ellipsoidal surface and u and v are the parameters. Of the three eigenvalues, λ³ is chosen to be the largest, representing the longest axis of the ellipsoid.

Figure 8 shows three representative configurations of the ellipsoid. In the optically thick limit, the ellipsoid is reduced to a sphere for the vanishing matter velocity (Figure 8(a)). Recall that we are working in the laboratory frame. The Eddington tensor has a single threefold degenerate eigenvalue and three mutually orthogonal eigenvectors, the orientation of which is arbitrary. In this case, the Jacobi method returns the eigenvectors as

$$\begin{eqnarray}\left(\begin{array}{ccc}{V}^{\theta 1} & {V}^{\theta 2} & {V}^{\theta 3}\\ {V}^{\phi 1} & {V}^{\phi 2} & {V}^{\phi 3}\\ {V}^{r1} & {V}^{r2} & {V}^{r3}\end{array}\right)=\left(\begin{array}{ccc}1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\end{array}\right).\end{eqnarray} \tag{ 38 }$$

Note, however, that the completely isotropic distribution is never realized in the supernova core owing to matter inhomogeneity, which would induce diffusion of neutrinos even in the optically thick region. In the M1 closure approximation applied in the fluid rest frame, the situation is the same, with the dominant isotropic term of γ^ij/3 in Equation (B12) being modified by the minor term of ${F}^{i}{F}^{j}/| F{| }^{2}$ in Equation (B11). It is noted that in this approximation F gives the eigenvector corresponding to the largest eigenvalue in almost the entire region for r ≳ 10 km at t = 10 ms postbounce in this study (see Appendix D). In the optically thin limit, on the other hand, the ellipsoid is reduced to a line (Figure 8(c)). There is only one nonvanishing eigenvalue. The corresponding eigenvector gives the direction of the neutrino energy flux also in this case. It is obvious that there is a twofold degenerate vanishing eigenvalue in addition to the unique nonvanishing eigenvalue. Finally, in the intermediate regime between the optically thick and thin limits, the ellipsoid is triaxial in general (Figure 8(b)). The longest principal axis is denoted by L hereafter. The directions of the neutrino energy flux F and the matter velocity V are also represented by the arrows with white and black heads, respectively, in the figure. The angle between F and L is designated as θ_FL, and its cosine ${\mu }_{{FL}}=\cos {\theta }_{{FL}}$ is used for the later analysis.

Zoom In Zoom Out Reset image size

Figure 9 shows the distribution functions f of ν_e in momentum space (top panels) and the corresponding ellipsoids of the Eddington tensors obtained in the Boltzmann simulation (middle panels) and those evaluated in the M1 approximation (bottom panels) at three radii from 10 to 40 km for the neutrino energy ${\epsilon }_{\mathrm{FR}}=\langle {\epsilon }_{\mathrm{FR}}\rangle$, i.e., the average energy. At r = 10 km, the neutrino distribution in momentum space is almost isotropic at all energies in the fluid rest frame owing to the tight coupling with matter, and as a result, F is nearly parallel to V in the laboratory frame (Figure 9(d)). The Eddington tensor that is evaluated in the laboratory frame is also slightly elongated in the direction of F. With the increasing radius, high-energy neutrinos are depleted (Figure 9(b)) and, more importantly, anisotropy becomes pronounced particularly for low-energy neutrinos (Figure 9(c)); it is observed at the same time (Figures 9(e) and (f)) that F and V start to be misaligned with each other in the convective region (r ∼ 40 km) as the interactions of neutrinos with matter get weaker and neutrinos diffuse out according to the local gradient of the neutrino number density. It is noted that the ellipsoids for the Boltzmann simulation (Figures 9(d)–(f)) agree well with those for the M1 approximation (Figures 9(g)–(i)), with the longest axis of the ellipsoid L being almost aligned with the direction of the energy flux F.

Zoom In Zoom Out Reset image size

Figure 10 shows the same quantities as Figure 9 does except for the spatial positions considered: 60–80 km. At these larger distances from the center of the core, the neutrinos are more concentrated toward the +P_Z axis in momentum space (Figures 10(a)–(c)). The energy flux F and the matter velocity V tend to the +P_Z and −P_Z axes, respectively, with increasing radius (Figures 10(d)–(f)). The ellipsoid changes its shape from a near sphere (Figure 10(d)) to a horizontally elongated ellipsoid (Figure 10(e)), and finally to a vertically extended ellipsoid (Figure 10(f)), where the terms "horizontal" and "vertical" mean "perpendicular" and "parallel" to the direction of F, respectively. The horizontally elongated configuration observed above is consistent with the fact that k^rr (k^θθ, k^ϕϕ) is lower (higher) than 1/3 there (Figures 7(a)–(c)), which can also be seen in the results of the neutrino radiation transport simulations solving the Boltzmann equation with the S_N method (Smit et al. 2000) and with the Monte Carlo approach (Janka et al. 1992; Murchikova et al. 2017). It occurs not only in 3D but also in 1D and 2D simulations. In the bottom panels of Figure 10, we exhibit the ellipsoids given by the M1 closure approximation. At 60 km, the directions of the longest principal axis L are different between the Boltzmann simulation and the M1 approximation, although the ellipsoids are almost spherical in both cases (Figures 10(d) and (g)). At r = 70 km, however, the difference is more substantial. In fact, the ellipsoid in the M1 approximation is extended vertically, whereas it is wider horizontally in the Boltzmann simulation (Figures 10(e) and (h)). The difference remains even at r = 80 km, where the M1 approximation gives a more elongated ellipsoid, although L in the Boltzmann simulation is now aligned with F (Figures 10(f) and (i)): ${k}_{M1}^{{rr}}$ > k^rr, ${k}_{M1}^{\theta \theta }\lt {k}^{\theta \theta }$, and ${k}_{M1}^{\phi \phi }\lt {k}^{\phi \phi }$ at r > 70 km (Figures 7(a)–(c)), which is also common to the 1D, 2D, and 3D simulations. Note that to what extent these are the case is rather sensitive to the numerical resolution (see Appendix A). We stress again, however, that L is always almost aligned with F in the M1 approximation (Figures 10(g), (h), and (i)), while it is not necessarily the case in reality (Figures 10(d) and (e)).

Zoom In Zoom Out Reset image size

Figure 11 shows the radial profile of $\langle | {\mu }_{{FL}}| {\rangle }_{{\rm{\Omega }}}$ for ν_e (left panels), ${\bar{\nu }}_{e}$ (middle panels), and ν_x (right panels), where $\langle | {\mu }_{{FL}}| {\rangle }_{{\rm{\Omega }}}$ denotes the absolute value of μ_FL averaged over the whole solid angle Ω at each radius. It is unity (zero) if L is parallel (perpendicular) to F. The vertical dashed line in the figure indicates the radial position of the shock wave. The prompt convection prevails from r ∼ 10 to 60 km while a laminar shocked flow exists from 60 to 70 km, and there is a supersonic accretion from 70 to 200 km. The results for ν_e in the Boltzmann simulation show that in the convective region $\langle | {\mu }_{{FL}}| {\rangle }_{{\rm{\Omega }}}$ gets small, ∼0.5, for low-energy neutrinos (reddish lines), while for the average-energy neutrinos (yellow and green lines) it is deviated from unity substantially in the laminar region; in the supersonic accretion flow high-energy neutrinos (bluish lines) are affected. The results for ${\bar{\nu }}_{e}$ and ν_x are more or less similar to those for ν_e, although the dip in the laminar region is less remarkable for the average-energy neutrinos (Figures 11(b) and (c)). The results for the M1 closure approximation are quite different, on the other hand. The longest principal axis of the ellipsoid obtained from ${k}_{M1}^{{ij}}$ is roughly parallel to the neutrino energy flux over the entire region, except just above the shock wave (Figures 11(d)–(f)). This is due to the fact that the terms of FⁱF^j/F² in Equation (B11) have a dominant role to determine the orientation of the ellipsoid in the M1 closure approximation even in the laboratory frame. As a result, the M1 prescription Equation (B9) cannot reproduce the horizontally wide ellipsoids, which are realized in different regions depending on the neutrino energy.

Zoom In Zoom Out Reset image size

Figure 12 shows the radial profile of Λ for ν_e (left panels), ${\bar{\nu }}_{e}$ (middle panels), and ν_x (right panels), where Λ is the mean free path normalized by the radius. The shock position is indicated by the vertical thick dashed line, whereas the horizontal thin dashed line corresponds to Λ = 1. In the bottom row of Figure 12, we show the radial profiles of the contributions to Λ of various interactions for the average neutrino energy (Figures 12(d)–(f)). In both of the convective and laminar shocked regions, the electron capture (ecp) gives the smallest Λ for ν_e, while the scattering with free nucleons (nsc) is dominant for ${\bar{\nu }}_{e}$ and ν_x. In the supersonic accretion flow, on the other hand, the coherent scattering with heavy nuclei (csc) determines Λ for all neutrino species. The region with Λ ≲ 1 is the transition region, and that is exactly where the approximation for neutrino transport should be validated by the Boltzmann simulation. For low-energy neutrinos (reddish lines), Λ ≲ 1 in the prompt convection region from 10 to 60 km. For the neutrinos around the average energy (yellow and green lines) the region with Λ ≲ 1 extends itself over both the prompt convection region and the laminar shocked region. The value of Λ for high-energy neutrinos (bluish lines) is below unity in the entire region up to 200 km, including the supersonic region. It is hence clear that the misalignment of F and L occurs in their own transition regions for different energies of neutrinos.

Zoom In Zoom Out Reset image size

One finds from Figure 10(b) for ${\epsilon }_{\mathrm{FR}}=\langle {\epsilon }_{\mathrm{FR}}\rangle \sim 12.6\,\mathrm{MeV}$ just behind the shock wave that the neutrino angular distribution becomes hemispheric, i.e., almost isotropic only in one hemisphere, and this is associated with the horizontally elongated ellipsoids in these regions discussed earlier. The hemispheric distribution is understood intuitively as follows: outgoing neutrinos with θ_ν < π/2 are mostly generated in the deeper and hence optically thicker region, while those going inward with θ_ν > π/2 are emitted or back-scattered in the outer optically thinner region; as the reaction rates are certainly larger for the former, we obtain the angular distribution with the latter neutrinos highly depleted. As the mean free path increases, the angular distribution of f in the neutrino-rich hemisphere becomes forward-peaked while f in the opposite hemisphere gets more diminished. The hemispheric distribution emerges typically at Λ ∼ 0.1 at t = 10 ms. It is noted that the maximum entropy distributions for Fermi–Dirac radiation (e.g., Janka et al. 1992; Cernohorsky & Bludman 1994), making the normalized radiation pressure smaller than 1/3 depending on the occupation density and flux factor, can give a similar distribution to the hemispheric one, where the boundary between neutrino-rich and neutrino-poor is located at μ_ν < 0. The relation of our results to the maximum entropy distributions for Fermi–Dirac radiation will be discussed in an upcoming paper. It should be stressed that in multidimensional cases these hemispheres do not exactly correspond to the outgoing or ingoing directions. At the early postbounce period we consider here, this is particularly clear for low-energy neutrinos. In Figure 6 we find the typical hemispheric distribution as the orange surface for  ≲ 10 MeV at r = 30 km. Its shape is axisymmetric with respect to the local radial direction in 1D (Figure 6(a)). In 2D and 3D (Figures 6(b) and (c)), it is distorted by the prompt convection, and, as mentioned, the two hemispheres are inclined from the local radial direction, although it remains plane-symmetric with respect to the P_Z–P_X plane in 2D. It is repeated that these hemispheric distributions are the main origin of the horizontally elongated ellipsoid, for which L and F are orthogonal to each other, and the M1 approximation failed to reproduce the Eddington tensor (Figure 10(e)).

In the multidimensional cases, the region in which L is perpendicular to F appears in space in a complex way. Figure 13 shows the color maps of unaveraged $| {\mu }_{{FL}}|$ on the equatorial plane for the Boltzmann simulation (first row) and the M1 closure approximation (second row) as well as Λ (third row) and $| F| /E$ (fourth row) for ν_e at t = 10 ms. The left, middle, and right panels are the results for _FR = 3.3, 6.4, and 12.6 MeV, respectively. In the upper two rows for $| {\mu }_{{FL}}|$, the thick colored region is the place where L and F are orthogonal to each other and the ellipsoid becomes horizontally wide. As was observed in Figure 11, the radius of the thick colored region gets larger with increasing energy in the Boltzmann simulation (Figures 13(a)–(c)), while there is no such region for the M1 closure approximation (Figures 13(d)–(f)). At low neutrino energies, _FR = 3.3 and 6.4 MeV, these regions are distributed in a patchy fashion in the prompt convection zone (Figures 13(a) and (b)) and the angle average of $| {\mu }_{{FL}}|$ becomes $\langle | {\mu }_{{FL}}| {\rangle }_{{\rm{\Omega }}}\gtrsim 0.3$ there (Figure 11(a)). For the average neutrino energy _FR ∼ 12.6 MeV (Figure 13(c)), on the other hand, it is extended uniformly in the laminar shocked flow and, as a result, $\langle | {\mu }_{{FL}}| {\rangle }_{{\rm{\Omega }}}$ is close to zero even just behind the shock wave (Figure 11(a)). However, there are also pockets of regions in the prompt convection zone, in which L and F are not aligned (Figure 13(c)). They appear correlated not with the color map of Λ (Figure 13(i)) but with $| F| /E$ (Figure 13(l)).

Zoom In Zoom Out Reset image size

Figures 14 and 15 are the Mollweide projection maps of $| {\mu }_{{FL}}|$ (left panels) and $| F| /E$ (right panels) at various radii for _FR = 12.6 and 3.3 MeV, respectively, at t = 10 ms. In Figure 14 for _FR = 12.6 MeV, we can confirm the correlation between the $| {\mu }_{{FL}}|$ and $| F| /E$, which we found in Figures 13(c) and (l) above, at r = 20–50 km, where Λ is ≳0.01. It is also observed that the regions of small $| {\mu }_{{FL}}|$ ellipsoids roughly agree with those of local minima of $| F| /E$. It turns out that in these regions, where Λ ≳ 0.01 for _FR = 12.6 MeV, the ellipsoids are almost spherical and slightly elongated in the direction perpendicular to F (see Figures 9(d)–(f)). For _FR = 3. 3 MeV, the pattern in $| {\mu }_{{FL}}|$ also looks very similar to that in $| F| /E$ at r = 40 km (Figures 15(c), (h)), where Λ is ≲1 (Figure 12(a)). Although the neutrino distribution starts to become forward peaked, it is still hemispheric and the ellipsoid is horizontally wide in the regions of $| {\mu }_{{FL}}| \sim 0$, which again roughly agree with the regions of local minima of $| F| /E$. This is not true, however, at smaller radii (Figures 15(a), (b), (f), (g)). Note that there are still many regions where $| {\mu }_{{FL}}|$ is much smaller than 1 and the ellipsoid is horizontally elongated substantially. They are located at Λ ∼ 0.1. It seems, however, that the convection makes situations much more complicated in 3D. In fact, the neutrino distribution function at _FR = 3.3 MeV is severely distorted by complex matter motions in these regions (see the orange-colored isosurfaces in Figures 6(c) and (f)), which makes the neutrino flux less correlated with the shape of the ellipsoid. If convective motions at Λ ∼ 0.1 are responsible indeed, it may not be easy to predict the place where the horizontally elongated ellipsoids are produced and make an appropriate modification to the M1 prescription. That will be a future work.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have done a radiation-hydrodynamical simulation of core-collapse supernovae for a 11.2M_⊙ progenitor model in 3D space with the full Boltzmann neutrino transport until 20 ms after bounce. The time-dependent six-dimensional Boltzmann equations for three species of neutrinos and the 3D hydrodynamic equations with the monopole Newtonian self-gravity have been solved with Furusawa and Togashi's EOS. What we have done in this paper is summarized as follows:

• 1.  
  We have investigated the neutrino distributions in the 3D space at this early postbounce phase. Multiple round vortices are generated in the prompt convection that sets in at ∼10 ms and produce, in turn, 3D structures of density, entropy, and electron fraction in space. We have observed that neutrinos move along with matter in the optically thick region, start to decouple from matter in the transition region between the optically thick and thin limits, and finally propagate outward in the optically thin region.
• 2.  
  We have confirmed that the neutrino angular distributions in momentum space have no symmetry for 3D, while they have axisymmetry with respect to the radial direction and reflection symmetry with respect to the P_Z–P_X plane for 1D and 2D, respectively; all the off-diagonal components of Eddington tensor are nonvanishing for 3D; and there are some differences between the Eddington tensors obtained from the Boltzmann simulation and from the M1 closure approximation.
• 3.  
  We have applied a new analysis based on the principal axis transformation for the Eddington tensor to better understand these differences. We visualize the Eddington tensor as the ellipsoid whose principal axes are parallel to its eigenvectors, having lengths proportional to the corresponding eigenvalues. The ellipsoid is reduced to a sphere in the optically thick limit and to a line parallel to the energy flux in the optically thin limit. In between it is a triaxial ellipsoid in general. We have found that the ellipsoid obtained directly from the Boltzmann simulation sometimes becomes horizontally wide, that is, elongated in the direction perpendicular to the energy flux, in the transition regime between the optically thick and thin limits. This is in sharp contrast to the M1 closure approximation, in which the longest principal axis of the ellipsoid is always almost parallel to the direction of the energy flux.
• 4.  
  The horizontally wide Eddington tensor emerges with high probabilities in the transition region between the optically thick and thin limits, where the mean free path divided by the radius is Λ ∼ 0.1. As a matter of fact, in the convective region, the horizontally wide ellipsoid tends to occur at places where Λ is a bit higher or lower than 0.1 and, in addition, the neutrino velocity $| F| /E$ takes locally minimum values. We have observed, however, that convective matter motions make the situation more complicated. This may make it difficult to improve the M1 approximation in these regions.

This paper is the very first step in our project of 3D Boltzmann radiation-hydrodynamics simulations of core-collapse supernovae. In fact, we have focused only on the prompt convection phase at t ∼ 10 ms postbounce and paid attention to the neutrino distributions in momentum space. In this very early phase, the entire postshock region is optically thick with Λ < 1 for neutrinos with ≳10 MeV. The convection region is hence located deep in the optically thick region for the average neutrinos. At much later times, the neutrino-driven convection or the SASI and even the lepton-driven convection occur in more optically thin regions. We hence should examine how the Eddington tensor is affected by these multidimensional matter motions: in particular, we are interested in whether the horizontally elongated ellipsoid occurs in a similar fashion and, more importantly, if it has an impact on the shock revival. If it does, it will be important to try seriously to improve the M1 closure approximation. Not to mention, our simulation is hardly perfect. The numerical resolutions in space and momentum space, the monopolar Newtonian gravity, and the 1D treatment of the innermost region at r < 8 km are certainly concerns, to mention a few. We are planning to perform longer simulations with better resolutions on the Fugaku computer, the next-generation flagship supercomputer of Japan whose capability is expected to be more than 10 times as large as K-computer's, and will address these issues there.

This research used high-performance computing resources of the K-computer/the supercomputer Fugaku provided by RIKEN, the FX10 provided by Tokyo University, the FX100 provided by the Nagoya University, the Grand Chariot provided by Hokkaido University, and the Oakforest-PACS provided by JCAHPC through the HPCI System Research Project (Project ID: hp130025, 140211, 150225, 150262, 160071, 160211, 170031, 170230, 170304, 180111, 180179, 180239, 190100, 190160, 200102, 200124), SR16000 and XC40 at YITP of Kyoto University, NEC SX Aurora Tsubasa at KEK, Research Center for Nuclear Physics (RCNP) at Osaka University, and the XC50 and the general common-use computer system provided by CfCA in the National Astronomical Observatory of Japan (NAOJ). Large-scale storage of numerical data is supported by JLDG constructed over SINET4 of NII. This work is supported in part by Grants-in-Aid for Scientific Research (26104006, 15K05093, 19K03837, 20H01905) and Grant-in-Aid for Scientific Research on Innovative areas "Gravitational wave physics and astronomy: Genesis" (17H06357, 17H06365) and "Unraveling the History of the Universe and Matter Evolution with Underground Physics" (19H05802) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. This work was also partly supported by research programs at K-computer of the RIKEN R-CCS, HPCI Strategic Program of Japanese MEXT, Priority Issue on Post-K-computer (Elucidation of the Fundamental Laws and Evolution of the Universe), Joint Institute for Computational Fundamental Sciences (JICFus) and MEXT as "Program for Promoting Researches on the Supercomputer Fugaku" (Toward a unified view of the universe: from large scale structures to planets). A.H. is supported in part by Japan Society for the Promotion of Science (JSPS) Grant-in-Aid for Research Activity Start-up (19K23435). S.F. is supported by JSPS KAKENHI (19K14723). S.Y. is supported by Institute for Advanced Theoretical and Experimental Physics, and Waseda University and the Waseda University Grant for Special Research Projects (project No. 2020-C273). K.S. acknowledges the Particle, Nuclear and Astro Physics Simulation Program (2019-002, 2020-004) at KEK.

In this appendix we report some of the results of resolution tests we have conducted that we think are of relevance for this paper. We refer readers also to our earlier publications for more tests (Richers et al. 2017; Nagakura et al. 2018; Harada et al. 2019).

The canonical number of numerical grid points deployed in this paper is N_r × N_θ × N_ϕ × N × N_θν × N_ϕν = 256 × 48 × 96 × 16 × 6 × 6, where N_r, N_θ, N_ϕ, N, ${N}_{{\theta }_{\nu }}$, and ${N}_{{\phi }_{\nu }}$ are the grid numbers of radial, polar, and azimuthal angles in space, neutrino energy, and polar and azimuthal angles in momentum space, respectively. Note that the grid spacing of the radial mesh is the same as in Nagakura et al. (2018), while the computational domain is confined to r ≤ 200 km; the number of numerical grid points in momentum space is reduced from ${N}_{\epsilon }\times {N}_{{\theta }_{\nu }}\,\times {N}_{{\phi }_{\nu }}=20\times 10\times 6$ in Nagakura et al. (2018) to 16 × 6 × 6 in this paper.

We first conduct the resolution test for the spatial angular mesh in 2D. Comparing the results among N_θ = 48, 64, and 96, we find that at t = 10 ms the prompt convection grows with three coherent vortex rings for N_θ = 48, whereas the number increases to four for N_θ = 64. At N_θ = 96 we observe three highly deformed and less coherent rings. We also observe that modes with smaller sizes grow more rapidly, reaching the nonlinear phase earlier, as the spatial resolution increases. The differences are not so large as to affect the conclusions in this paper, though. We hence set N_θ = 48 throughout the paper. The number of the azimuthal grid points N_ϕ is chosen to be twice N_θ.

Next, we perform 1D simulations to examine effects of the energy and angular resolutions in momentum space. Richers et al. (2017) conducted similar resolution tests, comparing the results from the Boltzmann simulations with those obtained with their Monte Carlo code. They employed three different grids: $({N}_{\epsilon },{N}_{{\theta }_{\nu }})$ = (20, 10), (40, 20), and (80, 40). Their results indicate that even $({N}_{\epsilon },{N}_{{\theta }_{\nu }})=(80,40)$ may not be large enough in the discrete ordinates code (Figure 8 in Richers et al. 2017). Since they changed both N and ${N}_{{\theta }_{\nu }}$ at the same time, however, it is not clear which resolution is more important for their results. Here we vary each of the two numbers individually. Figure 16 shows the energy resolution dependence of (a) k^rr and (b) ${k}_{M1}^{{rr}}$, in which N = 8, 16, and 32. The radial profiles of k^rr for N = 16 are in good agreement with those for N = 32, with deviations of the order of O(10⁻³). This is also the case for ${k}_{M1}^{{rr}}$. We hence conclude that N = 16 is good enough for our analysis in this paper. Note also that this is still larger than the number of energy grid points in most supernova simulations.

Zoom In Zoom Out Reset image size

On the other hand, Figure 17 shows the angular resolution dependence of (a) k^rr and (b) ${k}_{M1}^{{rr}}$, in which we deploy ${N}_{{\theta }_{\nu }}=6$, 12, and 48. In the optically thin region for ${\epsilon }_{\mathrm{FR}}=\langle {\epsilon }_{\mathrm{FR}}\rangle$ at r > 70 km, k^rr becomes larger with increasing ${N}_{{\theta }_{\nu }};$ it is obvious that ${N}_{{\theta }_{\nu }}=6$ is not large enough, and it is interesting that ${k}_{M1}^{{rr}}$ is smaller than k^rr for ${N}_{{\theta }_{\nu }}=12$ and 48 while it is larger for ${N}_{{\theta }_{\nu }}=6$. It should be noted, however, that the Eddington tensor evaluated with the M1 closure approximation is based on the same numerical data obtained by the Boltzmann simulation. Deeper inside at r = 30–60 km, the radial profiles of k^rr and ${k}_{M1}^{{rr}}$ for different values of ${N}_{{\theta }_{\nu }}$ both indicate again that ${N}_{{\theta }_{\nu }}=6$ is not large enough to get numerical convergence. However, the important feature that there exists a region at r ∼ 55–70 km, in which k^rr < 1/3 while ${k}_{M1}^{{rr}}\gt 1/3$, is common to the three cases with ${N}_{{\theta }_{\nu }}=6$, 12, and 48. We hence consider that ${N}_{{\theta }_{\nu }}=6$ is good enough for the principal axis analysis of the Eddington tensor in this paper.

Zoom In Zoom Out Reset image size

Now we turn to the 3D case and consider the angular resolution in momentum space. The higher-resolution runs with $({N}_{{\theta }_{\nu }},{N}_{{\phi }_{\nu }})=(6,12)$ and (12,6) are initiated from t = 9.5 ms postbounce, interpolating the result for the original resolution with $({N}_{{\theta }_{\nu }},{N}_{{\phi }_{\nu }})=(6,6)$. Figure 18 shows the radial profiles of diagonal (top panels) and off-diagonal (bottom panels) components of the Eddington tensor k^ij (left panels) and ${k}_{M1}^{{ij}}$ (right panels) at ${\epsilon }_{\mathrm{FR}}=\langle {\epsilon }_{\mathrm{FR}}\rangle$ for ν_e at t = 10 ms. We find for the diagonal components that the larger ${N}_{{\theta }_{\nu }}$ raises k^rr and lowers ${k}^{\theta \theta }$ and ${k}^{\phi \phi }$ at r > 70 km, which is essentially the same as what we saw above in the 1D tests, while the higher value of ${N}_{{\phi }_{\nu }}$ does not induce much change for them in the whole region (Figure 18(a)). As we mentioned in the previous paragraph for the 1D resolution tests, the region at r > 70 km is optically thin for ${\epsilon }_{\mathrm{FR}}=\langle {\epsilon }_{\mathrm{FR}}\rangle$ and the ellipsoids are elongated vertically. Those are the features better captured with higher resolutions. At smaller radii r < 70 km the value of ${N}_{{\theta }_{\nu }}$ less affects the diagonal components of the Eddington tensor just as in the 1D resolution tests. In fact, the important feature that the horizontally wide ellipsoids occur is unaffected by the resolution. This is also the case for ${k}_{M1}^{{ij}}$ in the sense that it produces almost always vertically long ones irrespective of the angular resolutions. The off-diagonal components, on the other hand, are also affected by the angular resolutions, particularly at r > 70 km for the same reason (Figure 18(c)). Note that they are much smaller than the diagonal components (see the difference in the scales of the vertical axes in the left panels). It is interesting that ${k}^{\theta \phi }$ is somewhat more sensitive to the resolution in ϕ_ν direction. We also point out that it is insensitive to ${N}_{{\theta }_{\nu }}$ and ${N}_{{\phi }_{\nu }}$ that ${k}_{M1}^{{rr}}$ (${k}_{M1}^{\theta \theta }$ and ${k}_{M1}^{\phi \phi }$) increases (decrease) monotonically just behind the shock wave (Figure 18(b)) and that ${k}_{M1}^{\theta \phi }$ is much smaller than ${k}_{M1}^{r\theta }$ and ${k}_{M1}^{r\phi }$ at r > 30 km (Figure 18(d)).

Zoom In Zoom Out Reset image size

Now we investigate how the resolution dependence observed in Figure 18 is reflected in the principal axis analysis for the Eddington tensor. Figure 19 shows the radial profiles of $\langle | {\mu }_{{FL}}| {\rangle }_{{\rm{\Omega }}}$ for k^ij and ${k}_{M1}^{{ij}}$ at various energies. It is recalled that $| {\mu }_{{FL}}| =0$ and 1 mean, respectively, that the longest principal axis of the Eddington tensor is perpendicular and parallel to the energy flux. We find that the results for the higher resolutions agree well with those for the original resolution below the shock wave at every energy (Figures 19(a), (b), and (c)). The results for the M1 closure approximation share the same feature of $\langle | {\mu }_{{FL}}| {\rangle }_{{\rm{\Omega }}}\sim 1$ among the different resolutions (Figures 19(d), (e), and (f)). Figures 20 and 21 present the color maps of $| {\mu }_{{FL}}|$ at _FR = 12.6 and 3.3 MeV, respectively, on the equatorial plane for ν_e at t = 10 ms in the Boltzmann simulations. The 1D results are also shown for reference. The thick colors indicate the regions where the longest principal axis of the ellipsoid is perpendicular to the energy flux. In the 1D case for _FR = 12. 6 MeV, there are two such regions (Figure 20(a)). In the inner region at Λ ∼ 0.01 the ellipsoids are almost spherical with aspect ratios of ≲1.004, while in the outer region at Λ ∼ 0.1 they have higher aspect ratios of ≲1.147. In the 1D simulation for _FR = 3.3 MeV, on the other hand, there exists one thick colored region at 0.1 ≲ Λ ≲ 1 (Figure 21(a)), where the ellipsoids have aspect ratios of ≲1.042. The appearance patterns of the thick colored regions in the 3D runs agree well among all the resolutions: $({N}_{{\theta }_{\nu }},{N}_{{\phi }_{\nu }})=(6,6),(12,6)$, and (6, 12) for both _FR = 12.6 (Figures 20(b), (c), and (d)) and 3.3 MeV (Figures 21(b), (c), and (d)). The outer circular region is what we found below the shock front commonly irrespective of the dimension of the simulations. The inner zone observed in the 1D case was shredded into smaller pieces by the convective motions of matter in 3D. The important thing, as we mentioned, is that the pattern is essentially the same among the different resolutions. We focus in this paper only on these features that are rather insensitive to the grid resolution. Note that the Eddington tensor corresponds to the second angular moments of the distribution function, that is, rather low order, and its ellipsoid in particular reflects not very fine but global features in the angular distribution in momentum space. It is hence not surprising that even the angular grid with $({N}_{{\theta }_{\nu }},{N}_{{\phi }_{\nu }})=(6,6)$ can capture those low-order features rather well except in the optically thin region. It is true that we have not yet seen the numerical convergence as observed in Figure 18, but we believe that the remaining errors do not change the conclusions in this paper. However, they are based on the single 3D simulation, and we certainly need to confirm them with more simulations for different postbounce times, other progenitor models, and different spatial and momentum resolutions in the future.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this appendix we give some definitions of the angular moments of the distribution function in phase space, as well as of the Eddington tensor, following the moment formalism described in Thorne (1981) and Shibata et al. (2011). An unprojected second moment is first defined in an arbitrary frame as

$$\begin{eqnarray}&&{M}^{\alpha \beta }({\epsilon }_{\mathrm{FR}})\equiv \int f({p}^{\alpha })\delta \left(\displaystyle \frac{{\epsilon }_{\mathrm{FR}}^{3}}{3}-\displaystyle \frac{{\left(-{u}_{\alpha }{p}^{\alpha }\right)}^{3}}{3}\right){p}^{\alpha }{p}^{\beta }{{dV}}_{p},\end{eqnarray} \tag{ B1 }$$

where f, p^α, u^α, and dV_p are the neutrino distribution function, four-momentum of neutrinos, four-velocity of medium, and invariant integration element, respectively. Note that _FR is a parameter to be specified, for example, as the neutrino energy measured in the fluid rest frame as adopted here. The Greek indices α, β, and γ run over 0, 1, 2, and 3, representing the time and space components. The second angular moment can be written as

$$\begin{eqnarray}\begin{array}{rcl}{M}^{\alpha \beta }({\epsilon }_{\mathrm{FR}}) & = & J({\epsilon }_{\mathrm{FR}}){u}^{\alpha }{u}^{\beta }+{H}^{\alpha }({\epsilon }_{\mathrm{FR}}){u}^{\beta }\\ & & +\ {H}^{\beta }({\epsilon }_{\mathrm{FR}}){u}^{\alpha }+{L}^{\alpha \beta }({\epsilon }_{\mathrm{FR}}),\end{array}\end{eqnarray} \tag{ B2 }$$

where the energy density J, energy flux H, and radiation pressure tensor L are the variables projected on to the fluid rest flame. In this study, in which general relativity is ignored, they are calculated as follows:

$$\begin{eqnarray}&&J({\epsilon }_{\mathrm{FR}})={J}_{\mathrm{FR}}({\epsilon }_{\mathrm{FR}}),\end{eqnarray} \tag{ B3 }$$

$$\begin{eqnarray}&&{H}^{\alpha }({\epsilon }_{\mathrm{FR}})={{\rm{\Lambda }}}_{\beta }^{\alpha }{H}_{\mathrm{FR}}^{\beta }({\epsilon }_{\mathrm{FR}}),\end{eqnarray} \tag{ B4 }$$

$$\begin{eqnarray}&&{L}^{\alpha \beta }({\epsilon }_{\mathrm{FR}})={{\rm{\Lambda }}}_{\gamma }^{\alpha }{{\rm{\Lambda }}}_{\delta }^{\beta }{L}_{\mathrm{FR}}^{\gamma \delta }({\epsilon }_{\mathrm{FR}}),\end{eqnarray} \tag{ B5 }$$

where ${{\rm{\Lambda }}}_{\beta }^{\alpha }$ is an appropriate Lorentz transformation; the energy density J_FR, energy flux H_FR, and radiation pressure tensor L_FR on the right-hand side are derived from the numerical integration of the distribution function obtained in the Boltzmann simulation. The same moment can be also expressed as

$$\begin{eqnarray}\begin{array}{rcl}{M}^{\alpha \beta }({\epsilon }_{\mathrm{FR}}) & = & E({\epsilon }_{\mathrm{FR}}){n}^{\alpha }{n}^{\beta }+{F}^{\alpha }({\epsilon }_{\mathrm{FR}}){n}^{\beta }\\ & & +\ {F}^{\beta }({\epsilon }_{\mathrm{FR}}){n}^{\alpha }+{P}^{\alpha \beta }({\epsilon }_{\mathrm{FR}}),\end{array}\end{eqnarray} \tag{ B6 }$$

where n^α is a unit timelike vector orthogonal to the hypersurface of constant coordinate time. The energy density E, energy flux F, and radiation pressure tensor P on the right-hand side of this equation are the variables projected onto the laboratory frame. The Eddington tensor k^ij is then defined as

$$\begin{eqnarray}&&{k}^{{ij}}({\epsilon }_{\mathrm{FR}})=\displaystyle \frac{{P}^{{ij}}({\epsilon }_{\mathrm{FR}})}{E({\epsilon }_{\mathrm{FR}})},\end{eqnarray} \tag{ B7 }$$

where the indices i and j run over 1, 2, and 3, corresponding to the space components.

In the truncated moment scheme, the time evolution equations of E and F are normally solved with the algebraic closure relation imposed by hand. In this paper, we investigate the M1 closure relation, the most favorite choice these days, for which the radiation pressure tensor is expressed as

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

where ζ is the variable Eddington factor given, for example (Levermore 1984), as

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

where the flux factor is defined, following Shibata et al. (2011), as

$$\begin{eqnarray}&&\bar{F}({\epsilon }_{\mathrm{FR}})=\sqrt{\displaystyle \frac{{h}_{\alpha \beta }{H}^{\alpha }({\epsilon }_{\mathrm{FR}}){H}^{\beta }({\epsilon }_{\mathrm{FR}})}{J{\left({\epsilon }_{\mathrm{FR}}\right)}^{2}}}.\end{eqnarray} \tag{ B10 }$$

The optically thin limit of P^ij is

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

whereas the optically thick limit is

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

where Vⁱ and γ^ij = g^ij + nⁱn^j stand for the 3D fluid velocity and the spatial metric on the hypersurface of constant coordinate time, respectively. The Eddington tensor in the M1 prescription is defined as

$$\begin{eqnarray}&&{k}_{M1}^{{ij}}({\epsilon }_{\mathrm{FR}})=\displaystyle \frac{{P}_{M1}^{{ij}}({\epsilon }_{\mathrm{FR}})}{E({\epsilon }_{\mathrm{FR}})}.\end{eqnarray} \tag{ B13 }$$

Jacobi's approach is one of the methods to obtain eigenvalues and eigenvectors of real symmetric matrices (e.g., Press et al. 1992). Let us consider diagonalizing a real symmetric matrix  A. If the maximum absolute value of the off-diagonal components of  A is given by the pq component, the matrix  A is transformed to

$$\begin{eqnarray}&&{\boldsymbol{A}}^{\prime} ={{\boldsymbol{P}}}_{{pq}}^{t}\cdot {\boldsymbol{A}}\cdot {{\boldsymbol{P}}}_{{pq}},\end{eqnarray} \tag{ C1 }$$

where  P_pq is the basic Jacobi rotation matrix,

$$\begin{eqnarray}{{\boldsymbol{P}}}_{{pq}}=\left(\begin{array}{ccccccc}1 & \cdots & 0 & \cdots & 0 & \cdots & 0\\ \vdots & \ddots & \vdots & & \vdots & & \vdots \\ 0 & \cdots & \cos \theta & \cdots & \sin \theta & \cdots & 0\\ \vdots & & \vdots & \ddots & \vdots & & \vdots \\ 0 & \cdots & -\sin \theta & \cdots & \cos \theta & \cdots & 0\\ \vdots & & \vdots & & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & \cdots & 0 & \cdots & 1\end{array}\right).\end{eqnarray} \tag{ C2 }$$

Note that only the components on the pth and qth rows and columns are transformed as follows:

$$\begin{eqnarray}{\boldsymbol{A}}^{\prime} =\left(\begin{array}{ccccccc} & \cdots & {a}_{1p}^{{\prime} } & \cdots & {a}_{1q}^{{\prime} } & \cdots & \\ \vdots & \ddots & \vdots & & \vdots & & \vdots \\ {a}_{p1}^{{\prime} } & \cdots & {a}_{{pp}}^{{\prime} } & \cdots & {a}_{{pq}}^{{\prime} } & \cdots & {a}_{{pn}}^{{\prime} }\\ \vdots & & \vdots & \ddots & \vdots & & \vdots \\ {a}_{q1}^{{\prime} } & \cdots & {a}_{{qp}}^{{\prime} } & \cdots & {a}_{{qq}}^{{\prime} } & \cdots & {a}_{{qn}}^{{\prime} }\\ \vdots & & \vdots & & \vdots & \ddots & \vdots \\ & \cdots & {a}_{{np}}^{{\prime} } & \cdots & {a}_{{nq}}^{{\prime} } & \cdots & \end{array}\right),\end{eqnarray} \tag{ C3 }$$

$$\begin{eqnarray}&&\begin{array}{l}{a}_{{rp}}^{{\prime} }={a}_{{rp}}\cos \theta -{a}_{{rq}}\sin \theta ,\\ \quad \ (r\ne p,r\ne q,r=1,2,\cdots ,n),\\ {a}_{{rq}}^{{\prime} }={a}_{{rq}}\cos \theta +{a}_{{rp}}\sin \theta ,\\ \quad \ (r\ne p,r\ne q,r=1,2,\cdots ,n),\\ {a}_{{pp}}^{{\prime} }={a}_{{pp}}{\cos }^{2}\theta +{a}_{{qq}}{\sin }^{2}\theta -2{a}_{{pq}}\sin \theta \cos \theta ,\\ {a}_{{qq}}^{{\prime} }={a}_{{qq}}{\cos }^{2}\theta +{a}_{{pp}}{\sin }^{2}\theta +2{a}_{{pq}}\sin \theta \cos \theta ,\\ {a}_{{pq}}^{{\prime} }={a}_{{pq}}({\cos }^{2}\theta -{\sin }^{2}\theta )+({a}_{{pp}}-{a}_{{qq}})\sin \theta \cos \theta .\end{array}\end{eqnarray} \tag{ C4 }$$

In the Jacobi method, each rotation is chosen so that the off-diagonal pq component should be zero, ${a}_{{pq}}^{{\prime} }=0$. Then, we have the following relations:

$$\begin{eqnarray}\begin{array}{rcl}\cos \theta & = & \sqrt{\displaystyle \frac{1+\alpha }{2}},\\ \sin \theta & = & \sqrt{\displaystyle \frac{1-\alpha }{2}}\mathrm{sgn}(\beta \gamma ),\\ \alpha & = & \displaystyle \frac{| \beta | }{\sqrt{{\beta }^{2}+{\gamma }^{2}}},\ \ \beta =\displaystyle \frac{1}{2}({a}_{{pp}}-{a}_{{qq}}),\ \ \gamma =-{a}_{{pq}}.\end{array}\end{eqnarray} \tag{ C5 }$$

Repeating this operation until all the off-diagonal components vanish (actually fall below 1.0 × 10⁻⁸), we eventually obtain the diagonal matrix  D that satisfies

$$\begin{eqnarray}&&{\boldsymbol{D}}={{\boldsymbol{V}}}^{t}\cdot {\boldsymbol{A}}\cdot {\boldsymbol{V}},\end{eqnarray} \tag{ C6 }$$

where  V is the rotation matrix given as

$$\begin{eqnarray}&&{\boldsymbol{V}}={{\boldsymbol{P}}}_{1}\cdot {{\boldsymbol{P}}}_{2}\cdot {{\boldsymbol{P}}}_{3}\cdots ,\end{eqnarray} \tag{ C7 }$$

with  P_i being one of the successive Jacobi rotation matrices. The diagonal components of  D finally give the eigenvalues of the original matrix  A, and the columns of  V are the corresponding eigenvectors.

In this paper, we analyze the Eddington tensor defined in the laboratory frame (Equation (B7)). We can do the same in the fluid rest frame, though. It is in fact easier in the optically thick regime. Figure 22 shows the radial profiles of $\langle | {\mu }_{{FL}}| {\rangle }_{{\rm{\Omega }}}$ in the fluid rest frame for ν_e at t = 10 ms postbounce, where $\langle | {\mu }_{{FL}}| {\rangle }_{{\rm{\Omega }}}$ is the absolute value of μ_FL, the cosine of the angle between F and L, averaged over the whole solid angle Ω. It should be mentioned first that if the angular distribution is completely isotropic, the Eddington tensor is proportional to the unit tensor and the longest principal axis is arbitrary but set to be radial in our code; then, μ_FL becomes the cosine of the angle between F and r. This never happens in the supernova core in reality, since neutrinos diffuse out of fluid elements to the adjacent ones according to the gradient of energy density even in the optically thick region. This is demonstrated in Figure 22(a), in which we show the cosine of the angle between F and r averaged over the solid angle to be compared with the $\langle | {\mu }_{\mathrm{FL}}| {\rangle }_{{\rm{\Omega }}}$ presented in Figure 22(c) in the same figure. In the fluid rest frame, the pressure tensor in the M1 closure approximation has just two terms of γ^ij and FⁱF^j, in which case L should be simply parallel to F. This is indicated indeed at r > 14 km in Figure 22(b), in which we show the result for the M1 approximation. The small deviation from unity at r < 14 km comes from ζ − 1/3 ≲ O(10⁻⁸), where ζ is the variable Eddington factor in Equation (B9). Then, the first term of Equation (B8) divided by E(_FR), which determines the direction of L in the fluid rest frame, becomes lower than O(10⁻⁸) at r < 14 km. This means that the analysis can be precisely available for r > 14 km at t = 10 ms postbounce since the accuracy of the Jacobi method is O(10⁻⁸). It is further pointed out that the radial profiles of $\langle | {\mu }_{{FL}}| {\rangle }_{{\rm{\Omega }}}$ in the fluid rest frame for the Boltzmann simulation (Figure 22(c)) are qualitatively different from those of $\langle | {\mu }_{{FL}}| {\rangle }_{{\rm{\Omega }}}$ in Figure 22(a) in the transition regions. This is simply because L is not radial there. The radius at which $\langle | {\mu }_{{FL}}| {\rangle }_{{\rm{\Omega }}}$ gets minimum in the fluid rest frame (Figure 22(c)) gets larger with the neutrino energy, the feature also observed in the laboratory frame (Figure 11(a)). The minimum values in the fluid rest frame, on the other hand, tend to be lower than those in the laboratory frame. Regardless of the reference frame, the message here is that there are regions in which the longest principal axis of the ellipsoid derived from the Eddington tensor is perpendicular to the energy flux in the Boltzmann simulation.

Zoom In Zoom Out Reset image size