$\nu\texttt{bhlight}$: Radiation GRMHD for Neutrino-Driven Accretion Flows

TL;DR

The paper introduces $ u$bhlight, a new computational scheme that models neutrino transport in general relativistic magnetohydrodynamics to better understand neutron star mergers and their associated phenomena.

Contribution

It presents $ u$bhlight$, a novel Monte Carlo-based scheme for energy-dependent neutrino transport in 3+1D GRMHD simulations of neutron star merger accretion disks.

Findings

The scheme accurately models neutrino effects in relativistic astrophysical flows.

Tests demonstrate the scheme's effectiveness and necessity for realistic simulations.

Sample calculations show the impact on accretion disk composition and dynamics.

Abstract

The 2017 detection of the in-spiral and merger of two neutron stars was a landmark discovery in astrophysics. We now know that such mergers are central engines of short gamma ray bursts and sites of r-process nucleosynthesis, where the heaviest elements in our universe are formed. In the coming years, we expect many more such mergers. Modeling such systems presents a significant computational challenge along with the observational one. To meet this challenge, we present $\nu\texttt{bhlight}$, a scheme for solving general relativistic magnetohydrodynamics with energy-dependent neutrino transport in full (3+1)-dimensions, facilitated by Monte Carlo methods. We present a suite of tests demonstrating the accuracy, efficacy, and necessity of our scheme. We demonstrate the potential of our scheme by running a sample calculation in a domain of interest---the dynamics and composition of the accretion disk formed by a binary neutron star merger.

Tables

Table 1. Table 1 Emission and Absorption Processes used in ν ​ bhlight 𝜈 bhlight \nu\texttt{bhlight} .

Type	Processes	Charged/Neutral	Corrections/Approximations
Abs./Emis. on Neutrons	${\nu }_e+n\leftrightarrow e^{-}+p$ ${\nu }_{\mu }+n\leftrightarrow {\mu }^{-}+p$	Charged	Blocking/Stimulated Abs. Weak Magnetism Recoil
Abs./Emis. on Protons	${\overline{\nu }}_e+p\leftrightarrow e^{+}+n$ ${\overline{\nu }}_{\mu }+p\leftrightarrow {\mu }^{+}+n$	Charged	Blocking/Stimulated Abs. Weak Magnetism Recoil
Abs./Emis. on Ions	${\nu }_{e}A\leftrightarrow A^{\prime }{e}^{-}$	Charged	Blocking/Stimulated Abs. Recoil
Electron Capture on Ions	$e^{-}+A\leftrightarrow A^{\prime }+{\nu }_e$	Charged	Blocking/Stimulated Abs. Recoil
$e^{+}-e^{-}$ Annihilation	$e^{+}{e}^{-}\leftrightarrow {\nu }_i{\overline{\nu }}_i$	Charged + Neutral	single-$\nu$ Blocking Recoil
$n_i$-$n_i$ Brehmsstrahlung	$n_i^1+n_i^2\to n_i^3+n_i^4+{\nu }_i{\overline{\nu }}_i$	Neutral	single-$\nu$ Blocking Recoil

Table 2. Table 2 Parameters of the initial data for our black-hole-disk system.

Parameter	Value	Units
$M_{BH}$	2.80	$M_{\odot }$
$a_{BH}$	0.80	$J_{BH}/M_{BH}^2$
$s$	4.00	$k_b/$baryon
$Y_e$	0.10	$n_e/\left(n_e+n_p\right)$
$M_{\text{disk}}$	0.05	$M_{\odot }$
$R_{\text{in}}$	3.70	$M_{BH}$
$R_{\text{max}}$	9.03	$M_{BH}$

Table 3. Table 3 Numerical parameters used for the post-merger disk calculation.

Parameter	Value
$N_1$	192
$N_2$	128
$N_3$	66
$N_{\nu }/$cell	$5$
$N_t$	$1.6\times {10}^6$
$R_{\text{out}}$	${10}^3{M}_{BH}$

Equations

$${\partial }_t\left(\sqrt{-g}{\rho }_0{u}^t\right)+{\partial }_i\left(\sqrt{-g}{\rho }_0{u}^i\right)$$

$${\partial }_t\left[\sqrt{-g}\left(T_{\nu }^t+{\rho }_0{u}^t{\delta }_{\nu }^t\right)\right]+{\partial }_i\left[\sqrt{-g}\left(T_{\nu }^i+{\rho }_0{u}^i{\delta }_{\nu }^t\right)\right]$$

$${\partial }_t\left(\sqrt{-g}{B}^i\right)-{\partial }_j\left[\sqrt{-g}\left(b^j{u}^i-b^i{u}^j\right)\right]$$

$${\partial }_t\left(\sqrt{-g}{\rho }_0{Y}_e{u}^t\right)+{\partial }_i\left(\sqrt{-g}{\rho }_0{Y}_e{u}^i\right)$$

$$T_{\nu }^{\mu }=\left({\rho }_0+u+P+b^2\right)u^{\mu }{u}_{\nu }+\left(P+\frac{1}{2}{b}^2\right){\delta }_{\nu }^{\mu }-b^{\mu }{b}_{\nu }$$

$$B^i{=}^{\ast }{F}^{it}$$

$${}^{\ast }{F}^{\mu \nu }=b^{\mu }{u}^{\nu }-b^{\nu }{u}^{\mu }.$$

$$P=P(\rho ,u,Y_e).$$

$$P$$

$$\frac{u}{\rho }:=\epsilon$$

$$1-Y_e$$

$${\nu }_i+p$$

$${\nu }_i+n$$

$${\nu }_i+A$$

$${\nu }_i+\alpha$$

$$-k^{\mu }{\eta }_{\mu }=ϵ=h\nu ,$$

$$\frac{D}{d\lambda }\left(\frac{h^3{I}_{ϵ,f}}{{ϵ}^3}\right)=\left(\frac{h^2{\eta }_{ϵ,f}}{{ϵ}^2}\right)-\left(\frac{ϵ{\chi }_{ϵ,f}}{h}\right)\left(\frac{h^3{I}_{ϵ,f}}{{ϵ}^3}\right),$$

$${\chi }_{ϵ,f}={\alpha }_{ϵ,f}+{\sigma }_{ϵ,f}^a$$

$${\eta }_{ϵ,f}=j_{ϵ,f}+{\eta }_{ϵ,f}^s(I_{ϵ,f})$$

$$e_{(a)}^{\mu }$$

$$e_{(a)}^{\mu }{e}_{\nu }^{(b)}={\eta }_{\nu }^{\mu }{\eta }_{(b)}^{(a)}$$

$$e_{(t)}^{\mu }=u^{\mu },$$

$$G_{(a)}=\frac{1}{h}\int dϵd\Omega \left({\chi }_{ϵ,f}{I}_{ϵ,f}-{\eta }_{ϵ,f}\right)n_{(a)},$$

$$G^{\mu }=e_{(a)}^{\mu }{G}^{(a)}.$$

$$G_{ye}=\frac{m_p}{h}\text{sign}(f)\int \frac{{\chi }_{ϵ,f}{I}_{ϵ,f}-{\eta }_{ϵ,f}}{ϵ}d\Omega dϵ,$$

$$\text{sign}(f)=\{\begin{array}{ll}1& \text{if }f={\nu }_e\\ -1& \text{if }f={\bar{\nu }}_e\\ 0& \text{if }f={\nu }_x\end{array}$$

$$\frac{1}{\sqrt{-g}}\frac{d{N}_s}{d^{3}xdtd\nu d\Omega }=\frac{1}{w\sqrt{-g}}\frac{dN}{d^{3}xdtd\nu d\Omega }=\frac{1}{\omega }\frac{j_{ϵ,f}}{h\nu },$$

$$N_{s,tot}=\Delta t\sum _{f\in \{{\nu }_e,{\nu }_{\bar{e}},{\nu }_x\}}\int \sqrt{-g}{d}^{3}xd\nu d\Omega \frac{1}{w}\frac{j_{ϵ,f}}{h\nu },$$

$$N_{s,i,f}=\Delta t{\Delta }^{3}x\int \sqrt{-g}d\nu d\Omega \frac{1}{w}\frac{j_{ϵ,f}}{h\nu }.$$

$$w=\frac{C}{\nu }$$

$$E_s=wh\nu =hC;$$

$$C=\frac{\Delta t}{h{N}_{target}}\sum _{f\in \{{\nu }_e,{\nu }_{\bar{e}},{\nu }_x\}}\int \sqrt{-g}{d}^{3}xd\nu d\Omega j_{ϵ,f}.$$

$$\Delta {\tau }_a(\nu ,f)=\nu {\alpha }_{ϵ,f}\Delta \lambda ,$$

$$\Delta {\tau }_a(\nu ,f)>-\ln (r_a),$$

Full text

1. Introduction

We now know that the in-spiral and merger of two neutron stars is a central engine of short gamma ray bursts Soares-Santos et al. (2017); Eichler et al. (1989); Narayan et al. (1992) and a site of r-process nucleosynthesis Abbott et al. (2017a), where the heaviest elements in our universe are formed Lattimer and Schramm (1976); Lattimer et al. (1977); Côté et al. (2018). In the coming years, many more such mergers are expected Abbott et al. (2017b).

This breakthrough poses a number of questions. What are the dynamics driving the gamma ray burst? Is the relativistic burst of material out the poles driven by neutrino annihilation Jaroszynski (1996) or magnetic fields Blandford and Znajek (1977)? What fraction of these jets escapes and what fraction is slowed down by ambient material Mooley et al. (2018)? What fraction of the r-process nucleosynthetic yields comes from material in the tidal tails of the merging stars and what fraction comes from wind driven off of material accreting onto the central remnant Tanvir et al. (2017)?

This last question is of particular importance for understanding the spectrum of the optical and infrared afterglow of the merger event Tanvir et al. (2017). The heavy elements produced via r-process nucleosynthesis radioactively decay, producing this afterglow—the so-called macronova or kilonova . The remnant accretion disk ejects mass as a wind, which may be thermal, magnetically-driven, or neutrino-driven.111Fernández and Metzger (2016) and references therein offer a nice summary of these processes. Along with the tidal tails in the merger, this wind may be one site of r-process nucleosynthesis.

The dynamics of the r-process in the wind depend on its composition, which depends on lepton number, and thus neutrino processes and transport. The mass and morphology of the wind depends on the dynamics of the disk, which depends on magnetically-driven turbulence via the magneto-rotational instability Balbus and Hawley (1991), neutrino transport, and general relativistic effects such as frame dragging Wald (2010). Therefore, accurately computing the nucleosynthetic yields produced—and thus the spectrum of the kilonova—depends sensitively on the interplay of gravity, plasma physics, and neutrino radiation transport. In other words, they are well-modeled by general relativistic radiation magnetohydrodynamics (GRRMHD).

Although black hole accretion disk physics is a large and well-explored topic,222See Abramowicz and Fragile (2013) and references therein for a review. very few three-dimensional (3D) calculations of accretion disks formed by a compact binary merger including neutrino physics have been performed, and those only recently. (Indeed, few GRRMHD simulations of disks have been performed in any context.) Foucart et al. (2015) and Hossein Nouri et al. (2018) used a moment method to treat the radiation in a disk formed by the merger of a black hole and a neutron star. Siegel and Metzger (2018) modeled a disk formed by the merger of two neutron stars with general relativistic magnetohydrodynamics on a Cartesian grid and a leakage scheme for the neutrinos. Siegel et al. (2018) use a similar calculation to argue that r-process nucleosynthesis can occur in disks formed by the collapse of massive stars. And Fernández et al. (2018) performed a suite of studies of the disk outflow with a cooling function treatment for the neutrino physics.

These groundbreaking efforts, although heroic, make significant approximations in the treatment of the radiation transport. Realistic modeling of neutrino transport requires solving the $(6+1)$-dimensional kinetic Boltzmann equation for each neutrino species, which is computationally expensive and numerically challenging. In the limit of infinite optical depth, the radiation field may be treated with diffusion physics—see, e.g., Miralles et al. (1993). For vanishing optical depth, a cooling or leakage scheme, where neutrinos are allowed to freely stream through material is appropriate. Indeed this is the approach taken by Siegel et al. (2018) and Fernández et al. (2018). For intermediate optical depths, the transport equations must be solved directly. Moment-based schemes, where the continuum limit of the radiation field is taken and a set of hydrodynamic-like equations are attained, side-step this requirement by imposing strong assumptions on the radiation field in order to close the system of equations, with poorly understood consequences on modeling accuracy. This approach was used by Foucart et al. (2015) and Hossein Nouri et al. (2018).333Much progress has also been made in the postprocessing of simulations of accretion disks with realistic transport. See, for example, Richers et al. (2015) and Foucart (2018).

We seek to the remedy this gap. We present $\nu\texttt{bhlight}$, a new GRRMHD code with accurate neutrino transport via Monte Carlo methods. Monte Carlo methods solve the full kinetic Boltzmann equation, discretized with particles, each of which represents a packet of radiation. $\nu\texttt{bhlight}$ is built on the successful photon GRRMHD code bhlight Ryan et al. (2015) and is designed specifically to tackle the post-merger disk problem.

In section 2, we describe in detail the system of equations $\nu\texttt{bhlight}$ is designed to solve. In section 3, we describe the methods used. In section 4 we describe code tests used to verify $\nu\texttt{bhlight}$. In section 5, we demonstrate our new capabilities with an example calculation of a post-merger disk in full 3D with realistic neutrino transport. Finally, in section 6, we offer some concluding thoughts.

2. System

We solve the equations of relativistic ideal MHD coupled to neutrino radiation. We use a formulation almost identical to that presented in Gammie et al. (2003); Dolence et al. (2009); Ryan et al. (2015). However, there are a few key differences. We evolve the conserved lepton number density (encapsulated in the electron fraction $Y_{e}$). Unlike Gammie et al. (2003), this necessitates a realistic, tabulated equation of state. Unlike in Dolence et al. (2009); Ryan et al. (2015), our radiation is relativistic neutrinos, not photons. Our radiation sector carries conserved lepton number as well as energy and momentum. Moreover, while there is only one “type” of photon, there are several flavors of neutrinos, which we bundle into three types: electrons, anti-electrons, and heavies. We discuss the details of these differences below.

2.1. Fluid

The fluid sector consists of the following system of equations.

{widetext}

[TABLE]

where the energy-momentum tensor $T^{\mu}_{\ \nu}$ is assumed to be

[TABLE]

for metric $g_{\mu\nu}$, rest energy $\rho_{0}$ fluid four-velocity $u^{\mu}$, internal energy density $u$, pressure $P$, and Christoffel connection $\Gamma^{\alpha}_{\beta\gamma}$.

Equation (1) represents conservation of baryon number. Equation (2) represents conservation of energy-momentum, subject to the radiation four-force $G_{\nu}$ (not to be confused with the Einstein tensor). Note that we have added a multiple of equation (1) to the $\nu=0$ index of the canonical form of the energy-momentum equation to arrive at equation (2). This is equivalent to the canonical form, but removes rest-energy from the energy conservation law. We have found this approach to be more numerically favorable. See Martí et al. (1991) for one influential work that uses this trick.

Equation (3) describes the evolution of magnetic fields, where

[TABLE]

comprise the magnetic field components of the Maxwell tensor $F_{\mu\nu}$ and $b^{\mu}$ is the magnetic field four-vector

[TABLE]

Finally, equation (4) describes the conservation of lepton number. $G_{ye}$ is a source term describing the rate at which lepton density is transferred between the fluid and the radiation field. It will be described in more detail below.

The system is closed by an equation of state, which relates the pressure $P$ to the density $\rho$, internal energy $u$, and and electron fraction $Y_{e}$:

[TABLE]

We use an equivalent, temperature-dependent formulation of the equation of state, which relates the pressure $P$ and specific internal energy $\varepsilon=u/\rho$ to the density $\rho$, temperature $T$, and electron fraction $Y_{e}$:

[TABLE]

We invert equation (10) to find the temperature and then calculate the pressure using equation (9).

2.2. Neutrino Physics

We are interested in r-process nucleosynthesis, which depends on the fraction of free neutrons in our gas, or

[TABLE]

The electron fraction $Y_{e}$ is affected by the emission or absorption of electron neutrinos (denoted $\nu_{e}$) and electron antineutrinos (denoted $\bar{\nu}_{e}$) but unaffected by the emission and absorption of all other neutrinos. We therefore dub the other neutrinos, which do not modify electron number heavy neutrinos and denote them $\nu_{x}$. If neutrino species does not matter (or we wish to iterate over species, depending on context), we denote the neutrinos as $\nu_{i}$.

We include many different interactions of neutrinos with matter. We categorize these processes as “absorption or emission” and as “scattering” processes. We list the absorption and emission processes in table 1. Those that involve the absorption or emission of an electron neutrino or antineutrino can change the electron fraction and therefore the number of free neutrons in the gas. We include the elastic scattering processes listed below,

[TABLE]

where $n$ represent neutrons, $p$ protons, $\nu_{i}$ neutrinos of arbitrary type, $A$ heavy ions, and $\alpha$ alpha particles.

There are several effects which we are neglecting, mainly inelastic scattering of neutrinos off of electrons Bruenn (1985); neutrino-neutrino annihilation, which may help drive the gamma ray burst Eichler et al. (1989); and neutrino oscillations Duan et al. (2011). We also neglect ion screening, electron polarization, and form factor corrections to neutrino-heavy ion scattering (13). Although this effect is subdominant in core-collapse supernovae Bruenn and Mezzacappa (1997), we do not know how important it is for the disk problem. We also note that the pair processes, i.e., nucleon-nucleon bremsstrahlung and particle-antiparticle annihilation, do not impose any conditions on pairs of Monte Carlo radiation packets. Rather, they are approximated as isotropic processes and incorporated into our emissivities and absorption opacities. On large scales, i.e., $GM_{BH}/c^{2}$ for a black hole of mass $M_{BH}$, we believe this is a good approximation. Since this work is a “first pass” at accurately tracking neutrino physics in neutrino driven accretion flows, we believe ignoring these effects initially is justified. We will incorporate and study them in future work.

Neutrino interactions with matter have a long history in astrophysics Freedman (1974); Tubbs and Schramm (1975); Fuller et al. (1982); Bruenn (1985); Leinson et al. (1988); Aufderheide et al. (1994); Horowitz (1997). We borrow these results to produce our emissivities, opacities, and cross sections. Our emissivities and opacities in particular are drawn from tabulated data first presented in Burrows et al. (2006), which also accounts for subdominant high-density many-body effects. Scattering is treated on an interaction-by-interaction basis and we use analytic single-particle cross-sections. We use cross-sections as summarized in Burrows et al. (2006).

2.3. Treatment of the Neutrino Radiation Field

We assume our neutrinos are massless and travel on null geodesics and obey a light-like dispersion relation

[TABLE]

where $h$ is Planck’s constant, $\epsilon$ is the energy of a neutrino with wavevector $k^{\mu}$ as measured by an observer traveling along a timelike Killing vector $\eta^{\mu}$. Here $\nu$ is the frequency of the neutrino. However, to avoid notational confusion, we will usually use $\epsilon$ rather than $\nu$ when referring to neutrino energies and frequencies, which are interchangeable via a factor of Planck’s constant. Since the neutrino mass is both small and unknown—far smaller than the many-MeV energies neutrinos attain in post-merger disks—we believe this is a reasonable approximation.

We thus recast our neutrino transport as the standard radiative transfer equation

[TABLE]

where $D/d\lambda$ is a derivative along a neutrino trajectory in phase space, $I_{\epsilon,f}$ is the intensity of the neutrino field of flavor $f\in\{\nu_{e},\bar{\nu}_{e},\nu_{x}\}$,

[TABLE]

is the extinction coefficient that combines absorption coefficient $\alpha_{\epsilon,f}$ and scattering extinction $\sigma_{\epsilon,f}^{a}$ for scattering interaction $a$ and

[TABLE]

is the emissivity combining fluid emissivity $j_{\epsilon,f}$ and emission due to scattering from $\eta_{\epsilon,f}^{s}$. Note that every neutrino flavor has its own radiation field and interactions with matter. (Equivalently, the radiation field has an extra, discrete index specifying neutrino flavor.) Each of the quantities in expression (16) is invariant.

2.4. Radiation-Fluid Interactions

We define an orthonormal tetrad444Here roman Greek indices indicate a lab frame and latin indices in parentheses indicate a comoving frame.

[TABLE]

with

[TABLE]

so that

[TABLE]

i.e., so that it is comoving with the fluid. In this frame, the radiation four-force is

[TABLE]

where $n_{(a)}=p_{(a)}/\epsilon$. A coordinate transformation then maps the comoving radiation four-force into the lab frame:

[TABLE]

The scalar source term $G_{ye}$ for lepton conservation (4) is similar. Evaluated in the fluid frame, it is given by

[TABLE]

where $m_{p}$ is the mass of a proton and

[TABLE]

determines the sign of the contribution.

3. Methods

3.1. Fluid Integration

We evolve our fluid via a standard second-order conservative high-resolution shock capturing finite volume method. We base our implementation in this sector on HARM Gammie et al. (2003) and use the same set of primitive and conserved variables as described in Gammie et al. (2003) and Ryan et al. (2015). The lone exception being the electron fraction. We describe our implementation of the electron fraction in more detail in section 3.3.

We use a local Lax-Friedrichs (LLF) approximate Riemann solver Harten et al. (1983). For reconstructions, we use either a linear reconstruction with a monotized central slope limiter Toro (2013) or a fifth-order WENO reconstruction Liu et al. (1994). The form we use is the variant described in Tchekhovskoy et al. (2007), although we do not manually reduce the order of reconstruction near discontinuities.

We treat our magnetic fields via a constrained transport method described by Toth Tóth (2000). This version of constrained transport uses cell-centered magnetic fields. We use a special second-order derivative operator which ensures that a discrete, corner-centered divergence of the magnetic field vanishes. Since this scheme uses centered-differencing, it neglects emf upwinding, which can be important for flux loop advection. For more details, see Gammie et al. (2003).

In general relativity, the conversion between conserved variables and primitive variables is not known analytically and involves the numerical root finding of a complex algebraic function. We use the procedure described by Mignone and McKinney in Mignone and McKinney (2007).

3.2. Radiation Transport

We have modified the way bhlight   performs radiation transport. We transport three types of radiation packet, each one corresponding to electron neutrinos, electron antineutrinos, or heavy neutrinos. Each type of neutrino has separate emissivities and opacities. This implies that the probability that a given radiation packet is emitted, scattered, or absorbed depends on the neutrino type. There are two ways we could account for this:

Treat each neutrino type as a separate class of radiation object and draw probabilities from completely separate probability distributions, one for each neutrino type. 2. 2.

Draw probabilities from a joint probability distribution, which depends on neutrino type.

The former allows for more fine-grained control over how phase space is sampled, while the latter has the advantage of being simpler to implement. Because of its simplicity, we have implemented option 2. We now describe this approach in more detail. Our treatment closely follows that in Dolence et al. (2009) with a few modifications. Here we emphasize the differences between our algorithm and that described in Dolence et al. (2009) and Ryan et al. (2015).

3.2.1 Emissivity

The probability distribution of emitted radiation packets is given by

[TABLE]

where $N_{s}$ is the number of “superneutrinos,” or radiation packets with $w$ physical neutrinos per packet, $N$ is the number of physical neutrinos, $j_{\epsilon,f}$ is the emissivity (in the plasma frame) of neutrinos of species $i$ with de Broglie frequency $\nu$ and flavor $f\in\{\nu_{e},\bar{\nu}_{e},\nu_{x}\}$, and $h$ is Planck’s constant.

This implies that the number of superneutrinos created in a time interval $\Delta t$ is

[TABLE]

and the number of superneutrinos of flavor $f$ created in a finite volume cell $i$ of volume $\Delta^{3}x$ is given by

[TABLE]

We approximately fix the total number of superneutrinos created per timestep by setting the weights $w$.

We set the weight to

[TABLE]

where $C$ is a constant.

This ensures that a superneutrino of frequency $\nu$ and $w(\nu)$ has a total energy

[TABLE]

superneutrino energy is independent of frequency.

We calculate the constant $C$ by fixing the total number of superneutrinos created to be $N_{target}$ and inverting equation (26) for $C$. We decide $N_{target}$ by trying to keep the total number of superneutrinos constant over time. This means $N_{target}$ is chosen so that the superneutrinos created and scattered replace those lost to absorption or that leave the domain. This results in the integral quantity

[TABLE]

To summarize, to produce superneutrinos, we sample them from a species dependent probability distribution which has a weight calculated by integrating over the total probability distribution for all species.

3.2.2 Absorption

Our treatment of absorption is identical to that described in Ryan et al. (2015), except that absorption extinction coefficients are now evaluated per neutrino species. Absorption is treated probabilistically. If a radiation packet of neutrino flavor $f$ travels an affine distance $\Delta\lambda$, it passes through an incremental optical depth to absorption

[TABLE]

where $\alpha_{\epsilon,f}$ is the absorption extinction coefficient for neutrino radiation of flavor $f$ and frequency $\nu$. An absorption event occurs if

[TABLE]

where $r_{a}$ is a random variable sampled uniformly from the interval $[0,1)$.

3.2.3 Scattering

Like in bhlight, scattering in $\nu\texttt{bhlight}$ is treated probabilistically. We generalize the approach in bhlight to treat scattering of radiation off of multiple scattering particles, each with their own cross-section. We allow our neutrinos to scatter elastically off of protons (11), neutrons (12), heavy nuclei (13), and alpha particles (14). We calculate the individual number densities of the constituent particles via the appropriate mass fraction, which is tabulated in our equations of state.

For each neutrino flavor $f$ and type $p$ of gas particle off of which a neutrino can scatter, we construct a scattering extinction coefficient $\alpha_{s}(\nu,f,p)$. Then a superneutrino of flavor $f$ scatters off of a particle of species $p$ if

[TABLE]

where $\Delta\tau_{p}(f)$ is the scattering optical depth due to an interaction between the neutrino and the particle, constructed analogously to the absorption optical depth (31), and $b_{s}(\nu,f,p)$ is a bias parameter that enhances the probability of scattering. To ensure that the biased process reflects nature, we reset the weight of the scattered superneutrino to $w/b$ for a conservative process with incident superneutrino of weight $w$. For more details, see Dolence et al. (2009) and Ryan et al. (2015).

How do we sample multiple different interactions? After all, the neutrino should be subject to absorption and scattering against all kinds of particles. When a superneutrino travels an affine distance of $\Delta\lambda$, we construct all optical depths

[TABLE]

and biases

[TABLE]

for absorption and scattering against all particles. Then we sample a uniform random variable

[TABLE]

for each type of interaction. The interaction that occurs is the one for which the ratio

[TABLE]

is smallest for all $i\in\{a,p\}$.

Unlike in Dolence et al. (2009) and Ryan et al. (2015), we allow our bias parameters to depend individually on both the neutrino flavor $f$ and the species of interacting particle $p$. In particular, we demand that the parameter

[TABLE]

be approximately equal for all scattering processes. This ensures that all scattering processes are equally well sampled. In section 4.5, we provide evidence that this procedure is both necessary and effective.

3.2.4 Sampling The Scattered Superneutrino

To generate a new superneutrino from an incident one with wavevector $k^{\mu}$, we follow a modified version of the procedure presented in Dolence et al. (2009); Ryan et al. (2015):

We boost into the rest frame of the plasma. 2. 2.

We sample the four-momentum $p^{\mu}$ of the particle off of which the superneutrino scatters from a thermal relativistic Maxwell distribution using the procedure described in Canfield et al. (1987). Note that this requires the differential single-particle cross section for the interaction of the scattering particle with the neutrino. 3. 3.

We boost into the rest frame of the scattering particle. 4. 4.

We sample the wavevector $k^{\mu}_{s}$ of the scattered superneutrino from the differential single-particle cross section. 5. 5.

We transform $k^{\mu}_{s}$ back into the lab frame.

This is why we must perform scattering on a per scatterer basis. Otherwise, the differential cross section is inaccessible.

3.2.5 Radiation Force

We calculate the radiation four-force on the fluid by conserving four-momentum, in a manner identical to that described in Ryan et al. (2015). The only complication emerges from the tracking of lepton number.

3.2.6 Tracking Lepton Number

When a superneutrino is emitted or absorbed, it can modify the electron fraction of the gas. We couple this contribution to the electron fraction evolution via an operator split update analogous to the radiation four-force update in Ryan et al. (2015). The emission of a neutrino radiation packet of flavor $f$ and weight $w$ provides a discrete contribution to source term (23) of magnitude

[TABLE]

where $u^{0}$ is the time component of the fluid four-velocity, $m_{b}$ is unit mass per baryon in the gas, and $\sqrt{-g}\Delta^{4}x$ is the invariant four-volume of a discrete finite volumes cell (including time step). The contribution to this source term for absorption is equal and opposite.

3.2.7 A Note on Timesteps

We pause briefly to note limits on the size of our timestep. Our method is fully explicit, both in the radiation and fluid sectors. We feel comfortable applying a fully explicit approach because the systems we are interested in have modest optical depths, and cooling times are long. Moreover, our approach to scattering and absorption requires that a superneutrino not travel more than one cell distance in one timestep. This timestep restriction means that the restriction on timestep due to using a fully explicit approach is not so severe.

We therefore insist our timestep is smaller than the following quantities:

The light crossing time within any cell 2. 2.

The cooling time due to emissivity $u/\int d\nu d\Omega j_{\nu}$

The emissivity condition (2) isn’t a guarantee for stability. Rather, it is a guarantee that the system will converge to a stable solution in the limit of large superneutrino number. In practice we find that the light crossing time within a cell is almost always the smallest of these quantities and the cooling time restriction is not severe.

3.3. Advection and Passive Scalars

As discussed in section 2, the electron fraction $Y_{e}$ evolves via equation (4). We have implemented a generic framework for evolving variables which are “passively” advected by the fluid, so called passive variables. Our framework allows for two methods of advection, which differ by what variable is considered the primitive variable. We dub these two approaches advect intrinsics and advect numbers.

3.3.1 Advecting Intrinsics

An extrinsic thermodynamic quantity is one which grows linearly with volume. Extrinsic quantities include total energy and entropy. An intrinsic quantity is the corresponding extrinsic quantity per volume. This is in contrast to, say, specific quantities, which are extrinsic quantities per mass.

Consider an intrinsic variable $\phi$ that is passively advected by the fluid. It obeys the differential equation

[TABLE]

where $u^{\mu}$ is the fluid four-velocity. If we perform a $(3+1)$ split and translate this into the language of finite volumes, $\phi$ is our primitive variable, $\phi u^{0}$ is our conserved variable, and $\phi u^{i}$ is our flux in the $i^{th}$ direction. If

[TABLE]

then we recover equation (4) for conservation of lepton number.

3.3.2 Advecting Numbers

Consider a “number” quantity $X$, which is neither intrinsic nor extrinsic. The “density” $X\rho$ will be intrinsic. This fact suggests an alternative form of equation (37):

[TABLE]

where the conserved variable is now $X\rho u^{0}$ and the $i^{th}$ flux is $X\rho u^{i}$. However, there is a degeneracy in the primitive variable. If we treat $X\rho$ as the primitive, we recover the formalism in section 3.3.1. If we treat $X$ as the primitive, we recover a mathematically equivalent but numerically distinct approach. This is the advect numbers scheme.

If $X=Y_{e}$, then equation (39) becomes equation (4) for the conservation of lepton number. In practice, we have found that advecting $Y_{e}$ as a number as per equation (39) to be more numerically robust—particularly in the atmosphere—than advecting the conserved proton mass density as per equation (37) and this is the approach we use.

3.4. Equation of State

We have implemented realistic, tabulated nuclear equations of state. We use the tables as generated and described in O’Connor and Ott (2010, 10). Our equation of state tables provide thermodynamic quantities in terms of the log of the density, $\log_{10}{\rho}$, the log of the temperature $\log_{10}{T}$, and the election fraction $Y_{e}$.

Since temperature is not one of our primitive (or conserved) variables, we must solve for it via one-dimensional root finding. In our implementation, we use Newton’s method but default to bisection if Newton’s method fails. To do so, we invert the relation

[TABLE]

with a given specific internal energy $\varepsilon$ or the relation

[TABLE]

with a given enthalpy by volume $w$ for $\log_{10}{T}$, where $\log_{10}{\rho}$ and $Y_{e}$ are given by the primitive state. We can then extract thermodynamic quantities such as pressure and sound speed. (We use equation (3.4) when solving for our primitive variables from our conserved variables and (40) everywhere else.)

3.5. Atmosphere Treatment

Equations (1), (2), and (4) are valid only for non-vanishing density $\rho$. Moreover, only some values of internal energy $u$ and electron fraction $Y_{e}$ are physically valid. Therefore, we must impose floors on these quantities to keep them in a physically valid regime.

An additional complication is that tabulated thermodynamic values are available only for a finite range of temperature, pressure, and electron fraction:

[TABLE]

Therefore these limits must be accounted for in some way. The physically allowed values of electron fraction sit well within the range given by equation (3.5), so we simply set floors and ceilings for the electron fraction given by $(Y_{e})_{\text{min}}$ and $(Y_{e})_{\text{max}}$.

For density in black hole metrics, we demand that

[TABLE]

where we choose $\rho_{0}\approx 10^{-5}$ for our disk simulations.555This choice is problem dependent. This implies that near the black hole, our floor is about $10^{-5}$ in code units.666Typically, approximately $10^{5}$ g$/$cm3 in physical units. However, the floor decays in radius, so it is much smaller at large radii. This treatment is designed to ensure that the atmosphere does not interfere with diffuse winds blown off the disk.

For these large radii, $\rho_{\text{flr}}<\log{\rho}_{\text{min}}$. We therefore analytically extend the table with a cold, polytropic equation of state (which depends on electron fraction),

[TABLE]

where $K$ and $\Gamma$ are chosen so that $P_{\text{poly}}$ and $(\partial P_{\text{poly}}/\partial\rho)_{s}$ match the table at $(\log{\rho}_{\text{min}},\log{T}_{\text{min}},Y_{e})$.

For internal energy $u$, we demand that

[TABLE]

where the specific internal energy $\varepsilon$ can contain contributions from binding energy and thus be negative. For consistency with the cold nature of the polytropic equation of state (44), we set $u=u_{\text{flr}}$ when $\rho$ is less than some threshold value, typically a $\log{\rho}_{\text{min}}/10^{5}$ or less. In magnetically dominated regions, these floors are imposed in the rest frame of the fluid. However, in matter or kinetic energy dominated regions, they are imposed in the lab frame.

We note that our treatment of the floors is not thermodynamically consistent. Moreover, any application of density floors is unphysical and can, if care is not taken, change the results of a simulation. Unfortunately, within an Eulerian framework, we have no choice but to apply density floors. Fortunately, this inconsistency affects only very low density regions and should therefore not change the results of a simulation if used judiciously. We have experimented with floor values and found the results of our simulations to be insensitive to these choices.

Electron fraction has no meaning in the atmosphere, but numerically, we must set it to something. For simplicity, we set the atmosphere to have an electron fraction of $Y_{e}=0.5$ at the initial time. When density floors are enforced, $Y_{e}$ is not reset. Rather, the bounds on $Y_{e}$ are enforced independently.

3.6. Tracer Particles

We have added tracer particles to $\nu\texttt{bhlight}$. A tracer particle is a numerical representation of a Lagrangian fluid packet, which is passively advected with the fluid. In the (3+1)-split of general relativity, tracers obey the equation of motion

[TABLE]

where $x^{i}$ are the spatial components of the tracer’s position vector, $\alpha$ is the lapse, $\beta^{i}$ the components of the shift, $v^{i}$ the three-velocity of the fluid, and $u^{\mu}$ the four-velocity of the fluid Foucart et al. (2014).

We interpolate the velocities in equation (46) via second-order Lagrange interpolation. The tracers are integrated in time via a second-order explicit Runge-Kutta scheme and, because we integrate them in lockstep with the fluid, the coupling between tracers and fluid is fully second-order. We utilize the particle infrastructure already implemented in $\nu\texttt{bhlight}$ to treat tracers with the same shared and distributed memory parallelism as for superneutrinos.

In this work, we follow Bovard and Rezzolla (2017) and roughly uniformly sample our tracers everywhere where there is physical fluid (i.e., everywhere that is not atmosphere). When we set up a disk, for each cell containing disk material, we calculate the total mass within the cell, and equally distribute it between $N$ tracer particles that are placed in the cell. The tracer positions within the cell are randomly sampled from the uniform distribution.

4. Code Verification

Many aspects of $\nu\texttt{bhlight}$ have been tested rigorously in previous works, such as Gammie et al. (2003); Dolence et al. (2009); Ryan et al. (2015). Here we discuss tests of the added functionality required to study neutrino transport.

4.1. Advection Tests

To test the passive scalar framework described in section 3.3, we study the advection of an intrinsic (in the thermodynamic sense) scalar field by a constant (in space and time) fluid flow in flat Minkowski space. Under these conditions and in one spatial dimension, equation (37) reduces to777Although here we describe the problem setup for the advection of intrinsic variables, the test (and results) for advection of numbers is identical.

[TABLE]

where $u^{x}$ is the velocity of the fluid in the $x$-direction in Minkowski coordinates.

We solve equation (47) as an initial value problem using the techniques described in section 3.1 with initial conditions

[TABLE]

for scalar fields $\phi_{0}$ and $\phi_{1}$, both obeying equation (47), with periodic boundary conditions on the domain $x\in[-1,1]$. We plot a solution to this initial-boundary value problem at $t=2\sqrt{1-(u^{x})^{2}}/u^{x}$ in figure 1.

We use (49) to check for convergence and (48) to test the handling of discontinuities. A three-dimensional version of this test can be constructed by rotating initial data (49) by 45 degrees about the $y$ and $z$ axes. We plot a two-dimensional slice of $\phi_{1}$ in the three-dimensional test in figure 2 and a two-dimensional slice of the pointwise error in figure 3. Finally, we plot the convergence of infinity norm of the error in $\phi_{1}$ in three dimensions in Figure 4. As expected for a second-order Godunov-type method, our solution converges at second-order.

4.2. Tracer Particles

We test the tracer particle infrastructure described in section 3.6 by constructing known fluid flows and watching the tracers advect with the fluid.

4.2.1 Advection

One simple known flow is that given by equation (47) in section 4.1. In this case—where the velocity field is uniform and constant in time—there is no truncation error in the spatial discretization. Therefore, even for very small resolutions, errors on the order of machine precision can easily be achieved. Indeed, we run this test with a mere 16 cells and achieve machine precision accuracy in the positions of the tracer particles.

4.2.2 Equilibrium Torus

Another simple flow is a torus in hydrostatic equilibrium about a black hole.888The equilibrium torus is in an unstable equilibrium. Eventually a Papaloizou-Pringle instability will form. Fortunately the growth time is long compared to the simulation time presented here Papaloizou and Pringle (1984). For this test, we use the initial conditions as described in section 5.0.1. Briefly, we assume constant entropy and specific angular momentum with a tabulated SFHo equation of state (see section 5.0.1 for details). We use a relatively coarse grid—$96\times 96\times 64$ cells and $\sim 657,000$ tracer particles. We evolve the system for $500$ gravitational times ($GM_{BH}/c^{3}$) with no seed magnetic field.

The continuum initial data is in equilibrium, but the numerical initial data is not. We allow the disk to relax towards numerical equilibrium for $200$ gravitational times, then select tracer particles within $1\times GM_{BH}/c^{2}$ of the midplane of the disk. Figure 5 shows projections of tracks of a random selection of these tracer particles onto the $xy$-plane. The dynamical time in the inner region is shorter than in the outer region, so the innermost trace covers many orbits, while the outermost trace covers only one. Similarly, the outermost region is not yet in numerical equilibrium, hence why the track does not close.999Since integration is not symplectic, we do not expect orbits to completely close, even with a perfectly relaxed disk. However, compared to the effect of disk relaxation, this effect is negligible and it is not visible here.

4.3. Fake Table Tests

To test our tabulated EOS reader, we tabulate the ideal gas law

[TABLE]

where $\Gamma$ is the ratio of specific heats. Using this “fake” table, we can repeat tests presented in Ryan et al. (2015) in the absence of radiation. We perform the non-relativistic linear waves and shock tube tests presented in Ryan et al. (2015).

The tables treat all quantities on a logarithmic scale. For equation (50), the log of all thermodynamic quantities except sound speed is linear and the interpolation is exact, yielding identical results to an analytic EOS. The logarithm of the sound speed is linear at low velocities, but not at relativistic velocities. However, the tests we reproduce from Ryan et al. (2015) are non-relativistic and so these nonlinearities are not present. Therefore, we expect agreement up to machine precision. And indeed, we find this to be the case.

4.4. Artificial Neutrino Cooling

As a basic test of the coupling of neutrinos to matter, we study optically thin neutrino cooling in a simplified context. For this test, we choose to emit either only electron neutrinos or only electron antineutrinos. In each case, we define an emissivity of the form

[TABLE]

where $C$ ensures the units and scale are appropriate,

[TABLE]

is the selection function, and $y_{f}(Y_{e})$ is given by

[TABLE]

Assuming a homogeneous and isotropic fluid at rest, this implies that the electron fraction $Y_{e}$ and internal energy density $u$ obey ordinary differential equations of the form

[TABLE]

where

[TABLE]

Equation (54) has a solution

[TABLE]

where $(Y_{e})_{0}=Y_{e}(t=0)$. This implies that the electron fraction either exponentially decays to zero or exponentially approaches $1/2$, depending on whether we emit electron neutrinos or electron antineutrinos. With this solution in hand, we can solve equation (55) to find that

[TABLE]

where $u_{0}=u(t=0)$.

In our tests, we choose

[TABLE]

so that equation (58) reduces to

[TABLE]

and equation (59) reduces to

[TABLE]

and $u(t)$ asymptotes to $u_{0}-B_{C}/(2A_{C})$. Figure 6 shows the electron fraction as a function of time for a gas cooled by electron neutrinos and figure 7 shows the analogous quantity for a gas cooled by electron antineutrinos. The energy density as a function of time looks identical whether we cool by electron neutrinos or electron antineutrinos. We plot this in figure 8. In all cases, the agreement is good. A small deviation appears at late times, but since the plot is on a log scale, this deviation is extremely small.

4.5. Artificial Single Scattering Events

Our procedure for biasing scattering on a per-interaction basis is a novel part of $\nu\texttt{bhlight}$. Therefore it is worth checking both that it works well and that it is worth the added complexity. We seek to test this here. In this test, we take a homogeneous, zero temperature gas and propagate a steady stream of heavy neutrinos traveling in the $z$-direction through it so that each neutrino traverses a scattering optical depth for the most likely scattering interaction of roughly $\Delta\tau=1$. The optical depth for the less likely scattering interactions will be less than unity, the least likely interaction significantly so. If the resolution is chosen so that the most likely scattering interaction is just barely well-sampled, a naive biasing algorithm will undersample these less likely interactions.

After propagating the neutrinos, we investigate the directions of the neutrinos that have scattered exactly once. The probability distribution of directions traveled by these scattered neutrinos should match the total probability distribution formed by the sum of all differential scattering cross sections

[TABLE]

where $N(\theta,\phi)$ is the total number of heavy neutrinos traveling in the $(\theta,\phi)$ direction and $d\sigma_{p}/d\Omega$ is the differential cross section for a heavy neutrino scattering off of a gas particle of species $p$. We measure $N$ integrated over the azimuthal direction:

[TABLE]

For the purpose of this test, we introduce three fake particles, with three fake, anisotropic, elastic scattering kernels:

[TABLE]

where $\mu=\cos(\theta)$. These interactions have different total cross sections and thus different probabilities that an individual neutrino will scatter via a given process. However, if all processes are well sampled, we should be able to measure a probability distribution that matches equation (62).

We perform this experiment in two ways. First, we use a single global bias which modifies scattering probability uniformly across all interactions. For a resolution which marginally well-samples the most likely interaction ($i=2$), the less likely interactions ($i=0,1$) will be undersampled. Second, we bias each scattering interaction individually, as described in section 3.2.3. This second approach should more evenly sample all interactions for a given resolution. In both cases, we use the same number of superneutrinos (roughly $10^{6}$) and set the biases such that the same number of unscattered superneutrinos scatters each timestep (roughly $4\times 10^{4}$).

We show our results in figure 9. The area of the blue shaded region is the integral $\int(dN/d\Omega)d\Omega$ of equation (61), meaning the boundary of the shaded region is given by equation (63). The green dashed curve is the probability distribution of superneutrinos measured when the experiment is performed with a single global bias. The red solid curve is the probability distribution measured when the experiment is performed with per-interaction biases.

When global biases are used $d\sigma_{0}/d\Omega$ and $d\sigma_{1}/d\Omega$ are undersampled with respect to $d\sigma_{2}/d\Omega$. But when per-interaction biases are used, the agreement with equation (62) is quite good. This indicates both the necessity and efficacy of per-interaction biases.

4.6. Two-Dimensional Lepton Transport

A major motivation for treating neutrino radiation accurately is the fact that neutrinos can carry lepton number from one place to another. In the context of an accretion disk, this means leptons—and thus electron fraction—can be transported from one part of the disk to another.

To demonstrate this capability in $\nu\texttt{bhlight}$, we perform the following simple two-dimensional test. Consider a gas in a two-dimensional, periodic box,

[TABLE]

The gas is at constant density and temperature

[TABLE]

and piecewise-constant electron fraction defined by

[TABLE]

so that there is one “hot spot” of $Y_{e}$ and one “cold spot.” The hot and cold spot regions are separated from the rest of the gas by membranes which are impermeable by the gas but through which neutrinos can travel freely. In this way, the gas does not evolve due to pressure gradients.

Over time, as neutrinos are emitted and absorbed, the electron fractions in the hot spot and the cold spot will come to equilibrium with each other. Indeed, the gas itself will come into equilibrium with the radiation field. Figure 10 shows the electron fraction as a function of space for three times. Figure 11 shows the evolution of the electron fraction as a function of time. The average electron fraction experiences an early transient as leptons are carried into the radiation field but then remains stable. The hot spot and cold spot converge to the average exponentially with time. The final electron fraction is not the average $Y_{e}$ in the initial condition.

4.7. Code Comparisons

The tests described in sections 4.4 and 4.5 use artificial emissivities and scattering cross sections. This has the advantage of permitting an analytic solution against which we can compare. However, it has the significant disadvantage of being unphysical. We would also like to test the performance of our physical emissivities and absorption opacities.

To this end, we compare our code to the supernova code FORNAX Skinner et al. (2018) in two simple test cases. The two codes are designed for different scenarios and use significantly different methods. $\nu\texttt{bhlight}$ is fully general relativistic, whereas FORNAX uses non-relativistic dynamics and an approximate treatment for gravity. $\nu\texttt{bhlight}$ uses Monte Carlo transport for neutrinos, while FORNAX uses a multi-group moment formalism with the M1 closure model Castor (2004). Some care is thus required to choose test cases where both codes converge to the same, physically correct, solution.

We therefore use a simple zero-dimensional setup. We use a homogenous and isotropic gas at rest on a periodic domain in Minkowski space. We also use the same equation of state—SFHo by Steiner et al. Steiner et al. (2013). This eliminates discrepancies due to treatment of the gas. We use the same emission and absorption opacities, as presented in Burrows et al. (2006); Skinner et al. (2018) and provided by Burrows (2018). Since the scattering cross-sections are different between the codes, we disable scattering for these tests. Moreover, by studying only homogeneous, isotropic radiation, we enter a regime where the M1 closure model is valid.

In both our comparison tests, we use the following initial conditions for the gas:

[TABLE]

which roughly mimic conditions one might encounter in a neutrino-driven accretion flow. We run each calculation for total duration of 0.5 seconds. We assume that there is no radiation at the initial time. In both cases we run FORNAX with 200 energy groups and energies ranging from 1 to 300 MeV. We run $\nu\texttt{bhlight}$ with a target number of $10^{5}$ superneutrinos which can have energies in the same range as in FORNAX.

4.7.1 Optically Thin Cooling

In this test, we enable emissivity but disable absorption and scattering. Traces of electron fraction and temperature are shown for both FORNAX and $\nu\texttt{bhlight}$ in figure 12. The cooling rate is a steep function of temperature. As the gas cools rapidly, the electron fraction changes rapidly before reaching equilibrium.

The agreement between FORNAXand $\nu\texttt{bhlight}$is at the percent level, as shown in figure 13. Given that the codes use dramatically different methods, this agreement is quite satisfactory for the problem of interest.

4.7.2 Thermalization

In this test, we enable both emission and absorption but disable scattering. The goal is to watch as the gas and the radiation reach thermal equilibrium. We plot the electron fraction and the temperature for both FORNAX and $\nu\texttt{bhlight}$ in figure 14. The electron fraction changes very rapidly, but the cooling rate is subdued thanks to absorption.

$\nu\texttt{bhlight}$ and FORNAX again agree within about a percent, as shown in figure 15. Again, given that the codes use dramatically different methods, this agreement is quite good.

5. Post-Merger Disk

As a demonstration of $\nu\texttt{bhlight}$’s capabilities, we perform a fully three-dimensional neutrino radiation GRMHD simulation of a representative accretion disk that formed from a compact binary merger. Although we do not perform a detailed analysis, we believe our qualitative results compellingly demonstrate both the capabilities of our code and a need for those capabilities.

5.0.1 Disk Setup

We set up an axisymmetric disk in hydrostatic equilibrium on a Kerr background.101010We perform the setup in Boyer-Lindquist coordinates but transition to the modified Kerr-Schild coordinates described in Gammie et al. (2003) for the evolution. For a detailed discussion of these various coordinate systems, see Poisson (2004) and references therein. We demand our disk have constant fixed specific entropy $s$, specific angular momentum $l$, and electron fraction $Y_{e}$. Under these conditions, the relativistic Euler equations can be written as an exterior differential system, which can be integrated along characteristics for the specific enthalpy $h$ as derived in Fishbone and Moncrief (1976).

Figure 16 shows the specific enthalpy as a function of $\log_{10}{\rho}$ and $\log_{10}{T}$ for the Hempel SFHo Hempel et al. (2012) equation of state with fixed electron fraction $Y_{e}=0.1$. Overlayed on top of the heatmap are contours of constant entropy. To construct a constant entropy disk, we find one of these contours and move along it. Each contour represents a relationship between $\log_{10}{\rho}$ and $\log_{10}{T}$.

The exterior differential system for the enthalpy introduces several constants of integration. These are set by the innermost radius of the disk, $R_{\text{in}}$, the radius of maximum pressure $R_{\text{max}}$, and the limit

[TABLE]

for a given entropy $s$ and electron fraction $Y_{e}$. For ideal gasses, $h_{0}\geq c$. However, this is not the case for more general equations of state. The only constraint is that $h\geq 0$, as required by the weak energy condition. We plot $h-c^{2}$ vs. pressure for $s=4$ $k_{b}/$baryon and $Y_{e}=0.1$ in figure 17. Note the offset along the y-axis.

We initialize the disk with a uniform, weak, purely poloidal magnetic field with a minimum ratio of gas to magnetic pressure

[TABLE]

which acts as the seed field for the magneto-rotational instability Velikhov (1959). We summarize our parameter choices for our disk in table 2. We summarize the numerical parameters used in table 3.

5.0.2 Results

We run our simulation for 10,000 $GM_{BH}/c^{3}$, or $\sim 138$ milliseconds. Figures 18 and 19 show the density of the disk at late time.. Neutrinos can carry lepton number and vary the electron fraction as a function of space and time. Figures 20 and 21 show the electron fraction of the disk-wind system at late times. The core of the disk still has very low electron fraction—close to $Y_{e}\sim 0.1$. However, the outflow has a composition that varies significantly in space. Material in the equatorial plane still has low electron fraction, near $Y_{e}\sim 0.2$. However, material far from the midplane has electron fraction as large as $Y_{e}\sim 0.3$.

We make these qualitative observations more quantitative by examining tracer data. For this analysis, we use tracers that have reached radii greater than a fixed extraction radius $r_{min}=100$ at time $10,000GM_{BH}/c^{3}$. We plot density, temperature, and electron fraction as a function of time for several characteristic tracer particles both in the midplane and off-plane in figure 22. Although we do not calculate yields here, these kinds of tracers may be used as input into a nuclear reaction network to calculate yields. Although not conclusive, we believe these results are highly suggestive that realistic neutrino transport is required to accurately model these systems. We will pursue this detailed modeling in future work.

6. Concluding Thoughts

We have developed the capacity to accurately study neutrino driven accretion flows via Monte Carlo methods. We have validated the accuracy of our approach via a variety of code tests, as discussed in section 4. Moreover, we have demonstrated this capability by performing a fully three-dimensional general relativistic neutrino-radiation-magnetohydrodynamics calculation of a representative black-hole-accretion-disk system formed by a compact binary merger. With realistic neutrino transport active, we observe a rich phenomenology of the wind morphology and composition.

Since we will observe many more compact binary mergers in the coming years, accurately modeling these systems, and their remnants, is of critical importance. We believe $\nu\texttt{bhlight}$ represents a key tool in this modeling effort. In the future, we will use $\nu\texttt{bhlight}$ to investigate the morphology of the disk-wind system in the context of multimessenger observables.

7. Acknowledgements

The authors thank Chris Fryer, Francois Foucart, Daniel Siegel, Oleg Korobkin, Jonas Lippuner, and Roseanne Cheng for many valuable discussions. We are especially grateful to Adam Burrows for providing us with neutrino opacity tables and explaining how these tables are produced in Burrows (2018).

We acknowledge support from the U.S. Department of Energy Office of Science and the Office of Advanced Scientific Computing Research via the Scientific Discovery through Advanced Computing (SciDAC4) program and Grant DE-SC0018297

This work was supported by the US Department of Energy through the Los Alamos National Laboratory. Additional funding was provided by the Laboratory Directed Research and Development Program and the Center for Nonlinear Studies at Los Alamos National Laboratory under project number 20170508DR. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001).

This research used resources provided by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the U.S. Department of Energy National Nuclear Security Administration under Contract No. 89233218CNA000001.

This article is cleared for unlimited release, LA-UR-19-20336.

We are grateful to the countless developers contributing to open source projects on which we relied in this work, including Python Rossum (1995), numpy and scipy van der Walt et al. (2011); Jones et al. (01), and Matplotlib Hunter (2007).