The r-process nucleosynthesis in the outflows from short GRB accretion disks

TL;DR

This study uses relativistic magneto-hydrodynamic simulations to analyze the formation, composition, and nucleosynthesis of outflows from accretion disks around black holes in short gamma-ray bursts, explaining kilonova observations.

Contribution

It provides new insights into the properties of outflows from post-merger accretion disks, including their composition, velocity, and contribution to kilonova lightcurves, using realistic nuclear equations of state.

Findings

Fast wind outflows with velocities 0.11-0.23c are launched.

Outflows have electron fractions Y_e of 0.1-0.4, explaining kilonova components.

Mass loss from the disk ranges from 2% to 17% of initial disk mass.

Abstract

Short gamma-ray bursts require a rotating black hole, surrounded by a magnetized relativistic accretion disk, such as the one formed by coalescing binary neutron stars or neutron star - black hole systems. The accretion onto a Kerr black hole is the mechanism of launching a baryon-free relativistic jet. An additional uncollimated outflow, consisting of sub-relativistic neutron-rich material which becomes unbound by thermal, magnetic and viscous forces, is responsible for blue and red kilonova. We explore the formation, composition and geometry of the secondary outflow by means of simulating accretion disks with relativistic magneto-hydrodynamics and employing realistic nuclear equation of state. We calculate the nucleosynthetic r-process yields by sampling the outflow with a dense set of tracer particles. Nuclear heating from the residual r-process radioactivities in the freshly synthesized nuclei is expected to power a red kilonova, contributing independently from the dynamical ejecta component, launched at the time of merger, and neutron-poor broad polar outflow, launched from the surface of the hypermassive neutron star by neutrino wind. Our simulations show that both magnetisation of the disk and high black hole spin are able to launch fast wind outflows ($v/c\sim 0.11-0.23$) with a broad range of electron fraction $Y_{\rm e}\sim 0.1-0.4$, and help explain the multiple components observed in the kilonova lightcurves. The total mass loss from the post-merger disk via unbound outflows is between 2\% and 17\% of the initial disk mass.

Tables

Table 1. Table 1: Summary of the models

Model	Torus mass
$M_{\odot }$ units
	BH spin
	$r_{\mathrm{in}}$
	$r_{\max }$
	${\beta }_0$	${\dot{M}}_{out}$
	$\Delta M_{out}$	$M_{\odot }$ s-1 units
	$M_{\odot }$ units
HS-Therm	0.1031	0.0848	0.9	3.8	9.75	100	$7.64\cdot {10}^{-2}$	$4.42\cdot {10}^{-3}$
HS-Magn	0.1031	0.0741	0.9	3.8	9.75	10	$8.26\cdot {10}^{-2}$	$1.72\cdot {10}^{-2}$
LS-Therm	0.1104	0.0895	0.6	4.5	9.1	100	$8.76\cdot {10}^{-2}$	$2.29\cdot {10}^{-3}$
LS-Magn	0.1104	0.0682	0.6	4.5	9.1	10	$1.20\cdot {10}^{-1}$	$1.81\cdot {10}^{-2}$

Table 2. Table 2: Thermodynamic properties of the outflows

Model	Unbound outflow mass [$M_{\odot }$]	average $Y_e$
	average $s$ [$k_B/b$]	at $T=5GK$
	average $v$ [$c$]
	at $800r_g$
HS-Therm	0.00043	0.34	17.57	0.151
HS-Magn	0.00386	0.28	14.56	0.229
LS-Therm	0.00012	0.26	12.83	0.110
LS-Magn	0.00315	0.22	12.89	0.177

Equations

$${T_{\left(m\right)}}^{\mu \nu }=\rho \xi u^{\mu }{u}^{\nu }+p{g}^{\mu \nu }$$

$${T_{\left(em\right)}}^{\mu \nu }=b^{\kappa }{b}_{\kappa }{u}^{\mu }{u}^{\nu }+\frac{1}{2}{b}^{\kappa }{b}_{\kappa }{g}^{\mu \nu }-b^{\mu }{b}^{\nu }$$

$$T^{\mu \nu }={T_{\left(m\right)}}^{\mu \nu }+{T_{\left(em\right)}}^{\mu \nu }$$

$$(\rho u_{\mu }{)}_{;\nu }=0$$

$$T_{\nu ;\mu }^{\mu }=0.$$

$$r=R_0+\exp \left[x^{(1)}\right]$$

$$\theta =\frac{\pi }{2}\left(1+x^{(2)}\right)+\frac{1-h}{2}\sin \left[\pi \left(1+x^{(2)}\right)\right]$$

$$A_{\phi }=\frac{\bar{\rho }}{{\rho }_{\max }}-{\rho }_0$$

$$L_{\mathrm{unit}}=\frac{GM}{c^2}=1.48\cdot 10^5\frac{M}{M_{\odot }}\mathrm{cm}$$

$$T_{\mathrm{unit}}=\frac{r_g}{c}=4.9\cdot 10^{-6}\frac{M}{M_{\odot }}\mathrm{s}$$

$$n_{\mathrm{p}}({\Gamma }_{p+e^{-}\to n+{\nu }_e}+{\Gamma }_{p+{\bar{\nu }}_e\to n+e^{+}}+{\Gamma }_{p+e^{-}+{\bar{\nu }}_e\to n})=n_{\mathrm{n}}({\Gamma }_{n+e^{+}\to p+{\bar{\nu }}_e}+{\Gamma }_{n+{\nu }_e\to p+e^{-}}+{\Gamma }_{n\to p+e^{-}+{\bar{\nu }}_e}).$$

$$n_{\mathrm{i}}=\frac{\sqrt{2}}{{\pi }^2}\frac{(m_i{c}^2{)}^3}{(\hslash c{)}^3}{\beta }_i^{3/2}\left[F_{1/2}({\eta }_{\mathrm{i}},{\beta }_{\mathrm{i}})+\frac{1}{2}{\beta }_{\mathrm{i}}{F}_{3/2}({\eta }_{\mathrm{i}},{\beta }_{\mathrm{i}})\right]$$

$$Y_{\mathrm{e}}=\frac{n_{{\mathrm{e}}^{-}}-n_{{\mathrm{e}}^{+}}}{n_{\mathrm{b}}}.$$

$$P=P_{\mathrm{nucl}}+P_{\mathrm{He}}+P_{\mathrm{rad}}+P_{\nu }$$

$$P_{\mathrm{i}}=\frac{2\sqrt{2}}{3{\pi }^2}\frac{(m_i{c}^2{)}^4}{(\hslash c{)}^3}{\beta }_i^{5/2}\left[F_{3/2}({\eta }_{\mathrm{i}},{\beta }_{\mathrm{i}})+\frac{1}{2}{\beta }_{\mathrm{i}}{F}_{5/2}({\eta }_{\mathrm{i}},{\beta }_{\mathrm{i}})\right]$$

$$P_{\nu }=\frac{7}{8}\frac{{\pi }^2}{15}\frac{(kT{)}^4}{3(\hslash c{)}^3}\sum _{i=e,\mu ,\tau }\frac{0.5({\tau }_{a,i}+{\tau }_s)+\frac{1}{\sqrt{3}}}{0.5({\tau }_{a,i}+{\tau }_s)+\frac{1}{\sqrt{3}}+\frac{1}{3{\tau }_{a,i}}}$$

$$Q_{\nu }=\frac{(7/8)\sigma T^4}{(3/4)}\sum _{i=e,\mu ,\tau }\frac{1}{0.5({\tau }_{a,i}+{\tau }_s)+\frac{1}{\sqrt{3}}+\frac{1}{3{\tau }_{a,i}}}$$

$${\lambda }_{\mathrm{MRI}}=\frac{2\pi }{\Omega }\frac{b}{\sqrt{4\pi \rho h+b^2}}$$

Full text

accretion, accretion disks – black hole physics – gamma-ray burst: general – winds, outflows – MHD –nuclear reactions

1 Introduction

The compact object binaries, namely the Black Hole-Neutron Star (BHNS) and the Neutron Star-Neutron Star (NSNS) binaries (Eichler et al., 1989; Paczynski, 1991; Narayan et al., 1992), are favorable progenitors for the Short Gamma Ray Bursts (sGRB) (see however Berger (2011) and references therein for alternative explanations, including magnetars formed through a massive star core-collapse, or binary white dwarf mergers, or the white dwarf accretion-induced collapse, or accretion-induced collapse of neutron stars). The complex nature of macroscopic and microphysical properties of the central engines requires a still ongoing effort to identify crucial aspects of their operation. Some of the fundamental requirements for the mechanism responsible for the outflow launching were already stated years ago. Apart from the collimated jets (Sari et al., 1999; Rhoads, 1999), these engines are also supposed to produce the uncollimated, equatorial wind outflows, mediated by magnetic fields and centrifugal forces (Narayan et al., 2012; Fernández et al., 2015).

A major breakthrough occurred recently with the multi-messenger detection of the GRB 170817A through the VIRGO/LIGO and Fermi-GBM observatories. The waveform of the emitted gravitational waves is consistent with the merging of an NSNS system (Abbott et al., 2017). The event was followed by a sGRB after $\sim 1.7s$, as given by the Fermi Gamma-ray Space Telescope trigger. The associated GRB was a few orders of magnitude fainter than a typical short burst, and its energy $E_{\rm iso}=\left(3.1\pm 0.7\right)\cdot 10^{46}$ ergs was explained as the off jet-axis observation and the line of sight passing through the surrounding cocoon (Lazzati et al., 2017), or an intrinsic property of the specific NSNS merging event (Murguia-Berthier et al., 2017; Zhang et al., 2017). The optical counterpart of the GRB was discovered and identification of NGC 4993 as its host was reported by Coulter et al. (2017).

On the outcome of the merging process, Granot et al. (2017) gave a comprehensive discussion. Among the potential remnants, a massive NS and the long lasting, differentially rotating supra-massive neutron star (SMNS) do not lead to the formation of a BH-torus system that powers the bursts (Shibata et al., 2000; Margalit et al., 2015). However, a hyper-massive neutron star (HMNS) may also form a BH-torus engine. In case of the GW 170817 event, a delayed collapse of a HMNS might be the reason for the observed time delay between the merger and GRB phenomena (Granot et al., 2017). In the same event, the early optical and infrared emission due to r-process is also in favor of the HMNS scenario (Margalit & Metzger, 2017). Such emission has been confirmed by a number of ground based observations, e.g., Smartt et al. (2017).

On the theoretical ground, it was proposed already e.g. by Li & Paczyński (1998) that compact binary mergers eject a small fraction of matter with a subrelativistic velocities. This medium condenses into neutron-rich nuclei, most of which are radioactive and provide a long-term heat source for the expanding envelope. The first tentative signal of this kind was reported in 2013, when the ground-based optical and Hubble Space Telescope optical and near-IR observations of the short-hard GRB 130603B revealed the presence of near-IR emission. It was explained as the effect of an r-process powered transient (Berger et al., 2013; Tanvir et al., 2013).

This emission may originate from the long tidal tails of coalescing compact binaries, as proposed by the hydrodynamic and nuclear reaction network calculations (see, e.g., Roberts et al. (2011)). In addition, the r-process nucleosynthesis may be contributed by the black hole accretion disk outflows. Such studies were performed in the past, e.g. in the frame of semi-analytic, time-dependent evolution models (Fujimoto et al., 2004; Wanajo & Janka, 2012), as well as in numerical simulations in the pseudo-Newtonian gravity (Just et al., 2015; Wu et al., 2016).

In the present work, we present a numerical simulation of the short GRB central engine as composed of the rotating black hole, surrounded by a remnant torus which has already formed after the NSNS binary coalescence. We use the full general relativistic framework and we present an axisymmetric simulation in the fixed Kerr background metric. The equation of state (EOS) of the dense matter in the torus describes the Fermi gas where the gas pressure is contributed by partially degenerate nucleons, electron positron pairs, and Helium nuclei. The weak interactions are controlled by the nuclear equilibrium conditions and establish the neutronisation level in the torus plane, as well as in the outflows launched from its surface. The torus matter is magnetized and its accurate evolution is followed when the MRI turbulence (Balbus & Hawley, 1991) drives the mass inflow. In addition, the magnetic field is responsible for launching the uncollimated outflows of the plasma, while the rotation of the black hole powers the polar jets. Our simulations include also the neutrino emission, as in the framework of the so-called neutrino-dominated accretion flows, NDAF (Kohri et al., 2005).

Our numerical scheme is based on the GRMHD code HARM (Gammie et al., 2003; Noble et al., 2006), but equipped with the non-adiabatic EOS. It is embedded in the relativistic conservative scheme and interpolation over the temperatures and densities spans several orders of magnitude, as first proposed by Janiuk (2017). The main focus, and novelty of the present study is the imprint of the torus properties on the amount and chemical composition of the emerging ejecta. We follow the wind outflow, and compute the synthesis of subsequent r-process elements on the trajectories, where mass is ejected in sub-relativistic particles. The effectiveness of the adopted method is that it distributes tracers uniformly in rest-mass density in flow, which is a non-trivial task in numerical relativity (Bovard & Rezzolla, 2017).

The current paper is an extension of our previous works (Janiuk et al., 2013; Janiuk, 2017), where we studied the neutrino cooling and microphysical properties of the GRB central engine. In the former studies, the accretion was also modeled with general relativistic MHD simulations, but we did not follow the dynamics of the outflows, and we implemented only the nuclear reaction networks producing heavy ions under the nuclear statistical equilibrium (Hix & Meyer, 2006). In the current approach, we follow the neutron rich ejecta, where the synthesis of isotopes proceeds faster than the equilibrium timescale, and we are able to obtain the abundance patters reaching the third peak ($A\sim 195$). The effective postprocessing is made with the modular reaction network library (Lippuner & Roberts, 2017).

The article is organized as follows. In Section 2 we present the simulation scheme and the initial configuration of our model. We also describe the properties of the EOS (Section 2.2), and the concept of tracer particles (Section 2.3). The results showing the general properties of the flow, and the nuclear reaction network simulation outcome, appear in Section 3. Discussion and conclusions are the subject of Section 4.

2 The Simulations Setup

We use the general relativistic magneto-hydrodynamic code, HARM (Gammie et al., 2003; Noble et al., 2006), to integrate our model under a fixed Kerr metric, i.e. neglecting effects like the self gravity of the disrupted material, or the BH spin changes. The HARM code is a finite volume, shock capturing scheme that solves the hyperbolic system of the partial differential equations of GR MHD. The numerical scheme is based on GR MHD equations with the plasma energy-momentum tensor, $T^{\mu\nu}$, with contributions from gas and electromagnetic field

[TABLE]

where $u^{\mu}$ is the four-velocity of gas, $u$ denotes internal energy density, p is pressure, $b^{\mu}$ is magnetic four-vector, and $\xi$ is the fluid specific enthalpy, $\xi=(\rho+p+u)/\rho$. The continuity and momentum conservation equation reads:

[TABLE]

They are brought in conservative form, by implementing the Harten, Lax, van Leer (HLL) solver to calculate numerically the corresponding fluxes.

In terms of the Boyer-Lindquist coordinates, $\left(r,\theta,\phi\right)$, the black hole is located at $0<r\leq r_{\rm h}$, where $r_{\rm h}=\left(1+\sqrt{1-a^{2}}\right)r_{\rm g}$ is the horizon radius of a rotating black hole with mass $M$ and angular momentum $J$ in geometrized units, $r_{\rm g}=GM/c^{2}$, and $a$ is the dimensionless Kerr parameter, $a=J/(Mc),0\leq a\leq 1$. In our simulations we investigate the rotating black holes, $a=0.6-0.9$.

The HARM code doesn’t perform the integration in the Boyer-Lindquist coordinates, but instead in the so called Modified Kerr-Schild ones: $\left(t,x^{(1)},x^{(2)},\phi\right)$ (Noble et al., 2006). The transformation between the coordinate systems is given by:

[TABLE]

where $R_{0}$ is the innermost radial distance of the grid, $0\leq x^{(2)}\leq 1$, and $h$ is a parameter that determines the concentration of points at the mid-plane. In our models we use $h=0.3$ (notice that for $h=1$ and a uniform grid on $x^{(2)}$ we obtain an equally spaced grid on $\theta$, while for $h=1$ the points concentrate on the mid plane). The exponential resolution in the $r$-direction leads to higher resolution and it is adjusted to resolve the initial propagation of the outflow. Our grid resolution is $256\times 256$.

2.1 The Torus Initial Configuration

The accreting material is modeled following Fishbone & Moncrief (1976) (hereafter FM) who provided an analytic solution of a constant specific angular momentum, in a steady state configuration of a pressure-supported ideal fluid in the Kerr black hole potential. Other similar configurations like Chakrabarti (1985) with a power law radial evolution of the angular momentum or of independently varying the Bernoulli parameter (sum of the kinetic and potential energy, and enthalpy of the gas) and disk thickness are also possible (Penna et al., 2013).

In the FM model, the position of the material reservoir is determined by the radial distance of the innermost cusp of the torus, $r_{\rm in}$, and the distance where the maximum pressure occurs, $r_{\rm max}$. Because of its geometry, the relative difference of the two radii determines also the dimension of the torus, with higher differences resulting to extended cross section. Subsequently the $r_{\rm in}$ and $r_{\rm max}$ determine also the angular momentum value and the distribution of the angular velocity along the torus.

The initial torus is embedded in a poloidal magnetic field, prescribed with the vector potential of

[TABLE]

where $\bar{\rho}$ is the average density in the torus, $\rho_{\rm max}$ is the density maximum, and we use $\rho_{0}=0.2$. As a consequence, the magnetized flows considered here are not in equilibrium, and the angular momentum is transported. Regardless of the specific magnitude of the plasma $\beta-$parameter the flow slowly relaxes its initial configuration, becomes geometrically thinner, and launches the outflows.

We examine two sets of models with different initial $\beta-$parameter, defined as the ratio of the fluid’s thermal to the magnetic pressure, $\beta\equiv p_{g}/p_{\rm mag}$ for the torus configurations, while every set includes models differing with the black hole spin.

Following Janiuk et al. (2013), we use the value of the total initial mass of the torus to scale the density over the integration space, and we base our simulations on the physical units. We use

[TABLE]

as the spatial and time units, respectively. Now, the density scale is related to the spatial unit as $D_{\rm unit}=M_{\rm scale}/L_{\rm unit}^{3}$. We conveniently adopt the scaling factor of $M_{\rm scale}=1.5\times 10^{-5}M_{\odot}$, so that the total mass contained within the simulation volume (mainly in the FM torus) is about 0.1 $M_{\odot}$, for a black hole mass of 3 $M_{\odot}$. The $r_{\rm in}$ and $r_{\rm max}$ radii of the torus are having the fiducial values of $r_{\rm in}=(3.5-4.5)\,r_{g}$ and $r_{\rm max}=(9-11)\,r_{g}$ (see Table 3).

To sum up, the GRB engines are modeled with the following global parameters: mass of the black hole $M_{\rm BH}$, the torus mass, and black hole dimensionless spin, $a$. Typical parameters are $M_{\rm BH}=3\,M_{\odot}$, and $a=0.6-0.9$, while the torus mass is about 0.1 $M_{\odot}$, as resulting from its size in geometrical units, and adopted density scaling. The accretion rate onto the black hole, measured as the mass flux transported through the horizon, is varying with time. Its mean value expressed in physical units is on the order of 0.1 $M_{\odot}$s*-1*, as converted from the density scaling and time unit.

2.2 Equation of state and neutrino cooling

We consider the torus composed of free protons, neutrons, electron-positron pairs, and Helium nuclei. For a given baryon number density and temperature, the equilibrium condition is assumed between the reactions of electron-positron capture on nucleons, and neutron decays (Reddy et al., 1998).

Namely, the ratio of free protons, is determined from the equilibrium between the transition reactions from neutrons to protons, and from protons to neutrons: $p+e^{-}\rightarrow n+\nu_{e}$, $p+\bar{\nu}_{e}\rightarrow n+e^{+}$, $p+e^{-}+\bar{\nu}_{e}\rightarrow n$, $n+e^{+}\rightarrow p+\bar{\nu}_{e}$, $n+\nu_{e}\rightarrow p+e^{-}$, and $n\rightarrow p+e^{-}+\bar{\nu}_{e}$. The balance equation has a form:

[TABLE]

The above reaction rates are the sum of forward and backward rates and the appropriate formulae for $\Gamma$’s are given in Kohri et al. (2005).

The number density of fermions under arbitrary degeneracy is determined by

[TABLE]

with $F_{\rm k}$ being the Fermi-Dirac integrals of the order $k$, and $\beta_{i}=kT_{i}/(m_{i}c^{2})$.

This is supplemented with the charge neutrality condition, $n_{e^{-}}=n_{e^{+}}+n_{p}+n^{0}$, where $n^{0}=2n_{\rm He}=(1-X_{\rm nuc})n_{b}/2$ is the number of protons in Helium, and the conservation of baryon number, $n_{\rm p}+n_{\rm n}=X_{\rm nuc}n_{\rm b}$. The fraction of free nuclei is a strong function of density and temperature, and is adopted from the fitting formula after Qian & Woosley (1996).

We can define the proton-to-baryon number density ratio $Y_{\rm p}$, and the electron fraction:

[TABLE]

The free species may have an arbitrary degeneracy level. The total pressure is contributed by free nucleons, pairs, radiation, alpha particles, and also trapped neutrinos (Yuan, 2005; Janiuk et al., 2007), and is given by:

[TABLE]

where $P_{\rm nucl}=P_{e^{-}}+P_{e^{+}}+P_{n}+P_{p}$** ** and

[TABLE]

where $\eta_{\rm e}$, $\eta_{\rm p}$ and $\eta_{\rm n}$ are the reduced chemical potentials. Here, $\eta_{i}=\mu_{i}/kT$ is the degeneracy parameter (where $\mu_{i}$ is the standard chemical potential). Reduced chemical potential of positrons is $\eta_{\rm e+}=-\eta_{\rm e}-2/\beta_{\rm e}$.

The total pressure of subnuclear matter is mainly contributed by electrons, and therefore it is influenced by the changes in the electron fraction. Reduced electron chemical potentials are somewhat larger than unity, because at these accretion rates which we consider here, electrons are slightly degenerate. The pressure of Helium is taken to be of ideal gas, $P_{\rm He}=(n_{b}/4)(1-X_{\rm nuc})kT$. The radiation pressure, $P_{\rm rad}=(4\sigma)/(3c)T^{4}$, is in GRB disks a few orders of magnitude smaller than other components.

2.3 Neutrino cooling

The reactions of the electron and positron capture on nucleons (a.k.a. URCA reactions, see above), and also the electron-positron pair annihillation, $e^{+}+e^{-}\rightarrow\nu_{i}+\bar{\nu}_{i}$, nucleon bremsstrahlung, $n+n\rightarrow n+n+\nu_{i}+\bar{\nu}_{i}$, and plasmon decay, $\tilde{\gamma}\rightarrow\nu_{e}+\bar{\nu}_{e}$, are producing neutrinos, which act as a cooling mechanism for the plasma.

The cooling rates due to the bremsstrahlung and plasmon decay, are $q_{\rm brems}=3.35\cdot 10^{27}\rho_{10}^{2}T_{11}^{5.5}$, and $q_{\rm plasm}=1.5\cdot 10^{32}T_{11}^{9}\gamma_{p}^{6}e^{-\gamma_{p}}(1+\gamma_{p})(2+\gamma_{p}^{2}/(1+\gamma_{p})$ where $\gamma_{p}=5.56\cdot 10^{-2}\sqrt{(\pi^{2}+3\eta_{e}^{2})/3}$ (Ruffert et al., 1996).

The cooling rates due to URCA processes, $q_{\rm urca}$, and pair annihillation, $q_{\rm pair}$, reactions are having more complex forms, and can be found i.e. in Janiuk et al. (2007) (see equations A7-A15 therein). They involve the distribution functions, and the blocking factors, which describe the extent on which neutrinos are trapped (see below). In the GRB accretion disk, we consider the neutrino transparent and opaque regions, and transition between the two. In terms of the blocking factors, $0\leq b_{i}\leq 1$, the neutrinos of $\nu_{i}$ flavor might be freely escaping or trapped, consistently with the two-stream approximation (Di Matteo et al., 2002). The blocking factor of trapped neutrinos is however used only for the URCA emissivities. For the annihillation reaction it is not used, because these emissivities are much smaller, and do not change the electron fraction.

Trapped neutrinos give a contribution to the pressure with a component

[TABLE]

where $\tau_{a,i}$ and $\tau_{s}$ denote absorption and scattering. Here $\tau_{a,i}=(H/((7/8)\sigma T^{4}))q_{a,\nu_{i}}$, where $q_{a,\nu_{e}}=q_{pair}+q_{urca}+q_{plasm}+{1\over 3}q_{brems}$ and $q_{a,\nu_{\mu}}=q_{pair}+{1\over 3}q_{brems}$. The scale $H$ is the local thickness of the disk given by the pressure scaleheight. The free escape of neutrinos is further limited by their scattering on neutrons and protons, for which we use a formula $\tau_{\rm s}=24.3\cdot 10^{-5}((kT/m_{e}c^{2}))^{2}H(C_{s,p}n_{p}+C_{s,n}n_{n})$, with $(C_{s,p}=[4(C_{V}-1)^{2}+5\alpha^{2}]/24$ and $(C_{s,n}=[1+5\alpha^{2}]/24$, $C_{V}=1/2+2\sin^{2}\theta_{C}$, and $\sin^{2}\theta_{\rm C}=0.23$. The ’blocking factor’ for electron neutrino is then defined as $b_{e}=((\tau_{a,e}+\tau_{s})/2+1/\sqrt{3})/(\tau_{a,e}/2+1/\sqrt{3}+1/(3\tau_{a,e})$.

Neutrino cooling rate is finally computed as

[TABLE]

in the optically thick regime. In the optically thin regime, it will be given by $Q_{\nu}=H(q_{\rm pair}+q_{\rm urca}+q_{\rm plasm}+q_{\rm brems})$.

We store the neutrino cooling rate as a function of the rest mass density and temperature, and we also store the corresponding electron fraction and total pressure values, that result from the EOS. The numerically computed EOS is tabulated and the pressure that is used in the MHD evolution, comes out from interpolation over the density and temperature over the grid, by means of the spline method (Akima, 1970). The GR MHD computation in this case is non-trivial, and also technically demanding, as the numerical scheme solves at every time-step for the inversion between the so called ’primitive’ and ’conserved’ variables, the latter being the total energy and comoving density (see e.g., Noble et al. (2006)). Our modified HARM scheme, introduced previously in Janiuk (2017), takes into account the total pressure as given by Eq. 6, and also its derivatives over density and energy, when computing the sound speed square. Hence the magneto-sonic velocity, defined as in (Ibáñez et al., 2015) is then properly adopted by the MHD stress tensor and source terms in the equations of motion. Notice that in the pioneering works dealing with GR MHD disks (McKinney et al., 2012), the adiabatic equation of state was used in the form of $p=(\gamma_{\rm ad}-1)u$ (where $\gamma_{\rm ad}$ is an adiabatic index) during the dynamical simulations because they were addressed to much lower temperature and density systems, such as AGN disks.

2.4 The use of tracer particles

We define the tracers already while initializing the GR MHD simulation. In every grid cell we check first, if the density is larger than some minimum value (e.g., $10^{-4}$ of the maximum density in the torus), hence we pick all the particles which are embedded in the densest parts of the accreting flow. During the evolution, we update the positions of tracers according to their new (contravariant) velocity components. The linear interpolation of the four-velocities is performed, so that we are always in the grid cell centers.

While tracking the moving particle, we record their coordinates, density, temperature, and electron fraction, which will be used later to compute the nucleosynthesis. This is done always, whenever the ultimate fate of the trajectory is determined. The particles can escape from the computational domain either through the inner boundary, $R_{\rm in}$, or through the outer boundary, $R_{\rm out}$. While the inflow corresponds to the black hole accretion, in this study we are interested mostly in the outflow properties. Figure 1 shows the trajectories of tracer particles which left through $R_{\rm out}=1000~{}r_{\rm g}$ as computed during the evolution of the torus until time $t_{f}=20000~{}M$. We present four models, that are listed in Table 3. In Fig. 1, we plotted only the uncollimated outflows. The collimated ones (i.e. the “jets”), and ignored. These traces are defined as having the polar angles less than a minimum value (here: $\theta_{\rm min}=0.02\pi$), as measured from either of the polar semi-axes.

The trajectories recorded in such way are further post-processed to compute the $r$-process nucleosynthesis. Given the evolution of density, temperature, and electron fraction, they provide input for the thermonuclear reaction network.

3 Results