A practical guide to a moment approach for neutrino transport in numerical relativity

The development of a neutrino moment based radiative-transfer code to simulate binary neutron-star mergers can easily become an obstacle path because of the numerous ways in which the solution of the equations may fail. We describe the implementation of the grey M1 scheme in our fully general-relativistic magnetohydrodynamics code and detail those choices and strategies that could lead either to a robust scheme or to a series of failures. In addition, we present new tests designed to show the consistency and accuracy of our code in conditions that are similar to realistic merging conditions and introduce a new, publicly available, benchmark based on the head-on collision of two neutron stars. This test, which is computationally less expensive than a complete merging binary but has all the potential pitfalls of the full scenario, can be used to compare future implementations of M1 schemes with the one presented here.


INTRODUCTION
The merger of binary neutron-star (BNS) system offers a window into some of the most extreme physics in the universe (Baiotti & Rezzolla 2017;Paschalidis 2017).They are linked to short gammaray bursts (Rezzolla et al. 2011;Palenzuela et al. 2013;Kiuchi et al. 2015a), to the synthesis of heavy elements through rapid neutron capture processes (Metzger et al. 2010) and, of course, to the emission of gravitational waves (GWs).As such, understanding in detail the dynamics of these systems is a key challenge that the scientific community has been addressing with increasing levels of sophistication for over a decade.The first observation of GWs from a BNS merger by the LIGO-VIRGO collaboration (The LIGO Scientific Collaboration & The Virgo Collaboration 2017), together with a rich set of counterparts in the electromagnetic spectrum (The LIGO Scientific Collaboration et al. 2017;LIGO Scientific Collaboration et al. 2017), has provided constraints on gravity (Abbott et al. 2019) as well as precious information about the properties of matter at the extreme density regimes realised in the cores of neutron stars (see, e.g., Abbott et al. 2018;Bauswein et al. 2017;Margalit & Metzger 2017;Rezzolla et al. 2018;Most et al. 2018).Future observations of GWs from late BNS inspirals and early post-merger emission by the hypermassive neutron star (HMNS) remnant by current and futuregeneration detectors (Punturo et al. 2010) will further constrain the equation of state (EOS) of matter above nuclear saturation density, thus probing a regime where the properties of the QCD phase diagram cannot be easily explored by standard techniques.However, GWs from the inspiral and post-merger of a BNS system alone cannot provide us with the all the details of this complex picture.Particu-★ E-mail: musolno@itp.uni-frankfurt.de(KTS) larly important are the explanation and accurate modelling of key phenomena associated with these events, such as the launching of a relativistic jet (Rezzolla et al. 2011;Just et al. 2016;Ciolfi 2020;Hayashi et al. 2022), the composition and properties of the ejecta leading to the emission of a kilonova (Bovard et al. 2017;Papenfort et al. 2018;Combi & Siegel 2022), and the long-timescale evolution of the system beyond the gravity-dominated regime.All of these questions need to be addressed both for their astrophysical interest as well as for the complementary information they may provide about the EOS of neutron-star matter (Combi & Siegel 2022;Papenfort et al. 2022).
Modelling the system self-consistently and over timescales that are hundreds or thousand times the dynamical timescale proves to be a formidable challenge both from the theoretical as well as from the computational perspective: if the inspiral is mostly dictated by gravity and the strong interactions between the constituent particles of the NS fluid, considerably more complex physical phenomena become relevant at the merger.First, it was shown that even small magnetic fields in the two constituent NSs will amplify due to turbulence and Kelvin-Helmholtz instabilities (Kiuchi et al. 2015b;Aguilera-Miret et al. 2022;Chabanov et al. 2023) to reach almost equipartition with the fluid kinetic energy, making electromagnetic forces extremely relevant on the long timescale in which the remnant evolves.Moreover, weak interactions play a crucial role in the dynamics of the system after merger for a number reasons.First, they contribute to reprocessing the ejected material increasing its electron fraction (Combi & Siegel 2022;Radice et al. 2016;Sekiguchi et al. 2016).Second, they provide a pressure contribution that will lead to the shedding of additional reprocessed material in the form of a neutrino wind (Dessart et al. 2009;Perego et al. 2014;Fujibayashi et al. 2017).Third, on the timescale of the diffusion of neutrinos through the optically-thick remnant, they cool the system very efficiently, thus counteracting the heating due to magnetic-field turbulence.Finally, it has long been suggested that neutrinos may have a role in the powering of a collimated polar outflow from the merger remnant (Eichler et al. 1989;Ruffert et al. 1997;Just et al. 2016).
Several approaches are available in the literature for the inclusion of the effects of neutrinos in general-relativistic hydrodynamical or magnetohydrodynamical (GRMHD) simulations of BNS mergers.These range from very simple and computationally efficient "leakage-type" schemes (Ruffert et al. 1997;Galeazzi et al. 2013;Most et al. 2019), where the local heating/cooling rates are directly estimated from the reaction cross-sections corrected with a diffusion prescription, over to the so-called "moment schemes", where a varying number of moments of the Boltzmann equation for neutrinos is solved (Rezzolla & Miller 1994;Foucart et al. 2015;Just et al. 2015;Kuroda et al. 2016;Skinner et al. 2019;Melon Fuksman & Mignone 2019;Weih et al. 2020a;Radice et al. 2022;Sun et al. 2022;Izquierdo et al. 2022).The most advanced approaches even consider the direct solution of the radiative transfer equation via Montecarlo or other methods (Radice et al. 2013;Foucart et al. 2020;Roth et al. 2022).
In this paper we present FIL-M1: a new implementation of a twomoment (M1) scheme for general-relativistic radiative transfer applied to BNS mergers.Obviously, this is not the first paper discussing such an implementation and papers with a similar spirit have been published very recently (see, e.g., Weih et al. 2020a;Radice et al. 2022).Here, however, besides describing in detail our implementation of the grey M1 scheme in a fully general-relativistic MHD, we also concentrate on highlighting the possible pitfalls one may encounter in the process of developing a similar code.In addition, we present new tests designed to show the consistency and accuracy of our code in conditions that are similar to realistic merging conditions: the head-on collision of two neutron stars.The neutrino luminosities obtained in this way and the evolution of the maximum temperature and density can serve therefore as an effective benchmark against which test new implementations.
The content of the paper is organised as follows: in Sec. 2 we review the theoretical background of moment-based methods for radiative transfer and the basic equations involved.In Sec. 3 we describe in detail the numerical scheme implemented in our code, paying special attention to troublesome aspects of the algorithm which are key to a successful application of the methods.Finally, in Sec. 4 we present results from standard implementation tests, as well as from a set of realistic tests designed to probe the accuracy and convergence properties of the code when presented with a realistic neutron-star simulation.

M1 SCHEME FOR GENERAL RELATIVISTIC RADIATION TRANSFER
In the limit where neutrinos are considered to be massless particles, their evolution can be expressed in terms of the radiation intensity, which in turn has to satisfy the Boltzmann equation (see, e.g., Rezzolla & Zanotti 2013).Since this equation is intrinsically 3 + 3 + 1 dimensional, its direct solution is computationally too expensive (see however Weih et al. 2020b, for some efficient methods to discretize the Boltzmann equation) to be feasible in conjunction with the solution of the coupled Einstein and GRMHD system.Following the moment formalism by Thorne (Thorne 1981), the M1 scheme consists in evolving the zeroth and first moments of the radiation intensity whilst providing an approximate closure for the second.In the fol-lowing, we will assume that spacetime is coupled with an "ordinary" perfect fluid described by the equations of GRMHD which couples to neutrinos by weak processes to be defined below.To this end, we consider the second moment of the radiation transfer equation which reads (see Shibata et al. 2011, for a detailed derivation) where   (  ) is the radiation intensity integrated twice over the (null) direction of propagation,   is the fluid four-velocity,  is the frequency of the neutrino, and   (  ) is the moment of the source term.For simplicity, we will here consider a so-called "grey scheme", i.e., one where we evolve the moments integrated over the neutrino frequency (energy)  .As a result, the second term in Eq. ( 1) vanishes out and we are left with a conservation equation for the radiation energy momentum (2) Obviously, the total energy momentum tensor, which is the sum of the energy-momentum tensor of the ordinary fluid   flu and of the radiation fluid   rad will have to be conserved, so that (3) In the expression above, and for a magnetised fluid in the ideal-MHD limit, the energy-momentum tensor reads where  and ℎ are the rest-mass density and specific enthalpy,   is the fluid's four-velocity,   the four-metric and   the magnetic field four-vector in the frame comoving with the fluid, with  2 :=     .On the other hand, the energy-momentum tensor of neutrinos can be expressed in terms of the first three moments of the radiation intensity in a frame comoving with the fluid The M1 scheme has the advantage of providing immediate physical interpretation of the quantities at hand: in analogy with the theory of ideal fluids we identify  as the energy density,   as the energy flux density and   as the pressure tensor of the radiation as measured in the fluid comoving frame.It is, however, worth noting that interpreting the neutrino field as a ideal fluid only makes sense in completely optically-thick media where the radiation is in local thermodynamic equilibrium with the background matter (see discussion in Rezzolla & Zanotti 2013).To obtain a set of conservation equations suitable for numerical integration, and as is customary for ordinary GRMHD, we can make use of a 3 + 1 foliation of spacetime to project the radiation energy momentum tensor onto (and in the direction orthogonal to) a spacelike hypersurface, which yields Here,  is the radiation energy density and, by construction,   is a spatial vector and    is a symmetric rank-two tensor lying in the tangent space to the spatial hypersurface.The equations for the grey M1 scheme finally read (Shibata et al. 2011) where, as customary,   is the shift vector,  the lapse,    the spatial three-metric, and  its determinant.We also define the Lorentz factor  associated with the fluid four-velocity   as  :=   .Hence, in addition to the standard GRMHD equations in our code (Most et al. 2019), we evolve equations ( 7) and ( 8) for three effective neutrino species: the electron neutrinos   , the electron antineutrinos ν , and an effective heavy neutrino species   .While the approximation of evolving the energy integrated moments (grey) is crucial to attain the desired computational efficiency, it foregoes any information regarding the neutrino energy spectrum, which might be extremely important to correctly capture the composition of matter in the aftermath of a BNS merger.Moreover, the set of Eqs. ( 7), (8) does not guarantee the conservation of the total lepton number.To partially overcome these issues and following Foucart et al. (2016); Radice et al. (2022), we evolve an additional equation for the neutrino number density.In particular, we introduce an effective neutrino number current   and assume that the flow of neutrino number comoves with the neutrino energy flux, i.e., where we have defined the neutrino number density in the fluid frame as ñ.Defining then with the last equality following from the fact that   is orthogonal to   (    = 0), we can write the evolution equation for the neutrino number density as where the source term N will be presented below [see Eq. ( 15)] and   is the projection of   in the tangent space of the three-hypersurface Note also that in writing Eq. ( 11) we have made the implicit assumption that the neutrinos have the same average energy in the advection term as in the energy-flux term.This is a valid assumption if all the neutrinos have the same energy and it is not strictly true in general; however, it represents a reasonable approximation in most of the cases we are interested in.
We now turn our attention to the source terms, which provide the coupling between the radiation and the ordinary fluid.These terms will model phenomena of isotropic, elastic scattering of neutrinos with the fluid, absorption and the emission of radiation.In particular, we can write the source terms as so that, in the Eulerian frame, this yields which will serve as coupling to the energy and momentum conservation equations.The RHS of Eq. ( 11), which will be coupled to the fluid's electron fraction, reads In Eqs. ( 14), (15) above,  er ( nr ) refer to the total energy (number) emissivity of a given neutrino species, while  er  ( nr  ) are the opacities due to energy (number) absorption and  er  ( nr  ) the opacities due to scattering; we additionally indicate the total opacity by  nr,er :=  nr,er  +  nr,er  . (16) In the following sections we will often forego the superscript and from the context it will be understood whether we are referring to energy or the number rates.Note that, as written, Eqs. ( 7) and ( 8) are exact, but the corresponding system cannot be solved because it lacks an evolution equation of the rank-two tensor    .This is not unusual in a moment-based scheme where the the solution of the lower-order moments requires the knowledge of higher-order ones and the system is complete only with an infinite expansion.The standard solution to this otherwise unavoidable problem is via the introduction of a "closure relation", that is, an equation relating higher-order moments with the lower ones.In our case this amounts to prescribing a mathematically well behaved and physically well-motivated relation for the pressure tensor of the type In deciding for such a closure we are aided by a large literature that has experimented on prescriptions that allow us to smoothly express the pressure tensor in the limit of an optically-thick medium, where the radiation is in thermodynamic equilibrium with the fluid, and in the optically-thin limit for radiation emitted by a point-like source.
We discuss below the details of the closures implemented in FIL-M1.

Closure Relations
As customary in implementations of the M1 scheme, we here adopt the following ansatz for computing the neutrino pressure tensor which is expressed as a linear weighted average between the expression it assumes in the optically-thin and optically-thick regimes, namely where FIL-M1 employs Minerbo's (Minerbo 1978) closure by default, for which with  being the variable Eddington factor defined as which is zero in the optically-thick limit and one in the optically-thin one.Note that ( = 0) = 1/3 and ( = 1) = 1 and that Eq. ( 20) requires the knowledge of the fluid frame moments of the radiation field since   and  appear in Eq. ( 20) as opposed to   and .
This is due to the fact that   is not guaranteed to be small in the optically-thick limit (Shibata et al. 2011).
Because physical quantities related to the ordinary fluid or to the radiation fluid are meaningfully expressed in the locally comoving frame, but the actual evolution equations are expressed in terms of an Eulerian frame, it is important to express the Eulerian-frame moments (we recall that these are the radiation energy density  and its flux   ) in terms of the fluid-frame ones Using again the fact that   is orthogonal to   , we can contract Eq. ( 22) with  and find from which we can obtain an expression for the projection of   along the unit normal Finally, we can compute the pressure tensor as Specialising now to the optically-thick case, for which we can assume the pressure tensor to be isotropic in the fluid rest-frame we then use Eqs.( 21) and ( 22) to write and Eq. ( 24) to express     in Eq. ( 27) to obtain We then employ Eqs. ( 28) and ( 29) to find from which, using Eqs.( 25) and ( 26), we finally obtain the desired expression for the radiation-pressure tensor in the optically-thick limit Similarly, when dealing with the optically-thin case, we can express the pressure tensor in the Eulerian frame for a point-like emitter as where we introduced the normalised energy flux Using expression (32), we can easily compute the fluid-frame moments and We are now in a position to write down the fluid frame moments in terms of the Eulerian frame ones.From Eqs. ( 29), (34) we then obtain and from Eqs. ( 30), (35), and ( 24) where and 0  := − , (47) We now have all the information needed to compute the closure.As can be seen from Eqs. ( 36) and (37), the computation of  is a nonlinear problem, where the knowledge of the evolved moments ,   is not sufficient to find the closure and a numerical root-finding approach is necessary.Following Foucart et al. (2015); Weih et al. (2020a), in FIL-M1 this is done by looking for the root of via a Newton-Raphson method in conjunction with a Brent method as a fallback in case the first approach fails.It is worth pointing out here that the choice of a closure function in the form of a relation  = () ∈ [1/3, 1] is not unique and indeed it represents one of the main sources of uncertainty in the M1 scheme and more in general in any truncated formalism approach (Rezzolla & Miller 1994).The choice of the Minerbo closure, which is the classical form of the maximum-entropy closure, is made following Murchikova et al. (2017), where it was identified as the closure which performed best on average in a series of tests among a series of possible analytic closures.FIL-M1 , however, leaves the option of implementing additional closures and comes equipped with the the Levermore closure (Levermore 1984).

NUMERICAL SCHEME
In this section we provide a detailed description of the numerical algorithm employed by FIL-M1.As a first step, it is worth remarking that the full set of the GRMHD and radiative-transfer equations represents a system of nonlinear partial differential equations that can be cast in a conservative form (Rezzolla & Zanotti 2013) where  is the state vector and contains all of the conserved quantities that are evolved in time, while   and  are the fluxes and source terms, respectively.For compactness, we will concentrate here only on the subset of variables that are relevant for the radiative-transfer portion of the system and remind the interested reader to the various papers where the corresponding GRMHD part is presented (see, e.g., Fambri et al. 2018).In this case, the system of equations can still be cast in the conservative form ( 51), but the state vector will be just a part of the full state vector and we will distinguish the source term  into a part that is related to the radiation quantities, and that we indicate with G, and another one that contains instead also information on the spacetime metric and extrinsic curvature, that we refer to as S. As a result, the conservative form of the radiativetransfer equations reads where and [see Eq. ( 14) for a definition of with and The importance of casting the system in the conservative formulation ( 53) is that can use standard high-resolution shock-capturing (HRSC) methods to numerically solve Eqs.(52).In particular: FIL-M1 discretizes the fluxes in Eqs. ( 52) by standard second-order accurate finite-volume techniques.Special care has to be taken when dealing with the M1 system for two reasons: (i) when the mean free path of the particles tends to zero, the equations tend asymptotically to diffusion equations for the energy and number density of neutrinos.Standard techniques for hyperbolic equations do not reproduce this limit accurately introducing instead significant amounts of numerical dissipation in the solution; (ii) in hot and dense regions as those normally occurring in BNS simulations, the collisional source terms G can become stiff leading to inaccurate solutions.
We discuss below how to address these issues in practice.

Flux discretization
As anticipated, FIL-M1 discretizes the equations via second-order accurate, HRSC finite-volume methods.Given a conserved variable    discretized at time   and at a spatial point   , we write the semidiscrete equivalent of the system (52) as where the dots indicate the collisional sources whose treatment will be detailed in sec.3.2.The numerical fluxes F  are computed at cell interfaces based on the reconstructed values of the primitive variables for the M1 scheme.In FIL-M1 we reconstruct (/, ,   /) to ensure the causality of the energy fluxes using the standard secondorder total-variation diminishing (TVD) monotonized central reconstruction scheme (Foucart et al. 2015).The code has the option of also using a simple MinMod reconstruction or a third-order accurate WENO reconstruction method (see, e.g., Rezzolla & Zanotti 2013, for details).To compute the fluxes, we utilise the two-wave Harten-Lax-van Leer-Einfeldt (HLLE) approximate Riemann solver (Harten et al. 1983;Einfeldt 1988) with eigenspeeds (Foucart et al. 2015;Shibata et al. 2011;Weih et al. 2020a) which read where and As mentioned above, when computing the numerical fluxes at cell interfaces we need to ensure that the scheme is asymptotically preserving: i.e., we need to correct the fluxes in a way that will make them suitable to solve a diffusion-type equation in the opticallythick limit.Several approaches to this problem have been proposed in the literature.One possibility, employed for instance by Foucart et al. (2015) and Weih et al. (2020a), is to interpolate between the HLLE fluxes in the optically-thin regime and a diffusive flux obtained from a finite-difference approximation of the derivative of  in the optically-thick one.This approach, while shown to yield good results, is cumbersome and requires the numerical differentiation of an expression which itself contains a derivative.On the other hand, Radice et al. (2022) correct the fluxes by using centered differences in the optically-thick regime and employing a flux limiter to hybridise the high-order centered flux with a diffusive flux near shocks or extrema of the solution.While this second approach is computationally less expensive, it neglects altogether the causal information coming from the eigenvalues of the system.For these reasons, in FIL-M1 we apply the following correction to the HLLE fluxes for energy and number densities (see also Kuroda et al. 2016;Skinner et al. 2019;Cheong et al. 2023, for a similar approach) where the diffusion-limiter factor  is defined as and represents effectively is the mean-free-path of neutrinos computed with an opacity  which is taken to be the average among the two neighbouring cells at the interface.For the energy fluxes   we instead apply the simpler correction given by Expression ( 64) effectively switches off the diffusive term in the HLLE solver at high opacities, turning the fluxes into a form similar to a centered-differencing scheme (where however the fluxes are computed at the left and right side of the interface as opposed to the cell centres of neighbouring cells).

Collisional sources
As anticipated above, the reaction rates (emissivity, opacities) involved in the computation of the collisional sources, and which serve as coupling between the radiation and the fluid, can become very large in the hot and dense matter found in BNS mergers (Endrizzi et al. 2020).As a result, the collisional sources in Eq. ( 52) are not only a nonlinear function of the evolved variables but can also become very stiff.For this reason it is common for M1 schemes (Weih et 2022) and evolve the system in time according to where the the upper index ( * ) is employed to indicate an intermediate step in the time evolution and where fluid is updated at the end of the second substep.We then need to solve the nonlinear system where we indicate by X the vector of partially updated variables (i.e., where the flux and geometric source terms have already been added).Since the reaction rates are kept fixed through the neutrino update, the source term for the neutrino number density decouples from the rest of the system and can be inverted analytically.Following Radice et al. (2022), we solve this by a globally convergent Newton-Raphson procedure where the full Jacobian of the system is evaluated.In particular, the equation for the collisional source term of the neutrino number density, i.e., Eq. ( 15), is linear in the evolved variable and its inversion is trivial.We therefore only need to focus on the subsystem that comprises the sources for  and   , and Eq. ( 14), that is, a system of four coupled nonlinear equations in the evolved radiation moments.The Jacobian can be formally written down as with Ū being defined as the sub-vector of evolved radiation variables where the number density has been excluded.We thus have Note that the derivatives can be computed from Eqs. ( 36) and ( 37), and that the term in   proportional to the hypersurface normal  does not contribute to the Jacobian, since only     enters the collisional source terms.The relevant derivatives are (see also Radice et al. 2022;Izquierdo et al. 2022) As an initial guess, we also use the first-order in  approximation to the sources, which in the fluid frame reads where  and   are evaluated after the explicit update.To transform back to the evolved variables, the closure is assumed to be optically thick (Radice et al. 2022), Eqs. ( 27) and ( 28) are used and the projection  •  can be found via the orthogonality to  and is simply given by In particular, the process of inverting the implicit sources in FIL-M1 proceeds as follows (i) The closure is updated using the Eulerian-frame moments already updated according to the explicit part of the scheme.At this point, the moments are also limited to ensure that the energy density is larger than a small but positive "atmosphere" value and the energy flux is limited to ensure that the flux factor (the ratio | |/ of the energy flux to the energy density) is not greater than unity (which would otherwise break causality); (ii) The solver is initialised according to the initial guess in Eq. ( 77); (iii) One iteration of the inversion method is performed; (iv) The closure is updated again and the moments checked for consistency with the conditions of point (i), namely that the energy density is above atmosphere and the fluxes respect causality; (v) Upon convergence, the source term of the neutrino number density is computed based on the updated neutrino moments.

Transport Rates
What remains to be fixed at this point are the transport rates dictating the coupling of the radiation field to the fluid, i.e., the absorption opacity  er  , the scattering opacity  er  , the emissivity  and the corresponding number rates.Since FIL-M1 is a grey M1 scheme, the reaction rates need to be appropriately averaged over the neutrino energies.In particular, from the energy-integrated radiation transport equation where F represents the radiation distribution function.
Table 1 offers a list of all the reactions employed in FIL-M1 and the energy-integrated rates for all these processes can be estimated analytically as a function of the local thermodynamic state of the fluid (explicit expressions will be given in Appendix A).Furthermore, the chemical potential of neutrinos as a function of (, ,   ) is calculated, as is customary for radiative-transfer codes, assuming cold beta-equilibrium

𝜇
eq where    indicates the nucleon mass difference and all the chemical potentials  * include the rest mass contribution.While this is appropriate for cold neutron-star matter, significant corrections should be applied in high temperature regions of the EOS table which in principle are probed during BNS merger simulations (see, e.g., the discussion in Hammond et al. 2021).Many leakage codes in the literature employ a corrected version of this equilibrium chemical potential which is suppressed by a factor (1 − )) where  is the optical depth.This correction ensures that the chemical potential goes to zero as matter becomes optically thin.Analogously to (Foucart et al. 2015) we find that the choice between  eq and the "leakage" prescription has very little impact on the simulations performed with M1, and we employ the corrected chemical potential by default.
Since the dominant charge current absorption rates are proportional to the square average energy of neutrinos, we can employ the improved knowledge about the neutrino spectrum coming from our effective evolution of the neutrino number density to correct the free rates.Indeed, when neutrinos decouple from the plasma and become free streaming, their temperature no longer coincides with the local fluid temperature and remains roughly constant.Because of this, computing the absorption and scattering rates due to charge current reactions according to the local thermodynamic state of the fluid leads to a large underestimation of the neutrino transport rates in optically-thin regions.For this reason, and in analogy with Foucart An anlogous correction is applied to all the scattering opacities.The neutrino temperature   is defined as where F  () is the Fermi integral of order  and is defined as and the average neutrino energy is extracted from the evolved variables according to After this temperature correction has been applied, and in order to guarantee that the neutrinos are in thermal equilibrium with the fluid in the optically trapped regions ( 1), we follow (Foucart et al. 2015;Radice et al. 2022) and recompute the emission of electronflavour neutrinos according to an approximate form of Kirchhoff's theorem given by where  is the neutrino fugacity and  is the black-body spectrum with the statistical weight of the neutrino flavour   = 1,   = 4.In particular, we set  er,nr (  , ν ) =  er,nr  * ,  . (90) Note that the values of  * and  *  appearing in Eqs. ( 89) and ( 90) are set depending on the timescale of beta equilibration.More precisely, we estimate the beta-equilibration timescale as , and compare it with the simulation time step.If  − > Δ, then we simply set ( * ,  *  ) = (,   ).On the other hand, if  − < Δ/2, then the equilibration process cannot be resolved by the simulation.Since the rates are updated only once per timestep, and remain constant throughout the nonlinear solver iterations used to determine the collisional source terms, using the local state of the fluid can cause overshooting the equilibrium state and thus lead to spurious oscillations that can be especially troublesome near black-hole horizons or neutron-star surfaces.For this reason, following Radice et al. (2022), we utilise instead of (,   ) their value at beta equilibrium.Crucially, neutrino-less beta equilibrium is not appropriate in this instance since the fluid is optically thick ( − ∼ 1/).We therefore need to perform a nonlinear root-finding iteration to find ( * ,  *  ) = ( eq ,  eq  ).
In particular, we solve the following equations (Perego et al. 2019) =   eq , where (,   ) = (1+ ) (, ,   ), with  being the specific internal energy computed from the EOS table, and ,   are defined as It is worth noting that this inversion requires interpolating the table for  at each iteration.Finally, in the intermediate the case in which 1/2 <  − /Δ < 1, we determine ( * ,  *  ) via a linear interpolation between the evolved and equilibrium values.

IMPLEMENTATION TESTS
FIL-M1 is written to be used in conjunction with the EinsteinToolkit (Loeffler et al. 2012;Zlochower & et al. 2022), exploiting the Carpet AMR driver (Schnetter et al. 2006) and the evolution code-suite developed in Frankfurt consisting of the FIL code for the higher-order finite-difference solution of the GRMHD equations and of the Antelope spacetime solver (Most et al. 2019) for the evolution of the constraint damping formulation of the Z4 formulation of the Einstein equations (Bernuzzi & Hilditch 2010;Alic et al. 2012).Most of the analysis of the data presented in the following sections is performed with the kuibit Python library (Bozzola 2021).
In what follows we will present a series of standard and new implementation tests for a two-moment radiation-transport schemes in general relativity.The tests aim at asserting the efficacy of the various aspects of the numerical scheme employed by FIL-M1.It is worth pointing out that all the tests were performed on a 3D Cartesian grid even when underlying symmetries are present in the test and that, unless otherwise stated, all the tests employ a fixed CFL factor of 0.25.

Straight Beam
The first test presented is the so-called "straight-beam" test and consists in setting up a beam of radiation oriented along a single spatial direction.We therefore set up a grid covering [0, 0.75] × [−0.2, 0.2] × [−0.1, 0.1] ∈ R 3 with 300 × 160 × 80 equidistant points.A radiation beam parallel to the −axis is injected from the left side of the domain by setting  = 1 =   in the outer boundary ghostzones of the domain with  = 0.
Since the spacetime is flat and the background fluid is perfectly transparent (i.e., we set  er =  er = 0), the beam is expected to travel straight towards the positive end of the −axis and at the speed of light.The outcome of this test can be observed in the left panel of Fig. 1, which reports the distribution of the normalised radiation energy density at  = 0.4.Note that FIL-M1 maintains a reasonably focused straight beam, which propagates at the correct speed as indicated by the expected position of the beams' front, which is marked with a vertical red line.Naturally, a small amount of diffusion is present both at the edges of the beam and at its leading front.This diffusion converges away with resolution and is the result of the interaction with the surrounding artificial atmosphere.

Curved Beam
Since FIL-M1 is meant to be employed in general-relativistic simulations, we next consider a test that involves curved spacetimes and hence non-trivial values for the metric and curvature terms on the right-hand-side of the radiative-transfer equations.For this reason, we consider another classical test where a beam of radiation is injected in a Schwarzschild spacetime to validate the capability of the code of captured correctly the bending of rays in a curved spacetime (see also Foucart et al. 2015;Radice et al. 2022;Weih et al. 2020a).We therefore set up a grid spanning [0, 5] × [−0.1, 0.1] × [0, 4] ∈ R 3 with a uniform spacing of Δ = 0.05 in all directions.We then employ a background fluid that is completely transparent and static around a unit-mass Schwarzschild black hole place at the origin of the Cartesian grid.The beam is injected in the region where  = 0 and 3.0 ≤  ≤ 3.5 with  = 1 and the energy flux set up so that     =  2 ,  , = − , /, so that the beam is initially parallel to the −axis (in this test the CFL factor is set to 0.2 as required by a stable Runge-Kutta-3 evolution scheme).
The results of this curved-beam test are reported in the right panel of Fig. 1 and show that FIL-M1 is able to reproduce the bending of radiation geodesics fairly accurately, albeit with the presence of some non-negligible dissipation.In turn, this leads to a significant broadening of the flux and consequent reduction of the intensity of the beam's section when it reaches the −axis at  = 0. (the integral of the energy density at  = 0 is reduced by 50% when compared with the initial value).Also shown in the figure with dashed red lines are the two geodesics bounding the initial beam of radiation.As can be readily observed, the vast majority of the beam's energy remains enclosed within the analytically correct trajectories.Overall, this behaviour is in line with what is seen from other similar codes in the literature when performing this test (Foucart et al. 2015;Radice et al. 2022;Weih et al. 2020a).

Shadow casting
In the next two tests we consider the interaction of a radiation beam with an obstacle whose absorption properties are extremely large so that the radiation will be fully absorbed when interacting with the obstacle.In turn, this will imply that the radiation field downstream A practical guide to a moment approach in NR 9 of the obstacle will be very small and thus the obstacle will effectively cast a shadow.Both of the tests are performed in a flat spacetime.We first consider a straight beam as the one presented in Sec.4.1 entering the Cartesian domain from the left boundary.The grid spans [−4, 4] × [−2.5, 2.5] × [−2.5, 2.5] ∈ R 3 and is covered by 640 × 400 × 400 equidistant points.The spacetime is again a flat one and the domain is optically thin except for a cylindrical region A centered at  A = (−1.5,0, 0) with a radius of  A = 0.5, where the opacity to absorption is set to a very large number, i.e.,  er  ( ∈ A) = 10 10 .The results of the test are presented in Fig. 2 and clearly show that FIL-M1 is able to produce a sharp shadow with only little diffusion behind or inside of the high-absorption region.
Next, in order to determine the extent to which the Cartesian grid affects the beam propagation and absorption onto the obstacle, consider a scenario in which the radiation field has a spherical geometry.Therefore, following Just et al. (2015) and Kuroda et al. (2016), we set up a Cartesian grid spanning [−7.5, 7.5] × [−7.5, 7.5] × [−7.5, 7.5] ∈ R 3 and covered by 196 3 equidistant points.We furthermore set up a refinement level centered at the origin whose radius we set to 2.5 and whose refinement factor is taken to be two.This allows us to check for a potential interference among refinement levels with the two-moment scheme.
The geometry is again flat and the radiation energy density is set to a small but positive value everywhere of  0 = 10 −15 .The radiation energy flux is initially set to zero on the whole grid.The domain is optically thin except for a central spherical emitting region E of radius  E = 1.5 and an absorbing sphere A of radius  A = 1 centered at  A = (3.5, 0).We then set the scattering opacity to be zero across the whole domain, while the absorption opacity is set to follow the prescription The emissivity is set to have an equilibrium energy density  eq = 10 −1 for  ∈ E and to understand the rationale behind this choice we recall that sources for  er  = 0 and   = 0 are [see Eqs. ( 14) and ( 54)] from which it is apparent that for zero flux  eq =  er / er  = const..As can be seen from Fig. 3, FIL-M1 handles this test-case almost as well as its counterpart where the radiation flux was aligned with one of the grid principal directions, without excessive diffusion near the absorbing region and without artefacts spoiling the spherical symmetry of the simulation due to the Cartesian grid or the presence of a refinement level (only small oscillations can be seen at  3.5).It is worth pointing out that the flux factor is limited to be smaller than 0.999 in FIL-M1.x t = 4.0 .Note that also in this case only a small diffusion is present downstream of the sphere despite the propagation is not along one of the main coordinates of the numerical grid.

Diffusion in a stationary background medium
The next test probes the effectiveness of FIL-M1 in capturing the asymptotic diffusion limit of the M1 equations.We follow the standard approach in this type of test (see, e.g., Pons et al. 2000;Kuroda et al. 2016;Weih et al. 2020a;Radice et al. 2022) and track the evolution of a spherical radiation wave whose profile is initially a Dirac −function at the origin of the Cartesian grid.Hence, we set up a stationary medium in flat spacetime with  er  =  er = 0 and  er  = (10 2 , 10 5 ).Within this scenario, the M1 system reduces to a diffusion equation for the energy density whose analytic solution is given by (Pons et al. 2000) where  is the number of dimensions.We therefore layout a 3D Cartesian domain [−0.5, 0.5] × [−0.5, 0.5] × [−0.5, 0.5] ∈ R 3 covered by 100 3 points in the case of  er  = 10 2 and 32 3 points with two nested refinement levels at  = ±0.25 and  = ±0.15 in the case of  er  = 10 5 .As a result, the diffusion-limiter factor  appearing in Eq. ( 64) is simply given by  = 1 [see Eq. ( 62)], so that the equations are being solved with the full HLLE solver.In the case of  er  = 10 5 , on the other hand, the Peclet number, defined as the opacity normalised to the grid spacing, i.e., Pe :=  er  Δ, is Pe 10 3 at the center of the grid and changes with the refinement level.As a result, the equations are solved with the diffusion part almost completely switched off in the Riemann solver.
Figure 4 shows the results obtained by FIL-M1 in this test and reports the very good agreement with the analytic solution, thus validating the ability to properly capture the asymptotically diffusive limit of the M1 equations.

Diffusion in a moving background medium
Despite the fact that in some of the tests considered so far (in particular the ones shown in Sec.4.3 and 4.4) the code was confronted with stiff collisional source terms, none of these has really put the implicit source-term solver to the test.Indeed, whenever the background is set to be static (  = 0), the fluid frame moments obviously coincide with the Eulerian ones, and the source terms become linear functions of the evolved variables.In particular, the initial guess used by the solver Eq. ( 77) are already the solution.We therefore set out to perform the test first introduced by Radice et al. (2022), where a diffusing radiation wave is evolved on a moving background medium.In particular, we set the scattering opacity  er  = 10 3 , a background velocity   = 0.5 and initialise the radiation energy density in the Eulerian frame to simply be We then use Eq. ( 28) with   = 0 to initialise the flux assuming a completely optically-thick medium.As fiducial solution to compare to, we employ Eq. ( 95) translated at a velocity   = 0.5, with  = 1 due to the slab geometry of the problem and with a diffusion coefficient appropriately scaled to get Eq.( 97) as initial condition.
For this test we employ a grid covering  = [−3, 3] in the −axis and spanning only a few points in the other two directions.The grid is covered by 600 points on the -axis, leading to a spacing of Δ = 0.01 code units.The numerical solution obtained by FIL-M1 is shown in Fig. 5, and reproduces the fiducial solution within the expected error level.The importance of this test is that it probes the capability of the code of dealing with stiff sources in a moving medium where the nonlinear solver is strictly required in order to converge to the true solution.Moreover, it provides an additional validation of the flux formulation in the diffusive limit.

Uniformly radiating-sphere test
The final test we consider in a flat spacetime is represented by the so-called uniformly radiating-sphere test and consists in a spherical region E of radius  E centered around the origin of the coordinate system and that is allowed to emit and absorb radiation.The sphere is surrounded by vacuum and the opacity and emissivity are therefore defined as While this test can be thought of as a very crude model for an astrophysical spherical source of radiation, its important lays in serving as a test of the collisional source terms when the absorption coefficient is nonzero.This problem admits an analytic solution where the distribution function is given by (Pons et al. 2000;Murchikova et al. 2017) with For this test we extend the domain to 4  E in all directions with a resolution of Δ = 0.0125 and employ reflection symmetries across all coordinate planes to reduce computational costs.
The results obtained by FIL-M1 are shown in Fig. 6 using three different closures: the Minerbo closure (Minerbo 1978), which is the default for FIL-M1, the Levermore closure (Levermore 1984), which consists in always assuming an optically-thick regime in the fluid frame, and with the fit to the analytic closure proposed by Murchikova et al. (2017).As can be seen, the Minerbo closure reproduces the result quite well although it does not converge to the analytic solution.On the other hand, the fit to the analytic closure of Murchikova et al. (2017) is accurate to within 1% for this specific problem and FIL-M1 can be seen to accurately reproduce the analytic result with this choice.These results emphasise how the choice of a closure that is most appropriate for the problem at hand can considerably improve the accuracy of the M1 approximation.It is therefore desirable when developing an M1 code to allow by construction for the use of arbitrary closures, as done in FIL-M1 .

NEUTRON-STAR SIMULATIONS
The next series of tests is rather different as it meant to assess the ability and self-convergence properties of FIL-M1 in providing accurate solutions under the physical conditions that are expected to be present in realistic simulations involving neutron stars and where neutrino transport plays an important role.For this purpose, we consider three different tests.First, and following previous literature (Foucart et al. 2015;Just et al. 2015;Foucart et al. 2020), we compare the results obtained by FIL-M1 with those provided by the open-source energy-dependent M1 code GR1D corresponding to a proto neutron star resulting from a core-collapse supernova (Sec.5.1).This test will show that the results obtained by our code in a realistic context and when including all source terms and coupling to the astrophysical plasma are consistent with well-tested implementations of similar schemes.Secondly, we probe the self-convergence of FIL-M1 in the context of a simple oscillating nonrotating star (Sec.5.2).Finally, we discuss in detail the simulation of the head-on collision of two neutron stars (Sec.5.3), whose results will be made publicly available so as to providing a new multidimensional and non-trivial test for future M1 implementations.

Core-collapse supernova test
For this test we simulate the collapse of a 20 progenitor (Woosley & Heger 2006) employing the tabulated DD2 EOS (Hempel & Schaffner-Bielich 2010) with the publicly available GR1D code (O'Connor 2015).In particular, we import the (1D) radial solution obtained by GR1D after core bounce into the 3D Cartesian grid employed my FIL-M1 with reflection symmetries across all coordinate planes.Furthermore, we add five static mesh refinement boxes at radii RL = {40, 80, 120, 160}  , so that the outermost box has an outer boundary at  max = 208  is covered by 26 × 26 × 26 points leading to a finest spacing of Δ 5 = 0.5  740 m.On the other hand, the (radial) grid in GR1D covers the domain [0, 13]  with a uniform spacing Δ  = 0.2  and extends outwards with a logarithmically spaced grid.The import is performed with a simple linear interpolation.The weak rates in GR1D are chosen to be as similar as possible as those in FIL-M1.In particular, we enable the re-computation of the emission coefficients from the absorption opacities using of Kirchhoff's law and we deactivate in FIL-M1 the calculation plasmon-decay emissivity since this is not supported by the the neutrino-interaction library employed by GR1D, namely, NuLib (Sullivan et al. 2015).
The systems of radiative-transfer equations are then evolved in both GR1D and FIL-M1 from  ini.−  bounce 0.15 s, and up to  fin.8 ms, keeping the spacetime fixed in both codes and updating only the hydrodynamical variables in FIL-M1 through the backreaction of neutrinos.This is done for two different reasons.First, because the 3D Cartesian grid cannot cover the whole domain that would allow a self-consistent hydrodynamical evolution without significant computational overhead.Second, because we are only interested in comparing the radiation-transfer implementations across the two codes.
The comparison of the results obtained by the two codes is shown in Fig. 7, with the left panel referring to the temperature and the right one to the electron fraction.Clearly, the agreement in the temperatures is remarkable, while the electron fraction suffers from larger uncertainties.This is to be expected since the solution from FIL-M1 is energy integrated (grey) while that in GR1D takes into account twelve energy bins.Overall, the results reported in Fig. 7 are comparable with those shown in similar code comparisons (Foucart et al. 2015).The overall consistency between the two radiative-transfer solutions is confirmed also by the panels in the top row of Fig. 8, which show the profile of the neutrinos average energies for the three neutrino species.The profiles show very good agreements in the case of electron-flavour neutrinos, while a systematic overestimation characterises the effective average energy of heavy neutrinos and was already reported by Foucart et al. (2016) (the case that can be directly compared to FIL-M1 is the one reported as  = ∞ in Tab.(1) of the reference).It is worth noting that the agreement shown in the average energies is indeed remarkable since all information about the neutrino energy spectrum is lost in a grey M1 scheme.We believe that the reason why FIL-M1 is able to reproduce the results from the energydependent transport code is that it evolves an additional equation for the number density of neutrinos, namely Eq. ( 11).In particular, the correction to absorption and scattering rates detailed in Eq. ( 82) was already shown in Foucart et al. (2015) to be crucial in allowing a grey M1 scheme to capture the features of the temperature profile in this test, and the average energy of neutrinos contained in this correction can be estimated more accurately through the additional information provided by the evolved neutrino number density.Finally, we report in the bottom row of Fig. 8 the neutrino luminosities defined as where  * stands for either   , ν , or   , Σ is a coordinate 2-sphere placed of radius  and Σ  its oriented unit normal.Also in this case, the results presented in lower panels of Fig. 8 are comparable with those reported in the literature and validate the ability of FIL-M1 to describe a non-trivial and dynamical radiative-transfer problem.

Self-convergence of a hot nonrotating star
We next consider a very important test that has never been performed in previous implementations of an M1 scheme, namely establishing the self-convergence properties of FIL-M1 in a fully coupled system in which both the Einstein and the GRMHD equations are evolved.
While this test should be straightforward in principle, it is extremely challenging in practice.First, applying such a test to a full binary system would be computationally prohibitive as the costs associated with a series of simulations, two of which need to be at resolutions sufficiently high to show a convergent postmerger behaviour, would easily exhaust the computational resources typically available in research groups.Second, these simulations would require sufficiently long times after the merger to allow the codes to recover from the low convergence order produced by the large shocks at merger (see  The curves refer to the last available snapshot at  −  ini.8 ms, while the dashed line refers to the initial solution at  ini.−  bounce 150 ms.Note that agreement in the temperature is within the expected error levels and, more importantly, the "gain region" , where neutrino absorption exceed emission and net heating occurs, starts at  100km and is properly described by the grey scheme of FIL-M1. .Top row: Average neutrino-energy profiles for the three neutrino species in the evolution of a core-collapse supernova remnant as computed by FIL-M1 (red solid line) and by the GR1D code (red solid line); the curves refer to the last available snapshot at  fin.8 ms.Note that electron-flavour neutrinos show an excellent agreement across the whole domain, whereas heavy-lepton neutrinos suffer from a systematic error in the optically-thin region.Bottom row: The same as above but for the neutrino luminosities profiles.
this, it is not surprising that rigorous self-convergence tests for simulations of binary mergers have not been performed before in previous implementations of an M1 scheme (however, see Radice et al. 2022, for self-convergence tests on simpler setups).
In view of these considerations we validate the convergence properties of a simpler scenario, namely, the neutrino luminosity from a spherically symmetric hot neutron star described by a temperature and composition-dependent EOS.It should be stressed that while simpler, this test remains challenging for at least three different reasons.First, being hot, the star does not have a background equilibrium which it can converge to as the resolution is increased (Font et al. 2002); as a result, although in a convergent regime, the star will exhibit a secular behaviour as it cools and expands.Second, the inevitable oscillations that are triggered in the star at the initial time will lead to periodic weak processes and neutrino emissions that will introduce short-timescale oscillations in the overall behaviour.Finally, the quantity whose convergence properties we want to show, namely, the neutrino luminosity, is a derived quantity and not a conserved one.As such, it will be affected by a number of other numerical operations (most importantly root-findings and interpolations) that inevitably will impact the convergence order.Notwithstanding these considerations, we will show below that FIL-M1 is able to provide a converged solution for this scenario.
For our test we consider a nonrotating stellar model constructed from the isentropic slice with entropy per baryon  = 1  B of the DD2 EOS (Hempel & Schaffner-Bielich 2010) and with a central density of   = 7.41 × 10 14 g/cm 3 ; the star has an ADM mass of TOV = 1.98  and a central temperature of the order of ∼ 30 MeV.The star is then evolved with three different resolutions differing from each other by factors of two and respectively of ΔLR 531.6 m (low resolution), ΔMR 265.8 m (medium resolution), and with ΔHR 132.9 m (high resolution).In all cases, we consider a fixed mesh four refinement levels, where the coarsest grid extends to ∼ 380 km in each direction.The neutrino luminosities are extracted at a coordinate sphere of radius  150 km, which is large enough to guarantee that no spurious diffused material will ever reach the detector.
The luminosities of each neutrino species for each resolution is shown in Fig. 9, where it is already clear that the code displays a very consistent behaviour across resolutions for all of the three neutrino species, although the strong dependence of the location of the neutrinosphere induces changes of up to a factor five between the low and high-resolution simulations (see discussion below).We then compute the convergence order for a given observable  computed at resolutions ℎ, ℓ,  as (Rezzolla & Zanotti 2013) where ℎ < ℓ <  and ℎ = ℓ/, ℓ = /.In using expression (100) it is of course assumed that the truncation error on the observable  can be expressed as a simple power-law with index , which is only true for grid spacings that are small enough to make the sub-leading terms in the error expansion negligible (i.e., in the so-called "convergence regime").The convergence order is reported in Fig. 10 for each neutrino species, showing with dashed lines the instantaneous convergence order and with solid lines of the same colour the timeaverage corresponding order employing a uniform moving filter of 0.2 ms.It is worth remarking that we compute all convergence orders using Eq. ( 100) and that in the case of neutrino luminosities  is intended as the timeseries of the observable, whereas for the lapse and rest-mass density  indicates the pointwise value of the field evaluated at the center of the grid.
The oscillating behaviour of the unfiltered convergence order is most likely due to the presence of additional, high-frequency oscillation modes that are captured in the star evolved with the highest resolution.Since the short timescale variability of the neutrino luminosity is mostly dictated by the expansion/contraction of the neutrinosphere, that is defined as the location where the optical depth   * = 1 and which itself is regulated by the stellar oscillations (Galeazzi et al. 2013;Perego et al. 2014), the luminosities from simulations with higher resolutions inevitably contain higher frequency components, which add up when the errors are computed.Indeed, an analogous oscillatory behaviour can also be observed when looking at the locations of the neutrinospheres as reported in Fig. 11.In particular, the figure offers a spacetime diagram of the temperature evolution for the simulation at the highest simulation.Also reported with coloured solid lines are the worldlines of the neutrinospheres in the three species, while the white dashed line is the worldline of the putative stellar surface, which we set to be where the rest-mass density reaches a specific value, i.e.,  = 10 12 g cm −3 .
The challenging nature of this self-convergence test is clear also when considering the convergence order computed on the pointwise values of the central rest-mass density and lapse function reported in Fig. 10 with coloured solid lines.Also for these quantities, that are genuinely part of the vectors of numerically evolved variables, a degradation of the expected convergence rate, i.e., an order of three (Most et al. 2019), is measured.This is not particularly surprising and is due to the presence of the stellar surface, where small shocks inevitably appear and lead to a degradation of the convergence order (see Radice & Rezzolla 2012;Radice et al. 2014a, for a discussion).On the other hand, the convergence order for the neutrino luminosities reported in Fig. 10 is of the order of 1 − 1.5, which is smaller than expected second order implemented in FIL-M1 , but aligned with what already observed for the bulk spacetime and hydrodynamical variables.We also note that the convergence order is somewhat different for the various neutrino species.Since the radiative-transfer equations solved are independent of the neutrino species, the reason for different convergence behaviours should be sought in the different hydrodynamic and thermodynamic conditions experienced by the different neutrino species.More specifically, the different convergence behaviours follow from the various depths within the star where the various neutrinos decouple from the fluid and become free streaming.As can be clearly seen in Fig. 11 through the worldlines of the neutrinospheres, heavy lepton neutrinos systematically decouple deeper within the star than other neutrino species, whereas electron neutrinos are strongly coupled with matter essentially up to the stellar surface.Because the position of the latter is in general in a lowerorder convergent regime (large gradients and small shocks develop near the surface), it is natural to expect that also electron neutrino luminosities are those with the lowest convergence order.Moreover, the oscillatory behaviour of the convergence order in Fig. 10 is not the directly related to the eigenfrequencies of the underlining star, but rather produced by the slightly different oscillation periods at different resolutions.Overall, the results presented in section provide evidence that the solution of the full set of the Einstein, GRMHD, and radiative-transfer equations leads to an overall convergence order which 3 for the hydrodynamical and spacetime variables, and 2 for quantities related to neutrinos.The thick lines are calculated as a moving average with a window of width of 0.2 ms, whereas the unfiltered data is shown with the dashed lines of the same colour.Also reported are the convergence orders on the pointwise values of the central rest-mass density (black solid line) and lapse function (orange solid line).The lower convergence order of the luminosities is due to the intrinsic lower order of the numerical scheme and the strong dependence of the neutrinosphere on the stellar temperature near the stellar surface.

Head-on collision of two neutron stars
The final test of the robustness of our code is also the most involved but hopefully also the most useful as it can be used in the future as an effective, comparatively inexpensive but all-round benchmark test for the implementation of an M1 scheme under realistic conditions of spacetime curvature but also of hydrodynamical and thermodynamical states.In particular, we consider the head-on collision of two equal-mass neutron stars described by the temperature dependent DD2 EOS (Hempel & Schaffner-Bielich 2010), each having a mass of 0.91  and a central density of   = 4.6 × 10 14 g/cm 3 .To maintain the setup as simple as possible, we do not seek a constraint- Figure 11.Spacetime diagram of the temperature evolution for the hot nonrotating star in the simulation with the highest resolution.Also reported with coloured solid lines are the worldlines of the neutrinospheres in the three species considered here, while a white dashed line marks the worldline of the putative stellar surface, which is set to be where  = 10 12 g cm −3 .satisfying initial solution but simply consider the linear combination of the spacetimes corresponding to two stars in isolation and with a boost of | in | = ±0.02along the direction of the collision.Given this setup, the resulting object is not expected to produce a black hole but a remnant star (Rezzolla & Takami 2013;Koeppel et al. 2019).The grid is set to have a total extent of  max = 256  containing five nested fixed refinement levels at positions RL = {32, 64, 80, 100, 160}  .The outermost box has an outer boundary at and is covered by 40 × 40 × 20 points leading to a finest spacing of Δ 5 = 310.08 m.
Reflection symmetry across the −plane is employed to save computational resources.
Figure 12 shows three representative snapshots of the evolution of the colliding stars, reporting with colourcodes the temperature distribution (left part of each panel) and the rest-mass density distribution (right part of each panel).In the leftmost panel the stars are about to collide and the temperature is very low everywhere except near the stellar surfaces, where spurious heating takes place due to interactions between the stars and the artificial atmosphere.The central panel refers instead to an instant during the "dynamical" phase where, the remnant star is undergoing violent oscillations and internal shock heating.Finally, the rightmost panel shows the state of the system at the end of the simulation, consisting of a collision remnant that is secularly stable and that reaches a roughly spherically symmetric state with a ring-like structure in the temperature, as also encountered in the aftermath of BNS mergers (see, e.g., Kastaun & Galeazzi 2015;Hanauske et al. 2017).The origin of this temperature distribution has to be found in the different manner in which the thermal energy is produced and redistributed within the remnant.While immediately after the violent collision of the two stars (central panel of Fig. 12) the temperature is redistributed via strong shocks and hot material reaches the deepest regions of the remnant's core, during the subsequent oscillations, the temperature is redistributed mostly via conduction and the transport of thermal energy of takes place along isopycnic levels and not across them.As a result, the thermal energy will not be able to reach the dense inner regions of the remnant, but will instead accumulate on an almost spherical ring.It is also worth remarking that unlike a BNS merger, where the stellar material is heated by the shearing of the two grazing stars at merger, the heating here comes from the very strong shocks produced by the collinear collision.As as a result, the temperature in the remnant's core is almost a factor of ten larger than what is typically be produced in a BNS merger remnant (see e.g., Tootle et al. 2022).
Figure 13 shows instead the evolution of the neutrino luminosities for the three species (left panel), of the maximum temperature (middle panel) and of the rest-mass density (right panel).Note that the variations in the neutrino luminosities are not quite simultaneous as they are produced by different neutrinospheres but are strongly correlated with the corresponding variations in temperature and restmass density which, in turn, are due to the violent oscillations of the remnant system as it attains a new equilibrium.The strongest peaks in the temperature are reached at  = 1.3 ms, when the two stellar cores bounce after the initial collision at  = 0.5 ms, after which the oscillations in the temperature are much smaller and the system enters a slow, neutrino driven, cooling stage.Interestingly, we record a slight delay between the oscillations in the rest-mass density and those in the neutrino luminosity, which are likely due to the time needed by the neutrinos to stream out of the remnant.After about 3 ms most of the oscillations have been damped and the evolution enters a quasi-stationary phase in which the the neutrino luminosities reach almost constant values and the temperature experiences very slow exponential decay due to cooling of the remnant star.Finally, note in this new equilibrium the central rest-mass density of the remnant stabilises around a value that is larger than the initial one (horizontal dotted line) but also slightly smaller than the one of a zero-temperature remnant with the same mass (horizontal dashed line).This is due to the additional internal energy gained through the collision that reduces the central density; as the remnant cools down, the zero-temperature solution is progressively reached.The most important data corresponding to this simulation, in terms of spacetime, hydrodynamical, and thermodynamical quantities, will be made accessible freely online for comparison with future M1-scheme implementations.

CONCLUSIONS
We have presented FIL-M1 , a new implementation of a two-moment radiative-transfer scheme for neutrino transport in numericalrelativity simulations of neutron stars.FIL-M1 incorporates many recent algorithmic developments in radiation transport that contribute to its stability and accuracy.In what follows we list briefly those steps in the implementations that should be followed and the pitfalls that, on the contrary, should be avoided.
• The collisional source terms of the M1 system become extremely stiff when the fluid velocity is very large.In these cases, resorting to simpler linearized expressions of the source terms is not advisable as it invariably leads to inaccurate or unstable evolutions of the system.On the other hand, the inversion of the full set of nonlinear equations (Eqs ( 67)) together with implicit-explicit time-stepping, albeit more involved, has shown to lead to stable and accurate solutions.
• Some radiative-transfer codes iterate the full state vector of the MHD and radiation variables to find the correct value for the implicit sources (Melon Fuksman & Mignone 2019).While mathematically correct, this approach is computationally unfeasible in simulations of neutron stars, where the conservative-to-primitive inversion procedure contributes to most of the computational costs.To circumvent this problem, the fluid state can be kept constant during the iterations of the nonlinear root-finding solver for the reaction source terms (which implies that the reaction rates are also kept fixed during said iteration).This approach can lead to oscillations near beta equilibrium that can amplify and lead to unstable evolutions near the stellar surface or the black-hole's apparent horizon.In these cases, fixing the reaction rates with the beta-equilibrated values of  and   [Eq.( 91)] has proven to be essential to obtain stable evolutions.
• When computing the numerical fluxes at cell interfaces it is important that the scheme is asymptotically preserving: i.e., that the optically-thick limit the fluxes are those of a diffusion-type equation.The formulation employed by FIL-M1 [see Eqs. ( 64) and ( 62)] is asymptotically preserving and has two desirable features: -it reduces to the three-state HLLE Riemann solver with the correct eigenspeeds of the M1 system in the low Peclet-number limit, so that the causal structure of the underlying partial differential equation is properly captured by the scheme.
-it circumvents the need of explicitly replacing the numerical flux with the correct asymptotic form, thus avoiding potential problems with the energy density in the fluid frame, as well as with the velocity-dependent terms.
Following these strategies, we have reported the results of a number of tests for the validation of the solution of the radiative-transfer equations, all of which FIL-M1 passes successfully.In addition to these standard tests in flat spacetimes, we also consider the validation of the code across non-trivial scenarios of curved spacetimes such as that involving the comparison with energy dependent but spherically symmetric radiative-transfer calculations of core-collapse supernovae.Finally, we introduced two novel tests which involve the solution of the full set of the Einstein, GRMHD and radiative-transfer equations.These are the study of the neutrino emission from an isolated hot neutron star and from the head-on collision of two equalmass neutron stars.While in the first test we show that FIL-M1 is able Three strong peaks appear and are obviously correlated but not simultaneous and are produced by the initial violent oscillations of the remnant.As the are damped, the luminosities reach a very slow exponential decay due to the cooling of the remnant star.Middle panel: Evolution of the maximum temperature, with the strongest peak occurring when the two disrupted cores bounce after the initial collision.Subsequently, the oscillations in the temperature are damped and the system enters a slow, neutrino-driven cooling stage.Right panel: Evolution of the maximum rest-mass density, which clearly shows the behaviour typical of a damped oscillator.Note that the final central density is larger than the initial one (horizontal dotted line) but also slightly smaller than the one of a zero-temperature remnant with the same mass (horizontal dashed line), to which it progressively converges.
to provide solutions in a convergent regime, the second test is particularly useful since the most important data in terms of spacetime, hydrodynamical, and thermodynamical quantities, is freely available and can be used as a benchmark test in new implementations of the M1 schemes.
The production rates for electron-positron pair annihilation relative to electron-flavour neutrinos are given by where  * = 1/137.036is the fine-structure constant,  = 5.565 × 10 −2 √︃ 1/3( 2 + 3 2  ), and the mean energy of the produced neutrinos in plasmon decay we take to be Finally, we include production of heavy lepton neutrinos via nucleon-nucleon bremsstrahlung as (Burrows et al. 2006)  =1   ,brems = 2.08 × 10 2   2  2  +  2  + 28 3      5.5 , (A21) where we set  = 0.5 (Burrows et al. 2006).It is important to note that the reported rate is for all four neutrino species combined.We then follow Ardevol-Pulpillo et al.This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .Figure 2 .
Figure 1.Left: Distribution at  = 0.4 of the radiation energy density normalised to the maximum value for a radiation beam propagating in a Minkowski spacetime (straight-beam test).The expected position of the beam front is shown as a vertical red line.Right: The same as in the left but for a radiation beam propagating near a Schwarzschild black hole.Red dashed lines indicate the exact geodesics at the edges of the injected beam.

Figure 4 .
Figure 4. Radial profiles of the energy density (top panels) and of the energy flux density (bottom panels) at different times of sphere of radiation diffusing in a Minkowski spacetime (diffusion-wave test).The left column represents the solution for  er  = 10 2 , whereas the right column for  er  = 10 5 .The numerical (analytical) solutions are shown with crosses (solid lines) and show a very good agreement between the two.Also shown with a dashed gray line in the top-right panel is the Plecet number, which is clearly always much larger than unity.

Figure 5 .
Figure 5. Radial profiles at different retarded times  −  0 of the energy density in the frame comoving with the fluid  (moving-fluid diffusion test).The numerical (analytical) solutions are shown with crosses (solid lines) and show a very good agreement between the two.

Figure 6 .
Figure 6.Radial profiles of the energy density  for a uniform sphere of radiation (uniform-sphere test).The results for different closures (Murchikova+, Levermore and Minerbo) are shown with symbols, while the exact solution is shown with a black solid line.

Figure 7 .
Figure 7. Radial profiles of the temperature (left panel) and of the electron fraction (right panel) in the evolution of a core-collapse supernova remnant as computed by FIL-M1 (black solid line) and by the GR1D code (red solid line).The curves refer to the last available snapshot at  −  ini.8 ms, while the dashed line refers to the initial solution at  ini.−  bounce 150 ms.Note that agreement in the temperature is within the expected error levels and, more importantly, the "gain region" , where neutrino absorption exceed emission and net heating occurs, starts at  100km and is properly described by the grey scheme of FIL-M1.
Figure8.Top row: Average neutrino-energy profiles for the three neutrino species in the evolution of a core-collapse supernova remnant as computed by FIL-M1 (red solid line) and by the GR1D code (red solid line); the curves refer to the last available snapshot at  fin.8 ms.Note that electron-flavour neutrinos show an excellent agreement across the whole domain, whereas heavy-lepton neutrinos suffer from a systematic error in the optically-thin region.Bottom row: The same as above but for the neutrino luminosities profiles.

Figure 9 .Figure 10 .
Figure 9. Neutrino luminosities of the three neutrino species considered and produce by a hot and radiating nonrotating star.The luminosities are shown as a function of the retarded time at the detector and for the three different resolutions (high, medium and low) employed in the self-convergence test.Note that the strong dependence of the location of the neutrinosphere induces changes of up to a factor five between the low and high-resolution simulations.

Figure 12 .]Figure 13 .
Figure 12.Representative snapshots of the head-on collision of two equal-mass neutron stars.All panels show the temperature (left half) and the rest-mass (right side).Note the appearance of a ring-like structure in the temperature of the relaxed post-collision remnant.