Black hole accretion in gamma ray bursts
Agnieszka Janiuk Michal Bejger Petra Sukova Szymon Charzynski








§ 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_⊙ 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 β 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 α-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.




§ ACCRETION FLOW IN THE GRB ENGINE








 §.§ Equation of state


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


    P = P_ nucl+P_ He+P_ rad+P_ν,


where P_ nucl includes free neutrons, protons,
and the electron-positrons:


    P_ nucl=P_ e-+P_ e++P_ n+P_ p


with


    P_ i = 2 √(2) 3π^2(m_ic^2)^4 (ħ c)^3β_i^5/2[F_3/2(η_ i,β_ i)+1 2β_ iF_5/2(η_ i,β_ i)].



Here, F_ k are the Fermi-Dirac integrals of the order k, and
η_ e, η_ p and η_ n are the reduced chemical
potentials, η_i = μ_i/kT
is the degeneracy
parameter, μ_i denoting the standard chemical
potential. Reduced chemical potential of positrons is
η_ e+=-η_ e-2/β_ e.
Relativity parameters are defined as β_ i=kT/m_ ic^2.
The term P_ rad, describing the pressure of radiation,
scales with the temperature as T^4, and is added here together with the
pressure of neutrinos, P_ν, 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).



 §.§ Neutrino cooling of the accretion flow in GRBs


In the hot and dense torus, with temperature of ∼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:


    p + e^-→ n + ν_ e
    
    p + ν̅_ e→ n + e^+
    
    p + e^- + ν̅_e→ n 
    
    n + e^+→ p + ν̅_ e
    
    n → p + e^- + ν̅_ e
    
    n + ν_ e→ p + e^-,

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:


    e^-+e^+→ν_ i+ν̅_ i


    n+n → n+n+ν_ i+ν̅_ i


    γ̃→ν_ e+ν̅_ e .

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):


    Q^-_ν = 7  8σ T^43  4∑_i=e,μ,τ 1 τ_ a, ν_i + τ_ s 2
    + 1 √(3) +
    1  3τ_ a, ν_i×1  H  [ erg  s^-1 cm^-3] .





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θ_ 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.



 §.§ 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:

    (ρ u^μ)_;μ = 0


    T^μ_ν;μ = 0


    P = Kρ^γ = (γ-1) u

where:

    T^μν = T^μν_gaz + T^μν_EM


    T^μν_gaz = ρ h u^μ u^ν + pg^μν =(ρ + u + p) u^μ u^ν + pg^μν


    T^μν_EM = b^2 u^μ u^ν + 1/2b^2 g^μν - b^μ b^ν ;    b^μ = u_ν^μν

where u^μ is the four-velocity of gas, u is the internal energy density,
b^μ = 1/2ϵ^μνρσu_ν F_ρσ is the magnetic four-vector,
and F is the electromagnetic stress tensor. The metric tensor is denoted by g^μν.
Within the force-free approximation, we have E_ν=u_μF^μν=0.

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


    _t () = -_i ^i() + 𝐒(),



where  denotes a vector of “conserved” variables (momentum, energy density,
number density, taken in the coordinate frame),  is a vector of
'primitive' variables (rest mass density, internal energy),
 ^i are the fluxes, and 𝐒 is a vector of source terms.


In contrast to non-relativistic MHD, where → and  → have a closed-form solution, in relativistic MHD,

() is a complicated, nonlinear relation. Inversion
() 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+ρ denotes the enthalpy, W=γ^2 w, D=√(-g)ρ u^t, and v^2=g^μνu^μu^ν.
In our scheme, we have tabulated (p,u)(ρ,T) values and then interpolate over this table.
The Jacobian ∂ /∂ 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 θ directions.



 §.§ Example simulation results


Parameters used in our computations are: BH mass, with a fiducial value of M=10 M_⊙, BH spin a=0.6-0.98, and accretion disk mass M_d≈ 1 M_⊙, 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-θ, 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_ϕ= (ρ/ρ_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.







§ 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, α, γ).
 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:

    Ẏ_ i = ∑_j N_j^iλ_jY_j +
    ∑_j,k N_j,k^iρ N_ A Y_jY_k +
    ∑_j,k,l N_j,k,l^iρ^2 N^2_ AY_jY_kY_l

Abundances of the isotopes are calculated under the assumption of nucleon number and charge conservation for a given density, temperature and electron fraction (T ≤ 1 MeV). 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., ^45Ti, ^57Mn,
or ^57Co (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).






§ 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.




 §.§ Merger of binary BHs
 

The GW150914 event was interpreted to be a merger of two BHs of the
masses of 36^+5_-4 M_⊙ and 29^+4_-4 (Abbott et al. 2016). The final BH parameters were estimated
to be of 62^+4_-4 M_⊙ 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^∘ and 135^∘ 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 ψ_4 . The final BH spin obtained in this simulation was 0.68.

Curvature 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 Ψ_0, ... , Ψ_4. These are different components of curvature tensor, which have different physical interpretation.
Ψ_4 describes the outgoing radiation for an asymptotically flat spacetime.

 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×48×48 M. In AMR we use 7 levels of refinement, and the
resolution ranges from Δ
x=2.0M for the coarsest grid to Δ 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 ∼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.



 §.§ 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 × 10^-7 erg cm^-2.
The implied source luminosity in gamma rays equals to 1.8^+1.5_-1.0× 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=62 M_⊙, the BH spin of a=0.7, and disk mass M_d=15 M_⊙ 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_ν, 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_ BZ, we need
the electromagnetic part of the stress tensor,

    T^μν_ EM =
    b^2 u^μ u^ν + b^22 g^μν - b^μ b^ν,

where the four-vector b^μ with b^t ≡
g_iμ B^i u^μ and b^i ≡ (B^i + u^i b^t)/u^t.
We then evaluate the radial energy flux as the power of the Blandford-Znajek process:

    L_ BZ≡Ė = 2π∫_0^π  dθ √(-g)F_E

where F_E ≡ -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_ 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 <ref>), 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.






§ 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 ∼ 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.








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.





mdpi






999

ref-journal Abbott, B.P., et al. Observation of Gravitational Waves from a Binary Black Hole Merger. Phys. Rev. Lett. 2016a, 116, 061102

ref-journal Abbott, B.P., et al. Astrophysical Implications of the Binary Black-hole Merger GW150914. ApJL 2016b, 818, L22 

ref-journal Abramowicz, M.; Jaroszynski, M.; Sikora, M. Relativistic, accreting disks. A&A 1978, 63, 221

ref-journal Barkov, M.V.; Baushev, A.N. Accretion of a massive magnetized torus on a rotating black hole. New Astr. 2011, 16, 46

ref-journal Barkov, M.V., Komissarov, S.S. Close binary progenitors of gamma-ray bursts. MNRAS 2010, 401, 1644

ref-journal Blandford, R.D., Znajek, R.L. Electromagnetic extraction of energy from Kerr black holes. MNRAS 1977 179, 433

ref-journal Chen, W.X.; Beloborodov, A.M. Neutrino-cooled Accretion Disks around Spinning Black Holes. ApJ 2007, 657, 383

ref-journal Connaughton, V., et al. Fermi GBM Observations of LIGO Gravitational-wave Event GW150914. ApJ Lett. 2016, 826, 6

ref-journal Di Matteo, T.; Perna, R.; Narayan, R.. Neutrino Trapping and Accretion Models for Gamma-Ray Bursts. ApJ 2002, 579, 706
	
ref-journal Fishbone, L. G.; Moncrief, V. Relativistic fluid disks in orbit around Kerr black holes. ApJ 1976, 207, 962
	
ref-journal Gammie, C. F.; McKinney, J.C.; Tóth, G. HARM: A Numerical Scheme for General Relativistic Magnetohydrodynamics. ApJ 2003, 589, 444

ref-journal Greiner, J; Burgess, J. M.; Savchenko, V.; Yu, H.-F. On the Fermi-GBM Event 0.4 s after GW150914. ApJ Lett. 2016, 827, 38

ref-journal Hix, W.R., Meyer, B.S., Thermonuclear kinetics in astrophysics. Nucl Phys., 2006 777, 188

ref-journal Janiuk, A.; Bejger, M.; Charzynski, S.; Sukova, P. On the possible gamma-ray burst-gravitational wave association in GW150914. New Astr. 2017, 51, 7

ref-journal Janiuk, A. Nucleosynthesis of heavy elements in gamma ray burst engines. A&A 2014, 568, 105



ref-journal Janiuk, A.; Mioduszewski, P.; Moscibrodzka, M. Accretion and Outflow from a Magnetized, Neutrino Cooled Torus around the Gamma-Ray Burst Central Engine. ApJ 2013, 776, 105

ref-journal Janiuk, A.; Charzynski, S.; Bejger, M. Long gamma ray bursts from binary black holes. A&A 2013, 560, 25

ref-journal Janiuk, A.; Proga, D.; Moderski, R. On the Duration of Long GRBs: Effects of Black Hole Spin.  ApJ 2008, 687, 433

ref-journal Janiuk, A.; Yuan, Y.F. The role of black hole spin and magnetic field threading the unstable neutrino disk in gamma ray bursts. A&A, 2010, 509, 55

ref-journal Janiuk, A.; Yuan, Y.F., Perna, R.; Di Matteo, T. Instabilities in the Time-Dependent Neutrino Disk in Gamma-Ray Bursts. ApJ 2007, 664, 1011

ref-journal Janiuk, A.; Perna, R., Di Matteo, T.; Czerny, B. Evolution of a neutrino-cooled disc in gamma-ray bursts. MNRAS, 2004, 355, 950

ref-journal Kocsis, B.; Haiman, Z.; Loeb, A. Gas pile-up, gap overflow and Type 1.5 migration in circumbinary discs: general theory. MNRAS 2012, 427, 2660

ref-journal Kohri, K.; Mineshige, S. Can Neutrino-cooled Accretion Disks Be an Origin of Gamma-Ray Bursts? ApJ 2002, 577, 311


ref-journal Kohri, K.; Narayan, R.; Piran, T. Neutrino-dominated Accretion and Supernovae. ApJ 2005, 629, 341

ref-journal Löffler, F.; et al. The Einstein Toolkit: a community computational infrastructure for relativistic astrophysics. Class. Quantum Grav. 2012, 29, 115001

ref-journal Li, L.X., Paczynski, B.P. Transient Events from Neutron Star Mergers, ApJ 1998, 507, 59

ref-journal Loeb, A. Electromagnetic Counterparts to Black Hole Mergers Detected by LIGO. ApJL 2016, 819, 21

ref-journal Mc Kinney, J.C.; Gammie, C.F. A Measurement of the Electromagnetic Luminosity of a Kerr Black Hole. ApJ 2004, 611, 977

ref-journal Noble, S.C.; Gammie, C.F.; McKinney, J.C.; Del Zanna, L.
Primitive Variable Solvers for Conservative General Relativistic Magnetohydrodynamics. ApJ 2006, 641, 626

ref-journal Perna, R.; Lazzati, D.; Giacomazzo, B.; Short Gamma-Ray Bursts from the Merger of Two Black Holes. ApJL 2016, 821, 18

ref-journal Popham, R.; Woosley S.E.; Fryer, C. Hyperaccreting Black Holes and Gamma-Ray Bursts. ApJ 1999, 518, 356

ref-journal Reynoso, M. M.; Romero, G. E.; Sampayo, O. A. Precession of neutrino-cooled accretion disks in gamma-ray burst engines. A& A 2006, 454, 11

ref-journal Spera, M.; Mapelli, M.; & Bressan, A. The mass spectrum of compact remnants from the PARSEC stellar evolution tracks. MNRAS 2015, 451, 4086

ref-journal Yuan, Y.-F. Electron-positron capture rates and a steady state equilibrium condition for an electron-positron plasma with nucleons. Phys. Rev. D 2005, 72, 013007

ref-journal Wallerstein, G. et al. Synthesis of the elements in stars: forty years of progress.  Rev. Mod. Phys. 1997, 69, 995

ref-journal Zhang; B. Mergers of Charged Black Holes: Gravitational-wave Events, Short Gamma-Ray Bursts, and Fast Radio Bursts. ApJL 2016, 827, 31


















