# Monte-Carlo closure for moment-based transport schemes in general relativistic radiation hydrodynamics simulations

## Abstract

General relativistic radiation hydrodynamics simulations are necessary to accurately model a number of astrophysical systems involving black holes and neutron stars. Photon transport plays a crucial role in radiatively dominated accretion disks, while neutrino transport is critical to core-collapse supernovae and to the modeling of electromagnetic transients and nucleosynthesis in neutron star mergers. However, evolving the full Boltzmann equations of radiative transport is extremely expensive. Here, we describe the implementation in the general relativistic SpEC code of a cheaper radiation hydrodynamics method which theoretically converges to a solution of Boltzmann’s equation in the limit of infinite numerical resources. The algorithm is based on a gray two-moment scheme, in which we evolve the energy density and momentum density of the radiation. Two-moment schemes require a closure which fills in missing information about the energy spectrum and higher-order moments of the radiation. Instead of the approximate analytical closure currently used in core-collapse and merger simulations, we complement the two-moment scheme with a low-accuracy Monte-Carlo evolution. The Monte-Carlo results can provide any or all of the missing information in the evolution of the moments, as desired by the user. As a first test of our methods, we study a set of idealized problems demonstrating that our algorithm performs significantly better than existing analytical closures. We also discuss the current limitations of our method, in particular open questions regarding the stability of the fully coupled scheme.

###### keywords:

keyword1 – keyword2 – keyword3^{1}

^{2}

## 1 Introduction

Neutrinos and photons play a critical role in numerical studies of astrophysical systems. For example, general relativistic photon transport is required to study radiatively dominated accretion disks, while neutrino-matter interactions are crucial to the explosion mechanism of core-collapse supernovae. In neutron star mergers, neutrinos do not directly affect the dynamics of the merger, yet they are the main source of cooling of the accretion disks and massive neutron stars formed as a result of many mergers. Neutrino-matter interactions also determine the evolution of the composition of the matter ejected by these mergers. These neutron-rich outflows undergo rapid neutron capture (r-process) nucleosynthesis, making mergers promising candidates as the site of production of many heavy elements (e.g. gold, platinum, uranium, Korobkin et al. 2012; Wanajo et al. 2014). Radioactive decay of the ashes of the r-process also powers bright electromagnetic counterparts to the gravitational waves emitted by mergers: the optical/infrared transients called kilonovae (Lattimer & Schramm, 1976; Li & Paczynski, 1998; Roberts et al., 2011; Kasen et al., 2013). A kilonova was for example observed in the afterglow of the first neutron star merger detected through gravitational waves (Abbott et al., 2017; Cowperthwaite et al., 2017).

The mass and composition of merger outflows are the main determinant of the color, duration, and luminosity of kilonovae (Barnes & Kasen, 2013), while the composition of the outflows largely set the relative yields of different heavy elements as a result of r-process nucleosynthesis (Korobkin et al., 2012; Wanajo et al., 2014; Lippuner & Roberts, 2015). Accordingly, neutrino transport is an important component in any effort to model kilonovae and r-process nucleosynthesis in mergers. Neutrino-antineutrino annihilations can also deposit a significant amount of energy in low-density regions above the remnant (Perego et al., 2017). While, on their own, annihilations probably do not deposit enough energy to power short gamma-ray bursts (SGRBs), they may help by clearing the polar regions of baryonic matter (Fujibayashi et al., 2017).

The main objective of a radiation transport scheme is to evolve the distribution function of neutrinos or photons, where are the spacetime coordinates and are the components of the 4-momentum of the neutrinos/photons, with some affine parameter. Here we neglect neutrino masses, and the 4-momentum is thus a null vector, i.e . The distribution function of photons and of each species of neutrinos evolves according to Boltzmann’s equation

(1) |

where are the Christoffel symbols, and the right-hand side includes all collisional
processes. Solving Boltzmann’s equation thus requires the evolution
in time of a 6-dimensional function, a very steep computational challenge. The main objective of this work is to provide a
relatively cheap algorithm for general relativistic radiation transport which, while approximate for the numerical resolutions
that we can currently afford, asymptotes to Boltzmann’s equation in the limit of infinite computational resources. While our
algorithm can theoretically be used in various problems involving general relativistic photon or neutrino transport, we implement and test it with the problem of neutrino transport in merger simulations in mind. We implement the algorithm in the SpEC code^{3}

In the context of neutron star-neutron star (NSNS) and neutron star-black hole (NSBH) mergers, neutrinos have so far been modeled using approximate transport schemes, with sophistication ranging from order-of-magnitude accurate leakage schemes (Deaton et al., 2013; Neilsen et al., 2014), to various moment schemes in which only the lowest or two-lowest moments of in momentum space are evolved (Wanajo et al., 2014; Foucart et al., 2015; Radice et al., 2016). Neutrinos are also generally evolved in the gray approximation, in which information about the energy spectrum of the neutrinos is either unavailable or, in the case of our most recent moment scheme, limited to the knowledge of the average energy of the neutrinos (Foucart et al., 2016b). Spectral (energy-dependent) leakage and moment schemes have been used in Newtonian simulations of post-merger accretion disks (Just et al., 2015; Perego et al., 2016), or in the context of core-collapse supernovae (Roberts et al., 2016), but not for NSNS/NSBH mergers. This is particularly problematic because neutrino cross-sections have a strong dependence in the energy of the neutrinos. In NSNS mergers, this leads to significant uncertainties in the computation of the composition of the ejecta (Foucart et al., 2016a; Foucart et al., 2016b).

Moment schemes also require information about moments of that are not evolved in the simulation. For example, in the two-moment formalism, the energy and flux density of the neutrinos are evolved, and the evolution equations require information about the pressure tensor of the neutrinos. Semi-arbitrary choices thus have to be made in order to close the system of equations. The most common choice in recent simulations has been the Minerbo (M1) closure (Minerbo, 1978). The M1 closure, however, is only guaranteed to be correct in two regimes: in the optically thick limit, and for a single beam of free-streaming neutrinos. Otherwise, the closure is inaccurate. The best known consequence of this inaccurate closure is the presence of radiation shocks whenever neutrino beams cross or converge.

To go beyond the moment formalism, at least two directions can be considered. The first is to discretize in momentum space. The cost associated with a 6-dimensional grid is however prohibitive. Recently, a first axisymmetric simulation with a full discretization of has been performed in the context of core-collapse supernovae (Nagakura et al., 2017a). While such a scheme may one day become affordable without reducing the dimensionality of the problem, it remains at this point well beyond our reach. A second possibility is to rely on Monte-Carlo (MC) methods, randomly sampling the 6-dimensional space of . In MC codes, we evolve neutrino packets representing large numbers of neutrinos that propagate through the numerical grid. The distribution function can be reconstructed at any point from these packets, within statistical errors due to the finite number of samples/packets. Recently, the MC code Sedonu has been used to perform Newtonian evolutions of neutrinos in time-independent snapshots of supernovae simulations (Richers et al., 2017), while the general relativistic code bhlight has been developed for simulations of accretion disks (Ryan et al., 2015). While 6-dimensional general relativistic simulations of neutrino transport coupled to general relativistic hydrodynamics simulations of supernovae or neutron star mergers with a MC code appear closer than with a grid-based code, they still require a very significant investment of computational resources.

In this paper, we present the first implementation in a general relativistic hydrodynamics code of a somewhat cheaper option: a two-moment scheme where both the unknown moments of and any spectral information needed by the code are obtained from MC evolution of the neutrinos. Our method is developed in the spirit of the Variable Eddington Tensor algorithm implemented in the Athena code by Davis et al. (2012) for Newtonian simulations of accretion disks, except that we use the noisier but cheaper Monte-Carlo methods to provide the closure instead of the short-characteristic method used by Davis et al. (2012)

The fact that we have to evolve neutrinos with a MC scheme may give the impression that our algorithm does not provide any advantage over a pure MC evolution. However, using MC only as a way to close the equations for the evolution of the moments allows us to take two important, cost-saving shortcuts. First, we use MC to compute time-averaged moments of . At fixed number of neutrino packets on the grid, this allows us to reduce the statistical error in the computation of the moments without increasing the cost of the simulation. This comes, however, at the cost of smoothing over time the value of these moments. Second, we can completely side-step a common issue in MC simulations: the high cost of evolving optically thick regions, where many packets are constantly created and reabsorbed. MC codes can partially avoid this issue by switching to a different treatment of the neutrinos in optically thick regions, e.g. by using MC methods to approximate a diffusion equation through optically thick regions. We instead simply rely on the moment scheme in optically thick regions, where it is very accurate, and only evolve MC packets below a certain optical depth.

The use of a two-moment scheme offers another practical advantage: it is generally much simpler to satisfy the Hamiltonian and Momentum constraints of Einstein’s equations, as well as conservation laws, in a grid-based moment scheme than in a pure MC evolution. Finally, and most importantly, as the MC evolution asymptotes to the exact solution to Boltzmann’s equation, and as the evolution of the moments can take from that MC evolution all of the missing information about higher moments of the distribution function and about neutrino spectra that would otherwise be approximated using analytical prescriptions, the algorithm satisfies our requirement that it theoretically converges to a true solution of the transport problem as more computational resources become available.

At this point, however, we should mention an important caveat: while the test results obtained with our code, and described in this manuscript, are so far encouraging, we have not rigorously proven that the coupled MC-Moment system is numerically stable. Whether any additional work is required to guarantee that the coupled MC-moment equations are well-behaved for realistic astrophysical simulations is an important open questions, that will require additional investigation. In practice, in the near future, we plan to use our algorithm in two distinct ways. First, we will use the MC evolution to compute the neutrino spectra, the absorption and scattering opacities, and the rate of annihilation in low-density regions above a neutron star merger remnant, and to obtain detailed information about the distribution function of neutrinos in neutron star mergers. This provides a significant improvement over a gray two-moment scheme, but not the desired outcome of an algorithm converging to a true solution to Boltzmann’s equation (as, e.g., the neutrino pressure tensor would still be approximated using the M1 closure). In a second stage, we will use the fully coupled system presented in this manuscript, possibly improved to handle any additional stability issues encountered in realistic astrophysical systems. We note that our approach to this problem makes it easy to rely on either the MC evolution or analytical approximation to provide any of the missing information in the evolution of the moments of . This makes switching from a “partially coupled” MC-Moment scheme to a “fully coupled” MC-Moment scheme, as appropriate for any given project, a fairly easy task.

In the rest of this paper, we assume that . For subscripts and superscripts, greek letters are spacetime indices going from 0 to 3, and roman letters are spatial indices going from 1 to 3. When a black hole is involved in the simulation its mass is . Otherwise, we work with an arbitrary mass unit (code units), or with . A table summarizing the meaning of the various symbols used throughout this work is provided at the end of the manuscript (Table 1). We begin with a description of our numerical method (Sec. 2), then demonstrate the performance of our algorithm in a few idealized test cases (Sec. 3), and finally conclude with a discussion of the strengths and limitations of our algorithm in its current form, and of potential applications (Sec. 4).

## 2 Numerical Methods

The SpEC code can be used to evolve Einstein’s equations of general relativity coupled to the general relativistic equations of hydrodynamics. In this paper, however, we focus on the development of a new neutrino transport scheme for SpEC. For the code tests presented here, we do not evolve Einstein’s equations, while the fluid is at most evolved through its coupling to the neutrinos. From our experience developing a two-moment neutrino transport scheme in SpEC, we do not expect the full coupling to Einstein’s equations and the equations of hydrodynamics to create new problems in our evolutions, although whether this remains true for the coupled MC-moment scheme developed here still has to be tested in practice. In the following sections, we first define the reference frames in which we solve the radiation transport problem (Sec. 2.1), then discuss the implementation of the two-moment (Sec. 2.2) and Monte-Carlo (Sec. 2.3) transport schemes, the methods used to couple the two schemes (Sec. 2.4), and finally how the various pieces of our radiation transport algorithm fit within the evolution of the general relativistic radiation hydrodynamics equations in SpEC (Sec. 2.5).

### 2.1 Definitions and reference frames

We assume a 3+1 decomposition of spacetime. The 3+1 decomposition relies on a foliation of spacetime into slices of constant time coordinate , with timelike unit normal and line element

(2) | |||||

where is the 4-metric, the 3-metric on a slice of constant , the lapse, and the shift vector. Our numerical grid is discretized in the spatial coordinates . We will refer to the coordinates as the grid frame.

The fluid is described by its baryon density , temperature , electron fraction , and 4-velocity . Two special observers will play an important role in the description of our neutrino transport algorithm: inertial observers, whose timeline is tangent to , and comoving observers, whose timeline is tangent to . We also define the coordinates of the fluid rest frame , which are defined at a point so that

(3) |

with the Minkowski metric, and . We construct these local coordinates from an orthonormal tetrad , with , and . The three other components of the tetrad are obtained by applying Gramm-Schmidt’s algorithm to the three vectors (). The orthonormal tetrad and the corresponding one-forms are precomputed and stored for each grid cell, and can be used to easily perform transformations from the fluid rest frame coordinates to grid coordinates (and vice-versa) by simple matrix-vector multiplications. We note that their are many alternative methods to choose a convenient tetrad for radiation transport, e.g. Lindquist (1966); Cardall et al. (2013); Shibata et al. (2014); Nagakura et al. (2017b). Because we do not assume any symmetry when constructing our tetrad, there is no obvious geometrical meaning to this tetrad, besides the fact that its first element is tangent to the world line of an observer comoving with the fluid. In special relativistic problems and for a fluid velocity aligned with a coordinate axis, our coordinate transformation from the grid frame to the fluid rest frame is a standard Lorentz boost. For fluid velocities of arbitrary orientation, it already differs from the standard choice by a spatial rotation in the fluid rest frame.

Unless otherwise specified, vector and tensor components are expressed in the grid coordinates. The fluid rest frame is mostly used to compute neutrino-matter interactions. Quantities expressed in the fluid rest frame will be primed, e.g. for the momentum of the neutrinos in the fluid rest frame.

### 2.2 Two-moment transport

The first layer of our transport scheme is based on the two-moment formalism. A two-moment scheme was already implemented in the SpEC code, and has been used for the study of NSBH (Foucart et al., 2015) and NSNS (Foucart et al., 2016a; Foucart et al., 2016b) mergers. The general idea of the scheme, proposed by Thorne (1981) and Shibata et al. (2011), is to evolve moments of the distribution function . These can be defined from the stress-energy tensor of the neutrinos,

(4) |

Here, is the energy density of the neutrinos measured by an inertial observer. The energy flux and pressure tensor are both normal to , i.e. . Evolution equations for , , with the determinant of the 3-metric , are obtained by integrating Boltzmann’s equation:

where is the extrinsic curvature of the current spatial slice. To close this system of equations, we require two pieces of information: the pressure and the source terms . We define

(7) |

being the Eddington tensor. In this manuscript, we also write the source terms

(8) |

Here, is the absorption opacity, the scattering opacity, and the emissivity. The Eddington tensor and opacities are quantities to be provided by the MC code, while the emissivity is taken from tabulated values for neutrino-matter interactions (see Secs. 2.3.1 and 2.3.2). We note that this form for , chosen for its convenience in the tests used in this manuscript, excludes potentially important processes, e.g. pair annihilations (which couple neutrinos and anti-neutrinos) or inelastic scatterings. Energy and momentum deposition from these processes can however be computed from the MC evolution, and added as source terms to the evolution of the moments as needed. It is fairly easy to do so as long as the additional source terms are not too large, and can thus be treated explicitly

The moments and are, respectively, the energy density and flux density measured by an observer comoving with the fluid (but with vector components in grid coordinates). They can be defined from the stress-energy tensor of the neutrinos

(9) |

with the pressure tensor measured by a comoving observer, and . and can thus be computed as functions of , , and , by taking projections of .

An equivalent, and more intuitive definition of the moments can be obtained if we start from the comoving moments

(10) | |||||

(11) | |||||

(12) |

with the neutrino energy in the fluid rest frame, integrals over solid angle on a unit sphere in momentum space, and

(13) |

the 4-momentum of neutrinos, with .

We use high-order finite volume methods to evolve . A locally implicit time stepping allows us to handle stiff neutrino-matter interaction terms, while the flux terms are computed explicitly. The only terms to be treated implicitly in the equations are those containing the source terms . Thanks to the linearity of the equations in , we can solve the implicit problem exactly by inverting a 4x4 matrix for each neutrino species at each point.

In practice, we compute the fluxes and on cell faces using the HLL approximate Riemann solver (Harten et al., 1983), and a fifth-order WENO reconstruction of , , on cell faces from their cell-centered values (Liu et al., 1994; Jiang & Shu, 1996). The time discretization uses an implicit second-order Runge-Kutta method: to evolve the system by a time step , we first take a test step , with the fluxes and explicit source terms computed at the beginning of the time step. We then take a full step , with the fluxes and explicit source terms evaluated from the result of the half step.

To compute the HLL fluxes, we need the characteristic speeds of the system. For the linear system considered here, the minimum and maximum speeds across a face in the direction are

(14) |

Finally, we note that the two-moment scheme is corrected in very optically thick cells (i.e. cells in which , with the grid spacing). Without such a correction, the diffusion of neutrinos in optically thick regions is set by numerical viscosity. This correction effectively transitions between a two-moment scheme and a one-moment scheme, with in optically thick regions set by its known value in the diffusion limit. In practice, we follow the method of Jin & Levermore (1996). Details of its implementation in the SpEC code are discussed in Foucart et al. (2015).

### 2.3 MC transport

The second and most expensive layer of our transport scheme is a MC algorithm, which we implement in SpEC for this work. Our MC methods are largely inspired from the general relativistic radiation hydrodynamics code bhlight (Ryan et al., 2015). We refer the reader to that work for more details on the derivation of a general relativistic MC scheme, and focus here on a summary of the algorithm, as well as on the additional work required to use a MC algorithm to close the two-moment equations.

In a MC scheme, we evolve neutrino packets that sample the distribution function of neutrinos, . The ensemble of packets in the simulation at time , with packet representing neutrinos at the spatial coordinates and with 4-momentum , provides an estimate of the distribution function

(15) |

We note that the use of lower indices in and upper indices in is required for this equation to be a relativistic invariant. To couple the MC evolution to the moment formalism, we also rely on the MC estimate of the neutrino stress-tensor

(16) |

In a MC transport scheme, Boltzmann’s equation for can be translated into prescriptions for the creation, annihilation, scattering and propagation of the neutrino packets sampling , which we discuss in the following sections. In SpEC, the state of a neutrino packet is entirely determined by its grid coordinates , momentum , and neutrino species, as well as the number of neutrinos represented by the packet. As we neglect the masses of neutrinos, the fourth component of the momentum can easily be computed from , e.g. . The spatial components of the grid coordinates and momenta are evolved in time, while and the neutrino species remain constant during the evolution of a packet. We also regularly need the fluid rest frame energy of the neutrinos in the packet,

(17) |

where is the Lorentz factor of the fluid with respect to an inertial observer.

For high-accuracy evolution of neutrino packets in the MC framework, high-order interpolation of the fluid and metric variables to the position of a neutrino packet would be desirable. However, this would significantly increase the cost of the MC scheme. As limiting the cost of the MC algorithm is our main concern at this point, we make a simpler, less accurate choice: we use cell-centered values of the fluid variables, metric variables, and derivatives of the metric, which are already computed during the evolution of the general relativistic equations of hydrodynamics. We will see in Sec. 3 that this is unlikely to be an important contribution to the error budget of our simulations.

#### Tabulated neutrino-matter interactions

In this manuscript, we ignore inelastic scatterings and reactions involving two or more neutrinos. Neutrino-matter interactions are described by a neutrino emissivity , absorption opacity , and elastic scattering opacity , for each species of neutrinos and antineutrinos. The exact interpretation of these variables in the context of a MC algorithm will be discussed in more detail in the following sections.

We use tabulated values of these quantities produced by the open-source NuLib library (O’Connor & Ott, 2010), which provides a flexible framework to include a range of neutrino-matter interactions of importance to the merger problem, and let us choose the density of the table in the fluid variables . We use linear interpolation in to interpolate between tabulated values. We discretize the neutrino spectrum into energy bins, with bounds . NuLib provides us with values of the opacities at the center of each energy bin. We obtain the opacities at other energies by interpolating linearly in . NuLib also provides us with the emissivity per unit volume, integrated over each energy bin. The tabulated values guarantee that and satisfy Kirchoff’s Law, i.e. that for each energy bin is the energy density of neutrinos in equilibrium with the fluid, integrated over that bin. By enforcing Kirchoff’s law in that manner, we make sure that the neutrinos reach the desired equilibrium distribution in optically thick regions. For testing purposes, we also implement an alternative framework in which are provided as functions of the spacetime coordinates.

We note that incorporating inelastic scattering of neutrinos by nucleons, electrons, and nuclei into the MC evolution is not particularly complex from an algorithmic point of view, but requires very large tables for the inelastic scattering cross-sections, coupling all energy bins. Neutrino-antineutrino annihilations into pairs are more challenging if one wants to account for the blocking factor of the electrons. Ignoring that factor (which is probably acceptable in the low-density regions where pair annihilation may play an important role), one can define the energy deposition due to pair annihilation from the moments of the distribution function of each neutrino species (Fujibayashi et al., 2017), which we already have at our disposal. We plan to incorporate annihilations in that approximation, once we begin to use our code to study compact binary mergers.

#### Emission

To sample the emission of neutrinos, we compute for each cell of volume and time interval the number of neutrinos packets emitted with a fluid rest frame energy within a given energy bin . Taking advantage of the invariance of the 4-volume , and defining as the tabulated emissivity per unit volume in that energy bin, we find

(18) |

where is the desired energy of each neutrino packet in the fluid frame and . is provided as an input to the code, and effectively sets the number of packets evolved by the MC algorithm. If is the largest integer smaller than , then we create neutrinos packets, with a probability of creating one additional packet.

All packets are created with a fluid rest frame energy , and represent a number of neutrinos . The choice to initialize all neutrinos with the energy of the center of the bin was proposed by Richers et al. (2017) and has the advantage to be consistent with the way Kirchoff’s law is enforced in the NuLib tables (i.e. by balancing the emissivity integrated over the entire energy bin with the opacity at the center of the bin). The time of creation and initial position of the packets are drawn from uniform distributions in within the volume and time interval . We note that the choice of uniform sampling in the coordinates could be a significant approximation when using curvilinear coordinates, or if the metric varies significantly over the length of a cell. It may be necessary to use sampling methods which take into account variations in the proper volume element within a cell for applications requiring higher accuracy and/or for non-cartesian grid structures.

The initial direction of propagation of the neutrino packets is drawn so that the packets are isotropically distributed in the fluid rest frame. The 4-momentum of a neutrino in a given packet is, in that frame and using the orthonormal tetrad defined in Sec. 2.1,

(19) |

We draw from a uniform distribution in the range and from a uniform distribution in the range . The 4-momentum of neutrinos in grid coordinates can then be computed using that same orthonormal tetrad:

(20) | |||||

(21) |

The angles , used in this sampling process are convenient to obtain an isotropic distribution in the fluid rest frame. As opposed to angles defined with respect to the radial direction of a spherical grid, however, they have no clean interpretation. The “polar” angle is the angle between the direction of propagation of the neutrinos and the last element of the orthonormal tetrad constructed by applying Gramm-Schmidt’s algorithm to the coordinate axes of of the grid on which we evolve the moment’s equations. In curved spacetime, this is not a particularly meaningful vector.

We note that, after that transformation and summing over all bins, the total energy of the neutrinos created within a cell is, on average and as measured by an inertial observer, . By comparing this result with the definition of the energy integrated emissivity used in the evolution of the moments equations, we find that, unsurprisingly, . The emissivity in the gray two-moment transport can thus be obtained directly from the tabulated emissivities used by the MC algorithm.

#### Initialization of a cell

The MC algorithm also needs a prescription for the initialization of packets within a cell at the beginning of a time step. Such initialization is required at the beginning of the first time step, as well as in some optically thick regions discussed in Sec. 2.3.5. In optically thick regions, we want these packets to sample the equilibrium distribution of neutrinos. In optically thin regions, we do not have any particularly good guess to rely on, and thus do not create any neutrino packets at the initial time. We assume that the number density of neutrinos, on the initial slice and within the energy bin is, in the fluid rest frame,

(22) |

with . The constant is semi-arbitrary, and sets the initial energy density of neutrinos in regions of moderate optical depth. At and for the problems presented here, we choose , where is a length scale comparable to the typical length scale for variations of . In our tests, we use , with (in code units) for idealized tests and for our test using a core-collapse fluid profile.

A slight complication when sampling the neutrino distribution function on a spatial slice, also discussed in Ryan et al. (2015), is that the 3-volume is not a relativistic invariant. Accordingly, one has to be careful when sampling in the inertial frame a distribution function which is known only in the fluid rest frame. We can rely on the invariance of and the fact that and in the fluid rest frame to derive the ratio between the volume element in the fluid rest frame, , and the volume element in the grid frame, : . Using this result, we sample the distribution of neutrinos by creating neutrino packets

(23) |

Each packet represents neutrinos with fluid rest frame energy . The non-integer part of is treated as a probability to emit one more packet, the initial position of each packet is drawn from a uniform distribution in , and the momentum of the neutrinos is drawn from an isotropic distribution in the fluid frame, as for the main emission procedure. Each packet represents neutrinos, and not neutrinos, to account for the ratio .

#### Propagation, absorption, and scattering

In our code, the evolution of MC packets currently involves three types of operation: the propagation of packets along null geodesics, as well as absorption and scattering of packets sampling neutrino-matter interactions. To evolve a neutrino packet by a time interval , we first determine whether the packet is free-streaming, or whether it is absorbed or scattered by the fluid. Absorption and scattering probabilities can be computed from the infinitesimal optical depth along a geodesic, , with the increment in the affine parameter (). The time interval before the first absorption/scattering is then

(24) |

with drawn from a uniform distribution in . We then determine the smallest of the three time intervals . If is the smallest time interval, the packet is propagated by without interacting with the fluid. If is the smallest interval, the packet is propagated by and then absorbed (i.e. removed from the simulation). Finally, if is the smallest interval, the packet is propagated by , then scattered by the fluid. After scattering we begin a new time step with . Scattering is performed in the fluid rest frame. As we only consider isotropic elastic scattering, we simply redraw the 4-momentum (at constant fluid rest frame energy ), from the same istropic distribution as during packet creation.

Propagation of a packet along null geodesics is performed following the prescription of Hughes et al. (1994),

(25) | |||||

(26) |

We use the same second-order Runge-Kutta time stepping as for the two-moment algorithm, with the caveat that we use metric quantities evaluated at the cell center closest to the initial position of the packet. We only switch the cell from which we gather the fluid and metric variables at the end of a time step. This naturally imposes a limit on the timestep , where reflects our tolerance for how far packets potentially move into a neighboring cell before we use updated values of the metric, and is the maximum grid-coordinate value of the speed of light. In most tests, we choose , which is comparable to the Courant condition for the evolution of the fluid and of the moments .

#### Optically thick regions

In theory, we now have at our disposal all of the pieces required for a basic MC scheme: initialization of the simulation, emission of neutrino packets, and interactions with the fluid. We could move on to the computation of the backreaction on the evolution of the moments. However, this would be prohibitively expensive for at least two reasons. First, the scheme would continuously create and absorb a large number of packets in optically thick regions, where a very small fraction of the emitted packets survive any single time step. Second, in regions in which the scattering opacity , the code would spend too much time propagating neutrino packets over a large number of time intervals . Both of these issues can, however, be handled through much cheaper methods.

To do so, we define a 4-dimensional grid in , with the spatial discretization being provided by the grid on which we evolve the moments, and the energy discretization provided by the binning used when generating the NuLib table for neutrino-matter interactions. We define three types of cells on this grid. Cells with are optically thick cells. There, we assume that the neutrinos are in thermal equilibrium with the fluid. Cells which are not optically thick and satisfy (with the time step in the fluid rest frame) are called high-scattering cells. There, we assume that the diffusion equation is a good approximation to the evolution of the energy density of neutrinos. Finally, cells which are neither optically thick nor high-scattering use the standard MC algorithm outlined in the previous sections. The parameters and can be specified at run time. It is important to note that for our algorithm to formally converge to a solution of Boltzmann’s equations, we should take and , in addition to increasing the number of neutrino packets and the accuracy of the evolution of the moments. For the grid spacings and number of packets that we can currently afford, however, the errors introduced by the use of approximate methods in optically thick / high-scattering cells are subdominant.

As our algorithm couples the MC evolution to the evolution of the moments, we have a fairly simple solution to the costly evolution of optically thick cells. We simply skip the MC algorithm in optically thick regions, and rely solely on the two-moment scheme in that regime. The two-moment scheme is known to perform well in the diffusion regime, as long as some corrections are applied to the fluxes in order to limit the impact of numerical dissipation. A cell is masked if . Any neutrino packet which ends a time step in a masked cell is removed from the simulation. When coupling the MC and moment schemes, we assume that neutrinos are in thermal equilibrium with the fluid in masked cells. The proper choice for depends on the typical length scale over which vary. We note that we only assume a thermal distribution of neutrinos for the computation of high-order moments and of . The evolution of the moments computes the diffusion of neutrinos through optically thick regions without requiring exact thermal equilibrium.

Entirely neglecting masked cells may create errors close to the boundary of the masked region. To provide an effective boundary condition to the MC algorithm, we evolve any masked cell that shares some boundary with an evolved cell (in 4-dimensional space (), and including 1D, 2D, and 3D interfaces), using a modified MC algorithm. All packets that finish a time step in a masked boundary cell are removed from the simulation. Packets sampling an equilibrium distribution of neutrinos are then redrawn at the beginning of each time step, as during the initialization of the MC scheme (but setting the constant ). The emission of neutrino packets during a time step proceeds as in evolved cells. However, we modify so that no more than packets are created in a single energy bin of a single boundary cell. In this paper, we use . While destroying and re-creating packets in these cells at each time step may seem highly inefficient, we note that this generally happens in cells where the number of packets required to sample the equilibrium distribution of neutrinos is comparable to the number of packets created during a time step. The additional use of a strict limit on the number of packets created within each of these cells guarantees that our effective boundary conditions does not become an important cost in the MC algorithm.

#### Approximate treatment of high-scattering regions

In high scattering regions, instead of treating each scattering event individually, we rely on the expected diffusion velocity of packets away from their original location in the fluid rest frame. Solving the diffusion equation in spherical symmetry, we get that the probability distribution for a packet to move a distance from its original location in the fluid rest frame, over a time interval in that frame, is

(27) | |||||

(28) |

We could tabulate the function , defined implicitly by

(29) |

and then randomly draw from a uniform distribution in to obtain . In fact, we find that we can obtain a better approximation of the diffusion limit for moderate values of if we set

(30) |

and set whenever that formula predicts , to avoid violating causality. Cases in which are interpreted as packets that do not experience any scattering events, and thus behave as free-streaming neutrinos. Eq. 30 is obtained by renormalizing Eq. 29, using the constraint that the probability for a packet to avoid any scattering (i.e. to get ) is . With that normalization, we find good agreement between the approximate distribution of and the distribution obtained by treating scattering events individually, for optical depths as low as (see Fig. 1).

In practice, when evolving a neutrino packet, we always evolve the packet exactly up to its first interaction. If that interaction is a scattering event, and the remaining time in the time step is such that , we continue the time step using the diffusion approximation. We then draw a new absorption time and evolve the packet by as follows. We randomly draw from a uniform distribution in to obtain . We interpret as follow: we transport the neutrino packet with the fluid for a time , and then propagate the packet for a time in a random direction drawn from an isotropic distribution in the fluid rest frame. If the packet is absorbed, we then remove it from the simulation.

The advection of a packet by the fluid proceeds as follows. We define a momentum , setting and by requiring that and . Then, we propagate the packet along the null geodesic defined by for a time . This propagates the packet over exactly the coordinate distance covered by the fluid during the time interval , but along a null geodesic instead of a timelike curve. Practically,

(31) | |||||

(32) | |||||

(33) |

As these equations are singular for a fluid with , we only go through this procedure when . Otherwise, advection with the fluid motion simply keeps a packet at its current location. This algorithm does not exactly reproduce the diffusion of a packet in curved spacetime, but guarantees that packets experience the proper gravitational redshift, and are advected by the fluid. The choice to evolve along a null geodesic instead of along is done so that we can evolve the momentum of the packet according to Eq. 26, which assumes . We then propagate the packet for along . The vector is timelike, not null, but we make the approximation that the 4-momentum of neutrinos in the packet is constant during that step. In all of the above, the conversion between ‘primed’ time intervals and ‘unprimed’ time intervals is .

After this ‘advection’ step, we move to the ‘free-streaming’ step. We randomly draw a direction of propagation from an isotropic distribution in the fluid frame, as when emitting new packets, keeping the fluid frame energy constant. This generates a new fluid frame momentum for the packet, . There is now a subtlety in the determination of the time interval over which we let the packet free stream in the direction of . It is not simply , because packets propagating in different directions should be propagated for the same time interval in the fluid rest frame. The correct time step for the free propagation of a packet is . With this weighting of the time step, the average propagation velocity of freely propagating packets is , as desired. It is worth noting that because of this weighting, a given packet may be propagated to a time larger or smaller than initially requested. In our MC algorithm, there is however no particular requirement for the propagation of all packets to end at the same time in any given call to the MC algorithm, so this does not create any issue.

The last step in our approximate treatment of regions of high is the most complex. Once a packet has been propagated, we need to determine what its momentum will be at the beginning of the next time step. We define the angles determining the orientation of that momentum vector in the fluid frame, in a spherical polar coordinate system whose polar axis is parallel to the momentum . The momentum at the beginning of the next time step is then, in the fluid frame,

with , the angles defining the orientation of in the fluid-frame orthonormal tetrad defined in Sec. 2.1.

By symmetry, we can always draw from a uniform distribution in . The simplest options for would be to either draw from a uniform distribution in , which effectively means that is isotropically distributed in the fluid frame, or to set and use the average direction of propagation of the packets away from the motion of the fluid as their final direction of propagation. We will see that both options fail to reproduce solutions obtained without using any approximation in regions of high . By inspection of three sets of scattering experiments in which we scattered packets through optical depths , however, we find that to high accuracy the distribution of only depends on the ratio . Furthermore, a good fit to the result of these experiments can be obtained by setting

(34) |

with drawn from a uniform distribution in . To capture the fact that should be isotropic for and that for , we want and . In practice, we fit to match our scattering experiments, and tabulate the result of those fits. We then interpolate linearly in . We use with

(35) | |||||

and . This approximation is tested in Sec. 3.4, and shows very good agreement with a full treatment of individual scatterings. Nevertheless, it is only an approximation. When using this prescription, finite errors remain even in the limit of an infinite number of particles. Thus, despite the fact that the accuracy of this scheme is a small source of error compared to numerical errors at the grid resolution / number of MC packets that we can currently afford, we should keep in mind that the coupled Moment-MC scheme only formally converges to an exact solution of Boltzmann’s equation in the limit (i.e. when this approximate method is used nowhere).

### 2.4 MC-Moments Coupling

The role of the MC evolution in our algorithm is to provide missing information about the distribution function of neutrinos to the evolution of the moments. In particular, we have shown in Sec. 2.2 that the moment scheme needs the Eddington tensor , and the energy-averaged absorption and scattering opacities . In order to compute the evolution of the electron fraction of the fluid, we also need the number emission and absorption opacity, and . Indeed, the evolution equation for the composition of the fluid is

(36) |

where the sum is taken over all neutrino species, and is for , for , and for heavy-lepton (anti)neutrinos.

We have already seen that the emissivity is simply , with the integrated emissivity within an energy bin. Similarly, . The other quantities will be derived from ratios of the MC moments , , , , , , and . All of these moments are computed following the same strategy. Any MC moments is divided into three contributions

(37) |

is the contribution from high opacity cells, as defined in Sec. 2.3.5, i.e. cells where we do not evolve MC packets, and simply assume that neutrinos are in equilibrium with the fluid for the purpose of computing . is the contribution from packets propagating without any approximation, and from the portion of the time step that a packet in a high-scattering region (i.e. treated using the approximation of Sec. 2.3.6) spends ‘free-streaming’. Finally, is the contribution from the portion of the time step that those same packets spend being explicitly advected with the fluid. This decomposition is necessary because MC packets are treated very differently (or entirely ignored) in regions of high opacity. It is not, however, equivalent to a more traditional separation between ‘trapped’ and ’free-streaming’ neutrinos. While neutrinos contributing to are trapped, neutrinos contributing to the other two components could be either trapped or free-streaming. In the limit of infinite resources, all packets should contribute to .

Each of these three components relies on a different estimate of the stress-energy tensor, adapted to the treatment of neutrino packets in that region. For , each energy bin which is not evolved by the MC algorithm is assumed to contribute

(38) |

as appropriate for neutrinos in equilibrium with the fluid. Here, is the absorption opacity of the fluid for neutrinos in the b energy bin (computed at the center of the bin). For each propagating packet, we have

(39) |

while for packets advected with the fluid, we use the approximation

(40) |

which assumes an isotropic distribution of momenta in the fluid rest frame. The full moment is then obtained by integrating over the proper volume of the cell (for the optically thick limit), or the propagating time of a packet. For example, for the energy density ,

The first sum is taken over the energy bins in which the cell is assumed to be optically thick, while the other two are taken over the paths of packets either propagating along null geodesics or advected with the fluid. We note that there are several abuses of notation in the previous formula. A single packet undergoing scattering interactions may contribute to both the free-streaming and advected components of the stress-energy tensor during a single time step, with different , and may contribute multiple time to the free-streaming component of the stress-energy tensor with different and . Practically, the calculation is performed by letting each packet contribute to the computation of the moments whenever it is evolved in time. Packets evolved without approximations contribute once to for each time interval spent free-streaming in-between scattering or absorption events, while packets in high-scattering regions evolved using the approximation from Sec. 2.3.6 contribute to both (for the period spent approximately free-streaming) and (for the period spent being advected by the fluid). The optically thick regions are taken into account during the emission step, when we determine which cells are not evolved using the MC algorithm, and thus contribute to .

The other moments can be obtained by taking different projections of the stress tensor, except for which is

With these moments at hand, we can compute , and . We also define the average number of packets per cell (multiplied by the proper time interval spent in that cell) as

To time-average the MC moments over multiple time steps, and thus allow us to use a lower number of packets in the simulation, we continuously add to all of the above moments. Instead of resetting the moments to zero at the beginning of a time step, we damp all moments (including ) using

(44) |

Thus, is a minimum damping timescale in the fluid rest frame, while is the desired number of packets over which we average the neutrino distribution function, assuming that we can accumulate packets over a timescale shorter than . The distance provides a rough estimate of the propagation time across a cell. A large implies small statistical errors in the MC moments, but longer averaging timescales. Our standard choice in this manuscript is , which should lead to relative errors in the MC moments. The actual accuracy of the method will however be better than this in optically thick regions, where more than 100 packets live in a cell at any given time. The impact of varying on the noise in the MC moments is explored in the tests of Sec. 3.

We also need a backup prescription for the computation of the moments when the MC algorithm does not have enough packets to provide any reliable information. Whenever we have a cell with , all opacities are set to zero and the Eddington tensor is computed using the M1 closure, . During the moment evolution, we also use at cell faces if both neighboring cells have . When that is the case, the characteristic speeds of the system are estimated assuming the nonlinear M1 closure. In this work, we use . We note that this prescription is necessary in many idealized problems in which neutrinos may be confined to a small portion of the grid, but probably less important in simulations of neutron star mergers, where neutrinos are present everywhere.

The most natural interpretation of this ‘time-averaging’ of the moments is that we damp the contribution of a packet to over the time scale needed for packets to cross a given grid cell, or over the user-prescribed time scale (whichever is smaller). The damping time scale thus varies from grid cell to grid cell, and over time. A disadvantage of this prescription is that we do not have absolute measurements of the moments, as cannot be easily normalized to provide a true average of . In our algorithm, only ratios of the moments are meaningful quantities to use. Thus, we can use the Eddington tensor , but not the energy density itself. We have found, however, that a fixed damping timescale (which would allow us to properly normalize ) can lead to instabilities in optically thick regions (if the time scale is too long), or to a highly inaccurate closure far away from the sources (if the time scale is too short). As any given cell can switch from the former to the latter over the course of a simulation, we need some adaptivity in the choice of the damping time scale.

In SpEC, we evolve the MC packets before evolving the moments, and after evolving the fluid. When information from the MC evolution is required for the evolution of the moments, we always use the last-computed moments . In regions where more than packets cross a cell over a given time step, this means that is simply the average of over all packets evolved during the last MC time step. In other regions, is an average of over the recent past of the system. This is a rather simple choice which, while appropriate for the type of conditions existing in neutron star merger simulations, may be problematic in any system where the neutrino-fluid and/or moment-MC couplings are stiff enough that they have to be treated implicitly (e.g., possibly, in rapidly cooling optically thin regions existing in radiation-dominated accretion disks).

Finally, we will show in the code tests that, at finite resolution, the moment and MC algorithms can become inconsistent. Whether this will be a significant issue for any given application remains an open question, as is the best method to avoid this issue. To illustrate a potential way around this issue, we implement a fairly simple-minded method which attempts to switch to a pure moment evolution when the moment and MC schemes become so inconsistent that using the MC moments may be meaningless. We define a deviation between the MC and moment algorithms as

(45) | |||||

(46) |

The Eddington tensor is then defined as