Neutrino transport in general relativistic neutron star merger simulations

Numerical simulations of neutron star–neutron star and neutron star–black hole binaries play an important role in our ability to model gravitational-wave and electromagnetic signals powered by these systems. These simulations have to take into account a wide range of physical processes including general relativity, magnetohydrodynamics, and neutrino radiation transport. The latter is particularly important in order to understand the properties of the matter ejected by many mergers, the optical/infrared signals powered by nuclear reactions in the ejecta, and the contribution of that ejecta to astrophysical nucleosynthesis. However, accurate evolutions of the neutrino transport equations that include all relevant physical processes remain beyond our current reach. In this review, I will discuss the current state of neutrino modeling in general relativistic simulations of neutron star mergers and of their post-merger remnants. I will focus on the three main types of algorithms used in simulations so far: leakage, moments, and Monte-Carlo scheme. I will review the advantages and limitations of each scheme, as well as the various neutrino–matter interactions that should be included in simulations. We will see that the quality of the treatment of neutrinos in merger simulations has greatly increased over the last decade, but also that many potentially important interactions remain difficult to take into account in simulations (pair annihilation, oscillations, inelastic scattering).


Introduction
Over the last decade, the study of merging compact objects has made tremendous progress. Recently observed astrophysical events provide us with some of the most reliable information currently at our disposal regarding the population of stellar mass black holes in the nearby Universe. Rarer events that include neutron stars also inform us about the mass distribution of neutron stars, the equation of state of dense matter, and the origin of heavy elements formed through rapid neutron capture nucleosynthesis (r-process). Our ability to study these systems has largely grown in tandem with the sensitivity of the LIGO and Virgo gravitational wave detectors. Gravitational wave observatories have now detected dozens of binary black hole (BBH) mergers, as well as two likely binary neutron star (BNS) mergers and at least two likely neutron star-black hole (NSBH) mergers (see Sect. 2.2 for a more detailed discussion of these events). An overview of these events can be found in the three GWTC catalogues (Abbott et al, 2019(Abbott et al, , 2021a. While BNS and NSBH mergers are not as commonly observed as BBH mergers, they do have important advantages for nuclear astrophysics. The presence of a neutron star means that these systems can potentially be used to constrain the equation of state of cold, neutron rich dense matter (Abbott et al, 2018) -a crucial source of information about many-nucleon interactions and, potentially, the high-density states of quantum chromodynamics. Additionally, some mergers and post-merger remnants eject material that undergoes r-process nucleosynthesis. The radioactive decay of the ashes of the r-process can then power optical/infrared emission days to weeks after the merger: a kilonova (Lattimer and Schramm, 1976;Li and Paczynski, 1998;Metzger et al, 2010;Roberts et al, 2011;Kasen et al, 2013). The production site(s) of r-process elements remain(s) very uncertain today, and the observation of neutron star mergers and associated kilonovae may help us solve the long-standing problem of their astrophysical origin. Additionally, some post-merger remnants likely produce collimated relativistic outflows (jets) that are currently believed to be the source of short-hard gamma-ray bursts (SGRBs) (Eichler et al, 1989;Nakar, 2007;Fong and Berger, 2013). The exact process powering SGRBs is however not well understood, and further observations of neutron star mergers could help us ellucidate how these high-energy events occur in practice. Finally, joint observations of neutron star mergers using both gravitational and electromagnetic waves may also provide additional information about the properties of the merging compact objects, the position of the merging binary, and even the value of the Hubble constant (Holz and Hughes, 2005;Nissanke et al, 2010;Abbott et al, 2017a;Hotokezaka et al, 2019).
Neutron star mergers involve a wide range of nonlinear physical processes, preventing us from providing quantitative theoretical predictions for the result of a merger using purely analytical methods. As a result, numerical simulations are an important tool in current attempts to model the gravitational wave and electromagnetic signals powered by compact binary mergers. Gravity, fluid dynamics, magnetic fields and neutrinos all play major roles during and after neutron star mergers, with out-of-equlibrium nuclear reactions also becoming important on longer time scales (∼ seconds). In theory, merger simulations thus need to solve Boltzmann's equations of radiation transport coupled to the relativistic equations of magnetohydrodynamics and Einstein's equation of general relativity. However, no simulation can do this with the desired level of realism at this point. Two major roadblocks to this modeling efforts are our inability to properly resolve magnetohydrodynamical instabilities during merger (and thus the dynamo process that may follow the growth of magnetic fields due to these instabilities) (Kiuchi et al, 2015), as well as the difficulty of properly solving Boltzmann's equation of radiation transport for the evolution of neutrinos . In this review, we focus on the second problem. The role of magnetic fields in merger simulations is discussed in more detail, for example, in Baiotti and Rezzolla (2017); Paschalidis (2017); Burns (2020).
Neutrinos play a number of roles in neutron star mergers, with particularly noticeable impacts on the production of r-process elements and the properties of kilonovae. However, properly accounting for neutrino-matter interactions in neutron star mergers remains a difficult problem because, within a merger remnant, neutrinos transition from being in equilibrium with the fluid (in dense hot regions) to mostly free-streaming through the ejected material (far away). In the intermediate regions, neutrino-matter interactions play an important role in the evolution of the temperature and composition of the fluid, but neutrinos cannot be assumed to be in equilibrium with the fluid. Numerical methods that properly capture both regimes are technically challenging and/or computationally expensive. As a result, most merger simulations use approximate neutrino transport algorithms that introduce potentially significant and often hard to quantify errors in our predictions for the nuclei produced during r-process nucleosynthesis and for the properties of kilonovae.
The main objective of this review is to provide an overview of the various algorithms currently used in general relativistic simulations of neutron star mergers and of their post-merger remnants. These can be broadly classified into three groups: leakage methods, which do not explicitly transport neutrinos; moment schemes, which evolve a truncated expansion of the transport equations in momentum space with methods highly similar to those used to evolve the equations of relativistic magnetohydrodynamics; and Monte-Carlo methods, which sample the distribution of neutrinos with packets (or superparticles) propagating through numerical simulations. These are discussed in detail in Sect. 4. Section 2 aims to provide some scientific background about merging neutron stars, while Sect. 3 provides an overview of neutrino physics in neutron star mergers, and of the important neutrino-matter interactions that are currently included or neglected in simulations. Finally, Sect. 5 discusses what existing simulations can tell us about the ways in which our choice of algorithm impacts our numerical results. We note that the objective here is not to review all results in the study of neutron star mergers with neutrinos, but rather to focus on the numerical methods used to perform general relativistic radiation transport. We will thus focus on comparisons of different numerical methods, rather that provide an extensive review of existing simulations that make use of neutrino transport.
Conventions: In this manuscript, latin letters are used for the indices of spatial 3-dimensional vectors/tensors, while greek letters are used for the indices of 4-dimensional vectors / tensors. Sections discussing numerical methods will often use units such that h = c = G = 1, but we explicitly keep physical constants in our expressions when discussing interaction rates.
2 Scientific background 2.1 Overview of neutron star mergers physics Before delving deeper into the topic of radiation transport in neutron star mergers, it is worth reviewing how we currently understand the evolution of these systems, as well as when different physical processes are expected to play an important role. When discussing neutron star merger simulations, we are typically concerned with the evolution of a binary from tens of milliseconds before merger to a few seconds after merger, i.e., from the moment standard post-Newtonian methods can no longer accurately model the gravitational wave signal to the moment when the accretion disk formed during a merger has lost most of its mass. In the late inspiral (O(10) orbits before merger), the tidal distortion of a neutron star by its binary companion has a potentially measurable impact on the gravitational wave signal, which can be used to put constraints on the equation of state of neutron stars (Flanagan and Hinderer, 2008;Abbott et al, 2018). The main role of numerical simulations in that regime is to help test and calibrate analytical waveform models used in the analysis of gravitational wave events (e.g., Bernuzzi Thompson et al (2020); Matas et al (2020) for NSBH mergers). General relativity, fluid dynamics, and the choice of equation of state are important at that stage, but magnetic fields only impact relatively weak pre-merger electromagnetic signals and neutrinos have practically no impact on the evolution of the system.
For NSBH binaries, the same remains true during the merger itself, i.e., the few milliseconds during which the neutron star is either tidally disrupted by its black hole companion, or absorbed whole by the black hole. The outcome of the merger is determined by the masses and spins of the compact objects, the equation of state of dense matter (Lattimer and Schramm, 1976;Pannarale et al, 2011;Foucart, 2012), and the eccentricity of the orbit (East et al, 2015). Numerical simulations of low-eccentricity binaries have shown that only low mass and/or high spin black holes disrupt their neutron star companions (M BH 5 M for non-spinning compact objects and circular orbits), a prerequisite to the production of any post-merger electromagnetic signal. If the neutron star is tidally disrupted, a few percents of a solar mass of very neutron rich, cold matter is typically ejected, and tenths of a solar mass remain in a bound accretion disk and/or tidal tail around the black hole (see e.g., Foucart In disrupting NSBH systems, most of the matter is rapidly accreted onto the black hole, while the rest forms an accretion disk and extended tidal tail. Low-mass NSNS binaries form a massive neutron star remnant surrounded by a bound disk, with a smaller amount of material ejected in the tidal tail. The right panel is reproduced from Foucart et al (2017); the left panel visualizes a simulation from Foucart et al (2016a). 2020; Kyutoku et al 2021 for recent reviews, and Fig. 1). In eccentric binaries, neutron stars are typically easier to disrupt, and eject more mass in their tidal tails.
For BNS systems, on the other hand, other physical processes become important once the neutron stars collide. First, the shear region that is naturally created between the merging neutron stars is unstable to the Kelvin-Helmoltz instability, leading to the rapid growth of small scale turbulence (Kiuchi et al, 2015). Magnetic fields are quickly amplified to B ∼ 10 16 G as a result, and start to play an important role in the evolution of the system. Whether a dynamo process can generate a large scale magnetic field from this turbulent state is an important open questions that simulations have not so far been able to answer. The collision of the two neutron stars also creates hot regions where neutrino emission and absorption can no longer be ignored. BNS mergers eject relatively small amounts of cold tidal ejecta ( 0.01 M ), as well as hotter material coming from the regions where the cores of the neutron stars collide. We will see that neutrinos play an important role in the evolution of that hot ejecta. Depending on the equation of state and on the mass of the system, the remnant may immediately collapse to a black hole (on milliseconds time scales), remain temporarily supported by rotation and/or thermal pressure, or form a long-lived neutron star (as on Fig. 1). In all cases, that remnant is surrounded by a hot accretion disk -with more asymmetric systems producing more massive disks (see e.g., Baiotti and Rezzolla 2017;Burns 2020; Radice et al 2020 for recent reviews).
After merger, neutrino emission is the main source of cooling for the accretion disk and remnant neutron star (if there is one), and neutrino-matter interactions drive changes in the composition of the disk material and of the outflows. Initially, the efficiency of neutrinos in cooling the disk lies in between the radiatively efficient (thin disks) and radiatively inefficient (thick disks) regimes observed in AGNs. NSBH and BNS simulations including radiation transport show a disk aspect ratio H/R ∼ (0.2 − 0.3) (with H the scale height of the disk and R its radius) (Foucart et al, 2015;Fujibayashi et al, 2018). Hydrodynamical shocks and/or fluid instabilities and then turbulence driven by the magnetorotational instability (MRI) lead to angular momentum transport and heating in the disk, and drive accretion onto the compact object. If a large scale poloidal magnetic field threads the disk, magnetically driven outflows are likely to unbind ∼ 20% of the mass of the disk (Siegel and Metzger, 2017;Fernández et al, 2019) -but this is not a given considering uncertainties about the large scale structure of the magnetic field in post-merger accretion disks. Indeed, while it is possible to grow such a large scale field after merger (Christie et al, 2019), this takes too long to efficiently contribute to the production of winds. A large scale magnetic field generated during merger appears to be required for these winds to exist.
After O(100 ms), the density of the disk decreases enough that neutrino cooling becomes inneficient (Fernández and Metzger, 2013;De and Siegel, 2021), while the MRI remains active. The disk becomes advection dominated. It puffs up to H/R ∼ 1, and viscous spreading of the disk leads to the ejection of 5% − 25% of the disk mass (viscous outflows) (Fernández and Metzger, 2013). Neutrino-matter interactions directly impact the properties of magnetically driven outflows, and indirectly impact the properties of viscous outflows (due to neutrino-matter interactions during the early evolution of the disk, before weak-interaction freeze-out).
The post-merger evolution is also impacted by the presence and life time of a massive neutron star remnant. A hot neutron star remnant is a bright source of neutrinos that can accelerate changes to the composition of matter outflows in the polar regions. How efficiently matter can accrete onto the neutron star remains uncertain. Axisymmetric simulations treating the neutron star surface as a hard boundary predict the eventual ejection of most of the remnant disk (Metzger and Fernández, 2014); whether this would remain true for more realistic boundary conditions is unclear, but it is at least likely that a larger fraction of the disk is eventually unbound for neutron star remnants than for black hole remnants. The neutron star remnants themselves are initially differentially rotating, and simulations generally find rotation profiles that are stable to the MRI in most of the star (the angular velocity increases with radius). Some other angular momentum transport mechanism is thus required to bring these remnants to uniform rotation, e.g., convection and/or the Spruit-Taylor dynamo (Margalit et al, 2022). The exact impact of the interaction between the neutron star remnant, its external magnetic field, and the surrounding accretion disk on the evolution of the system remains very uncertain. Examples of post-merger remnants are shown in Fig. 2

Observables and existing observations
The main signals observed so far in neutron star mergers include gravitational wave emission during the late inspiral of the binary towards mergers, SGRBs Fig. 2 Post-merger remnant a few milliseconds after a BNS merger (Left), and 0.3 s after a NSBH merger (Right). The BNS system forms a massive, differentially rotating neutron star surrounded by a low-mass accretion disk, with shocked spiral arms visible in the disk. The NSBH system forms an extended accretion disk around the remnant black hole, with collimated magnetic fields in the polar region. The right panel is reproduced from ; the left panel visualizes a simulation from Foucart et al (2016a).
(and their multi-wavelength afterglows) likely due to relativistic jets powered by the post-merger remnant, and kilonovae. For a system with component masses m 1 , m 2 , the gravitational waves provide us with a very accurate measurement of the chirp mass M c = (m 1 m 2 ) 0.6 /(m 1 + m 2 ) 0.2 , as well as, for sufficiently loud signals, less accurate information about the mass ratio (and thus the component masses), the spins of the compact objects, the equation of state of neutron stars (through their tidal deformability), as well as the distance, orientation, and sky localization of the source (especially for multidetector observations). We will not discuss the gravitational wave signal in much more detail here, as it is not meaningfully impacted by neutrinos. Outflows generated during and after the merger (see previous section) will be the main source of post-merger electromagnetic signals. Relativistic collimated outflows power SGRBs detectable by observers located along the spin axis of the remnant. As the jet material becomes less relativistic, SGRBs are followed by longer wavelength afterglows detectable by off-axis observers (Fong et al, 2015). The gamma-ray emission is very short lived ( 2 s for a typical SGRB), but radio afterglows can still be observed a year after the merger (Mooley et al, 2018). The exact mechanism powering the relativistic jet remains unknown. The most commonly discussed model requires the formation of a large scale poloidal magnetic field threading a black hole remnant, with energy extraction from the black hole's rotation though a Blandford-Znajek-like process (Blandford and Znajek, 1977). Some SGRB models are however powered by neutrino-antineutrino pair annihilations in the polar regions. Explaining the most energetic SGRBs through this mechanism is difficult given what is currently known of the neutrino luminosity of post-merger remnants and the efficiency of the pair annihilation process (Just et al, 2016), yet even in a magnetically-powered SGRB, energy deposition due to neutrino pair annihilation or baryon loading of the polar regions due to neutrino-driven winds could impact the formation of a jet (Fujibayashi et al, 2017).
The properties of kilonovae and the role of neutron star mergers in astrophysical nucleosynthesis are likely to be much more significantly impacted by neutrinos than gravitational waves or even SGRBs. Absorption and emission of electron-type neutrinos (ν e ) and antineutrinos (ν e ) modifies the relative number of neutrons and protons in the fluid. This is usually expressed through the lepton fraction with n e ± , n n , n p , n νe , nν e the number density of electrons, positrons, neutrons, protons, ν e andν e respectively. Many simulations use the net electron fraction Y e instead of the lepton fraction, and assume that charge neutrality requires n e − − n e + = n p 1 , so that Y e = n p n p + n n . ( The electron fraction is a crucial determinant of the outcome of r-process nucleosynthesis in merger outflows. Low Y e outflows (roughly Y e 0.25) produce heavier r-process elements, while higher Y e outflows produce lighter r-process elements (Lippuner and Roberts, 2015). In particular, for the conditions typically observed in merger outflows, there is not much production of elements above the "2nd peak" of the r-process (at atomic number A ∼ 130) for high-Y e outflows, and an under-production of elements below the 2nd peak for neutronrich (low Y e ) outflows. Cold outflows that do not interact much with neutrinos are typically neutron-rich, but hotter outflows can end up with Y e ∼ 0.4 − 0.5 due to neutrino-matter interactions. Cooling from neutrino emission and heating from neutrino absorption are also important to the thermodynamics of the remnant and of the outflows, and neutrino absorption in the disk corona and close to the neutron star surface can lead to the production of neutrino-driven winds (Dessart et al, 2009). It is thus clear that neutrino-matter interactions should be properly understood if we aim to model the role of neutron star mergers in the production of r-process elements.
The impact of neutrinos on kilonovae is less direct but no less important. Most of the r-process occurs within a few seconds of the merger, after which the outflows are mainly composed of radioactively unstable heavy nuclei. Radioactive decays of these nuclei will continue to release energy over much longer timescales. Initially, the outflows are opaque to most photons, and decay products are thermalized -except for neutrinos, which immediately escape the outflows. As the density of the outflows decrease, however, they will eventually become optically thin to optical/infrared photons. When this transition happens depends on the composition of the outflows. Lanthanides and actinides, which are among the heavier r-process elements that are only produced by neutron-rich outflows, have much higher opacities than other nuclei produced during the r-process. As a result, neutron-rich outflows become optically thin later than higher Y e outflows (∼ 10 days vs ∼ 1 day), and the corresponding kilonova signal is redder (peaks in the infrared, instead of in the optical). Overall, the duration, color, and magnitude of a kilonova tell us about the mass of the outflows, their composition, and their velocity . For a given binary merger, it will also depend on the relative orientation of the binary and the observer, as different types of outflows have different geometry.
Other electromagnetic counterparts to neutron star mergers have been proposed, with no confirmed observations so far. This include bursts of radiation before merger (Tsang, 2013), continuous emission from magnetosphere interactions (Palenzuela et al, 2013), coherent emission from magnetosphere interactions (Most and Philippov, 2022), and months to decades-long synchrotron radio emission from the mildly relativistic ejecta as it interacts with the interstellar medium (Hotokezaka et al, 2016). Neutrinos have no impact on the first three, however, and only a minor impact on the third (as neutrinomatter interactions may slightly change the mass / velocity of the outflows). More detailed discussions of the range of electromagnetic transients that may follow a merger can be found, e.g., in Fernández and Metzger (2016); Burns (2020) Electromagnetic emission from neutron star mergers has likely been observed for decades now in the form of SGRBs, and a first kilonova may have been observed in the afterglow of GRB130603B as early as 2013 (Tanvir et al, 2013;Berger et al, 2013). However, our current understanding of the engine powering SGRBs is not sufficient to provide us with much information about the parameters of the binary system that created the burst -or even to differentiate between a BNS and NSBH merger. Gravitational wave observations provide more direct information about the properties of the compact objects. So far, two systems have been observed with component masses most easily explained by the merger of two neutron stars: GW170817 (Abbott et al, 2017b) and GW1902425 (Abbott et al, 2020). The former is a relatively low mass system, whose observation was followed by a weak SGRB (most likely observed off-axis), radio emission most likely associated with a relativistic jet, and a clear kilonova signal most easily explain by a combination of at least two outflow components -one that led to strong r-process nucleosynthesis, and one that did not. The exact process that produced these outflows remain a subject of research today. GW190425 has a higher total mass (3.4 M ). There was no observed electromagnetic counterpart to that signal, a relatively unsurprising result considering the large uncertainty in the location of the source and the high likelihood that such a system did not eject a significant amount of matter (Barbieri et al, 2021;Raaijmakers et al, 2021;Dudi et al, 2021;Camilletti et al, 2022). At least two NSBH mergers were observed in 2020 (Abbott et al, 2021c), with more candidates also available in the latest gravitational wave catalogue (Abbott et al, 2021b). None of these systems was however expected to lead to the disruption of their neutron star, and thus their lack of electromagnetic counterpart was unsurprising.
Overall, we note that the analysis of current and future observations of neutron star mergers would benefit from accurate models of kilonova signals, as well as from an improved understanding of the engine behind gamma-ray bursts. In that respect, it is particularly important to understand the role of neutrinos in setting the composition of the outflows powering kilonovae, and possibly their impact on the production of relativistic jets. In the rest of this review, we will mainly focus on these issues, and on the methods available to evolve neutrinos in merger simulations.

Neutrinos in mergers 3.1 Definitions
When solving the general relativistic equations of radiation transport, we would ideally evolve Boltzmann's equation, or the quantum kinetics equations (QKE, when accounting for neutrino oscillations). Classically, we evolve the distribution function of neutrinos f ν (t, x i , p j ), defined such that is the number of neutrinos within a 6D volume of phase space V . Here, x i are the spatial coordinates and p j the spatial components of the 4-momentum one-form p µ , while h is Planck's constant. When using the classical equations of radiation transport, we usually neglect neutrino masses and assume p µ p µ = 0. Boltzmann's equation is then with τ the proper time in the fluid frame, ν the neutrino energy in the fluid frame, and Γ α βγ the Christoffel symbols. The left-hand side simply implies that neutrinos follow null geodesics, while the right-hand side includes all neutrino-matter and neutrino-neutrino interactions, and thus hides most of the complexity in these equations. We note that we should evolve a separate f ν for each type of neutrinos (ν e , ν µ , ν τ ) and antineutrinos (ν e ,ν µ ,ν τ ); and that these distributions functions may be coupled through the collision terms. As neutrinos are fermions, we have 0 ≤ f ν ≤ 1.
The spacial coordinate volume d 3 x = dxdydz and momentum volume d 3 p = dp x dp y dp z are not invariant under coordinate transformations, but d 3 xp t √ −g and d 3 p(p t √ −g) −1 are, with g the determinant of the spacetime metric g µν . Thus d 3 xd 3 p is invariant under coordinate transformations. The stress-energy tensor of neutrinos at (t, x i ) is In general relativistic merger simulations, we often use the 3+1 decomposition of the metric with α the lapse, β i the shift, and γ ij the 3-metric on a slice of constant time t.
The unit normal one-form to such a slice is then n µ = (−α, 0, 0, 0), and the 4vector n µ = g µν n ν can be interpreted as the 4-velocity of an observer moving along that normal -which we will call normal observer from now on. From there, we can deduce that = −p µ n µ = αp t is the energy of a neutrino of 4momentum p µ as measured by a normal observer. More generally, the energy of a neutrino measured by an observer with 4-velocity u µ is ν = −p µ u µ . Here, we will generally reserve the symbol for the energy measured by normal observers, and ν for the energy measured in the fluid rest frame, i.e., when u µ is the 4-velocity of the fluid.

Equilibrium distribution
We will often make use of the equilibrium distribution of neutrinos. For neutrinos in equilibrium with a fluid at temperature T moving with 4-velocity u µ , that is the Fermi-Dirac distribution with µ the chemical potential of neutrinos, and k B Boltzmann's constant. We note that in an orthonormal frame (t,x i ) the energy density of neutrinos is withˆ = pt the energy of neutrinos as measured by a stationary observer in the orthonormal frame. We thus see that we recover the expected results for the equilibrium energy of a fermion gas in the fluid frame, where in the last expression we used the special relativistic result ν = p c. This is more easily expressed in terms of the Fermi integrals F n , which we will use extensively in this section: From this definition, we see that Similarly, the equilibrium number density of neutrinos is and the average energy of neutrinos in equilibrium with the fluid (which asymptotes to 3.15k B T at low densities, when µ k B T ).

Commonly considered reactions
Let us now discuss the various neutrino-matter interactions that are commonly considered in neutron star merger simulations. Our objective here is not to provide detailed derivations of all interaction rates, but rather to review the reactions that may be taken into consideration and to get reasonable estimates of the scaling of reaction rates with the fluid properties. This will allow us to estimate when different reactions become important to the evolution of the system. Accordingly, for the sake of brevity, the cross-sections and reaction rates presented here sometimes make stronger approximations than what is done in merger simulations. However, for each reaction we provide references to more detailed discussions of these cross-sections. We will also make use of our discussion of the p + e − ↔ n + ν e and e + e − ↔ νν reactions to illustrate a number of issues that arise when attempting to include collision terms in the radiation transport equations, and thus discuss these reactions in more detail than the others. Given the significant overlap between reactions important to neutron star merger simulations and reactions important to core-collapse supernova simulations, a number of expressions in this section are slight modifications of the interaction rates presented in the review of neutrino reactions in core-collapse supernovae of Burrows et al (2006), though for numerical estimates of interaction rates we focus on the conditions most commonly found in neutron star mergers and post-merger remnants.

Charged-current reactions
The reactions with the strongest impact on the observable properties of neutron star mergers involve absorption and emission of ν e andν e . Indeed, these reactions are often (but not always) the main source of cooling in the system, and they are the only reactions that lead to changes in the electron fraction Y e of the fluid. In the hot, dense remnant of a BNS or NSBH merger, this mostly occurs through the reactions p + e − ↔ n + ν e ; n + e + ↔ p +ν e which are typically included at least approximately in all merger simulations that attempt to account for neutrino-matter interactions. Self-consistently calculating the forward and backward reaction rates can be difficult. Final state blocking means that these reactions depend on the distribution functions of p, n, e + , e − , ν e ,ν e . While we can typically assume equilibrium distributions at the fluid temperature and composition for n, p, e + , e − in neutron star mergers, at least in regions where neutrino-matter interactions are important, the neutrinos may be far out of equilibrium -and many approximate schemes used in simulations today do not contain enough information about the neutrino distribution function to fully account for the value of f ν in all reactions.
To illustrate these issues, and some of the ways in which they are handled in existing simulations, let us consider the cross-section per baryon for the reaction n + ν e → p + e − , the dominant absorption process in merger outflows, derived by Bruenn (1985). Following the notation of Burrows et al (2006), we get ν νe the fluid frame neutrino energy, ∆ np = (m n − m p )c 2 = 1.293 MeV the difference in rest mass energy between neutrons and protons, m e the mass of an electron, and W M a small correction for weak magnetism and recoil (2.5% for 20 MeV neutrinos) (Vogel, 1984). Neutrinos in BNS and NSBH mergers have typical energies ν 10 MeV, significantly larger than the rest mass energy of an electron. Thus, to a reasonably good approximation (for the purpose of our qualitative discussion here at least), This dependence of neutrino cross-sections on the square of the neutrino energies is found in many reactions relevant to neutron star mergers, and is going to be a significant source of uncertainty in our simulations, as many approximate transport algorithms do not provide detailed information about the neutrino spectrum.
The opacity for the absorption of ν e on n is then with f n , f p the distribution functions of neutrons and protons, and E ≈ p 2 /2m the kinetic energy of the baryons (ignoring the difference in mass between protons and neutrons and momentum transfer onto the proton). In the last expression, which ignores the final state blocking factor of the protons, n n is the neutron number density. That expression would be very inaccurate in the densest region of a star (where f p cannot be neglected), but is quite accurate in the lower-density regions where neutrinos decouple from the fluid.
To gain a more intuitive understanding of the rate of these interactions, let us assume that the typical length scale within a neutron star is ∼ 1 km. We can see from this expression that for a 20 MeV neutrino, we expect κ a = 1 km −1 for n n ∼ 10 −3 fm −3 , i.e., for a neutron mass density of ∼ 10 12 g/cm 3 . As the center of a neutron star has density ρ c ∼ 10 15 g/cm 3 , we see that neutrinos inside the neutron star have a mean free path much shorter than the size of the star, and decouple from the matter as they move through the crust of the neutron star.
Similar scalings apply to the p +ν e → n + e + reaction, as and κ a ≈ n p σν e p for the absorption ofν e on protons, under the same assumptions as for absorption onto neutrons. The correction WM is more significant than W M (∼ 15% at 20 MeV) (Vogel, 1984;Horowitz, 2002), though still not large enough to impact our order of magnitude estimates. As n p < n n in most regions of a neutron star merger remnant, the absorption opacity forν e is smaller than for ν e . It is also possible to include in simulations the impact of ν e and/orν e absorption on atomic nuclei. This is typically more important in the corecollapse context than in mergers, as in mergers most of the matter is in the form of free nucleons in regions where neutrino-matter interactions are significant. Additionally, simulations do not keep track of the abundances of individual nuclei, and equations of state for the fluid do not always contain that information, complicating any estimate of the absorption cross-section for this process. Cross-sections for the absorption of ν e onto nuclei can be found in Bruenn (1985). In high-density, low-temperature, neutron-rich regions inside of merging neutron stars, the modified URCA processes (Yakovlev et al, 2001;Alford et al, 2021) N + n → N + p + e − +ν e ; N + p + e − → N + n + ν (with N a spectator nucleon) may also play a role in the evolution of the system through the creation of an effective bulk viscosity in the post-merger remnant (Alford et al, 2018).
In the expressions derived so far for neutrino absorption, we have generally ignored final state blocking factors. These can however be approximately calculated if we rely on the fact that the fluid particles are in statistical equilibrium at a given temperature. Final state blocking factors for neutrino emission are slightly more complex to take into account. For neutrinos of a given energy and momentum, the neutrino emission rate will generally be of the form η = η * (1 − f ν ), where the (1 − f ν ) term captures Pauli blocking for neutrinos in the final state. This is not a form that is practical to use in simulations, as we would like the emission rate and opacities to depend solely on the properties of the fluid, without any dependency on f ν . Burrows et al (2006) show that a convenient redefinition of the emissivity and absorption opacity can solve this problem. If we directly use η * as our emission rate (without neutrino blocking factor), and define κ * a = κ a /(1−f eq ν ) as our absorption opacity (with f eq ν taken from Eq. (7)), then the collision term for charged-current reactions in Boltzmann's equation can be written in the two equivalent ways with η the emissivity per unit of solid angle and neutrino energy. Importantly, in the first expression η ν depends of f ν , but in the second η * ν does not. Accordingly, most simulations use η * and κ * a to parametrize neutrino-matter interactions. In our discussion of numerical algorithms for neutrino transport, emissivity and absorption opacity will generally refer to these corrected values.
We also note that when all reactions are accounted for, η * = cκ * a f eq ν (Kirchoff's law). This allows us to calculate only one of (η * , κ * a ), then set the other to make sure that the equilibrium energy density of neutrinos has the desired physical value. This is particularly useful in dense, hot regions, where neutrinos quickly reach equilibrium with the fluid. In that regime, the exact emission and absorption rate can be more difficult to calculate (due to blocking factors), but they are also fairly unimportant: what matters is that neutrinos quickly reach their equilibrium density, and then diffuse through the dense regions. This is guaranteed when using Kirchoff's law, even if η * and κ * a are not extremely accurate.
The total emission rate of neutrinos due to a given reaction can be calculated by integrating η * over both solid angle and neutrino energy. In terms of the absorption opacity, we get For comparison with results for other reactions, we can estimate this emission rate for ν e , ignoring the final state blocking factor of protons in the inverse reaction and using W M ∼ 1. We then get for the emission of electron neutrinos due to electron capture on protons (energy per unit volume) with η = µ/(k B T ). Similarly, the number of neutrinos emitted per unit volume is simply and the average energy of emitted neutrinos For η ν 1, ν ∼ 5.1k B T . We note that this is higher than the average energy of neutrinos in equilibrium with the fluid. This will generally be true whenever neutrinos are allowed to directly escape from an emission region instead of thermalizing with the fluid first. A more explicit expression for Q pe − is We see that the emission rate of neutrinos has a strong dependence in the fluid temperature, with Q ∝ T 6 , and a linear dependence in the fluid density (ignoring the Fermi integral term). The emission rate ofν e can be computed in the exact same manner, In this expression, we made use of the fact that η eq νe = −η eq νe . The dependence of these emission rates on n n and n p may seem counterintuitive, as Q pe − involves absorption of electrons on protons, yet is proportional to n n . This is however a natural result of using Kirchoff's law; the complete dependence of Q pe − in the density of all fluid particles is practically hidden in the Fermi integral term F 5 (η eq νe ), and the assumption of statistical equilibrium in the fluid. In particular, as F n (η) monotonically increase with η, and neutrino emission in post-merger remnants comes from regions of the fluid where η νe < 0 (more neutron-rich than in equilibrium), we generally get Q ne + > Q pe − even though n p < n n .

Pair processes
After charged current reactions, the most commonly considered processes for the emission and asborption of neutrinos are the pair processes i.e., electron-positron annihilation, plasmon decay, and nucleon-nucleon Bremsstrahlung. Here, each pair can be ν eνe , ν µνµ , or ν τντ . Pair processes will be the dominant source of neutrino emission for muon and tau neutrinos, as charged-current reactions involving muons and taus are significantly less common than charged current reactions involving electrons in the merger context (the mass of a muon is 105 MeV, while most of the post-merger remnant has temperature T 50 MeV, and the neutrinospheres and optically thin regions are at even lower temperatures). Pair processes are however harder to accurately include in simulations due to their nonlinear dependencies in the neutrino distribution functions. The reaction rates for the νν pair productions (forward reactions) depend on the distribution function of both neutrinos and antineutrinos through blocking factors, which are typically difficult to estimate accurately with existing transport algorithms. Worse, the reaction rates for pair annihilations (inverse reactions) are directly propoprtional to the product of the distribution functions of neutrinos and antineutrinos. Let us consider for example the reactions νν → e + e − , for neutrinos of energy significantly higher than m e c 2 (as appropriate in neutron star mergers). We can slightly adapt the results of Salmonson and Wilson (1999), based on the Newtonian rate calculations of Cooperstein et al (1986); Goodman et al (1987), to find the rate of momentum deposition per unit volume (30) with G F = 5.29 × 10 −44 cm 2 MeV −2 and D = 2.34 for electron type neutrinos, while D = 0.50 for muon or tau neutrinos. We have here chosen to rewrite the results of Salmonson and Wilson (1999) into a manifestly covariant expression more appropriate for general relativistic simulations. From this expression, we can see that the probability that a given neutrino is annihilated will depend on both the momentum of that neutrino and the distribution function of its antiparticle.
To limit the computational cost of this calculation, it is often convenient to make some assumptions regarding the distribution function of neutrinos, e.g. ignoring neutrino blocking factors (for the forward reactions), assuming equilibrium distributions of neutrinos (for either direction), or, in moment schemes, using approximate moments of the distribution functions (for the backward reactions). The most common strategy in existing merger simulations has been to compute the forward reaction rates assuming equilibrium distributions of neutrinos or ignoring blocking factors, either for all neutrinos or only for the muon and tau neutrinos. The inverse reaction rates are then computed using Kirchhoff's law, even though that law is not necessarily valid for pair processes (O'Connor, 2015). These approximations are generally reasonable for heavy lepton neutrinos close to the neutrinosphere, i.e., where most of the neutrinos that leave the remnant are emitted, because as long as charged-current reactions including muons and taus are negligible, the distribution functions of ν µ ,ν µ , ν τ ,ν τ are all identical and close to equilibrium. They are however very unreliable for electron-type neutrinos and for calculations of the rate of νν annihilation in regions where the neutrinos are not in equilibrium with the fluid (e.g. in the polar regions). We will consider some of the ways in which the latter process has been studied in our discussion of specific radiation transport algorithm and simulations. Now that we have established the difficulty of properly treating pair processes, let us estimate their importance in the merger context. We begin again with e + e − annihilation. Ignoring neutrino blocking factors, Burrows et al (2006) integrate the reactions rate of Dicus (1972) to find a total emissivity in νν pair of with Q 0 = 9.76×10 24 ergs cm −3 s −1 for ν eνe , and Q 0 = 4.17×10 24 ergs cm −3 s −1 for all other neutrinos combined. The average energy of the emitted neutrinos in the fluid frame is (e.g., ν ≈ 4.1T when η e = 0). Approximate expressions for the energy spectrum of the neutrinos are also found in Burrows et al (2006), following the work of Bruenn (1985). If we compare this result to Q pe − and Q ne + from the previous section, we see that for ν eνe , pair processes will only dominate over charge current reactions in very hot and/or low density regions of the fluid, where neutrinos will either rapidly reach their equilibrium distribution or rapidly cool the fluid. In such cases, getting exact reaction rates is not overly important as long as we obtain the correct equilibrium distribution and have sufficiently high emission rates. Even in dense regions of the fluid, neglecting ν eνe production is thus not a particularly strong approximation (but neglecting pair annihilation in low-density regions might be, as we will see).
What about the heavy-lepton neutrinos? The equilibrium energy density of neutrinos is At T = 1 MeV, the timescale for neutrinos to reach that equilibrium density solely through e + e − emission is thus O(10 s), but at T = 10 MeV, it is O(0.1 ms), i.e., much shorter than the dynamical timescale of a neutron star merger. In hot regions, heavy-lepton neutrinos (muons and taus) will thus reach their expected equilibrium density, and the neutrino luminosity of ν µνµ ν τντ will be set by the diffusion timescale of neutrinos through the hot, dense remnant. For heavy-lepton neutrinos, ignoring pair processes (and missing the associated cooling of the remnant) would be significantly worse than incuding approximate reaction rates, as long as those rates properly recover the equilibrium energy of neutrinos in dense regions.
Let us now briefly consider other pair processes. For nucleon-nucleon Bremsstrahlung, Burrows et al (2006) (building on results by Brinkmann and Turner (1988); Hannestad and Raffelt (1998)) find the total neutrino emissivity per species to be Q nb ≈ (1.5 × 10 26 ergs cm −3 s −1 ) n n 10 36 cm −3 We see that Bremsstrahlung will dominate over e + e − annihilation in denser, colder regions. In the densest region of a post-merger accretion disk (typically n n ∼ 10 35−37 cm −3 , T ∼ (1 − 10) MeV), we see that the process dominating the production of heavy-lepton neutrinos may thus vary, and we can neglect neither Bremsstrahlung nor e + e − pair production/annihilation. Approximate formula for the total energy emission from plasmon decays and for the average energy of the neutrinos emitted through that process can be found in Ruffert et al (1996). They are equivalent to with Q 0,pl = (6 × 10 23 ergs cm −3 s −1 ) per species for electron-type neutrinos and Q 0,pl = (10 21 ergs cm −3 s −1 ) per species for other neutrinos. The blocking factor B = (1 − f ν ) (1 − fν) can only be evaluated assuming a specific neutrino distribution function, while γ ≈ 0.056 (π 2 + 3η 2 e )/3 is a parameter with strong sensitivity to the electron degeneracy parameter η e . We see that we need relatively fine-tuned conditions for plasmon decay to dominate over pair annihilation, especially for heavy-lepton neutrinos -as a result, this reaction is often ignored in merger simulations.
From these estimates of the emissivity of pair processes, we can understand one additional difficulty in the use of these processes in simulations. Both e + e − creation/annihilation and plasmon decays have Q ∝ T 9 , with no explicit dependence in the fluid density (at least in regions where blocking factors are negligible). This can prove problematic in merger simulations, where numerical errors can lead to the creation of hot low-density regions whose properties are not necessarily well modeled by equations of state built to capture the properties of dense matter. As a result, some simulations ignore pair processes below an ad-hoc density threshold.

Neutrino scattering
Scattering of neutrinos on protons, neutrons, nuclei and electrons plays an important role in setting the diffusion timescale of neutrinos through the densest regions of merger remnants. The total cross-sections per baryon for the nearly elastic scattering of neutrinos onto protons and neutrons are (Yueh and Buchler, 1976;Bruenn, 1985;Burrows et al, 2006) and the scattering opacity for these two processes combined is thus We see that these quasi-elastic scatterings are about as likely as chargedcurrent absorption for ν e andν e , and will often be a dominant contribution to the total opacity of the fluid for heavy-lepton neutrinos. We also note that the differential cross-sections are with µ = cos θ and θ the scattering angle. As δ p ∼ −0.2 and δ n ∼ −0.1 at the most relevant neutrino energies (Burrows et al, 2006), back-scattering is favored. However, most merger simulations assume isotropic elastic scatterings in the fluid frame (δ n,p = 0), an approximation whose impact on simulation results has not been tested so far. Similar calculations can be made for scattering on atomic nuclei. For example, elastric scattering on α particles has a total cross-section per nucleus of (Yueh and Buchler, 1976;Burrows et al, 2006) We note that in the neutron star merger context, we typically have n α n n in regions where neutrino scattering is important, and similar results apply to heavier nuclei. As many equations of state used in merger simulations do not provide detailed information about the abundances of individual atomic nuclei, the contribution of nuclei to the total scattering opacity is often only approximately taken into account (e.g. considering only α particles, or α particles and some 'representative' nucleus of fixed proton number Z and atomic number A), or completely ignored.
Including inelastric scattering of neutrinos on electrons is a more difficult problem, and as a result inelastic scattering has not so far been taken into account in merger simulations. To understand these issues, we can look at the methods used to treat inelastic scattering in core-collapse supernovae (Bruenn, 1985;Burrows et al, 2006). The relevant part of Boltzmann's equation can be written (40) Note that f ν is the distribution function for neutrinos with energy ν and momentum p µ , while f ν is the distribution function for neutrinos with energy ν and momentum p µ . R in,out are the scattering kernels to scatter into/out of the energy bin ν from/to ν . Even ignoring the blocking factors (1 − f ν ) and (1 − f ν ), the collision terms clearly depend on the distribution function of neutrinos, and couple the values of f ν at all neutrino momenta. One possible approximation is to use a truncated expansion of the kernels in cos θ: with Φ 0,1 known functions of the incoming and outgoing neutrino energies. The integrals over p i are then similarly truncated using moments of the distribution function f ν . While this makes the evolution of f ν slightly more tractable numerically, we still end up with numerically stiff terms coupling every pair of neutrino energies, which makes these reactions expensive to include in simulations. The scattering kernels have complex dependencies in the incoming and outgoing neutrino energies, and the temperature of the fluid (which sets the electron distribution function). As a very rough order of magnitude estimate, and assuming that the neutrinos have energies larger than or comparable to the electrons, we have with ν the typical energy of neutrinos at the current point, and e the typical energy of electrons. This leads to an effective opacity (i.e., the inverse of the mean free path of neutrinos with respect to scattering on electrons) We thus see that at the densities at which neutrino-matter interactions are most important in neutron star mergers, inelastic scattering on electrons has a significantly lower opacity than elastic scattering on nucleons or chargedcurrent reactions, but not necessarily smaller than absorption opacities for pair processes. Accordingly, its direct impact on ν e andν e is likely subdominant, but it could be important to the thermalization of heavy-lepton neutrinos. Finally, we note that scattering on nucleons is not perfectly elastic. The typical exchange of energy between neutrinos and the fluid is much smaller for nucleon scatterings than for electron scatterings, but as seen above, the crosssections for nucleon scatterings are larger in regions where neutrino-matter interactions are important. In core-collapse supernovae, Wang and Burrows (2020) showed that the smaller energy transfer during each scattering can be used to treat inelastic scattering on nucleons as a diffusion process in energy space, leading to much cheaper calculations than when using scattering kernels: there is no need to couple all energy bins through numerically stiff interaction terms. In Wang and Burrows (2020), the impact of neutrino-nucleon scattering on the thermalization of heavy-lepton nucleons was also shown to be comparable to the impact of neutrino-electron scattering. Accounting for inelastic scattering on nucleons could thus provide an avenue to partially account for the thermalization effect of scattering events without an implementation of inelastic scattering on electrons.

Discussion
From the previous sections, we see that the reactions currently used in our most advanced merger simulations can, if properly included in a transport algorithm, capture the dominant processes for emission, absorption, diffusion, and thermalization of ν e andν e in most of a post-merger remnant. Without even getting into the complications of approximate transport methods, however, we see that the situation is already more complex for other species of neutrinos. Emission of heavy-lepton neutrinos is dominated by pair processes which are poorly modeled as soon as neutrinos are out of equilibrium with the fluid. Thermalization of these neutrinos is likely impacted by inelastic scattering, which current simulations do not take into account. Finally, pair annihilation of all types of neutrinos in low-density regions is difficult to include, but possibly important to jet formation. It is thus worth noting that uncertainties in transport schemes are not the only potential sources of errors in our modeling of neutrinos today; the choice of physical processes included in the simulations, and the accuracy to which they are modeled, remains an area where significant improvements are possible.

Quantum kinetics and neutrino oscillations
So far, we have considered neutrinos as particles in well-defined flavor states (electron, muon, tau). However, we know that this is only an approximation. Even in vacuum, the fact that the mass eigenstates of neutrinos are different from their flavor eigenstates leads to oscillations between flavors. Vacuum oscillations occur on length scales too long to impact the evolution of a postmerger remnant, though if neutrinos from a neutron star merger were ever to be observed, oscillations between the source and the Earth would certainly be significant. There are however other processes that lead to flavor transformation with more relevance to the merger problem. Generally, any process that transform electron type neutrinos into heavy-lepton neutrinos (or viceversa) close enough to the merger remnant that neutrino-matter interactions are still impacting the composition of the outflows has the potential to change the properties of kilonovae and the outcome of nucleosynthesis in neutron star mergers.
One way to study neutrino oscillations is through the quantum kinetic equations (QKE). In that formalism, neutrinos are described by the 3 × 3 density matrix ρ(t, x i , p µ ). The diagonal terms of this matrix can be understood as equivalent to the distribution functions f νe , f νµ , f ντ , while the off-diagonal terms encode quantum coherence between flavors. A second matrixρ contains information about antineutrinos. The density matrix evolves according to (Vlasenko et al, 2014) where the left-hand side is a total time derivative in phase-space, and the two terms on the right-hand side are responsible for, respectively, oscillations and collisions. The Hamiltonian H can be decomposed as with H vac responsible for vacuum oscillations, H mat for interactions between neutrinos and the matter potential, and H νν for neutrino self-interactions. At least two types of oscillations have been found to be potentially improtant in the merger context. The Matter-Neutrino Resonance (MNR) occurs when the matter potential is equal to the neutrino self-interaction potential, and can impact the luminosity of ν e andν e within a few radii of the post-merger remnant (Caballero et al, 2014;Zhu et al, 2016). The flast-flavor instability (FFI), on the other hand, is due solely to the neutrino self-potential, and seems to occur in regions where the sign of the net lepton flux (number flux of ν e minus number flux ofν e ) changes between different directions of propagation of the neutrinos (Banerjee et al, 2011;Wu et al, 2017;Grohs et al, 2022). The FFI occurs on very short timescales (∼ns, i.e., cm length scales), and is likely active in many regions close to the post-merger remnant (Grohs et al, 2022). How much flavor transformation occurs as a result of the FFI remains uncertain, but recent studies using simplified prescriptions for where the FFI occurs and how much flavor transformation happens as a result have shown that it could plausibly lead to significant changes in the composition of matter outflows (Li and Siegel, 2021;Fernández et al, 2022). As quantum kinetics is not at this point studied as part of general relativistic radiation transport algorithms coupled to merger simulations, but rather evaluated using either simple approximations or specialized zoomed-in simulations, we do not discuss it in more detail here. We do however emphasize that these oscillations could very well have an impact on the composition of the matter outflows produced in mergers. The fact that they are not included directly within simulations is due largely to the additional technical difficulty of evolving the quantum kinetic equations and to the very short timescales involved in the FFI, rather than to a certainty that oscillations are not important to astrophysical results. Obtaining better models for the role of oscillations in merger simulations is certainly an important open problem in merger simulations today.

Radiation transport algorithms
Having discussed the reactions that we would like to take into account in neutron star merger simulations, we can now turn to a discussion of the various methods used so far to treat neutrino transport and neutrinomatter interactions. These can be broadly classified into quasi-local leakage schemes, approximate transport schemes based on the moment formalism, and Monte-Carlo evolution of Boltzmann's equation. Multiple simulations have also considered mixed leakage-moment schemes, while algorithms mixing Monte-Carlo methods with a moment scheme have been considered but not successfully used in merger simulations. For most of this section, we attempt to keep the discussion focused to the methods used for general relativistic radiation hydrodynamics simulations, either in the context of neutron star merger simulations or for the evolution of their post-merger remnant. We will however discuss along the way a number of techniques that were first developed for Newtonian simulations or for simulations using quasi-Newtonian potentials that have either been ported to general relativistic simulations, or are likely to be used in that context in the near future. This is particularly true for the more advanced leakage schemes, which were first developed in non-relativistic codes but are currently being integrated in general relativistic simulations. We also note that neutrino radiation transport algorithms used in simulations of core-collapse supernovae are often more advanced than any of the algorithms used in merger simulations (e.g. Mezzacappa and Bruenn (1993) (2019)). In fact, many algorithms used in merger simulations today are directly inspired from work done in the core-collapse community. Accordingly, while we do not attempt to review the algorithms used in the core-collapse context, we will occcasionally refer to methods developed for core-collapse simulations if they have been used in the merger context. More advanced methods have also been proposed, but not yet applied to the merger problem; e.g. methods for a fully covariant evolution of the radiative transport equations appropriate for a direct discretization of Boltzmann's equation have recently been studied in Davis and Gammie (2020), Lattice-Boltzmann methods have been implemented and used on test problems in Weih et al (2020a), and the MOCMC (Method of Characteristics Moment Closure) method has shown that it is possible to combine particle and moment formalisms to improve on the convergence properties of a pure Monte-Carlo radiation transport code (Ryan and Dolence, 2020).
In this review, we focus particularly on moment methods, as they have been used in the majority of the most advanced radiation hydrodynamics simulations of mergers to-date. Most general relativistic simulations using leakage schemes use methods that are at best order-of-magnitude accurate, while very few simulations have been performed with the recently developed Monte-Carlo algorithms. Accordingly, moment schemes remain at the moment our best source of information about the role of neutrinos in neutron star mergers.
We note that while most simulations consider 3 species of neutrinos and 3 species of antineutrinos, it is fairly common for simulations to assume that the distribution function of all heavy-lepton neutrinos ν µ , ν τ ,ν µ ,ν τ are identical, and thus to replace the evolution of those 4 species by the evolution of a single species ν x that represent them all; we will use the notation ν x to represent all heavy-lepton neutrinos here as well.

Overview
Leakage algorithms are the simplest methods used to treat neutrinos in neutron star merger and post-merger simulations. In their most basic form, they can capture the cooling of the post-merger remnant at the order-of-magnitude level, but not the evolution of the composition of the outflows. More advanced leakage schemes have however been developed for post-merger simulations, and mixed leakage-moment schemes have been used in general relativistic merger simulations. Those advanced schemes can at least approximately capture absorption within the outflows of neutrinos emitted by a post-merger accretion disk or by a neutron star remnant.
The leakage schemes used in merger simulations today were first developed by Ruffert et al (1996); Rosswog and Liebendörfer (2003). They generally rely on a local computation of the neutrino energy and number emission rates per unit volume, Q ν,free and R ν,free , for each species of neutrinos. In addition, leakage schemes compute an estimate of the optical depth between a given grid cell and the outer boundary of the computational domain, in order to estimate the diffusion timescale of trapped neutrinos through the remnant.
The energy emission rate Q ν,free is calculated as described in Sect. 3.2. In the merger context, simulations have usually considered the total (energyintegrated) emission rate, but energy-dependent leakage schemes have been developed for post-merger simulations . In the former case, R ν,free = Q ν,free / ν , with ν the average energy of emitted neutrinos. In the latter case, R ν,free for each energy bin is just Q ν,free / ν , with ν the energy at the center of the bin. In optically thin regions, this is sufficient to calculate the cooling rate and composition changes of the fluid. In optically thick regions, however, the rates at which neutrinos carry away energy and lepton number are much lower than the free emission rates. In those regions, we expect neutrinos to quickly reach their equilibrium distribution function f eq ν , and to slowly diffuse out of the remnant over a time scale t diff . The rate at which neutrinos carry energy away from a given cell is then approximately given, for neutrinos of a given energy ν, by with E eq ν the equilibrium energy density of neutrinos. Most leakage schemes used in merger simulations implement the diffusion time scale prescription of Rosswog and Liebendörfer (2003) where κ tot is the total opacity at the current point (including all absorption and scattering processes considered in the simulation), and τ ν is the estimated optical depth between that point and the domain boundary. The parameter α diff is calibrated to the result of transport simulations; Rosswog and Liebendörfer (2003) use α diff = 3, but this choice is not unique (e.g. O'Connor and Ott (2010) argue for an increase of α diff by a factor of two). The optical depth is defined as with the minimum taken over all possible paths Γ starting from the current point and ending at the boundary of the computational domain. The optical depth τ ν is energy dependent, but the reactions that dominate the calculation of κ tot all have κ ∝ ν 2 (charged-current reactions and elastic scatterings).
Calculating τ ν at a single energy and then assuming τ ν ∝ ν 2 is thus a reasonable approximation. We discuss different methods to estimate τ ν later in this section. More complex estimates of t diff are also possible; for example, in the Improved leakage-equilibration-absorption scheme (ILEAS) of Ardevol-Pulpillo et al (2019), separate diffusion timescales are calculated for the number and energy emission rate from the local gradient of the number density and energy density. Such a local calculation also has the advantage of allowing for emission rates that match the expected diffusion limit in optically thick regions, which is not possible using the simpler dimensional analysis of earlier schemes.
In an energy integrated leakage scheme, one then needs to integrate Q diff (ν) and R diff (ν) separately over ν: or, if we assume t diff ∝ ν 2 , and similarly for the number diffusion rate The average energy of escaping neutrinos in the diffusion regime is then We note that the average energy of diffusing neutrinos is significantly lower than the average energy of a thermal spectrum, reflecting the fact that low energy neutrinos diffuse faster than high energy neutrinos.
The actual rate at which neutrinos leave a given region of the fluid is then given by an interpolation between the estimates valid at low and high optical depth, effectively considering that neutrino transport is limited by the lowest of those two rates. Ruffert et al (1996) uses Alternatively, Sekiguchi (2010) considers an exponential transition between the two regimes These results can then be coupled to the evolution of the fluid equations using conservation of energy-momentum and lepton number The neutrino pressure term can for example be computed assuming a relativistic gas of neutrinos in equilibrium with the fluid, i.e., P νi = E eq νi /3 (O'Connor and Ott, 2010), potentially suppressed in regions where neutrinos are not trapped. Alternatively, we will see below that more advanced leakage schemes have been coupled to evolution equations for the trapped neutrinos -in which case the neutrino pressure is calculated directly from the estimated energy density of trapped neutrinos.

Leakage in general relativistic merger simulations
In general relativistic merger simulations, the first published leakage scheme was developed by Sekiguchi (2010). This algorithm is a mixed moment-leakage scheme with significantly more complexity that the simple scheme described above. The algorithm assumes evolution equations for the stress-energy tensor T αβ fl+tr of the fluid and trapped neutrinos combined, T αβ st of streaming neutrinos, and for the lepton fraction Y l . The streaming neutrinos are evolved using a moment scheme (see Sect. 4.2), which allowed later iterations of this algorithm to realitively easily take into account reabsorption of the streaming neutrinos in low density regions. We can see from the evolution equations that this algorithm has the advantage of guaranteeing exact conservation of energy-momentum. Additionally, the algorithm evolves the fractions Y e , Y νe , Yν e and Y νx of electrons and neutrinos, assuming that these fractions reach their equilibrium value (at given temperature, density and lepton fraction) in regions where neutrino-matter interactions are fast compared to the numerical time step. This allows for relatively simple estimate of the contribution of neutrinos to the fluid pressure, assuming that the neutrinos are a relativistic gas.
An algorithm closer to the original methods of Ruffert et al (1996); Rosswog and Liebendörfer (2003) was first used in general relativistic simulations by Deaton et al (2013). In that work, the only contribution of neutrinos to the evolution of the system is the source terms of Eqs. (55)-(56). In Deaton et al (2013), the minimum optical depth was calculated by considering lines along the coordinate directionsx,ŷ,ẑ of a cartesian grid, as well as along the diagonals of a cube in the same coordinates, a method similar to that previously used by Ruffert and Janka (1999). This algorithm however requires global communications between all points of the numerical grid whenever τ ν is computed, and creates preferred directions along the axis of the cartesian coordinates. An improved method to calculate τ ν was later proposed by Neilsen et al (2014), and is now the most commonly adopted algorithm in numerical relativity simulations. Their method relies on finding the path of shortest optical depth linking neighboring cell centers on a grid. We can indeed discretize our equation for τ ν at pointx as with the minimum taken over all neighboring pointsx n . Here, ∆s n is the distance betweenx andx n . Given an initial guess τ 0 ν for the optical depth at each point, we can solve this equation iteratively, using < at all points for some small constant . Because the optical depth evolves slowly over time, a single iteration initialized with the value of τ ν at the previous time step is generally sufficient to maintain a good estimate of τ ν everywhere, except when computing τ ν for the first time. This method is now widely used in neutron star merger simulations (Foucart et al, 2014;Radice et al, 2018;Mösta et al, 2020;Cipolletta et al, 2021;Most and Raithel, 2021). A conceptually similar algorithm that does not rely on the existence of an underlying cartesian grid has also been developed in Perego et al (2014a), allowing for the easy use of this method in grid-less simulations (e.g. SPH), while an improved numerical methods to solve for τ ν by solving the eikonal equation has been proposed by Palenzuela et al (2022).

Leakage limitations and improved leakage schemes
The accuracy of a leakage scheme can be very problem dependent. The free parameters in most leakage schemes are calibrated to spherically symmetric transport problems, and tend to perform best in that context -while neutron star mergers and their post-merger remnants are very asymmetric. Even ignoring symmetry issues, however, the standard scheme discussed above has a number of important limitations. We have already mentioned the fact that the simplest leakage schemes do not accurately capture the local diffusion rate of neutrinos in the high optical depth limit; in this section we consider a few additional notable issues.
First, in regions where the total optical depth is high but the absorption and inelastic scattering optical depths are low, the assumption that neutrinos reach their equilibrium distribution function can be inaccurate. This is particularly problematic for heavy-lepton neutrinos, which typically have much lower absorption opacity than scattering opacity. One way to solve this issue is to keep track of the energy density of neutrinos, rather than assuming an equilibrium energy density. In Newtonian simulations of post-merger disks, the Advanced Spectal Leakage (ASL) scheme of Perego et al (2016), for example, calculates the energy and number density of trapped neutrinos assuming that the distribution function of trapped neutrinos f tr ν satisfies the equation for the estimated production time scale t ν,prod and diffusion time scale t ν,diff . Here, ∆t is the time step of the evolution. We see that in this scheme, if the time scale to produce neutrinos is long, the trapped neutrinos do not reach equilibrium with the fluid. Perego et al (2016) further assume with τ en the optical depth ignoring elasctic scatterings, and γ a function of position only. This allows for the use of the single unknown γ to represent the function f tr ν . The ASL scheme has been adapted for use in merger simulations within a smoothed particle hydrodynamics (SPH) code (Gizzi et al, 2021), which in conjunction with the development of a first general relativistic SPH code (Rosswog et al, 2022) should allow for the use of ASL in general relativsitic merger simulations in the near future. We note that the general relativistic algorithm of Sekiguchi (2010) also keeps track of the energy density of trapped neutrinos when, in the notation of Perego et al (2016), t diff ≤ t prod . However, as Sekiguchi (2010) calculates the neutrino diffusion rate from the equilibrium density of neutrinos rather than from the energy density evolved within the simulation, the scheme only partially correct for the difference between those two estimates of the energy. The ILEAS scheme (Ardevol-Pulpillo et al, 2019) similarly tracks trapped neutrinos through an "equilibration-advection" step where neutrinos in optically thick regions are assumed to be in equilibrium with the fluid, and "trapped" neutrinos are otherwise advected through the computational grid to guarantee exact conservation of lepton numbers as well as to approximately account for the impact of trapped neutrinos on the properties of the fluid.
A second common issue with leakage schemes is that the simple interpolation between free emission and diffusion is naturally going to be inexact in semi-transparent regions. Corrections to the emission rate can be calibrated to specific systems (e.g., Perego et al 2016; Gizzi et al 2021 for accretion disks and mergers), but the region where most of the neutrinos leaving the system are emitted is, by definition, particularly difficult to model accurately in a leakage scheme.
A third issue is that the energy of neutrinos diffusing through the system is not constant; inelastic scatterings, as well as absorption and re-emission, are both taken into account in the energy diffusion rate, and these processes tend to bring the neutrinos closer to thermal equilibrium with the fluid. On the other hand, the integrated values of Q diff and R diff computed above assume that neutrinos random walk through the fluid at constant energy. Diffusing neutrinos have lower average energies than the emitted neutrinos, but only because the opacity of the system to high-energy neutrinos is larger; the scheme does not account for thermalization of the neutrinos. The ASL scheme attempts to partially account for this thermalization by suppresing R ν by a factor of e −τen/10 , then renormalizing the energy-intergrated R ν to keep the total number of neutrinos constant while reducing their average energy.
The last important issue is that, in its simplest form, a leakage scheme does not account for the energy, momentum, and lepton number deposited in the fluid by reabsorption of neutrinos emitted in different regions of the system. Neutrinos either leave or do not leave the simulation domain, but the interactions between neutrinos and the fluid after emission do not feedback onto the evolution of the fluid. This is particularly problematic when studying the composition of matter outflows, which is known to be significantly impacted by absorption of ν e (Wanajo et al, 2014). In general relativistic merger simulations, multiple leakage schemes have attempted to approximately include these absorption effects. In Mösta et al (2020); Cipolletta et al (2021), absorption is taken into account by calculating the emission and absorption of neutrinos along radial rays in the post merger remnant. Along each ray, the neutrino emissivity is estimated from the energy and number density emission rates predicted by a leakage scheme. The heating rate and change in composition of the fluid can then be estimated from the neutrino luminosity integrated over each ray, and from a local estimate of κ a . Radice et al (2018) instead consider approximate transport along each ray, evolving the energy density of the free-streaming neutrinos (zeroth moment, see Sect. 4.2). In all of these works, absorption of free streaming neutrinos is only included in optically thin regions, with an exponential cutoff ∝ e −τν applied to the absorption rate in optically thick regions. In the ILEAS scheme (Ardevol-Pulpillo et al, 2019), neutrino emission (as predicted by the leakage) is similarly followed along "rays" and potentially reabsorbed. However, no specific geometry is assumed for these rays; instead, neutrinos follow the gradient of the optical depth. This leads to many more rays that may jointly contribute to absorption in any given cell, thus complicating the algorithm, but also to a propagation of neutrinos that should better match the geometry of the system.
We note that whenever an attempt is made at taking into account neutrino absorption, it is important to calculate the absorption opacities for the energy of the free-streaming neutrinos, rather than for the energy of neutrinos locally in equilibrium with the fluid. Typically, neutrinos in merger simulations have ν ∼ (10 − 20) MeV, while in the outflows T ∼ 1 MeV. As κ a ∝ ν 2 , assuming thermal equilibrium of the neutrinos with the fluid can lead to large underestimates of the absorption rate (Foucart et al, 2016b;Radice et al, 2018). This is an issue not only for leakage schemes, but also for energy-integrated moment schemes, as we will see.
In Newtonian (or pseudo-Newtonian) post-merger simulations, simple lightbulb prescriptions have instead been used to model the spatial distribution of emitted neutrinos (Fernández and Metzger, 2013;Metzger and Fernández, 2014). In a lightbulb model, the total luminosity of neutrinos is calculated by integrating over the entire simulation, then that luminosity is assumed to come from a specific region -in the case of post-merger systems, an annulus around the densest region of the remnant accretion disk and/or the surface of the remnant neutron star (if present). Alternatively, simulations using the ASL scheme have combined that scheme with a more advanced approximate transport algorithm (see Perego et al 2014b;Gizzi et al 2019Gizzi et al , 2021. In that algorithm, neutrinos emitted in optically thick regions (τ ν > 2/3) are assumed to diffuse to the neutrinopshere (τ ν = 2/3 surface) along the gradient of τ ν (i.e., along the path of smallest optical depth), before free-streaming from all points of the neutrinosphere. Neutrinos emitted in optically thin regions just freestream from their point of emission. The number density of neutrinos outside of the neutrinosphere can then be combined with calculations of the absorption and scattering optical depth to determine energy and momentum deposition in optically thin regions. As the calculation of the number density of neutrinos requires first the determination of a map linking each point at τ ν > 2/3 to a point on the neutrinosphere, and then for each point at τ ν < 2/3 an integral over the entire neutrinosphere and all optically thin regions, this algorithm is significantly more expensive than other leakage schemes, even in Newtonian physics where free-streaming neutrinos can be assumed to propagate along straight lines.

Discussion
Overall, we thus see that leakage schemes provide us with a simple, costeffective method to approximately incorporate neutrino cooling in merger simulations. Using a basic leakage scheme typically has no noticeable effect on the cost of a merger simulations, but the results are only order-of-magnitude accurate and do not account for the important role of neutrino absorption in matter outflows. More advanced leakage methods can be developed through coupling to radiation transport algorithms, or by attempting to predict where the leaking neutrinos will go as they leave the system. These more advanced algorithms will naturally be more costly, but they have been shown to provide a better match to the result of full radiation transport simulations, at least on test problems. Estimating the error of such an algorithm without comparison to a more advanced simulation is however always difficult.

Moments-based radiation transport
To go beyond leakage schemes and attempt to transport neutrinos along geodesics through a numerical simulation, multiple general relativistic merger codes have now implemented "tuncated moments" schemes, i.e., algorithms evolving moments of the distribution function of neutrinos. The formalism for general relativistic moments schemes was derived by Thorne (1981), while the development of methods for the coupled evolution of the fluid and moment equations largely build on work performed for photon transport in Newtonian simulations (Stone et al, 1992;Audit et al, 2002;Vaytet et al, 2011;Skinner and Ostriker, 2013). For general relativistic transport, it saw early uses in spherical symmetry (Rezzolla and Miller, 1994), before being adapted to the methods used in 3D neutron star merger simulations (Shibata et al, 2011). Our understanding of the moment equations in relativistic systems also significantly benefited from work done in the context of photon transport in accretion disks (Sadowski et al, 2013). The moment formalism leads to evolution equations that are typically more complex and costly to evolve than the leakage scheme, but automatically include emission, transport, and reabsorption of neutrinos. We will see that moment schemes have been very successful in providing us with an improved understanding of the role of neutrinos in neutron star mergers, but also that the truncated expansion of the distribution function used in these schemes will force us to choose approximate analytical closure for unknown higher moments of f ν , leading to evolution equations that do not converge to the true solution of Boltzmann's equation, with hard-to-quantify errors in the results.
We organize our discussion of moment schemes as follow. In Sect. 4.2.1, we discuss the basic ideas behind the moment formalism and the derivation of the moment equations. In Sect. 4.2.2, we review evolution equations for the moments in the frame of a numerical simulation as well as prescriptions to account for energy and momentum transfer between neutrinos and the fluid. In Sect. 4.2.3, we discuss lepton number exchange and the option for energyintegrated moment schemes to evolve the neutrino number density to guarantee exact lepton number conservation. In Sect. 4.2.4, we provide more details on the calculation of the coupling terms between neutrinos and matter, with a focus on issues that arise when the energy spectrum of neutrinos is not known. In Sect. 4.2.5 we discuss the approximate analytical prescriptions used to close the moment equations, while Sect. 4.2.6 focuses on the numerical implementation of the moment formalism in general relativistic merger codes. In Sect. 4.2.7 we discuss specific issues with the use of moment schemes in high-density regions, and ways to recover the proper diffusion rate of neutrinos through those regions. Finally, Sect. 4.2.8 discusses approximate implementations of neutrino-antineutrino pair annihilation in low-density regions.
In this section, we choose units such that h = c = 1.

Truncated moments formalism
To understand the basic idea behind the moment formalism, let us write the neutrino 4-momentum as witht µ a timelike unit vector, and l µ a spacelike unit vector orthogonal tot µ . We see that is the neutrino energy for an observer with 4-velocityt µ , and that by construction p µ p µ = 0. The n th moment of the distribution function f ν according to our observer with 4-velocityt µ is defined as with dΩ an integral over solid angle on the unit sphere in momentum space. For example, the 1 st moment is while the energy density of neutrinos (0 th moment) for that same observer is We thus see that only the spatial components of M α are independent of the 0 th moment; these are usually denoted as the flux density with, by construction, F αt α = 0, i.e., F α is a purely spatial vector according to our chosen observer. Combining these results, we get Similarly, for the second moment, M αβt α = −M β and M αβt β = −M α , and so the only components of the second moment that cannot be directly reconstructed from E and F α are the spatial components of the pressure tensor with P αβt α = P αβt β = 0. With those definitions, we can write the second moment as M αβ = Et αtβ + F αtβ + F βtα + P αβ .
We note that these moments are well-defined tensors, but they do generally require the choice of a specific frame in which neutrino energies are measured. The one exception is the energy integrated second moment of f ν , which can be written as As d dΩ = d 3 p/(pt √ −g) is the invariant integration volume in momentum space, written in an orthonormal frame with timelike coordinatet α , that expression is independent of the coordinates in which we measure . In fact, comparing with Eq. (5), we see that M αβ tot is the stress-energy tensor for the chosen species of neutrinos.
The first moment equation can be derived by taking the covariant divergence of the first moment with neutrino energies measured in the fluid rest frame (i.e., takingt µ = u µ ), and combining it with Boltzmann's equation. We get (Thorne, 1981) with ν the neutrino energy in the fluid rest frame. The right-hand side is defined from with S α the first moment of S in the fluid rest frame, i.e.
In merger simulations using the moment formalism, the source term has so far been limited to isotropic emission, absorption, and isotropic elastic scattering of neutrinos. Then with η the emissivity, κ a,s the absorption and scattering opacities, and J, H α the energy density and momentum flux in the fluid frame, i.e.
for our choice oft µ = u µ . The second moment equation is (Thorne, 1981) and we typically decompose M αβ in the fluid frame as with L αβ the pressure tensor in that reference frame. For an observer using an orthonormal tetrad in the fluid rest frame, we can interpret the second moment equation as an evolution equation for J and Hî. Indeed, However, the evolution equations are not closed: they depend on the next two moments of f ν , L αβ and M αβγ . To obtain a well defined system of equations, one thus needs a closure relation. That closure should provide higher moments of f ν that are not evolved in our equations, as a function of the evolved variables. This is practically similar to the need for an equation of state P (ρ, T, Y e ) in the evolution of the fluid equations, except that in a fluid P can be calculated assuming statistical equilibrium of the fluid particles. For out-of-equilibrium particles, we do not have an analytical expression for the closure; our choice of closure will thus introduce an error in our solution.
We note that everything in our derivation so far has assumed that the moments are functions of the neutrino energy in the fluid frame. These functions can be discretized in energy space to obtain a spectral or energydependent moment scheme. Such a scheme has been used for the evolution of post-merger remnants in Newtonian simulations (Just et al, 2015) and in general relativistic simulations of core-collapse supernovae (Kuroda et al, 2016;Roberts et al, 2016), but not so far in general relativistic merger simulations. As our focus here is mainly on the latter, we will mainly discuss the cheaper alternative: evolving energy-intergrated moments of f ν using a grey moment scheme. Grey moment schemes require different closure relations. On the one hand, we can see from Eq. (79) that, under the physical assumption that f ν drops to zero faster than ν −1 as ν → ∞, the term involving M αβγ in that equation does not contribute to the energy-intergrated moment equation. On the other hand, given the strong dependence of κ a,s on the energy of the neutrinos, the caculation of an energy-integrated S α requires a choice of neutrino spectrum -a new closure relation that will significantly impact the assumed cross-sections of neutrino-matter interactions. We will come back to this choice later in this review.

Moment equations in the simulation frame
One of the main advantage of the moment formalisms is that the evolution equations for the moments can be put into a "conservative" form very similar to that commonly used for the evolution of the fluid equations (see below). To do so, however, it is useful to move away from moments computed in the fluid rest frame. Consider instead the decomposition of the second moment with n α the unit normal to a constant-time slice in the simulation, and F α n α = P αβ n α = P αβ n β = 0. We can then interpret E, F α , P αβ as the energy density, flux density, and pressure tensor of the neutrinos of a given species as measured by a normal observer. Remembering that the energy integrated second moment is independent of the reference frame in which we measure the neutrino energy, this expression makes it easy to calculate E, F α , P αβ as functions of J, H α , L αβ (and vice-versa) using projections of the second moment. Let us define the Lorentz factor W amd the fluid 3-velocity v µ such that u µ = W (n µ + v µ ) and v µ n µ = 0, and the 3-metric γ αβ = g αβ + n α n β . We then get E = M αβ n α n β = JW 2 − 2W (H α n α ) + L αβ n α n β (84) Alternatively, using h αβ = g αβ + u α u β , we get We will not need this last expression, and thus do not provide a more explicit expansion. We note that n µ = (−α, 0, 0, 0), with α the lapse function, which simplifies many expressions when calculating E, F α , P αβ from J, H α , L αβ . We also note that F α , P αβ , v α are all purely spatial tensors, e.g. F t = 0. Indices for these tensors can be raised and lowered with the 4-metric using all indices, or with the 3-metric using only spatial indices (e.g. v i = γ ij v j ). The second moment equation can now be recast as evolution equations for the energy integrated moments weighted by √ γ, with γ the determinant of the 3-metric γ ij . DefiningX = X √ γ for any tensor X, and considering now energy-integrated moments (i.e., all moments are integrated from ν = 0 to ν = ∞), we get the equations for a grey two-moment scheme in the simulation frame: (91) with K ij the extrinsic curvature of the underlying spacetime. These equations are particularly convenient to use in general relativistic simulations using finitedifference or finite-volume conservative methods for the evolution of the fluid equations. Indeed, the moment equations are now nearly identical to the ideal fluid equations without neutrino-matter interactions (in a fluid,P ij = √ γP γ ij , with P the fluid pressure). The only new terms are the source terms involving S µ . For well chosen closures (see below), these equations also form a wellposed system of hypebolic, causal equations, an important property in order to obtain stable numerical evolutions (Pons et al, 2000). We note that we have a separate system of equations for each species of neutrinos and antineutrinos.
The back-reaction of the neutrinos onto the fluid is easily computed: energymomentum conservation requires that the source terms transfering energy and momentum from the fluid to the neutrinos exactly cancel the source terms transfering energy and momentum from the neutrinos to the fluid. Our evolution equations are, for each neutrino species ν i , and thus with the sum being over all species of neutrinos and antineutrinos. General relativisitic merger simulations relying on the evolution of the moments of f ν have evolved either (Ẽ,F i ) while providing closures forP ij and the energy spectrum of the neutrinos; orẼ with a closure forF i ,P ij , and the energy spectrum. At this point, it is worth commenting on the use of grey schemes in general relativistic merger simulations so far. As we will see over the course of this section, the formalism for an energy-dependent moment scheme in general relativity has been fully developed, and has been used e.g. in core-collapse simulations (Kuroda et al, 2016;Roberts et al, 2016). Why is it not used in merger simulations? One issue with an energy-dependent scheme is of course that the moments of f ν within each energy bin have to be evolved, and thus the cost of the neutrino evolution scales at least linearly with the number of energy groups evolved. This would be a steep price to pay, even for simulations with relatively coarse energy resolution. More importantly, however, neutron star mergers and their post-merger remnants include regions where the fluid is moving at relativistic speed with respect to an observer at rest in the simulation frame, as well as steep velocity gradients. This creates two important issues. First, a neutrino propagating through the remnant may rapidly change energy group, if the energy discretization is done in the fluid frame. If the discretization is done in the simulation frame, on the other hand, transforming to fluid frame energies to calculate the source terms is non trivial, as that transformation depends on the unknown direction of propagation of the neutrinos (and neutrinos with the same energy in one frame may have very different energies in the other). Rapid variations of the neutrino energies can additionally lead to significant numerical diffusion in energy space and to a smoothing of the neutrino energy spectrum. Second, the flux in energy space F α ν = (νM αβγ ∇ γ u β ) can be large enough that explicit time stepping becomes unstable, at least for the time steps otherwise used for the evolution of the equations of general relativistic hydrodynamics. This means that either the time step has to be decreased (potentially drastically), or the energy flux has to be treated implicitly, thus coupling in an implicit time step the evolution of all energy groups. Either choice introduces a significant additional computational cost. This is not to say that an energy-dependent scheme in merger simulations is impossible. However, the development of a cost-effective and stable evolution scheme for an energy-dependent moment algorithm applicable to general relativistic merger simulations remains an unsolved problem, and thus any discussion of which scheme would be pratical is currently based on conjecture only.

Number density evolution and lepton number conservation
If we want more information about neutrino energies without going all the way to an energy-dependent scheme, a potentially useful extension of the standard grey two-moment approach described in the previous section is to define number-weighted moments in addition to the energy-weighted moments discussed so far, i.e., moments In this case, the first moment is independent of the choice of time vectort µ , as From the definition of N α , we see that in any given orthonormal frame, Nt is the number density of neutrinos, while Nî is the number flux of neutrinos. Taking the covariant divergence of N α and combining with Boltzmann's equation, we get This is simply expressing the conservation of neutrino number up to the contribution of the source term S N . If we write N α = N n α + F α N with F α N n α = 0, this is equivalent to Adding this evolution equation to the evolution ofẼ and, possibly,F i has two important advantages. The first is that one can get some information about the average energy of neutrinos from the local values ofẼ,F i ,Ñ . The other is that this equation allows for explicit conservation of all lepton numbers. For example, for the electron lepton number, Transforming to the variables typically used in fluid simulations, the electron fraction Y e and weighted energy density ρ * Y e = n e − − n e + n p + n n ; ρ * = m b (n p + n n )W √ γ with n i the number density of species i in the fluid frame and m b an arbitrary reference baryon mass, we get with v i T = α −1 v i − β i the transport velocity. This couples the evolution of the composition of the fluid to the evolution of electron type neutrinos in a way that guarantees conservation of electron lepton number. The main disadvantage, besides the cost of evolving an additional variable (which is in practice minimal), is that we now need a closure relation for the number fluxF i N , and still have to make semi-arbitrary assumption for the shape of the neutrino spectrum (Foucart et al, 2016b).
We note that whether the number densityÑ is evolved or not, we need to calculateS N . Indeed, capturing the evolution of Y e due to neutrino-matter interactions is one of the main objective of merger simulations including neutrino transport. Evaluating the source terms in Eq. (100) is thus a necessity in any transport scheme.

Source terms
As already mentioned, existing general relativistic simulations of neutron star mergers typically assume a source term of the form for neutrinos of energy ν in the fluid frame, with η and κ a calculated so that This assumes that emission is isotropic in the fluid frame, and guarantees that neutrinos reach their equilibrium energy density in optically thick regions. We note that this is not the most general form that the source terms can take, and in fact implicitly makes a number of important assumptions. Besides isotropic emission, this choice assumes isotropic elastic scattering, and it is not an accurate model for pair processes. We will discuss alternative methods to account for pair processes in optically thin regions in Sect. 4.2.8. In optically thick regions, we tend to rely on the fact that as long as neutrinos reach their equilbrium density on a time scale short compared to the dynamical timescale of the system, imposing Eq. (102) will be sufficient to approximately recover the correct physical behavior. We discuss the limits of this approach in Sect. 4.2.8 as well.
The coefficients η, κ a , κ s can be tabulated as functions of the fluid properties (ρ 0 , T, Y e ) and the neutrino energy ν (as in e.g. the NuLib library (O'Connor and Ott, 2010)). For grey schemes, we need instead the integrated emissivity and average opacities, defined such that with J tot and H tot the energy integrated energy density and flux. The calculation of the total emissivity is trivial: The calculation of the average opacities is however more complex. We would ideally want However, while we know κ a (ν), we only evolve J tot , not J(ν). Without information about the neutrino energy spectrum, the simplest choice aims to guarantee that the equilibrium neutrino energy density takes the expected value J eq tot = η tot c κ eq a (106) Unfortunately, this is only correct for the energy spectrum of neutrinos in equilibrium with the fluid. As it is quite common for neutrinos in low-density regions to have much higher energy than neutrinos in equilibrium with the fluid, and as κ a ∝ ν 2 for many reactions, this can lead to severe underestimates of κ a .
Alternatively, given an average neutrino energy ν , we may take advantage of the approximate scaling of opacities with ν to write κ a = κ eq a ν 2 ν eq 2 , with ν eq the average energy of neutrinos in equilibrium with the fluid, Estimates for the average neutrino energy have been taken from leakage predictions (Foucart et al, 2015), or from the evolution of N . For example, in Foucart et al (2016b) we used This last equation is an approximation (it ignores differences between the energy-weighted average energy of neutrinos and their flux-weighted average energy), yet it should provide a significantly better estimate of ν than the assumption of equilibrium with the fluid. The scattering opacity can use the same rescaling: with κ eq s calculated assuming an equilibrium spectrum of neutrinos. We note however that in using the same scaling for κ s and κ a , we again implicitly assume the same energy spectrum for J and H α . This is generally not true in the diffusion regime (see Sec.4.2.7).
We additionally need the source terms entering the evolution ofÑ and/or Y e . For neutrinos of a given energy ν, that source term is simply and, if we know η(ν), we can define an integrated number emissivity As before, we would like to define κ N such that Computing κ N requires us to once more guess at the neutrino energy spectrum. In Foucart et al (2016b) we chose κ N so that neutrinos properly thermalize to an equilibrium spectrum at the fluid temperature if the optical depth is large enough, i.e.
with ν calculated as before, and the last term included to take into account the energy dependence of the cross-sections. This choice is however far from unique. One can e.g. also choose with ν being estimated either assuming equilibrium between the neutrinos in the fluid (at the risk of strongly underestimating its value) or from a separate estimate of neutrino energies (e.g., from the energy of neutrinos according to a leakage scheme (Foucart et al, 2015)). We see that getting an accurate estimate of the source terms in a grey moment scheme is thus quite difficult, particularly when estimating the source terms entering into the evolution of Y e . As that is also one of the most important parameters to be impacted by neutrino-matter interactions, this may certainly be a significant issue limiting the accuracy of current merger simulations (see Foucart et al 2016b and the discussion of merger simulations here for the impact of inaccurate energy estimates).

Closures
The main approximations made by moment schemes are their choice of analytical closure. Grey schemes need an energy closure specifying the neutrino spectrum. Simulations evolving onlyẼ (1-moment schemes) additionally need a closure forF i andP ij , while simulations evolvingẼ andF i (2-moment schemes) need a closure forP ij . Energy-dependent schemes also have to specifỹ M αβγ .
There is no unique prescriptions for these closures that work in all possible regions of a simulation, but there are regimes in which they can be fairly easily calculated -mainly regions where neutrinos are in thermal equilibrium with the fluid (optically thick regime), as well as regions far away from a localized source of neutrinos, where all neutrinos approximately propagate away from that source.
In the first case, we can assume that for an orthonormal tetrad in the fluid rest frame Lîĵ ≈ 1 3 δîĵJ, as appropriate for a relativistic gas of neutrinos. In covariant form, this is Going back to an orthonormal tetrad in the fluid frame, we can write the moments equations in the optically thick regime as ∂tHî + 1 3 ∂îJ = −(κ a + κ s )Hî In this regime, neutrinos reach a quasi-equilibrium state on a time scale short compared to the time step of the simulation, and thus This gives us closure relations appropriate for a 1-moment scheme (Eq. (120)) and for a 2-moments scheme (Eq. (116)). We note however that these equations are provided in the fluid frame. The more useful relations providing F i or P ij as functions of lower order moments can be recovered using the relationship between fluid-frame and simulation-frame moments (see e.g., Shibata et al 2011;Cardall et al 2013;Foucart et al 2015). Applying this optically thick closure for all values of the opacity corresponds to the diffusion approximation (Flick's law). Such a choice is however clearly problematic for low value of the opacities, as the flux becomes infinite when (κ a + κ s ) → 0. For free-streaming neutrinos that all propagate in the same direction (as expected far from a post-merger remnant), we expect instead withf i the unit vector pointing away from the remnant, and Any form of Eq. (122) is a valid closure for a 2-moment scheme. Eq. (121) could on the other hand only be used in a 1-moment scheme if we choose a direction of propagation for the neutrinos, e.g., in a ray-by-ray transport scheme. For example, the M0 scheme of Radice et al (2016), which is only used to evolve free-streaming neutrinos in their mixed moment-leakage scheme, takesf i to be the unit vector in the t − r plane orthogonal to the fluid velocity. For 1-moment schemes, one way to combine these two limits is flux-limited diffusion (see e.g. Levermore and Pomraning (1981)), which sets for some function λ(R) which asymptotes to 1/3 in optically thick regions (small R) and 1 in optically thin regions (large R). For two-moment schemes, a number of proposals have been made for the "optimal" choice of closure, e.g. Wilson et al (1975); Minerbo (1978); Levermore and Pomraning (1981); Smit et al (2000), usually under the assumption of a preferred direction for the propogation of neutrinos (spherical or planar symmetry). If we write F i = f En i , the tensorial closure above can then be replaced by a prescription for the scalar Eddington factor p = P ij n j /E . Existing general relativistic simulations use the maximum entropy closure for p derived by Minerbo (1978) for photon transport, and updated by Cernohorsky et al (1989) for neutrino transport. A closed-form expression for this closure was derived by Cernohorsky and Bludman (1994). The full pressure tensor can be recovered from p using the prescription of Levermore (1984); Dubroca and Feugeas (1999). Overall, this results in the closure relation using the optically thick limit described earlier in this section (which is expressed in the fluid frame, as discussed below), as well as and We see that in optically thick regions (ξ = 0) we recover the optically thick closure P ij thick , while in free-streaming regions (ξ = 1), we get P ij thin . It is however worth noting that while the former is asymptotically correct for large opacities, the latter is not correct for free-streaming neutrinos: neutrinos in vacuum are generally not all propagating in the same direction, and thus do not satisfy Eq. (122). When making this choice, the interpolation between two well-defined asymptotic regime is thus not the only source of error: we have to come to terms with the fact that the optically thin regime uses a closure that is inaccurate in most regions where neutrino-matter interactions are important. A well known consequence of that choice is the creation of artificial radiation shocks whenever neutrino beams cross (see Fig. 4, and similar results in the two-beam test performed by Sadowski et al (2013)). In neutron star mergers, this also leads to an overdensity of neutrinos in the polar regions (see discussion of simulations).
We also note that in the optically thin limit, we listed a number of closures that were theoretically equivalent as long as F k F k = E 2 , a condition satisfied if all neutrinos propagate in the same direction. These expressions do however differ in regimes when F k F k = E 2 . Shibata et al (2011) show that the choice has the advantage to guarantee P k k = E, but can lead to superluminal characteristic speeds for the moment equations if used directly as a closure, i.e. for leads to P k k = E when F k F k = E 2 . Their recommendation to prefer Eq. (126), consistent with Levermore (1984); Dubroca and Feugeas (1999), has been followed in merger simulations so far, and appears to lead to a system of equations that is both hyperbolic and causal (Dubroca and Feugeas, 1999;Shibata et al, 2011).
As ξ is a function of the fluid-frame moments (J, H α ), which are themselves functions of the evolved simulation frame moments (E, F i , P ij ); and as P ij depends on ξ, we see that our prescription for ξ is an implicit equation. We can solve for ξ by searching for roots of the function with the known physical bounds ξ ∈ [0, 1]. As ξ changes relatively slowly, the value of ξ at the end of a recent time step is usually a good initial guess for the current value of ξ, allowing for rapid convergence of many root finding algorithms. For grey schemes evolving (Ẽ,F i ), the above closure combined with a choice of neutrino energy spectrum is sufficient to close the system of equations. When evolvingÑ as well, the situation is a little more complex, as we also need to estimate the number fluxF i N . For monoenergetic neutrinos, by definition, However, after integrating over neutrino energies, we should get The first term in this flux is important to capture advection of neutrinos with the fluid, while the second is important to capture neutrino propagation in the fluid frame. The factor ν J is an energy-weighted average neutrino energy that can reasonably be estimated using Eq. (109). On the other hand, ν H is a flux-weighted average energy that, in dense regions, is likely smaller than ν J . Indeed, low-energy neutrinos diffuse faster than high-energy neutrinos, and thus H α (ν) has a softer energy spectrum than J(ν). In fact, there is no particular reason to expect ν H to be the same for every component of H α ! Foucart et al (2016b) developed a rather complex procedure to estimate ν H during a simulation, which involves a new evolution equation for the spectral index of each neutrino species. While that procedure was calibrated to match the result of energy-dependent radiation transport in specific test cases, and aims to capture the transition from a soft spectrum for H α in the diffusion regime to a thermal spectrum close to the neutrinosphere, its accuracy for more complex physical configurations is unknown. Radice et al (2022), who also evolveÑ , make the simpler approximation ν H = ν J = ν . Regardless of the choice made, the difficulty of properly capturing the evolution ofÑ 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 in the context of general relativistic simulations, it is worth mentionning that Shibata et al (2011) and Cardall et al (2013) have provided closures for the third moment M αβγ . One method assumes that the expansion of the distribution function in l α 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 αβγ . These closures are not needed in existing grey moment schemes, but would be necessary for energy-dependent schemes.

Numerical implementation
In the previous sections, we covered most of the ingredients needed to evolve moments of f ν : the moment equations in the simulation frame, their coupling to the fluid equations, and the required analytical closures for higher-order moments of f ν (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 with U the vector of evolved variables (e.g. U = (Ẽ,F i ,Ñ ) 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 ∆t α CF L ∆x, with ∆x the grid spacing and α CF L 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 where S imp + S exp = S, and S 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 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, 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 Modern numerical simulations often use higher-order, less dissipative methods. In those advanced methods, values of F, U at (i, i + 1) are replaced with leftbiased and right-biased stencils used to interpolate (F, U ) on cell faces (van Leer, 1977;Jiang and Shu, 1996;Borges et al, 2008), 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 and Rezzolla 2012 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 (2011), and details of the calculation of F * for twomoment algorithms used in merger simulations are available in Shibata et al (2011);Foucart et al (2015); Weih et al (2020b); Radice et al (2022); Sun et al (2022). The only difficulty specific to the radiation transport equations is the treatment of high-density regions, discussed in more detail in Sect. 4.2.7. The treatment of the implicit terms has steadily improved over the years, starting from the hybrid leakage-moment schemes that do not require any implicit treatment of the source terms (Shibata et al, 2011), to first approximate (Foucart et al, 2015), and then full implicit-explicit time-stepping (Weih et al, 2020b) linearizing the source terms around the zero-state (E = F i = N = 0), to most recently a linearization of the problem around an arbitrary state for the neutrino radiation field (Radice et al, 2022). We review here the latest methods of Radice et al (2022). In that work, the Jacobian matrix is calculated explicitly for the Minerbo closure around an arbitrary state U . The implicit equations for U n+1 can then be solved using standard iterative methods for the determination of the roots of a multi-dimensional function f (U ), given f , ∂f /∂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 = (Ẽ,F x ,F y ,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 [J(U ) is independent ofÑ if U = (Ẽ,F i )]. As J(U ) is quite complex, we refer the reader to Radice et al (2022) for its exact form. Many older two-moment schemes however use more approximate values of the Jacobian matrix 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 fl and U rad are the vectors of evolved variables for the fluid and neutrinos respectively, the coupled problem is evolved using equations of the form 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 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 et al 2016b). 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 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 (Ẽ,F i ) and a single linear implicit equation forÑ . A full implicit solve of the fluid and radiation equations, on the other hand, would require solving a system of 6 + 5N ν non-linear implicit equations for the N ν neutrino species and the fluid (assuming 6 fluid variables and 5 components of the moments for a grey two-moment scheme). As long as the simpler method does not lead to numerical instabilities, its significantly lower computational cost makes it much more appealing.

Diffusion regime
In regions where (κ a + κ s )L 1, with L a typical lengthscale of our system, we expect the neutrino momentum density in the fluid frame to be ∝ −(κ a + κ s ) −1 ∇J. This leads to an evolution equation for the energy density in the fluid frame i.e., a diffusion equation with diffusion coefficient D = [3 (κ a + κ s )] −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. (134), and ignoring spatial variations in the opacities, the discretized equation becomes (Audit et al, 2002) ∂tJ − 1 3(κ a + κ s ) We clearly see that in any region where (κ a + κ s )∆x 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 (2002) suggest to effectively use a one-moment scheme in regions where (κ a + κ s )∆x 1, i.e., to replace the numerical flux in the evolution ofẼ by its value assuming that This leads to neutrinos being advected with the fluid, with an additional slow diffusion provided by H α . A more detailed discussion of this method in the merger context can be found in Foucart et al (2015). As pointed out in Radice et al (2022), 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 (2022) 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 (κ a +κ s )∆x 1. In both algorithms, one transitions smoothly between the standard shock-capturing fluxes and the modified fluxes around (κ a + κ s )∆x ∼ 1. For example, Radice et al (2022) use F * = (1 − a)F HO + aF LO (142) with F LO given by Eq. (134) (replacing c by the value of the speed of light in the simulation coordinates) and F HO using the simple second-order accurate prescription The transition coefficient is and the opacities are estimated using the average values of neighboring grid points. Outside of the merger context, a conceptually similar correction limiting the use of the dissipative fluxes in high-opacity regions had been used by Sadowski et al (2013). Radice et al (2022) show that the approximate linearization of the source terms used in many existing two-moment schemes and 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 κ s ∆x 1, v/c ∼ 1, and κ a ∆x 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 postmerger 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.

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 Sect. 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 ν µνµ pairs and ν τντ 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. Deep into the neutron star, heavy lepton neutrinos are nonetheless in thermal equilibrium with the fluid, at least as long as κ a t diff ∼ 1. Due to the slow diffusion rate through the dense matter, this remains true even in regions where their absorption optical depth τ a < 1, as their scattering optical depth τ s 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, we eventually reach regions where κ a t diff < 1 and τ s > 1. There, heavy lepton neutrinos slowly diffuse out while out of thermal equilibrium with the fluid (typically, neutrinos will be hotter than if they were in equilibrium at the fluid temperature). 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, ν e andν 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: νν 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 Eq. (30), we get that the appropriate energy-integrated source term for the moment equations is, for annihilation into e + e − pairs and for the energy deposited by the neutrinos only, The integral over p (ν) , 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 et al (2017) propose the approximation with (ν) some appropriate average of the energy of annihilated neutrinos. Given the scaling in the source term, the choice would be reasonable for a quasi-thermal spectrum at an estimated temperature T (ν) . This remains a significant approximation. Indeed, by using S α ∝ u α , 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  indicate that S α 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.

Discussion
As the main algorithms used so far for radiation transport in general relativistic 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 et al, 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 for estimating neutrino energies provide wildly different answers for the emission of heavylepton neutrinos ) -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 et al, 2022). 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 et al, 2016b). 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 pressure closure chosen in existing simulations is also a source of concern. Factor of ∼ 2 errors in the energy density of neutrinos in the polar regions observed in comparisons with Monte-Carlo simulations ) 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 of a 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 et al, 2017;Foucart et al, 2018) (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 , with disagreements at the 10% − 20% level in most global quantities, including the ν e andν 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).

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  and post-merger (Miller et al, 2019b) simulations.
As for moment schemes, the use of Monte-Carlo methods for the evolution of radiation coupled to a fluid has a long history outside of the context of neutron star mergers and general relativistic radiation hydrodynamics. Monte Carlo methods have been developed to, among other applications, study stellar profiles (Lucy, 1999), stellar formation (Ercolano et al, 2003;Haworth and Harries, 2012), black hole accretion (Ryan et al, 2015), post-merger accretion disks (Richers et al, 2015), supernova ejecta (Kasen et al, 2006;Noebauer et al, 2012;Wollaeger et al, 2013), or core-collapse supernovae (Janka, 1991;Abdikamalov et al, 2012). In this section, we will first review the formalism of Monte-Carlo radiation transport in general relativistic simulations (Sect. 4.3.1). We will then discuss technical issues related to the coupling of the neutrinos to the fluid (Sect. 4.3.2), the treatment of optically thick regions (Sect. 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. We base our discussion of general relativistic Monte-Carlo transport in large part on the methods developed for photon transport in black hole accretion disk by Ryan et al (2015), with modifications made for applications to the merger problem in Foucart et al (2021). We also comment on the recent development of another Monte-Carlo code aimed at axisymmetric post-merger simulations by Kawaguchi et al (2022). That code generally makes less simplifying assumptions and aims for higher order methods than Foucart et al (2021), taking advantage of the expected lower cost of axisymmetric simulations.
As in the previous section, we assume h = c = 1 unless noted otherwise.

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.
with P the ensemble of all packets, x i k (t) the position of packet k as a function of time, and p k i (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 (ν) .
The finite number of packets in a Monte-Carlo simulation will inevitably lead to sampling noise in f (ν) , and in any variable derived from f (ν) . We typically expect relative errors in quantities obtained by summing over N packets to be ∼ N −1/2 . If one wanted ∼ 1% errors for the energy density of neutrinos within each cell of a computational simulation at any given time, one would thus need ∼ 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 each individual packet in a given grid cell of coordinate volume ∆V represents a total neutrino energy E p , and that the energy-integrated emissivity is η tot , we create packets over a time step ∆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 their energy from the distribution f (ν) = η(ν)/η tot . The initial time and position of a packet can be drawn from an homogeneous distribution within the 4-volume (∆V ∆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 ∝ 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 η(ν) within specific energy bins, rather than a known continuous function η(ν). 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 values of κ a at the center of the same energy bins, the equilibrium energy density of neutrinos will match its expeceted physical value to numerical roundoff (Richers et al, 2015).
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 et al (2021). We note however that many other choices of E p are possible, and that the optimal choice is problem dependent. Kawaguchi et al (2022) instead propose to choose E p as a fraction of the energy of neutrinos in thermal equilibrium with the fluid in a cell, with a prescription to destroy and resample low energy packets at the end of each step to avoid the continuous evolution of a large number of low-energy packets in optically thin regions. Simulation of accretion disks around black holes can rely on simpler prescriptions for E p (e.g. constant packet energy), as they do not need to deal with the dense, hot regions observed in neutron star remnants.
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 (1994) 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: The first two lines are evolution equations for x i , p i , while the third comes from the constraint that p µ p µ = 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 et al 2021).
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 λ i , or opacity κ i = λ −1 i , then for each process we can randomly draw the time to the next interaction from a Poisson distribution, e.g.
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 ∆τ is in the reference frame in which we provide κ i . In simulations, this is usually the fluid frame. Transforming to the coordinate time, we get (Ryan et al, 2015) ∆t The time to the next event is then ∆t next = min (∆t i , ∆t step ) with ∆t step the time step, and with the minimum taken over all possible processes. If the minimum is ∆t step , a packet is simply propagated for that time. Otherwise, whichever process provided the minimum ∆t occurs after that time interval. Existing simulations have considered an absorption opacity κ a and an isotropic elastic scattering opacity κ 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 results 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 ∆t i is only valid for a constant κ 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 κ i within each grid cell , which is numerically convenient but may need to be improved. 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 et al (2021), 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 are desirable, but complicate calculations of the times-to-interaction ∆t i and, to be consistent, also require the use of higher-order estimates of the emissivities, and thus inhomogeneous production of packets within a single grid cell. A second-order accurate scheme has recently been developed by Kawaguchi et al (2022), for applications to axisymmetric systems, where the additional cost of a higher-order scheme is more manageable than in full 3D merger simulations.
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 (Sect. 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 κ a ∆t step 1), and/or slowly diffuse through the fluid (if κ s ∆t step 1). We discuss how one may deal with these issues in Sect. 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 ν : Before getting into these issues, a useful additional piece of formalism to discuss is the computation of moments of f (ν) 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 δ in the definition of f (ν) can be practically problematic when considering moments of f (ν) , and it is thus often preferable to consider the average moments within a volume V (e.g. a grid cell), defined as For the second moment, this is and similary for the number flux 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 with ν k = −p µ k u µ . 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.

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 some 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 ∆p µ k = p µ after − p µ before , then the fluid should undergo a change of 4-momentum −N k ∆p µ k as a reaction. While this change is generally instantaneous, we can distribute it over a 4-volume element ∆V ∆t to get the 4-force density (Ryan et al, 2015) 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 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 ∆t are The most commonly evolved fluid variables in merger simulations are the internal energy densityτ =T µν n µ n ν − ρ * and momentum densityS i = −T µ i n µ . For these variables, we get Similarly, for the evolution of the electron fraction, with s k = 0 for muon and tau neutrinos and for any scattering event, s k = 1 for emission of ν e and absorption ofν e , and s k = −1 for emission ofν e and absorption of ν 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. It is best suited to simulations with fairly homegeneous packet energies, or a sufficiently large number of packets to avoid shot noise issues. This is the method adopted e.g. in Miller et al (2019a). 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 et al (2021), 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, As when using exact conservation laws, we then need to integrate this source term over a time step ∆t and a small volume ∆V , as well as over the entire neutrino energy spectrum. We then get In this expression, the sum is taken over all packets k, and all time intervals ∆t k,j during which packet k is propagating along a geodesic while inside of the volume ∆V and time interval ∆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 ∆t, we get For the electron fraction, we get instead with s k = 1 for ν e , s k = −1 forν e , and s k = 0 otherwise. With these choices, the source terms in low-density regions (κ∆t 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 et al, 2021). 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 ∆t k,j .
An intermediate scheme, used e.g. by Kawaguchi et al (2022), is to absorb full packets in optically thick regions, but to damp N k in optically thin regions, following the method of Dolence et al (2009). The number of neutrinos in a packet then evolves as with τ the time in the fluid frame. In this case, the energy / lepton number absorbed over a time step can easily be added to the fluid while guaranteeing that we exactly satisfy the relevant conservation laws, but one may end up with a larger number of packets as individual packets are no longer destroyed by default (an issue solved in Kawaguchi et al (2022) by resampling packets in optically thin regions at the end of a step anyway). 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, and 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.

Optically thick regime
The optically thick regime poses a more fundamental problem for Monte-Carlo methods. One issue is that if (κ a +κ s )∆t = N , we expect ∼ N neutrino-matter interactions per time-step. In dense or hot regions of mergers, we can easily get N 1, making evolutions costly. A related issue is that if we calculate η, κ a , κ 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 (2019a) 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 (2020), and handled with approximate methods inspired by the implicit Monte-Carlo method of Fleck and Cummings (1971) for regions of high absorption optical depth, and the random walk method of Fleck and Canfield (1984) for regions of high scattering optical depth. We briefly summarize their methods here, as well as the proposed algorithm of Kawaguchi et al (2022), 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 κ s ∆t 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 (2017) andFoucart (2018) propose similar methods to treat this regime, based on the work of Fleck and Canfield (1984). Both assume that in regions of sufficiently high κ s ∆t, neutrinos are advected with the fluid, while undergoing a slow random walk away from the fluid's motion. When determining the outcome of that random walk, Richers et al (2017) draw from the distribution of times needed for a neutrino to diffuse a certain distance in the fluid frame. Foucart (2018) draws instead 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 et al (2017);Foucart (2018) 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 ∆p µ of each packet between the beginning and end of a step, but should in theory correct that result for the natural change in p µ due to the evolution of a packet along a geodesic in curved spacetime. If using expectation values of the interaction rates instead, Foucart (2018) uses the approximation that p µ = νu µ for a fraction f adv of the time step, and p µ = ν(u µ + l µ ) for the remaining of the time step, with f adv and l µ 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 κ s ∆t ∼ 3 (Foucart, 2018).
Regions with κ a ∆t 1 create a more difficult challenge. In these regions, we expect neutrinos to reach statistical equilibrium with the fluid on a time scale κ −1 a ∆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 η/κ a , and their diffusion timescale to be set by (κ a + κ s ) −1 . Then, the evolution of the system can be approximately captured if we accept that we cannot resolve what happens on timescales much shorter than ∆t (Foucart et al, 2021), and rely instead on changes to the source terms inspired by Implicit Monte-Carlo methods (Fleck and Cummings, 1971). In Foucart et al (2021), we propose the transformation which guarantees that η /κ a = η/κ a and κ a + κ s = κ a + κ s , but modifies the equilibration time scale of neutrinos from κ −1 a to (aκ a ) −1 . We then choose a such that κ a ∆t 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 subpercent errors in the neutrino luminosities when comparing this method to the solution of Boltzmann's equation (Foucart et al, 2021). 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 and Cummings (1971), 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 and potentially reduce the errors that it creates.
In their proposed Monte-Carlo algorithm, Kawaguchi et al (2022) adopt an algorithm much closer to that of Fleck and Cummings (1971). A single parameter a is chosen for all energy bins of a given neutrino species, and the additional scattering opacity induced by the algorithm is assumed to be inelastic: the total energy of scattered packets does not change, but the energy of individual neutrinos within the packet samples the equilibrium energy distribution of neutrinos. With such an algorithm, Implicit Monte-Carlo can be seen as just a specific discretization of the original transport equations (Fleck and Cummings, 1971), a significant theoretical advantage. This comes at the cost of slightly less flexibility in our ability to reduce absorption and reemission of packets, and more complexity in the treatment of scattering events.

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 ν iνi → e + e − reaction. The absorption cross section for neutrinos of 4-momentum p α could be written, under the same assumptions as Eq. (30) and assuming that each packet of antineutrinos represent a spatially uniform distribution of neutrinos within a grid cell of volume V , as withp α k 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 ∆t νν = κ −1 νν ln r with r ∈ [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 withT αβ the stress-energy tensor of the antineutrinos. This allows for faster calculations for large numbers of antineutrinos if we precomputeT αβ (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 approximateT αβ 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 contractions with the metric, a code that stores p t and p i would ideally express this opacity as and thus specifically precomputeĒ,F i ,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.

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 ν . Monte-Carlo simulations have already allowed important tests 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 a 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 combined 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 ν . 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 ν 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 ν e minus flux ofν 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.

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.

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 erg/s. That emission however rapidly decays on ∼ 10 ms timescale (Deaton et al, 2013;Sekiguchi et al, 2016;Radice et al, 2022). 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 (Fernández et al, 2020;Shibata et al, 2021;Fujibayashi et al, 2022). It may remain as high as ∼ 10 52−53 erg/s for O(100 ms). This emission is not sufficient to create a thin disk (the disk aspect ratio remains ∼ 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Ṁ 10 −3 M /s (De and Siegel, 2021)), 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 (2020) find neutrino luminosities of ∼ 10 52−53 erg/s more than 5 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 ms) (see Fig. 6).
In both cases, simulations generally agree on an energy hierarchy νx > νe > νe , with average energy of ∼ 10 MeV for ν e but potentially above 20 MeV for ν x . This is simply due to the higher absoption opacity of the fluid to ν e , which puts the ν e neutrinosphere farther out in the remnant than theν e neutrinosphere (and even more the ν 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ν e emission over ν 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. Long simulations of accretion disks with mild electron degeneracy using a leakage scheme have found that weak interactions can lead to self-regulation of the composition of the disk midplane to Y e ∼ 0.1 (Siegel and Metzger, 2018), with a lower density outer disk and hotter at significantly higher Y e (Fig. 7). Simulations of accretion disks with moment transport find that neutrino absorption in the disk can lead to higher electron fractions Y e ∼ 0.15 − 0.2 Foucart et al (2015); Just et al (2021). These values can also be significantly impacted by irradiation of the disk by a central neutron star (Fig. 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 (Fernández et al, 2017), 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 α viscosity models as a subgrid model to capture angular momentum transport and heating from magneticallydriven turbulence (Shakura and Sunyaev, 1973), 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. MHD simulations starting with different magnetic field configurations can also produce very different amounts of outflow,indicating that predictions for matter outflows in post-merger remnants are likely to depend on the unknown large scale properties of magnetic fields in the post-merger remnant (Christie et al, 2019;. 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 ∼ (0.2 − 0.4)Wanajo et al (2014) (2022), 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 0.05), even in simulations that include neutrino absorption (Sekiguchi et al, 2016;Foucart et al, 2017;Radice et al, 2018;Camilletti et al, 2022). 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 , but fairly low mass M outflow ≈ 10 −(3−4) M (Just et al, 2015;Foucart et al, 2015) (see e.g. polar regions of Fig. 7). They are likely to be dominated by magneticallydriven winds (Siegel and Metzger, 2018;Fernández et al, 2019), if those winds are significant (see above), and by viscous outflows in the advection-dominated phase (Fernández and Metzger, 2013). In the absence of strong magneticallydriven winds, neutrino-driven winds could however impact kilonova signals, as they might be geometrically separated from both the dynamical ejecta and the viscously-driven winds. Magnetically-driven outflows are faster (∼ 0.1c), and have less time to absorb neutrinos, while viscous outflows are slower (∼ 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 ∼ 0.15 − 0.25, starting from very neutron-rich initial conditions (Siegel and Metzger, 2018). Shorter simulations using an advanced Monte-Carlo transport scheme and including magnetic fields (Miller et al, 2019b), with similar initial conditions, find a broader distribution of Y e peaking just below Y e ∼ 0.2 and extending up to Y e ∼ 0.4. Similarly, long simulations of a NSBH merger remnant with a moment scheme and magnetic fields , initialized from the outcome of a merger simulation, find a broad Y e ∼ 0.15−0.4 distribution peaking just above Y e ∼ 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 (Fernández et al, 2020), the lifetime of the massive neutron star remnant (if present) (Metzger and Fernández, 2014), and the initial composition used in the simulation (Fernández et al, 2017). 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 et al, 2019) and a neutron star remnant (de Haas et al, 2022). 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 νν pair annihilation in relativistic merger simulations. Fujibayashi et al (2017), using the approximate moment method described in Sect. 4.2.8, finds that pair annihilation can accelerate the matter in the polar regions to mildly relativistic speeds (Lorentz factor Γ ∼ 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). In NSBH remnants with very little matter in the polar regions, Just et al (2016) even find relativistic outflows powered by pair annihilation, but not with enough energy to power the brightest SGRBs. Foucart et al (2018) find that the energy deposition predicted by the approximate moment scheme is likely accurate within a factor of ∼ 2-3. More advanced transport calculations had earlier been performed on a prescribed fluid background (Popham et al, 1999;Kneller et al, 2006;Zalamea and Beloborodov, 2009;Dessart et al, 2009;Perego et al, 2017). 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.

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 twomoment scheme (with Minerbo closure, evolvingẼ andF i ) was performed in the context of a low-mass neutron star merger in Foucart et al (2016a). In that simulation, leakage and moment schemes agreed quite well on theν e luminosity, while other luminosities varied by factors of 2-3 in the first 10 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Ṁ ∼ 10 −2 M /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 ( Y e = 0.2 vs Y e = 0.1, and s = 20k B vs s = 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 et al (2022), 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 ∼ 2 with respect to the moment simulations, 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 ν eν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 ∼ 30% less mass ejected than the hybrid scheme, and significantly higher electron fractions (∆Y e ∼ 0.1). The use of the Eddington or Minerbo closure had again a much smaller effect (∼ 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.
The impact of the chosen energy closure in grey two-moment schemes was investigated in Foucart et al (2016b), for the same physical configuration as in Foucart et al (2016a). 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 ν = max (T leak , T ), with T 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 ν , 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 Results are shown for three algorithms: a hybrid leakage-moment scheme (blue), a twomoment scheme with Minerbo closure (orange), and a two-moment scheme with Eddington closure (green). Figures reproduced from Radice et al (2022).
∼ 50% at the end of the evolution (8 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 et al (2022), 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 ∆Y e ∼ (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 when opacities are computed 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 et al (2022) 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 twomoment scheme evolving the number density were performed without backreaction of the Monte-Carlo code to the fluid (i.e., using the same fluid evolution as in the moment evolution)  as well as with coupling to the fluid  in a neutron star merger simulation. In the fully coupled simulation, the luminosity of ν e andν e and the energy of neutrinos was in better agreement (10%-20% differences) than in other comparisons between transport methods discussed so far, with the exception of the comparison between different choices of closure relations performed by Radice et al (2022). The heavy-lepton luminosity showed ∼ 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 than in other comparisons discussed in this section, with the average Y e changing by ∆Y e ∼ 0.03, and the maximum Y e by ∆Y e ∼ 0.05 -not necessarily negligible changes, but a significantly better agreement than in other comparisons (see Fig. 9). 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 ν 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 ∼ 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 Sect. 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. 10). Overall, we thus note that there has been a significant improvement in the magnitude of those errors, with the more modern moment (Foucart et al, 2016b;Radice et al, 2022) and Monte-Carlo  schemes getting estimated relative errors in the ∼ 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).

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 et al, 2022), and dependent on the algorithms used for neutrino evolution in the simulations used to calibrate them (Nedora et al, 2022). 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.