The r-process Nucleosynthesis in the Outflows from Short GRB Accretion Disks

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_in, and the distance where the maximum pressure occurs, r_max. Because of its geometry, the relative difference of the two radii determines also the dimension of the torus, with higher differences resulting in an extended cross section. Subsequently, the r_in and r_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

$$\begin{eqnarray}&&{A}_{\varphi }=\displaystyle \frac{\bar{\rho }}{{\rho }_{\max }}-{\rho }_{0},\end{eqnarray} \tag{ 1 }$$

where $\bar{\rho }$ is the average density in the torus, ρ_max is the density maximum, and we use ρ₀ = 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 β-parameter the flow slowly relaxes its initial configuration, becomes geometrically thinner, and launches the outflows.

We examine two sets of models with different initial β-parameter, defined as the ratio of the fluid's thermal to the magnetic pressure, $\beta \equiv {p}_{g}/{p}_{\mathrm{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

$$\begin{eqnarray*}\begin{array}{rcl}{L}_{\mathrm{unit}} & = & \displaystyle \frac{{GM}}{{c}^{2}}=1.48\cdot {10}^{5}\displaystyle \frac{M}{{M}_{\odot }}\,\,\mathrm{cm}\\ {T}_{\mathrm{unit}} & = & \displaystyle \frac{{r}_{g}}{c}=4.9\cdot {10}^{-6}\displaystyle \frac{M}{{M}_{\odot }}\,\,{\rm{s}}\end{array}\end{eqnarray*}$$

as the spatial and time units, respectively. Now, the density scale is related to the spatial unit as ${D}_{\mathrm{unit}}={M}_{\mathrm{scale}}/{L}_{\mathrm{unit}}^{3}$. We conveniently adopt the scaling factor of M_scale = 1.5 × 10⁻⁵ M_⊙, so that the total mass contained within the simulation volume (mainly in the FM torus) is about 0.1 M_⊙, for a black hole mass of 3 M_⊙. The r_in and r_max radii of the torus have fiducial values of r_in = (3.5–4.5)r_g and r_max = (9–11)r_g (see Table 1).

Table 1.  Summary of the Models

Model	Torus Mass		BH Spin	rin	rmax	β0	${\dot{M}}_{\mathrm{out}}$	ΔMout
	M⊙ units		a	rg units	rg units		M⊙ s−1 units	M⊙ units
	M0	Mend
HS-Therm	0.1031	0.0848	0.9	3.8	9.75	100	7.64 · 10−2	4.42 · 10−3
HS-Magn	0.1031	0.0741	0.9	3.8	9.75	10	8.26 · 10−2	1.72 · 10−2
LS-Therm	0.1104	0.0895	0.6	4.5	9.1	100	8.76 · 10−2	2.29 · 10−3
LS-Magn	0.1104	0.0682	0.6	4.5	9.1	10	1.20 · 10−1	1.81 · 10−2

Note. The first two models refer to a highly spinning black hole in a short GRB engine, and the last two models represent a moderately spinning black hole. The inner radius of the torus, r_in, the radius of pressure maximum, r_max, and the plasma-β are given as the initial state parameters. The last two columns give the mass loss rate through the outer boundary, averaged over the simulation time, and the cumulative mass lost in the outflow.

Download table as:  ASCIITypeset image

To sum up, the GRB engines are modeled with the following global parameters: mass of the black hole M_BH, the torus mass, and black hole dimensionless spin, a. Typical parameters are M_BH = 3 M_⊙, and a = 0.6–0.9, while the torus mass is about 0.1 M_⊙, 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_⊙ s⁻¹, as converted from the density scaling and time unit.

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}^{-}\to n+{\nu }_{e}$, $p+{\bar{\nu }}_{e}\to n+{e}^{+}$, $p+{e}^{-}+{\bar{\nu }}_{e}\to n$, $n+{e}^{+}\to p+{\bar{\nu }}_{e}$, $n+{\nu }_{e}\to p+{e}^{-}$, and $n\to p+{e}^{-}+{\bar{\nu }}_{e}$. The balance equation has a form:

$$\begin{eqnarray}&&\begin{array}{l}{n}_{{\rm{p}}}({{\rm{\Gamma }}}_{p+{e}^{-}\to n+{\nu }_{e}}+{{\rm{\Gamma }}}_{p+{\bar{\nu }}_{e}\to n+{e}^{+}}+{{\rm{\Gamma }}}_{p+{e}^{-}+{\bar{\nu }}_{e}\to n})\\ \quad =\,{n}_{{\rm{n}}}({{\rm{\Gamma }}}_{n+{e}^{+}\to p+{\bar{\nu }}_{e}}+{{\rm{\Gamma }}}_{n+{\nu }_{e}\to p+{e}^{-}}+{{\rm{\Gamma }}}_{n\to p+{e}^{-}+{\bar{\nu }}_{e}}).\end{array}\end{eqnarray} \tag{ 2 }$$

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

The number density of fermions under arbitrary degeneracy is determined by

$$\begin{eqnarray}&&{n}_{{\rm{i}}}=\displaystyle \frac{\sqrt{2}}{{\pi }^{2}}\displaystyle \frac{{({m}_{i}{c}^{2})}^{3}}{{({\hslash }c)}^{3}}{\beta }_{i}^{3/2}\left[{F}_{1/2}({\eta }_{{\rm{i}}},{\beta }_{{\rm{i}}})+\displaystyle \frac{1}{2}{\beta }_{{\rm{i}}}{F}_{3/2}({\eta }_{{\rm{i}}},{\beta }_{{\rm{i}}})\right],\end{eqnarray} \tag{ 3 }$$

with F_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}\,=2{n}_{\mathrm{He}}=(1-{X}_{\mathrm{nuc}}){n}_{b}/2$ is the number of protons in helium, and the conservation of baryon number, n_p + n_n = X_nucn_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_p, and the electron fraction:

$$\begin{eqnarray}&&{Y}_{{\rm{e}}}=\displaystyle \frac{{n}_{{{\rm{e}}}^{-}}-{n}_{{{\rm{e}}}^{+}}}{{n}_{{\rm{b}}}}.\end{eqnarray} \tag{ 4 }$$

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:

$$\begin{eqnarray}&&P={P}_{\mathrm{nucl}}+{P}_{\mathrm{He}}+{P}_{\mathrm{rad}}+{P}_{\nu }\end{eqnarray} \tag{ 5 }$$

where ${P}_{\mathrm{nucl}}={P}_{{e}^{-}}+{P}_{{e}^{+}}+{P}_{n}+{P}_{p}$ and

$$\begin{eqnarray}&&{P}_{{\rm{i}}}=\displaystyle \frac{2\sqrt{2}}{3{\pi }^{2}}\displaystyle \frac{{({m}_{i}{c}^{2})}^{4}}{{({\hslash }c)}^{3}}{\beta }_{i}^{5/2}\left[{F}_{3/2}({\eta }_{{\rm{i}}},{\beta }_{{\rm{i}}})+\displaystyle \frac{1}{2}{\beta }_{{\rm{i}}}{F}_{5/2}({\eta }_{{\rm{i}}},{\beta }_{{\rm{i}}})\right],\end{eqnarray} \tag{ 6 }$$

where η_e, η_p, and η_n are the reduced chemical potentials. Here, η_i = μ_i/kT is the degeneracy parameter (where μ_i is the standard chemical potential). The 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 changes in the electron fraction. Reduced electron chemical potentials are somewhat larger than unity, because at the accretion rates that we consider here electrons are slightly degenerate. The pressure of helium is taken to be of ideal gas, ${P}_{\mathrm{He}}=({n}_{b}/4)(1-{X}_{\mathrm{nuc}}){kT}$. The radiation pressure, ${P}_{\mathrm{rad}}=(4\sigma )/(3c){T}^{4}$, is in GRB disks a few orders of magnitude smaller than other components.

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}^{-}\to {\nu }_{i}+{\bar{\nu }}_{i}$, nucleon bremsstrahlung, $n+n\to n+n+{\nu }_{i}+{\bar{\nu }}_{i}$, and plasmon decay, $\tilde{\gamma }\to {\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}_{\mathrm{brems}}=3.35\cdot {10}^{27}{\rho }_{10}^{2}{T}_{11}^{5.5}$, and ${q}_{\mathrm{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_urca, and pair annihillation, q_pair, reactions have 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 to 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 ≤ b_i ≤ 1, the neutrinos of ν_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

$$\begin{eqnarray}&&{P}_{\nu }=\displaystyle \frac{7}{8}\displaystyle \frac{{\pi }^{2}}{15}\displaystyle \frac{{({kT})}^{4}}{3{({\hslash }c)}^{3}}\displaystyle \sum _{i=e,\mu ,\tau }\displaystyle \frac{0.5({\tau }_{a,i}+{\tau }_{s})+\tfrac{1}{\sqrt{3}}}{0.5({\tau }_{a,i}+{\tau }_{s})+\tfrac{1}{\sqrt{3}}+\tfrac{1}{3{\tau }_{a,i}}},\end{eqnarray} \tag{ 7 }$$

where τ_a,i and τ_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}_{\mathrm{pair}}+{q}_{\mathrm{urca}}\,+{q}_{\mathrm{plasm}}+\displaystyle \frac{1}{3}{q}_{\mathrm{brems}}$ and ${q}_{a,{\nu }_{\mu }}={q}_{\mathrm{pair}}+\displaystyle \frac{1}{3}{q}_{\mathrm{brems}}$. The scale H is the local thickness of the disk given by the pressure scale height. The free escape of neutrinos is further limited by their scattering on neutrons and protons, for which we use the 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 the 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})$.

The neutrino cooling rate is finally computed as

$$\begin{eqnarray}&&{Q}_{\nu }=\displaystyle \frac{(7/8)\sigma {T}^{4}}{(3/4)}\displaystyle \sum _{i=e,\mu ,\tau }\displaystyle \frac{1}{0.5({\tau }_{a,i}+{\tau }_{s})+\tfrac{1}{\sqrt{3}}+\tfrac{1}{3{\tau }_{a,i}}}\end{eqnarray} \tag{ 8 }$$

in the optically thick regime. In the optically thin regime, it will be given by Q_ν = H(q_pair + q_urca + q_plasm + q_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 is revealed 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 nontrivial, 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 Equation (6), and also its derivatives over density and energy, when computing the sound speed square. Hence the magnetosonic 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 EOS was used in the form of $p=({\gamma }_{\mathrm{ad}}-1)u$ (where γ_ad is an adiabatic index) during the dynamical simulations because they were addressed to much lower temperature and density systems, such as AGN disks.

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⁻⁴ of the maximum density in the torus), hence we pick all the particles that 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 always done when the ultimate fate of the trajectory is determined. The particles can escape from the computational domain either through the inner boundary, R_in, or through the outer boundary, R_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 that left through R_out = 1000 r_g as computed during the evolution of the torus until time t_f = 20,000 M. We present four models, that are listed in Table 1. In Figure 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 θ_min = 0.02π), as measured from either of the polar semiaxes.

Zoom In Zoom Out Reset image size

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

Our simulations are divided into two sets with the Magn class of models referring to a weakly magnetized torus and the Therm class to torus with significantly (10 times) higher plasma β-parameter. Table 1 presents a summary of the models and their initial parameters, the total mass of the disk at the initial state, and the averaged mass loss rate through the outer boundary. The physical conditions in the flow, namely its density, temperature, and electron fraction, are governed by the global parameters: accretion rate, black hole mass and spin, and magnetic field normalization.

The overall evolution follows a similar pattern for all models that have been investigated. The FM torus initial state is close to a steady state retaining its structure for a relatively long integration time. Its innermost, densest part is enclosed within ∼50r_g, and forms a narrow cusp through which the matter sinks under the black hole horizon. The surface layers of the torus, which extend to about 200 r_g, form a kind of corona. This region is the base of subrelativistic outflows that start being launched from the center as soon as the initial condition of the simulation is relaxed (this is at about 2000 M, which is 0.03 s for the adopted black hole mass). The models were evolved until the time t_f = 20,000 M (geometrical units, M = G M_BH/c³, make ${t}_{f}=0.3\,{\rm{s}}$ for M_BH = 3M_⊙).

The magnetorotational turbulence is resolved in our simulations with the proper scaling of the cell sizes. We check this by computing the ratio of the wavelength of the fastest growing mode, as given by:

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

We find that our grid always provides at least 10 cells per wavelength, with the MRI being at least marginally resolved inside the torus, and well resolved in the regions of the outflow. In Figure 2 we plot the ratio of the λ_MRI with respect to the local grid resolution, i.e., ${Q}_{\mathrm{MRI}}=\displaystyle \frac{{\lambda }_{\mathrm{MRI}}}{{\rm{\Delta }}\theta }$ (see e.g., Siegel & Metzger 2018), for three representative times in simulation LS-Magn.

Zoom In Zoom Out Reset image size

In Figure 3 we show the profiles of density, magnetic field, neutrino emissivity, and electron fraction, in the r−θ plane, as obtained at the end of each simulation. The models, from top to bottom, represent the cases of LS-Therm (low spin a = 0.6, large gas to magnetic pressure ratio, β = 100), LS-Magn (low spin a = 0.6, smaller gas to magnetic pressure ratio, β = 10), HS-Therm (high spin a = 0.9, large gas to magnetic pressure ratio, β = 100), and HS-Magn (high spin a = 0.9, smaller gas to magnetic pressure ratio, β = 10). As can be seen, the more magnetized models have an effect on the larger extension of the magnetically driven and neutrino-cooled winds at the equatorial plane and intermediate latitudes. These winds are moderately dense (ρ ∼ (10⁵−5) × 10⁸ g cm⁻³) and hot (T ∼ (10⁹−3) × 10¹⁰ K). In addition, there are hot luminous regions near the black hole rotation axis that have a very low baryon density. Their neutrino emissivity is correlated with the black hole spin value (see, e.g., also Caballero et al. 2016 for the neutrino emissivity of the disk dependent on the BH spin).

Zoom In Zoom Out Reset image size

The Y_e value is determined in the grid at every time step (from Equation (4)). The electron fraction is very low (Y_e ∼ 0.1–0.2) only in the very central, innermost parts of the accretion torus. The wind material follows the trajectories that have this initial value of Y_e, and reach the outer boundary during the dynamical simulation. Notice that in the electron fraction maps shown in Figure 3 we use the logarithmic scale in radius, as most of the neutronized matter is located below 100 r_g. In the outflows that reach outer boundary, the thin filaments (visible in the maps of magnetic field distribution, made in the linear scale) drive outwards the neutronized matter. However, the bulk of background material (with very low density) has the final electron fraction of Y_e ∼ 0.5.

The electron fraction in the outflow rises with time (see Figure 3). However, some earlier outflows must have much lower Ye in order to produce the second and third peak elements. The time dependence of the electron fraction is presented in Figure 4, as measured along the tracers and averaged over all angles. As shown in the plot, during the first second the average electron fraction of the ejecta is between 0.2 and 0.35. It then rises up to above Y_e = 0.45, as the outflows expand. The trends observed for the four models considered are such that the more magnetized outflows are on average more neutron rich, for the same black hole spin. Also, this angle-averaged distribution implies that the low spinning black hole produces outflows that are on average more neutron rich.

Zoom In Zoom Out Reset image size

The temperature maps are not shown, but the distribution of temperature is directly followed by the neutrino emissivity (the cooling rates for the neutrino emission by URCA, electron–positron anihillation, bremsstrahlung, and plasmon process, scale with temperature in the power 6, 9, 5.5, and 9, respectively, see Janiuk et al. 2004).

The wind density can be visualized also by means of the outflow tracer's distribution, as was shown in Figure 1. We notice that despite the fact that all these tracers were uniformly distributed in the initial state, their final distribution strongly depends on the model parameters. For a higher BH spin value, much fewer tracer particles are detected in the wind outflows than for a low BH spin, if the flow is more magnetized. In contrast, for the thermally dominated models, the higher BH spin helps launch denser wind outflows. This fact can be understood in terms of competitive action of the MHD turbulence for the wind acceleration, and the Blandford–Znajek process. The latter is strongly dependent on the BH spin value so that for almost maximally rotating black holes the B-Z process will overcome on the uncollimated outflows (see Sapountzis & Janiuk 2019).

In Figure 5 we show the time dependence of the mass loss rate through the outer boundary in the function of time during our simulations. As shown in the figure, the mass outflow is huge at the beginning of the simulation, when the initial condition of the pressure equilibrium torus is being relaxed. However, after some 0.03 s (which corresponds to about 2000 M), the outflow rate saturates at a small value. At the end of the simulation, the outflow rate is below 0.05 M_⊙ s⁻¹. The models with more magnetic pressure result in higher time-averaged mass loss rates (see Table 1). Also, the total mass lost from the outer boundary follows the same trend with magnetization. For large BH spin of a = 0.9, the mass outflow rate in the second half of the simulation is larger in the less magnetized model, while initially it was smaller. The effect seen in the figure might be partially an artifact of the initial condition, and the MRI turbulence decay at late times. We notice that the total mass that is lost from the disk during our simulations is between 17% and 38% of the initial disk mass (see Table 1). This number, however, takes into account both the unbound outflows and mass accreted through the BH horizon, as well as the polar jets (albeit these are of a very low density). The mass lost through the outer boundary appears to be in the range of 2%–16% only, in contrast to the results of Fernández et al. (2019). The instantaneous mass of the unbound outflowing material, as estimated via sampling the tracers, is also consistent with the above results and with the intuitive prediction that the larger magnetic field helps drive larger mass of the outflows (see Table 2). Finally, we checked for the amount of unbound matter with −hu_t > 1, i.e., the condition corresponding to a positive Bernoulli parameter in Newtonian gravity. It results in the potential mass of the outflows in the range between 10⁻³ and 10⁻² M_⊙, varying with time in the simulation, with higher values obtained for more magnetized models.

Zoom In Zoom Out Reset image size

We use the data from our simulations, and the particle trajectories, as an input to the nuclear reaction network. These computations are performed using the code SkyNet (Lippuner & Roberts 2017), in the postprocessing simulation. The code is capable of tracing nucleosynthesis in the rapid neutron capture process and involves a large database of over a thousand isotopes. It takes into account the fission reactions and electron screening. In our SkyNet runs, the Helmholtz EOS is used by this code, and it calls for the weak and strong reaction libraries, and for the spontaneous and symmetric fission. Self-heating is taken into account, while the screening is neglected.

The results of MHD simulations that are stored in the tracer particles data are the density, temperature, and electron fraction. In the postprocessing, the density distribution is followed along each of the particle trajectories, and the piecewise linear function is used for interpolation of the densities. When the time is exceeding the original simulation timescale, the density profile is extrapolated with a power-law continuation of t^−3.0, in good agreement with the homologous expansion of the outflow. Because the nuclear self-heating is taken into account, we take only the initial value of temperature on each trajectory, and then the temperature adjusts itself to the reactions balance. We checked that the average temperature in the outflows rises during the first ∼10 s, and then drops to below 10⁶ K after several hundred seconds from the outflow launching. The electron fraction value is also read from each of the trajectories initially, to establish the conditions for nuclear statistical equilibrium and initial chemical abundance pattern. When following the outflow, the electron fraction is updated according to the nuclear reactions balance. The values less than Y_e ∼ 0.05–0.3 are found in the tracers below 1 s of expansion. After a hundred seconds, the electron fraction level in the wind saturates, and stays around Y_e ∼ 0.4–0.5, depending on the angle.

To probe the nucleosynthesis process in our simulations, we investigated the thermodynamic properties of the outflowing ejecta. In Figure 6 we present histograms of the electron fraction, entropy, and velocity, as distributed according to the mass carried in the outflows, showing how much mass within the outflows carries these quantities of certain values. The distributions are plotted at the time when the outflow temperature is still large, and equal to 5 GK. As shown in the first panel of this figure, the largest mass of the outflow with smallest electron fraction, Y_e < 0.2, is launched in the simulations with stronger magnetic field. Here, the comparison between HS-Magn, and LS-Magn histograms, reveals the additional influence of the black hole spin on the results. The Y_e distribution of highly neutronized outflows is narrower for the spin a = 0.6, than for a = 0.9. It means that quickly rotating black holes tend to launch a slightly smaller amount of the most neutron-rich outflows, while the overall mass of the outflow with Y_e > 0.2 remains similar. These outflows have also a broader distribution of entropy and velocity (middle and right panels).

Zoom In Zoom Out Reset image size

For simulations with the weak magnetization, both the total mass of the outflow and the fraction of mass with Y_e < 0.2 are very small, and they are smaller for the lower black hole spin. The specific entropy in these outflows is concentrated around 10 k_B/b, and their velocity does not exceed 0.3 c at the distance 800 r_g. The mean values of the outflow velocity, electron fraction, and entropy, are summarized in Table 2.

In Figure 7 we show the final results of r-process nucleosynthesis in our simulated black hole accretion disk outflows. Shown are four models, for two values of the black hole spin, a = 0.6, and two values of the gas-to-magnetic pressure ratio, β = 100 (left panel) and β = 10 (right panel). The results for the relative abundances of created isotopes are obtained here after extrapolation of the wind outflow up to the time of 1 Myr, and compared with the solar abundance data taken from Arnould et al. (2007).

Zoom In Zoom Out Reset image size

We note that the second and third peaks of the r-process element distribution are well reproduced in our simulated profiles. The relative abundances in the first peak are, however, somewhat smaller than observed. This is found at most trajectories (gray lines) and in the profiles averaged over all trajectories (red line). This would suggest that the iron group isotopes can be synthesized in the short GRB accretion disk outflows, but only at moderate rates.

In the network simulations, we can trace the angular distribution of the outflow trajectories, on which the elements of subsequent abundance peaks are produced. We found that the most efficient production of the second and third peak elements, A ∼ 130 and A ∼ 200, is concentrated on the polar angles 70° < θ < 120°, so roughly around the equator. For the first peak, the most abundant production of the elements with A ∼ 60 is found at trajectories below θ ∼ 70° and above θ ∼ 120°, so close to the polar axes. Because these trajectories give only a moderate contribution to the total mass outflow, the averaged abundance pattern is under-representing the first peak magnitude, relatively to the second and third one, as seen in Figure 7. The geometrical dependence of the abundance pattern could potentially reveal a higher magnitude of the first peak, in some specific cases of GRB-related kilonovae as observed closer to the polar axis.

As for the relative distributions of the elements with mass number A ∼ 130 (xenon group), their profiles resulting from our simulations match the solar data more closely for the short GRB models with higher-magnetized winds, and a less spinning black hole in their engine. This might tentatively suggest that such black holes should be more frequently produced via the binary neutron star mergers, which is in agreement with the moderate spins of black holes resulting from the simulation of the hypermassive neutron star that undergoes delayed collapse (Ruiz et al. 2016). We intend to investigate further the problem of the black hole spin influence on the kilonova signal in future work, implementing the three-dimensional scheme for the MHD turbulence to disentangle the effects of magnetic driving and BH rotation.