the choice made, the difficulty of properly capturing the evolution of
$\tilde{N}$ in the diffusion regime, and thus the diffusion of $Y_{e}$ through
dense regions, is a limitation of existing moment schemes.

Finally, for completeness on the topic of closures, it is worth mentionning
that Shibata al shibata:11  and Cardall et al PhysRevD.87.103004  have
provided closures for the third moment $M^{\alpha\beta\gamma}$. One method
assumes that the expansion of the distribution function in $l^{\alpha}$ is
truncated at second order, in which case the third moment is an explicit
function of the first moment. An alternative is to use the Minerbo closure to
again interpolate between optically thick and free-streaming estimates of
$L^{\alpha\beta\gamma}$. These closures are not needed in existing grey moment
schemes, but would be necessary for energy-dependent schemes.

####  4.2.6 Numerical implementation

In the previous sections, we covered most of the ingredients needed to evolve
moments of $f_{\nu}$: the moment equations in the simulation frame, their
coupling to the fluid equations, and the required analytical closures for
higher-order moments of $f_{\nu}$ (and for the neutrino energy spectrum in
grey schemes). These equations are very similar to the equations of
hydrodynamics in conservative form, i.e. they can be expressed as

$\partial_{t}U+\partial_{i}F^{i}(U)=S(U)$ (129)

with $U$ the vector of evolved variables (e.g.
$U=(\tilde{E},\tilde{F}_{i},\tilde{N})$ in a grey M1 scheme that evolves the
number density), $F^{i}$ the fluxes, and $S$ the local source terms. These
equations can be evolved using the same high-order shock-capturing methods
developed for the equations of fluid dynamics, guaranteeing that conservation
laws are satisfied to round-off accuracy in our evolutions. The main
complications will come from the source terms, which can be extremely large
for radiation transport, thus requiring the right-hand side of this equation
to be treated implicitly. The flux terms, on the other hand, can be treated
explicitly as long as the timestep satisfies the usual Courant condition
$\Delta t\lesssim\alpha_{CFL}\Delta x$, with $\Delta x$ the grid spacing and
$\alpha_{CFL}$ a constant of order unity that depends on the exact numerical
methods used to evolve these equations. In practice, it is thus useful to
consider split implicit-explicit time evolutions

$\partial_{t}U=\left(S_{\rm exp}(U)-\partial_{i}F^{i}(U)\right)+S_{\rm
imp}(U)$ (130)

with $S_{\rm imp}+S_{\rm exp}=S$, and $S_{\rm imp}$ contains all terms that we
choose to treat implicitly (typically, neutrino-matter interactions). A first-
order in time discretization would then be the implicit equation

$U^{n+1}=U^{n}+\Delta t\left(S_{\rm
exp}(U^{n})-\partial_{i}F^{i}(U^{n})\right)+(\Delta t)S_{\rm imp}(U^{n+1}).$
(131)

where upper indices $(n,n+1)$ refer to the beginning/end of the time step. In
a conservative scheme, the spatial discretization requires us to consider
values of the fields at grid points, and halfway between grid points. For
example, in 1D,

$U^{n+1}_{i}=U^{n}_{i}+\Delta t\left(S_{\rm
exp}(U^{n}_{i})-\frac{F^{*}(U^{n}_{i+1/2})-F^{*}(U^{n}_{i-1/2})}{\Delta
x}\right)+(\Delta t)S_{\rm imp}(U^{n+1}_{i})$

with the subscripts referring to grid points / cell centers (integer values)
and half grid points / cell faces (half-integer values). Many possible choices
can be made for the calculation of the numerical fluxes $F^{*}$, which will
offer different orders of convergence and shock-capturing capabilities. A
simple, very dissipative choice would for example be the Lax-Friedrich flux

$F^{*}_{i+1/2}=\frac{1}{2}\left(F_{i}+F_{i+1}\right)-\frac{c}{2}\left(U_{i+1}-U_{i}\right).$
(132)

Modern numerical simulations however often use higher-order, less dissipative
methods. In those advanced methods, values of $F,U$ at $(i,i+1)$ are replaced
with left-biased and right-biased stencils used to interpolate $(F,U)$ on cell
faces MC ; Jiang1996202 ; Borges , while the speed of light $c$ is replaced by
the maximum characteristic speed of the system of equations. The general idea,
however, remains: the numerical fluxes include a first term estimating the
flux at the mid-point, and a second term providing numerical dissipation at
shocks to smooth the solution and avoid instabilities. Alternatively, the
evolution equations can be projected onto their characteristic fields,
allowing for the use of even less dissipative methods (see e.g. Radice:2012cu
for general relativistic fluid dynamics). We will not review here the
extensive literature discussing choices for these numerical fluxes, but note
that the characteristic speeds for the two-moments algorithm can for example
be found in Shibata et al shibata:11 , and details of the calculation of
$F^{*}$ for two-moment algorithms used in merger simulations are available in
shibata:11 ; FoucartM1:2015 ; Radice:2021jtw ; Sun:2022vri . The only
difficulty specific to the radiation transport equations is the treatment of
high-density regions, discussed in more detail in Sec. 4.2.7.

For the treatment of the implicit terms, the most up-to-date methods were
derived by Radice et al Radice:2021jtw . In that work, the Jacobian matrix

$J(U)=\frac{\partial S_{\rm imp}(U)}{\partial U}$ (133)

is calculated explicitly for the Minerbo closure. This allows them to solve
the implicit equations for $U^{n+1}$ using standard iterative methods for the
determination of the roots of a multi-dimensional function $f(U)$, given $f$,
$\partial f/\partial U$, and a reasonable initial guess $U_{g}$ for the
solution (e.g. the value of $U$ at the beginning of a time step, or the
equilibrium value of $U$). We note that in this case,
$U=(\tilde{E},\tilde{F}_{x},\tilde{F}_{y},\tilde{F}_{z}$). Each species of
neutrinos is treated separately, and the evolution of the number equation (if
included) can be performed after the evolution of $U$ (as the implicit source
terms for $U$ are independent of $N$). As $J(U)$ is quite complex, we refer
the reader to Radice:2021jtw  for its exact form. Many older two-moment
schemes however use approximate values of the Jacobian matrix, e.g.
linearizing the problem around $E=F_{i}=0$ and/or fixing some variables
entering the calculation of $J$ to their value for $U=U^{n}$ instead of using
$U=U^{n+1}$, in order to simplify the calculation of $J$ at the cost of some
accuracy in the implicit solve.

Finally, we note that most moment algorithms use a split operator method to
couple the fluid and neutrino evolution. Assuming that $U_{\rm fl}$ and
$U_{\rm rad}$ are the vectors of evolved variables for the fluid and neutrinos
respectively, the coupled problem is evolved using equations of the form

$\displaystyle U^{n+1,*}_{\rm fl}$ $\displaystyle=$ $\displaystyle f_{\rm
fl}(U^{n}_{\rm fl},\Delta t)$ (134) $\displaystyle U^{n+1}_{\rm rad}$
$\displaystyle=$ $\displaystyle f_{\rm rad}(U^{n+1,*}_{\rm fl},U^{n}_{\rm
rad},U^{n+1}_{\rm rad},\Delta t)$ (135) $\displaystyle U^{n+1}_{\rm fl}$
$\displaystyle=$ $\displaystyle U^{n+1,*}_{\rm fl}+S_{\rm fl/rad}\Delta t.$
(136)

The first line represents the usual evolution of the fluid equations without
coupling to the neutrinos. The second line represents the mixed implicit-
explicit evolution of the radiation equations on a given fluid background.
Finally, the third line represents backreaction of the radiation onto the
fluid. It is also possible to improve on this scheme by using a guess
$U^{n+1,g}_{\rm fl}$ for the fluid variables on the second line. This is
particularly useful in dense regions, where we might use the expected state of
the fluid once neutrinos and matter equilibrate (see e.g. Foucart:2016rxm ).
Using such a guess is sometimes necessary to avoid numerical instabilities in
regions where the coupling between neutrinos and matter leads to stiff source
terms in the evolution of the fluid equations (i.e. in the $S_{\rm fl/rad}$
term), in addition to the stiff source terms that nearly always exist in the
evolution of the neutrino moments.

A more self-consistent method would be of course to use an implicit-explicit
solver to evolve jointly the fluid equations and the moment equations for all
species of neutrinos – but this would make the implicit part of the solver
significantly more costly. Indeed, the standard scheme solves, for each
neutrino species, a system of 4 coupled non-linear implicit equations for
$(\tilde{E},\tilde{F}_{i})$ and a single linear implicit equation for
$\tilde{N}$. A full implicit solve of the fluid and radiation equations, on
the other hand, would require solving a system of $6+5N_{\nu}$ non-linear
implicit equations for the $N_{\nu}$ neutrino species and the fluid (assuming
$6$ fluid variables and 5 components of the moments). As long as the simpler
method does not lead to numerical instabilities, its significantly lower
computational cost makes it much more appealing.

####  4.2.7 Diffusion regime

In regions where $(\kappa_{a}+\kappa_{s})L\gg 1$, with $L$ a typical
lengthscale of our system, we expect the neutrino momentum density in the
fluid frame to be $\propto-(\kappa_{a}+\kappa_{s})^{-1}\nabla J$. This leads
to an evolution equation for the energy density in the fluid frame

$\partial_{\hat{t}}J-\delta^{\hat{i}\hat{j}}\partial_{\hat{i}}\left(\frac{1}{3(\kappa_{a}+\kappa_{s})}\partial_{\hat{j}}J\right)=\kappa_{a}(J^{\rm
eq}-J).$ (137)

i.e. a diffusion equation with diffusion coefficient
$D=\left[3\left(\kappa_{a}+\kappa_{s}\right)\right]^{-1}$. However, when
evolving the two-moment equations with shock-capturing methods, the
dissipative terms in the numerical fluxes modify this optically thick
solution. For example, with a low-order numerical flux like Eq. 132, and
ignoring spatial variations in the opacities, the discretized equation becomes
Audit:2002hr

$\partial_{\hat{t}}J-\left(\frac{1}{3(\kappa_{a}+\kappa_{s})}+\frac{c\Delta
x}{2}\right)\delta^{\hat{i}\hat{j}}\partial^{2}_{\hat{i}\hat{j}}J=\kappa_{a}(J^{\rm
eq}-J).$ (138)

We clearly see that in any region where $(\kappa_{a}+\kappa_{s})\Delta
x\gtrsim 1$, numerical diffusion will be larger than physical diffusion.
Higher-order fluxes will be more forgiving, but even high-order shock
capturing methods behave similarly to their low-order counterparts near shocks
or in underresolved regions. As $J$ may vary rapidly in the hot, dense regions
of a merger remnant, the diffusion rate of neutrinos through that remnant
could plausibly be impacted by numerical dissipation.

To avoid this issue, two methods have been proposed so far. Audit et al
Audit:2002hr  suggest to effectively use a one-moment scheme in regions where
$(\kappa_{a}+\kappa_{s})\Delta x\gtrsim 1$, i.e. to replace the numerical flux
in the evolution of $\tilde{E}$ by its value assuming that

$H^{\alpha}=-\frac{h^{\alpha\beta}}{3(\kappa_{a}+\kappa_{s})}\nabla_{\beta}J\mathchar
24635\relax\;\,\,L^{\alpha\beta}=\frac{h^{\alpha\beta}}{3}J.$ (139)

This leads to neutrinos being advected with the fluid, with an additional slow
diffusion provided by $H^{\alpha}$. A more detailed discussion of this method
in the merger context can be found in FoucartM1:2015 . As pointed out in
Radice:2021jtw , this method does however transform our evolution equations
into a single diffusion equation, which is known to be acausal and potentially
unstable in general relativity. Radice et al Radice:2021jtw  thus propose an
alternative, numerically simpler method: they use a high-order flux without
numerical dissipation (i.e. using finite difference methods) in regions where
$(\kappa_{a}+\kappa_{s})\Delta x\gtrsim 1$. In both algorithms, one
transitions smoothly between the standard shock-capturing fluxes and the
modified fluxes around $(\kappa_{a}+\kappa_{s})\Delta x\sim 1$. For example,
Radice et al Radice:2021jtw  use

$F^{*}=(1-a)F^{\rm HO}+aF^{\rm LO}$ (140)

with $F^{\rm LO}$ given by Eq. 132 (replacing $c$ by the value of the speed of
light in the simulation coordinates) and $F^{\rm HO}$ using the simple second-
order accurate prescription

$F^{\rm HO}_{i+1/2}=\frac{F_{i}+F_{i+1}}{2}.$ (141)

The transition coefficient is

$a=\min{\left(1,\frac{1}{\Delta x(\kappa_{a}+\kappa_{s})}\right)}$ (142)

and the opacities are estimated using the average values of neighboring grid
points.

Figure 4: Advection of trapped neutrinos in a rapidly moving fluid ($v=0.5c$)
with high scattering opacity, and no absorption or emission. The different
simulaitons use (i) the latest implementation of the Jacobian matrix and high-
order numerical fluxes from Radice:2021jtw  (blue), (ii) approximate sources
and a diffusion equation (ZelmaniM1 code, black), or (iii) only first order in
$v/c$ terms (green). Figure reproduced from Radice al (2022) Radice:2021jtw .

Radice et al Radice:2021jtw  show that the approximate linearization of the
source terms used in many existing two-moment schemes combined with their
reliance on the solution of a diffusion equation leads to inaccuracies in the
evolution of the transport equations at the very least in regions where
$\kappa_{s}\Delta x\gg 1$, $v/c\sim 1$, and $\kappa_{a}\Delta x\ll 1$ (i.e.
where advection and diffusion are important, and the evolution is not driven
to the equilibrium value of the moments by a high absorption rate). This
regime is likely to be relevant in the outer regions of a rotating post-merger
neutron star remnant for heavy-lepton neutrinos. Additionally, these
approximations may lead to other sources of errors in yet untested regimes. On
the other hand, the use of non-dissipative fluxes could, in theory at least,
lead to issues at shocks in dense region – though this does not seem to have
been observed in simulations using these newer methods so far.

####  4.2.8 Pair annihilation

Pair processes pose a particular challenge for moment-based transport
algorithms – and, in fact, for many other types of radiation transport schemes
– because they naturally lead to non-linear source terms in the transport
equations that couple the distribution function of neutrinos and antineutrinos
(see Sec. 3.2.2). On the other hand, fully ignoring pair production and
annihilation is not an option for muon and tau neutrinos. In grey moment
schemes, we thus typically calculate an emission rate of
$\nu_{\mu}\bar{\nu}_{\mu}$ pairs and $\nu_{\tau}\bar{\nu}_{\tau}$ pairs with
blocking factors computed assuming an equilibrium density of neutrinos, and
model the inverse reaction (pair annihilation) through an absorption opacity
calculated using Kirchhoff’s law.

For a typical merger profile, we expect heavy lepton neutrinos to experience
higher scattering opacities than absorption opacities. Accordingly, those
neutrinos remain in thermal equilibrium with the fluid while in regions where
their absorption optical depth $\tau_{a}<1$, but their scattering optical
depth $\tau_{s}\gg 1$. In those regions, our assumptions for pair production
and annihilation may be inaccurate, but they should nonetheless practically
capture the physics of neutrino diffusion quite well. Closer to the surface of
the remnant, however, heavy lepton neutrinos slowly diffuse out while out of
thermal equilibrium with the fluid (typically, neutrinos will be hotter than
expected). In these regions, our assumptions likely underestimate the rate of
pair annihilation. The impact of this approximation has not been studied so
far.

For electron type neutrinos, $\nu_{e}$ and $\bar{\nu}_{e}$ decouple from the
fluid in different regions, and the above assumptions would be even more
problematic. However, in optically thick regions, pair production and pair
annihilation are expected to be subdominant – except maybe in the hottest
regions, where neutrinos will be in equilibrium with the fluid regardless of
the reactions included in a calculation. Accordingly, ignoring pair production
for electron-type neutrinos has often been considered safer than including it
in a very approximate manner.

This leaves us with one important issue, however: $\nu\bar{\nu}$ annihilation
is expected to be an important process for energy deposition in low-density
regions around the rotation axis of a post-merger remnant, and may also
deposit energy and momentum in low-density matter outflows elsewhere. From
equation Eq. 29, we get that the appropriate energy-integrated source term for
the moment equations is, for annihilation into $e^{+}e^{-}$ pairs,

$\displaystyle S^{\alpha}_{(\nu)}$ $\displaystyle=$
$\displaystyle-\int\int\frac{d^{3}p_{(\nu)}}{\sqrt{-g}p^{t}_{(\nu)}h^{3}}\frac{d^{3}p_{(\bar{\nu})}}{\sqrt{-g}p^{t}_{(\bar{\nu})}h^{3}}f_{(\nu)}f_{(\bar{\nu})}(p^{\alpha}_{(\nu)})\frac{DG_{F}^{2}}{3\pi}\left(-p_{(\nu)}^{\beta}p_{\beta(\bar{\nu})}c^{2}\right)^{2}$
(143) $\displaystyle=$
$\displaystyle-\frac{c^{4}DG_{F}^{2}}{3\pi}M_{\beta\gamma,(\bar{\nu})}\int\frac{d^{3}p_{(\nu)}}{\sqrt{-g}p^{t}_{(\nu)}h^{3}}f_{(\nu)}p^{\alpha}_{(\nu)}p^{\beta}_{(\nu)}p^{\gamma}_{(\nu)}.$

The integral over $p_{(\nu)}$, however, does not match the moments used in our
evolution equations (it includes an extra power of the neutrino energy). A
similar term should also be included for the antineutrinos. Fujibayashi al
Fujibayashi:2017xsz  propose the approximation

$S^{\alpha}_{(\nu)}=-\langle\epsilon_{(\nu)}\rangle
u^{\alpha}\frac{c^{4}DG_{F}^{2}}{3\pi}M_{\beta\gamma,(\bar{\nu})}M^{\beta\gamma}_{(\nu)}$
(144)

with $\langle\epsilon_{(\nu)}\rangle$ some appropriate average of the energy
of annihilated neutrinos. Given the scaling in the source term, the choice

$\langle\epsilon_{(\nu)}\rangle=\frac{F_{4}(\eta_{(\nu)})}{F_{3}(\eta_{(\nu)})}k_{B}T_{(\nu)}$
(145)

would be reasonable for a quasi-thermal spectrum at an estimated temperature
$T_{(\nu)}$.

Figure 5: Density, velocity, and heating rate from neutrino absorption
$20\,{\rm ms}$ after a neutron star merger in a simulation using pair
annihilation (Left), and a similar simulation ignoring that effect (Right).
Note the significant change in the velocity of the outflows, which persists up
to the end of these simulations ($300\,{\rm ms}$ post-merger). Figure
reproduced from Fujibayashi et al (2017) Fujibayashi:2017xsz .

This remains a significant approximation. Indeed, by using $S^{\alpha}\propto
u^{\alpha}$, the above formula ignores momentum deposition in the fluid frame,
even though neutrinos leaving the remnant definitely have a preferred
direction of propagation. Additionally, the source term depends on the
neutrino pressure tensor, which is estimated through the chosen analytical
closure. Polar regions where pair annihilations are expected to be important
(see Fig. 5) are also where that closure creates the largest errors. Despite
these limitations, comparisons with Monte-Carlo transport results
Foucart:2018gis  indicate that $S^{\alpha}$ is at least order-of-magnitude
accurate when using this approximation – which is about the best that one can
hope for in the context of a grey moment scheme.

####  4.2.9 Discussion

The main algorithms for general relativistic radiation transport implemented
so far in merger simulations, moment schemes have allowed us to greatly
increase our understanding of the role of neutrinos in mergers. In particular,
moment simulations showed that neutrino-matter interactions have a fairly
dramatic impact on the composition of hot matter outflows, with important
consequences for r-process nucleosynthesis Wanajo:2014 . The cost and
complexity of moment schemes is however significantly larger than those of the
simplest leakage schemes: the moment equations without source terms are
comparable in complexity to the evolution of the fluid equations, and if large
regions of the computational domain require the use of implicit timestepping
(as in most merger simulations that do not result in the formation of a black
hole), radiation transport can easily become the main computational cost of a
simulation. Mixed leakage-moment schemes can do away with that last issue, but
at the cost of a more approximate treatment of neutrino diffusion in dense
regions.

The most up-to-date moment schemes can likely capture very well the diffusion
of single-energy neutrinos in dense regions, yet as moment schemes used in
general relativistic merger simulations do not keep detailed information about
the neutrino energy spectrum, and the diffusion of neutrinos through a dense
medium is strongly energy dependent, their ability to properly capture the
diffusion of both energy and lepton number consistently is more difficult to
assess. We will see that simulations with different methods provide wildly
different answers for the emission of heavy-lepton neutrinos Foucart:2020qjb
– which may be in part because those neutrinos spend a significant amount of
time in regions where they are out of thermal equilibrium with the fluid, yet
still experiencing high scattering opacities. Approximate treatment of the
source terms and the use of a solution to the acausal diffusion equations in
many moment schemes may also be a source of uncertainty in that regime
Radice:2021jtw . In the semi-transparent and optically thin regimes, moment
algorithms offer a reasonable approximation to the qualitative evolution of
the composition and temperature of the outflows. However, simulations have
shown that the assumed energy spectrum of neutrinos in those regions has a
non-negligible impact on the final composition of the ejecta Foucart:2016rxm .
Energy-dependent moment schemes could do away with most of those issues, but
would be computationally expensive, even if difficulties related to the choice
of reference frame in which the neutrino energy is discretized were to be
solved.

In optically thin regimes, the closure chosen in existing simulations is also
a source of concern. Factor of $\sim 2$ errors in the energy density of
neutrinos in the polar region observed in comparisons with Monte-Carlo
simulations Foucart:2018gis  are very likely due to this approximate closure.
Combined with the lack of information about the exact distribution function of
neutrinos, this limits the ability of moment schemes to take into account the
role of pair annihilation in the development of a baryon-free region and
relativistic jet along the remnant’s rotation axis – even though existing
simulations that are expected to capture the importance of pair annihilation
at least at the qualitative level indicate that pair annihilation may
significantly impact the properties of polar regions Fujibayashi:2017xsz ;
Foucart:2018gis  (see also Fig. 5).

Despite these limitations, comparisons of moment schemes with simulations
using Monte-Carlo transport offer some reassurances regarding the validity of
the moment results Foucart:2020qjb , with disagreements at the $10\%-20\%$
level in most global quantities, including the $\nu_{e}$ and $\bar{\nu}_{e}$
luminosities, the average energy of neutrinos, and the mass and composition of
the outflows. Given the difficulty of estimating errors directly in moment
simulations, further studies of the uncertainties associated with the many
different moment schemes currenlty used in the merger community are certainly
required, but the overall understanding of the most important neutrino
processes in merger simulations gathered from these simulations is likely
accurate (except for the potential role of neutrino oscillations).

###  4.3 Monte-Carlo radiation transport

Grey moment schemes are at this point the most commonly used algorithms to
approximately evolve Boltzmann’s equation in merger simulations. We have seen
however that one of their limitations is the strong assumptions that need to
be made about the energy spectrum of the neutrinos, and the impact of these
assumptions on neutrino-matter interaction rates. The use of an approximate
closure for the pressure also introduces an important approximation in
simulations. As a result, even in the limit of infinite resolution, a two-
moment scheme does not converge to the correct solution to Boltzmann’s
equations.

Direct discretizations of Boltzmann’s equation in both position space and
momentum space using finite difference, finite volume, or spectral methods
have not yet been used in general relativistic merger simulations, and are
likely beyond our current computational capabilities. While we wait for these
methods to become practical in the merger context, however, an interesting
alternative is the use of Monte-Carlo radiation transport. Monte-Carlo methods
are particularly efficient for low-cost evolutions of high-dimensional, highly
inhomogeneous problems. Their scaling with increasing computational resources
is significantly worse than the other methods mentioned so far, and they will
thus nearly certainly become less efficient than those other algorithms at
some point in the future – but at the moment, they are the only algorithms
going beyond the moment formalism to have been implemented in general
relativistic merger Foucart:2020qjb  and post-merger Miller:2019dpt
simulations.

In this section, we will first review the formalism of Monte-Carlo radiation
transport in general relativistic simulations (Sec. 4.3.1). We will then
discuss technical issues related to the coupling of the neutrinos to the fluid
(Sec. 4.3.2), the treatment of optically thick regions (Sec. 4.3.3), and the
inclusion of non-linear source terms such as pair processes (Sec 4.3.4). As
the use of Monte-Carlo methods in neutron star merger simulations is
significantly less mature than the use of leakage or moment schemes, and the
exact methods are likely to change rapidly in the next few years, we will keep
our discussion more general than in the previous section, focusing on the
important components of an algorithm and the main difficulties that have been
encoutered so far.

As in the previous section, we assume $h=c=1$ unless noted otherwise.

####  4.3.1 Formalism

The general idea behind Monte-Carlo methods for radiation transport is to
discretize the distribution function of neutrinos using Monte-Carlo packets
(or superparticles) each representing a large number of neutrinos, i.e.

$f_{(\nu)}=\sum_{k\in
P}N_{k}\delta^{3}\left(x^{i}-x^{i}_{k}(t)\right)\delta^{3}\left(p_{i}-p_{i}^{k}(t)\right)$
(146)

with $P$ the ensemble of all packets, $x^{i}_{k}(t)$ the position of packet
$k$ as a function of time, and $p_{i}^{k}(t)$ the spatial components of the
momentum one-form of the neutrinos in that packet. A Monte-Carlo algorithm
needs prescriptions to create packets (emission), propagate them on the grid,
destroy them (absorption), and handle non-destructive interactions with the
fluid or other neutrinos (e.g. scattering events). The propagation of
neutrinos can simply be performed by moving packets along geodesics. All other
processes have to be treated probabilistically, in such a way that the
ensemble of packets remains a good sampling of the true distribution function
$f_{(\nu)}$.

The finite number of packets in a Monte-Carlo simulation will inevitably lead
to sampling noise in $f_{(\nu)}$, and in any variable derived from
$f_{(\nu)}$. We typically expect relative errors in quantities obtained by
summing over $N$ packets to be $\sim N^{-1/2}$. If one wanted $\sim 1\%$
errors for the energy density of neutrinos within each cell of a computational
simulation at any given time, one would thus need $\sim 10^{4}$ packets per
cell… and scaling to smaller error bars is extremely expensive. A saving grace
for Monte-Carlo simulations in the merger context, however, is that neutrino-
matter interactions are generally not dynamically important, and tend to
change the properties of the fluid over relatively long time scales (with a
few exceptions). In practice, this means that we will be able to perform
reasonably accurate Monte-Carlo simulations by relying on averaging neutrino-
matter interactions over time scales significantly longer than our numerical
time step in the vast majority of the computational domain. Practically, many
regions can in fact be evolved stably and reliably with less than a packet per
grid cell. This should be contrasted with astrophysical systems where
radiation is dynamically important (e.g. radiation-dominated accretion disks).
In those systems, sampling noise in e.g. the radiation pressure can be
problematic. The price to pay for the use of a small number of packets,
however, is that we lose the ability to get instantaneous estimates of the
neutrino distribution function.

Emission: The creation of Monte-Carlo packets is theoretically simple to
handle, as long as we know the emission rate of neutrinos as a function of
momentum, and are provided with a prescription to choose the number of
neutrinos represented by each packet. For example, if we assume that all
packets within a given grid cell of coordinate volume $\Delta V$ represent a
total neutrino energy $E_{p}$, and that the energy-integrated emissivity is
$\eta_{\rm tot}$, we create

$N=\frac{\eta_{\rm tot}\sqrt{-g}\Delta V\Delta t}{E_{p}}$ (147)

packets over a time step $\Delta t$ in that cell. Non-integer values of $N$
can be treated statistically (e.g. $N=0.2$ implies a $20\%$ chance of creating
a single packet). The direction of propagation of the packets is drawn from an
isotropic distribution in the fluid frame and from the energy distribution
$f(\nu)=\eta(\nu)/\eta_{\rm tot}$. The initial time and position of a packet
can be drawn from an homogeneous distribution within the 4-volume $(\Delta
V\Delta t)$, or possibly from a spatial distribution more adapted to the
problem (e.g. in axisymmetry, it is beneficial to assume that the volume
element is $\propto rdr$, making it more likely that packets are created in
the outer region of a cell than in its inner region).

In many simulations, we have at our disposal tabulated values of $\eta(\nu)$
within specific energy bins, rather than a known continuous function
$\eta(\nu)$. In that case, a common strategy is to calculate within each bin
the number of emitted packets $N_{b}$. Packets can then be created in the
fluid frame with the energy of the center of their energy bin. Giving all
packets the energy of the center of the bin guarantees that if we have
tabulated value of $\kappa_{a}$ at the center of the same energy bins, the
equilibrium energy density of neutrinos can match its expeceted physical value
to numerical roundoff richers:15 .

We see that the emission step is thus conceptually simple. However, it
involves a choice crucial to the efficiency of Monte-Carlo methods: that of
the energy of a packet $E_{p}$. That energy can vary from cell to cell, over
time, and from energy bin to energy bin, and is one of the main tool available
to distribute computational resources efficiently through a simulation. In
neutron star mergers, for example, it is convenient to choose $E_{p}$ in
optically thick regions in order to get a desired number of packets per cell
(and thus a desired statistical noise), while $E_{p}$ in optically thin
regions can be chosen to obtain a desired number of packets over the entire
simulation (setting the overall cost of the simulation). Such an algorithm is
detailed in Foucart:2021mcb . We note however that many other choices of
$E_{p}$ are possible, and that the optimal choice is problem dependent.

Propagation: Neglecting the finite mass of neutrinos, we expect packets to
propagate along null geodesics in between interactions with the fluid and
other neutrinos. Hughes et al Hughes1994  showed that to do so, it is
convenient to evolve the position $x^{i}$ and spatial components of the
momentum one-form $p_{i}$. The corresponding geodesic equations were initially
developed for ray-tracing and photon transport, but are also appropriate for
neutrinos:

$\displaystyle\frac{dx^{i}}{dt}$ $\displaystyle=$
$\displaystyle\gamma^{ij}\frac{p_{j}}{p^{t}}-\beta^{i}$ (148)
$\displaystyle\frac{dp_{i}}{dt}$ $\displaystyle=$ $\displaystyle-\alpha
p^{t}\partial_{i}\alpha+p_{j}\partial_{i}\beta^{j}-\frac{1}{2}p_{j}p_{k}\partial_{i}\gamma^{jk}$
(149) $\displaystyle p^{t}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{\gamma^{ij}p_{i}p_{j}}}{\alpha}.$ (150)

The first two lines are evolution equations for $x^{i}$, $p_{i}$, while the
third comes from the constraint that $p^{\mu}p_{\mu}=0$. Numerically, the main
choices to make here are the time stepping algorithm, and a method for the
computation of the metric and its derivatives at the location of the packet.
We note that as packets are either slowly diffusing through the system (in
which case we do not directly use these equations; see below) or rapidly free-
streaming out of the computational domain, even low-order methods perform well
enough that the propagation step is not a leading source of error in
simulations (see e.g. Foucart:2021mcb ).

Additional interactions: Finally, we have to allow packets to interact with
the fluid and/or other neutrinos. If we imagine that we have a set of
processes with known mean free path $\lambda_{i}$, or opacity
$\kappa_{i}=\lambda_{i}^{-1}$, then for each process we can randomly draw the
time to the next interaction from a Poisson distribution, e.g.

$\Delta\tau_{i}=-\lambda_{i}\ln r_{i}=-\kappa^{-1}_{i}\ln r_{i}$ (151)

with $r_{i}$ a random number drawn from an homogeneous distribution in $[0,1)$
(the repeated index $i$ does not imply summation here). We note that
$\Delta\tau$ is in the reference frame in which we provide $\kappa_{i}$. In
simulations, this is usually the fluid frame. Transforming to the coordinate
time, we get Ryan2015

$\Delta t_{i}=-\kappa^{-1}_{i}\frac{p^{t}}{\nu}\ln r_{i}.$ (152)

The time to the next event is then $\Delta t_{\rm next}=\min{(\Delta
t_{i},\Delta t_{\rm step})}$ with $\Delta t_{\rm step}$ the time step, and
with the minimum taken over all possible processes. If the minimum is $\Delta
t_{\rm step}$, a packet is simply propagated for that time. Otherwise,
whichever process provided the minimum $\Delta t$ occurs after that time
interval. Existing simulations have considered an absorption opacity
$\kappa_{a}$ and an isotropic elastic scattering opacity $\kappa_{s}$, as in
moment schemes, but this method can be generalized to a larger number of
interactions. Absorption simply results in the packet being removed from the
simulation, while isotropic elastic scattering result in redrawing the
direction of propagation of the packet from an isotropic distribution in the
fluid frame, under the constraint that the fluid frame energy is conserved.

As for the propagation of packets, an important step here is how to
interpolate the opacities to the position of the packets and, if opacities are
tabulated, in energy space. Additionally, our estimate for $\Delta t_{i}$ is
only valid for a constant $\kappa_{i}$. If opacities are changing rapidly,
taking too large of a time step can introduce significant errors. The existing
simulation of a neutron star merger with Monte-Carlo transport used constant
$\kappa_{i}$ within each grid cell Foucart:2020qjb , which is numerically
convenient but may need to be improved upon in the future. In addition to
being low-order, this choice leads to an algorithm that is sensitive to the
time at which we determine a packet leaves a cell, and thus recompute the
opacities. As discussed in Foucart:2021mcb , recovering accurate diffusion
rates through high-opacity regions is then only possible if the packets take
fairly small time steps when close to a cell boundary. Higher-order methods
would be desirable, but would complicate calculations of the times-to-
interaction $\Delta t_{i}$ and, to be consistent, would also require the use
of higher-order estimates of the emissivities, and thus inhomogeneous
production of packets within a single grid cell.

Potential Issues: We see that the main steps of a Monte-Carlo algorithm are
fairly straightforward, at least when using low-order methods. They also
naturally follow the expected behavior of individual neutrinos. Unfortunately,
this simple algorithm runs into a number of important roadblocks in practice.
The simplest one to solve is that the evolution still needs to be coupled to
the fluid. Ideally, this would be done in a way that satisfies conservation
laws and avoids unnecessary shot noise in the fluid evolution. We will see
that in practice, there is a trade-off between these two objectives (Sec.
4.3.2). A more significant issue is that the algorithm as proposed would be
extremely inefficient and potentially unstable in optically thick regions.
This is because we then expect rapid creation, annihilation, and scattering of
individual packets, requiring expensive calculations and potentially creating
stiff source terms in the evolution of the fluid equations. This is despite
the fact that the outcome of the evolution of neutrinos in that regime is
known: they simply get into statistical equilibrium with the fluid (if
$\kappa_{a}\Delta t_{\rm step}\gg 1$), and/or slowly diffuse through the fluid
(if $\kappa_{s}\Delta t_{\rm step}\gg 1$). We discuss how one may deal with
these issues in Sec. 4.3.3. Finally, neutrino packets are very inhomogeneously
distributed through the computational grid. This is a desirable effect for
optimal use of our computational resources, but it means that we cannot simply
divide our computational domain into regions with roughly the same number of
grid cells, and then distribute those regions onto different processors.
Parallelization is not a fundamental roadblock to the use of Monte-Carlo
transport, as packets can be distributed to processors in such a way that the
number of packets on each processor is well balanced. Doing so does however
require significant reorganization in merger codes that were not built with
radiation transport in mind, and thus tend to assume that the computational
cost of evolving a given grid cell is roughly the same regardless of the
chosen cell.

Moments of $f_{\nu}$: Before getting into these issues, a useful additional
piece of formalism to discuss is the computation of moments of $f_{(\nu)}$ in
a Monte-Carlo code. Moments can be useful to compute neutrino-matter
interactions, and of course to compare Monte-Carlo results with moment
transport algorithms. The Dirac $\delta$ in the definition of $f_{(\nu)}$ can
be practically problematic when considering moments of $f_{(\nu)}$, and it is
thus often preferable to consider the average moments within a volume $V$
(e.g. a grid cell), defined as

$\bar{M}^{\alpha_{1}...\alpha_{n}}(t,x^{i})=\int_{V}\frac{dV\epsilon^{3}}{V}\int
d\Omega
f_{\nu}(t,x^{i},\epsilon,\Omega)(\hat{t}^{\alpha_{1}}+l^{\alpha_{1}})...(\hat{t}^{\alpha_{n}}+l^{\alpha_{n}}).$
(153)

For the second moment, this is

$\bar{M}^{\alpha\beta}=\sum_{k\in
V}N_{k}\frac{p^{\alpha}_{k}p^{\beta}_{k}}{\sqrt{-g}Vp^{t}_{k}}$ (154)

and similary for the number flux

$\bar{N}^{\alpha}=\sum_{k\in
V}N_{k}\frac{p^{\alpha}_{k}}{\sqrt{-g}Vp^{t}_{k}}.$ (155)

In both cases, the sum is over all packets within the volume $V$. From these
expressions, calculating moments only requires tensor projections. For
example, the energy density in the fluid frame is

$\bar{J}=\sum_{k\in V}N_{k}\frac{\nu_{k}^{2}}{\sqrt{-g}Vp^{t}_{k}}$ (156)

with $\nu_{k}=-p^{\mu}_{k}u_{\mu}$. We note that while this expression will
turn out to be convenient for numerical simulations, its theoretical
interpretation is a little convoluted, as it represents a moment in the fluid
frame, averaged over a volume element in the simulation frame at a constant
simulation time.

####  4.3.2 Coupling to the fluid

In the moment formalism, the back-reaction of the neutrinos onto the fluid was
relatively easy to compute, because the moment equations are directly
equivalent to the equations of conservation of energy, momentum, and, if
evolving the number density, lepton number. A coupling scheme that explicitly
conserves these quantities can also be designed for Monte-Carlo methods, with
a few additional calculations; but we’ll see that it comes with significant
pitfalls.

To get a conservative scheme, we can keep track of changes in the momentum of
neutrinos due to emission, absorption, and scattering. If the $N_{k}$
neutrinos in a packet undergo a change of 4-momentum $\Delta
p_{k}^{\mu}=p^{\mu}_{\rm after}-p^{\mu}_{\rm before}$, then the fluid should
undergo a change of 4-momentum $-N_{k}\Delta p_{k}^{\mu}$ as a reaction. While
this change is generally instantaneous, we can distribute it over a 4-volume
element $\Delta V\Delta t$ to get the 4-force density Ryan2015

$G^{\mu}=-\sum_{\rm events}\frac{N_{k}\Delta p_{k}^{\mu}}{\sqrt{-g}\Delta
V\Delta t},$ (157)

with the sum being over all events changing a packet’s 4-momentum within the
given 4-volume. We note that a single packet may be subject to no event, or to
many events, and that a packet may interact with the fluid in different volume
elements over the course of a time step; thus, properly computing that sum
requires careful bookkeeping. The fluid equations become

$\nabla_{\mu}T_{\rm fl}^{\mu\nu}=G^{\nu}.$ (158)

When using operator splitting, one might instead want to apply the total
change of momentum within the 4-volume to the evolved fluids variables as a
postprocessing step, after evolving the fluid and neutrinos. This requires a
time integration of this equation. The relevant changes for the evolved fluid
variables within a volume $V$ after a time step $\Delta t$ are

$\Delta(\tilde{T}^{\mu\nu}n_{\nu})=\sum_{\rm events}\frac{N_{k}\Delta
p_{k}^{\mu}}{\Delta V}.$ (159)

The most commonly evolved fluid variables in merger simulations are the
internal energy density
$\tilde{\tau}=\tilde{T}^{\mu\nu}n_{\mu}n_{\nu}-\rho_{*}$ and momentum density
$\tilde{S}_{i}=-\tilde{T}^{\mu}_{i}n_{\mu}$. For these variables, we get

$\displaystyle\Delta\tilde{\tau}$ $\displaystyle=$ $\displaystyle\sum_{\rm
events}\frac{N_{k}\Delta p_{k}^{\mu}n_{\mu}}{\Delta V}$ (160)
$\displaystyle\Delta\tilde{S}_{i}$ $\displaystyle=$ $\displaystyle-\sum_{\rm
events}\frac{N_{k}\Delta p^{k}_{i}}{\Delta V}.$ (161)

Similarly, for the evolution of the electron fraction,

$\Delta(\rho_{*}Y_{e})=-m_{b}\sum_{\rm events}\frac{N_{k}s_{k}}{\Delta V}$
(162)

with $s_{k}=0$ for muon and tau neutrinos and for any scattering event,
$s_{k}=1$ for emission of $\nu_{e}$ and absorption of $\bar{\nu}_{e}$, and
$s_{k}=-1$ for emission of $\bar{\nu}_{e}$ and absorption of $\nu_{e}$. The
mass $m_{b}$ is again the reference baryon mass entering the definition of the
rest mass density. This method has the advantage of imposing exact
conservation laws: whatever the neutrinos gain, the fluid looses, and vice-
versa. Its disadvantage is to be fairly sensitive to shot noise. If a packet
is emitted, absorbed, or scattered in a very low density regions, that event
would lead to extremely large changes in the temperature, momentum and
composition of the fluid. These changes may even lead to unphysical values of
the fluid variables after coupling to the neutrinos.

One alternative is to give up on exact conservation laws, counting instead on
conservation laws being statistically satisfied over many interactions. In
such an algorithm, described for example in Foucart:2021mcb , one may write
the source terms for the fluid equations as in the moment equations, i.e. for
neutrinos of a given energy and considering only emission, absorption, and
elastic scattering,

$G^{\mu}=-(\eta-\kappa_{a}J)u^{\mu}+(\kappa_{a}+\kappa_{s})H^{\mu}.$ (163)

As when using exact conservation laws, we then need to integrate this source
term over a time step $\Delta t$ and a small volume $\Delta V$, as well as
over the entire neutrino energy spectrum. We then get

$\int d\nu\int dV\int dtG^{\mu}=-\eta_{\rm tot}\Delta V\Delta
tu^{\mu}+\sum_{k,j}\frac{\Delta
t_{k,j}N_{k}\nu_{k}}{\sqrt{-g}p^{t}_{k}}\left(\kappa_{a,k}p^{\mu}+\kappa_{s,k}p_{k}^{\alpha}h^{\mu}_{\alpha}\right).$

In this expression, the sum is taken over all packets $k$, and all time
intervals $\Delta t_{k,j}$ during which packet $k$ is propagating along a
geodesic while inside of the volume $\Delta V$ and time interval $\Delta t$
(due to scattering events, a single packet may appear multiple times in this
sum). Converting to the change per unit volume in the evolved variables over a
time $\Delta t$, we get

$\displaystyle\Delta\tilde{\tau}$ $\displaystyle=$
$\displaystyle-\tilde{\eta}_{\rm tot}\alpha W\Delta t+\sum_{k,j}\frac{\Delta
t_{k,j}N_{k}\nu_{k}}{\Delta
Vp^{t}_{k}}\left((\kappa_{a,k}+\kappa_{s,k})\epsilon_{k}-\kappa_{s,k}W\nu\right)$
(164) $\displaystyle\Delta\tilde{S}_{i}$ $\displaystyle=$
$\displaystyle-\tilde{\eta}_{\rm tot}\alpha u_{i}\Delta
t+\sum_{k,j}\frac{\Delta t_{k,j}N_{k}\nu_{k}}{\Delta
Vp^{t}_{k}}\left((\kappa_{a,k}+\kappa_{s,k})p_{i}^{k}-\kappa_{s,k}\nu
u_{i}\right).$ (165)

For the electron fraction, we get instead

$\Delta(\rho_{*}Y_{e})=m_{b}\alpha\left(\tilde{\eta}^{(\bar{\nu}_{e})}_{N,\rm
tot}-\tilde{\eta}^{(\nu_{e})}_{N,\rm
tot}\right)+m_{b}\sum_{k,j}s_{k}\kappa_{a,k}\Delta
t_{k,j}N_{k}\frac{\nu_{k}}{\Delta Vp^{t}_{k}}$ (166)

with $s_{k}=1$ for $\nu_{e}$, $s_{k}=-1$ for $\bar{\nu}_{e}$, and $s_{k}=0$
otherwise. With these choices, the source terms in low-density regions
($\kappa\Delta t\ll 1$) are a lot less noisy. Indeed, every packet passing
through such a region contributes to the source terms exactly the expectation
value of these source terms, instead of creating large source terms when a
packet actually interacts with the fluid and none otherwise. Importantly, the
time-averaged value of the source terms is the same as in the previous method
Foucart:2021mcb . We decrease the shot noise in the source terms at the cost
of loosing exact conservation laws. We note that in this latter technique, an
important decision is when to switch a packet from contributing to one grid
cell to contributing to its neighbor. The obvious answer would be to do so at
the exact time a packet crosses a cell boundary, but this is not necessarily
trivial to predict in high-opacity and/or high-curvature regimes. Any
algorithm using this method thus has to consider a balance between the cost of
accurately determining the time at which a packet crosses a cell boundary, and
the cost of using approximate values for the $\Delta t_{k,j}$.

As more Monte-Carlo simulations are performed in the future, other methods are
likely to be developed. At the moment, exact conservation appears preferable
if the neutrinos contain a significant fraction of the total energy /momentum
of the system, or are dynamically important to the evolution of the system,
i.e. if a lot of packets are used in the simulations. In mergers and post-
merger remnants, where neutrinos are not dynamically important and we tend to
use few packets in low-density regions, the second method may be preferable.

####  4.3.3 Optically thick regime

The optically thick regime poses a more fundamental problem for Monte-Carlo
methods. One issue is that if $(\kappa_{a}+\kappa_{s})\Delta t=N$, we expect
$\sim N$ neutrino-matter interactions per time-step. In dense or hot regions
of mergers, we can easily get $N\gg 1$, making evolutions costly. A related
issue is that if we calculate $\eta,\kappa_{a},\kappa_{s}$ using the
properties of the fluid at the beginning of a time step, and the source terms
lead to large changes in the fluid variables, the evolution can become
numerically unstable. To avoid taking extremely small time steps, we would
then need to at the very least obtain a good guess for the fluid variables at
the end of the time step, or even iteratively solve for these values. The
latter would be costly for Monte-Carlo simulations, as all packets have to be
reevolved for every new guess of the fluid variables. We note that in the
merger context, these issues are mostly due to the presence of a dense, hot
neutron star remnant. For their simulations of post-merger accretion disks,
Miller et al Miller:2019gig  simply require a time step smaller than the
cooling timescale of any given cell, and find that condition to be generally
less constraining than the Courant condition. Accordingly, the problem of very
optically thick regions has only really been encountered in the merger
simulations of Foucart et al Foucart:2020qjb , and handled with approximate
methods inspired by work in other astrophysical contexts. We briefly summarize
their methods here, but note that given the relative novelty of the problem,
these methods remain poorly tested and understood.

A first change that can be made to the standard Monte-Carlo algorithm in dense
region, appropriate when $\kappa_{s}\Delta t\gg 1$, is to once more rely on
the idea that, in that regime, neutrinos are slowly diffusing in the reference
frame of the fluid. Richers et al Richers:2017awc  and Foucart Foucart:2017mbt
propose similar methods to treat this regime. Both assume that in regions of
sufficiently high $\kappa_{s}\Delta t$, neutrinos are advected with the fluid,
while undergoing a slow random walk away from the fluid’s motion. Richers et
al Richers:2017awc  draw from the distribution of times needed for a neutrino
to diffuse a certain distance in the fluid frame. Foucart Foucart:2017mbt
draws from the distribution of distances that the neutrinos move in the fluid
frame after a fixed time, and additionally draws the final momentum of the
diffusing neutrinos from a distribution function calibrated to a solution of
the full Boltzmann equation. Both algorithms also correct the solution of the
diffusion equation so that neutrinos cannot move faster than the speed of
light (i.e. effectively assume that if the diffusion equation predicts
superluminal motion, no scattering occurs and neutrinos just propagate at the
speed of light along a geodesic). We refer the reader to Richers:2017awc ;
Foucart:2017mbt  for the exact choices of distribution functions. When using
such an approximation, calculating the source terms for coupling to the fluid
can be slightly more involved. If using the actual value of the momentum
transfer, we can still compute $\Delta p^{\mu}$ of each packet between the
beginning and end of a step, but should in theory correct that result for the
natural change in $p^{\mu}$ due to the evolution of a packet along a geodesic
in curved spacetime. If using expectation values of the interaction rates
instead, Foucart Foucart:2017mbt  uses the approximation that $\langle
p^{\mu}\rangle=\nu u^{\mu}$ for a fraction $f_{\rm adv}$ of the time step, and
$p^{\mu}=\nu(u^{\mu}+l^{\mu})$ for the remaining of the time step, with
$f_{\rm adv}$ and $l^{\mu}$ chosen so that the packet ends at the desired
location. The source terms can then be computed as in the previous section,
using these approximate values of the neutrino 4-momentum. Overall, this
process has been demonstrated to work well even down to $\kappa_{s}\Delta
t\sim 3$ Foucart:2017mbt .

Regions with $\kappa_{a}\Delta t\gg 1$ create a more difficult challenge. In
these regions, we expect neutrinos to reach statistical equilibrium with the
fluid on a time scale $\kappa_{a}^{-1}\ll\Delta t$, and the temperature and
composition of the fluid themselves may quickly change as more neutrinos or
antineutrinos are emitted, absorbed, or diffuse away. In these regions, we
expect the energy density of neutrinos of a given energy to be set by
$\eta/\kappa_{a}$, and their diffusion timescale to be set by
$(\kappa_{a}+\kappa_{s})$. Then, the evolution of the system can be
approximately captured if we accept that we cannot resolve what happens on
timescales much shorter than $\Delta t$ Foucart:2021mcb , and rely instead on
changes to the source terms inspired by Implicit Monte-Carlo methods
1971JCoPh…8..313F . In Foucart:2021mcb , we propose the transformation

$\displaystyle\eta^{\prime}$ $\displaystyle=$ $\displaystyle a\eta$ (167)
$\displaystyle\kappa_{a}^{\prime}$ $\displaystyle=$ $\displaystyle
a\kappa_{a}$ (168) $\displaystyle\kappa_{s}^{\prime}$ $\displaystyle=$
$\displaystyle\kappa_{s}+(1-a)\kappa_{a}$ (169)

which guarantees that $\eta^{\prime}/\kappa_{a}^{\prime}=\eta/\kappa_{a}$ and
$\kappa_{a}+\kappa_{s}=\kappa_{a}^{\prime}+\kappa_{s}^{\prime}$, but modifies
the equilibration time scale of neutrinos from $\kappa_{a}^{-1}$ to
$(a\kappa_{a})^{-1}$. We then choose $a$ such that $\kappa_{a}^{\prime}\Delta
t\lesssim 1$, making all relevant time scales longer than the time step. We
note that a different value of $a$ may be used for each energy bin and each
neutrino species. This method has the advantage of limiting the number of
emissions and absorptions of neutrinos to roughly what is needed to get to
statistical equilibrium in a few time steps, while avoiding stiff source terms
in the fluid evolution equations. For neutron star merger remnants, we find
sub-percent errors in the neutrino luminosities when comparing this method to
the solution of Boltzmann’s equation Foucart:2021mcb . We note however that
this method could easily impact the diffusion rate of neutrinos within the
densest regions if the neutrino spectrum changes on scales smaller than the
grid scale, and that its accuracy in the presence of significant inelastic
scattering is also uncertain. As opposed to the more complex method used by
Fleck & Cummings 1971JCoPh…8..313F , the method proposed here is not a direct
discretization of the coupled fluid-radiation equations, though deviations
from those equations are all related to sub grid scale and sub time step
features of the solution. There is no doubt that further improvements and/or
tests of the method would be benefitial to properly understand the errors it
creates.

####  4.3.4 Pair processes and other non-linear source terms

In this last section discussing Monte-Carlo methods, we turn to pair processes
as an example of issues that may arise when non-linear terms in the neutrino
distribution functions come into play. In theory, Monte-Carlo methods provide
us with a direct discretization of the neutrino distribution functions in 6D.
One could thus evaluate blocking factors and non-linear source terms
explicitly. However, in the presence of a large number of packets, this is an
expensive computation. If on ther other hand only a small number of packets
are present, sampling noise in the distribution function may lead to large
errors in the resulting reaction rates.

Consider for example the $\nu_{i}\bar{\nu}_{i}\rightarrow e^{+}e^{-}$
reaction. The absorption cross section for neutrinos of 4-momentum
$p^{\alpha}$ could be written, under the same assumptions as Eq. 29 and
assuming that each packet of antineutrinos represent a spatially uniform
distribution of neutrinos within a grid cell of volume $V$

$\kappa_{\nu\bar{\nu}}=\frac{DG_{F}^{2}}{3\pi}\frac{1}{\sqrt{-g}Vp^{t}}\sum_{k,\bar{\nu}_{i}\in
V}N_{k}\frac{(-p^{\alpha}\bar{p}_{\alpha,k})^{2}}{\bar{p}^{t}_{k}}$ (170)

with $\bar{p}_{k}^{\alpha}$ the 4-momentum of antineutrinos in packet $k$, and
the sum covering all antineutrino packets of the correct species in the chosen
grid cell. We note that as opposed to previous opacities, which were computed
in the fluid frame, this opacity is computed in the simulation frame, i.e. if
we draw the time to the next annihilation event, we should use $\Delta
t_{\nu\bar{\nu}}=\kappa^{-1}_{\nu\bar{\nu}}\ln{r}$ with $r\in[0,1)$. While the
mathematical expression is relatively simple, this requires for each packet
iteration over all packets of antineutrinos in the same cell. This is costly
if the number of packets is large, and sensitive to shot noise if the number
of packets is small. We note however that according to our earlier expression
for the stress-energy tensor in Monte-Carlo, this is equivalent to

$\kappa_{\nu\bar{\nu}}=\frac{DG_{F}^{2}}{3\pi}\frac{p^{\alpha}p^{\beta}}{p^{t}}\bar{T}_{\alpha\beta}$
(171)

with $\bar{T}_{\alpha\beta}$ the stress-energy tensor of the antineutrinos.
This allows for faster calculations for large numbers of antineutrinos if we
precompute $\bar{T}_{\alpha\beta}$ (the algorithm will scale linearly with the
number of packets, instead of scaling with the number of packets squared).
Additionally, this method allows us to approximate $\bar{T}_{\alpha\beta}$ in
other ways in regions where few packets are available, e.g. using an average
over a larger spatial region, or a longer time interval, in order to gather
information from a larger number of packets at the cost of a smoothing of the
neutrino distribution function. Quite importantly however, as the opacity
depends explicitly on the 4-momentum of the neutrinos, we have to recompute it
for each packet, at each time. To avoid costly contraction with the metric, a
code that stores $p^{t}$ and $p_{i}$ would ideally express this opacity as

$\kappa_{\nu\bar{\nu}}=\frac{DG_{F}^{2}}{3\pi
p^{t}}\left([p^{t}\alpha]^{2}\bar{E}-2[p^{t}\alpha][p_{i}\bar{F}^{i}]+p_{i}p_{j}\bar{P}^{ij}\right)$
(172)

and thus specifically precompute $\bar{E},\bar{F}^{i},\bar{P}^{ij}$ to
minimize computations. These expressions have however not been used in general
relativistic merger simulations so far.

A potential complication is that if pair processes are included in the
emission rate of neutrinos (as is typically done for muon and tau neutrinos),
and the absorption opacity is calculated according to Kirchoff’s law, then the
opacity for pair annihilation is already included in the simulations under the
assumption that neutrinos and antineutrinos are both in equilibrium with the
fluid. So simply adding the opacity calculated here to a simulation would
double count pair annihilation in optically thick regions. On the other hand,
we know that when assuming equilibrium we underestimate the rate of pair
annihilations by many orders of magnitude in low-density regions above the
remnant. Without a fully self-consistent calculation of pairs everywhere
(which would require the inclusion of blocking factors in the emissivity, and
thus make the emission step dependent on the current distribution of
neutrinos), the method to calculate pair annihilation outlined here should
only be used in optically thin regions.

####  4.3.5 Discussion

The use of Monte-Carlo algorithms in merger and post-merger simulations is a
relatively novel development, with few simulations published so far. Early
results however indicate that Monte-Carlo methods can be used at a
surprisingly low cost, comparable to or lower than that of the most complex
moment schemes, while automatically taking into account the energy dependence
of the distribution function $f_{\nu}$. Monte-Carlo simulations have already
allowed important test of the accuracy of moment schemes, and will certainly
continue to remain a valuable tool until the evolution of Boltzmann’s
equations with higher-order methods becomes a realistic prospect.
Nevertheless, these methods have their own drawbacks. For monoenergetic
neutrinos, they are certainly less accurate than moment schemes in high-
density regimes, were approximations are currently made to avoid the use of
implicite time stepping and the calculations of large number of interactions.
For more realistic neutrino spectra, existing tests of these approximations
indicate that they are probably subdominant sources of error in current merger
simulations, but only a few of these tests have been performed so far. Further
improvements to the behavior of Monte-Carlo algorithms in dense regions,
possibly coupled to the use of implicit methods for their coupling to the
fluid, may be desirable in the future.

Additionally, one of the main supposed advantage of an evolution of the full
Boltzmann equations is the availability of $f_{\nu}$. For existing Monte-Carlo
scheme, this availability is doubtful. Current simulations use a very small
number of packets in low-density regions, so that $f_{\nu}$ can only be
recovered with reasonably low statistical noise if one averages over long time
scales and/or spatial scales. We have noted the importance of this problem for
pair annihilation, but this could also be an issue for calculations of
oscillations due to fast flavor instabilities. The FFI is typically triggered
due to changes in the sign of the net lepton flux (flux of $\nu_{e}$ minus
flux of $\bar{\nu}_{e}$), and calculations of that net flux from existing
Monte-Carlo simulations would be entirely dominated by statistical noise.
While Monte-Carlo simulations are certainly an important step forward in our
modeling of radiation transport, they are thus far from a one-size-fit-all
solution to the problem of radiation transport in merger simulations.

##  5 General relativistic merger simulations

In the previous sections, we provided a detailed discussion of the three broad
classes of algorithms used in general relativistic simulations of neutron star
binary mergers so far. These algorithms have been used in a large number of
simulations, and discussing each of these results goes beyond the scope of
this review. However, in order to put these algorithms in context, and provide
examples of their limitations and known sources of errors, it is useful to
discuss at least two aspects of these simulations: what they broadly tell us
about neutrino physics in neutron star mergers, and how different algorithms
compare when used to simulate the same physical configuration.

###  5.1 Neutrinos in neutron star merger simulations

Neutrinos play two major roles in the evolution of the properties of post-
merger remnants: cooling the remnant disk and, if present, the remnant neutron
star, and modifying the composition of the fluid. There is broad agreement
that, for systems that rapidly collapse to a black hole, or for black hole-
neutron star systems, shock heating during the merger and the circularization
of the accretion disk will lead to bright neutrino emission peaking at
$10^{53-54}\,{\rm erg/s}$. That emission however rapidly decays on $\sim
10\,{\rm ms}$ timescale Deaton:2013sla ; Sekiguchi:2016bjd ; Radice:2021jtw .
Late time emission depends on the relatively poorly constrained heating due to
magnetically-driven turbulence in the post-merger disk, as well as on the mass
of the disk itself Fernandez:2020oow ; Shibata:2021xmo ; Fujibayashi:2022ftg .
It may remain as high as $\sim 10^{52-53}\,{\rm erg/s}$ for $O(100\,{\rm
ms})$. This emission is not sufficient to create a thin disk (the disk aspect
ratio remains $\sim 0.2$), but simulations that do include cooling find
remnant disks that are significantly more compact than simulations without
cooling, indicating a loss of gravitational binding energy. Once the mass
accretion rate drops sufficiently (to roughly $\dot{M}\lesssim
10^{-3}M_{\odot}\,{\rm/s}$ De:2020jdt ), neutrino emission is no longer
sufficient to cool the disk, which becomes an advection dominated accretion
flow.

For systems with a massive neutron star remnant, the peak emission is at a
level comparable to the black hole-disk remnant, but emission from the neutron
star can continue over much longer timescales. For example, the long
axisymmetric simulations of Fujibayashi et al Fujibayashi:2020dvr  find
neutrino luminosities of $\sim 10^{52-53}\,{\rm erg/s}$ more than $5\,{\rm
ms}$ post-merger, with no sign of the luminosity decreasing on those
timescales. The neutron star remnant will thus be the dominant source of
neutrinos after $O(100\,{\rm ms})$ (see Figure 6).

Figure 6: Total neutrino luminosity from the neutron star (solid curves) and
accretion disks (dashed curve) in three long simulations of NSNS merger
remnants in which the central object remains a neutron star. We observe both
the decay of the disk emission on $O(100\,{\rm ms})$ timescales and the
sustained emission from the neutron star. Figure reproduced from Fujibayashi
et al (2020) Fujibayashi:2020dvr .

In both cases, simulations generally agree on an energy hierarchy
$\epsilon_{\nu_{x}}>\epsilon_{\bar{\nu}_{e}}>\epsilon_{\nu_{e}}$, with average
energy of $\sim 10\,{\rm MeV}$ for $\nu_{e}$ but potentially above $20\,{\rm
MeV}$ for $\nu_{x}$. This is simply due to the higher absoption opacity of the
fluid to $\nu_{e}$, which puts the $\nu_{e}$ neutrinosphere farther out in the
remnant than the $\bar{\nu}_{e}$ neutrinosphere (and even more the $\nu_{x}$
neutrinosphere). As the temperature of the remnant close to the neutrinosphere
increases with density, heavy-lepton neutrinos have higher energies.

Weak reactions during the merger and early post-merger evolution are typically
well out of equilibrium, and existing simulations find an initial
overabundance of $\bar{\nu}_{e}$ emission over $\nu_{e}$ emission as the
remnant’s electron fraction rapidly increases - for both NSBH and BNS mergers,
and regardless of whether a remnant black hole is formed or not. After
formation of an accretion disk, however, different regimes can be found. Early
in the evolution of a massive accretion disk, when electrons in the inner disk
are still degenerate and the disk is optically thick to neutrinos, the
composition of the inner disk self-regulates to a relatively neutron-rich
value $Y_{e}\sim 0.1$ 2003ApJ…588..931B ; 2018ApJ…858…52S , while the lower
density outer disk and the hotter corona have significantly higher $Y_{e}$
(Figure 7). For optically thin disks lacking degenerate matter, on the other
hand, the inner disk immediately evolves to higher $Y_{e}\sim 0.15-0.2$
FoucartM1:2015 ; Just:2021cls . These values can also be significantly
impacted by irradiation of the disk by a central neutron star (Figure 7).
Farther out in the disk, or once neutrino emission becomes inefficient and the
disk becomes advection dominated, the electron fraction largely freezes out.
The initial density and mass accretion rate of the disk thus plays a role in
determining the composition of post-merger outflows Fernandez:2016sbf , and
that composition is also affected by the location from which the matter is
ejected (midplane vs corona). In that respect, it is worth noting that
outflows produced from self-consistent MHD simulations and outflows produced
by simulations using $\alpha$ viscosity models, even when they agree on the
amount of matter ejected, make very different predictions for the history of
the fluid elements ejected from an accretion disk – which may impact
predictions of the final $Y_{e}$ of the outflows.

Figure 7: Electron fraction in post-merger remnants. We show a vertical slice
through a NSBH merger (Left) and a BNS merger (Right). The BNS merger is
evolved using two different choices for the energy closure in a two-moment
scheme. Figures reproduced from Foucart et al (2015) FoucartM1:2015  and
Foucart et al (2016) Foucart:2016rxm .

All of the above features are found in simulations using leakage, moments, or
Monte-Carlo transport, although disagreements between methods can be found on
the exact neutrino luminosity and composition of the remnant (see below).
Simulations that include reabsorption of neutrinos in the matter outflows also
find that hot outflows originating from the colliding cores of two neutron
stars or the hot corona of an accretion disk rapidly evolve to electron
fractions $Y_{e}\sim(0.2-0.4)$Wanajo:2014 ; FoucartM1:2015 ; Sekiguchi:2016bjd
; Radice:2018pdn ; Foucart:2020qjb ; Radice:2021jtw ; Camilletti:2022jms ,
though the exact value of $Y_{e}$ can vary significantly depending on the
chosen numerical algorithm (see below). Faster and colder outflows associated
with the tidal disruption of a low-mass neutron star by a more massive
companion (in either NSBH or BNS mergers) do not capture enough neutrinos to
undergo sigificant changes of composition, and thus remain very neutron-rich
($Y_{e}\lesssim 0.05$), even in simulations that include neutrino absorption
Sekiguchi:2016bjd ; Foucart:2016vxd ; Radice:2018pdn ; Camilletti:2022jms .
There is higher uncertainty on the composition of the post-merger outflows.
Neutrino-driven outflows observed in post-merger simulations that do not
include magnetic fields tend to have very high $Y_{e}$ 2015MNRAS.448..541J ;
FoucartM1:2015  (see e.g. polar regions of Fig. 7), but matter outflows in the
early post-merger phase are likely dominated by magnetically-driven winds
2018ApJ…858…52S ; Fernandez:2018kax , and viscous outflows in the advection-
dominated phase Fernandez:2013tya . Magnetically-driven outflows are faster
($\sim 0.1c$), and have less time to absorb neutrinos, while viscous outflows
are slower ($\sim 0.03c$) and have a composition set by the electron fraction
of the fluid at the time at which weak interactions freeze-out. Long 3D
simulations including magnetic fields and approximtely accounting for neutrino
absorption in post-processing find outflows with $Y_{e}\sim 0.15-0.25$,
starting from very neutron-rich initial conditions 2018ApJ…858…52S . Shorter
simulations using an advanced Monte-Carlo transport scheme and including
magnetic fields Miller:2019dpt , with similar initial conditions, find a
broader distribution of $Y_{e}$ peaking just below $Y_{e}\sim 0.2$ and
extending up to $Y_{e}\sim 0.4$. Similarly, long simulations of a NSBH merger
remnant with a moment scheme and magnetic fields Hayashi:2021oxy , initialized
from the outcome of a merger simulation, find a broad $Y_{e}\sim 0.15-0.4$
distribution peaking just above $Y_{e}\sim 0.2$. Very few long 3D simulations
including both magnetic fields and neutrinos are however available. Parameter
space exploration with axisymmetric simulations (using artificial viscosity
instead of evolving the magnetic field) find a significant sensitivity of the
results on the compactness of the disk Fernandez:2020oow , the lifetime of the
massive neutron star remnant (if present) Metzger:2014ila , and the initial
composition used in the simulation Fernandez:2016sbf . 3D simulations have
also demonstrated the importance of the large scale structure of the post-
merger magnetic field both with a black hole remnant Christie:2019lim  and a
neutron star remnant deHaas:2022ytm . Given our limited understanding of the
properties of post-merger remnant as a function of the initial binary
configuration, it is thus fairly difficult at this point to build reliable
models of post-merger outflows.

Finally, in the previous sections we already emphasized the difficulties of
properly implementing $\nu\bar{\nu}$ pair annihilation in relativistic merger
simulations. Fujibayashi et al Fujibayashi:2017xsz , using the approximate
moment method described in Sec. 4.2.8, finds that pair annihilation can
accelerate the matter in the polar regions to mildly relativistic speeds
(Lorentz factor $\Gamma\sim 2$) – sufficient to be important to the dynamics
of the fluid in those low-density regions, but not enough to power SGRBs (see
Fig. 5). Foucart et al Foucart:2018gis  find that the energy deposition
predicted by the approximate moment scheme is likely accurate within a factor
of $\sim 2-3$. More advanced transport calculations had earlier been performed
on a prescribed fluid background 1999ApJ…518..356P ; 2006JPhG…32..443K ;
Zalamea:2008dq ; Dessart:2008zd ; just:16 ; Perego:2017fho . These simulations
find that at most a few percents of the neutrino luminosity is deposited in
the polar region through neutrino annihilation, with a rapid drop in
efficiency as the accretion rate decreases. These two sets of results are
qualitatively consistent, indicating that while neglecting pair annihilation
may be acceptable when considering only the low-velocity outflows powering
kilonovae, it may be important to take into account when attempting to resolve
the evolution of the polar regions, and in particular the formation of jets in
that region.

###  5.2 Comparisons of numerical algorithms

To finish our discussion of general relativistic transport methods for neutron
star mergers, let us look into a few direct comparisons of numerical
algorithms to gauge their accuracy.

A direct comparison of a leakage scheme (whithout absorption) and a two-moment
scheme (with Minerbo closure, evolving $\tilde{E}$ and $\tilde{F}_{i}$) was
performed in the context of a low-mass neutron star merger in Foucart:2015gaa
. In that simulation, leakage and moment schemes agreed quite well on the
$\bar{\nu}_{e}$ luminosity, while other luminosities varied by factors of
$2-3$ in the first $10\,{\rm ms}$ following the merger. The inclusion of
neutrino absorption led to the production of a neutrino-driven wind that did
not exist in the leakage simulation, but the mass outflow rate was only
$\dot{M}\sim 10^{-2}M_{\odot}/s$, i.e. less than what one might expect at that
time from magnetically-driven winds. Outflows in the moment simulation were
also significantly hotter and less neutron-rich than in the leakage simulation
($\langle Y_{e}\rangle=0.2$ vs $\langle Y_{e}\rangle=0.1$, and $\langle
s\rangle=20k_{B}$ vs $\langle s\rangle=10k_{B}$ per baryon). This confirms
that neutrino luminosities are only order-of-magnitude accurate in simple
leakage schemes, and the crucial impact of neutrino absorption on the
properties of matter outflows.

In Radice:2021jtw , the authors consider two neutron star mergers, one for
which the remnant collapses to a black hole a few milliseconds after contact,
and one forming a long-lived neutron star remnant. They compare results using
a hybrid moment-leakage scheme, the standard two-moment scheme with Minerbo
closure, as well as a two-moment scheme using the Eddington closure (i.e. the
optically thick closure everywhere). The hybrid schemes overestimates neutrino
luminosities by a factor of $\sim 2$, while the two-moment simulations with
different closures are in much better agreement ($10\%-30\%$ differences, see
Fig. 8). All three schemes agree very well on the average energy of escaping
neutrinos for $\nu_{e}\bar{\nu}_{e}$, while early emission of higher energy
heavy-lepton neutrinos is predicted by the two-moment scheme but not the
hybrid scheme. The two-moment scheme also predicts $\sim 30\%$ less mass
ejected than the hybrid scheme, and significantly higher electron fractions
($\Delta Y_{e}\sim 0.1$). The use of the Eddington or Minerbo closure had
again a much smaller effect ($\sim 10\%$ relative error in mass, shifts of a
few percents in $Y_{e}$). Differences between the hybrid and two-moment scheme
are thus slightly lower than between a pure leakage and a two-moment scheme,
but nonetheless significant. The fact that the Minerbo and Eddington closure
are in much closer agreement is however quite encouraging, as a reasonable
proxy for the error due to the use of an approximate closure in semi-
transparent regions.

Figure 8: Neutrino luminosities for $\nu_{e}$, $\bar{\nu}_{e}$ and $\nu_{x}$
in a BNS merger with a long-lived remnant. Results are shown for three
algorithms: a hybrid leakage-moment scheme (blue), a two-moment scheme with
Minerbo closure (orange), and a two-moment scheme with Eddington closure
(green). Figures reproduced from Radice et al (2022) Radice:2021jtw .

The impact of the chosen energy closure in grey two-moment schemes was
investigated in Foucart:2016rxm , for the same physical configuration as in
Foucart:2015gaa . In that manuscript, two energy closures are considered. In
the first, the average energy of neutrinos is taken from a black-body
distribution at temperature $T_{\nu}=\max{(T_{\rm leak},T)}$, with $T_{\rm
leak}$ the temperature predicted by a leakage scheme (globally) and $T$ the
temperature of the fluid. In the second, the neutrino number density is
evolved to obtain a local estimate of $T_{\nu}$, and the neutrino spectral
index is evolved to attempt to accounting for softer spectra in the regions
optically thick to scattering but optically thin to absorption. The luminosity
of neutrinos initially is in reasonable agreement in both schemes ($20\%-30\%$
differences), but differences increase over time to $\sim 50\%$ at the end of
the evolution ($8\,{\rm ms}$ post-merger). This is less than the difference
between leakage and two-moment schemes, but still quite significant. One
likely source of error here is divergence in the evolution of $Y_{e}$ in the
remnant, due in part to the fact that the scheme that does not evolve the
number density cannot explicitly conserve the total lepton number of the
system. As in Radice:2021jtw , the average energy of electron-type neutrinos
is reasonably close in both schemes, but calculations based on the leakage
scheme underestimate the initial energy of heavy-lepton neutrinos. Both
schemes have similar outflow masses, but difference in electron fraction
$\Delta Y_{e}\sim(0.05-0.1)$ (see Fig. 7). We note that this is despite the
fact that the estimated average neutrino energies are similar in both
simulations, and mostly due to the fact that polar neutrinos (which interact
with hot outflows) have higher energy than equatorial neutrinos (which do
not), and are thus more strongly absorbed than one would think using a global
estimate of the average energy. Using the local fluid temperature to estimate
the neutrino energies would lead to significantly larger errors: factors of a
few changes in the neutrino energies instead of tens of percent.

One can also take a broader view of these comparisons between transport
schemes. Instead of directly comparing numerical simulations, Nedora al
Nedora:2020qtd  compare datasets of simulations using different microphysical
inputs, and provide numerical fits for the outcome of these simulations. Their
results are consistent with the direct comparisons discussed above, and show
that the choice of neutrino transport algorithm in merger simulations remain
an important source of error in outflow modeling.

Finally, comparisons of a Monte-Carlo transport scheme with a two-moment
scheme evolving the number density were performed without back-reaction of the
Monte-Carlo code to the fluid (i.e. using the same fluid evolution as in the
moment evolution) Foucart:2018gis  as well as with coupling to the fluid
Foucart:2020qjb  in a neutron star merger simulation. In the fully coupled
simulation, the luminosity of $\nu_{e}$ and $\bar{\nu}_{e}$ and the energy of
neutrinos was in better agreement than in most previous comparisons, excepting
the impact of the pressure closure ($10\%-20\%$ differences). The heavy-lepton
luminosity showed $\sim 50\%$ changes, showing once more the difficulty of
properly capturing the evolution of those neutrinos. Similarly, the
composition of the outflows was in much better agreement, with the average
$Y_{e}$ changing by $\Delta Y_{e}\sim 0.03$, and the maximum $Y_{e}$ by
$\Delta Y_{e}\sim 0.05$ – not necessarily negligible changes, but a
significantly better agreement than in other comparisons (see Fig. 9).

Figure 9: Comparison of the electron fraction $Y_{e}$ in a vertical slice of a
post-merger remnant for simulations using Monte-Carlo methods (Left) and a
two-moment scheme (Middle). The right panel shows the $Y_{e}$ distribution of
matter outflows in both simulations. Figures reproduced from Foucart et al
(2020) Foucart:2020qjb .

The simulation that does not fully couple the Monte-Carlo algorithm to the
fluid evolution, thus avoiding any drift of the fluid variables over time due
to small differences between the algorithms, shows similar differences in the
neutrino energies, but better agreement for the $\nu_{x}$ luminosity. A more
detailed study of the spatial distribution of neutrinos however indicates that
the moment scheme greatly overestimates the density of neutrinos close to the
pole (by $\sim 50\%-100\%$), and underestimates their density farther out – an
important consequence of the choice of closure made for the pressure tensor.
That simulation also computed the rate of neutrino pair annihilation using the
moment scheme and the Monte-Carlo methods. We note that, as shown in Sec.
4.2.8, the moment calculation requires the use of both an approximate average
energy for the pairs and of the approximate pressure closure in a regime in
which it is inaccurate. Additionally, it is naturally impacted by errors in
the estimated energy density of neutrinos close to the poles. Interestingly,
for the specific numerical algorithms studied here, those errors partially
cancelled out, leaving us with factors of $2-3$ errors in the actual
annihilation rate for the moment scheme. There is however now guarantee that
this would also be true for a different binary configuration, or with
different estimates for the neutrino energy, as individually some of these
errors can modify the annihilation rate by close to an order of magnitude (see
Fig. LABEL:MCM1Lum).

Figure 10: Left: Energy spectrum of the $\nu_{x}$ neutrinos $14\,{\rm ms}$
after a BNS merger, according to a Monte-Carlo simulation. Vertical lines show
the average energy esimated using Monte-Carlo (Solid line) or a two-moment
scheme (dashed line). Right: Angular distribution of the neutrinos in the same
simulations, using a Monte-Carlo evolution (black) or a two-moment evolution
(red). Figures reproduced from Foucart et al (2018) Foucart:2018gis .

Overall, we thus note that there has been a significant improvement in the
magnitude of those errors, with the more modern moment Foucart:2016rxm ;
Radice:2021jtw  and Monte-Carlo Foucart:2020qjb  schemes getting estimated
relative errors in the $\sim 10\%-20\%$ range for the variables that most
significantly impact astrophysical observables – by which point neutrino
transport is certainly a subdominant source of error when compared to the
underresolved evolution of magnetic fields, the nuclear physics uncetainties
in kilonova modeling, or maybe even the impact of neglected / poorly modeled
neutrino physics (oscillations, inelastic scattering, pair process).

##  6 Conclusions

The inclusion of radiation transport algorithm in neutron star merger
simulations has taken significant step forward over the last decade. The
development of improved two-moment schemes and Monte-Carlo algorithms, in
particular, allows for reasonably accurate evolution of the transport
equations for relatively simple neutrino physics.

This leaves us however with a few important challenges. First, very few
simulations have made use of these new methods, and thus efforts to model the
observable counterparts to neutron star mergers still heavily rely on results
obtained with simpler microphysics. As a result, current model are often
unreliable Henkel:2022naw , and dependent on the algorithms used for neutrino
evolution in the simulations used to calibrate them Nedora:2020qtd . Second,
we know that a number of potentially important processes are not included in
existing simulations, or are poorly modeled in those simulations. This include
at least neutrino oscillations, pair annihilation, inelastic scattering, and
potentially direct and modified URCA processes. Finally, neutrino transport is
only one part of the problem when modeling neutron star mergers. Magnetic
fields are crucial to the evolution of neutron star mergers and their post-
merger remnants. The growth of large scale magnetic fields due to MHD
instabilities is not sufficiently resolved even in simulations that do not
include neutrinos, or that include them very approximately. Combining high-
resolution MHD simulations with our most advanced neutrino transport schemes
over the seconds time scales needed to follow the evolution of a post-merger
remnants remains an extremely difficult problem that will likely remain an
important source of uncertainty in our modeling of electromagnetic signals
from neutron star mergers for years to come.

Acknowledgments

The author is grateful to Adam Burrows, Gail McLaughlin, Philipp Moesta, Elias
Most and Sherwood Richers for valuable discussions on various aspects of this
review during its preparation, and to Sho Fujibayashi, Kota Hayashi, and David
Radice for allowing reproduction of their results in this manuscript. The
author acknowledges support from the DOE Office of Science, Office of Nuclear
Physics, under contract number DE-AC02-05CH11231, from NASA through grant
80NSSC18K0565, and from the NSF through grant AST-2107932. Part of this work
was performed at the Aspen Center for Physics, which is supported by National
Science Foundation grant PHY-1607611.

## References

  *   * (1) B.P. Abbott, et al., GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs.  Phys. Rev. X 9(3), 031,040 (2019).  10.1103/PhysRevX.9.031040.  arXiv:1811.12907 [astro-ph.HE]
  * (2) R. Abbott, et al., GWTC-2.1: Deep Extended Catalog of Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run (2021).  arXiv:2108.01045 [gr-qc]
  * (3) R. Abbott, et al., GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run (2021).  arXiv:2111.03606 [gr-qc]
  * (4) B. Abbott, et al., GW170817: Measurements of neutron star radii and equation of state.  Phys. Rev. Lett. 121(16), 161,101 (2018).  10.1103/PhysRevLett.121.161101.  arXiv:1805.11581 [gr-qc]
  * (5) J.M. Lattimer, D.N. Schramm, The tidal disruption of neutron stars by black holes in close binaries.  Astroph.J. 210, 549–567 (1976).  10.1086/154860
  * (6) L.F. Roberts, D. Kasen, W.H. Lee, E. Ramirez-Ruiz, Electromagnetic transients powered by nuclear decay in the tidal tails of coalescing compact binaries.  Astrophys. J. Lett. 736, L21 (2011).  10.1088/2041-8205/736/1/L21.  arXiv:1104.5504 [astro-ph.HE]
  * (7) D. Kasen, N.R. Badnell, J. Barnes, Opacities and Spectra of the r-process Ejecta from Neutron Star Mergers.  Astroph.J. 774, 25 (2013).  10.1088/0004-637X/774/1/25.  arXiv:1303.5788 [astro-ph.HE]
  * (8) L.X. Li, B. Paczynski, Transient events from neutron star mergers.  Astroph.J. 507, L59 (1998).  10.1086/311680.  arXiv:astro-ph/9807272 [astro-ph]
  * (9) W. Fong, E. Berger, The Locations of Short Gamma-Ray Bursts as Evidence for Compact Object Binary Progenitors.  Astroph.J. 776, 18 (2013).  10.1088/0004-637X/776/1/18.  arXiv:1307.0819 [astro-ph.HE]
  * (10) D.E. Holz, S.A. Hughes, Using Gravitational-Wave Standard Sirens.  Astroph.J. 629(1), 15–22 (2005).  10.1086/431341.  arXiv:astro-ph/0504616 [astro-ph]
  * (11) S. Nissanke, D.E. Holz, S.A. Hughes, N. Dalal, J.L. Sievers, Exploring Short Gamma-ray Bursts as Gravitational-wave Standard Sirens.  Astroph.J. 725(1), 496–514 (2010).  10.1088/0004-637X/725/1/496.  arXiv:0904.1017 [astro-ph.CO]
  * (12) B.P. Abbott, et al., A gravitational-wave standard siren measurement of the Hubble constant.  Nature 551(7678), 85–88 (2017).  10.1038/nature24471.  arXiv:1710.05835 [astro-ph.CO]
  * (13) K. Hotokezaka, E. Nakar, O. Gottlieb, S. Nissanke, K. Masuda, G. Hallinan, K.P. Mooley, A.T. Deller, A Hubble constant measurement from superluminal motion of the jet in GW170817.  Nature Astron. 3(10), 940–944 (2019).  10.1038/s41550-019-0820-1.  arXiv:1806.10596 [astro-ph.CO]
  * (14) K. Kiuchi, P. Cerdá-Durán, K. Kyutoku, Y. Sekiguchi, M. Shibata, Efficient magnetic-field amplification due to the Kelvin-Helmholtz instability in binary neutron star<EMAIL_ADDRESS>92(12), 124034 (2015).  10.1103/PhysRevD.92.124034.  arXiv:1509.09205 [astro-ph.HE]
  * (15) F. Foucart, M.D. Duez, L.E. Kidder, R. Nguyen, H.P. Pfeiffer, M.A. Scheel, Evaluating radiation transport errors in merger simulations using a Monte Carlo algorithm.  Phys. Rev. D 98(6), 063,007 (2018).  10.1103/PhysRevD.98.063007.  arXiv:1806.02349 [astro-ph.HE]
  * (16) L. Baiotti, L. Rezzolla, Binary neutron star mergers: a review of Einstein’s richest laboratory.  Rept. Prog. Phys. 80(9), 096,901 (2017).  10.1088/1361-6633/aa67bb.  arXiv:1607.03540 [gr-qc]
  * (17) V. Paschalidis, General relativistic simulations of compact binary mergers as engines for short gamma-ray bursts.  Class.Quant.Grav. 34(8), 084002 (2017).  10.1088/1361-6382/aa61ce.  arXiv:1611.01519 [astro-ph.HE]
  * (18) E. Burns, Neutron Star Mergers and How to Study Them.  Living Rev. Rel. 23(1), 4 (2020).  10.1007/s41114-020-00028-7.  arXiv:1909.06085 [astro-ph.HE]
  * (19) É.É. Flanagan, T. Hinderer, Constraining neutron-star tidal Love numbers with gravitational-wave detectors.  Phys.Rev.D 77(2), 021502 (2008).  10.1103/PhysRevD.77.021502.  arXiv:0709.1915
  * (20) S. Bernuzzi, A. Nagar, M. Thierfelder, B. Brügmann, Tidal effects in binary neutron star coalescence.  Phys.Rev.D. 86(4), 044030 (2012).  10.1103/PhysRevD.86.044030.  arXiv:1205.3403 [gr-qc]
  * (21) T. Hinderer, A. Taracchini, F. Foucart, A. Buonanno, J. Steinhoff, M. Duez, L.E. Kidder, H.P. Pfeiffer, M.A. Scheel, B. Szilágyi, K. Hotokezaka, K. Kyutoku, M. Shibata, C.W. Carpenter, Effects of neutron-star dynamic tides on gravitational waveforms within the effective-one-body approach.  Phys. Rev. Lett. 116(18), 181101 (2016).  10.1103/PhysRevLett.116.181101.  arXiv:1602.00599 [gr-qc]
  * (22) S. Akcay, S. Bernuzzi, F. Messina, A. Nagar, N. Ortiz, P. Ret̃tegno, Effective-one-body multipolar waveform for tidally interacting binary neutron stars up to merger.  Phys. Rev. D 99(4), 044,051 (2019).  10.1103/PhysRevD.99.044051.  arXiv:1812.02744 [gr-qc]
  * (23) J.E. Thompson, E. Fauchon-Jones, S. Khan, E. Nitoglia, F.e. Pannarale, T. Dietrich, M. Hannam, Modeling the gravitational wave signature of neutron star black hole coalescences: PhenomNSBH.  Phys. Rev. D 101(12), 124,059 (2020).  10.1103/PhysRevD.101.124059.  arXiv:2002.08383 [gr-qc]
  * (24) A. Matas, et al., Aligned-spin neutron-star–black-hole waveform model based on the effective-one-body approach and numerical-relativity simulations.  Phys. Rev. D 102(4), 043,023 (2020).  10.1103/PhysRevD.102.043023.  arXiv:2004.10001 [gr-qc]
  * (25) F. Foucart, D. Desai, W. Brege, M.D. Duez, D. Kasen, D.A. Hemberger, L.E. Kidder, H.P. Pfeiffer, M.A. Scheel, Dynamical ejecta from precessing neutron star-black hole mergers with a hot, nuclear-theory based equation of state.  Class. Quant. Grav. 34(4), 044,002 (2017).  10.1088/1361-6382/aa573b.  arXiv:1611.01159 [astro-ph.HE]
  * (26) F. Pannarale, A. Tonita, L. Rezzolla, Black hole-neutron star mergers and short GRBs: a relativistic toy model to estimate the mass of the torus.  Astrophys.J. 727, 95 (2011).  10.1088/0004-637X/727/2/95.  arXiv:1007.4160 [astro-ph.HE]
  * (27) F. Foucart, Black-hole-neutron-star mergers: Disk mass predictions.  Phys.Rev.D 86(12), 124007 (2012).  10.1103/PhysRevD.86.124007.  arXiv:1207.6304 [astro-ph.HE]
  * (28) W.E. East, V. Paschalidis, F. Pretorius, Eccentric Mergers of Black Holes with Spinning Neutron Stars.  ApJ Letters 807, L3 (2015).  10.1088/2041-8205/807/1/L3.  arXiv:1503.07171 [astro-ph.HE]
  * (29) F. Foucart, A brief overview of black hole-neutron star mergers.  Front. Astron. Space Sci. 7, 46 (2020).  10.3389/fspas.2020.00046.  arXiv:2006.10570 [astro-ph.HE]
  * (30) K. Kyutoku, M. Shibata, K. Taniguchi, Coalescence of black hole–neutron star binaries.  Living Rev. Rel. 24(1), 5 (2021).  10.1007/s41114-021-00033-4.  arXiv:2110.06218 [astro-ph.HE]
  * (31) D. Radice, S. Bernuzzi, A. Perego, The Dynamics of Binary Neutron Star Mergers and GW170817.  Ann. Rev. Nucl. Part. Sci. 70, 95–119 (2020).  10.1146/annurev-nucl-013120-114541.  arXiv:2002.03863 [astro-ph.HE]
  * (32) D.M. Siegel, B.D. Metzger, Three-dimensional general-relativistic magnetohydrodynamic simulations of remnant accretion disks from neutron star mergers: Outflows and $r$-process nucleosynthesis.  Phys. Rev. Lett. 119(23), 231,102 (2017).  10.1103/PhysRevLett.119.231102.  arXiv:1705.05473 [astro-ph.HE]
  * (33) R. Fernández, A. Tchekhovskoy, E. Quataert, F. Foucart, D. Kasen, Long-term GRMHD simulations of neutron star merger accretion discs: implications for electromagnetic counterparts\̇lx@bibnewblockMon. Not. Roy. Astron. Soc. 482(3), 3373–3393 (2019).  10.1093/mnras/sty2932.  arXiv:1808.00461 [astro-ph.HE]
  * (34) I.M. Christie, A. Lalakos, A. Tchekhovskoy, R. Fernández, F. Foucart, E. Quataert, D. Kaseñ, The Role of Magnetic Field Geometry in the Evolution of Neutron Star Merger Accretion Discs.  Mon. Not. Roy. Astron. Soc. 490(4), 4811–4825 (2019).  10.1093/mnras/stz2552.  arXiv:1907.02079 [astro-ph.HE]
  * (35) R. Fernández, B.D. Metzger, Delayed outflows from black hole accretion tori following neutron star binary coalescence.  Mon. Not. Roy. Astron. Soc. 435, 502 (2013).  10.1093/mnras/stt1312.  arXiv:1304.6720 [astro-ph.HE]
  * (36) S. De, D.M. Siegel, Igniting Weak Interactions in Neutron Star Postmerger Accretion Disks.  Astrophys. J. 921(1), 94 (2021).  10.3847/1538-4357/ac110b.  arXiv:2011.07176 [astro-ph.HE]
  * (37) B.D. Metzger, R. Fernández, Red or blue? A potential kilonova imprint of the delay until black hole formation following a neutron star merger.  Mon. Not. Roy. Astron. Soc. 441, 3444–3453 (2014).  10.1093/mnras/stu802.  arXiv:1402.4803 [astro-ph.HE]
  * (38) B. Margalit, A.S. Jermyn, B.D. Metzger, L.F. Roberts, E. Quataert, Angular Momentum Transport in Proto-Neutron Stars and the Fate of Neutron Star Merger Remnants (2022).  arXiv:2206.10645 [astro-ph.HE]
  * (39) K. Hayashi, S. Fujibayashi, K. Kiuchi, K. Kyutoku, Y. Sekiguchi, M. Shibata, General-relativistic neutrino-radiation magnetohydrodynamic simulation of seconds-long black hole-neutron star mergers.  Phys. Rev. D 106(2), 023,008 (2022).  10.1103/PhysRevD.106.023008.  arXiv:2111.04621 [astro-ph.HE]
  * (40) W. Fong, E. Berger, R. Margutti, B.A. Zauderer, A Decade of Short-duration Gamma-Ray Burst Broadband Afterglows: Energetics, Circumburst Densities, and Jet O pening Angles.  Astroph.J. 815, 102 (2015).  10.1088/0004-637X/815/2/102.  arXiv:1509.02922 [astro-ph.HE]
  * (41) K.P. Mooley, et al., A Strong Jet Signature in the Late-time Light Curve of GW170817.  Astrophys. J. Lett. 868(1), L11 (2018).  10.3847/2041-8213/aaeda7.  arXiv:1810.12927 [astro-ph.HE]
  * (42) R.D. Blandford, R.L. Znajek, Electromagnetic extraction of energy from Kerr black holes.  Mon.Not.Roy.Astron.Soc. 179, 433–456 (1977)
  * (43) O. Just, M. Obergaulinger, H.T. Janka, A. Bauswein, N. Schwarz, Neutron-star Merger Ejecta as Obstacles to Neutrino-powered Jets of Gamma-Ray Bursts.  Astroph.J.Lett. 816, L30 (2016).  10.3847/2041-8205/816/2/L30.  arXiv:1510.04288 [astro-ph.HE]
  * (44) S. Fujibayashi, Y. Sekiguchi, K. Kiuchi, M. Shibata, Properties of Neutrino-driven Ejecta from the Remnant of a Binary Neutron Star Merger: Pure Radiation Hydrodynamics Case.  Astrophys. J. 846(2), 114 (2017).  10.3847/1538-4357/aa8039.  arXiv:1703.10191 [astro-ph.HE]
  * (45) J. Lippuner, L.F. Roberts, r-process Lanthanide Production and Heating Rates in Kilonovae.  Astroph.J. 815, 82 (2015).  10.1088/0004-637X/815/2/82.  arXiv:1508.03133 [astro-ph.HE]
  * (46) L. Dessart, C.D. Ott, A. Burrows, S. Rosswog, E. Livne, Neutrino signatures and the neutrino-driven wind in Binary Neutron Star Mergers.  Astrophys. J. 690, 1681 (2009).  10.1088/0004-637X/690/2/1681.  arXiv:0806.4380 [astro-ph]
  * (47) J. Barnes, D. Kasen, Effect of a High Opacity on the Light Curves of Radioactively Powered Transients from Compact Object<EMAIL_ADDRESS>775, 18 (2013).  10.1088/0004-637X/775/1/18.  arXiv:1303.5787 [astro-ph.HE]
  * (48) D. Tsang, Shattering flares during close encounters of neutron stars.  The Astrophysical Journal 777(2), 103 (2013).  URL http://stacks.iop.org/0004-637X/777/i=2/a=103
  * (49) C. Palenzuela, L. Lehner, M. Ponce, S.L. Liebling, M. Anderson, D. Neilsen, P. Motl, Electromagnetic and Gravitational Outputs from Binary-Neutron-Star Coalescence.  Phys. Rev. Lett. 111(6), 061105 (2013).  10.1103/PhysRevLett.111.061105.  arXiv:1301.7074 [gr-qc]
  * (50) E.R. Most, A.A. Philippov, Reconnection-powered fast radio transients from coalescing neutron star binaries (2022).  arXiv:2207.14435 [astro-ph.HE]
  * (51) K. Hotokezaka, S. Nissanke, G. Hallinan, T.J.W. Lazio, E. Nakar, T. Piran, Radio Counterparts of Compact Binary Mergers detectable in Gravitational Waves: A Simulation for an Optimized Survey.  Astrophys. J. 831(2), 190 (2016).  10.3847/0004-637X/831/2/190.  arXiv:1605.09395 [astro-ph.HE]
  * (52) R. Fernández, B.D. Metzger, Electromagnetic Signatures of Neutron Star Mergers in the Advanced LIGO Era.  Ann. Rev. Nucl. Part. Sci. 66, 23–45 (2016).  10.1146/annurev-nucl-102115-044819.  arXiv:1512.05435 [astro-ph.HE]
  * (53) N.R. Tanvir, A.J. Levan, A.S. Fruchter, J. Hjorth, R.A. Hounsell, K. Wiersema, R.L. Tunnicliffe, A ‘kilonova’ associated with the short-duration $\gamma$-ray burst GRB130603B.  Nature 500, 547–549 (2013).  10.1038/nature12505.  arXiv:1306.4971 [astro-ph.HE]
  * (54) E. Berger, W. Fong, R. Chornock, An r-process Kilonova Associated with the Short-hard GRB 130603B.  Ap. J. Lett. 774, L23 (2013).  10.1088/2041-8205/774/2/L23.  arXiv:1306.3960 [astro-ph.HE]
  * (55) B.P. Abbott, et al., GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral.  Phys. Rev. Lett. 119(16), 161,101 (2017).  10.1103/PhysRevLett.119.161101.  arXiv:1710.05832 [gr-qc]
  * (56) B. Abbott, et al., GW190425: Observation of a Compact Binary Coalescence with Total Mass $\sim 3.4M_{\odot}$.  Astrophys. J. Lett. 892(1), L3 (2020).  10.3847/2041-8213/ab75f5.  arXiv:2001.01761 [astro-ph.HE]
  * (57) C. Barbieri, O.S. Salafia, M. Colpi, G. Ghirlanda, A. Perego, Exploring the nature of ambiguous merging systems: GW190425 in low latency.  Astron. Astrophys. 654, A12 (2021).  10.1051/0004-6361/202037778.  arXiv:2002.09395 [astro-ph.HE]
  * (58) G. Raaijmakers, et al., The Challenges Ahead for Multimessenger Analyses of Gravitational Waves and Kilonova: A Case Study on GW190425.  Astrophys. J. 922(2), 269 (2021).  10.3847/1538-4357/ac222d.  arXiv:2102.11569 [astro-ph.HE]
  * (59) R. Dudi, A. Adhikari, B. Brügmann, T. Dietrich, K. Hayashi, K. Kawaguchi, K. Kiuchi, K. Kyutoku, M. Shibata, W. Tichy, Investigating GW190425 with numerical-relativity simulations (2021).  arXiv:2109.04063 [astro-ph.HE]
  * (60) A. Camilletti, L. Chiesa, G. Ricigliano, A. Perego, L.C. Lippold, S. Padamata, S. Bernuzzi, D. Radice, D. Logoteta, F.M. Guercilena, Numerical relativity simulations of the neutron star merger GW190425: microphysics and mass ratio effects (2022).  arXiv:2204.05336 [astro-ph.HE]
  * (61) R. Abbott, et al., Observation of Gravitational Waves from Two Neutron Star–Black Hole Coalescences.  Astrophys. J. Lett. 915(1), L5 (2021).  10.3847/2041-8213/ac082e.  arXiv:2106.15163 [astro-ph.HE]
  * (62) A. Burrows, S. Reddy, T.A. Thompson, Neutrino opacities in nuclear matter.  Nucl. Phys. A 777, 356–394 (2006).  10.1016/j.nuclphysa.2004.06.012.  arXiv:astro-ph/0404432
  * (63) P. Vogel, Analysis of the antineutrino capture on protons.  Phys. Rev. D 29, 1918–1922 (1984).  10.1103/PhysRevD.29.1918.  URL https://link.aps.org/doi/10.1103/PhysRevD.29.1918
  * (64) S.W. Bruenn, Stellar core collapse - Numerical model and infall epoch.  Astroph.J.Suppl. 58, 771–841 (1985).  10.1086/191056
  * (65) D.G. Yakovlev, A.D. Kaminker, O.Y. Gnedin, P. Haensel, Neutrino emission from neutron stars.  Phys. Rept. 354, 1 (2001).  10.1016/S0370-1573(00)00131-9.  arXiv:astro-ph/0012122
  * (66) M.G. Alford, A. Haber, S.P. Harris, Z. Zhang, Beta Equilibrium Under Neutron Star Merger Conditions.  Universe 7(11), 399 (2021).  10.3390/universe7110399.  arXiv:2108.03324 [nucl-th]
  * (67) E.R. Most, A. Haber, S.P. Harris, Z. Zhang, M.G. Alford, J. Noronha, Emergence of microphysical viscosity in binary neutron star post-merger dynamics (2022).  arXiv:2207.00442 [astro-ph.HE]
  * (68) J.D. Salmonson, J.R. Wilson, General Relativistic Augmentation of Neutrino Pair Annihilation Energy Deposition nea r Neutron Stars.  Astroph.J. 517, 859–865 (1999).  10.1086/307232.  astro-ph/9908017
  * (69) M. Ruffert, H.T. Janka, G. Schaefer, Coalescing neutron stars - a step towards physical models. I. Hydrodynamic evolution and gravitational-wave emission.  Astron.Astroph. 311, 532–566 (1996).  astro-ph/9509006
  * (70) T. Wang, A. Burrows, Generalized Kompaneets formalism for inelastic neutrino-nucleon scattering in supernova simulations.  Phys. Rev. D 102(2), 023,017 (2020).  10.1103/PhysRevD.102.023017.  arXiv:2006.12240 [astro-ph.HE]
  * (71) A. Vlasenko, G.M. Fuller, V. Cirigliano, Neutrino Quantum Kinetics.  Phys. Rev. D 89(10), 105,004 (2014).  10.1103/PhysRevD.89.105004.  arXiv:1309.2628 [hep-ph]
  * (72) O.L. Caballero, A.C. Malkus, G.C. McLaughlin, R.A. Surman, The influence of neutrinos on the nucleosynthesis of accretion disc outflows.  Journal of Physics G Nuclear Physics 41(4), 044004 (2014).  10.1088/0954-3899/41/4/044004
  * (73) Y.L. Zhu, A. Perego, G.C. McLaughlin, Matter Neutrino Resonance Transitions above a Neutron Star Merger Remnant.  Phys. Rev. D 94(10), 105,006 (2016).  10.1103/PhysRevD.94.105006.  arXiv:1607.04671 [hep-ph]
  * (74) A. Banerjee, A. Dighe, G. Raffelt, Linearized flavor-stability analysis of dense neutrino streams.  Phys. Rev. D 84, 053,013 (2011).  10.1103/PhysRevD.84.053013.  URL https://link.aps.org/doi/10.1103/PhysRevD.84.053013
  * (75) M.R. Wu, I. Tamborra, O. Just, H.T. Janka, Imprints of neutrino-pair flavor conversions on nucleosynthesis in ejecta from neutron-star merger remnants.  Phys. Rev. D 96(12), 123,015 (2017).  10.1103/PhysRevD.96.123015.  arXiv:1711.00477 [astro-ph.HE]
  * (76) E. Grohs, S. Richers, S.M. Couch, F. Foucart, J.P. Kneller, G.C. McLaughlin, Neutrino Fast Flavor Instability in three dimensions for a Neutron Star Merger (2022).  arXiv:2207.02214 [hep-ph]
  * (77) X. Li, D.M. Siegel, Neutrino Fast Flavor Conversions in Neutron-Star Postmerger Accretion Disks.  Phys. Rev. Lett. 126(25), 251,101 (2021).  10.1103/PhysRevLett.126.251101.  arXiv:2103.02616 [astro-ph.HE]
  * (78) R. Fernández, S. Richers, N. Mulyk, S. Fahlman, The Fast Flavor Instability in Hypermassive Neutron Star Disk Outflows (2022).  arXiv:2207.10680 [astro-ph.HE]
  * (79) S. Rosswog, M. Liebendörfer, High resolution calculations of merging neutron stars. 2: Neutrino emission.  Mon.Not.Roy.Astron.Soc. 342, 673 (2003).  10.1046/j.1365-8711.2003.06579.x.  arXiv:astro-ph/0302301 [astro-ph]
  * (80) A. Perego, R. Cabezón, R. Käppeli, An advanced leakage scheme for neutrino treatment in astrophysical simulations.  Astrophys. J. Suppl. 223(2), 22 (2016).  10.3847/0067-0049/223/2/22.  arXiv:1511.08519 [astro-ph.IM]
  * (81) E. O’Connor, C.D. Ott, A New Open-Source Code for Spherically-Symmetric Stellar Collapse to Neutron Stars and Black Holes.  Class. Quant. Grav. 27, 114,103 (2010).  10.1088/0264-9381/27/11/114103.  arXiv:0912.2393 [astro-ph.HE]
  * (82) R. Ardevol-Pulpillo, H.T. Janka, O. Just, A. Bauswein, Improved Leakage-Equilibration-Absorption Scheme (ILEAS) for Neutrino Physics in Compact Object Mergers.  Mon. Not. Roy. Astron. Soc. 485(4), 4754–4789 (2019).  10.1093/mnras/stz613.  arXiv:1808.00006 [astro-ph.HE]
  * (83) Y. Sekiguchi, An implementation of the microphysics in full general relativity: A general relativistic neutrino leakage scheme.  Class. Quant. Grav. 27, 114,107 (2010).  10.1088/0264-9381/27/11/114107.  arXiv:1009.3358 [astro-ph.HE]
  * (84) M.B. Deaton, M.D. Duez, F. Foucart, E. O’Connor, C.D. Ott, L.E. Kidder, C.D. Muhlberger, M.A. Scheel, B. Szilagyi, Black Hole-Neutron Star Mergers with a Hot Nuclear Equation of State: Outflow and Neutrino-Cooled Disk for a Low-Mass, High-Spin Case.  Astrophys. J. 776, 47 (2013).  10.1088/0004-637X/776/1/47.  arXiv:1304.3384 [astro-ph.HE]
  * (85) F. Galeazzi, W. Kastaun, L. Rezzolla, J.A. Font, Implementation of a simplified approach to radiative transfer in general relativity.  Phys. Rev. D 88, 064,009 (2013).  10.1103/PhysRevD.88.064009.  arXiv:1306.4953 [gr-qc]
  * (86) D. Neilsen, S.L. Liebling, M. Anderson, L. Lehner, E. O’Connor, C. Palenzuela, Magnetized Neutron Stars With Realistic Equations of State and Neutrino Cooling.  Phys. Rev. D 89(10), 104,029 (2014).  10.1103/PhysRevD.89.104029.  arXiv:1403.3680 [gr-qc]
  * (87) F. Foucart, M.B. Deaton, M.D. Duez, E. O’Connor, C.D. Ott, R. Haas, L.E. Kidder, H.P. Pfeiffer, M.A. Scheel, B. Szilagyi, Neutron star-black hole mergers with a nuclear equation of state and neutrino cooling: Dependence in the binary parameters.  Phys. Rev. D 90, 024,026 (2014).  10.1103/PhysRevD.90.024026.  arXiv:1405.1121 [astro-ph.HE]
  * (88) D. Radice, A. Perego, K. Hotokezaka, S.A. Fromm, S. Bernuzzi, L.F. Roberts, Binary Neutron Star Mergers: Mass Ejection, Electromagnetic Counterparts and Nucleosynthesis.  Astrophys. J. 869(2), 130 (2018).  10.3847/1538-4357/aaf054.  arXiv:1809.11161 [astro-ph.HE]
  * (89) P. Mösta, D. Radice, R. Haas, E. Schnetter, S. Bernuzzi, A magnetar engine for short GRBs and kilonovae.  Astrophys. J. Lett. 901, L37 (2020).  10.3847/2041-8213/abb6ef.  arXiv:2003.06043 [astro-ph.HE]
  * (90) F. Cipolletta, J.V. Kalinani, E. Giangrandi, B. Giacomazzo, R. Ciolfi, L. Sala, B. Giudici, Spritz: General Relativistic Magnetohydrodynamics with Neutrinos.  Class. Quant. Grav. 38(8), 085,021 (2021).  10.1088/1361-6382/abebb7.  arXiv:2012.10174 [astro-ph.HE]
  * (91) E.R. Most, C.A. Raithel, Impact of the nuclear symmetry energy on the post-merger phase of a binary neutron star coalescence.  Phys. Rev. D 104(12), 124,012 (2021).  10.1103/PhysRevD.104.124012.  arXiv:2107.06804 [astro-ph.HE]
  * (92) A. Perego, E. Gafton, R. Cabezón, S. Rosswog, M. Liebendörfer, MODA: a new algorithm to compute optical depths in multidimensional hydrodynamic simulations.  Astron. Astrophys. 568, A11 (2014).  10.1051/0004-6361/201423755.  arXiv:1403.1297 [astro-ph.HE]
  * (93) D. Gizzi, C. Lundman, E. O’Connor, S. Rosswog, A. Perego, Calibration of the Advanced Spectral Leakage scheme for neutron star merger simulations, and extension to smoothed-particle hydrodynamics.  Mon. Not. Roy. Astron. Soc. 505(2), 2575–2593 (2021).  10.1093/mnras/stab1432.  arXiv:2102.08882 [astro-ph.HE]
  * (94) S. Rosswog, P. Diener, F. Torsello, Thinking Outside the Box: Numerical Relativity with Particles.  Symmetry 14(6), 1280 (2022).  10.3390/sym14061280.  arXiv:2205.08130 [gr-qc]
  * (95) S. Wanajo, Y. Sekiguchi, N. Nishimura, K. Kiuchi, K. Kyutoku, M. Shibata, Production of All the r-process Nuclides in the Dynamical Ejecta of Neutron Star Mergers.  Astroph.J.Lett. 789(2), L39 (2014).  10.1088/2041-8205/789/2/L39.  arXiv:1402.7317 [astro-ph.SR]
  * (96) F. Foucart, E. O’Connor, L. Roberts, L.E. Kidder, H.P. Pfeiffer, M.A. Scheel, Impact of an improved neutrino energy estimate on outflows in neutron star merger simulations.  Phys. Rev. D 94(12), 123,016 (2016).  10.1103/PhysRevD.94.123016.  arXiv:1607.07450 [astro-ph.HE]
  * (97) A. Perego, S. Rosswog, R.M. Cabezón, O. Korobkin, R. Käppeli, A. Arcones, M. Liebendörfer, Neutrino-driven winds from neutron star merger remnants.  Mon. Not. Roy. Astron. Soc. 443(4), 3134–3156 (2014).  10.1093/mnras/stu1352.  arXiv:1405.6730 [astro-ph.HE]
  * (98) D. Gizzi, E. O’Connor, S. Rosswog, A. Perego, R. Cabezón, L. Nativi, A multidimensional implementation of the Advanced Spectral neutrino Leakage scheme.  Mon. Not. Roy. Astron. Soc. 490(3), 4211–4229 (2019).  10.1093/mnras/stz2911.  arXiv:1906.11494 [astro-ph.HE]
  * (99) K.S. Thorne, Relativistic radiative transfer - Moment formalisms.  Mon.Not.Roy.Astro.Soc. 194, 439–473 (1981)
  * (100) M. Shibata, K. Kiuchi, Y. Sekiguchi, Y. Suwa, Truncated Moment Formalism for Radiation Hydrodynamics in Numerical Relativity.  Progress of Theoretical Physics 125, 1255 (2011).  arXiv:1104.3937 [astro-ph.HE]
  * (101) O. Just, A. Bauswein, R. Ardevol Pulpillo, S. Goriely, H.T. Janka, Comprehensive nucleosynthesis analysis for ejecta of compact binary mergers.  Mon.Not.Roy.Astron.Soc. 448(1), 541–567 (2015).  10.1093/mnras/stv009.  arXiv:1406.2687 [astro-ph.SR]
  * (102) T. Kuroda, T. Takiwaki, K. Kotake, A New Multi-Energy Neutrino Radiation-Hydrodynamics Code in Full General Relativity and Its Application to Gravitational Collapse of Massive Stars.  Astrophys. J. Suppl. 222(2), 20 (2016).  10.3847/0067-0049/222/2/20.  arXiv:1501.06330 [astro-ph.HE]
  * (103) L.F. Roberts, C.D. Ott, R. Haas, E.P. O’Connor, P. Diener, E. Schnetter, General Relativistic Three-Dimensional Multi-Group Neutrino Radiation-Hydrodynamics Simulations of Core-Collapse Supernovae.  Astrophys. J. 831, 98 (2016).  10.3847/0004-637X/831/1/98.  arXiv:1604.07848 [astro-ph.HE]
  * (104) F. Foucart, E. O’Connor, L. Roberts, M.D. Duez, R. Haas, L.E. Kidder, C.D. Ott, H.P. Pfeiffer, M.A. Scheel, B. Szilágyi, Post-merger evolution of a neutron star-black hole binary with neutrino transport.  Phys.Rev.D 91(12), 124021 (2015).  10.1103/PhysRevD.91.124021.  arXiv:1502.04146 [astro-ph.HE]
  * (105) C.Y. Cardall, E. Endeve, A. Mezzacappa, Conservative $3\mathbf{+}1$ general relativistic variable eddington tensor radiation transport equations.  Phys. Rev. D 87, 103,004 (2013).  10.1103/PhysRevD.87.103004.  URL https://link.aps.org/doi/10.1103/PhysRevD.87.103004
  * (106) D. Radice, F. Galeazzi, J. Lippuner, L.F. Roberts, C.D. Ott, L. Rezzolla, Dynamical Mass Ejection from Binary Neutron Star Mergers.  Mon. Not. Roy. Astron. Soc. 460(3), 3255–3271 (2016).  10.1093/mnras/stw1227.  arXiv:1601.02426 [astro-ph.HE]
  * (107) C.D. Levermore, G.C. Pomraning, A flux-limited diffusion theory.  Astroph.J. 248, 321–334 (1981).  10.1086/159157
  * (108) G.N. Minerbo, Maximum entropy Eddington factors.  J.Quant.Spec.Radiat.Transf. 20, 541–545 (1978).  10.1016/0022-4073(78)90024-9
  * (109) F. Foucart, Monte Carlo closure for moment-based transport schemes in general relativistic radiation hydrodynamic simulations.  Mon. Not. Roy. Astron. Soc. 475(3), 4186–4207 (2018).  10.1093/mnras/sty108.  arXiv:1708.08452 [astro-ph.HE]
  * (110) D. Radice, S. Bernuzzi, A. Perego, R. Haas, A new moment-based general-relativistic neutrino-radiation transport code: Methods and first applications to neutron star mergers.  Mon. Not. Roy. Astron. Soc. 512(1), 1499–1521 (2022).  10.1093/mnras/stac589.  arXiv:2111.14858 [astro-ph.HE]
  * (111) B. van Leer, Towards the Ultimate Conservative Difference Scheme. IV. A New Approach to Numerical Convection.  J. of Comp.Phys. 23, 276 (1977).  10.1016/0021-9991(77)90095-X
  * (112) G.S. Jiang, C.W. Shu, Efficient implementation of weighted eno schemes.  J. of Comp.Phys. 126(1), 202 – 228 (1996).  10.1006/jcph.1996.0130
  * (113) R. Borges, M. Carmona, B. Costa, W.S. Don, An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws.  J. of Comp.Phys. 227(6), 3191 – 3211 (2008).  10.1016/j.jcp.2007.11.038
  * (114) D. Radice, L. Rezzolla, THC: a new high-order finite-difference high-resolution shock-capturing code for special-relativistic hydrodynamics.  Astron.Astrophys. 547, A26 (2012).  arXiv:1206.6502 [astro-ph.IM]
  * (115) L. Sun, M. Ruiz, S.L. Shapiro, A. Tsokaros, Jet launching from binary neutron star mergers: Incorporating neutrino transport and magnetic fields.  Phys. Rev. D 105(10), 104,028 (2022).  10.1103/PhysRevD.105.104028.  arXiv:2202.12901 [astro-ph.HE]
  * (116) E. Audit, P. Charrier, J.P. Chieze, B. Dubroca, A radiation-hydrodynamics scheme valid from the transport to the diffusion limit (2002).  arXiv:astro-ph/0206281
  * (117) F. Foucart, M.D. Duez, F. Hebert, L.E. Kidder, H.P. Pfeiffer, M.A. Scheel, Monte-Carlo neutrino transport in neutron star merger simulations.  Astrophys. J. Lett. 902, L27 (2020).  10.3847/2041-8213/abbb87.  arXiv:2008.08089 [astro-ph.HE]
  * (118) J.M. Miller, B.R. Ryan, J.C. Dolence, A. Burrows, C.J. Fontes, C.L. Fryer, O. Korobkin, J. Lippuner, M.R. Mumpower, R.T. Wollaeger, Full Transport Model of GW170817-Like Disk Produces a Blue Kilonova.  Phys. Rev. D 100(2), 023,008 (2019).  10.1103/PhysRevD.100.023008.  arXiv:1905.07477 [astro-ph.HE]
  * (119) S. Richers, D. Kasen, E. O’Connor, R. Fernández, C.D. Ott, Monte Carlo Neutrino Transport through Remnant Disks from Neutron Star Mergers.  Astroph.J. 813, 38 (2015).  10.1088/0004-637X/813/1/38.  arXiv:1507.03606 [astro-ph.HE]
  * (120) F. Foucart, M.D. Duez, F. Hebert, L.E. Kidder, P. Kovarik, H.P. Pfeiffer, M.A. Scheel, Implementation of Monte Carlo Transport in the General Relativistic SpEC Code.  Astrophys. J. 920(2), 82 (2021).  10.3847/1538-4357/ac1737.  arXiv:2103.16588 [astro-ph.HE]
  * (121) S.A. Hughes, C.R. Keeton, II, P. Walker, K.T. Walsh, S.L. Shapiro, S.A. Teukolsky, Finding black holes in numerical spacetimes.  Phys.Rev.D 49, 4004–4015 (1994).  10.1103/PhysRevD.49.4004
  * (122) B.R. Ryan, J.C. Dolence, C.F. Gammie, bhlight: General Relativistic Radiation Magnetohydrodynamics with Monte Carlo Transport.  Astroph.J. 807, 31 (2015).  10.1088/0004-637X/807/1/31.  arXiv:1505.05119 [astro-ph.HE]
  * (123) J.M. Miller, B.R. Ryan, J.C. Dolence, $\nu\texttt{bhlight}$: Radiation GRMHD for Neutrino-Driven Accretion Flows.  Astrophys. J. Suppl. 241(2), 30 (2019).  10.3847/1538-4365/ab09fc.  arXiv:1903.09273 [astro-ph.IM]
  * (124) S. Richers, H. Nagakura, C.D. Ott, J. Dolence, K. Sumiyoshi, S. Yamada, A Detailed Comparison of Multidimensional Boltzmann Neutrino Transport Methods in Core-collapse Supernovae.  Astrophys. J. 847(2), 133 (2017).  10.3847/1538-4357/aa8bb2.  arXiv:1706.06187 [astro-ph.HE]
  * (125) J. Fleck, J. A., J.D. Cummings, An implicit Monte Carlo scheme for calculating time and frequency dependent nonlinear radiation transport.  Journal of Computational Physics 8(3), 313–342 (1971).  10.1016/0021-9991(71)90015-5
  * (126) Y. Sekiguchi, K. Kiuchi, K. Kyutoku, M. Shibata, K. Taniguchi, Dynamical mass ejection from the merger of asymmetric binary neutron stars: Radiation-hydrodynamics study in general relativity.  Phys. Rev. D 93(12), 124,046 (2016).  10.1103/PhysRevD.93.124046.  arXiv:1603.01918 [astro-ph.HE]
  * (127) R. Fernández, F. Foucart, J. Lippuner, The landscape of disc outflows from black hole–neutron star mergers.  Mon. Not. Roy. Astron. Soc. 497(3), 3221–3233 (2020).  10.1093/mnras/staa2209.  arXiv:2005.14208 [astro-ph.HE]
  * (128) M. Shibata, S. Fujibayashi, Y. Sekiguchi, Long-term evolution of neutron-star merger remnants in general relativistic resistive magnetohydrodynamics with a mean-field dynamo term.  Phys. Rev. D 104(6), 063,026 (2021).  10.1103/PhysRevD.104.063026.  arXiv:2109.08732 [astro-ph.HE]
  * (129) S. Fujibayashi, K. Kiuchi, S. Wanajo, K. Kyutoku, Y. Sekiguchi, M. Shibata, Comprehensive study on the mass ejection and nucleosynthesis in the binary neutron star mergers leaving short-lived massive neutron stars (2022).  arXiv:2205.05557 [astro-ph.HE]
  * (130) S. Fujibayashi, S. Wanajo, K. Kiuchi, K. Kyutoku, Y. Sekiguchi, M. Shibata, Postmerger Mass Ejection of Low-mass Binary Neutron Stars.  Astrophys. J. 901(2), 122 (2020).  10.3847/1538-4357/abafc2.  arXiv:2007.00474 [astro-ph.HE]
  * (131) A.M. Beloborodov, Nuclear Composition of Gamma-Ray Burst Fireballs.  Astroph.J 588(2), 931–944 (2003).  10.1086/374217.  arXiv:astro-ph/0210522 [astro-ph]
  * (132) D.M. Siegel, B.D. Metzger, Three-dimensional GRMHD Simulations of Neutrino-cooled Accretion Disks from Neutron Star Mergers.  Astroph.J. 858(1), 52 (2018).  10.3847/1538-4357/aabaec.  arXiv:1711.00868 [astro-ph.HE]
  * (133) O. Just, S. Goriely, H.T. Janka, S. Nagataki, A. Bauswein, Neutrino absorption and other physics dependencies in neutrino-cooled black hole accretion discs.  Mon. Not. Roy. Astron. Soc. 509(1), 1377–1412 (2021).  10.1093/mnras/stab2861.  arXiv:2102.08387 [astro-ph.HE]
  * (134) R. Fernández, F. Foucart, D. Kasen, J. Lippuner, D. Desai, L.F. Roberts, Dynamics, nucleosynthesis, and kilonova signature of black hole—neutron star merger ejecta.  Class. Quant. Grav. 34(15), 154,001 (2017).  10.1088/1361-6382/aa7a77.  arXiv:1612.04829 [astro-ph.HE]
  * (135) S. de Haas, P. Bosch, P. Mösta, S. Curtis, N. Schut, Magnetic field effects on nucleosynthesis and kilonovae from neutron star merger remnants (2022).  arXiv:2208.05330 [astro-ph.HE]
  * (136) R. Popham, S.E. Woosley, C. Fryer, Hyperaccreting Black Holes and Gamma-Ray Bursts.  Astroph.J. 518(1), 356–374 (1999).  10.1086/307259.  arXiv:astro-ph/9807028 [astro-ph]
  * (137) J.P. Kneller, G.C. McLaughlin, R.A. Surman, Neutrino scattering, absorption and annihilation above the accretion discs of gamma ray bursts.  Journal of Physics G Nuclear Physics 32(4), 443–462 (2006).  10.1088/0954-3899/32/4/004.  arXiv:astro-ph/0410397 [astro-ph]
  * (138) I. Zalamea, A.M. Beloborodov, Efficiency of Neutrino Annihilation around Spinning Black Holes.  AIP Conf. Proc. 1133(1), 121–123 (2009).  10.1063/1.3155863.  arXiv:0812.4041 [astro-ph]
  * (139) A. Perego, H. Yasin, A. Arcones, Neutrino pair annihilation above merger remnants: implications of a long-lived massive neutron star.  J. Phys. G 44(8), 084,007 (2017).  10.1088/1361-6471/aa7bdc.  arXiv:1701.02017 [astro-ph.HE]
  * (140) F. Foucart, R. Haas, M.D. Duez, E. O’Connor, C.D. Ott, L. Roberts, L.E. Kidder, J. Lippuner, H.P. Pfeiffer, M.A. Scheel, Low mass binary neutron star mergers : gravitational waves and neutrino emission.  Phys. Rev. D 93(4), 044,019 (2016).  10.1103/PhysRevD.93.044019.  arXiv:1510.06398 [astro-ph.HE]
  * (141) V. Nedora, F. Schianchi, S. Bernuzzi, D. Radice, B. Daszuta, A. Endrizzi, A. Perego, A. Prakash, F. Zappa, Mapping dynamical ejecta and disk masses from numerical relativity simulations of neutron star mergers.  Class. Quant. Grav. 39(1), 015,008 (2022).  10.1088/1361-6382/ac35a8.  arXiv:2011.11110 [astro-ph.HE]
  * (142) A. Henkel, F. Foucart, G. Raaijmakers, S. Nissanke, A study of the agreement between binary neutron star ejecta models derived from numerical relativity simulations (2022).  arXiv:2207.07658 [astro-ph.HE]