Three-dimensional GRMHD Simulations of Neutrino-cooled Accretion Disks from Neutron Star Mergers

The equations of ideal GRMHD with weak interactions include energy and momentum conservation, baryon number conservation, lepton number conservation, and Maxwell's equations,

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }{T}^{\mu \nu }={{Qu}}^{\nu },\end{eqnarray} \tag{ 1 }$$

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

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }({n}_{{\rm{e}}}{u}^{\mu })=R,\end{eqnarray} \tag{ 3 }$$

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

where

$$\begin{eqnarray}&&{T}^{\mu \nu }=(\rho h+{b}^{2}){u}^{\mu }{u}^{\nu }+\left(p+\displaystyle \frac{{b}^{2}}{2}\right){g}^{\mu \nu }-{b}^{\mu }{b}^{\nu }\end{eqnarray} \tag{ 5 }$$

is the energy–momentum tensor, u^μ is the four-velocity, ${n}_{{\rm{b}}}$ is the baryon number density, ${n}_{{\rm{e}}}$ is the electron number density, and ${F}^{* \mu \nu }$ is the dual of the Faraday electromagnetic tensor. Furthermore, p is the pressure; $h=1+\epsilon +p/\rho$ denotes the specific enthalpy, with being the specific internal energy; ${b}^{\mu }\equiv {(4\pi )}^{-1/2}{F}^{* \mu \nu }{u}_{\nu }$ is the magnetic field vector in the frame comoving with the fluid; b² ≡ b^μb_μ; and g_μν is the spacetime metric.⁶ We assume that the thermodynamic properties of matter can be described by a finite-temperature, composition-dependent (three-parameter) EOS formulated as a function of density $\rho ={n}_{{\rm{b}}}{m}_{{\rm{b}}}$, where m_b denotes the baryon mass; temperature T; and electron fraction ${Y}_{{\rm{e}}}={n}_{{\rm{e}}}{n}_{{\rm{b}}}^{-1}$. The evolution of Y_e is described by Equation (3). The source terms Qu^ν and R on the right-hand side of Equations (1) and (3) account for the evolution of Y_e due to weak interactions, which create neutrinos and antineutrinos that carry away energy and momentum from the system.

For numerical evolution, Equations (1)–(4) can essentially be transformed into a set of conservation equations in flat space by adopting a 3 + 1 split of spacetime into nonintersecting space-like hypersurfaces of constant coordinate time t (Lichnerowicz 1944; Arnowitt et al. 2008), in which case, the line element can be written as

$$\begin{eqnarray}&&{{ds}}^{2}=-{\alpha }^{2}{{dt}}^{2}+{\gamma }_{{ij}}({{dx}}^{i}+{\beta }^{i}{dt})({{dx}}^{j}+{\beta }^{j}{dt}),\end{eqnarray} \tag{ 6 }$$

where α denotes the lapse function, βⁱ is the shift vector, and γ_ij is the metric induced on every spatial hypersurface. The hypersurfaces are characterized by the time-like unit normal n^μ = (α⁻¹, −α⁻¹βⁱ) (${n}_{\mu }=(-\alpha ,0,0,0)$), which also defines the Eulerian observer, i.e., the observer moving through spacetime with four-velocity n^μ perpendicular to the hypersurfaces. Equations (1)–(4) can then be written as

$$\begin{eqnarray}&&{\partial }_{t}(\sqrt{\gamma }{\boldsymbol{q}})+{\partial }_{i}[\alpha \sqrt{\gamma }{{\boldsymbol{f}}}^{(i)}({\boldsymbol{p}},{\boldsymbol{q}})]=\alpha \sqrt{\gamma }{\boldsymbol{s}}({\boldsymbol{p}}),\end{eqnarray} \tag{ 7 }$$

where γ is the determinant of the spatial metric γ_ij and

$$\begin{eqnarray}&&{\boldsymbol{q}}\equiv [D,{S}_{i},\tau ,{B}^{i},{{DY}}_{{\rm{e}}}]\end{eqnarray} \tag{ 8 }$$

denotes the vector of conserved variables. The latter is composed of the conserved density, momenta, and energy, defined as

$$\begin{eqnarray}&&D\equiv \rho W,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{S}_{i}\equiv -{n}_{\mu }{T}_{i}^{\mu }=\alpha {T}_{i}^{0}=(\rho h+{b}^{2}){W}^{2}{v}_{i}-\alpha {b}^{0}{b}_{i},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\tau \equiv {n}_{\mu }{n}_{\nu }{T}^{\mu \nu }-D\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&=\,(\rho h+{b}^{2}){W}^{2}-\left(p+\displaystyle \frac{{b}^{2}}{2}\right)-{\alpha }^{2}{({b}^{0})}^{2}-D,\end{eqnarray} \tag{ 12 }$$

respectively, the three-vector components of the magnetic field ${B}^{\mu }\equiv {(4\pi )}^{-1/2}{F}^{* \mu \nu }{n}_{\nu }$ as measured by the Eulerian observer, as well as the conserved electron fraction DY_e. The Eulerian three-velocity is defined by

$$\begin{eqnarray}&&{v}^{i}\equiv \displaystyle \frac{{\gamma }_{\mu }^{i}{u}^{\mu }}{-{u}^{\mu }{n}_{\mu }}=\displaystyle \frac{{u}^{i}}{W}+\displaystyle \frac{{\beta }^{i}}{\alpha },{v}_{i}=\displaystyle \frac{{\gamma }_{i\mu }{u}^{\mu }}{-{u}^{\mu }{n}_{\mu }}=\displaystyle \frac{{u}_{i}}{W},\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}&&W\equiv -{u}^{\mu }{n}_{\mu }=\alpha {u}^{0}=\displaystyle \frac{1}{\sqrt{1-{v}^{2}}}\end{eqnarray} \tag{ 14 }$$

denotes the relative Lorentz factor between u^μ and n^μ, with ${v}^{2}\equiv {\gamma }_{{ij}}{v}^{i}{v}^{j}$. For completeness, the comoving and Eulerian magnetic field components are related by

$$\begin{eqnarray}&&{b}^{i}=\displaystyle \frac{{B}^{i}}{W}+{b}^{0}(\alpha {v}^{i}-{\beta }^{i}),{b}_{i}=\displaystyle \frac{{B}_{i}}{W}+\alpha {b}^{0}{v}_{i}\end{eqnarray} \tag{ 15 }$$

and

$$\begin{eqnarray}&&{b}^{0}=\displaystyle \frac{W}{\alpha }{B}^{i}{v}_{i},{b}^{2}={b}^{\mu }{b}_{\mu }=\displaystyle \frac{{B}^{2}+{(\alpha {b}^{0})}^{2}}{{W}^{2}},\end{eqnarray} \tag{ 16 }$$

where ${B}^{2}\equiv {B}^{i}{B}_{i}$. Furthermore,

$$\begin{eqnarray}&&{\boldsymbol{p}}\equiv [\rho ,{v}^{i},\epsilon ,{B}^{i},{Y}_{{\rm{e}}}]\end{eqnarray} \tag{ 17 }$$

summarizes the set of primitive variables. The fluxes are given by

$$\begin{eqnarray}&&{{\boldsymbol{f}}}^{(i)}({\boldsymbol{p}},{\boldsymbol{q}})\equiv \left[\begin{array}{c}D{\widetilde{v}}^{i}\\ {S}_{j}{\widetilde{v}}^{i}+\left(p+\displaystyle \frac{{b}^{2}}{2}\right){\delta }_{j}^{i}-\displaystyle \frac{{B}^{i}}{W}{b}_{j}\\ \tau {\widetilde{v}}^{i}+\left(p+\displaystyle \frac{{b}^{2}}{2}\right){v}^{i}-\alpha {b}^{0}\displaystyle \frac{{B}^{i}}{W}\\ {\widetilde{v}}^{i}{B}^{k}-{\widetilde{v}}^{k}{B}^{i}\\ {DY}_{{\rm{e}}}{\widetilde{v}}^{i}\end{array}\right]\end{eqnarray} \tag{ 18 }$$

and the sources by

$$\begin{eqnarray}&&{\boldsymbol{s}}({\boldsymbol{p}})\equiv \left[\begin{array}{c}0\\ {T}^{\mu \nu }({\partial }_{\mu }{g}_{j\nu }-{{\rm{\Gamma }}}_{\nu \mu }^{\delta }{g}_{\delta j})+{{WQv}}_{j}\\ \alpha ({T}^{0\mu }{\partial }_{\mu }\mathrm{ln}\alpha -{T}^{\mu \nu }{{\rm{\Gamma }}}_{\mu \nu }^{0})+{WQ}\\ {0}^{k}\\ {{Rm}}_{{\rm{b}}}\end{array}\right],\end{eqnarray} \tag{ 19 }$$

where ${\tilde{v}}^{i}\equiv {v}^{i}-{\beta }^{i}{\alpha }^{-1}$, and ${{\rm{\Gamma }}}_{\beta \gamma }^{\alpha }$ are the Christoffel symbols constructed from g_μν.

Weak interactions and neutrino transport determine the source terms on the right-hand side of Equations (1) and (3) and the terms apart from the geometrical source terms in Equation (7) (cf. Equation (19)). For the present simulations, we employ an energy-averaged (gray) leakage scheme, which we have newly implemented into GRHydro. Such leakage schemes are widely used in both core-collapse supernova and compact-binary merger simulations (e.g., van Riper & Lattimer 1981; Ruffert et al. 1996; Rosswog & Liebendörfer 2003; Sekiguchi et al. 2011; Ott et al. 2013; Perego et al. 2016; Radice et al. 2016). Our implementation closely follows the one by Radice et al. (2016), which is based on Galeazzi et al. (2013), which, in turn, builds on Ruffert et al. (1996) and Bruenn (1985). We follow the procedure discussed in Neilsen et al. (2014) to compute optical depths, which is well suited for aspherical and complex geometries (such as that of an accretion disk). In the following, we briefly outline some aspects of our leakage scheme.

We specify the net neutrino heating/cooling rate per unit volume in the rest frame of the fluid, Q, and the net lepton emission/absorption rate per unit volume in the rest frame of the fluid, R (cf. Equations (1), (3), and (19)) as a local balance of absorption and emission of free-streaming neutrinos,

$$\begin{eqnarray}&&R=\displaystyle \sum _{{\nu }_{i}}{\kappa }_{{\nu }_{i}}{n}_{{\nu }_{i}}-({R}_{{\nu }_{{\rm{e}}}}^{\mathrm{eff}}-{R}_{{\bar{\nu }}_{{\rm{e}}}}^{\mathrm{eff}})\end{eqnarray} \tag{ 20 }$$

and

$$\begin{eqnarray}&&Q=\displaystyle \sum _{{\nu }_{i}}{\kappa }_{{\nu }_{i}}{n}_{{\nu }_{i}}{E}_{{\nu }_{i}}-\displaystyle \sum _{{\nu }_{i}}{Q}_{{\nu }_{i}}^{\mathrm{eff}}.\end{eqnarray} \tag{ 21 }$$

Here ${\nu }_{i}=\{{\nu }_{{\rm{e}}},{\bar{\nu }}_{{\rm{e}}},{\nu }_{x}\}$, where ${\nu }_{{\rm{e}}}$ denotes electron neutrinos, ${\nu }_{{\rm{a}}}$ denotes electron antineutrinos, and the heavy-lepton neutrinos ${\nu }_{\mu }$ and ${\nu }_{\tau }$ are collectively labeled as ${\nu }_{x}$. Furthermore, ${\kappa }_{{\nu }_{i}}$, ${n}_{{\nu }_{i}}$, and ${E}_{{\nu }_{i}}$ denote the corresponding absorption opacities, number densities, and mean energies of the free-streaming neutrinos in the rest frame of the fluid, respectively. Finally, ${R}_{{\nu }_{{\rm{e}}}}^{\mathrm{eff}}$, ${R}_{{\bar{\nu }}_{{\rm{e}}}}^{\mathrm{eff}}$, and ${Q}_{{\nu }_{i}}^{\mathrm{eff}}$ denote the corresponding effective number and energy emissivities in the rest frame of the fluid. For the present simulations, we neglect neutrino absorption, as the accretion disk simulated here remains optically thin to all neutrino species at all times (cf. Siegel & Metzger 2017). Neutrino absorption is only expected to appreciably change the outflow and disk dynamics for significantly more massive accretion disks (Fernández & Metzger 2013).

The effective emission/cooling rates ${R}_{{\nu }_{i}}^{\mathrm{eff}}$ and ${Q}_{{\nu }_{i}}^{\mathrm{eff}}$ take effects of finite optical depth into account and are computed from the intrinsic (free) emission rates ${R}_{{\nu }_{i}}$ and ${Q}_{{\nu }_{i}}$ by (cf. Equations (B22) and (B23) of Ruffert et al. 1996)

$$\begin{eqnarray}&&{R}_{{\nu }_{i}}^{\mathrm{eff}}=\displaystyle \frac{{R}_{{\nu }_{i}}}{1+\tfrac{{t}_{\mathrm{diff},{\nu }_{i}}}{{t}_{{\nu }_{i}}^{\mathrm{em},{\rm{R}}}}},{Q}_{{\nu }_{i}}^{\mathrm{eff}}=\displaystyle \frac{{Q}_{{\nu }_{i}}}{1+\tfrac{{t}_{\mathrm{diff},{\nu }_{i}}}{{t}_{{\nu }_{i}}^{\mathrm{em},{\rm{Q}}}}}.\end{eqnarray} \tag{ 22 }$$

Here

$$\begin{eqnarray}&&{t}_{\mathrm{diff},{\nu }_{i}}={D}_{\mathrm{diff}}{\kappa }_{{\nu }_{i}}^{-1}{\tau }_{{\nu }_{i}}^{2}\end{eqnarray} \tag{ 23 }$$

denote the local diffusion timescales, where ${\tau }_{{\nu }_{i}}$ are the corresponding optical depths (see below), and D_diff is a diffusion normalization factor, which we set to D_diff = 6 (O'Connor & Ott 2010). Furthermore,

$$\begin{eqnarray}&&{t}_{{\nu }_{i}}^{\mathrm{em},R}=\displaystyle \frac{{R}_{{\nu }_{i}}}{{n}_{{\nu }_{i}}},{t}_{{\nu }_{i}}^{\mathrm{em},Q}=\displaystyle \frac{{Q}_{{\nu }_{i}}}{{e}_{{\nu }_{i}}}\end{eqnarray} \tag{ 24 }$$

are the local neutrino number and energy emission timescales, where ${e}_{{\nu }_{i}}$ refers to the neutrino energy densities and

$$\begin{eqnarray}&&{R}_{{\nu }_{i}}={\delta }_{{\nu }_{i},{\nu }_{{\rm{e}}}}{R}_{{\nu }_{{\rm{e}}}}^{\beta }+{\delta }_{{\nu }_{i},{\bar{\nu }}_{{\rm{e}}}}{R}_{{\bar{\nu }}_{{\rm{e}}}}^{\beta }+{R}_{{\nu }_{i}}^{\mathrm{ee}}+{R}_{{\nu }_{i}}^{\gamma },\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{Q}_{{\nu }_{i}}={\delta }_{{\nu }_{i},{\nu }_{{\rm{e}}}}{Q}_{{\nu }_{{\rm{e}}}}^{\beta }+{\delta }_{{\nu }_{i},{\bar{\nu }}_{{\rm{e}}}}{Q}_{{\bar{\nu }}_{{\rm{e}}}}^{\beta }+{Q}_{{\nu }_{i}}^{\mathrm{ee}}+{Q}_{{\nu }_{i}}^{\gamma }\end{eqnarray} \tag{ 26 }$$

(cf. Equations (B18)–(B21) of Ruffert et al. 1996). The emission rates ${R}_{{\nu }_{i}}^{\beta }$ and ${Q}_{{\nu }_{i}}^{\beta }$, ${R}_{{\nu }_{i}}^{\mathrm{ee}}$ and ${Q}_{{\nu }_{i}}^{\mathrm{ee}}$, and ${R}_{{\nu }_{i}}^{\gamma }$ and ${Q}_{{\nu }_{i}}^{\gamma }$ are computed as in Galeazzi et al. (2013) and reflect the contributing neutrino emission mechanisms we consider. These are, respectively,

• (i)  
  charged-current β-processes,

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

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

  the strongest neutrino emission mechanism in hot and dense nuclear matter;
• (ii)  
  electron–positron pair annihilation,

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

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

  which is most relevant in nondegenerate nuclear matter at low densities and high temperatures; and
• (iii)  
  plasmon decay,

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

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

  which is efficient at intermediate densities and high temperatures.

2.2.1. Calculation of Opacities

The neutrino opacities ${\kappa }_{{\nu }_{i}}$ introduced above may be subdivided into contributions from absorption and scattering,

$$\begin{eqnarray}&&{\kappa }_{{\nu }_{i}}={\kappa }_{{\nu }_{i},\mathrm{abs}}+{\kappa }_{{\nu }_{i},\mathrm{scat}},\end{eqnarray} \tag{ 33 }$$

where

• (i)  
  ${\kappa }_{{\nu }_{i},\mathrm{abs}}$ refers to absorption of electron and anti-electron neutrinos only,

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

  $$\begin{eqnarray}&&{\bar{\nu }}_{{\rm{e}}}+p\to n+{e}^{+};\end{eqnarray} \tag{ 35 }$$

  and
• (ii)  
  ${\kappa }_{{\nu }_{i},\mathrm{scat}}$ refers to coherent scattering on heavy nuclei A and scattering on free nucleons,

  $$\begin{eqnarray}&&{\nu }_{i}+A\to {\nu }_{i}+A,\end{eqnarray} \tag{ 36 }$$

  $$\begin{eqnarray}&&{\bar{\nu }}_{i}+A\to {\bar{\nu }}_{i}+A,\end{eqnarray} \tag{ 37 }$$

  $$\begin{eqnarray}&&{\nu }_{i}+[n,p]\to {\nu }_{i}+[n,p],\end{eqnarray} \tag{ 38 }$$

  $$\begin{eqnarray}&&{\bar{\nu }}_{i}+[n,p]\to {\bar{\nu }}_{i}+[n,p].\end{eqnarray} \tag{ 39 }$$

The absorption and scattering opacities for these processes are computed as in Galeazzi et al. (2013).

2.2.2. Calculation of Optical Depths

We reduce the nonlocal computation of optical depths ${\tau }_{{\nu }_{i}}$ to an effective local problem by applying the method described in Neilsen et al. (2014), which is well suited for aspherical geometries such as an accretion disk. Global integrations are avoided by decomposing the optical depth at a given grid point into the optical depth to any neighboring point plus the already computed optical depth ${\tau }_{{\nu }_{i},\mathrm{neigh}}$ at the neighboring point, which we compute as

$$\begin{eqnarray}&&{\tau }_{{\nu }_{i},\mathrm{neigh}}+{\bar{\kappa }}_{{\nu }_{i}}{({\bar{\gamma }}_{{ab}}{{dx}}^{a}{{dx}}^{b})}^{1/2},\end{eqnarray} \tag{ 40 }$$

where ${{dx}}^{a}$ is the spatial coordinate distance vector between the two points, and ${\bar{\kappa }}_{{\nu }_{i}}$ and ${\bar{\gamma }}_{{ab}}$ denote the opacities and components of the spatial metric averaged between the two neighboring points. We define the optical depth at a given grid point as the minimum over all expressions (Equation (40)) computed for all neighboring points.

Conservative GRMHD schemes evolve the conserved variables ${\boldsymbol{q}}$ (cf. Equation (7)). This involves computing the flux terms ${{\boldsymbol{f}}}^{(i)}({\boldsymbol{p}},{\boldsymbol{q}})$ and source terms ${\boldsymbol{s}}({\boldsymbol{p}})$ for a given ${\boldsymbol{q}}$, which requires us to obtain the primitive variables ${\boldsymbol{p}}$ from the conserved ones. While the conservative variables as a function of primitive variables, ${\boldsymbol{q}}={\boldsymbol{q}}({\boldsymbol{p}})$, are given in analytic form by Equations (9)–(16), the inverse relation, ${\boldsymbol{p}}={\boldsymbol{p}}({\boldsymbol{q}})$, i.e., the recovery of primitive variables from conservative ones, is not known in closed form; this instead requires numerical inversion of the aforementioned set of nonlinear equations.

We have implemented a new framework for the recovery of primitive variables in GRHydro that provides support for any composition-dependent, finite-temperature (three-parameter) EOS, as well as a recovery scheme based on a three-dimensional Newton–Raphson solver using Equations (21), (22), and (28) in Cerdá-Durán et al. (2008). We find that this scheme has particularly fast convergence properties as compared to other schemes, typically involving a minimum of EOS calls (Siegel et al. 2018; Siegel & Mösta 2018).⁷ The latter fact is of particular importance for a three-parameter EOS, as most such EOSs are provided in the form of multidimensional tables, and table lookups can become computationally expensive. Furthermore, its ability to recover strongly magnetized regions is important for evolving low-density magnetized disk winds, as in the present simulation.

We base the microphysical description of matter at the relatively low densities and temperatures of our present simulation on the Helmholtz EOS (Timmes & Arnett 1999; Timmes & Swesty 2000), which we have newly implemented into GRHydro. Nuclear-reaction networks such as SkyNet (Lippuner & Roberts 2017), which we employ for calculating r-process abundance yields, also use the Helmholtz EOS, which is how we minimize thermodynamical inconsistencies between the simulation and subsequent postprocessing to obtain nucleosynthesis abundance yields.

The Helmholtz EOS is formulated in terms of a Helmholtz free energy, which takes into account contributions from nuclei (treated as ideal gas) with Coulomb corrections, electrons and positrons with an arbitrary degree of relativity and degeneracy, and photons in local thermodynamic equilibrium. As nuclei in the present simulation, we consider free neutrons and protons, as well as α-particles. We have modified the Helmholtz EOS to include the nuclear binding energy release from α-particle formation. We compute the abundances of nuclei at given (ρ, T, Y_e) assuming nuclear statistical equilibrium (NSE), i.e., by numerically solving the Saha equation supplemented with baryon number and charge conservation,

$$\begin{eqnarray}&&{n}_{{\rm{p}}}^{2}{n}_{{\rm{n}}}^{2}=2{n}_{\alpha }{\left(\displaystyle \frac{{m}_{{\rm{b}}}{k}_{{\rm{B}}}T}{2\pi {{\hslash }}^{2}}\right)}^{9/2}\exp (-{Q}_{\alpha }/{k}_{{\rm{B}}}T),\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&{n}_{{\rm{b}}}={n}_{{\rm{n}}}+{n}_{{\rm{p}}}+4{n}_{\alpha },\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&{n}_{{\rm{b}}}{Y}_{{\rm{e}}}={n}_{{\rm{p}}}+2{n}_{\alpha }.\end{eqnarray} \tag{ 43 }$$

Here k_B is the Boltzmann constant, ℏ is the reduced Planck constant, ${Q}_{\alpha }\simeq 28.3\,\mathrm{MeV}$ is the nuclear binding energy of an α-particle, and n_n, n_p, and n_α denote the number densities of neutrons, protons, and α-particles, respectively. We also include additional terms to the thermodynamical derivatives that arise from compositional changes with respect to (ρ, T, Y_e), i.e., from the fact that $\partial {n}_{{\rm{n}}}/\partial \rho$, $\partial {n}_{{\rm{n}}}/\partial T$, ∂n_n/∂Y_e, etc. from Equations (41)–(43) are nonzero. These additional terms can be important to the evolution code, as, e.g., the Riemann solver can depend on thermodynamic derivatives through the sound speed.

Magnetic stresses generated by turbulence mediate angular momentum transport and energy dissipation in accretion disks around compact objects. Turbulence is thought to be generated in this context by the MRI, which refers to certain exponentially growing modes that can develop in differentially rotating magnetized fluids (e.g., Velikhov 1959; Chandrasekhar 1960; Balbus & Hawley 1991, 1998; Balbus 2003; Armitage 2011). The MRI is a local instability, the growth of which is dominated by a fastest-growing MRI mode; in GRMHD, its wavelength can be estimated by (Siegel et al. 2013; Kiuchi et al. 2015b, 2017)

$$\begin{eqnarray}&&{\lambda }_{\mathrm{MRI}}\simeq \displaystyle \frac{2\pi }{{\rm{\Omega }}}\displaystyle \frac{b}{\sqrt{4\pi \rho h+{b}^{2}}},\end{eqnarray} \tag{ 45 }$$

where ${\rm{\Omega }}={u}^{\phi }/{u}^{0}$ is the angular frequency and $b\equiv \sqrt{{b}^{2}}$. The MRI is typically well resolved when λ_MRI is numerically resolved by at least 10 grid points and partially resolved with more than ∼5 grid points (e.g., Siegel et al. 2013; Kiuchi et al. 2015b).

At t = 0 ms, λ_MRI is only resolved by ∼5 grid points in the high-density region of the initial torus (cf. Figure 1, top panel). Within ≈1 ms, however, by initial relaxation and magnetic winding, the high-density part (∼10¹⁰–10¹¹ g cm⁻³) of the torus rapidly enters a regime in which λ_MRI is resolved by 10 or more grid points (cf. Figure 1, center panel). Indeed, starting at ≈1 ms, we witness the onset of magnetic field amplification in the poloidal field at the expected rate for the MRI $\propto \exp (t/{\tau }_{\mathrm{MRI}})$, where (Siegel et al. 2013)

$$\begin{eqnarray}&&{\tau }_{{\rm{M}}{\rm{R}}{\rm{I}}}\simeq \displaystyle \frac{1}{{\rm{\Omega }}}\end{eqnarray} \tag{ 46 }$$

until saturation (cf. Figure 2, top panel); the onset of the instability leads to a total amplification by roughly 1.5 orders of magnitude for the maximum poloidal magnetic field strength.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As we start with a purely poloidal magnetic field configuration, the toroidal magnetic field component first needs to be amplified by magnetic winding in order for the grid setup to resolve the MRI in the toroidal field. For the maximum toroidal magnetic field strength, this initial amplification process by magnetic winding takes a few ms (Figure 2, center panel) and slightly longer for other parts of the disk that start with smaller poloidal field strengths. Combined amplification by winding and the MRI leads to an overall increase of almost two orders of magnitude in the maximum total magnetic field strength within the first ≈5–10 ms (Figure 2, bottom panel).

By t = 20 ms, the disk has reached a quasi-stationary state, in which λ_MRI is typically resolved by 10 or more grid points (Figure 1, bottom panel). The MRI remains resolved in this way throughout the torus for the rest of the simulation, although properly resolving the MRI very close to the BH is a challenging task with current computational resources; close to the BH, we do not resolve the MRI with >10 grid points at all times and spatial points. However, we do not expect that this appreciably affects our results for the quantity and composition of the disk outflows, since these are typically generated on larger spatial scales (see, e.g., Section 4.4).

The quasi-stationary state reached at t = 20 ms and depicted in Figure 1 (bottom panel) and Figure 3 is very similar to the very early state of accretion disks obtained in recent NS–NS merger simulations. In particular, the typical magnetic field strengths of up to ∼10¹⁵ G close to the BH and disk midplane, as well as the typical magnetic-to-fluid pressure ratios of ∼10⁻³–10⁻¹ (cf. Figure 3), were also obtained by Kiuchi et al. (2015b) and Ciolfi et al. (2017). This state at t = 20 ms serves as initial data for the rest of the simulation, and all matter accreted onto the BH or ejected from the disk during the relaxation phase t < 20 ms is discarded from all further analysis.

Zoom In Zoom Out Reset image size

Strong magnetic fields ∼10¹⁵–10¹⁶ G (cf. Figure 3) can potentially modify the EOS and the neutrino emission and absorption rates (Equations (27)–(30) and (34)–(35)) through the quantization of energy levels for electrons and positrons and their motion perpendicular to the magnetic field (Lai & Qian 1998; Duan & Qian 2004, 2005). Such effects of Landau-level quantization may become relevant for densities below a critical density (Harding & Lai 2006; Haensel et al. 2007; Kiuchi et al. 2015a)

$$\begin{eqnarray}&&{\rho }_{{\rm{B}}}=2.23\times {10}^{9}{\left(\displaystyle \frac{{Y}_{e}}{0.1}\right)}^{-1}{\left(\displaystyle \frac{B}{{10}^{15}{\rm{G}}}\right)}^{3/2}\,{\rm{g}}\,{\mathrm{cm}}^{-3}\end{eqnarray} \tag{ 47 }$$

and/or below a critical temperature T_B (Harding & Lai 2006)

$$\begin{eqnarray}{T}_{{\rm{B}}}=\left\{\begin{array}{cc}\displaystyle \frac{{m}_{{\rm{e}}}{c}^{2}}{{k}_{{\rm{B}}}}\left(\sqrt{\displaystyle \frac{2B}{{B}_{Q}}+1}-1\right), & \rho \leqslant {\rho }_{{\rm{B}}}\\ \displaystyle \frac{{\hslash }{\omega }_{{\rm{c}}}}{{k}_{{\rm{B}}}}{(1+{x}_{F}^{2})}^{-1/2}, & \rho \gg {\rho }_{{\rm{B}}}\end{array}\right..\end{eqnarray} \tag{ 48 }$$

Here m_e is the electron mass, c is the speed of light, ${\omega }_{{\rm{c}}}={eB}/{m}_{{\rm{e}}}c$ is the cyclotron frequency, ${x}_{F}={\hslash }{(3{\pi }^{2}{Y}_{{\rm{e}}}\rho )}^{1/3}$ is the normalized relativistic Fermi momentum, and B_Q =4.414 × 10¹³ G is the critical QED magnetic field strength.

Figure 4 shows that, typically, ρ ≫ ρ_B and T ≳ T_B in the disk. Consequently, many Landau levels are populated, and their thermal widths are larger than the level spacing, such that the magnetic field is nonquantizing. In the polar funnel, $\rho \ll {\rho }_{{\rm{B}}}$ but still T ≳ T_B, such that, again, the magnetic field has a nonquantizing effect. Since the disk remains in this state throughout the entire simulation, we conclude that the effects of Landau-level quantization are not important for the disk evolution.

Zoom In Zoom Out Reset image size

In the neutron-rich environment of the post-merger accretion disk, one might naively expect positron captures onto neutrons, ${e}^{+}+n\to p+{\bar{\nu }}_{e}$ (Equation (28)), to be favored over electron captures (Equation (27)), such that the disk matter would protonize over viscous timescales of hundreds of ms, raising the proton/electron fraction Y_e (e.g., Metzger et al. 2009). This effect is indeed evident from Figure 5 in some portions of the disk. However, a monotonic rise of Y_e in the disk midplane raises the question of how outflows from the disk can remain sufficiently neutron-rich to synthesize heavy r-process elements, even at late times in the disk evolution. As we now describe, the reason is the existence of a self-regulation mechanism in the inner parts of the disk, which keeps a reservoir of neutron-rich material that is continuously fed into the outflows.

Zoom In Zoom Out Reset image size

Once the disk has reached a quasi-stationary state (cf. Sections 4.1 and 4.4), it regulates itself to mild electron degeneracy, which, in the presence of optically thin neutrino cooling, results in a low Y_e state (Y_e ∼ 0.1).⁸ This mechanism has been noted in the context of 1D models of neutrino-cooled accretion disks on analytical and semi-analytical grounds (Chen & Beloborodov 2007; Kawanaka & Mineshige 2007; Metzger et al. 2009), and the first evidence of self-regulation in a full 3D GRMHD simulation has been presented in Siegel & Metzger (2017). Here we elaborate on these results and discuss the mechanism in somewhat more detail; the existence of this mechanism is important for the generation of neutron-rich outflows from the disk (Section 4.4), their r-process nucleosynthesis yields (Section 5), and the resulting thermal emission (KN).

In the hot and dense accretion disk, the number densities of electrons and positrons (e^±) in thermodynamic equilibrium with the baryonic matter are given by

$$\begin{eqnarray}&&{n}_{\pm }=\displaystyle \frac{{({m}_{{\rm{e}}}c)}^{3}}{{\pi }^{2}{{\hslash }}^{3}}{\int }_{1}^{\infty }{f}_{\pm }(E,T,\mu )E\sqrt{{E}^{2}-1}\,{dE},\end{eqnarray} \tag{ 49 }$$

where E is the relativistic particle energy in units of ${m}_{{\rm{e}}}{c}^{2}$. Here f_± is the Fermi–Dirac function,

$$\begin{eqnarray}&&{f}_{\pm }(E,T,\mu )=\displaystyle \frac{1}{\exp [(E\pm \mu )/{\rm{\Theta }}]+1},\end{eqnarray} \tag{ 50 }$$

where ${\rm{\Theta }}={k}_{{\rm{B}}}T/{m}_{{\rm{e}}}{c}^{2}$ and $\mu \equiv {\mu }_{-}=-{\mu }_{+}$ is the electron chemical potential in units of ${m}_{{\rm{e}}}{c}^{2}$. Charge neutrality requires that

$$\begin{eqnarray}&&{n}_{-}-{n}_{+}={Y}_{{\rm{e}}}{n}_{{\rm{b}}},\end{eqnarray} \tag{ 51 }$$

which, together with Equation (49), determines μ and n_± at a given thermodynamic state $(\rho ,T,{Y}_{{\rm{e}}})$. For degenerate relativistic matter ($\mu /{\rm{\Theta }}\gg 1$), using the Sommerfeld expansion of Equation (49) in terms of $\mu /{\rm{\Theta }}$, one can show that the temperature dependence of μ is approximately given by (see Appendix)

$$\begin{eqnarray}&&\sqrt{{\mu }^{2}-1}=\sqrt{{E}_{{\rm{F}}}^{2}-1}\left(1-\displaystyle \frac{{\pi }^{2}}{6}\displaystyle \frac{{{\rm{\Theta }}}^{2}}{{E}_{{\rm{F}}}^{2}-1}\right),\end{eqnarray} \tag{ 52 }$$

where ${E}_{{\rm{F}}}\equiv \mu (T=0)$ is the Fermi energy. Furthermore, for degenerate matter, free e^± pairs can only be obtained from around the Fermi edge E ≃ μ with width ${\rm{\Delta }}E\simeq 4{\rm{\Theta }}$, which is very narrow (${\rm{\Delta }}E/E\simeq 4{\rm{\Theta }}/\mu \ll 1$); from Equation (49), one finds that for

$$\begin{eqnarray}&&(\mu /{\rm{\Theta }}\gg 1,E\simeq \mu )\displaystyle \frac{{n}_{+}}{{n}_{-}}\propto \exp (-2\mu /{\rm{\Theta }}),\end{eqnarray} \tag{ 53 }$$

i.e., e^± creation is heavily suppressed. Higher electron degeneracy η ≡ μ/Θ results in less electrons and positrons (cf. Equations (49) and (53)). This decreases the neutrino emission via charged-current interactions and pair annihilation (cf. Equations (27)–(30)); i.e., it results in a lower cooling rate and higher temperatures. Higher temperatures, in turn, decrease μ (cf. Equation (52)) and thus increase the degeneracy, i.e., η. Because of this negative feedback loop, whenever the disk enters the (strongly) degenerate regime, it will tend to self-regulate its degeneracy and maintain a state of mild electron degeneracy η ∼  1. Indeed, as shown by Figure 5, soon after reaching the quasi-stationary state, the disk has regulated itself to mild degeneracy η ∼ 1 in the inner parts of the disk in which neutrino cooling is energetically important (r ≲ 60 km or r ≲ 14 gravitational radii) and qualitatively remains in this state throughout the remainder of the simulation.

In the hot and dense matter of the inner parts of the disk, electron and positron capture (cf. Equations (27) and (28)) are the dominant cooling reactions. The equilibrium Y_e that results from conditions of mild degeneracy in this neutrino-transparent matter is then determined by equal rates of e^± capture,

$$\begin{eqnarray}&&{\dot{n}}_{{e}^{-}p}={\dot{n}}_{{e}^{+}n};\end{eqnarray} \tag{ 54 }$$

i.e., Equations (49), (51), and (54) determine Y_e for a given ρ and T. For mild degeneracy η ≳ 1, one can show that from Equation (54), the equilibrium Y_e is approximately given by (Beloborodov 2003)

$$\begin{eqnarray}&&{Y}_{{\rm{e}}}=0.5+\displaystyle \frac{7{\pi }^{4}}{1350\zeta (5)}\left(\displaystyle \frac{Q}{2{\rm{\Theta }}}-\eta \right)\end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray}&&=0.5+0.487\left(\displaystyle \frac{1.2655}{{\rm{\Theta }}}-\eta \right),\end{eqnarray} \tag{ 56 }$$

where ζ is the Riemann ζ-function and $Q=({m}_{{\rm{n}}}-{m}_{{\rm{p}}})/{m}_{{\rm{e}}}\,=2.531$ is the neutron–proton mass difference in units of the electron mass. A very mild electron degeneracy η ⪆ 1 in hot matter Θ ≈ 1 is therefore sufficient to generate conditions of neutron richness Y_e < 0.5. For the hot Θ ≳ 1 and mildly degenerate conditions η ≳ 1 of the inner parts of the disk, the resulting neutron richness adjusts to an equilibrium value of typically Y_e ∼ 0.1 or lower (see Figure 5).

The presence of this self-regulation mechanism to mild electron degeneracy, which implies a low Y_e ∼ 0.1, is important to allow for the generation of neutron-rich outflows that can undergo r-process nucleosynthesis (Sections 4.4 and 5). It forces the disk to keep a reservoir of neutron-rich material despite the ongoing protonization process in the rest of the disk—neutron-rich material that is continuously fed into the outflows to keep the overall mean electron fraction ${\bar{Y}}_{{\rm{e}}}$ of the outflow rather low over the lifetime of the disk (${\bar{Y}}_{{\rm{e}}}\sim 0.2$; see Table II of Siegel & Metzger 2017 and Section 5.2). This results in the possibility of generating a robust second-to-third-peak r-process (cf. Section 5) and thus the production of a significant amount of lanthanide material in the outflow. Due to its high opacity, this material can then produce a red KN, as observed in the recent GW170817 event.

Magnetic stresses generated by MHD turbulence via the MRI mediate angular momentum transport and thus energy dissipation in the disk. Turbulence also dissipates magnetic energy, which, however, is regenerated through a dynamo (e.g., Parker 1955; Brandenburg et al. 1995). The balance of the two processes results in a saturated steady turbulent, quasi-equilibrium state, which is characterized by a roughly constant ratio of magnetic to internal energy in the disk.

Figure 6 shows the temporal evolution of the density-averaged ratio of electromagnetic to internal energy $\langle {e}_{\mathrm{EM}}/{e}_{\mathrm{int}}{\rangle }_{\hat{D}}$ and the magnetic-to-fluid pressure ratio $\langle {p}_{{\rm{B}}}/{p}_{{\rm{f}}}{\rangle }_{\hat{D}}$, which are indeed indicative of a disk in a steady turbulent state. We define the rest-mass density average of a quantity χ by

$$\begin{eqnarray}&&\langle \chi {\rangle }_{\hat{D}}\equiv \displaystyle \frac{\int \chi \hat{D}{d}^{3}x}{\int \hat{D}{d}^{3}x},\end{eqnarray} \tag{ 57 }$$

where $\hat{D}=\sqrt{\gamma }\rho W$ is the conserved rest-mass density (cf. Equations (7)–(9)).⁹ Following Duez et al. (2006), we define the total internal energy

$$\begin{eqnarray}&&{E}_{\mathrm{int}}\equiv \int \epsilon \rho W\sqrt{\gamma }{d}^{3}x\end{eqnarray} \tag{ 58 }$$

and the total electromagnetic energy

$$\begin{eqnarray}&&{E}_{\mathrm{EM}}\equiv \int {n}_{\mu }{n}_{\nu }{T}_{\mathrm{EM}}^{\mu \nu }\sqrt{\gamma }{d}^{3}x,\end{eqnarray} \tag{ 59 }$$

where ${T}_{\mathrm{EM}}^{\mu \nu }$ is the electromagnetic part of the energy–momentum tensor. We thus define the local ratio of electromagnetic to internal energy by

$$\begin{eqnarray}&&\displaystyle \frac{{e}_{\mathrm{EM}}}{{e}_{\mathrm{int}}}\equiv \displaystyle \frac{{n}_{\mu }{n}_{\nu }{T}_{\mathrm{EM}}^{\mu \nu }}{\epsilon \rho W}.\end{eqnarray} \tag{ 60 }$$

Figure 6 shows that for t > 20 ms, this ratio remains roughly constant in a time-averaged sense and thus indicates that a steady turbulent state of the disk is indeed achieved and maintained. Furthermore, Figure 6 shows that

$$\begin{eqnarray}&&{\left\langle \displaystyle \frac{{p}_{{\rm{B}}}}{{p}_{{\rm{f}}}}\right\rangle }_{\hat{D}}\simeq 0.1,\end{eqnarray} \tag{ 61 }$$

which is also characteristic of such a steady turbulent state (e.g., Jiang et al. 2014b; Sa̧dowski et al. 2015). This ratio in the nonlinear saturated state is much larger than the initial value of p_B/p_f < 5 × 10⁻³ (cf. Section 3 and Table 1).

Zoom In Zoom Out Reset image size

The 3D nature of our disk simulation is crucial for generating a steady turbulent state. Due to the antidynamo theorem (Cowling 1933), magnetic fields cannot be regenerated by dynamo action in axisymmetry, and a steady turbulent state cannot thus be maintained.

Direct evidence for dynamo action in our disk simulation is depicted in the top panel of Figure 7, which shows a spacetime diagram of the radially averaged y-component of the magnetic field in the x–z plane. This "butterfly" diagram clearly indicates the presence of magnetic cycles with a period of roughly ∼20 ms throughout the entire simulation time domain. In the disk midplane, magnetic fields of temporally alternating polarity are generated by MHD turbulence. These fields slowly migrate off the midplane by magnetic pressure gradients and buoyancy, where they are gradually dissipated into heat. This migration and dissipation of magnetic energy contributes to establishing a "hot" corona above and below the midplane, as indicated by the middle panel of Figure 7. This spacetime diagram of the specific entropy shows strongly increasing specific entropies off the midplane where magnetic field strengths decrease. We note that the temperature, however, decreases as a function of height off the midplane. Therefore, the production of high-energy nonthermal neutrinos in the corona by upscattering of thermal neutrinos emitted from the midplane (cf. bottom panel of Figure 7) is not expected.¹⁰

Zoom In Zoom Out Reset image size

In the hot corona, powerful outflows are generated. In these regions of lower density, viscous heating from MHD turbulence and dissipation of magnetic energy exceeds cooling by neutrino emission, which is strongest in the disk midplane (cf. Figure 7, bottom panel). This heating-cooling imbalance results in launching neutron-rich winds from the disk. Above and below the midplane, the neutrino emissivities decrease as functions of "height" $| z|$, and the weak interactions (and thus Y_e) essentially "freeze out"; however, further mixing in the (initially turbulent) outflows can still change Y_e.

The outflows are tracked by 10⁴ passive tracer particles that are advected with the plasma. These tracer particles are of equal mass, placed within the initial torus at t = 0 ms with a probability proportional to the conserved rest-mass density $\hat{D}=\sqrt{\gamma }\rho W$. We distinguish between total outflow, defined as the entity of all tracer particles that have reached a radial coordinate distance of 10³ km from the center of the BH by the end of the simulation, and unbound outflow, or ejecta, defined as the entity of tracer particles that are additionally unbound according to the Bernoulli criterion $-{{hu}}_{0}\gt 1$ (nonvanishing escape velocity at infinity).

Outflows are generated over a wide range of radii. This is illustrated by the top panel of Figure 8, which shows mass histograms of the outflow tracer particles in terms of their cylindrical coordinate radii $\varpi =\sqrt{{x}^{2}+{y}^{2}}$ at the time of ejection from the disk, ${\varpi }_{\mathrm{ej}}\equiv \varpi (t={t}_{\mathrm{ej}})$. We define the time of ejection from the disk or corona t = t_ej as the time after which the radial coordinate position of a tracer particle $r=\sqrt{{x}^{2}+{y}^{2}+{z}^{2}}$ only increases monotonically with time. The total outflow shows a broad distribution with significant mass being ejected between ${\varpi }_{\mathrm{ej}}\approx 20\,\mathrm{km}$ and ϖ_ej > 600 km from the BH. However, we find that mass ejection is most efficient in a narrower range of ejection radii, as indicated by the histogram of unbound matter, the latter being ejected essentially in the range ${\varpi }_{\mathrm{ej}}\approx 100\mbox{--}400\,\mathrm{km}$ from the BH.

Zoom In Zoom Out Reset image size

Matter is typically unbound by recombination into α-particles. The imbalance of heating and cooling in the hot corona, as mentioned above, lifts material in the BH potential but typically only leads to marginally bound or unbound outflows. Subsequent nuclear binding energy release from recombination of free nucleons into α-particles rapidly generates specific enthalpy as matter approaches the recombination temperature and immediately "unbinds" the material; this is shown in Figure 9 for a few representative tracer particles. A spike in the specific enthalpy h is created by internal energy that becomes available during the recombination process (7 MeV per baryon per α-particle produced) plus the resulting pressure increase in a low-density environment. For a stationary relativistic fluid flow (isentropic, constant specific angular momentum), hu₀ is constant along a fluid world line (Equation (44)). As the material moves away from the disk, the outflows cool ($h\to 1$) and specific enthalpy is converted into kinetic energy keeping hu₀ constant, which sets the asymptotic escape velocity.

Zoom In Zoom Out Reset image size

The bottom panel of Figure 8 shows the distribution of kinetic energy of the unbound and total outflows in terms of their outflow velocities. We characterize the outflow by two velocities: v_{1000 km}, the velocity at a coordinate distance r = 10³ km from the BH, and ${v}_{\infty }$, the corresponding asymptotic escape velocity when the conversion of internal energy to kinetic energy has been completed. Here ${v}_{\infty }$ is computed from the corresponding asymptotic Lorentz factor ${W}_{\infty }\equiv -{{hu}}_{0}$, where hu₀ is evaluated either when the tracer particle leaves the computational domain or at the final time of the simulation if it stays inside the computational domain for the entire simulation time. Unbound and total outflows have similar velocity distributions in the range v_{1000 km} ≈ (0.03–0.15)c. The kinetic energy–weighted mean outflow velocities ${\bar{v}}_{1000\mathrm{km}}\,\equiv \sqrt{2{E}_{\mathrm{kin},\mathrm{tot}}/{M}_{\mathrm{ej}}}$ are $0.063c$ and 0.058 for unbound and total outflow, respectively. Here E_kin,tot denotes the total kinetic energy in the outflow type, and M_ej is the total mass of the outflow type. The asymptotic kinetic energy distribution of the unbound outflow, however, shows ${v}_{\infty }\approx (0.04\mbox{--}0.25)c$, with a higher kinetic energy–weighted mean of ${\bar{v}}_{\infty }=0.094c\approx 0.1c$.

Though not included in our simulations, the outflows will receive additional nuclear heating from the r-process on larger radial scales of ≈2–3 MeV per nucleon (Metzger et al. 2010a), which will boost its speed by an additional ≈10%–20%. We note that ${\bar{v}}_{\infty }$ of the unbound outflow corresponds to the kinetic energy–averaged value ${v}_{\mathrm{KN}}\approx 0.1c$, similar to that required to explain the red KN component observed in the recent GW170817 event (e.g., Chornock et al. 2017; Villar et al. 2017).

The total unbound mass from the disk at the end of the simulation amounts to ≈20% of its initial value. However, the true total ejecta mass, including late times after the simulation has terminated, is likely to be roughly twice as great, as estimated in greater detail in the following subsection. Additional properties of the outflow are summarized in Siegel & Metzger (2017).

The global disk structure as characterized by the radial profile of the vertical density scale height is shown in Figure 10. We define the scale height according to

$$\begin{eqnarray}&&{z}_{H}(\varpi )\equiv \langle | z| {\rangle }_{\hat{D},\mathrm{cyl}},\end{eqnarray} \tag{ 62 }$$

where

$$\begin{eqnarray} &&\langle \chi { \rangle }_{\hat{D},{\rm{c}}{\rm{y}}{\rm{l}}}\equiv \displaystyle \frac{\int {\int }_{0}^{2\pi }\chi \hat{D}\varpi d\phi dz}{\int {\int }_{0}^{2\pi }\hat{D}\varpi d\phi dz}\end{eqnarray} \tag{ 63 }$$

is the rest-mass density average of a quantity χ over azimuthal angle ϕ and height z as a function of the cylindrical coordinate radius ϖ.

Zoom In Zoom Out Reset image size

At large radii, ϖ ≳ 250 km, the disk remains geometrically thick at all times, with a density scale height of ${z}_{H}/\varpi \gtrsim 0.4\mbox{--}1$. This is because neutrino cooling is always inefficient in these low-density regions, as illustrated by the radial profile of the density-averaged electron neutrino emission rate $\langle {Q}_{{\nu }_{{\rm{e}}}}^{\mathrm{eff}}{\rangle }_{\hat{D},\mathrm{cyl}}$ in Figure 11. At late times, t > 200 ms, the density scale height z_H/ϖ exceeds unity in the radial region ϖ ≈ 100–300 km, which is due to the outflows being efficiently generated at these radii (see Section 4.4, Figure 8). The thickening of the disk as the accretion drops and the concomitant generation of outflows was predicted by 1D (height-integrated) models (Metzger et al. 2008a, 2009).

Zoom In Zoom Out Reset image size

The disk becomes thinner at smaller radii, starting at the characteristic radius ϖ_α, where α-particles dissociate into free nucleons. The α-dissociation consumes 7 MeV per nucleon, which acts to cool the accretion flow and results in a geometrically thinner disk. This radius is initially at ϖ_α ≈170 km and decreases to ϖ_α ≈ 100 km by the end of the simulation, as indicated by the radial profile of the density-averaged α-particle mass fraction $\langle {X}_{\alpha }{\rangle }_{\hat{D},\mathrm{cyl}}$ (cf. Figure 10 and the bottom panel of Figure 11).

At yet smaller radii, the accretion flow becomes geometrically even thinner as the result of neutrino cooling, with the density scale height z_H/ϖ ∼ 0.1 close to the BH, $\varpi \lessapprox 70\,\mathrm{km}$ (cf. Figure 10). This efficient neutrino cooling begins interior to the so-called "ignition" radius ϖ_ign < ϖ_α, which is defined as the location where the neutrino-cooling timescale becomes less than the local accretion timescale (Chen & Beloborodov 2007). This radius typically coincides with the location at which the energies of electrons and positrons become comparable to the neutron–proton mass difference $({m}_{{\rm{n}}}-{m}_{{\rm{p}}}){c}^{2}$, triggering the onset of the efficient Urca cooling reactions (Equations (27) and (28); see Figure 11, top panel). The same weak interactions typically result in further reduction in the electron fraction Y_e due to the increased degeneracy of the matter, as discussed in the previous subsection (cf. Figure 11, middle panel).

By the end of the simulation, the BH has accreted ≈60% of the initial torus mass. The BH accretion rate as measured by the mass flux through spherical coordinate detector surfaces is shown in Figure 12 (top panel). It decreases from $\sim 1\,{M}_{\odot }\,{{\rm{s}}}^{-1}$ at early times to $\sim {10}^{-4}\,{M}_{\odot }\,{{\rm{s}}}^{-1}$ by the end of the simulation. This leads to an essentially converged total accreted mass onto the BH of ≈1.20 × 10⁻² M_⊙ or ≈0.59 M_t,in. Here M_t,in = 2.02 × 10⁻² M_⊙ is the initial disk mass at t = 20 ms, excluding all matter that is accreted onto the BH or ejected from the disk during the initial relaxation phase (cf. Section 4.1). As the accretion rate continues to decrease as the disk viscously spreads outward (see below), the total accreted disk mass is unlikely to increase by a significant amount during the subsequent evolution.

Zoom In Zoom Out Reset image size

The MHD turbulence mediates angular momentum transport in the disk, which leads to accretion onto the BH but also to viscous radial spreading of the disk. Evidence for the latter effect is reported in the bottom panel of Figure 12, which shows that the density-averaged cylindrical radius $\langle \varpi {\rangle }_{\hat{D}}$ of matter in the simulation domain is monotonically growing after the initial relaxation phase. The same result is obtained when the disk corona and winds are explicitly excluded from the integration, i.e., by only integrating up to the local density scale height z_H of the disk (Equation (62)). However, equatorial winds are not straightforward to distinguish from the disk itself and thus remain in the analysis either way.

About ≈40% of the initial disk mass is unbound in outflows, which undergo r-process nucleosynthesis (Section 5). By the end of the simulation, roughly ≈20% of the initial disk mass has already been ejected from the disk; i.e., it has reached >1000 km and is unbound (cf. Section 4.4 and Table II of Siegel & Metzger 2017). However, the disk is still producing steady winds by the end of the simulation, which means the total unbound mass is likely to become significantly higher. Even as the disk dilutes with time and neutrino cooling becomes less important, viscous heating will still continue to drive winds. Furthermore, as the disk viscously spreads, additional material is lifted out of the BH potential, aided by nuclear binding energy release from the formation of α-particles and heavier nuclei as the material cools. With the total accreted mass having already converged, it is thus reasonable to assume that the remaining disk mass by the end of the simulation will eventually be evaporated, leading to an estimated total ejected mass of ≲0.4 M_t,in.

The inner parts of the disk are sufficiently hot and dense that neutrino emission becomes energetically important (cf. Figure 11 and Section 4.5). In this section, we discuss the characteristics of the neutrino radiation from the disk, which will serve as input to our r-process nucleosynthesis calculations presented in the next section.

We define the total neutrino luminosity for each neutrino species ${\nu }_{i}\in \{{\nu }_{{\rm{e}}},{\bar{\nu }}_{{\rm{e}}},{\nu }_{x}\}$ according to (cf. Equations (19) and (21))

$$\begin{eqnarray}&&{L}_{{\nu }_{i}}=\int \alpha {WQ}_{{\nu }_{i}}^{{\rm{e}}{\rm{f}}{\rm{f}}}\alpha \sqrt{\gamma }{d}^{3}x,\end{eqnarray} \tag{ 64 }$$

where an additional factor α is included to correct for the gravitational redshift due to the BH potential. This definition takes into account the effects of finite optical depth; i.e., it is based on the effective energy emission rates, but it neglects reabsorption of emitted neutrinos by matter.

Neutrino emission is purely thermal, characterized by the local emission temperature T (the temperature of matter). We assign mean neutrino emission temperatures for the different neutrino species to the disk, defined as the neutrino energy emission rate averaged quantities

$$\begin{eqnarray}&&{\bar{T}}_{{\nu }_{i}}\equiv \langle T{\rangle }_{{Q}_{{\nu }_{i}}}.\end{eqnarray} \tag{ 65 }$$

Here we have defined the neutrino emission rate average of a quantity χ by

$$\begin{eqnarray}&&\langle \chi {\rangle }_{{Q}_{{\nu }_{i}}}\equiv \displaystyle \frac{\int \chi {Q}_{{\nu }_{i}}^{\mathrm{eff}}W\alpha \sqrt{\gamma }{d}^{3}x}{\int {Q}_{{\nu }_{i}}^{\mathrm{eff}}W\alpha \sqrt{\gamma }{d}^{3}x}.\end{eqnarray} \tag{ 66 }$$

Note that ${Q}_{{\nu }_{i}}^{\mathrm{eff}}W\alpha \sqrt{\gamma }$ corresponds to the energy emitted per unit time and coordinate volume through neutrinos of species ${\nu }_{i}$ as seen by the Eulerian observer (cf. Equations (19) and (21)). For further reference, we also define a corresponding spherical blackbody emission radius,

$$\begin{eqnarray}&&{r}_{{\nu }_{i}}={\left(\displaystyle \frac{{L}_{{\nu }_{i}}}{4\pi \tfrac{7}{16}\sigma {\bar{T}}_{{\nu }_{i}}^{4}}\right)}^{\tfrac{1}{2}},\end{eqnarray} \tag{ 67 }$$

where σ is the Stefan–Boltzmann constant and the actual characteristic neutrino emission radius

$$\begin{eqnarray}&&{R}_{\mathrm{em},{\nu }_{i}}\equiv \langle \varpi {\rangle }_{{Q}_{{\nu }_{i}}}.\end{eqnarray} \tag{ 68 }$$

Figure 13 shows the total neutrino luminosities, average neutrino emission temperatures, and blackbody as well as characteristic emission radii as extracted from our simulation data. We extrapolate these quantities beyond the end of the simulation at t = 381 ms by power laws fitted to the late-time simulation data.

Zoom In Zoom Out Reset image size

The neutrino luminosities are initially high, with ${L}_{\nu }\,\sim {10}^{52}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ for electron and anti-electron neutrinos and at least an order of magnitude lower for the heavier neutrino species, but they quickly fade over timescales of hundreds of ms. We note that these initial neutrino luminosities are very similar to the values found in the early post-merger accretion systems of recent hydrodynamic NS–NS and BH–NS merger simulations (e.g., Foucart et al. 2017; Radice et al. 2016; Sekiguchi et al. 2016). The total energy radiated in neutrinos by the disk in terms of the various neutrino species is given by ${E}_{{\nu }_{{\rm{e}}}},{E}_{{\bar{\nu }}_{{\rm{e}}}},{E}_{{\nu }_{x}}=(4.2,6.1,0.083)\times {10}^{50}\,\mathrm{erg}$. Despite the fact that the neutrino luminosities fade rapidly compared to the evolution timescale of the disk, irradiation by neutrinos during the early phase of the evolution can still have an appreciable effect on the composition of the disk outflows and thus on r-process nucleosynthesis. We discuss this effect in the following section.

Previous r-process nucleosynthesis analyses of disk outflows from 2D Newtonian α-disk simulations have noted an overproduction of A = 132 nuclei with respect to the second r-process peak (A = 128–130) when compared to observed solar system abundances (Wu et al. 2016). This was ascribed to late-time, low-temperature convection in the disk outflow, i.e., fluid elements, whose ejection time t_ej (cf. Section 4.4) from the disk is much greater than ${t}_{5\mathrm{GK}}$. We define t_5GK as the last time the temperature of a fluid element (tracer particle) decreased below 5 GK, which is the characteristic temperature for NSE to break down and the r-process to set in.

Although our 3D GRMHD setup is expected to show less large-scale, low-temperature convection than 2D viscous hydrodynamics (because of the inverse turbulent cascade in 2D), we still find an overproduction at A = 132, which is evident from Figure 15.

In contrast to Wu et al. (2016), we find that this anomaly in our 3D GRMHD setup is not predominantly due to tracers that undergo late-time low-temperature convection, i.e., for which ${t}_{\mathrm{ej}}\gg {t}_{5\mathrm{GK}}$. This is shown in Figure 14, which reports t_ej versus t_5GK for all unbound tracer particles. The dominant contributors to this anomaly all follow the main correlation between t_ej and t_5GK, and tracers with ${t}_{\mathrm{ej}}\gg {t}_{5\mathrm{GK}}$ are not among them. The origin of this anomaly remains inconclusive at this point. It may point to a nuclear origin, at least for our present calculations with SkyNet, which requires further investigation concerning the nuclear physics input.

Zoom In Zoom Out Reset image size

In order to explore the effects of neutrino absorption on r-process nucleosynthesis in the ejecta material, we "light-bulb" irradiate the ejecta by neutrinos from the disk in a postprocessing step, employing two different assumptions to bracket the uncertainties in the neutrino emission geometry.

Spherical blackbody. In a first approach, following Roberts et al. (2017), we assume that neutrinos are emitted with luminosity ${L}_{{\nu }_{i}}$ and temperature ${\bar{T}}_{{\nu }_{i}}$ from a single spherical surface centered on the BH of radius ${r}_{{\nu }_{i}}$ (cf. Equations (64), (65), and (67)) and that they follow a Fermi–Dirac distribution in energy space,

$$\begin{eqnarray}&&{f}_{\mathrm{FD}}(E,{\bar{T}}_{{\nu }_{i}})=\displaystyle \frac{1}{\exp (E/{k}_{{\rm{B}}}{\bar{T}}_{{\nu }_{i}})+1},\end{eqnarray} \tag{ 69 }$$

where E denotes the neutrino energy. The radii of the neutrinospheres ${r}_{{\nu }_{i}}$ are typically on the order of tens of km and are roughly comparable to or smaller than the actual radii ${R}_{\mathrm{em},{\nu }_{i}}$ of the peak neutrino emission within the disk (see Figure 13, bottom panel). The neutrino distribution function in energy space as a function of coordinate radius r for species ${\nu }_{i}$ is then given by

$$\begin{eqnarray}&&{f}_{{\nu }_{i}}(E,r;{\bar{T}}_{{\nu }_{i}},{L}_{{\nu }_{i}})=\displaystyle \frac{1}{2}\left(1-\sqrt{1-\displaystyle \frac{{r}_{{\nu }_{i}}^{2}}{{r}^{2}}}\right){f}_{\mathrm{FD}}(E,{\bar{T}}_{{\nu }_{i}}).\end{eqnarray} \tag{ 70 }$$

Ringlike blackbody. In a second approach, following the neutrino emission geometry of Fernández & Metzger (2013), we assume that neutrinos are emitted with luminosity ${L}_{{\nu }_{i}}$ and temperature ${\bar{T}}_{{\nu }_{i}}$ from a ring of radius ${R}_{\mathrm{em},{\nu }_{i}}$ in the equatorial plane around the BH (cf. Equations (64), (65), and (68)). This geometry more closely resembles neutrino emission from the disk, as most of the emission is confined to regions close to the midplane (cf. Figure 7, bottom panel) and the effective emission rates ${Q}_{{\nu }_{i}}^{\mathrm{eff}}$ are indeed sharply peaked around some characteristic emission radius $r\simeq {R}_{\mathrm{em},{\nu }_{i}}$ (cf. Figure 11, top panel). In analogy to Equation (70), the neutrino distribution function in this case is given by

$$\begin{eqnarray}&&{f}_{{\nu }_{i}}(E,r,\theta ;{\bar{T}}_{{\nu }_{i}},{L}_{{\nu }_{i}},{R}_{\mathrm{em},{\nu }_{i}})=\displaystyle \frac{1}{2}{N}_{{\nu }_{i}}{{ \mathcal I }}_{{\nu }_{i}}{f}_{\mathrm{FD}}(E,{\bar{T}}_{{\nu }_{i}}),\end{eqnarray} \tag{ 71 }$$

where

$$\begin{eqnarray}&&{N}_{{\nu }_{i}}=\displaystyle \frac{{L}_{{\nu }_{i}}}{4\pi {R}_{\mathrm{em},{\nu }_{i}}^{2}\tfrac{7}{16}\sigma {\bar{T}}_{{\nu }_{i}}^{4}}\end{eqnarray} \tag{ 72 }$$

and

$$\begin{eqnarray}&&{{ \mathcal I }}_{{\nu }_{i}}=\displaystyle \frac{1}{2\pi }{\left(\displaystyle \frac{{R}_{\mathrm{em},{\nu }_{i}}}{r}\right)}^{2}{\int }_{0}^{2\pi }\displaystyle \frac{d{\phi }_{{\rm{R}}}}{2D(r,\theta ,{R}_{\mathrm{em},{\nu }_{i}},{\phi }_{{\rm{R}}})/{r}^{2}}.\end{eqnarray} \tag{ 73 }$$

Here r and θ denote the radial coordinate and polar angle, respectively, and ${\phi }_{{\rm{R}}}$ denotes the azimuthal angle that parameterizes the neutrino emission ring. Furthermore,

$$\begin{eqnarray}&&D=r{\left[1+{\left(\displaystyle \frac{{R}_{\mathrm{em},{\nu }_{i}}}{r}\right)}^{2}-2\displaystyle \frac{{R}_{\mathrm{em},{\nu }_{i}}}{r}\sin \theta \cos {\phi }_{{\rm{R}}}\right]}^{1/2}\end{eqnarray} \tag{ 74 }$$

is the distance between a spatial point (r, θ) and the neutrino emission ring at position ϕ_R (cf. Figure B2 of Fernández & Metzger 2013).

Figure 15 reports detailed abundance yields, including neutrino absorption, computed with the two methods outlined above, in comparison to previous results obtained by neglecting neutrino absorption (Siegel & Metzger 2017). It is reassuring that these results do not depend on the method by which neutrino absorption is included; both approaches lead to essentially the same abundance yields. This is not surprising, given that the source of neutrino radiation with a diameter of essentially 60–80 km is sufficiently compact compared to the spatial size of the entire disk and outflows (cf. Section 4.5).

Zoom In Zoom Out Reset image size

With neutrino absorption included, the production of the entire range of r-process nuclei from the first to the third peak of the r-process can be explained. Including neutrino absorption dramatically improves the agreement between the abundance yields of the lighter nuclei from the first to the second r-process peak (A ∼ 80–120) compared to the observed solar system abundances. This is due to neutrinos irradiating part of the outflow and the outer parts of the disk, thereby raising Y_e in part of the outflow (see Figure 16), which enhances the production of lighter r-process nuclei. However, a strong second-to-third-peak r-process is still maintained. The fact that the outflow well reaches the production of third-peak elements at the required level to explain solar abundances, even in the presence of strong neutrino irradiation, is at least in part due to the self-regulation mechanism discussed in Section 4.3, which continuously releases very neutron-rich material into the outflow. The excellent agreement with observed abundances is also reflected in the bottom panel of Figure 15, which compares the abundance yields from our simulation including neutrino absorption with observed abundances in metal-poor stars in the halo of the Milky Way.

Zoom In Zoom Out Reset image size