Evaluating radiation transport errors in merger simulations using a Monte Carlo algorithm

Neutrino-matter interactions play an important role in the postmerger evolution of neutron star-neutron star and black hole-neutron star mergers. Most notably, they determine the properties of the bright optical/infrared transients observable after a merger. Unfortunately, Boltzmann’s equations of radiation transport remain too costly to be evolved directly in merger simulations. Simulations rely instead on approximate transport algorithms with unquantified modeling errors. In this paper, we use for the first time a time-dependent general relativistic Monte Carlo (MC) algorithm to solve Boltzmann’s equations and estimate important properties of the neutrino distribution function $\sim 10\mathrm{ms}$ after a neutron star merger that resulted in the formation of a massive neutron star surrounded by an accretion disk. We do not fully couple the MC algorithm to the fluid evolution, but use a short evolution of the merger remnant to critically assess errors in our approximate gray two-moment transport scheme. We demonstrate that the analytical closure used by the moment scheme is highly inaccurate in the polar regions, but performs well elsewhere. While the average energy of polar neutrinos is reasonably well captured by the two-moment scheme, estimates for the neutrino energy become less accurate at lower latitudes. The two-moment formalism also overestimates the density of neutrinos in the polar regions by $\sim 50\%$, and underestimates the neutrino pair-annihilation rate at the poles by factors of 2–3. Although the latter is significantly more accurate than one might have expected before this study, our results indicate that predictions for the properties of polar outflows and for the creation of a baryon-free region at the poles are likely to be affected by errors in the two-moment scheme, thus limiting our ability to reliably model kilonovae and gamma-ray bursts.

Figure 1

Difference in the pitch-angle distribution and energy spectrum of neutrinos between two simulations placing the boundary of the MC evolution at different optical depth. We look at neutrinos within a 3 km radius of a point on the polar axis ($z=45\mathrm{km}$), about 5 ms after the beginning of the simulation. Errors are normalized by the expected Poisson error in our reference simulation. Both simulations have near-identical Poisson noise.

Figure 2

Poloidal (top) and equatorial (bottom) slices through the merger remnant. In each figure, solid lines show density contours of ${\log }_{10}(\rho /[1\mathrm{g}/{\mathrm{cm}}^3])=(8,9,10,11,12,13,14)$. Color scales show the fluid temperature (left) and electron fraction (right).

Figure 3

One-dimensional profiles of the neutrino moments in the MC (solid lines) and M1 (dashed lines) code for each species of neutrinos (${\nu }_x$ values are for all heavy-lepton neutrinos combined), 14 ms after merger. We show the energy density (top), the parallel (radial) component of the normalized transport flux (average velocity of the neutrinos, middle), and the parallel component of the Eddington tensor (bottom) along the $z$ axis (left) and $x$ axis (right). M1 values are computed by interpolation, while MC moments are computed using all packets within a ball of radius of $300\mathrm{m}+0.05*R$ centered on the desired point, with $R$ the distance to the center of the remnant. All values are computed in the inertial frame, without projection into an orthonormal frame (hence the few points where $P_{\parallel }>E$).

Figure 4

Same as Fig. 3, but for the average energy of the neutrinos in the inertial frame.

Figure 5

Energy density of neutrinos for all species in the M1 and MC codes, along a line 25° from the polar axis.

Figure 6

Poloidal slice through the merger remnant. The left column shows the diagonal and off-diagonal components of the Eddington tensor $f_{zz}$ and $f_{xz}$ as measured with the MC code. Solid lines are contours of $f_{zz},f_{xz}$ separated by 0.1 (alternating red and black lines). This provides an order-of-magnitude estimate of the statistical noise, with $\mathrm{\Delta }{f}_{\mathrm{stat}}<0.05$. The right column shows the difference between the MC and M1 Eddington tensors. In the massive neutron star, the difference is negligible. In the disk, it is dominated by statistical MC noise. Outside of the remnant, large errors in the Minerbo closure dominate. In particular, $\mathrm{\Delta }{f}_{zz}\sim 0.1–0.3$ in the polar regions. All results are for ${\nu }_e$ neutrinos.

Figure 7

Same as Fig. 6, but for an equatorial slice through the remnant, and showing components $f_{xx},f_{xy}$ of the Eddington tensor. The Minerbo closure is significantly more accurate here than in the polar regions, although some regions of the shocked spiral arms show consistent biases in the Minerbo closure at a level $\mathrm{\Delta }{f}_{ij}\sim 0.1$.

Figure 8

Distribution probability of neutrinos as a function of the azimuthal angle $\varphi$ of their 4-momentum, at different points on the polar axis. We sample all packets within $\mathrm{\Delta }d=0.1z$ of the target point. We have 400–2000 packets within each region and for each type of neutrino.

Figure 9

Same as Fig. 8, but for the distribution probability with respect to $\cos (\theta )$, with $\theta$ the neutrino pitch angle. We show multiple species for the closest point to the remnant, as the angular distribution is more sensitive to the finite size of the emitting region close to the remnant.

Figure 10

Same as Fig. 8, but for the distribution probability with respect to the neutrino energy $\epsilon$. We use 12 energy bins to generate this figure, identical to the bins of the NuLib table.

Figure 11

Left: Energy distribution of the neutrinos leaving the computational domain (measured on the surface of a parallelepiped 6 grid spacings from the subdomain boundary, i.e., approximately a parallelepiped of size $[480\times 480\times 240]\mathrm{km}$ centered on the origin), for all 3 species of neutrinos. In each plot, the dashed vertical line shows the average neutrino energy estimated by the M1 scheme, and the solid vertical line the same quantity estimated by the MC scheme. The solid grey line shows our best fit to the spectrum. Right: Angular distribution of the neutrinos leaving the grid. Here, $\theta$ is the usual spherical-polar coordinate, not the pitch angle of the neutrinos. Grey histograms show the MC results, and red histograms the M1 results. Errors in the Minerbo closure lead to a large overestimate of the neutrino density in the polar regions. In all plots, we integrate the neutrino fluxes over a $50\mu \mathrm{s}$ interval 14 ms after merger.

Figure 12

Vertical slices through the merger remnant. Top: The specific energy deposition rate due to $\nu \overline{\nu }$ pair annihilation using the neutrino moments and average energy predicted by the MC code (right), and the ratio of the geometric factor $\kappa$ obtained from the MC and M1 codes (left). The latter provides an estimate of the error due solely to the use of the Minerbo closure in the computation of Eq. (24). Bottom: The geometric factor $\kappa$ in Eq. (24), for ${\nu }_x{\overline{\nu }}_x$ annihilation (left) and ${\nu }_e{\overline{\nu }}_e$ annihilation (right). Note that as the MC code only stores moments normalized to the energy density $E$ (i.e., $F^i/E$, $P_{ij}/E$), the computation of $Q_{\mathrm{MC}}$ still uses the energy density evolved by the M1 code.

Figure 13

Geometric factor $\kappa$ in Eq. (24), along the polar axis. We show results for electron and heavy-lepton neutrinos, as well as best-fit curves using Eq. (25). The fits ignore the high-density regions inside the remnant neutron star, where pair annihilation is a subdominant process.