HARM3D+NUC: A New Method for Simulating the Post-merger Phase of Binary Neutron Star Mergers with GRMHD, Tabulated EOS, and Neutrino Leakage

HARM3D (Gammie et al. 2003; Noble et al. 2006, 2009) solves the GRMHD equations in conservative form. HARM3D is a well-tested code that can handle arbitrary coordinate systems, which allows for less numerical diffusion and better conservation of angular momentum when using coordinate systems that more closely conform to local symmetries of the problem (Zilhão & Noble 2014). Below, we set G = c = 1. The GRMHD equations of motion include the baryon conservation equation,

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }\left({n}_{b}{u}^{\mu }\right)=0,\end{eqnarray} \tag{ 1 }$$

the energy–momentum conservation equations (with a heating/cooling source, neglecting momentum transfer)

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }{{T}^{\mu }}_{\nu }={ \mathcal Q }{u}_{\nu },\end{eqnarray} \tag{ 2 }$$

and Maxwell's equations

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\nu }^{* }{F}^{\mu \nu }=0,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\nu }{F}^{\mu \nu }={J}^{\mu },\end{eqnarray} \tag{ 4 }$$

where u^μ is the four-velocity of the fluid, ${ \mathcal Q }$ is the energy change rate per volume in the comoving fluid frame (due to neutrino heating/cooling), n_b is the number density of baryons, F^{μ ν} is the Faraday tensor times $1/\sqrt{4\pi }$, ^* F^{μ ν} is the dual of this tensor or the Maxwell tensor times $1/\sqrt{4\pi }$, and J^μ is the four-current. ²² In practice, we do not use Equation (4), since we work in the limit of ideal MHD. The change in the conservation of lepton number is

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }\left({n}_{e}{u}^{\mu }\right)={ \mathcal R }/{m}_{b},\end{eqnarray} \tag{ 5 }$$

where n_e is the number density of electrons and ${ \mathcal R }\,=-{{ \mathcal R }}_{{\nu }_{e}}+{{ \mathcal R }}_{{\bar{\nu }}_{e}}$ is the difference in the net rate of neutrino and antineutrino number per volume in the comoving fluid frame.

Note that the rest-mass density of the gas (mass per unit volume) is dominated by the baryon mass, ρ ≈ m_b n_b , where m_b is the baryon mass. The baryon number conservation equation can then be replaced by the regular continuity equation:

$$\begin{eqnarray}&&0={m}_{b}{{\rm{\nabla }}}_{\mu }\left({n}_{b}{u}^{\mu }\right)={{\rm{\nabla }}}_{\mu }\left({m}_{b}{n}_{b}{u}^{\mu }\right)={{\rm{\nabla }}}_{\mu }\left(\rho {u}^{\mu }\right).\end{eqnarray} \tag{ 6 }$$

Instead of using n_e and n_b , we may use the fluid density ρ and the electron fraction Y_e :

$$\begin{eqnarray}&&{Y}_{e}\equiv \displaystyle \frac{{n}_{e}}{{n}_{b}}=\displaystyle \frac{{n}_{e}}{\rho /{m}_{b}}=\displaystyle \frac{{m}_{b}{n}_{e}}{\rho }\end{eqnarray} \tag{ 7 }$$

or Y_e ρ = m_b n_e. We can therefore multiply Equation (5) by m_b to yield the electron fraction equation:

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }\left(\rho {Y}_{e}{u}^{\mu }\right)={ \mathcal R }.\end{eqnarray} \tag{ 8 }$$

The total stress-energy tensor is the sum of the fluid part,

$$\begin{eqnarray}&&{T}_{\mathrm{fluid}}^{\mu \nu }=\rho {{hu}}^{\mu }{u}^{\nu }+{{Pg}}^{\mu \nu },\end{eqnarray} \tag{ 9 }$$

and the electromagnetic part

$$\begin{eqnarray}&&{T}_{\mathrm{EM}}^{\mu \nu }={F}^{\mu \lambda }{{F}^{\nu }}_{\lambda }-\displaystyle \frac{1}{4}{g}^{\mu \nu }{F}^{\lambda \kappa }{F}_{\lambda \kappa }\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&=\ | | b| {| }^{2}{u}^{\mu }{u}^{\nu }+\displaystyle \frac{1}{2}| | b| {| }^{2}{g}^{\mu \nu }-{b}^{\mu }{b}^{\nu },\end{eqnarray} \tag{ 11 }$$

where we adopt the ideal MHD condition

$$\begin{eqnarray}&&{u}_{\lambda }{F}^{\lambda \kappa }=0,\end{eqnarray} \tag{ 12 }$$

and where g_{μ ν} is the metric, $h=\left(1+\epsilon +P/\rho \right)$ is the specific enthalpy, P is the pressure, is the specific internal energy density, b^μ = ^* F^{ν μ} u_ν is the magnetic field four-vector, and ∣∣b∣∣² ≡ b^μ b_μ is twice the magnetic pressure P_m .

Equations (2)–(6) can be expressed in flux conservative form

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{U}}\left({\boldsymbol{P}}\right)=-{\partial }_{i}{{\boldsymbol{F}}}^{i}\left({\boldsymbol{P}}\right)+{\boldsymbol{S}}\left({\boldsymbol{P}}\right)\,,\end{eqnarray} \tag{ 13 }$$

where U is a vector of "conserved" variables, F ⁱ are the fluxes, S is a vector of source terms, and P is the vector of primitive variables. Explicitly, these are

$$\begin{eqnarray}&&{\boldsymbol{P}}={\left[\rho ,{{ \mathcal B }}^{k},{\tilde{u}}^{i},{Y}_{e},T\right]}^{T}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{\boldsymbol{U}}\left({\boldsymbol{P}}\right)=\sqrt{-g}{\left[\rho {u}^{t},{{T}^{t}}_{t}+\rho {u}^{t},{{T}^{t}}_{j},{B}^{k},\rho {Y}_{e}{u}^{t}\right]}^{T}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{{\boldsymbol{F}}}^{i}\left({\boldsymbol{P}}\right)=\sqrt{-g}{\left[\rho {u}^{i},{{T}^{i}}_{t}+\rho {u}^{i},{{T}^{i}}_{j},\left({b}^{i}{u}^{k}-{b}^{k}{u}^{i}\right),\rho {Y}_{e}{u}^{i}\right]}^{T}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{\boldsymbol{S}}\left({\boldsymbol{P}}\right)=\sqrt{-g}{\left[0,{{T}^{\kappa }}_{\lambda }{{{\rm{\Gamma }}}^{\lambda }}_{t\kappa }+{ \mathcal Q }{u}_{t},{{T}^{\kappa }}_{\lambda }{{{\rm{\Gamma }}}^{\lambda }}_{j\kappa }+{ \mathcal Q }{u}_{i},0,{ \mathcal R }\right]}^{T},\end{eqnarray} \tag{ 17 }$$

where g is the determinant of the metric, ${{{\rm{\Gamma }}}^{\lambda }}_{\mu \kappa }$ is the metric's affine connection, T is the temperature, and ${B}^{i}={{ \mathcal B }}^{i}/\alpha {=}^{* }{F}^{{it}}$ is the magnetic field.

The primitive velocity is the flow's four-velocity projected into a frame moving orthogonal to the space-like hypersurface:

$$\begin{eqnarray}&&{\tilde{u}}^{\mu }=\left({{\delta }^{\mu }}_{\nu }+{n}^{\mu }{n}_{\nu }\right){u}^{\nu },\end{eqnarray} \tag{ 18 }$$

which only has spatial coefficients

$$\begin{eqnarray}&&{\tilde{u}}^{i}={u}^{i}+\alpha \gamma {g}^{{ti}},\end{eqnarray} \tag{ 19 }$$

where $\alpha =1/\sqrt{-{g}^{{tt}}}$ is the lapse function, βⁱ = − g^ti /g^tt is the shift function, γ = α u^t is the Lorentz factor, and n^μ is the four-velocity of the orthogonal frame: n_μ = [ − α, 0, 0, 0] and ${n}^{\mu }={[1/\alpha ,-{\beta }^{i}/\alpha ]}^{T}$. Defining a fluid three-velocity ${v}^{i}\,={\tilde{u}}^{i}/\gamma$, it can be shown that $\gamma =1/\sqrt{1-{v}^{2}}$, where v² = v_i vⁱ .

In the following section, we describe the implementation of a tabulated EOS in HARM3D+NUC.

The tables and routines for interpolating tabulated quantities are provided by ²³ O'Connor & Ott (2010) and Schneider et al. (2017). The finite-temperature tables give thermodynamic variables, including, for example, the sound speed, as well as the chemical potentials of the nucleons, electrons/positrons and neutrinos/antineutrinos, as a function of the temperature (T), the electron fraction (Y_e ), and the rest-mass density (ρ). The linear interpolation routines are provided by O'Connor & Ott (2010) and Schneider et al. (2017). The interpolation is done in $\mathrm{log}T$, $\mathrm{log}\rho$, Y_e space for $\mathrm{log}\epsilon$, $\mathrm{log}P$, and the rest of the thermodynamical variables.

The tables consider an interpolation between a single nucleus approximation (SNA) in the high-density regime and nuclear statistical equilibrium (NSE) of several nuclides in the low-density regime. The SNA is composed of free nucleons, electrons, positrons, α-particles, and photons. In the high-density regime, nuclei are included using the liquid drop model. The regimes are smoothly interpolated. Using the tables, we have the advantage that the nuclear binding energy release due to recombination energy from the α-particles is included.

There are three main calls to the EOS in HARM3D+NUC:

• 1.  
  We call the EOS when setting the characteristic velocity in order to solve the Riemann problem (Gammie et al. 2003). The wave velocities depend on the relativistic sound speed (Gammie et al. 2003), which can be interpolated directly from the tables.
• 2.  
  We replaced the primitive variable u = ρ with the temperature as a reconstructed variable, which makes the interpolation of the pressure faster because all independent variables are known and can be used to perform the interpolation immediately. This means that we call the EOS to obtain the primitive energy density u after we update ρ, T, and Y_e from the conservation equations.
• 3.  
  We call the EOS repeatedly when converting from conserved variables to primitive variables.

Our implementation of a tabulated EOS into the conserved to primitive variables routine in HARM3D+NUC follows Siegel et al. (2018).

2.2.1. Primary Recovery: 3D Routine

The primary recovery routine follows a three-parameter root-finding method similar to ones implemented in Cerdá-Durán et al. (2008) and Siegel et al. (2018). We call this routine the "3D" routine. For this routine, we reduce the GRMHD equations into three equations that have three unknowns, allowing us to solve the following system:

$$\begin{eqnarray}&&{\tilde{Q}}^{2}=\left(1-\displaystyle \frac{1}{{\gamma }^{2}}\right){\left({{ \mathcal B }}^{2}+W\right)}^{2}-\displaystyle \frac{{\left({Q}_{\mu }{{ \mathcal B }}^{\mu }\right)}^{2}\left({{ \mathcal B }}^{2}+2W\right)}{{W}^{2}}\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{Q}_{\mu }{n}^{\mu }=-\displaystyle \frac{{{ \mathcal B }}^{2}}{2}\left(2-\displaystyle \frac{1}{{\gamma }^{2}}\right)+\displaystyle \frac{{\left({Q}_{\mu }{{ \mathcal B }}^{\mu }\right)}^{2}}{2{W}^{2}}-W\,+\,P(\rho ,{Y}_{e},T)\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&\epsilon =\epsilon (\rho ,{Y}_{e},T).\end{eqnarray} \tag{ 22 }$$

Using these equations, we perform Newton–Raphson iterations until we obtain sufficiently accurate values for the independent variables γ, T, and W. Here, ${Q}_{\mu }=-{n}_{\nu }{T}_{\mu }^{\nu }=\alpha {T}_{\mu }^{t}$, and W is related to the specific enthalpy through W = h ρ γ².

${\tilde{Q}}^{\mu }={{j}^{\mu }}_{\nu }{Q}^{\nu }$, j_{μ ν} = g_{μ ν} + n_μ n_ν , and P is the pressure interpolated from tables.

2.2.2. Backup Recovery 1: 2D Routine

2.2.3. Backup Recovery 2: 2D "Safe-guess" Routine

If there is non-convergence for this backup routine, we include an initial "safe guess" as described in Cerdá-Durán et al. (2008). We call this routine the "2D safe guess" routine. In this scenario, we use the upper limits of the EOS table to obtain the maximum thermodynamical quantities:

$$\begin{eqnarray}&&{\rho }_{\max }=D,\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{T}_{\max }={T}_{\max ,\mathrm{tables}},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{P}_{\max }=P({\rho }_{\max },{Y}_{e},{T}_{\max }).\end{eqnarray} \tag{ 25 }$$

Here, D is the density measured in the orthogonal frame:

$$\begin{eqnarray}&&D\equiv -\rho {n}_{\mu }{u}^{\mu }=\gamma \rho .\end{eqnarray} \tag{ 26 }$$

Then we can estimate the initial "safe guess" for the root-finding procedure:

$$\begin{eqnarray}&&{\gamma }_{\mathrm{guess}}={\gamma }_{\max }=50,\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{W}_{\mathrm{guess}}={Q}_{\mu }{n}^{\mu }+{P}_{\max }-\displaystyle \frac{{{ \mathcal B }}^{2}}{2}.\end{eqnarray} \tag{ 28 }$$

2.2.4. Backup Recovery 3: 2D Dog Leg Routine

2.2.5. Backup Recovery 4: "Palenzuela" Routine

If all else fails, we use the routine described in Palenzuela et al. (2015). This routine solves a 1D equation using the Brent method. In this routine, called "Palenzuela," the independent variable is a rescaled variable

$$\begin{eqnarray}&&{x}_{\mathrm{pal}}\equiv \displaystyle \frac{\rho h{\gamma }^{2}}{\rho \gamma }.\end{eqnarray} \tag{ 29 }$$

We use the auxiliary rescaled variables:

$$\begin{eqnarray}&&{q}_{\mathrm{pal}}\equiv \displaystyle \frac{-({Q}_{\mu }{n}^{\mu }+D)}{D},{r}_{\mathrm{pal}}\equiv \displaystyle \frac{{\tilde{Q}}^{2}}{{D}^{2}},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{s}_{\mathrm{pal}}\equiv \displaystyle \frac{{{ \mathcal B }}^{2}}{D},{t}_{\mathrm{pal}}\equiv \displaystyle \frac{{Q}_{\mu }{{ \mathcal B }}^{\mu }}{{D}^{3/2}}.\end{eqnarray} \tag{ 31 }$$

The independent variable should be bracketed between

$$\begin{eqnarray}&&1+{q}_{\mathrm{pal}}-{s}_{\mathrm{pal}}\gt {x}_{\mathrm{pal}}\gt 2\,+\,2{q}_{\mathrm{pal}}-{s}_{\mathrm{pal}}.\end{eqnarray} \tag{ 32 }$$

The method uses an initial guess for x_pal from the previous time step and gets approximate quantities. Using them, it updates x_pal and iterates again until convergence is reached. The method is the following (where approximate quantities will be denoted by a hat).

We obtain an approximate Lorentz factor ${\hat{\gamma }}^{-2}$:

$$\begin{eqnarray}&&{\hat{\gamma }}^{-2}=1-\displaystyle \frac{{x}_{\mathrm{pal}}^{2}{r}_{\mathrm{pal}}+(2{x}_{\mathrm{pal}}+{s}_{\mathrm{pal}}){t}_{\mathrm{pal}}^{2}}{{x}_{\mathrm{pal}}^{2}{\left({x}_{\mathrm{pal}}+{s}_{\mathrm{pal}}\right)}^{2}}.\end{eqnarray} \tag{ 33 }$$

With that, we can estimate

$$\begin{eqnarray}&&\hat{\rho }=\displaystyle \frac{D}{\hat{\gamma }}\end{eqnarray} \tag{ 34 }$$

and an approximate specific energy:

$$\begin{eqnarray}&&\hat{\epsilon }=\hat{\gamma }-1+\displaystyle \frac{{x}_{\mathrm{pal}}}{\hat{\gamma }}(1-{\hat{\gamma }}^{2})+\hat{\gamma }\left({q}_{\mathrm{pal}}-{s}_{\mathrm{pal}}+\displaystyle \frac{{t}_{\mathrm{pal}}^{2}}{2{x}_{\mathrm{pal}}^{2}}+\displaystyle \frac{{s}_{\mathrm{pal}}}{2{\hat{\gamma }}^{2}}\right).\end{eqnarray} \tag{ 35 }$$

A call to the EOS will give the pressure $\hat{P}(\hat{\rho },\hat{\epsilon },{Y}_{e})$, and with all those approximate quantities, we can solve for x_pal using the Brent method by solving

$$\begin{eqnarray}&&0=f({x}_{\mathrm{pal}})={x}_{\mathrm{pal}}-\hat{\gamma }\left(1+\hat{\epsilon }+\displaystyle \frac{\hat{P}}{\hat{\rho }}\right).\end{eqnarray} \tag{ 36 }$$

We repeat the estimation of all the hat quantities until the solution for x_pal converges.

In the following section, we describe how we implemented a leakage scheme that takes into account the heating/cooling due to neutrinos, as well as how their emission and absorption affect the electron fraction. This leakage scheme is suited to describe the contribution of neutrinos to the composition and energy.

2.3.1. Rates

The scheme calculates the absorption/emission rate as well as the energy-loss rates due to neutrinos. We use these rates in the source terms of Equations (2) and (5). The scheme uses energy-averaged quantities.

Like Ruffert et al. (1996), Galeazzi et al. (2013), Siegel & Metzger (2018), we consider the following neutrino reactions, each with their own absorption/emission rate (which has units of cm⁻³ s⁻¹), and the energy-loss rate due to neutrinos (with units of erg cm⁻³ s⁻¹):

• 1.  
  Charged β-process with ${{ \mathcal R }}_{{\nu }_{i}}^{\beta }$ and ${{ \mathcal Q }}_{{\nu }_{i}}^{\beta }$:

  $$\begin{eqnarray}&&{e}^{-}+p\to n+{\nu }_{e}\end{eqnarray} \tag{ 37 }$$

  $$\begin{eqnarray}&&{e}^{+}+n\to p+{\bar{\nu }}_{e}\end{eqnarray} \tag{ 38 }$$
• 2.  
  Plasmon decay with ${{ \mathcal R }}_{{\nu }_{i}}^{\gamma }$ and ${{ \mathcal Q }}_{{\nu }_{i}}^{\gamma }$:

  $$\begin{eqnarray}&&\gamma \to {\nu }_{e}+{\bar{\nu }}_{e}\end{eqnarray} \tag{ 39 }$$

  $$\begin{eqnarray}&&\gamma \to {\nu }_{x}+{\bar{\nu }}_{x},\end{eqnarray} \tag{ 40 }$$

  where x is the muon and tauon, and in this case, γ corresponds to a photon.
• 3.  
  Electron–positron pair annihilation with ${{ \mathcal R }}_{{\nu }_{i}}^{{ee}}$ and ${{ \mathcal Q }}_{{\nu }_{i}}^{{ee}}$

  $$\begin{eqnarray}&&{e}^{-}+{e}^{+}\to {\nu }_{e}+{\bar{\nu }}_{e}\end{eqnarray} \tag{ 41 }$$

  $$\begin{eqnarray}&&{e}^{-}+{e}^{+}\to {\nu }_{x}+{\bar{\nu }}_{x}.\end{eqnarray} \tag{ 42 }$$

Using the above reactions, we calculate the total number emission in the optically thin regime from species i as (Ruffert et al. 1996)

$$\begin{eqnarray}&&{{ \mathcal R }}_{{\nu }_{i}}={{ \mathcal R }}_{{\nu }_{i}}^{\beta }+{{ \mathcal R }}_{{\nu }_{i}}^{\gamma }+{{ \mathcal R }}_{{\nu }_{i}}^{{ee}},\end{eqnarray} \tag{ 43 }$$

and the total energy-loss rate in the optically thin regime is

$$\begin{eqnarray}&&{{ \mathcal Q }}_{{\nu }_{i}}={{ \mathcal Q }}_{{\nu }_{i}}^{\beta }+{{ \mathcal Q }}_{{\nu }_{i}}^{\gamma }+{{ \mathcal Q }}_{{\nu }_{i}}^{{ee}},\end{eqnarray} \tag{ 44 }$$

where "i" denotes the different neutrino/antineutrino flavors: electron, or muon and tauon.

The total emission/absorption rates and the energy-loss rates are given by an interpolation between the diffusive optically thick regime and the transparent optically thin regime (Ruffert et al. 1996):

$$\begin{eqnarray}&&{{ \mathcal R }}_{{\nu }_{i}}^{\mathrm{eff}}={{ \mathcal R }}_{{\nu }_{i}}{\left(1+\displaystyle \frac{{t}_{\mathrm{diff}}}{{t}_{\mathrm{emission},{ \mathcal R }}}\right)}^{-1}\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}&&{{ \mathcal Q }}_{{\nu }_{i}}^{\mathrm{eff}}={{ \mathcal Q }}_{{\nu }_{i}}{\left(1+\displaystyle \frac{{t}_{\mathrm{diff}}}{{t}_{\mathrm{emission},{ \mathcal Q }}}\right)}^{-1}.\end{eqnarray} \tag{ 46 }$$

Here, the diffusion timescale is given by

$$\begin{eqnarray}&&{t}_{\mathrm{diff}}=\displaystyle \frac{{D}_{\mathrm{diff}}{\tau }^{2}}{c{\kappa }_{{\nu }_{i}}},\end{eqnarray} \tag{ 47 }$$

where D_diff = 6 (Rosswog & Liebendörfer 2003; O'Connor & Ott 2010; Siegel & Metzger 2018), τ is the optical depth, and ${\kappa }_{{\nu }_{i}}$ is the energy-averaged opacity (in units of cm⁻¹) of ν_i . The absorption/emission and energy loss timescales are ${t}_{\mathrm{emission},{ \mathcal R }}\,={{ \mathcal R }}_{{\nu }_{i}}/{n}_{{\nu }_{i}}$, with ${n}_{{\nu }_{i}}$ being the neutrino number density (at chemical equilibrium), and ${t}_{\mathrm{emission},{ \mathcal Q }}={{ \mathcal Q }}_{{\nu }_{i}}/{\varepsilon }_{{\nu }_{i}}$, with ${\varepsilon }_{{\nu }_{i}}$ being the neutrino energy density. In the optically thick regime, the neutrino-loss rate is less than the diffusion time, which results in ${{ \mathcal R }}_{{\nu }_{i}}^{\mathrm{eff}}={n}_{{\nu }_{i}}/{t}_{\mathrm{diff}}$ and ${{ \mathcal Q }}_{{\nu }_{i}}^{\mathrm{eff}}={\varepsilon }_{{\nu }_{i}}/{t}_{\mathrm{diff}}$, whereas in the optically thin regime, we recover the rates from Equations (43) and (44). The rates for the muon and tauon neutrinos/antineutrinos estimated in Ruffert et al. (1996) take into account all four of those species. We also note that several quantities, including the chemical potentials, are obtained from EOS table interpolation.

2.3.2. Optical Depth

The transition between the two regimes will be set by the optical depth ${\tau }_{{\nu }_{i}}$, which is also needed to obtain the diffusion timescale. In order to get the optical depth, we consider the following reactions as the source of neutrino opacity:

$$\begin{eqnarray}&&{\nu }_{e}+n\to p\,+\,{e}^{-}\end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray}&&{\bar{\nu }}_{e}+p\to n\,+\,{e}^{+}\ \end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&{\nu }_{i}+p\to {\nu }_{i}+p\end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}&&{\nu }_{i}+n\to {\nu }_{i}+n.\end{eqnarray} \tag{ 51 }$$

The opacities are obtained from Ruffert et al. (1996). Electron scattering is neglected.

The usual global approach to calculate the optical depth of a point in the flow would be to integrate the opacity over all directions and determine the path of minimal absorption. This approach assumes that the neutrino will follow a straight path. However, we follow Neilsen et al. (2014) and Siegel & Metzger (2018), where a local, iterative approach is used instead of a global calculation, and where crooked minimal paths are acceptable. The optical depth is calculated by obtaining the shortest path of the neutrino out of the star using its neighbors. For the first time step, we begin by initializing the optical depth grid to zero. Next, we perform the first iteration, where we estimate the optical depth at each cell as the minimum of the optical depth of its neighbor (${\tau }_{{\nu }_{i},\mathrm{neighbor}}$) plus the optical depth needed for the neutrino to reach that neighbor (${\bar{\kappa }}_{{\nu }_{i}}{({g}_{{kj}}{{dx}}^{k}{{dx}}^{j})}^{1/2}$):

$$\begin{eqnarray}&&{\tau }_{{\nu }_{i}}=\min \left({\tau }_{{\nu }_{i},\mathrm{neighbor}}+{\bar{\kappa }}_{{\nu }_{i}}{({g}_{{kj}}{{dx}}^{k}{{dx}}^{j})}^{1/2}\right),\end{eqnarray} \tag{ 52 }$$

where ${\tau }_{{\nu }_{i},\mathrm{neighbor}}$ is the optical depth of the neighboring cell, ${\bar{\kappa }}_{{\nu }_{i}}$ is the average opacity between the cell and its neighbor, and ${({g}_{{kj}}{{dx}}^{k}{{dx}}^{j})}^{1/2}$ is the distance to the neighboring cell calculated by taking the average value of g_kj between the local and neighboring cells. We minimize over all neighbors.

This essentially traces the path of least resistance of the neutrino to a neighbor. We update the entire grid and perform the next iteration, where again, we minimize over all the adjacent neighbors. The next iteration will show the path to the neighbor two cells away. As we do more iterations, we trace the path of least resistance that the neutrinos will take out of the star. This will lead us to the final optical depth. During the first time step, we initialize the optical depth by by performing $20{N}_{\max }$ iterations, where ${N}_{\max }$ is the maximum number of cells in each direction, independent of resolution. This is done to trace a path to the edge of the domain initially. After the initial calculation, which has a fixed number of iterations, we continue to do iterations to obtain the final optical depth; however, we impose a convergence criterion in order to minimize the number of iterations. In order to converge, we set conditions on the difference between iteration k − 1 and k:

$$\begin{eqnarray}&&{R}_{\mathrm{change},\tau }(k)\equiv \displaystyle \frac{\left|\sum {\tau }_{k-1}-\sum {\tau }_{k}\right|}{\sum {\tau }_{k-1}}\lt {\epsilon }_{1}\end{eqnarray} \tag{ 53 }$$

or

$$\begin{eqnarray}&&\displaystyle \frac{| {R}_{\mathrm{change},\tau }(k-1)-{R}_{\mathrm{change},\tau }(k)| }{{R}_{\mathrm{change},\tau }(k-1)}\lt {\epsilon }_{2},\end{eqnarray} \tag{ 54 }$$

where ∑τ_k is the sum of all the optical depths in the grid at iteration k, and ₁ and ₂ are parameters that we choose to be ₁ = 10⁻⁴ and ₂ = 10⁻³, respectively. Only a few iterations are needed for convergence after the initial guess.

In order to validate the routines that transform the conserved variables into primitive variables with tabulated EOS, we created primitive variables out of a grid of density and temperature values within the EOS table. The magnetic field was set randomly to be either aligned or antialigned with the velocity vector. The magnitude of the magnetic field was set to be such that ${b}^{2}/2=\left(\tfrac{{P}_{\mathrm{mag}}}{{P}_{\mathrm{gas}}}\right){P}_{\mathrm{gas}}$, where (P_mag/P_gas) is set as a parameter, P_mag is the magnetic pressure, and P_gas is the gas pressure. We then obtained a set of conserved variables based on these primitives. The true primitives were then varied by randomly adding or subtracting a 5% perturbation to each primitive. This test is based on Siegel et al. (2018).

We then used these primitives as initial guesses for the various routines that transform the conserved variables to primitive variables and compared the resultant solution to the original.

We show the error we obtained for all primitive variables in Figure 1. It can be seen that the recovery error is low. Additionally, the figure shows that the 3D method is less robust, but more accurate, which is the reason it is set as the primary routine. The different 2D methods and the "Palenzuela" routine are more robust, but less accurate (and slower) than the 3D method, so they serve better as backup routines.

Zoom In Zoom Out Reset image size

To test the EOS implementation, we simulated a nonmagnetized torus that is in hydrostatic equilibrium with no leakage scheme, following Fishbone & Moncrief (1976).

Figure 2 shows the 3D hydrodynamical evolution of a torus constructed to be in hydrostatic equilibrium with a tabulated EOS without neutrino cooling. There are perturbations, in particular near the BH, due to accretion onto the BH, but the density is low in those regions. As can be seen from the figure, the torus remains in hydrostatic equilibrium throughout the simulation.

Zoom In Zoom Out Reset image size

3.2.1. Initial Conditions inside the Torus

The specific enthalpy inside the torus is implemented via Equation (3.6) of Fishbone & Moncrief (1976), but adding $\mathrm{ln}{h}_{\min }$ to the integration constant (see Section 3.2.2). By construction, the torus is in hydrostatic equilibrium with the ambient atmosphere. We also set the torus to be isentropic and have uniform electron fraction. Given a specific entropy s_disk, a specific enthalpy given by Fishbone & Moncrief (1976), and an electron fraction Y_e,disk), the temperature and density of the disk are found by solving the following equations:

$$\begin{eqnarray}&&{s}_{\mathrm{disk}}=s(\rho ,T,{Y}_{e})\end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray}&&h=h(\rho ,T,{Y}_{e}),\end{eqnarray} \tag{ 56 }$$

where s is the specific entropy.

3.2.2. Atmosphere

In the classical torus, the boundaries of the torus are defined where h = 1. In the tabulated EOS, though, negative internal energy densities are allowed because the internal energy per nucleon is measured relative to the free neutron rest mass energy. In this case, the minimum specific enthalpy is not restricted to 1, but rather it can be $1\gt {h}_{\min }\gt 0$, where ${h}_{\min }$ is the specific enthalpy from the table, given the atmospheric density and the disk's electron fraction and specific entropy. Thus, we set the torus boundary to be where $h={h}_{\min }$. For the background atmosphere, we set the minimum atmospheric density ρ_atm as a parameter. We then find the minimum specific enthalpy by doing a table inversion and finding ${h}_{\min }\,=h({\rho }_{\mathrm{atm}},{s}_{\mathrm{disk}},{Y}_{e,\mathrm{disk}})$. We also find the atmospheric temperature by doing a table inversion T_atm = T(ρ_atm, s_disk, Y_e,disk).

The density in the background is set to

$$\begin{eqnarray}&&\max \left({\rho }_{\mathrm{atm}},\displaystyle \frac{{\rho }_{0}}{{r}^{2}}\right),\end{eqnarray} \tag{ 57 }$$

where we set ρ₀ as a parameter as well. The background atmosphere temperature is set to

$$\begin{eqnarray}&&\max \left({T}_{\mathrm{atm}},\displaystyle \frac{{T}_{0}}{r}\right),\end{eqnarray} \tag{ 58 }$$

where T₀ is a parameter. The power-law dependence is set to provide the background atmosphere with more pressure support so that it does not rapidly accrete onto the BH. This ultimately helps with robustness near the BH, as the low-density and low-temperature zones with high velocity are where the conserved to primitive routines tend to fail. We note that, in the region where there is a power-law dependence, the specific enthalpy is not a constant, whereas once the background atmosphere is set to be constant, everything is thermodynamically consistent because it was constructed with the tabulated EOS tables. We set the electron fraction of the atmosphere to a constant value found by assuming β-equilibrium (where the neutrino chemical potential is zero) at T_atm and ρ_atm.

The units are normalized so that the maximum density in the torus is set to ${\rho }_{\max }=1$ in code units, which in this case corresponds to ${\rho }_{\max }=5.4\times {10}^{8}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ in cgs units. In the simulation we performed, the torus has a constant electron fraction of Y_e = 0.1 and a specific entropy of 10k_B/baryon, where K_B is Boltzmann's constant. The background atmosphere is characterized by ρ_atm = 6000 g cm⁻³, ρ₀ = 3 × 10⁵ g cm⁻³, and T₀ = 0.4 MeV. We used the SLy4 table with NSE from Schneider et al. (2017), and with that table the minimum specific enthalpy for our parameters is set to ${h}_{\min }=0.9974$ (in code units), and T_atm = 0.0053 MeV. The electron fraction in the atmosphere, given by β-equilibrium, is set to Y_e,atm = 0.45. The boundary conditions are outflow in the outer radial boundary, reflective in the angular coordinate θ, and periodic in the angular coordinate ϕ. The metric is Kerr–Schild in spherical coordinates for a non-spinning BH.

Following Miller et al. (2019a), we tested the leakage scheme in an optically thin regime by considering an isotropic gas of constant density and temperature such that the gas is optically thin to neutrinos. We conducted separate tests of both reactions in the charged β-process, wherein we included only either the neutrinos or the antineutrinos.

In this case, the GRMHD equations reduce to

$$\begin{eqnarray}&&{\partial }_{t}{{T}^{t}}_{t}={ \mathcal Q }.,\end{eqnarray} \tag{ 59 }$$

$$\begin{eqnarray}&&{\partial }_{t}{Y}_{e}={ \mathcal R }/\rho ,\end{eqnarray} \tag{ 60 }$$

where ${ \mathcal R }$ and ${ \mathcal Q }$ are the respective emission/absorption and energy-loss rates due to neutrinos or antineutrinos of the reactions in the β-process. The rates need to be calculated semi-analytically, since they depend on interpolated quantities, such as the degeneracy parameters. We can then solve the equations semi-analytically with a set of initial conditions and compare to simulations. We chose the initial density and temperature such that the medium is optically thin to neutrinos and antineutrinos.

For the initial conditions, we used an initial density of 617714 g cm⁻³ and temperature of 1 MeV, chosen so that the medium is optically thin to neutrinos and antineutrinos. We used Y_e,0 = 0.5 and Y_e,0 = 0.005 for the electron neutrino and antineutrino tests, respectively. We used 2 × 2 × 1 cells in each direction using a Cartesian grid with Minkowski metric.

In Figure 3, we show the comparison between the semi-analytical solution and the simulation for the β-process both for neutrinos and antineutrinos. We compare the change in the electron fraction due to the absorption/emission rate and the change in temperature due to the heating/cooling rate. As can be seen from the figure, HARM3D+NUC is able to recreate the semi-analytical solution.

Zoom In Zoom Out Reset image size

4.2.1. Constant Density Circular Disk

In order to test the optical depth calculation, we simulated a circular disk with uniform density and temperature embedded in an optically thin medium of constant density and temperature. The advantage of this scenario is that we can calculate the opacity inside the circle and then calculate the optical depth analytically. This way, we can compare to the simulation. The simulations were performed in 2D, and the domain is 2r_g, where r_g = GM/c² is the gravitational radius. We used a Minkowski metric with spherical coordinates. There are outflow conditions on the radial boundaries. The optical depth in the outer radial boundary was set to zero so that the neutrinos and antineutrinos could escape the domain. We simulated an optically thick circular disk that has a constant density of 9.8 × 10¹³ g cm⁻³, an electron fraction of 0.1, and a temperature of 8 MeV, embedded in an optically thin medium, with a density of 6 × 10⁷ g cm⁻³, an electron fraction of 0.5, and a temperature of 0.01 MeV. Figure 4 shows the optical depth for both the electron neutrino and antineutrino for different resolutions. As can be seen from the figure, the initial guess for the optical depth is accurate and the convergence to the solution does not change with resolution. At smaller optical depths, the optical depth is slightly overestimated at lower resolutions, but as the optical depth increases, the solution does not depend noticeably on resolution.

Zoom In Zoom Out Reset image size

4.2.2. Stripes

We can also test the optical depth algorithm by simulating stripes of high-density material with low-density material in between. In this scenario, it is expected that a neutrino created in the region with high optical depth material will travel to the region with low optical depth and stream freely from the surface. For the simulation, we used 4096 × 96 × 1 cells. The simulations were performed in 2D, and the domain is 1r_g large in radial extent, where r_g = GM/c² is the gravitational radius. We used a Minkowski metric with spherical coordinates and outflow conditions at the radial boundaries. The optical depth at the outer radial boundary was set to zero so that the neutrinos and antineutrinos could escape the domain. We simulated three stripes of material with high optical depth: ρ = 9.8 × 10¹³ g cm⁻³, Y_e = 0.1, and T = 8 MeV. In between the stripes, the optically thin gas was initialized to ρ = 6 × 10⁷ g cm⁻³, Y_e = 0.5, and T = 0.01 MeV. The high-opacity stripes start at r = 0r_g and have a width of r = 0.1r_g. The next stripes are located in r = 0.2r_g and r = 0.4r_g.

We show the results from this setup in Figure 5, where we compare the results from the simulation with the analytical estimate (length units are in r_g):

$$\begin{eqnarray}{\tau }_{\mathrm{analytical}}=\left\{\begin{array}{ll}{\displaystyle \int }_{0}^{0.2}\kappa {dr} & r\leqslant 0.2\\ {\displaystyle \int }_{0.2}^{0.25}\kappa {dr} & 0.2\leqslant r\leqslant 0.25\\ {\displaystyle \int }_{0.25}^{0.35}\kappa {dr} & 0.25\leqslant r\leqslant 0.35\\ {\displaystyle \int }_{0.4}^{0.45}\kappa {dr} & 0.4\leqslant r\leqslant 0.45\\ {\displaystyle \int }_{0.45}^{0.55}\kappa {dr} & 0.45\leqslant r\leqslant 0.55\end{array}\right..\end{eqnarray} \tag{ 61 }$$

Zoom In Zoom Out Reset image size

The initial conditions inside the torus follow a setup similar to that of Section 3.2.1, but with the addition of a poloidal magnetic field. In order to start with a magnetic field devoid of magnetic monopoles, we first set the vector potential to a prescribed distribution and calculate its curl using a finite difference operator compatible with our constrained transport method (see Zilhão & Noble 2014 for further details). Our poloidal magnetic field distribution results from a vector potential with only one nonzero component,

$$\begin{eqnarray}&&{A}_{\phi }=\max \left(\overline{\rho }/{\rho }_{\max }-{\rho }_{0,\mathrm{mag}},0\right),\end{eqnarray} \tag{ 62 }$$

where $\overline{\rho }$ is the average density at that position and ${\rho }_{\max }\,=1.66\times {10}^{11}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ is the maximum density of the torus. We set ρ_0,mag = 0.2 in code units, which corresponds to ρ_0,mag = 3.33 × 10¹⁰ g cm⁻³. We then build the magnetic field with the vector potential and normalize its magnitude such that the ratio of the integrated gas pressure to integrated magnetic pressure is 100. Inside the disk, the matter is set to be neutron-rich, Y_e = 0.1. The treatment of the atmosphere is the same as in Section 3.2.2, except the density scales as r^−3/2. In the atmosphere, the electron fraction is set to its value in β-equilibrium, where the chemical potential of the neutrinos is set to zero. We show the parameters used in Table 1.

The simulations were performed in 3D on a grid designed to focus more cells about the equator and toward the black hole horizon. We use the same grid as defined in Noble et al. (2010) but with different parameters. The azimuthal grid spacing is uniform. The logarithmic radial grid is such that Δr/r is fixed and the ith cell center is located at:

$$\begin{eqnarray}&&{r}_{i}={r}_{\min }\exp \left[\left(i+1/2\right){\mathrm{log}}_{10}({r}_{\max }/{r}_{\min })/{N}_{r}\right],\end{eqnarray} \tag{ 63 }$$

with ${r}_{\min }=1.303{r}_{{\rm{g}}}$, ${r}_{\max }=2000{r}_{{\rm{g}}}$, and $i\in \left[0,{N}_{r}-1\right]$. The θ grid uses a high-order polynomial function to provide a nearly uniform grid spacing spacing near the equator:

$$\begin{eqnarray}&&{\theta }_{j}=\displaystyle \frac{\pi }{2}\left[1+\left(1-\xi \right)\left(2{x}_{j}^{(2)}-1\right)+\left(\xi -\displaystyle \frac{2{\theta }_{c}}{\pi }\right){\left(2{x}_{j}^{(2)}-1\right)}^{n}\right],\end{eqnarray} \tag{ 64 }$$

where ξ is a parameter controlling the severity of the focusing, n is the order of polynomial used in the transformation, θ_c is the opening angle of the polar regions we excise, ${x}_{j}^{(2)}\,\equiv \left(j+1/2\right)/{N}_{\theta }$, and $j\in \left[0,{N}_{\theta }-1\right]$. In our run, we used θ_c = π10⁻¹⁴, ξ = 0.65, and n = 7. The number of cells per dimension used was N_r × N_θ × N_ϕ = 1024 × 160 × 256.

We performed scaling tests for this run for three different numbers of processors: 5120, 2560, and 1280 processors. For this setup, the numbers of time steps in the code per second per processor were 0.000723, 0.000781, and 0.000868, respectively. The difference between 5120 and 1280 processors is around 17%. If we exclude the neutrino-leakage scheme and include only a tabulated EOS, for 2560 processors, the number of steps per second per processor is 0.001328, which makes the leakage 58% slower than only considering the tabulated EOS.

In order to confirm that we are adequately resolving magnetic turbulence, we display in Figure 6 the number of grid cells per wavelength of the fastest growing mode of the magnetorotational instability (MRI), defined as (Noble et al. 2010; Hawley et al. 2011, 2013; Sorathia et al. 2012)

$$\begin{eqnarray}&&{Q}_{\mathrm{mri},{\rm{x}}}=\displaystyle \frac{{\lambda }_{x,\mathrm{mri}}}{{{\rm{\Delta }}}_{x}},\end{eqnarray} \tag{ 65 }$$

where x = θ, ϕ, Δ_x is the cell size, and the wavelength of the fastest MRI growing mode is

$$\begin{eqnarray}&&{\lambda }_{x,\mathrm{mri}}=\displaystyle \frac{2\pi }{{\rm{\Omega }}}\displaystyle \frac{| {b}^{x}| }{\sqrt{\rho h+{b}^{2}}}.\end{eqnarray} \tag{ 66 }$$

Zoom In Zoom Out Reset image size

As can be seen in Figure 7, our grid satisfies the criterion of Sano et al. (2004) everywhere except for later times within r < 50 km. While our results fail to meet criteria for asymptotic MRI convergence set forth in Hawley et al. (2011), our disk does satisfy Q_mri,3 > 10 everywhere, and Q_mri,2 > 6 for r ≳ 50 km for most of the run. One reason why our simulation may not reach larger Q_mri values is because we used the same random perturbations across all MPI processes in the initial conditions. Because we used 16 subdomains in the azimuthal dimension, this means that the simulation is nearly periodic over Δϕ = π/8, and the azimuthal modes with m < 8 start off significantly weaker because they are seeded with perturbations at only the round-off error level.

Zoom In Zoom Out Reset image size

Magnetic stresses will transport angular momentum in the disk, heating the gas, which will produce a high-velocity outflow (Fernández & Metzger 2013a; Siegel & Metzger 2018). This outflow will be affected by the addition of neutrinos formed through weak reactions. In the midplane, neutrinos will carry significant amounts of energy, which will cool and make the disk geometrically thinner. Another outflow is also expected to occur in the outer regions of the disk, due to the release in nuclear binding energy when there is recombination of free nucleons into α-particles, which produces enthalpy and unbinds material (Lee et al. 2009; Fernández & Metzger 2013a). In this subsection, we show the impact of both the emission of neutrinos and the recombination of free nucleons.

In Figures 8 and 11, we display the outflows that result from our simulations of a neutrino-cooled magnetized disk at 114 ms. In Figures 9 and 10, we plot the electron fraction and the density, respectively, at t = 114 ms. Movies can be found here. ²⁴

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Neutrino cooling is expected to happen on the diffusion timescale, which is on the order of milliseconds—much shorter than our evolution timescale. The inner regions of the disk are very neutron-rich, confirming the self-regulating phase found in Siegel & Metzger (2017, 2018). In this phase, there is a balance between the neutrino cooling and the heating driven by MHD that self-regulates the electron degeneracy parameter, and the final state is a neutron-rich disk (Siegel & Metzger 2017, 2018). We note that, although this new code does not include neutrino absorption in the ejecta, absorption will modify the electron fraction in the outflow (Just et al. 2021).

Zoom In Zoom Out Reset image size

In Figure 12, we plot the geometrical thickness of the disk, or H/r. We estimated this thickness using the scale height H following Noble et al. (2012),

$$\begin{eqnarray}&&H=\displaystyle \frac{\langle \rho \sqrt{{g}_{\theta \theta }}\left|\theta -\pi /2\right|\rangle }{\langle \rho \rangle },\end{eqnarray} \tag{ 67 }$$

where 〈X〉 is the average of the quantity X over a spherical shell,

$$\begin{eqnarray}&&\langle X\rangle =\displaystyle \frac{\int X\sqrt{-g}d\theta d\phi }{\int \sqrt{-g}d\theta d\phi }.\end{eqnarray} \tag{ 68 }$$

In the top panel of Figure 13, we show the mass accretion rate through the innermost stable circular orbit as a function of time, and in the bottom panel, we show the accretion rate as a function of radius. The outflow can be clearly seen as a negative mass accretion rate at larger radii, as well as a settling of the mass accretion rate as time passes.

In the deepest regions of the disk, the heating due to MHD turbulence helps create neutrinos/antineutrinos, which escape, remove energy, and geometrically thin the disk. Recombination of free nucleons into α-particles releases binding energy, effectively increasing the enthalpy and unbinds material. The effect of recombination is less severe than the geometric thinning due to neutrino/antineutrino losses. This transition can be seen at around 150 km.

Zoom In Zoom Out Reset image size

We may obtain the amount of energy radiated by each species of neutrino and antineutrino as was done in Siegel & Metzger (2018):

$$\begin{eqnarray}&&{L}_{{\nu }_{i}}=\int \alpha \gamma {{ \mathcal Q }}_{{\nu }_{i}}^{\mathrm{eff}}\sqrt{-g}{d}^{3}x.\end{eqnarray} \tag{ 69 }$$

In Figure 14, we show the luminosity for each species. It can be seen that the electron neutrino (and antineutrino) dominate the emission over all of the other species of neutrino. The luminosity roughly follows the mass accretion rate as seen in Figure 13, as heating from the magnetic stresses ignites the creation of neutrinos/antineutrinos. This suggests the radiative efficiency of neutrino/antineutrinos emission remains relatively steady.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our initial conditions are similar (although not identical) to the initial conditions in Siegel & Metzger (2018). They performed 3D simulations of a post-merger accretion disk with a relatively higher specific entropy and lower spin than this simulation. They used Cartesian coordinates, a Helmholtz EOS for relatively low densities, and a neutrino-leakage scheme. They also evolved the disk for longer times (380 ms). Even though we use a different EOS (Sly4), the disk thickness is qualitatively similar. At the inner regions of the disk, neutrino cooling dominates, whereas at outer regions (at radius higher than around 100 km), recombination is responsible for making the disk geometrically thicker. The neutrino/antineutrino luminosities are comparable; Siegel & Metzger (2018) report a higher luminosity, but that could be attributed to the difference in the initial disk specific entropy.

As the outflow expands, it will cool and heavy elements will be created via the r-process. We will explore this nucleosynthesis in a future paper.