Properties of Neutrino Transfer in a Deformed Remnant of a Neutron Star Merger

Neutron star mergers have attracted a surge of interest as a target of multi-messenger observations to reveal the fate of compact objects and their diverse roles in the universe (Abbott et al. 2017a; Villar et al. 2017) (see also references in Radice et al. 2018; Shibata & Hotokezaka 2019) since the first detection of gravitational waves (Abbott et al. 2017b). The dynamics of the merger and the evolution of the merger remnant are essential to pin down several elusive problems: the origin of heavy elements through r-process nucleosynthesis (Lattimer & Schramm 1974; Symbalisty & Schramm 1982; Eichler et al. 1989; Meyer 1989; Freiburghaus et al. 1999; Wanajo et al. 2014), kilonovas/macronovas (Li & Paczyński 1998; Kulkarni 2005; Metzger et al. 2010), and the possible connection to short, hard gamma-ray bursts (Eichler et al. 1989), to name a few. Mass ejection from the merger remnant is one of the keys to determine the composition of elements, the formation of a jet in the heated region, and the associated electromagnetic radiation.

Neutrinos can play significant roles in the dynamics through energy transport and the conversion of composition in the ejected material. In core-collapse supernovae, it is well known that they are essential to drive the core collapse, bounce, and explosion, which provide the burst of neutrino emission and their influence on the explosive nucleosynthesis (for a review, see Janka 2012). In neutron star mergers, neutrinos are influential in providing the heating around the merger remnant and modifying the compositional balance between neutrons and protons, which are crucial for r-process nucleosynthesis. The prediction of neutrino emission from neutron star mergers is also interesting as an observational signature (Ruffert et al. 1997; Rosswog & Liebendörfer 2003; Sekiguchi et al. 2011; Albert et al. 2017; Abe et al. 2018; Kyutoku & Kashiyama 2018).

Obtaining knowledge of the transport and reactions of neutrinos is therefore mandatory for investigating quantitatively the dynamics and properties of neutron star mergers and core-collapse supernovae, although their role has different impacts. The complication arises in highly deformed structures of the merger remnant composed of a massive neutron star and a torus, which can be opaque to neutrinos. Resulting neutrino emission can be largely asymmetric and the propagation of neutrinos may become nontrivial. It is problematic to describe the neutrino emission and absorption in and around such a highly deformed remnant. Neutrino emission can be aspherical and provides vastly different angle and energy spectra. Pair annihilation is known to be sensitive to anisotropic neutrino fluxes (Birkl et al. 2007; Liu et al. 2010; Zalamea & Beloborodov 2011), and different treatments may provide different magnitudes of heating. Irradiation by neutrinos with different energy spectra can drive a shift of electron fraction, which is the key quantity in r-process nucleosynthesis.

While the importance of neutrinos is obvious, numerical simulations of neutron star mergers have been performed with approximate or even simple methods of neutrino transport. This is because numerical-relativity simulations in full spatial dimensions are highly demanding with microphysics even with a simplified neutrino treatment. Moreover, long-term simulations are necessary to study gravitational waves from the merger, mass ejection from the merger remnant, and so on. Therefore, it is meaningful to assess the influence of neutrino reactions and transport in order to make further quantitative studies.

There have been continuous efforts in the numerical studies of neutron star mergers to implement neutrino transfer with various levels of approximations. Studies at the early frontiers (Ruffert et al. 1996; Rosswog & Liebendörfer 2003) adopt a leakage scheme for neutrino emissions, which was originally introduced in supernova studies (van Riper & Lattimer 1981), to consider the timescale of emission depending on the environment. Many modern studies utilize advanced versions of the leakage scheme to capture detailed variations (O'Connor & Ott 2010; Sekiguchi 2010; Galeazzi et al. 2013; Foucart et al. 2016; Perego et al. 2016; Radice et al. 2016; Fujibayashi et al. 2017, 2018; Ardevol-Pulpillo et al. 2019). Sophisticated treatments with the moment method with formulae for the closure relation (Shibata et al. 2011a; Sekiguchi et al. 2012; Just et al. 2015) and the Monte Carlo method (Richers et al. 2015; Foucart 2018; Miller et al. 2019; Foucart et al. 2020) are used to explore detailed dynamics in recent numerical simulations. Evaluation of the advantages and disadvantages of these methods has been conducted recently by comparisons of approximations in recent studies (Endrizzi et al. 2020).

The Monte Carlo method is one of the promising methods to provide detailed information on neutrino transfer (Janka & Hillebrandt 1989). It describes the solution of the Boltzmann equation by a sampling approach and provides the neutrino distribution in full dimensions. General relativistic calculations have recently been made (Foucart 2018; Miller et al. 2019) (see also Akaho et al. 2020) and applied to numerical simulations of neutron star mergers through coupling with hydrodynamics (Foucart et al. 2020) (see also Abdikamalov et al. 2012). The distribution in space, angle, and energy obtained by tracking the particle transport with reactions can be used to provide elaborate information on observational signals of supernova neutrinos (Keil et al. 2003; Kato et al. 2020) and to validate the closure relations for high angle moments (Foucart 2018; Foucart et al. 2018). It is advantageous to obtain accurate angle distributions at large distances in contrast to the limited ability of the discrete ordinate methods (Yamada et al. 1999). It is also suitable to describe the crossing of two beams, which is important for merger remnants (see, for example, Foucart 2018). On the other hand, it requires a large number of samplings to reduce the random noise, careful prescriptions to describe frequent reactions in the diffusion limit, and has only restricted capability in the rapidly varying background matter. In these respects, the Monte Carlo method is a complementary approach to the discrete ordinate method for the solution of the Boltzmann equation.

The direct evaluation of neutrino transfer by the Boltzmann equation is a very demanding computation in multidimensional space and has been applied to the merger of neutron stars in a limited study (Dessart et al. 2009) based on Ott et al. (2008). While general relativistic neutrino-radiation hydrodynamics with the Boltzmann equation under spherical symmetry (Yamada 1997; Yamada et al. 1999; Liebendörfer et al. 2004) has been applied in studies of core-collapse supernovae (Liebendörfer et al. 2001; Sumiyoshi et al. 2005, 2007; Fischer et al. 2010, 2011), neutrino-radiation hydrodynamics by the discrete ordinate (S_n) method for the Boltzmann equation for an axially symmetric case has been performed in a study by Ott et al. (2008). Numerical simulations of core-collapse supernovae in multiple dimensions have been extensively performed using sophisticated approximations of neutrino transfer (Janka 2012; Burrows 2013; Janka et al. 2016; Müller 2016). In order to follow the long-term evolution over several hundred milliseconds to find the outcome of an explosion with limited computational resources, there has been steady progress in numerical methods of neutrino transfer from the leakage scheme (van Riper & Lattimer 1981) to the diffusion approximations (Burrows et al. 2006a; Liebendörfer et al. 2009), variable Eddington factor (Rampp & Janka 2002), and two-moment schemes (Kuroda et al. 2012; Just et al. 2015). In addition, ray-by-ray methods (Buras et al. 2006) to efficiently handle directional variations have often been adopted to describe the multidimensional neutrino transfer.

Recently, the solution of neutrino transfer by the Boltzmann equation in six dimensions (6D) has become possible (Sumiyoshi & Yamada 2012; Sumiyoshi et al. 2015). The solver of the Boltzmann equation in 6D is applied to neutrino-radiation hydrodynamics in core-collapse supernovae under axial symmetry (Nagakura et al. 2014, 2016, 2017, 2018; Harada et al. 2019) and in three spatial dimensions (3D) (Iwakami et al. 2020). In the analysis by the 6D Boltzmann equation (Sumiyoshi et al. 2015), the feature of neutrino transfer in the profiles of core-collapse supernovae in three dimensions has been revealed with detailed information on neutrino distributions in neutrino angles and energy in 3D space.

In this study, we apply the numerical code to solve the 6D Boltzmann equation to investigate neutrino transfer in the matter profile taken from the numerical simulations of a neutron star merger. We performed numerical simulations for the fixed background to obtain the stationary state of neutrino distribution functions by solving the Boltzmann equation. We obtain all information on neutrino distribution functions in 5D (2D for space and 3D for neutrinos) to reveal the properties of neutrino transport and reactions. The simulation covers all regions—diffusion, intermediate, and transparent regimes—so that it is possible to analyze the behavior in a seamless manner.

We examine the basic feature of neutrino transfer: neutrino number density, flux distributions, angle moments, neutrino emission regions such as neutrinospheres, and angular variations of luminosities. We explore the neutrino emission inside an oblate remnant neutron star with a geometrically thick torus and the neutrino heating rates around them. We reveal that the neutrino fluxes are highly aspherical and focused in the region above the neutron star due to the deformed neutrinosphere elongated along the equator. We analyze the contributions of each neutrino reaction to find the important reaction for cooling and heating. We show that pair annihilation of μ and τ neutrinos plays an important role in neutrino heating in the region above the merger remnant.

We examine the properties of the Eddington tensor through comparisons between evaluation by the neutrino distribution functions and a closure relation, which is used in the moment formalism. We find that the closure relation provides a good approximation in most of the region, while it can be erroneous for high-energy neutrinos. Deviations of the Eddington tensor from the closure relation are seen in the regions where neutrino emission from both the neutron star and the torus is important.

The information on neutrino transfer obtained by the Boltzmann equation will be helpful in examining the approximations used in other studies of neutron star mergers. The present analysis will be also used to validate the approximate methods and to develop a new closure relation in future. We plan to compare the neutrino quantities such as neutrino luminosity in the Boltzmann evaluation and dynamical simulations more intensively in a separate paper.

We organize this paper as follows. We explain the profiles of the merger remnant used for the simulation of neutrino transfer in Section 2. We briefly describe the numerical treatment of neutrino transfer by solving the Boltzmann equation in Section 3. We report the basic features of neutrino transfer in the deformed remnant in Section 4. In Section 5, we analyze the angle moments (Section 5.1) and the Eddington tensor and compare them with those from the closure relation to examine the validity of the approximation (Section 5.2). We report the general trend of neutrino transfer in the profiles for selected time slices taken from the numerical-relativity simulation in Section 6. We summarize the paper in Section 7 with some discussions. We provide the details of the Eddington tensor in Appendix A and the angular resolution in Appendix B.

We utilize the matter profiles from the numerical-relativity simulations of binary neutron star mergers, which are composed of a massive neutron star with a surrounding torus. We take four snapshots from the time evolution in a 2D numerical simulation of radiation hydrodynamics under axial symmetry (Fujibayashi et al. 2017). The initial profile (0 ms) in the 2D simulation (Fujibayashi et al. 2017) is set up based on a result of 3D simulations (Sekiguchi et al. 2015) at ∼50 ms after the onset of the merger for the equal-mass model with a total mass of 2.7 M_⊙. The axially symmetric profile is constructed through an average over azimuthal angle from the 3D profile, which is quasi-stationary and nearly axisymmetric (Fujibayashi et al. 2017). We adopt the initial profile (0 ms) and the subsequent profiles at 30, 65, and 135 ms from the 2D simulation. The dynamical simulations were performed using the DD2 equation of state (EOS) by Banik et al. (2014). We adopt the same DD2 EOS for simulations of neutrino transfer.

In the 2D simulation, Einstein's equation is solved with a version of the puncture–Baumgarte–Shapiro–Shibata–Nakamura formalism (Shibata & Nakamura 1995; Baumgarte & Shapiro 1999; Marronetti et al. 2008). The radiation transfer equation for neutrinos is solved by employing a leakage-based scheme. In this scheme, neutrinos are separated into "trapped" and "streaming" neutrinos. The trapped neutrinos are assumed to couple the fluid tightly and are advected as a part of the fluid. To solve the streaming neutrinos, we employ Thorne's truncated moment formalism (Thorne 1981; Shibata et al. 2011b) in an energy-integrated manner with a closure relation (Levermore 1984; González et al. 2007). The detailed description of these schemes is found in Sekiguchi (2010) and Fujibayashi et al. (2017). For the weak interaction rates, we adopt the rates in Fuller et al. (1985) for the electron and positron capture processes, those in Cooperstein et al. (1986) for pair-production processes, those in Ruffert et al. (1996) for plasmon decay, and those in Burrows et al. (2006b) for nucleon–nucleon bremsstrahlung.

We use the central part of the profile up to 250 km for the simulation of neutrino transfer by the Boltzmann equation. We remap the rest-mass density, temperature, and electron fraction at the cell center of 128 and 48 grids for the radial and polar coordinates, respectively. (See details on numerical mesh in Section 3.) Note that we cover 90° of polar angle, and reflection symmetry is imposed with respect to Z = 0 as in the original simulations.⁵

We show the profile of merger remnant at 0 ms (∼50 ms after the merger) in Figure 1. The remnant neutron star at the center has an oblate shape with an extended torus. The region of high matter density (≳ 10¹¹ g cm⁻³) extends over 60 km along the equator. This extended region leads to trapping of neutrinos with energy of ∼10 MeV and contributes asymmetric neutrino emissions as we will see below. The high-temperature region appears at off-center locations that originate in the dynamics of neutron star merger. The electron fraction, which is equal to the total proton fraction under charge neutrality, is below 0.1 in the remnant neutron star and the torus. The electron fraction is high in the region above the merger remnant. The entropy per baryon is low in the merger remnant and high in the region above it. These thermodynamical and chemical conditions are inherited from the neutron star merger and the subsequent evolution.

Zoom In Zoom Out Reset image size

The deformed structure of the merger remnant persists in a quasi-static manner beyond 100 ms. Figure 2 shows the profiles of rest-mass density and temperature at 30, 65, and 135 ms. The remnant neutron star remains approximately stationary while the torus gradually shrinks and becomes more compact because of the cooling by neutrino emission. The peak of temperature is located in the off-center region. The rest-mass density above the neutron star gradually decreases due to the launch of mass ejection along the z-axis. The change in hydrodynamical quantities of the neutron star and the torus is relatively slow, therefore stationary simulations of neutrino transfer made by fixing the matter profiles as the background in each snapshot can be a good approximation.

Zoom In Zoom Out Reset image size

We adopt the numerical code to solve the 6D Boltzmann equation (Sumiyoshi & Yamada 2012). The code solves the time evolution of neutrino distribution functions in 6D (3D in the spherical coordinate system and 3D for two angles and the energy of neutrinos). The 6D Boltzmann equation is solved directly by the S_n method for multiple energy bins. We obtain the stationary state of neutrino distributions for fixed matter profiles by following their evolution for sufficient time periods over ∼10 ms. The same code is applied to study the neutrino transfer in 3D supernova cores, and details of the procedures can be found in Sumiyoshi et al. (2015). The solver of the 6D Boltzmann equation is extended to the special relativistic version (Nagakura et al. 2014); however, the velocity-dependent terms are dropped in the current study for fixed backgrounds. The simulation of the Boltzmann equation is done in flat space though we take into account the geometric factor for evaluation by volume integration (see for extensions toward the general relativistic simulation, Nagakura et al. 2017). The neutrino-radiation hydrodynamics using the solver of the relativistic 6D Boltzmann equation is applied to dynamical simulations of core-collapse supernovae in 2D (Nagakura et al. 2018, 2019a, 2019b; Harada et al. 2019, 2020) and in 3D (Iwakami et al. 2020).

We treat three species of neutrinos: ν_e, ${\bar{\nu }}_{e}$, and ν_μ. The μ-type neutrino, ν_μ, is a representative of the group of four species of (anti)neutrinos of heavy flavors. The basic set of neutrino reactions for supernovae is implemented in the collision term with angle- and energy-dependent expressions. The standard form of Bruenn's reaction rates (Bruenn 1985) is employed together with extended rates for pair processes (Sumiyoshi et al. 2005; Sumiyoshi & Yamada 2012). We note that the pair process is calculated using the distribution of the counterpart neutrinos. In order to make the Boltzmann equation linear in neutrino distribution function and reduce the computational cost, the counterpart neutrino distribution function at the previous step is inserted with the dependence of neutrino direction on polar angle (θ_ν) and the average over azimuthal angle (ϕ_ν) with respect to the radial coordinate. The use of distribution functions for the counterpart neutrinos at the previous step is validated in the study of 3D supernova cores because of the environment at high matter density and temperature. The table of the DD2 EOS by Banik et al. (2014) is used to match the original simulations.

The numbers of grid points for space are 128 and 48 for radial and polar angle coordinates. The 14 energy grid points are placed logarithmically to cover neutrino energies up to 300 MeV. High angular resolution is necessary to capture the neutrino transfer in the deformed neutron star with a torus. We set the numbers of grid points for neutrino angle as 56 and 12 for polar (θ_ν) and azimuthal (ϕ_ν) angles with respect to the radial coordinate. The angular resolution is higher than that in the case of 2D/3D core-collapse supernovae (Sumiyoshi et al. 2015). This is important for the convergence of angle distributions in nonspherical situations and integrated quantities such as heating rates along the z-axis. We describe the detailed examinations in Appendix B.

4.1. Neutrino Density, Flux, and Average Energy

The neutrino distributions have novel features due to the highly deformed shape of the merger remnant. The combination of the elongated neutron star and the extended torus contributes to the unique characteristic features of neutrino transfer. We show in Figure 3 the number density and flux of neutrinos for three species. The remnant neutron star abundantly contains neutrinos not at its center but off-center, where the temperature is high. Among the species, the electron-type antineutrino, ${\bar{\nu }}_{e}$, is the most abundant and the μ-type neutrino, ν_μ, is next. These species are produced by the pair process in such high-temperature environment. Electron-type neutrinos, ν_e, on the other hand, are produced through electron captures in an environment of high matter density. The distribution of ν_e is extended in the geometrically thick torus along the equator. The neutrino emission from the neutron star contributes to large neutrino fluxes in the region above it. Neutrino fluxes toward the equatorial region are suppressed since the geometrically thick torus is opaque to neutrinos. Neutrino emission from the torus contributes to the nonradial fluxes and leads to the concentration of neutrino fluxes above the merger remnant. Such a configuration results in the enhancement of heating in the region above the merger remnant (see Section 4.3).

Zoom In Zoom Out Reset image size

We display in Figure 4 the average energies of neutrinos in the central region. We show the first moment of energy in these plots. The average energies of three species are above 100 MeV in the high-temperature region. Moreover, the enhancement of the average energy for ν_μ in the region above the neutron star is noticeable. In the left panel of Figure 5, the average energies of neutrinos evaluated at the radius of 250 km, which is the outer boundary, are plotted as a function of polar angle. The average energies of neutrinos emitted from the merger remnant have a definite hierarchy among species. The average energy of ν_μ is distinctively higher than those of ν_e and ${\bar{\nu }}_{e}$ and shows a strong dependence on the polar angle.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 shows the radial and polar components of the energy fluxes for the three species. The fluxes are highly asymmetric with nonzero polar components. The radial component is dominant in the region above the neutron star for all species. The magnitude of the polar component is comparable to or even larger than that of the radial one. The energy fluxes of ν_e and ${\bar{\nu }}_{e}$ are widely extended due to the emission by charged current reactions from both the neutron star and the torus. The energy flux of ν_μ, on the other hand, is more focused due to the emission by pair productions from the high-temperature region of the neutron star with minor contributions from the torus. The contribution of the torus as a peripheral source is most remarkable for ν_e, moderate for ${\bar{\nu }}_{e}$, and minor for ν_μ, as we will see in the emission rates in Figure 11 (Section 4.4). In addition, the torus plays the role of a shield toward the equator for ν_μ having only the central source.

Zoom In Zoom Out Reset image size

In the right panel of Figure 5, the radial energy fluxes of neutrinos evaluated at the radius of 250 km are plotted as a function of polar angle. The dependence of energy fluxes on the polar angle is strong, especially for ν_μ. Integration of the radial energy fluxes over solid angle provides the total luminosities for all directions of 3.2 × 10⁵², 4.9 × 10⁵², and 4.2 × 10⁵² erg s⁻¹ for ν_e, ${\bar{\nu }}_{e}$, and ν_μ, respectively. The luminosities of ${\bar{\nu }}_{e}$ and ν_μ are high due to the thermal processes, especially nucleon–nucleon bremsstrahlung. The luminosities in Fujibayashi et al. (2017) for each species are  ∼2 × 10⁵²,  ∼3 × 10⁵², and  ∼6 × 10⁵¹ erg s⁻¹, respectively. The luminosity is larger by several tens per cent for ν_e and ${\bar{\nu }}_{e}$, and remarkably by a factor of ∼7 for ν_μ. The large difference in the luminosity for ν_μ is also found in Foucart et al. (2020), in which the neutrino transfer with the Monte Carlo method is compared to the moment formalism in binary neutron star mergers. Thus, this difference would be because of the different level of the sophistication of neutrino transfer. We will investigate the difference closely in our future work.

4.2. Neutrinosphere

The region of neutrino emission has a deformed shape due to the deformed nature of the merger remnant. The locations of the neutrinospheres are shown in Figure 7 on top of the rest-mass density profiles. The neutrinosphere is defined by the location for which the optical depth is 2/3. Here, the optical depth is evaluated as the integral of the inverse mean free path along the radial coordinate from outside using the mean free path for forward angle.⁶

Zoom In Zoom Out Reset image size

The shapes of neutrinosphere for three species are highly deformed, reflecting the matter distribution. The neutrinosphere for a neutrino energy of 13 MeV roughly follows the contour of rest-mass density at around 10¹¹–10¹² g cm⁻³. This explains why the neutrino fluxes are focused in the region above the neutron star. It is remarkable that the neutrinosphere extends to the regions of low temperature and positive neutrino chemical potential.

The degree of deformation and the equatorial extension of a neutrinosphere depends strongly on the neutrino energy. Figure 8 displays the locations of neutrinospheres for different energies. The shapes of the neutrinospheres are more compact for 4.9 MeV and more extended for 34 and 89 MeV than those for 13 MeV in Figure 7. The positions of neutrinospheres above the remnant neutron star are rather close for different energies because of a steep matter density gradient. Fluxes of high-energy neutrinos tend to be focused along the z-axis as a consequence. The strong dependence of neutrino flux on the energy may influence the composition of ejecta and the resulting product of nucleosynthesis.

Zoom In Zoom Out Reset image size

4.3. Heating Rates

Energy transfer to matter through neutrino reactions plays an important role in the ejection of mass from the merger remnant and the resulting nucleosynthesis. We show a contour map of the specific heating rate in Figure 9. The heating is most efficient in the region above the remnant neutron star. The specific heating rate is larger than 10²³ erg g⁻¹ s⁻¹ at a radius around 10–40 km. The region of strong heating extends over 45°. The cooling proceeds in the limited region of the neutron star and torus. The total heating rate amounts to 1.2 × 10⁵² erg s⁻¹ according to the volume integral of local heating rates in the heating region. The corresponding value in Fujibayashi et al. (2017) is ∼3 × 10⁵¹ erg s⁻¹, and thus the heating rate in the simulation of the Boltzmann equation is ∼4 times larger than that in the simulation with an energy-integrated leakage-based scheme. This difference is due to the stronger pair annihilation heating because of the larger luminosity in the simulation of the Boltzmann equation (see Sections 4.1 and 4.4).

Zoom In Zoom Out Reset image size

We show in Figure 10 the specific heating rate for each neutrino reaction. The heating by the charged current reactions with neutrons and protons extends in a region above the neutron star as well as in the torus for electron-type (anti)neutrinos. Heating by the pair annihilation reaction is effective in a narrow region above the neutron star for which the matter density is low and neutrino fluxes are high for three species. The contribution of heating by nucleon–nucleon bremsstrahlung is minor (not shown here) though it is important for emissions (see below).

Zoom In Zoom Out Reset image size

4.4. Emission

The intense heating in the region above the neutron star is caused by the influence of neutrino emission from the merger remnant. We examine the neutrino emission rates for several reactions in Figure 11. For ν_e and ${\bar{\nu }}_{e}$, the charged current reactions with nucleons are dominant. The region of strong emissivity extends into the torus especially for ν_e due to the contribution of electron capture by protons. Neutrino emission by pair creation proceeds mainly through nucleon–nucleon bremsstrahlung for all species, which is more efficient than electron–positron pair annihilation. For μ- and τ-type neutrinos, the emission of pairs (ν_μ, ${\bar{\nu }}_{\mu }$, ν_τ, and ${\bar{\nu }}_{\tau }$) also proceeds dominantly through nucleon–nucleon bremsstrahlung. These pairs of neutrinos contribute to the heating through neutrino pair annihilation in the region of low matter density, while electron-type (anti)neutrinos additionally contribute to the heating through absorption by nucleons away from the region where the matter density is relatively high, as shown above.

Zoom In Zoom Out Reset image size

5.1. Angle Moments

It is advantageous that the direct solution of the Boltzmann equation can provide the angle distributions for multiple energies over the whole spatial grids. We examine the angle moments from the neutrino distribution function in five dimensions. As an example, Figure 12 displays the average value of the squared angle moment, $\langle {\mu }_{\nu }^{2}\rangle$, where the angle factor is defined as ${\mu }_{\nu }=\cos {\theta }_{\nu }$, for three species. The region with $\langle {\mu }_{\nu }^{2}\rangle =\tfrac{1}{3}$ (isotropic) is not spherical but largely deformed including the torus along the equatorial plane. The shape of contour lines follows the torus near the equator and becomes an oblate spheroid in the outer region. The squared angle moment increases toward 1 (forward peak) at large distances.

Zoom In Zoom Out Reset image size

Figure 13 displays the behavior of the radial distributions of angle moments, 〈μ_ν〉 and $\langle {\mu }_{\nu }^{2}\rangle$ along directions at three polar angles. It is interesting to see that these quantities behave in a non-monotonic manner depending on the direction. The angle moments along the equator (right panel) remain at their isotropic values due to the neutrino trapping in the merger remnant and simply increase to the forward peak value. There are some wiggles within 10 km due to inward and outward flows caused by neutrino diffusion near the shell-like high-temperature region inside the neutron star (see upper right panel of Figure 1). Above the neutron star (left panel), the behavior is similar, but the transition occurs around 10 km. Along the direction above the extended torus (middle panel), the angle moments show a complicated increase and decrease from isotropic to forward peak due to the shape of the torus. The angle moment 〈μ_ν〉 approaches 1 at the outer edge, 250 km. We note that the convergence of the forward peak value of the angle moments at large distances is achieved in the current simulation with high angular resolution. The convergence with various numbers of angle grids has been checked as described in Appendix B.

Zoom In Zoom Out Reset image size

5.2. Eddington Tensor

We analyze the Eddington tensor obtained from the simulation of the Boltzmann equation. Here, the Eddington tensor can be directly evaluated by integrating the neutrino distribution functions (hereafter, the direct evaluation). We compare it with the Eddington tensor from the closure relation in order to assess the validity of this approximation used in the moment formalism. In the closure relation, the pressure tensor is evaluated by assuming relations with the lower moments: energy flux and energy density. The definition of the Eddington tensor and the procedure to evaluate it using the closure relation can be found in Appendix A.

Figure 14 displays the rr-component of the Eddington tensor for a neutrino energy of 34 MeV obtained by the direct evaluation and by the closure relation, and the difference between them, for three species. The profile of the rr-component of the Eddington tensor from the direct evaluation (top) is deformed according to the shape of the merger remnant. The region with a value of $\tfrac{1}{3}$, which corresponds to the isotropic distribution, is roughly located inside the neutrinosphere (see middle panel of Figure 8). It increases to the limiting value of 1 (forward peak) away from the remnant. The profile obtained by the closure relation (middle) in general resembles the profile from the direct evaluation. The difference in rr-components between the closure relation and the direct evaluation is shown in the bottom panels. Note that the plots show the absolute difference, not the relative one, from the direct evaluation. Although the difference is small in most of the region, there is a sizable difference (∼0.1) in the region close to both the neutron star and the torus, where neutrinos from both are important.

Zoom In Zoom Out Reset image size

Figure 15 shows θθ-, ϕϕ-, and rθ-components of the Eddington tensor for ν_e. The profiles of the θθ- and ϕϕ-components are similar to each other and also to that of the rr-component. The θθ- and ϕϕ-components have an isotropic value of $\tfrac{1}{3}$ inside the merger remnant and decrease to 0 at large distances. It is remarkable that the rθ-component has nonzero finite values in rather outer regions. The difference between the closure relation and the direct evaluation is generally small but noticeable in the region above the torus, as in the case of the rr-component.

Zoom In Zoom Out Reset image size

The feature of the Eddington tensor with deformed shape substantially depends on the neutrino energy, as we have seen for the large deformation of the neutrinosphere for high energies in Figure 8. The validity of the closure relation accordingly depends on the energy range for each species. For low energies, the degree of deformation is rather small, reflecting the small deformation of the neutrinosphere in the remnant. Figure 16 shows the rr-component of the Eddington tensor for ν_e with neutrino energies of 4.9 and 13 MeV. For these energies, the deviation of the closure relation is small around the torus. The overestimation above the remnant neutron star is noticeable instead.

Zoom In Zoom Out Reset image size

Figure 17 shows the rr-component of the Eddington tensor for three species with a neutrino energy of 89 MeV, which is higher than those examined in Figures 14–16. The deformed profile of low values (isotropic, i.e., $\sim \tfrac{1}{3}$) is extended further and has a large and thick shape, inside which there are some regions of non-isotropic values. This is also a consequence of the more deformed shape of the neutrinosphere for higher energies. The shape also depends on the neutrino species, showing a smaller torus for μ-type neutrinos than for the others. While the rr-component for electron-type neutrinos rapidly increases to 1 along the z-axis, the increase is slow for μ-type neutrinos. Due to the large deformation, the deviation of the closure relation is apparently large for high energies as seen in the bottom panels. The region of underestimation extends around the edge of the torus. Further studies are necessary to clarify these differences around the torus.

Zoom In Zoom Out Reset image size

These trends of deformed shapes depending on the energy and species are well in accord with those of the shape of neutrinosphere seen in Figures 7 and 8. In this respect, it requires more caution in approximations of neutrino transfer for higher energies due to large anisotropies than for lower energies. This also suggests the importance of the numerical simulation for multiple energies, which treats the neutrino transport at different energies.

The radial distributions of the Eddington tensor at different polar angles behave in various ways depending on the traversing profile of the remnant. Figure 18 shows the components of the Eddington tensor as functions of radial distance at the three polar angles (θ =  9.8°, 48°, and 90°) for a neutrino energy of 89 MeV. While the tensor varies rapidly near the z-axis (left panel) from the isotropic value in the merger remnant⁷ to the forward-peaked value in the ambient dilute material, it behaves in a non-monotonic way at the equator (right) and around the edge of the torus (middle) because of the deformed profile. The polar component of the flux of neutrinos contributes to sizable non-diagonal rθ-components and affects the behavior of diagonal components. The components of the Eddington tensor obtained by the closure relation (thin lines) follows most of the trends in the components obtained by the direct evaluation (thick lines), but there are some regions of substantial difference as shown in the bottom panels. The deviation can amount to over 0.1 in the absolute value of the Eddington tensor, which is not negligible for the pressure evaluation of the moment scheme in the transitional regime between diffusive and transparent situations.

Zoom In Zoom Out Reset image size

We examine the properties of neutrino transfer in the time sequence of the remnant of the neutron star merger at 30, 65, and 135 ms (see also Figure 2). As we will see, the feature of neutrino transfer remains similar over ∼100 ms, but the asymmetric features of neutrino quantities seen in the initial profile becomes gradually milder due to the shrinkage of the torus, and in particular of its high-temperature region.

Figure 19 displays the number density and flux of neutrinos for ν_e (left), ${\bar{\nu }}_{e}$ (middle), and ν_μ (right) for the profiles at 30, 65, and 135 ms. The asymmetric emission of neutrinos continues up to 135 ms. However, the region with high neutrino number densities in the torus gradually shrinks according to the evolutionary change due to the cooling through neutrino emission (see Figure 2). The reduction of temperature induces this trend as we can see in the region in and around the torus in the late profiles, although the structure of the matter density profile changes slowly. (See the bottom panels in Figure 2.) The neutrino emission from the remnant is gradually reduced accordingly.

Zoom In Zoom Out Reset image size

The feature in the time evolution of neutrino flux depends on the neutrino species. The fluxes of ${\bar{\nu }}_{e}$ and ν_μ are focused above the neutron star and this tendency is enhanced in the late profiles. We show the radial component of the neutrino energy fluxes at 135 ms for three species as contour plots in Figure 20. Compared with the situation in Figure 6, the radial energy flux is focused in a narrower region away from the equator. This tendency is stronger in ${\bar{\nu }}_{e}$ and ν_μ than in ν_e because the emission of ${\bar{\nu }}_{e}$ and ν_μ is confined to the high-temperature region of the neutron star whereas the emission of ν_e continues also from the torus. However, the influence of the torus on the neutrino transfer in general becomes minor as the extension becomes small.

Zoom In Zoom Out Reset image size

Figure 21 shows the average energy and luminosity at the outer boundary for the four snapshots. The average energy is evaluated as the average over the solid angle. The luminosity is the integral over the solid angle and scaled for 4π coverage. The average energies for ν_μ and ${\bar{\nu }}_{e}$ gradually increase whereas that for ν_e remains almost the same. On the other hand, the luminosity for ν_μ decreases slowly and those for ν_e and ${\bar{\nu }}_{e}$ decrease rapidly.

Zoom In Zoom Out Reset image size

The evolutionary trend of neutrino fluxes is related with the shrinkage of the region of neutrino emission and the extended structure of the neutrinosphere. Figure 22 shows the evolution of the neutrinosphere. Its shape remains extended to the equator due to the presence of the torus, although the degree of the extension gradually shrinks as the torus shrinks. As the matter density above the neutron star becomes low and that around the equatorial region remains high even in the late stages, the flux above the neutron star can be enhanced due to the opaque condition in the equatorial directions. Since the production of thermal neutrinos (${\bar{\nu }}_{e}$ and ν_μ) continues in the high-temperature region, the neutrino emission region becomes almost confined to the neutron star and the neutrino flux is focused due to the hindrance by the extended neutrinosphere along the equatorial directions. In contrast, the production of electron neutrinos continues through charged current reactions in both the neutron star and the torus, and thus the neutrino flux is not so highly focused.

Zoom In Zoom Out Reset image size

It is interesting to examine the time evolution of the heating rate above the neutron star merger as a source of mass ejection (and associated radiation) in the region above the remnant. We show the evolution of heating and cooling rates as contour plots in Figure 23. It is remarkable that the heating in the region just above the neutron star persists over 100 ms. The total heating rate over the volume scaled for the 4π coverage is plotted as a function of time in Figure 24. The total heating rate amounts to over 5 × 10⁵¹ erg s⁻¹ and stays constant even in the late phase. This continuation of strong heating may influence the ejection of material along the z-axis. The corresponding heating rate in Fujibayashi et al. (2017) is ∼(6–8) × 10⁵⁰ erg s⁻¹ for t ∼ 0.05–0.15 s, and thus the heating rate in the late phase could be ∼7 times larger with a more sophisticated neutrino transfer method. Detailed comparison of the two simulations has to be made with considerations of different numerical schemes as well as reaction rates, and a separate study will be reported elsewhere.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We studied the properties of neutrino transfer in a remnant system of a binary neutron star merger consisting of a massive neutron star and a torus by solving the multidimensional Boltzmann equation. Adopting snapshots of the merger remnant obtained in numerical-relativity simulations (Sekiguchi et al. 2015; Fujibayashi et al. 2017) as the background profile of matter, we performed numerical simulations of neutrino propagation with reactions to obtain the neutrino distributions in full dimensions with energy and angle dependence. The solution of the Boltzmann equation enabled us to examine the properties of neutrino transfer in such a deformed structure in detail.

We revealed that neutrino transfer is highly aspherical in the profiles of a massive neutron star with an elongated torus. We show that neutrinos are trapped inside the deformed neutron star and the torus extends along the equator. The emission of neutrinos proceeds from both the neutron star and the extended torus. The neutrino fluxes are focused above the neutron star since the torus is geometrically thick along the equator.

We showed that the shape of the neutrinosphere is largely extended in the deformed neutron star with a geometrically thick torus. The location of the neutrinosphere is determined by the matter distribution of the merger remnant with a dependence on neutrino species and energy. The large deformation of the neutrinosphere for high-energy neutrinos enhances the focused flux above the neutron star near the z-axis.

The aspherical features of neutrino transfer differ depending on the neutrino species due to different origins. While the energy flux of ν_e is widely extended, the energy fluxes of ${\bar{\nu }}_{e}$ and ν_μ are more focused above the neutron star. This is because the emission of ν_e proceeds in both the neutron star and the torus whereas the emission of ${\bar{\nu }}_{e}$ and ν_μ occurs only in the high-temperature region. In addition, the energy flux is focused in the region above the neutron star that has the optically thick condition of the torus in the equatorial region. Therefore, the neutrino energy fluxes have different angular variations depending on the species.

We found that neutrino heating proceeds dominantly just above the remnant neutron star by the neutrino energy fluxes focused along the z-axis. The heating occurs through the charged current reactions with nucleons and the pair annihilation of neutrinos. The neutrinos for heating originate from the emission of the trapped neutrinos in the deformed remnant. While the charged current reactions are the main process for electron-type neutrinos, nucleon–nucleon bremsstrahlung is additionally important for the emission of neutrino pairs, especially for μ- and τ-type neutrinos.

Our simulations of the Boltzmann equation enable us to directly evaluate the angle moments and the Eddington tensors from the full information of space, angle, and energy distributions of neutrinos. We studied the angle moments and tensor components by the direct integration of the neutrino distribution functions to validate the approximate methods of neutrino transfer. We compared the Eddington tensor obtained by the closure relation with those obtained by the direct evaluation. We found that the difference between the two methods is generally small, but it becomes large for high-energy neutrinos due to the contributions of fluxes from both the neutron star and the torus.

We found that the general trend of neutrino transfer is preserved for ∼100 ms after the merger by examining the neutrino transfer and adopting a series of snapshots as the background. As the torus along the equator shrinks, the aspherical features of neutrino emission become less drastic. The neutrino heating above the neutron star persists for ∼100 ms in these snapshots.

The direct evaluation of the neutrino distribution functions from the Boltzmann equation provides new information for the treatment of neutrino transfer in numerical simulations of the remnant of neutron star mergers. The validation by the direct evaluation would be helpful to improve the approximate methods such as the closure relations to provide the angle-averaged quantities in a truncated moment scheme and to gauge parameters of neutrino emission and absorption in simplified schemes.

We examined the convergence due to the resolution of angle bins in neutrino transfer. We demonstrated that it is important to have high angular resolution, as adopted in the current study, for the convergence of heating rates due to fluxes of neutrino pairs.

It would be interesting to pursue the evolution of neutrino transfer in a longer evolution of the remnant of the neutron star merger. Continuous neutrino heating may affect the mass ejection from the merger remnant. In its long-term evolution, the magnetic field is amplified in the neutron star and torus by magnetohydrodynamical instabilities (e.g., Balbus & Hawley 1991; Price & Rosswog 2006). The effective viscosity originating from the instabilities transports the angular momentum in the system. The viscosity also heats up the material in the remnant. In this way, it may modify the matter density and temperature structures of the remnant and affect the neutrino emission.

It will also be important to study differences of neutrino transfer further in various methods to painstakingly solve the Boltzmann equation. The Monte Carlo method, for example, reveals the asymmetric neutrino emission and provides neutrino energy spectra and angle distributions (Richers et al. 2015; Foucart et al. 2020), which enables examinations of the closure relation used in the M1 scheme (Foucart 2018) as in the current study. While these studies reveal the luminosities and average energies of neutrino emission and the heating rates through pair annihilation, it is necessary to investigate quantitative differences that arise from the methods with different settings and microphysics.

We are preparing for the detailed comparison of neutrino transfer with hydrodynamical simulations so that we can validate and/or improve the approximate methods necessarily adopted for the long-term evolution of numerical-relativity simulations. These analyses will be separately reported elsewhere.

This work is supported by Grant-in-Aid for Scientific Research (15K05093, 16H06341, 16K17706, 18K03642, 19K03837, 20H00158, 20H01905) and Grant-in-Aid for Scientific Research on Innovative areas "Gravitational wave physics and astronomy:Genesis" (17H06357, 17H06361, 17H06363, 17H06365) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan.

For providing high-performance computing resources, Computing Research Center, KEK, JLDG of NII, Research Center for Nuclear Physics, Osaka University, Yukawa Institute of Theoretical Physics, Kyoto University, Information Technology Center, University of Tokyo, and Max Planck Computing and Data Facility are acknowledged.

This work was supported by MEXT as "Program for Promoting Researches on the Supercomputer Fugaku" (Toward a unified view of the universe: from large scale structures to planets, Simulation for basic science: from fundamental laws of particles to creation of nuclei), JICFuS, and the Particle, Nuclear and Astro Physics Simulation Program (No. 2019-002, 2020-004) of Institute of Particle and Nuclear Studies, High Energy Accelerator Research Organization (KEK).

We briefly describe here the definition of angle and energy moments and the Eddington tensor with associated approximations. The energy density, flux, and pressure tensor of neutrinos are evaluated by the neutrino distribution function, f(ε, Ω), for energy, ε, and angle, Ω, as

$$\begin{eqnarray}&&E=\int \displaystyle \frac{d\varepsilon {\varepsilon }^{2}}{{\left(2\pi \right)}^{3}}\int d{\rm{\Omega }}\varepsilon f(\varepsilon ,{\rm{\Omega }}),\end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray}&&{F}^{i}=\int \displaystyle \frac{d\varepsilon {\varepsilon }^{2}}{{\left(2\pi \right)}^{3}}\int d{\rm{\Omega }}\,\varepsilon {n}_{i}f(\varepsilon ,{\rm{\Omega }}),\end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray}&&{P}^{{ij}}=\int \displaystyle \frac{d\varepsilon {\varepsilon }^{2}}{{\left(2\pi \right)}^{3}}\int d{\rm{\Omega }}\,\varepsilon {n}_{i}{n}_{j}f(\varepsilon ,{\rm{\Omega }}),\end{eqnarray} \tag{ A3 }$$

respectively. Hereafter we use natural units with ℏ = c = 1. Here we suppress the spatial coordinate, r, θ, and ϕ, in the neutrino distribution function. n_i and n_j are unit vectors for the three components of r, θ, and ϕ directions, which are denoted by the subscripts i and j. The detailed definition of variables and the unit vectors in the spherical coordinate system can be found in Sumiyoshi & Yamada (2012).

It is essential to handle the moment equations with multiple energies, for example, for core-collapse supernovae. We analyze the moments for the energy zone in the discretized energy variable. The energy density for the energy zone, ε_k, with width Δε_k, is expressed as

$$\begin{eqnarray}&&E({\varepsilon }_{k})=\displaystyle \frac{{\rm{\Delta }}{\varepsilon }_{k}{\varepsilon }_{k}^{2}}{{\left(2\pi \right)}^{3}}\int d{\rm{\Omega }}\,{\varepsilon }_{k}f({\varepsilon }_{k},{\rm{\Omega }}).\end{eqnarray} \tag{ A4 }$$

The pressure tensor for the energy zone is similarly expressed as

$$\begin{eqnarray}&&{P}^{{ij}}({\varepsilon }_{k})=\displaystyle \frac{{\rm{\Delta }}{\varepsilon }_{k}\,{\varepsilon }_{k}^{2}}{{\left(2\pi \right)}^{3}}\int d{\rm{\Omega }}\,{\varepsilon }_{k}{n}_{i}{n}_{j}f({\varepsilon }_{k},{\rm{\Omega }}).\end{eqnarray} \tag{ A5 }$$

The Eddington tensor for each energy zone can be evaluated from

$$\begin{eqnarray}&&{k}^{{ij}}({\varepsilon }_{k})=\displaystyle \frac{\int d{\rm{\Omega }}{n}_{i}{n}_{j}f({\varepsilon }_{k},{\rm{\Omega }})}{\int d{\rm{\Omega }}\,f({\varepsilon }_{k},{\rm{\Omega }})},\end{eqnarray} \tag{ A6 }$$

from the neutrino distribution function. Since we describe the neutrino distribution function by the Boltzmann equation in 6D, we can directly provide these angle moments.

In the moment formalism with termination at a certain rank, the highest moment has to be determined by lower moments. For example, it is necessary to provide the second moment, i.e., the pressure tensor, through a closure relation by a functional form of the energy density and flux. In a classic form of the closure relation by Levermore (1984), the closure relation is given by flux vectors via a form of the tensor, t^ij(ε_k),

$$\begin{eqnarray}&&{P}^{{ij}}({\varepsilon }_{k})={t}^{{ij}}({\varepsilon }_{k})E({\varepsilon }_{k}),\end{eqnarray} \tag{ A7 }$$

to determine the pressure tensor. The tensor form of closure is given by

$$\begin{eqnarray}&&{t}^{{ij}}({\varepsilon }_{k})={\delta }^{{ij}}\displaystyle \frac{1-\chi }{2}+{u}^{i}{u}^{j}\displaystyle \frac{3\chi -1}{2}.\end{eqnarray} \tag{ A8 }$$

uⁱ is the unit vector defined by

$$\begin{eqnarray}&&{u}^{i}=\displaystyle \frac{{f}^{i}}{f},\end{eqnarray} \tag{ A9 }$$

using the normalized flux

$$\begin{eqnarray}&&{f}^{i}=\displaystyle \frac{{F}^{i}({\varepsilon }_{k})}{E({\varepsilon }_{k})},\end{eqnarray} \tag{ A10 }$$

and its norm, f. χ is the Eddington factor, which smoothly connects the two limiting cases from the diffusion regime (χ = 1/3) to the transparent regime (χ = 1). We adopt the expression as a function of f by Levermore (1984):

$$\begin{eqnarray}&&\chi =\displaystyle \frac{1}{3}+\displaystyle \frac{2{f}^{2}}{2+\sqrt{4-3{f}^{2}}},\end{eqnarray} \tag{ A11 }$$

in the current study. We examine the difference in the tensor form of closure from the Eddington tensor

$$\begin{eqnarray}&&{\rm{\Delta }}{t}^{{ij}}={t}^{{ij}}-{k}^{{ij}}\end{eqnarray} \tag{ A12 }$$

to check the quality of approximation by the closure relation.

In order to assess the convergence of the grid resolution in angle mesh for neutrinos, we perform a series of simulations for additional selected models by changing the number of grids for angles. We study the dependence of the behavior of moments, tensor, and heating rate on the angular resolution for the models of the central core in the gravitational collapse of massive stars. As we show below, the number of angle grids we adopted is sufficient for the current purpose to study the deformed neutron star with torus. We take the procedure to obtain the stationary state of the neutrino distribution function by solving the Boltzmann equation with the fixed background profile chosen from hydrodynamics. We utilize the same setting of neutrino reactions as described in Section 3 using the Shen EOS (Shen et al. 1998a, 1998b) matched with hydrodynamics. We adopt the same number of grids for space as well as neutrino energy and azimuthal angle as in the case of profiles from neutron star mergers, but change the number of neutrino polar angle grids, ${N}_{{\theta }_{\nu }}$, from 6, 12, 18, 24, 30, 36, 40, 44, 48, 52, to 56.

A set of profiles is taken from the 2D numerical simulation (Sekiguchi et al. 2012; Kotake et al. 2012) of rotating collapse of a 100 M_⊙ star by Umeda & Nomoto (2008). The core-collapse simulation revealed the dynamics of a collapsar: the formation of neutron star with a geometrically thick torus from the rotating massive star as a model of long gamma-ray bursts. We adopt a profile of the central part including a neutron star with the torus-like extension formed after the core bounce. Due to the rotation given in the initial condition, the massive proto-neutron star has an elongated shape with extended matter distribution. This situation resembles the remnant in a neutron star merger. Therefore, the examination of angular resolution is applicable to the situation in the current study.

We show, in Figure 25, the profiles of hydrodynamics quantities in the collapsar. At the center, the massive proto-neutron star is born with the extended torus at high matter density and temperature. The torus has a geometrically thick structure along the equator and the material is dilute along the z-axis. The detailed information about the dynamics and snapshots can be found in Sekiguchi et al. (2012) and in Section 3.3.3 of Kotake et al. (2012). We demonstrate the profiles of neutrino quantities obtained by the simulations of the Boltzmann equation with the highest angular resolution in Figure 26. The neutrino distributions and fluxes are nonspherical, having abundant neutrinos in the extended neutron star with enhanced fluxes along the z-axis. The neutrino distributions is deformed with a nearly isotropic region ($\langle {\mu }_{\nu }^{2}\rangle \sim \tfrac{1}{3}$) extended along the equator.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We show in Figure 27 radial distributions of angle moments, 〈μ_ν〉 and $\langle {\mu }_{\nu }^{2}\rangle$, at two polar angles. The angle moments increase smoothly along the z-axis (left) in the situation from the neutron star to the dilute material, while they behave in a non-monotonic way along the edge of the geometrically thick torus. The convergence of the curve of angle moments by increasing the number of angle grids is apparent. In order to show the variations of the angle moments in detail, we plot the values at the outermost radial position as a function of the number of grids in Figure 28. The angle moments change rapidly for a small number of grids, but converge adequately for the cases of over 40 grids. We have checked that the energy moments converge much faster than the angle moments.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We show in Figures 29 and 30 the radial distributions of the components of the Eddington tensor for three species along the radial directions at 69° from the z-axis (the edge of the geometrically thick torus) for different numbers of angle grids. The components of the Eddington tensor behave in a non-monotonic manner for ν_e and in a monotonic manner for ${\bar{\nu }}_{e}$ and ν_μ. The overall feature remains similar for ${\bar{\nu }}_{e}$ and ν_μ and is sensitive to the number of grids for ν_e. It converges well for the cases with 48 and 56 angle grids even in the case for ν_e. Relative differences in the values for 48 and for 56 angle grids are within 3% and 7% for diagonal and non-diagonal components, respectively. Hence, we believe that the number of angle grids in the current study of neutron star mergers is sufficient to examine the quantities of the Eddington tensor in detail.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The heating rate is important for the hydrodynamics such as mass ejection, which may contribute to nucleosynthesis as well as jet formation. The amount of heating through the pair annihilation of neutrinos and antineutrinos, which is a major contribution in this situation, is sensitive to the neutrino fluxes with detailed angular distributions. Therefore, the angular resolution is crucial to reliable evaluation. We show in Figure 31 the contour plot of heating rates for different numbers of angle grids. The heating is dominant in the region just above the neutron star along the z-axis in a similar situation to the case of a neutron star merger. It is noticeable that stronger heating proceeds in the region above the neutron star for the case of a small number of angle grids. The heating distribution converges for a large number of angle grids in the lower panels. To see the convergence in a quantitative manner, we show in Figure 32 the total heating rate for different numbers of angle grids. The total heating rate is larger for lower angular resolution and decreases for higher angular resolution. It nearly converges for the finest angular resolution in the current study at 56. The relative difference in the values with 48 and with 56 angle grids is 1.3%. This is related to the sensitivity of the pair-annihilation rate to the angle distribution. It is therefore advisable to provide enough angular resolution to determine the heating rate in a reliable manner for the situation of a neutron star with an extended torus such as a collapsar and a neutron star merger.

Zoom In Zoom Out Reset image size

Slow convergence of the heating rate is caused by the extended geometry of material as seen in this study. This is different from ordinary situations in core-collapse supernovae where moderate angular resolution is tolerated to determine the heating rate to some extent. Here we additionally demonstrate the dependence of the total heating rate for the supernova core on the angular resolution for comparison. Note that high angular resolution is necessary to obtain the forward peak at large distances for all applications in principle. The detailed study of angular resolution of neutrino transfer by the Boltzmann equation is reported in Richers et al. (2017) and Iwakami et al. (2020).

We adopt a profile of a supernova core at 150 ms after the core bounce taken from the core-collapse simulations in 2D (Takiwaki et al. 2012, 2014; Horiuchi et al. 2014) of an 11.2 M_⊙ star (Woosley et al. 2002). The case of an 11.2 M_⊙ star, which leads to explosions in the 2D and 3D simulations by Takiwaki et al. (2012), is an example of the configuration of a propagating shock wave with a largely elongated shape. We obtained the stationary neutrino distributions by solving the Boltzmann equation by the same procedure as in Sumiyoshi et al. (2015) and by the one above. The EOS by Lattimer & Swesty (1991) with the incompressibility of 180 MeV is used to match the original simulations. We adopt the number of polar angle grids of neutrinos from 6, 12, 18, 24, 30, to 36 with the same settings of 256 for radial and 64 for polar grids in space and of 14 for energy and 12 for azimuthal angle of neutrinos.

The basic properties of neutrino transfer in the 3D profiles of supernova cores for the cases of 11.2 M_⊙ and 27 M_⊙ are reported in Sumiyoshi et al. (2015). The properties of neutrino distributions in the 2D supernova core are investigated in Abbar et al. (2019, 2020) for applications to collective neutrino oscillations.

We show the total heating rates in Figure 33 for different numbers of angle grids. It is evident that the total heating rate does not change very much with increasing number of angle grids. Although there is a slight overestimation with the case of six angle grids, it remains flat for larger values. This is clearly different from the situation in a collapsar seen in Figure 32. Hence, it is not crucial to increase the number of angle grids for the total heating in core-collapse supernovae.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size