Black hole accretion in gamma ray bursts

TL;DR

This paper investigates the structure, evolution, and nucleosynthesis in hyperaccreting disks around black holes in gamma ray bursts, and explores their possible connection to gravitational wave events like GW150914.

Contribution

It introduces a detailed model of hyperaccreting disks in gamma ray bursts, including nucleosynthesis and potential observational signatures, linking them to gravitational wave detections.

Findings

Neutron-rich environments lead to heavier nuclei formation.

The model explains the timing and properties of the observed gamma ray burst related to GW150914.

Simulations reproduce low spin and kick velocities of merged black holes.

Abstract

We study the structure and evolution of the hyperaccreting disks and outflows in the gamma ray bursts central engines. The torus around a stellar mass black hole is composed of free nucleons, Helium, electron-positron pairs, and is cooled by neutrino emission. Accretion of matter powers the relativistic jets, responsible for the gamma ray prompt emission. The significant number density of neutrons in the disk and outflowing material will cause subsequent formation of heavier nuclei. We study the process of nucleosynthesis and its possible observational consequences. We also apply our scenario to the recent observation of the gravitational wave signal, detected on September 14th, 2015 by the two Advanced LIGO detectors, and related to an inspiral and merger of a binary black hole system. A gamma ray burst that could possibly be related with the GW150914 event was observed by the Fermi satellite. It had a duration of about 1 second and appeared about 0.4 seconds after the gravitational-wave signal. We propose that a collapsing massive star and a black hole in a close binary could lead to the event. The gamma ray burst was powered by a weak neutrino flux produced in the star remnant's matter. Low spin and kick velocity of the merged black hole are reproduced in our simulations. Coincident gravitational-wave emission originates from the merger of the collapsed core and the companion black hole.

Figures

Equations

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

$$P_{\mathrm{nucl}}=P_{\mathrm{e}-}+P_{\mathrm{e}+}+P_{\mathrm{n}}+P_{\mathrm{p}}$$

$$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+e^{-}\to n+{\nu }_{\mathrm{e}}$$

$$p+{\bar{\nu }}_{\mathrm{e}}\to n+e^{+}$$

$$p+e^{-}+{\bar{\nu }}_e\to n$$

$$n+e^{+}\to p+{\bar{\nu }}_{\mathrm{e}}$$

$$n\to p+e^{-}+{\bar{\nu }}_{\mathrm{e}}$$

$$n+{\nu }_{\mathrm{e}}\to p+e^{-},$$

$$e^{-}+e^{+}\to {\nu }_{\mathrm{i}}+{\bar{\nu }}_{\mathrm{i}}$$

$$n+n\to n+n+{\nu }_{\mathrm{i}}+{\bar{\nu }}_{\mathrm{i}}$$

$$\tilde{\gamma }\to {\nu }_{\mathrm{e}}+{\bar{\nu }}_{\mathrm{e}}.$$

$$Q_{\nu }^{-}=\frac{\frac{7}{8}\sigma T^4}{\frac{3}{4}}\sum _{i=e,\mu ,\tau }\frac{1}{\frac{{\tau }_{\mathrm{a},{\nu }_{\mathrm{i}}}+{\tau }_{\mathrm{s}}}{2}+\frac{1}{\sqrt{3}}+\frac{1}{3{\tau }_{\mathrm{a},{\nu }_{\mathrm{i}}}}}\times \frac{1}{H}[\mathrm{erg}{\mathrm{s}}^{-1}\mathrm{c}{\mathrm{m}}^{-3}].$$

$$(\rho u^{\mu }{)}_{;\mu }=0$$

$${T^{\mu }}_{\nu ;\mu }=0$$

$$P=K{\rho }^{\gamma }=(\gamma -1)u$$

$$T^{\mu \nu }=T_{gaz}^{\mu \nu }+T_{EM}^{\mu \nu }$$

$$T_{gaz}^{\mu \nu }=\rho h{u}^{\mu }{u}^{\nu }+p{g}^{\mu \nu }=(\rho +u+p)u^{\mu }{u}^{\nu }+p{g}^{\mu \nu }$$

$$T_{EM}^{\mu \nu }=b^2{u}^{\mu }{u}^{\nu }+\frac{1}{2}{b}^2{g}^{\mu \nu }-b^{\mu }{b}^{\nu };b^{\mu }=u_{\nu }{{}^{{}^{\ast }}F}^{\mu \nu }$$

$${\partial }_t\mathbf{U}(\mathbf{P})=-{\partial }_i{\mathbf{F}}^i(\mathbf{P})+\mathbf{S}(\mathbf{P}),$$

$${\dot{Y}}_{\mathrm{i}}=\sum _j{N}_j^i{\lambda }_j{Y}_j+\sum _{j,k}{N}_{j,k}^i\rho N_{\mathrm{A}}{Y}_j{Y}_k+\sum _{j,k,l}{N}_{j,k,l}^i{\rho }^2{N}_{\mathrm{A}}^2{Y}_j{Y}_k{Y}_l$$

$$T_{\mathrm{EM}}^{\mu \nu }=b^2{u}^{\mu }{u}^{\nu }+\frac{b^2}{2}{g}^{\mu \nu }-b^{\mu }{b}^{\nu },$$

$$L_{\mathrm{BZ}}\equiv \dot{E}=2\pi {\int }_0^{\pi }d\theta \sqrt{-g}{F}_E$$

Full text

\articlenumber

x \doinum10.3390/—— \pubvolumexx

\externaleditorAcademic Editor: name \historyReceived: date; Accepted: date; Published: date

\TitleBlack hole accretion in gamma ray bursts \AuthorAgnieszka Janiuk 1*, Michal Bejger 2, Petra Sukova 1, and Szymon Charzynski 1,3 \AuthorNamesAgnieszka Janiuk, Michal Bejger, Petra Sukova and Szymon Charzynski

\corresCorrespondence: agnieszka.janiuk@gmail.com

1 Introduction

Gamma-ray bursts (GRBs) are energetic, transient events observed in the sky at high energies. Their known cosmological origin requires a physical process that causes them to be a cosmic explosion of great power. Proposed mechanisms involve the creation of a black hole (BH) in a cataclysmic event. This may either result from the collapse of a massive rotating star, or via the merger of compact objects’ binary, e.g. neutron stars (NSs) or a BH and a neutron star. The so-called ‘central engine’ of this process is composed of a hot and dense accretion disk with a hyper-Eddington accretion rate (up to a few $M_{\odot}s^{-1}$) near a spinning BH.

In the hyperaccreting disks that are present around the BHs in the central engine, the densities and temperatures are so high, that the equation of state (EOS) can no longer be assumed to be of an ideal gas. The pressure and chemical balance is required by the nuclear reactions that take part between the free protons, neutrons, and electron-positron pairs abundant in the plasma. The charged particles must satisfy the total neutrality condition. In the $\beta$ reactions, neutrinos are produced in three flavors and are subject to absorption and scattering processes. The partially trapped neutrinos contribute to the pressure in the plasma, together with nucleons, electron-positron pairs, radiation, and synthesized heavier particles (Helium).

Numerical computations of the structure and evolution of the accretion flows in the gamma ray bursts engine were historically first carried out with the use of steady-state, and later with the time-dependent models. These models were axially and vertically averaged, and used the classical $\alpha$-parameter prescription for the viscosity (Popham et al. 1999; Di Matteo et al. 2002; Kohri et al. 2002, 2005; Chen & Beloborodov 2007; Reynoso et al. 2006; Janiuk et al. 2004; 2007; 2010). More recently, accretion flows are described by fully relativistic MHD computations (e.g., Barkov & Komissarov 2011; Janiuk et al. 2013). This method has also been applied recently to describe the central engine of a putative gamma ray burst which could be associated with the event GW150914 (Janiuk et al. 2017). Here we report on this work, and in addition, we also present new results on the nucleosynthesis of heavy elements in the accretion flow in this engine.

2 Accretion flow in the GRB engine

2.1 Equation of state

In the EOS, contribution to the pressure comes from free nuclei and $e^{+}-e^{-}$ pairs, helium, radiation and trapped neutrinos:

[TABLE]

where $P_{\rm nucl}$ includes free neutrons, protons, and the electron-positrons:

[TABLE]

with

[TABLE]

Here, $F_{\rm k}$ are the Fermi-Dirac integrals of the order $k$, and $\eta_{\rm e}$, $\eta_{\rm p}$ and $\eta_{\rm n}$ are the reduced chemical potentials, $\eta_{i}=\mu_{i}/kT$ is the degeneracy parameter, $\mu_{i}$ denoting the standard chemical potential. Reduced chemical potential of positrons is $\eta_{\rm e+}=-\eta_{\rm e}-2/\beta_{\rm e}$. Relativity parameters are defined as $\beta_{\rm i}=kT/m_{\rm i}c^{2}$. The term $P_{\rm rad}$, describing the pressure of radiation, scales with the temperature as $T^{4}$, and is added here together with the pressure of neutrinos, $P_{\nu}$, which is computed from the two-stream approximation (see details in Janiuk et al. 2007). The disk is fully opaque to photons, however the radiation pressure is still by several orders of magnitude smaller than the pressure of neutrinos. This EOS is computed numerically by solving the balance of nuclear reactions (Yuan 2005; Janiuk et al. 2007; Janiuk & Yuan 2010; see also Lattimer & Swesty 1991; Setiawan et al. 2004).

2.2 Neutrino cooling of the accretion flow in GRBs

In the hot and dense torus, with temperature of ${\sim}10^{11}$ K and density $>10^{10}$ g cm*-3*, neutrinos are efficiently produced. The main reactions that lead to their emission are the electron/positron capture on nucleons, as well as the neutron decay. Their nuclear equilibrium is described by the following reactions:

[TABLE]

and the forward and backward reaction rates are equal. The reaction rates are given by the appropriate integrals (Reddy, Prakash & Lattimer 1998; see also Appendix A in Janiuk et al. 2007).

Other neutrino emission processes are: electron-positron pair annihilation, bremsstrahlung and plasmon decay:

[TABLE]

[TABLE]

[TABLE]

We calculate their rates numerically, with proper integrals over the distribution function of relativistic, partially degenerate species.

The neutrino cooling rate is finally given by the two-stream approximation, and includes the scattering and absorptive optical depths for neutrinos of all three flavors (Di Matteo et al. 2002):

[TABLE]

The neutrino cooling is self-consistenlty computed from the balance of the nuclear reactions, listed in Section 2.2, whose rates govern the values of the absorptive optical depths for electron (all reactions) and muon or tau neutrinos (pair annihilation and bremsstrahlung reactions). Also, the neutrino scattering optical depth is computed, using the values of proton and neutron number densities, and Cabbibo angle of $\sin^{2}\theta_{\rm C}=0.23$ (see Yuan 2005, Janiuk et al. 2007). The emissivity is averaged over the disk height, instead of the integration over the neutrinosphere, whose shape is rather complex.

2.3 GR MHD scheme

Our simulations of the accretion flow dynamics in the GRB central engine were performed with a 2D code HARM (High Accuracy Relativistic Magnetohydrodynamics; Gammie et al. 2003). Our version of the code was extended to include the microphysics, described in the previous Section. The code provides a solver for the continuity and energy-momentum conservation equations:

[TABLE]

[TABLE]

[TABLE]

where:

[TABLE]

[TABLE]

[TABLE]

where $u^{\mu}$ is the four-velocity of gas, $u$ is the internal energy density, $b^{\mu}=\frac{1}{2}\varepsilon^{\mu\nu\rho\sigma}u_{\nu}F_{\rho\sigma}$ is the magnetic four-vector, and $F$ is the electromagnetic stress tensor. The metric tensor is denoted by $g^{\mu\nu}$. Within the force-free approximation, we have $E_{\nu}=u_{\mu}F^{\mu\nu}=0$.

HARM-2D uses a conservative numerical scheme to obtain solutions of equations of the following type:

[TABLE]

where ${\bf U}$ denotes a vector of “conserved” variables (momentum, energy density, number density, taken in the coordinate frame), ${\bf P}$ is a vector of ’primitive’ variables (rest mass density, internal energy), ${\bf F}^{i}$ are the fluxes, and $\mathbf{S}$ is a vector of source terms. In contrast to non-relativistic MHD, where ${\bf P}\rightarrow{\bf U}$ and ${\bf U}\rightarrow{\bf P}$ have a closed-form solution, in relativistic MHD, ${\bf U}({\bf P})$ is a complicated, nonlinear relation. Inversion ${\bf P}({\bf U})$ is calculated therefore numerically, at every time step. The inverse transformation requires a solution to 5 non-linear equations, done here by means of a multi-dimensional Newton-Raphson routine (see, e.g., Noble et al. 2006 for details).

The procedure is simple, if the pressure-density relation is adiabatic. However, for a general EOS, one must also compute $dp/dw$, $dp/dv^{2}$ and $p(W,v^{2},D)$, where $w=p+u+\rho$ denotes the enthalpy, $W=\gamma^{2}w$, $D=\sqrt{-g}\rho u^{t}$, and $v^{2}=g^{\mu\nu}u^{\mu}u^{\nu}$. In our scheme, we have tabulated $(p,u)(\rho,T)$ values and then interpolate over this table. The Jacobian $\partial{\bf U}/\partial{\bf P}$ is computed numerically. The conserved variables are then evolved in time.

HARM-2D has been MPI-parallelized for the hydro-evolution, and also implemented with the shared memory hyperthreading for the EOS-table interpolation. Our 2D simulations typically run up to t=2000-3000 M (within a couple weeks of computations on the local cluster, or a few days on Cray XC40 with 1k nodes). The computations are performed in the polar coordinate system, with resolution of 256 x 256 points in the $r$ and $\theta$ directions.

2.4 Example simulation results

Parameters used in our computations are: BH mass, with a fiducial value of $M=10M_{\odot}$, BH spin $a=0.6-0.98$, and accretion disk mass $M_{d}\approx 1M_{\odot}$, which gives the scaling for the density in physical units, needed for the EOS computations. To describe the rotationally-supported torus around a spinning BH we use the initial condition given by the equilibrium solution, first developed analytically by Fishbone & Moncrief (1976) and Abramowicz et al. (1978). We adopt the polar coordinate system, $r-\theta$, and use the Kerr-Schild coordinates to avoid the coordinate singularity on the BH horizon. Initially, a poloidal configuration of the field is assumed, with the vector potential $A_{\phi}=(\rho/\rho_{max})$, and initial normalization of $P_{gas}/P_{mag}=50$. Magnetic turbulences develop within the dense matter torus and the field is advected with the infalling gas under the BH horizon. For a rapidly-spinning BH, the open magnetic field lines form along the rotation axis, and the magnetically driven jets can be produced. Energy extracted from the rotating BH through the Blandford-Znajek process (Blandford & Znajek 1977) gives then a substantial contribution to the jet luminosity. In addition, annihilation of neutrino-antineutrino pairs produced in the torus is yet another source of power to the jets. In Figure 1, we show the result of an exemplary simulation, which shows a snapshot taken from evolved run, with the distribution of magnetic field lines, magnetic and gas pressure ratio, and neutrino emissivity, in the region close to a black hole.

3 Nucleosynthesis of heavy elements in the GRB engine

In the astrophysical plasma, thermonuclear fusion occurs due to the capture and release of particles ($n$, $p$, $\alpha$, $\gamma$). Reaction sequence produces further isotopes. The nuclear reactions may proceed with one (decays, electron-positron capture, photodissociacion), two (encounters), or three (triple alpha reactions) nuclei. Therefore, the change of abundance of the $i$-th isotope is in general given by:

[TABLE]

Abundances of the isotopes are calculated under the assumption of nucleon number and charge conservation for a given density, temperature and electron fraction ($T\leq 1MeV$). Integrated cross-sections depending on temperature $kT$ are determined with the Maxwell-Boltzmann or Planck statistics, and the background screening and degeneracy of nucleons must be taken into account. The set of resulting non-linear differential equations is solved using the Euler method (Wallerstein et al. 1997).

In our computations, we used the thermonuclear reaction network code, http://webnucleo.org, and we computed the nuclear statistical equilibria established for the fusion reactions, for which the data were taken from the JINA reaclib online database (Hix & Meyer 2006). This network is appropriate for temperatures below $1$ MeV, which is the case at the outer radii of accretion disks in GRB engines. The mass fraction of the isotopes is solved for converged profiles of density, temperature and electron fraction in the disk.

We find that the free nuclei are present mostly below 100-200 gravitational radii, and the plasma there is rather neutron rich. Then, most abundant heavy elements formed in the accretion tori in GRBs are Nickel, Iron and Cobalt, as well as Argonium, Titatnium, Cuprum, Zinc, Silicon, Sulphur, Clorium, Manganium, Titatnium, and Vanadium.

As a consequence of the heavy element formation, the enhanced emission in the lower energy bands, i.e., ultraviolet or optical, due to the decay of species, may accompany the GRB (the ’macronova’; e.g., Li & Paczynski 1998). Also, the radio flares, occurring months to years after the GRB can be observed (Piran et al. 2012). Furthermore, certain isotopes decay should be detectable via emission lines. For the NuSTAR satellite, sensitive in the $5-80$ keV range, it’s in principle possible to detect the photons from Ti, Co, Mn, Cu, Zn, Ga, Cr unstable isotopes decays. Additionally, the EPIC detector onboard the XMM-Newton, could be able to see lines below 15 keV, e.g., ${}^{45}Ti$, ${}^{57}Mn$, or ${}^{57}Co$ (Janiuk 2014). In Figure 3, we plot an example result of the nucleosynthesis computations, showing the volume integrated abundance distribution of elements, in the function of their mass number. The input was taken from the simulation of the GRB central engine, as designed to model the putative GRB that could possibly be associated with GW 150914 (see next Section).

4 Gravitational wave source GW150914

Data from Fermi Gamma-ray Burst Monitor (GBM) satellite suggested that the recently registered gravitational wave event GW150914 (Abbott et al. 2016), a coalescing binary BH system, is potentially related to a GRB. The electromagnetic radiation in high energy (above 50 keV) originated from a weak transient source and lasted for about 1 second (Connaughton et al. 2016). Its localization is broadly consistent with the sky location of GW150914.

We speculate here on the possible scenario for the formation of a GRB accompanied by the gravitational wave (GW) event. Even though the presence of the GRB was recently questioned by other instruments measurements (e.g., Greiner et al. 2016), we envisage the possibility of a GRB coincident with the GW event from the BH merger to be worth exploring, in the context of the future observations.

Our model invokes coalescence of a collapsing star’s nucleus that forms a BH, with its companion BH in a binary system (Janiuk et al. 2013, 2017). We find that the recoil velocity acquired by the final BH through the GW emission allows it to take only a small fraction of matter from the host star, provided specific configuration of the binary spin vectors with the orbital angular momentum. The GRB is produced on the cost of accretion of this remnant matter onto the final BH. The moderate spin of this BH accounts for the GRB jet to be powered by a weak neutrino emission rather than the Blandford-Znajek mechanism, and hence agrees with the low power observed in the Fermi GRB signal.

The semi-analytical treatment of the BH core collapse and the star’s envelope spin-up (see Janiuk et al. 2008, for details) is followed by the general-relativistic simulations of the binary BH merger. Then, the GRB central engine evolution is carried numerically, to find the energy output from the magnetized, neutrino-cooled torus that powers the electromagnetic burst.

4.1 Merger of binary BHs

The GW150914 event was interpreted to be a merger of two BHs of the masses of $36^{+5}_{-4}\ M_{\odot}$ and $29^{+4}_{-4}$ (Abbott et al. 2016). The final BH parameters were estimated to be of $62^{+4}_{-4}\ M_{\odot}$ and $0.67^{+0.05}_{-0.07}$ for its mass and spin. Probabilities that the angles between spins and the normal to the orbital plane are between $45^{\circ}$ and $135^{\circ}$ are about $0.8$ for each component BH. Spin magnitudes are constrained to be smaller than $0.7$ and $0.8$ at 90% probability. Assumption of a strict co-alignment of spins with the orbital angular momentum results in an upper limit of $0.2$ and $0.3$ for the spins. Distance to the source was that of $410^{+160}_{-180}$ Mpc, corresponding to a redshift of about $z=0.09$ (assuming standard cosmology).

We made several runs for BH mass and range of spins, constrained from the Advanced LIGO data. In Figure 3, we plot the example results from the run with BH mass ratio of $0.82$, and aligned spins of $0.2$ and $0.3$. The waveform is plotted using the Weyl scalar $\psi_{4}$ 111Curvature of spacetime is described by curvature tensor - the Riemman tensor (20 independent components), which can be decomposed into sum of Ricci tensor (10 independent components) and Weyl tensor (10 independent components). For vacuum solutions of Einstein equations Ricci tensor vanishes, so curvature is described by Weyl tensor only. In the Newman-Penrose formalism components of Weyl tensor are encoded as five complex Weyl scalars $\Psi_{0}$, … , $\Psi_{4}$. These are different components of curvature tensor, which have different physical interpretation. $\Psi_{4}$ describes the outgoing radiation for an asymptotically flat spacetime.. The final BH spin obtained in this simulation was $0.68$.

These runs were performed with the Einstein Toolkit computational package (Loeffler et al 2012). The numerical methods used here are based on finite difference computations on a gridded mesh and follow the inspiral, merger and ringdown phases of the binary black hole evolution, using the technique of the adaptive mesh refinement. The Cartesian grid covers the volume of $48\times 48\times 48M$. In AMR we use 7 levels of refinement, and the resolution ranges from $\Delta x=2.0M$ for the coarsest grid to $\Delta x=0.03125M$ for the finest grid. From the simulation runs, we constrained not only the spin of the merger black hole (it is then used to be a parameter in the GRB engine model), but also the velocity of the gravitational recoil. This turned out to be on the order of $\sim 700$ km/s. Following Kocsis et al. (2012), we argue that it is likely that accretion of the accumulated matter onto the final BH is triggered while it moves towards the circumbinary disk after the merger.

4.2 Weak GRB powered by neutrino emission

The Fermi GRB coincident with the GW150914 event had a duration of about 1 sec and appeared about 0.4 seconds after the GW signal. Within the limit of uncertainty of the Advanced LIGO and Fermi detectors’ capabilities it could also be associated spatially. The GRB fluence in the range 1 keV-10 MeV, is of $2.8\times 10^{-7}$ erg cm*-2*. The implied source luminosity in gamma rays equals to $1.8^{+1.5}_{-1.0}\times 10^{49}$ erg/s.

We modeled the GRB engine, constrained by the post-merger conditions. Using the code HARM-2D, with the numerical EOS and neutrino cooling implemented into MHD evolution, we estimated the neutrino and Blandford-Znajek luminosities available to power the GRB jet in this source. In particular, the mass of the accreting torus, created from the circumbinary matter, was tuned to produce adequately low neutrino luminosity for a weak GRB. Example parameters for the BH mass of $M=62M_{\odot}$, the BH spin of $a=0.7$, and disk mass $M_{d}=15M_{\odot}$ are consistent wit the weak GRB luminosity. Assuming low efficiencies of conversion between neutrino annihilation, and weak conversion between the jet kinetic and radiative power, this scenario is able to meet quantitatively the limits put by the Fermi data. The disk mass is here just an order of magnitude estimate. We first checked, that a much larger mass, on the order of the mass of the merged black hole, is too large to be consistent with the observed limitations for luminosity. On the other hand, the supernova explosion which might have left a 30 Solar mass black hole, should have originated from at least 80 Solar mass star on zero-age main sequence, in a low-metallicity environment (Abbott et al. 2016b, Spera et al. 2015). Even if a significant amount of the star’s mass was ejected during the evolution and explosion, some remnant of this order is plausible.

The luminosity $L_{\nu}$, emitted in neutrinos is simply equal to their volume integrated emissivity. To compute the energy flux through the horizon of the BH and the Blandford-Znajek luminosity, $L_{\rm BZ}$, we need the electromagnetic part of the stress tensor,

[TABLE]

where the four-vector $b^{\mu}$ with $b^{t}\equiv g_{i\mu}B^{i}u^{\mu}$ and $b^{i}\equiv(B^{i}+u^{i}b^{t})/u^{t}$. We then evaluate the radial energy flux as the power of the Blandford-Znajek process:

[TABLE]

where $F_{E}\equiv-{{T^{r}_{t}}}$. This can be subdivided into a matter ${F}^{(MA)}_{E}$ and electromagnetic ${F}^{(EM)}_{E}$ part, although in the force-free limit the matter part vanishes (McKinney & Gammie 2004).

In Figure 4, left panel, we plot the values of $L_{\rm BZ}$ as it varies with time for various simulations, depending on the black hole spin. In the right panel, we show the corresponding evolution of the neutrino luminosity. The maximum of the neutrino luminosity obtained in our simulations is reached about 0.4 seconds after an equilibrium torus, prescribed by our initial conditions, had formed. This may tentatively give the lower limit for the timescale when the jet appears after the BH merger. The Blandford-Znajek luminosity is on the order of $10^{50}-10^{51}$ ergs/s, and only occasionally non-zero for low spins (the jet is powered by magnetic field only episodically). We find therefore, that the GRB luminosity inferred from our simulations can be reconciled with the observational upper limits, for moderate spins of the final BH ($a=0.6-0.8$).

The connection of our computation to the particular event GW150914 and the corresponding Fermi GRB is arguable because of several facts. At first, the flare produced by our simulation is longer than the reported very short duration of the GRB (1s). However, the observed burst was close to the sensitivity limit of the instrument, so it is likely that it could observe just the very top of the flare and left the rest of the burst undetected. Moreover, if our line of sight is offset with respect to the jet axis and we can see only the edge of the jet, the detailed shape of the flare would be affected by the behaviour of the propagating jet. The time delay between the burst and gravitational signal can be significantly longer, if there is a massive star remnant, through which the jet has to propagate. The state of the remnant matter after the second BH collapse is very uncertain and needs to be studied in details in future works. However, we can speculate that the collapse and subsequent orbiting of the two BHs is violent enough to destroy the star, expel part of the matter away and form a dead circumbinary disc of the remnant matter, in which the accretion has stopped, around the BH binary. Then after the merger the new born black hole can receive a substantial kick towards the circumbinary disc (see our BH merger simulations in Section 4.1), can meet the matter in fraction of second and produce the burst. In that case, there is only little material left along the axis, through which the jet has to propagate, and so the time delay can be quite short, on the order of a second. However, this speculative scenario needs both more theoretical modeling and mainly more future observations of coincident gravitational and electromagnetic signals, which are well above the instrumental limits.

5 Conclusions

• •

We have carried out numerical GR MHD simulations of the magnetized torus accreting onto BH in the GRB central engine. The EOS is no longer assumed to be of an ideal gas, but accounts for the proper microphysics and neutrino cooling of the GRB engine.

• •

The neutrino emission and absorption processes account for an additional pressure component in the plasma, and the species are partially degenerate. The resulting density and temperature, as well as the electron fraction, are then taken into account when determining the nuclear equilibrium conditions.

• •

We determine the abundances of heavy elements (above Helium, up to the mass number of $\sim 100$), which are synthesized in the accretion disk in GRB engine, and in the winds that are magnetically launched from its surface. The observational signatures of radioactive decay of some of the unstable isotopes can be, for instance, the emission lines seen in the X-rays (in principle, in the NuSTAR range), and also in the faint emission in lower energy band continuum (i.e. ’kilonova’ or ’macronova’ scenario).

• •

We compute and compare the efficiencies of the GRB jet powered through the neutrino annihilation, and the Blandfrod-Znajek mechanism. We find that the constraints of a moderate BH spin and a small disk mass to the BH mass ratio, are in tentative agreement with the Advanced LIGO and Fermi data, if the latter was really coincident with the GW event.

• •

A definite answer and a test for our scenario will be possible if further searches for the GW events will provide more data, also for their electromagnetic counterparts, in both gamma rays and lower energies.

We claim, that even though the observations are still subject to debate, the possible GRB-BBH merger coincidence is worth to investigate. The timescales obtained in our scenario are possibly more plausible for a lonng gamma ray burst scenario than for a short one (see also Perna et al. 2016; Loeb 2016), unless the Fermi burst itself was much longer, but below the detction limit. On the other hand, the model proposed by Zhang (2016) can satisfy the duration and delay timsecales, but the intepretaion of a charged black hole and its observational signatures also needs further investigations.

Acknowledgements.

This work was supported in part by the grants no. DEC-2012/05/E/ST9/03914 and UMO-2014/14/M/ST9/00707 from the Polish National Science Center. We also acknowledge support from the Interdisciplinary Center for Mathematical Modeling of the Warsaw University, through the computational grant G53-5.