TY - JOUR
AU - Murase,, Kohta
AB - Abstract The first direct detections of gravitational waves (GWs) from black hole (BH) mergers, GW 150914, GW 151226 and LVT 151012, give a robust lower limit |${\sim } 70\,000^{+170\,000}_{-61\,000}$| on the number of merged, highly spinning BHs in our Galaxy. The total spin energy is comparable to all the kinetic energy of supernovae that ever happened in our Galaxy. The BHs release the spin energy to relativistic jets by accreting matter and magnetic fields from the interstellar medium (ISM). By considering the distributions of the ISM density, BH mass and velocity, we calculate the luminosity function of the BH jets, and find that they can potentially accelerate TeV–PeV cosmic ray particles in our Galaxy with total power ∼1037 ± 3 erg s−1 as PeVatrons, positron factories and/or unidentified TeV gamma-ray sources. Additional ∼300 BH jet nebulae could be detectable by Cherenkov Telescope Array. We also argue that the accretion from the ISM can evaporate and blow away cold material around the BH, which has profound implications for some scenarios to predict electromagnetic counterparts to BH mergers. acceleration of particles, gravitational waves, stars: black holes, cosmic rays, ISM: jets and outflows, gamma-rays: general 1 INTRODUCTION A century after Einstein predicted the existence of gravitational waves (GWs), the Laser Interferometer Gravitational-Wave Observatory (LIGO) observed the first direct GW signal GW 150914 from a merger of two black holes (BHs) with masses of |$36_{-4}^{+5}$| and |$29_{-4}^{+4} \,\mathrm{M}_{{\odot }}$| and radiated energy |$3_{-0.4}^{+0.5} \,\mathrm{M}_{{\odot }}\,c^2$| (Abbott et al. 2016a). This is also the first discovery of a binary BH. During Advanced LIGO’s first observing period (O1), 2015 September 12 to 2016 January 19 (The LIGO Scientific Collaboration et al. 2016),1 the second event GW 151226 with masses |$14.2_{-3.7}^{+8.3}$| and |$7.5_{-2.3}^{+2.3} \,\mathrm{M}_{{\odot }}$| and radiated energy |$1.0_{-0.2}^{+0.1} \,\mathrm{M}_{{\odot }}\,c^2$| (Abbott et al. 2016b) and a candidate LVT 151012 with |$23_{-6}^{+18}$|⁠, |$13_{-5}^{+4}$| and |$1.5_{-0.4}^{+0.3} \,\mathrm{M}_{{\odot }}\,c^2$| have been also detected, and the existence of a population of merging BHs has been established. These ∼2.5 events give a relatively certain estimate on the merger rate in the range |$\mathscr {R}_{\rm GW} \sim 9$|–240 Gpc−3 yr−1 (The LIGO Scientific Collaboration et al. 2016). A new era of GW astrophysics has been finally opened and will be driven by a network of LIGO, Virgo, KAGRA and Indian Initiative in Gravitational-wave Observations (IndiGO), and by Laser Interferometer Space Antenna (LISA) and Deci-Hertz Interferometer Gravitational wave Observatory (DECIGO) satellites in the future (Kyutoku & Seto 2016; Nakamura et al. 2016b; Sesana 2016). The binary BH mergers are the most luminous events in the Universe, even brighter than gamma-ray bursts (GRBs). The peak luminosities are |$3.6_{-0.4}^{+0.5} \times 10^{56}$|⁠, |$3.3_{-1.6}^{+0.8} \times 10^{56}$| and |$3.1_{-1.8}^{+0.8} \times 10^{56}$| erg s−1 for GW 150914, GW 151226 and LVT 151012, respectively (The LIGO Scientific Collaboration et al. 2016), which reach ∼0.1 per cent of the Planck luminosity c5/G = 3.6 × 1059 erg s−1 = 2.0 × 105 M⊙ c2 s−1. Merged BHs also retain huge energy in the spin. The spin energy is about \begin{eqnarray} E_{\rm spin} &=& \left(1-\sqrt{\frac{1+\sqrt{1-a_*^2}}{2}}\right)Mc^2\nonumber\\ & \sim& 1 \times 10^{54}\, {\rm erg} \left(\frac{M}{10 \,\mathrm{M}_{\odot }}\right), \end{eqnarray} (1) where the spin parameter is typically a* = a/M ∼ 0.7 after a merger (e.g. Zlochower & Lousto 2015). On the contrary, the spin parameter is usually unknown for the field BHs before a merger and could be very small a* ≪ 0.7 as a result of the mass loss from the progenitor stars. Post-merger spinning BHs should also exist in our Galaxy, having a lot of energy in the spin. The number of such BHs is estimated as \begin{eqnarray} N_{\rm BH} \sim \mathscr {R}_{\rm GW} n_{\rm gal}^{-1} H_0^{-1} \sim 7 \times 10^{4}\, {\rm galaxy}^{-1} \left(\frac{\mathscr {R}_{\rm GW}}{70\,{\rm Gpc}^{-3}\,{\rm yr}^{-1}}\right), \end{eqnarray} (2) where we use |$\mathscr {R}_{\rm GW} \sim 70$| Gpc−3 yr−1 (The LIGO Scientific Collaboration et al. 2016), ngal ∼ 0.01 Mpc−3 is the number density of galaxies and |$H_0^{-1} \sim 10^{10}$| yr is the Hubble time. This estimate is applicable unless the merger rate changes very rapidly in a time much shorter than the Hubble time. Note that, although the large mass in GW 150914 suggests a low-metallicity environment with Z ≲ Z⊙/2 (Abbott et al. 2016c), our Galaxy had a low-metallicity environment in the past, and also incorporated low-metallicity galaxies in the hierarchical structure formation. The total spin energy stored in the merged BHs in our Galaxy is \begin{eqnarray} E_{\rm tot}=N_{\rm BH} E_{\rm spin} \sim 9 \times 10^{58}\, {\rm erg} \sim 9 \times 10^{7} E_{\rm SN}, \end{eqnarray} (3) where ESN ∼ 1051 erg is the kinetic energy of a supernova (SN). This is comparable to the total energy of SNe that ever happened in our Galaxy, i.e. ∼108 SNe exploded during the Hubble time!2 This is a robust lower limit on the total spin energy, obtained by the GW observations for the first time, because the spin energy of the field BHs could be very small. A natural question arises: How much spin energy is extracted from the merged BHs in our Galaxy? The spin energy of a BH can be extracted by a large-scale poloidal magnetic field threading the BH, i.e. through Blandford–Znajek (BZ) effect (Blandford & Znajek 1977), which is thought to produce a BH jet. We show that a sufficient magnetic field is advected by the Bondi–Hoyle accretion from the interstellar medium (ISM) and the jet power becomes comparable to the accretion rate, which is larger than the radiative power of the accretion disc. By taking into account the distributions of the ISM density and the BH mass and velocity, we estimate the luminosity function and the total power of the BH jets. Similar kinds of jets would have been already observed in X-ray binary systems, particularly associated with the low-hard state (e.g. Gallo, Fender & Pooley 2003; Fender, Belloni & Gallo 2004; Begelman & Armitage 2014). Based on the estimate of the luminosities and the acceleration energy, we suggest that the BH jets are potentially the origin of high-energy particles in our Galaxy. There are enigmatic high-energy sources in our Galaxy, such as still unknown PeVatrons accelerating cosmic rays (CRs) up to the knee energy εknee ∼ 3 × 1015 eV and beyond, sources of TeV CR positrons and unidentified TeV sources (TeV unIDs) that are dominant in the very high energy gamma-ray sky. These sources require only a small fraction of the spin energy Etot and could be powered by the BH jets. Our examination of the BH accretion and jet also suggests that it is very difficult to detect an electromagnetic counterpart to a BH merger after a GW event. In particular, the report of a GRB around the time of GW 150914 by the Fermi Gamma-ray Burst Monitor (GBM) (Connaughton et al. 2016) is most likely irrelevant to the GW event. This is consistent with a large number of follow-up searches after GW 150914 (Abbott et al. 2016d; Abe et al. 2016; Ackermann et al. 2016; Adriani et al. 2016; Adrián-Martínez et al. 2016; Evans et al. 2016a,b; KamLAND Collaboration et al. 2016; Kasliwal et al. 2016; Morokuma et al. 2016; Palliyaguru et al. 2016; Tavani et al. 2016; The Pierre Auger Collaboration et al. 2016; Troja et al. 2016). The organization of this paper is as follows. In Section 2, we discuss the physical mechanism of energy extraction from a spinning BH. We find that the accretion disc typically results in the so-called magnetically arrested disc (MAD) state and the magnetic field extracts the spin energy with the maximum efficiency for producing a jet. In Section 3, we calculate the luminosity function of the BH jets by considering the distributions of the BH mass, the peculiar velocity, the GW recoil velocity and the ISM density. The luminosity function also gives the total power of the BH jets. In Section 4, we discuss the connections of the BH jets with high-energy sources in our Galaxy, such as PeVatrons, CR positron sources and TeV unIDs. In Section 5, we encompass the uncertainties of our estimate on the total power within a factor of 10±3 by taking into account various effects such as the initial spin, the BH formation scenario and the wind feedbacks. This is much better than before; the factor was almost 10±∞ before the GW detections. In Section 6, we show that BHs are difficult to keep accretion discs until the merger that are massive enough for making a detectable electromagnetic counterpart for GW 150914. Section 7 is devoted to the summary and discussions. In Appendix A, we clarify novel points of our work compared with previous studies. 2 EXTRACTING SPIN ENERGY OF GW 150914-LIKE GALACTIC BHs The spin energy of a BH can be extracted by a large-scale magnetic field threading the BH ergosphere. The BH spin twists the magnetic field and the twisted magnetic field carries energy outward as a Poynting jet. This is the so-called BZ effect (Blandford & Znajek 1977; Koide et al. 2002). Although the BH itself cannot keep the magnetic field because of the no-hair theorem, accretion on to the BH can maintain the magnetic field on the BH. In this section we consider a BH in the ISM and estimate the luminosity of a BZ jet powered by the BH spin. For typical parameters, we find that the luminosity of a BH jet is comparable to the accretion rate |$L_{{\rm j}} \approx {\dot{M}} c^2$|⁠, with the accretion disc in the state of the so-called MAD. 2.1 Bondi accretion from the ISM The accretion rate on to a BH from the ISM is given by the Bondi–Hoyle rate (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; Bondi 1952), \begin{eqnarray} \dot{M} &=& 4 \pi r_{\rm B}^2 V \rho =\frac{4\pi G^2 M^2 n \mu m_{\rm u}}{V^{3}} \nonumber \\ &\sim& 5 \times 10^{35} \, {\rm erg}\, {\rm s}^{-1} \frac{1}{c^2} \left(\frac{n}{10\,{\rm cm}^{-3}}\right) \left(\frac{M}{10 \,\mathrm{M}_{\odot }}\right)^2 \left(\frac{V}{10\,{\rm km}\,{\rm s}^{-1}}\right)^{-3}\nonumber\\ &\sim& 4\times 10^{-4} \frac{{L}_{\rm Edd}}{c^2}\left(\frac{n}{10\,{\rm cm}^{-3}}\right)\!\left(\frac{M}{10 \,\mathrm{M}_{\odot }}\right) \!\left(\frac{V}{10\,{\rm km}\,{\rm s}^{-1}}\right)^{-3}\!, \end{eqnarray} (4) where n is the number density of the ISM, mu is the unified atomic mass unit, the mean molecular weight is μ = 1.41 for the Milky Way abundance and μ = 2.82 for molecular clouds (e.g. Kauffmann et al. 2008), M is the mass of the merged BH, LEdd is the Eddington luminosity, \begin{eqnarray} r_{\rm B} = \frac{GM}{V^2} \sim 1 \times 10^{15}\, {\rm cm} \left(\frac{M}{10 \,\mathrm{M}_{\odot }}\right) \left(\frac{V}{10\,{\rm km}\,{\rm s}^{-1}}\right)^{-2} \end{eqnarray} (5) is the Bondi radius and \begin{eqnarray} V=\sqrt{c_{\rm s}^2+v^2+v_{\rm GW}^2} \end{eqnarray} (6) includes the (effective) sound speed cs of the ISM, the centre-of-mass velocity v of the BH before the merger in the local ISM and the recoil velocity vGW due to the GW emission at the merger. The accretion rate |${\dot{M}}$| is proportional to M2v−3n. The discovery of a massive BH with mass M ∼ 60 M⊙ in GW 150914 significantly increases the estimate of |${\dot{M}}$|⁠, while the GW recoil tends to reduce it. The ISM density spans many decades. Thus we have to consider the distributions of mass, velocity and density to estimate the total power in Section 3. 2.2 Formation of an accretion disc and ADAF The accreted matter forms an accretion disc for typical parameters (Fujita et al. 1998; Agol & Kamionkowski 2002). The ISM density has a turbulent fluctuation with a Kolmogorov spectrum δρ/ρ ∼ [L/(6 × 1018 cm)]1/3 down to ∼108 cm (Armstrong, Rickett & Spangler 1995; Draine 2011). As a BH travels in the ISM, the accreting matter acquires a net specific angular momentum, \begin{eqnarray} \ell \sim \frac{1}{4}\frac{\Delta \rho }{\rho } V r_{\rm B}, \end{eqnarray} (7) where |$\Delta \rho /\rho =\delta \rho /\rho |_{L=2r_{\rm B}}$| is the density difference across the accretion cylinder.3 By equating this with the Keplerian angular momentum |$\ell _{\rm K}=\sqrt{GM r_{\rm disc}}$|⁠, we obtain the radius of the resulting accretion disc, \begin{eqnarray} \frac{r_{\rm disc}}{r_{\rm S}} &\sim & \frac{1}{16}\left(\frac{2 GM/V^2}{6\times 10^{18}\,{\rm cm}}\right)^{2/3} \frac{c^2}{2V^2} \nonumber \\ &\sim & 2\times 10^{5} \left(\frac{M}{10 \,\mathrm{M}_{{\odot }}}\right)^{2/3} \left(\frac{V}{10\,{\rm km}\,{\rm s}^{-1}}\right)^{-10/3}, \end{eqnarray} (8) where rS = 2GM/c2 is the Schwarzschild radius. The disc radius could be decreased if the magnetic breaking is effective. The accretion disc most likely forms hot, geometrically thick accretion flow such as advection-dominated accretion flow (ADAF; Fujita et al. 1998). The accreted matter is heated and eventually ionized because the collisional ionization rate is larger than the accretion rate and the recombination rate for typical parameters (see also Section 6). The accretion rate is much lower than the Eddington rate as in equation (4) and hence the low density makes the cooling ineffective (Ichimaru 1977; Narayan & Yi 1994, 1995; Kato, Fukue & Mineshige 2008; Yuan & Narayan 2014). The radiated energy from ADAF is much less than the total generated energy and almost all energy is advected into the BH (see also Section 5.3). For example, the luminosity of bremsstrahlung emission from electrons is only |$L_{\rm brem}\sim (\alpha _{\rm QED}/\alpha ^2)(m_{\rm e}/m_{\rm u}) ({\dot{M}}c^2/L_{\rm Edd}) {\dot{M}} c^2 \ll {\dot{M}} c^2$|⁠, where αQED is the fine-structure constant, α is the viscous parameter and me is the electron mass. As shown below, this is much smaller than the jet luminosity. Thus we concentrate on the jet in this paper and consider the disc emission in the future papers (Matsumoto, Teraki & Ioka 2017; Matsumoto & Ioka, in preparation). A transition to a cold standard disc outside the hot disc is not expected for typical parameters, although this is common in BH X-ray binaries (e.g. Esin, McClintock & Narayan 1997; Kato et al. 2008). The reason is that at the initial radius in equation (8), the disc is already hot (ionized) and the maximum accretion rate of the ADAF solution (Abramowicz et al. 1995) is larger than the accretion rate in equation (4), i.e. cooling is ineffective. Then we do not also expect soft X-ray transients (or X-ray novae) caused by the accumulation of the accreted matter at some radius because the thermal instability due to recombination is absent for the ionized flow (e.g. Kato et al. 2008). 2.3 Blandford–Znajek jet from a MAD state The accretion of the ISM also drags magnetic fields into the BH (see Fig. 1). The magnetic fields are well frozen in the accreting fluid because the loss time of the magnetic flux in the ISM is much longer than the accretion time (Nakano, Nishi & Umebayashi 2002). The formed disc is also thick, being able to advect the magnetic flux inward (Beckwith, Hawley & Krolik 2009; Cao 2011), although a thin disc cannot conserve a magnetic flux (Lubow, Papaloizou & Pringle 1994). Such a picture of flux advection is also consistent with the observations of active galactic nuclei (AGNs; Sikora & Begelman 2013) and BH binaries (Begelman & Armitage 2014). The coherent length of the magnetic field in the ISM is much larger than the Bondi radius, approximately about the scale of energy injection by SNe and stellar winds ∼1–10 pc (Han, Ferriere & Manchester 2004). Then the magnetic flux conservation implies the magnetic field strength on the horizon, \begin{eqnarray} B_{\rm H} \sim \left(\frac{r_{\rm B}}{r_{\rm H}}\right)^2 B_{\rm ISM}, \end{eqnarray} (9) where BISM is the magnetic field strength in the ISM, and |$r_{\rm H}=\frac{1}{2}\left(1+\sqrt{1-a_*^2}\right) r_{\rm S}$| is the radius of the BH horizon. On the other hand, for a given accretion rate, there is a maximum strength of the magnetic field on the horizon, \begin{eqnarray} B_{\rm H} &\sim & \left.\sqrt{\frac{4GM{\dot{M}}}{r^3 v_{\rm r}}}\right|_{r=r_{\rm H}} \nonumber \\ &\sim & 4 \times 10^{7}\, {\rm G} \left(\frac{n}{10\,{\rm cm}^{-3}}\right)^{1/2} \left(\frac{V}{10\,{\rm km}\,{\rm s}^{-1}}\right)^{-3/2}, \end{eqnarray} (10) because the pressure of the magnetic field, \begin{eqnarray} p_B=\frac{B^2}{8\pi }, \end{eqnarray} (11) cannot exceed the ram pressure of the accreting matter, \begin{eqnarray} p_{{\rm a}} = \frac{GM\Sigma }{r^2} \sim \frac{GM {\dot{M}}}{2\pi r^3 v_{\rm r}}, \end{eqnarray} (12) where |$\Sigma ={\dot{M}}/2\pi r v_{\rm r}$| is the surface density, vr ≡ εvff is the radial velocity, |$v_{\rm ff}=\sqrt{3GM/4\pi r}$| is the free-fall time and ε ∼ 0.05 is suggested by the numerical simulations and observations (Tchekhovskoy, Narayan & McKinney 2011; Zamaninasab et al. 2014). Although accumulation of the magnetic flux with the same polarity makes a magnetic barrier (Bisnovatyi-Kogan & Ruzmaikin 1976), the accretion continues through a magnetic flux via magnetic interchange instability (e.g. Arons & Lea 1976; McKinney, Tchekhovskoy & Blandford 2012). Such a magnetically dominated state is the so-called MAD (Bisnovatyi-Kogan & Ruzmaikin 1976; Narayan, Igumenshchev & Abramowicz 2003; Tchekhovskoy et al. 2011). The MAD state is realized if BH in equation (9) is larger than that in equation (10), i.e. \begin{eqnarray} B_{\rm ISM} &>& B_{\rm crit} \equiv \left(\frac{r_{\rm H}}{r_{\rm B}}\right)^2 \left.\sqrt{\frac{4GM{\dot{M}}}{r^3 v_{\rm r}}}\right|_{r=r_{\rm H}} \nonumber \\ &\sim & 1 \times 10^{-10}\,{\rm G} \left(\frac{V}{10\,{\rm km}\,{\rm s}^{-1}}\right)^{5/2} \left(\frac{n}{10\,{\rm cm}^{-3}}\right)^{1/2}. \end{eqnarray} (13) This is usually satisfied for typical BISM ∼ 3 μG since the velocity after the merger is roughly v ∼ 100 km s−1 (see Section 3.5). Thus the formed disc is likely MAD. Figure 1. Open in new tabDownload slide Schematic picture of a Blandford–Znajek jet from a spinning BH that is accreting from the ISM. Figure 1. Open in new tabDownload slide Schematic picture of a Blandford–Znajek jet from a spinning BH that is accreting from the ISM. A spinning BH immersed in large-scale poloidal magnetic fields releases energy through the BZ effect with a Poynting luminosity: \begin{eqnarray} L_{{\rm j}} \approx \frac{\kappa }{4\pi c} \Omega _{\rm H}^2 \Psi _{\rm BH}^2, \end{eqnarray} (14) where κ ≈ 0.05 (Tchekhovskoy et al. 2011), ΩH = a*c/2rH is the angular frequency of the BH and |$\Psi _{\rm BH} \sim \pi r_{\rm H}^2 B_{\rm H}$| is a magnetic flux on the BH. For a MAD state pB ∼ pa, the BZ luminosity is calculated as \begin{eqnarray} L_{{\rm j}} \sim \left(\frac{\kappa }{\epsilon } \sqrt{\frac{\pi ^3}{12}\frac{r_{\rm S}}{2r_{\rm H}}}\right) a_*^2 {\dot{M}} c^2 \approx {\dot{M}} c^2 \end{eqnarray} (15) for a typical spin parameter after the merger a* ≈ 0.7. Note that the nearly 100 per cent efficiency is correct within a factor for the spin parameter a* ≈ 0.7 (Tchekhovskoy et al. 2011). The jet component is not a radiative but a magnetic (and kinetic) component, so that the sub-Eddington accretion system is still radiatively inefficient until the jet is dissipated. In the following, we will use |$L_{{\rm j}} \approx {\dot{M}} c^2$| to estimate the jet luminosity of the merged BHs.4 Note that the net angular momentum direction of the accretion flow changes on a time-scale of crossing the density fluctuation ta ∼ rB/v ∼ 40 yr(M/10 M⊙)(v/10 km s−1)−3. However, the angular momentum vector near the BH is forced to align with the BH spin by the Bardeen–Petterson effect (Bardeen & Petterson 1975). In addition, although the direction of the poloidal magnetic fields is generally different from the BH spin direction, the magneto-spin alignment is also realized by the frame-dragging effect (McKinney, Tchekhovskoy & Blandford 2013). Therefore, we can consider that the direction of the jet is the same as that of the BH spin. 3 LUMINOSITY FUNCTION OF GW 150914-LIKE GALACTIC BH JETS Since the accretion rate depends on nM2v−3 that spans many decades, we calculate the luminosity function of jets from GW 150914-like merged BHs in our Galaxy as \begin{eqnarray} &&\!\!\!\frac{{\rm d}N}{{\rm d}{\dot{M}}} = N_{\rm BH} \int {\rm d}m_1\, \frac{{\rm d}p(m_1)}{{\rm d}m_1} \int {\rm d}m_2\, \frac{{\rm d}p(m_2|m_1)}{{\rm d}m_2} \int {\rm d}v\, \frac{{\rm d}f(v)}{{\rm d}v} \nonumber \\ &&\ \times \int {\rm d}n\, \frac{{\rm d}\xi (n)}{{\rm d}n} h(m_1,m_2,v) \,\delta \left[\dot{M}(n, m_1, m_2, v)-\dot{M}\right]\!, \end{eqnarray} (16) where dp(m1)/dm1 and dp(m2|m1)/dm2 are the distributions of BH masses (Section 3.1), df(v)/dv is the distribution of the pre-merger velocity (Section 3.2), dξ(n)/dn is the distribution of the ISM density (Section 3.3) and h(m1, m2, v) is the correction factor due to the scale heights of the ISM phases and BH distributions (Section 3.4). First, the δ function can be integrated over v analytically as \begin{eqnarray} \frac{{\rm d}N}{{\rm d}{\dot{M}}} &=& N_{\rm BH} \int {\rm d}m_1\, \frac{{\rm d}p(m_1)}{{\rm d}m_1} \int {\rm d}m_2\, \frac{{\rm d}p(m_2|m_1)}{{\rm d}m_2} \int {\rm d}n\, \frac{{\rm d}\xi (n)}{{\rm d}n} \nonumber \\ &&\times \,\, h(m_1,m_2,v_0) \frac{{\rm d}f(v_0)}{{\rm d}v} \frac{V_{v=v_0}^2}{3 v_0 {\dot{M}}}, \end{eqnarray} (17) where |$v_0^2\equiv (4\pi G^2 M^2 n \mu m_{\rm u}/{\dot{M}})^{2/3}-c_{\rm s}^2-v_{\rm GW}^2$| should be positive, otherwise the integrant vanishes. The other integrals are computed numerically. We adopt NBH ∼ 7 × 104 BHs galaxy−1 as a fiducial value, corresponding to the GW event rate |$\mathscr {R}_{\rm GW} \sim 70$| Gpc−3 yr−1 (The LIGO Scientific Collaboration et al. 2016) in equation (2). 3.1 Mass function We assume a Salpeter-like mass function for the primary BH, \begin{eqnarray} \frac{{\rm d}p(m_1)}{{\rm d}m_1} = C m_1^{-\gamma }, \end{eqnarray} (18) with a uniform distribution of the secondary mass, \begin{eqnarray} \frac{{\rm d}p(m_2|m_1)}{{\rm d}m_2} = \frac{1}{m_1-M_{\min }}, \end{eqnarray} (19) where γ = 2.35, Mmin  ≤ m2 ≤ m1 ≤ Mmax , Mmin  = 5 M⊙, Mmax  = 50 M⊙ and |$C=(\gamma -1)/(M_{\min }^{1-\gamma }-M_{\max }^{1-\gamma })$|⁠. Such mass functions are inferred by the observations of massive stars (Sana et al. 2012; Kobulnicky et al. 2014). Similar mass functions5 are adopted by the analysis of LIGO O1 data (The LIGO Scientific Collaboration et al. 2016) and are consistent with the GW observations. Note that the total luminosity is dominated by heavy masses for γ < 3. In this respect, GW 150914 is crucial by raising the maximum mass Mmax  and hence the expected luminosity more than was previously thought. (Cf. Mmax  = 13 M⊙ was adopted in Agol & Kamionkowski 2002. Note also γ = 0.35 in Agol & Kamionkowski 2002.) 3.2 Velocity distribution before a merger The velocity distribution for GW 150914-like BHs before mergers is described by a Maxwell–Boltzmann distribution, \begin{eqnarray} \frac{{\rm d}f(v)}{{\rm d}v}=\sqrt{\frac{2}{\pi }} \frac{v^2}{\sigma _v^3} \exp \left(-\frac{v^2}{2 \sigma _v^2}\right), \end{eqnarray} (20) where an isotropic Gaussian approximation is enough for our order-of-magnitude estimates. As a fiducial value, we take the velocity dispersion σv = 40 km s−1 by considering the isolated binary formation scenario. From a theoretical point of view, massive star progenitors are born from molecular clouds and their velocity dispersion is initially low σv ∼ 10 km s−1 (Binney & Merrifield 1998). Unless the BH formation is associated with an exceptionally large kick due to, e.g. an asymmetric mass ejection, the resulting BHs have also low velocities. If the kick velocity is inversely proportional to the mass following the momentum conservation, the kick velocity of neutron stars implies σv ∼ 200 km s−1(1.4 M⊙/M) ∼ 30 km s−1(M/10 M⊙)−1. Older stars tend to have larger velocity dispersion and σv = 40 km s−1 is reasonable for progenitors with metallicity Z ≲ 0.5 Z⊙ (Binney & Merrifield 1998). From an observational point of view, the rms distance ∼410 pc from the Galactic plane for BH low-mass X-ray binaries, corresponding to a scale height of 290 pc, suggests a velocity dispersion of σv ∼ 40 km s−1 (White & van Paradijs 1996). Although there are exceptions such as GRO 1655−40 with a peculiar velocity v ∼ −114 km s−1 (Brandt, Podsiadlowski & Sigurdsson 1995; Mirabel et al. 2002) and XTE J1118+480 with v ∼ 145 km s−1 (Mirabel et al. 2001), two populations likely exist with low and high kick velocities, similar to neutron stars (Cordes & Chernoff 1998; Pfahl et al. 2002). On the other hand, these observations are not for high-mass systems. In addition, these estimates are subject to systematic errors in the distance (Repetto, Davies & Sigurdsson 2012). The most reliable estimate is based on the astrometric observations (Miller-Jones 2014). Although there is only one sample for a high-mass system, the BH high-mass X-ray binary Cygnus X-1 has a relatively low proper motion ∼20 km s−1 (Chevalier & Ilovaisky 1998; Mirabel & Rodrigues 2003; Reid et al. 2011). We discuss a high-velocity case σv = 200 km s−1 later in Section 5.2. 3.3 ISM density We consider five phases of the ISM as listed in Table 1; the molecular clouds consisting mostly of H2, the cold neutral medium consisting of H i clouds (cold H i), the warm neutral medium in thermally equilibrium with cold H i (warm H i), the warm ionized medium (warm H ii) and the hot ionized medium (hot H ii; Bland-Hawthorn & Reynolds 2000; Ferrière 2001; Heyer & Dame 2015; Inutsuka et al. 2015). For each phase, we use the probability distribution of the number density, \begin{eqnarray} \frac{{\rm d}\xi (n)}{{\rm d}n} = D \xi _0 n^{-\beta },\quad (n_1<n<n_2), \end{eqnarray} (21) where n1, n2 and β are given in Table 1 (Berkhuijsen 1999), |$D=(\beta -1)/(n_1^{1-\beta }-n_2^{1-\beta })$| and |$\xi _0 = \int \frac{{\rm d}\xi }{{\rm d}n} {\rm d}n$| is the volume filling fraction (Scoville & Sanders 1987; Clemens, Sanders & Scoville 1988; Agol & Kamionkowski 2002).6 Each phase has its scale height Hd in the Galactic disc. We assume that the hot H ii phase has a sound velocity cs = 150 km s−1 corresponding to a temperature T ∼ 106 K, while the other phases have effective cs ∼ 10 km s−1 (corresponding to T ∼ 2 × 104 K) because these phases (even with T < 2 × 104 K) have also a turbulent velocity ∼10 km s−1 in approximately pressure balance with each other. More precisely, although the Larson’s law suggests a smaller turbulent velocity in a molecular cloud than ∼10 km s−1, the relative velocity between molecular clouds (and also other phases) is ∼10 km s−1, and this is relevant to an accreting BH because the BH should have left the cradle molecular cloud. The parameters in Table 1 are similar to those in Agol & Kamionkowski (2002). Table 1. Five ISM phases. The density distribution from n1 to n2 with an index β, the volume filling fraction ξ0, the (effective) sound velocity cs (including the turbulent velocity for cold phases) and the scale height of the disc Hd are shown. Phase . n1 (cm−3) . n2 (cm−3) . β . ξ0 . cs (km s−1) . Hd . Molecular clouds 102 105 2.8 10−3 10 75 pc Cold H i 10 102 3.8 0.04 10 150 pc Warm H i 0.3 – – 0.35 10 0.5 kpc Warm H ii 0.15 – – 0.2 10 1 kpc Hot H ii 0.002 – – 0.4 150 3 kpc Phase . n1 (cm−3) . n2 (cm−3) . β . ξ0 . cs (km s−1) . Hd . Molecular clouds 102 105 2.8 10−3 10 75 pc Cold H i 10 102 3.8 0.04 10 150 pc Warm H i 0.3 – – 0.35 10 0.5 kpc Warm H ii 0.15 – – 0.2 10 1 kpc Hot H ii 0.002 – – 0.4 150 3 kpc Open in new tab Table 1. Five ISM phases. The density distribution from n1 to n2 with an index β, the volume filling fraction ξ0, the (effective) sound velocity cs (including the turbulent velocity for cold phases) and the scale height of the disc Hd are shown. Phase . n1 (cm−3) . n2 (cm−3) . β . ξ0 . cs (km s−1) . Hd . Molecular clouds 102 105 2.8 10−3 10 75 pc Cold H i 10 102 3.8 0.04 10 150 pc Warm H i 0.3 – – 0.35 10 0.5 kpc Warm H ii 0.15 – – 0.2 10 1 kpc Hot H ii 0.002 – – 0.4 150 3 kpc Phase . n1 (cm−3) . n2 (cm−3) . β . ξ0 . cs (km s−1) . Hd . Molecular clouds 102 105 2.8 10−3 10 75 pc Cold H i 10 102 3.8 0.04 10 150 pc Warm H i 0.3 – – 0.35 10 0.5 kpc Warm H ii 0.15 – – 0.2 10 1 kpc Hot H ii 0.002 – – 0.4 150 3 kpc Open in new tab 3.4 Scale height The BHs have their scale height H(vz) in the Galactic disc. Each phase of the ISM has also its own scale height Hd in Table 1. Then the number of BHs in each phase is corrected by a factor \begin{eqnarray} h(m_1,m_2,v)=\min \left[1,\frac{H_{\rm d}}{H(v_z)}\right]. \end{eqnarray} (22) For simplicity, we make a one-dimensional analysis of the vertical structure, neglecting the coupling of the vertical and horizontal motions. The scale height H(vz) is determined by the velocity in the z-direction, \begin{eqnarray} \frac{1}{2}v_z^2=\Phi _z\left[H(v_z)\right], \end{eqnarray} (23) where the z-component of the velocity is |$v_z^2=\frac{1}{3}(v^2+v_{\rm GW}^2)$| and the gravitational potential in the z-direction, \begin{eqnarray} \frac{\Phi _z(z)}{2\pi G}=K \big(\sqrt{z^2+Z^2}-Z\big)+Fz^2, \end{eqnarray} (24) where Z ∼ 180 pc, K = 48 M⊙ pc−2 and F = 0.01 M⊙ pc−3 (Kuijken & Gilmore 1989a,b). This simple model is sufficient for our order-of-magnitude estimates. 3.5 GW recoil velocity Merged BHs receive a recoil due to the anisotropic GW emission (Bonnor & Rotenberg 1961; Peres 1962; Bekenstein 1973). GWs carry linear momentum if two merging BHs have different masses and/or finite spins. Fitting formulas for the recoil velocity are obtained by using numerical simulations in the post-Newtonian-inspired forms (Baker et al. 2007; Campanelli et al. 2007; Zlochower & Lousto 2015). The recoil velocity for a merger of non-spinning BHs is well approximated by \begin{eqnarray} v_{\rm GW} &=& A \eta ^2 \sqrt{1-4\eta } \left(1+B \eta \right), \end{eqnarray} (25) where A = 1.20 × 104 km s−1, B = −0.93 (Fitchett 1983; Gonzalez et al. 2007) and η = m1m2/(m1 + m2)2 is the symmetric mass ratio. The maximum value is vGW ∼ 175 km s−1 for η = 0.195. GW 150914 has η ∼ 0.247 and hence vGW ∼ 61 km s−1, and GW 151226 has η ∼ 0.226 and hence vGW ∼ 150 km s−1. Although we have to extrapolate the formula to small η, this does not affect our result so much. We do not consider the spin-induced recoil because we are now assuming that the pre-merger spin is low to make a conservative estimate on the released energy from BHs. If the pre-merger spin is high, the pre-merger BHs, which are much more abundant than the merged BHs, can release energy through the BZ mechanism even before the merger without the GW recoil. This case will be discussed in Section 5.1. Current observations of GWs show that the primary BH has a spin of <0.7 at 90 per cent confidence with no evidence for spins being both large and strongly aligned. For GW 151226, the effective spin parameter is |$0.21^{+0.20}_{-0.10}$| and may imply spinning BHs (The LIGO Scientific Collaboration et al. 2016). For reference, if the spin parameters parallel to the orbital axis are a*2∥ ∼ 0.2 and a*1∥ ∼ 0 with the same masses η ∼ 1/4, the recoil velocity is about vGW ∼ 40 km s−1, while it is vGW ∼ 260 km s−1 if the in-plane spins are a*2⊥ ∼ 0.2 and a*1⊥ ∼ 0 with the same masses η ∼ 1/4 (Zlochower & Lousto 2015). 3.6 Results of the luminosity function Fig. 2 shows the luminosity function of the BH jets from accreting BHs in our Galaxy for the fiducial case (see Table 2 for the other cases), calculated from equation (16). Each line corresponds the ISM phase where the BHs reside. As the accretion rate is proportional to the ISM density in equation (4), the jet luminosity is brighter for BHs in denser ISM such as molecular clouds. On the other hand, brighter jets are rarer because the volume filling fraction of denser medium is smaller in the ISM as in Table 1. We can find that the brightest sources in our Galaxy (with the number |${\rm d}N/{\rm d}\,{\rm log}\,{\dot{M}} \sim 1$|⁠) have Lj ∼ 1036 erg s−1 mostly residing in the cold HI, while fainter sources are more abundant. Figure 2. Open in new tabDownload slide Luminosity function of the BH jets from GW 150914-like merged BHs in our Galaxy for the isolated binary formation scenario with the velocity dispersion σv = 40 km s−1, calculated from equation (16). Each line corresponds the ISM phase where the BHs reside (see Table 1). We can see that the most luminous source in our Galaxy has Lj ∼ 1036 erg s−1 for this scenario. Dotted lines are calculated by setting the GW recoil velocity vGW = 0. The upper horizontal axis is the maximum energy of accelerated particles allowed by the Hillas condition for a given luminosity with the charge of accelerated particles Z = 1 in equation (27). Figure 2. Open in new tabDownload slide Luminosity function of the BH jets from GW 150914-like merged BHs in our Galaxy for the isolated binary formation scenario with the velocity dispersion σv = 40 km s−1, calculated from equation (16). Each line corresponds the ISM phase where the BHs reside (see Table 1). We can see that the most luminous source in our Galaxy has Lj ∼ 1036 erg s−1 for this scenario. Dotted lines are calculated by setting the GW recoil velocity vGW = 0. The upper horizontal axis is the maximum energy of accelerated particles allowed by the Hillas condition for a given luminosity with the charge of accelerated particles Z = 1 in equation (27). Table 2. Possible uncertainties of the model parameters and the resulting total power Ptot are summarized. The isolated binary formation scenario is the fiducial case. The column with ‘–’ equals the fiducial value. The model parameters are the number of BH jets NBH, the dispersion of the velocity distribution σv in equation (20), the initial BH spin |$a_*^{\rm i}$|⁠, the accretion rate profile s in equation (31) and the duty cycle |$\mathscr {D}$| reduced by the feedback on the ISM by the disc outflow. . Number NBH . Velocity σv (km s−1) . Spin |$a_*^{\rm i}$| . Disc s . Duty cycle |$\mathscr {D}$| . Total power Ptot (erg s−1) . Isolated binary (fiducial) 7 × 104 40 0 0 1 ∼1037 High spin 108 – 0.2 – – ∼1037 × 103 Stellar cluster – 200 – – – ∼1037 × 10−2 Wind – – – 1 – ∼1037 × 10−2 Feedback – – – 1 10−1 ∼1037 × 10−3 . Number NBH . Velocity σv (km s−1) . Spin |$a_*^{\rm i}$| . Disc s . Duty cycle |$\mathscr {D}$| . Total power Ptot (erg s−1) . Isolated binary (fiducial) 7 × 104 40 0 0 1 ∼1037 High spin 108 – 0.2 – – ∼1037 × 103 Stellar cluster – 200 – – – ∼1037 × 10−2 Wind – – – 1 – ∼1037 × 10−2 Feedback – – – 1 10−1 ∼1037 × 10−3 Open in new tab Table 2. Possible uncertainties of the model parameters and the resulting total power Ptot are summarized. The isolated binary formation scenario is the fiducial case. The column with ‘–’ equals the fiducial value. The model parameters are the number of BH jets NBH, the dispersion of the velocity distribution σv in equation (20), the initial BH spin |$a_*^{\rm i}$|⁠, the accretion rate profile s in equation (31) and the duty cycle |$\mathscr {D}$| reduced by the feedback on the ISM by the disc outflow. . Number NBH . Velocity σv (km s−1) . Spin |$a_*^{\rm i}$| . Disc s . Duty cycle |$\mathscr {D}$| . Total power Ptot (erg s−1) . Isolated binary (fiducial) 7 × 104 40 0 0 1 ∼1037 High spin 108 – 0.2 – – ∼1037 × 103 Stellar cluster – 200 – – – ∼1037 × 10−2 Wind – – – 1 – ∼1037 × 10−2 Feedback – – – 1 10−1 ∼1037 × 10−3 . Number NBH . Velocity σv (km s−1) . Spin |$a_*^{\rm i}$| . Disc s . Duty cycle |$\mathscr {D}$| . Total power Ptot (erg s−1) . Isolated binary (fiducial) 7 × 104 40 0 0 1 ∼1037 High spin 108 – 0.2 – – ∼1037 × 103 Stellar cluster – 200 – – – ∼1037 × 10−2 Wind – – – 1 – ∼1037 × 10−2 Feedback – – – 1 10−1 ∼1037 × 10−3 Open in new tab The GW recoil effect reduces the luminosity by an order-of-magnitude as we can see from the dotted lines, which are calculated by setting vGW = 0. This reduction is approximately determined by the v−3 dependence of the accretion rate in equation (4) as ∼(vGW/σv)−3 ∼ (100/40 km s−1)−3 ∼ 0.06. Note also that the luminosity function for the hot H ii phase has a peak because this phase has a large cs and so v ∼ cs has little dispersion. Fig. 3 is the same as Fig. 2 but with the vertical axis multiplied by |$L_{{\rm j}} \sim {\dot{M}}c^2$|⁠. This makes it clear that the most energy is generated by BHs in the cold H i medium. We can read the total power, \begin{eqnarray} P_{\rm tot}=\int L_{{\rm j}} \frac{{\rm d}N}{{\rm d}{\dot{M}}} {\rm d}{\dot{M}} \sim 10^{37}\, {\rm erg}\, {\rm s}^{-1} \left(\frac{N_{\rm BH}}{7 \times 10^{4}}\right), \end{eqnarray} (26) which is very roughly derived by |$P_{\rm tot} \sim N_{{\rm BH}} \times \xi _0^{{\rm H}{\rm {\small {I}}}} \times {\dot{M}}(10\,{\rm cm}^{-3},\ 5 \,\mathrm{M}_{{\odot }},\ 5 \,\mathrm{M}_{{\odot }},\ 100\,{\rm km}\,{\rm s}^{-1}) c^2 \times (M_{\max }/M_{\min })^{3-\gamma } \sim 7 \times 10^{4} \times 0.04 \times 5 \times 10^{32}\,{\rm erg}\,{\rm s}^{-1} \times 4.5 \sim 0.6 \times 10^{37}\,{\rm erg}\,{\rm s}^{-1}$|⁠. Note that the velocity dependencies of |$\dot{M}$| and f(v) cancel each other and the low-velocity BHs have smaller scale height than the cold H i. The total power is approximately ∼3 × 10−5 of that of SN explosions ESN/100 yr ∼ 3 × 1041 erg s−1. This is small but comparable to the required power for some high-energy particles in our Galaxy. Based on these results, we will discuss observational implications in the next section. Figure 3. Open in new tabDownload slide Same as Fig. 2 but with the vertical axis multiplied by |$L_{{\rm j}} \sim {\dot{M}}c^2$|⁠. We can find that the most energy is generated by BHs in the cold H i medium. Dashed line plots the required power spectrum for supplying the observed CRs. We can see that the total power is comparable to that of CRs above the knee energy at around ∼3 PeV within the model uncertainties that are listed in Table 2. Figure 3. Open in new tabDownload slide Same as Fig. 2 but with the vertical axis multiplied by |$L_{{\rm j}} \sim {\dot{M}}c^2$|⁠. We can find that the most energy is generated by BHs in the cold H i medium. Dashed line plots the required power spectrum for supplying the observed CRs. We can see that the total power is comparable to that of CRs above the knee energy at around ∼3 PeV within the model uncertainties that are listed in Table 2. 4 OBSERVATIONAL IMPLICATIONS 4.1 PeVatrons Particle acceleration is ubiquitous in the BH jet system as manifested in AGNs, GRBs and X-ray binaries (Longair 2011). The maximum acceleration energy is limited by the source size, i.e. the so-called Hillas condition, εmax  = ZqBr/Γ, where Z is the charge of accelerated particles, Γ is the Lorentz factor of the acceleration region and B is the lab-frame magnetic field. This can be written in terms of the Poynting luminosity |$L_{{\rm j}} \sim 2 \pi \theta _{\rm j}^2 r^{2} c (B^2/8\pi )$|⁠, where the magnetic field carries an energy density B2/8π at a radius r with the jet opening angle θj (Norman, Melrose & Achterberg 1995; Blandford 2000; Waxman 2004). With the causality condition Γθj ≳ 1, we have \begin{eqnarray} \varepsilon _{\max }=Zq B \frac{r}{\Gamma }\sim\! \frac{2Zq}{\Gamma \theta _{\rm j}} \sqrt{\frac{L_{{\rm j}}}{c}}\! \lesssim 3.5\, {\rm PeV}\, Z\! \left(\frac{L_{{\rm j}}}{10^{36}\,{\rm erg}\,{\rm s}^{-1}}\right)^{\!1/2}\!\!. \end{eqnarray} (27) Therefore bright sources, i.e. massive BHs in dense ISM, are potential ‘PeVatrons’ accelerating particles beyond PeV energy (Barkov, Khangulyan & Popov 2012; Kotera & Silk 2016). We plot εmax  on the upper horizontal axis in Figs 2 and 3. Possible acceleration sites are discussed in Section 7. In Fig. 3, we also plot the required power spectrum for supplying the observed CRs by using |$L_0 (\varepsilon /\varepsilon _{\min }^{\rm CR})^{2-s} \frac{{\rm d}\varepsilon }{\varepsilon } =L_0 (\varepsilon /\varepsilon _{\min }^{\rm CR})^{2-s} \frac{{\rm d}{\dot{M}}}{2{\dot{M}}}$|⁠, where the spectral index is s = 2.34 below the knee (Genolini et al. 2015)7 and 0.3 higher above the knee (Blümer, Engel & Hörandel 2009; Gaisser, Stanev & Tilav 2013). SN remnants are commonly believed to supply most Galactic CRs from the peak energy |$\varepsilon _{\min }^{\rm SN}=1$| GeV to the knee |$\varepsilon _{\max }^{\rm SN}=3$| PeV. The normalization L0 is determined by the fact that a fraction εCR = 0.1 of the SN kinetic energy can yield CRs, \begin{eqnarray} \frac{\epsilon _{\rm CR} E_{\rm SN}}{100\,{\rm yr}} =\int _{\varepsilon _{\min }^{\rm SN}}^{\varepsilon _{\max }^{\rm SN}} L_0 \left(\frac{\varepsilon }{\varepsilon _{\min }^{\rm SN}}\right)^{2-s} \frac{{\rm d}\varepsilon }{\varepsilon }. \end{eqnarray} (28) From Fig. 3 we can see that the BH jets can produce comparable energy to that required for the observed CRs at and beyond the knee energies ≳3 PeV, taking the model uncertainties into account (see Section 5). The origin of these CRs is not known (e.g. Hillas 2005; Blümer et al. 2009; Gaisser et al. 2013). Currently known gamma-ray sources, including even SN remnants, do not show the characteristic PeVatron spectrum extending without a cut-off or break to tens of TeV (Aharonian 2013), with a possible exception of the Galactic Centre Sagittarius A* (HESS Collaboration et al. 2016). Even if the SN remnants are responsible for CRs up to the knee, the transition from Galactic to extragalactic CRs occurs between the knee and the ankle. Ultrahigh-energy CRs above the ankle are extragalactic because of the observed isotropy (Abreu et al. 2010; Abbasi et al. 2017).8 If the knee corresponds to the proton cut-off and the source composition is solar, the rigidity-dependent cut-offs extending beyond the knee are not sufficient to fill the observed all-particle flux (Hillas 2005). This implies a second (Galactic) component at energies between the knee and the ankle, sometimes called ‘component B’. Our results suggest that the BH jets might be PeVatrons and/or fill the gap between the knee and the ankle. An unnatural point of this possibility is that the BH jets are totally irrelevant to the SN remnants. It is just a coincidence that the CR fluxes from two kinds of sources are the same within a factor. There are also orders-of-magnitude uncertainties in the estimate of the total power of the BH jets (see Section 5). Furthermore, it is difficult to calculate the CR spectrum and the acceleration efficiency at present. Nevertheless the BH jets can potentially accelerate the CRs at and beyond PeV energy with the flux comparable to the observations. 4.2 Cosmic ray positrons and electrons The CR positron fraction (the ratio of positrons to electrons plus positrons) has been measured by the PAMELA satellite (Adriani et al. 2009) and more precisely by the AMS-02 experiment (Aguilar et al. 2013). The observed positron fraction rises from ∼10 GeV at least to ∼300 GeV, indicating the presence of nearby positron sources within ∼1 kpc. Although the dark matter annihilation or decay scenario is now severely constrained by other messengers, there are still many astrophysical candidates and the true origin is unclear (e.g. Fujita et al. 2009; Ioka 2010; Kashiyama, Ioka & Kawanaka 2011; Serpico 2012; Kohri et al. 2016). The BH jets could accelerate electrons and positrons preferentially if the jets are not contaminated by baryons (Barkov et al. 2012) but associated with the pair cascade (see also Section 7). The maximum energy of particle acceleration is enough for producing the observed positrons as in equation (27). The required total power for the positron excess is about ∼10−4 of that of SN explosions, i.e. about ∼3 × 1037 erg s−1. This is comparable to that of the BH jets in equation (26) within the model uncertainties (see Section 5). Therefore the BH jets are eligible to join the possible sources, although it is again difficult to estimate the spectrum and efficiency of the positron acceleration. A BH jet likely forms an extended nebula in the ISM, similarly to an AGN cocoon/lobe (Begelman & Cioffi 1989) and a GRB cocoon (Ramirez-Ruiz, Celotti & Rees 2002; Mizuta & Ioka 2013). A BH jet collides with the ISM at the jet head. The shocked matter goes sideways forming a cocoon. Although the cocoon pressure initially collimates the jet, the collimation radius expands and finally reaches the termination (reverse) shock. The maximum size of the termination shock is given by the condition that the jet pressure balances with the ram pressure of the ISM, |$L_{{\rm j}}/2 \pi \theta ^2 r_{\rm h}^2 c \sim n \mu m_{\rm u} V^2$|⁠, i.e. \begin{eqnarray} r_{\rm h} &\sim & \sqrt{\frac{L_{{\rm j}}}{2\pi \theta ^2 c n \mu m_{\rm u} V^2}} \sim 2\, {\rm pc} \left(\frac{L_{{\rm j}}}{10^{36}\,{\rm erg}\,{\rm s}^{-1}}\right)^{1/2} \nonumber \\ &&\times \left(\frac{n}{10\,{\rm cm}^{-3}}\right)^{-1/2} \left(\frac{V}{10\,{\rm km}\,{\rm s}^{-1}}\right)^{-1} \left(\frac{\theta }{0.1}\right)^{-1}, \end{eqnarray} (29) at which point the jet is completely bent by the ISM and dissipated into the cocoon. The cocoon is also extended along the direction of the proper motion, leading to a more or less spherical shape. The forward shock of the BH jet nebula expands with a velocity vc ∼ (Lj/nμmu)1/5t−2/5 and a size rc ∼ vct, slowing down to vc ∼ v at the maximum size rc,max  ∼ (Lj/nmuv3)1/2 ∼ 80 pc (Lj/1036 erg s−1)1/2(n/10 cm−3)−1/2(v/10 km s−1)−3/2. The BH jet nebula is similar to an old pulsar wind nebula (PWN). Then the CR electrons and positrons likely escape from the nebula to the ISM without adiabatic cooling (Kashiyama et al. 2011). Radiative cooling such as synchrotron emission limits the maximum energy of electrons and positrons, which depends on the propagation time and the magnetic field in the nebula (Kawanaka, Ioka & Nojiri 2010). Future observations beyond TeV energies by Calorimeteric Electron Telescope (CALET), Dark Matter Particle Explorer (DAMPE) and Cherenkov Telescope Array (CTA) will probe such leptonic PeVatrons (Kobayashi et al. 2004; Kawanaka et al. 2011). 4.3 Unidentified TeV gamma-ray sources The Galactic plane survey carried out by High Energy Stereoscopic System (H.E.S.S.) led to the discovery of dozens of TeV gamma-ray sources. Among these, the most abundant category is dark accelerators, so-called TeV unidentified sources (TeV unIDs), which have no clear counterpart at other wavelengths (Aharonian et al. 2005, 2006, 2008). They lie close to the Galactic plane, suggesting Galactic sources. Their power-law spectra with an index of 2.1–2.5 imply a connection with CR accelerators. They are extended ΔΘ ∼ 0|${^{\circ}_{.}}$|05–0|${^{\circ}_{.}}$|3, corresponding to a physical size of ∼3 pc (ΔΘ/0|${^{\circ}_{.}}$|2)(D/kpc) for an unknown distance D. Still, their unID nature prevents us to identify their origin (e.g. Yamazaki et al. 2006; de Jager et al. 2009; Ioka & Mészáros 2010). In Fig. 4, we plot the observed flux distribution of the TeV unIDs at energies εγ > 0.2 TeV in terms of the cumulative number of sources above a flux N(> F), i.e. log N–log F plot. In order to compare it with the BH jets, we calculate the flux distribution from the luminosity function in equation (16) by integrating the number of sources above a given (bolometric) flux |$F = L_{{\rm j}}/4\pi D^2 \sim {\dot{M}} c^2/4\pi D^2$| as \begin{eqnarray} \frac{{\rm d}N}{{\rm d}F}=\int \frac{{\rm d}N}{{\rm d}{\dot{M}}} \frac{4\pi D^2}{c^2} \frac{D{\rm d}D {\rm d}\theta }{\pi R_{\rm d}^2}, \end{eqnarray} (30) where we approximate the spatial distribution of the BH jets by a thin uniform disc with a radius Rd = 15 kpc and the distance of the Sun to the Galactic Centre R⊙ = 8 kpc. A thin approximation is applicable if the observed distance is larger than the scale height ∼300 pc for the fiducial case (see Table 2). Figure 4. Open in new tabDownload slide Flux distribution expressed by the cumulative number of sources above a given flux, i.e. log N(> F)–log F plot, with equations (30) and (16). This is compared with the observations of TeV unIDs. Both are comparable if the gamma-ray efficiency is εγ ∼ 10−2. CTA could detect additional ∼300 BH jets in the near future. Figure 4. Open in new tabDownload slide Flux distribution expressed by the cumulative number of sources above a given flux, i.e. log N(> F)–log F plot, with equations (30) and (16). This is compared with the observations of TeV unIDs. Both are comparable if the gamma-ray efficiency is εγ ∼ 10−2. CTA could detect additional ∼300 BH jets in the near future. Fig. 4 shows that the flux distribution is comparable with that of TeV unIDs if the gamma-ray efficiency is about εγ ∼ 10−2 for the fiducial parameters (see Table 2). Note that the inverse Compton cooling time of 10 TeV electrons is ∼105 yr. If the age is ∼105 yr, the TeV gamma-ray flux is ∼0.1–0.02 Le, implying that εe ∼ 0.1–0.5. This is comparable with values considered in GRB jets and PWNe. If this is the case, the CTA observatory will increase the number of TeV unIDs up to ∼300 by improving the sensitivity by about an order of magnitude in the near future (Acharya et al. 2013). Note that the flux distribution follows N(> F) ∝ D2 ∝ F−1 if the spatial distribution is disc-like, which is different from N(> F) ∝ F−1.5 for the 3D Euclidian space. The uniform disc approximation is acceptable for the current observations, which have not reached the Galactic Centre yet. For future observations, we have to consider the high-density region near the Galactic Centre. The nebular size in equation (29) is also consistent with the extended nature of TeV unIDs. The BH jet nebula also evades strong upper limits in X-rays with a TeV to X-ray flux ratio up to ≳50 (Bamba et al. 2007, 2009; Matsumoto et al. 2007, 2008; Fujinaga et al. 2011; Sakai, Yajima & Matsumoto 2011). This is because the physical parameters such as the energy density and the magnetic field are similar to those of an old PWN. Their emission spectra have the unID nature thanks to the old age (Yamazaki et al. 2006; de Jager et al. 2009; Ioka & Mészáros 2010). In addition, the ADAF disc is radiatively inefficient. The X-ray flux of the ADAF disc is about |$F_{\nu } \sim (\alpha _{\rm QED}/\alpha ^2) ({m_{\rm e}}/{m_{\rm u}}) ({{\dot{M}}c^2}/{L_{\rm Edd}}) {\dot{M}}/4\pi D^2 m_{\rm e} \sim 1 \times 10^{-18}\,{\rm erg}\,{\rm s}^{-1}\,{\rm cm}^{-2}\,{\rm keV}^{-1}$| (α/0.1)−2 (n/10 cm−3)2 (M/10 M⊙)3 (v/10 km s−1)−6 (D/kpc)−2, below the current limit. 5 MODEL UNCERTAINTIES Although the GW observations significantly narrow down the possible parameter space, in particular putting a lower bound on the number of spinning BHs in equation (2), there are still large uncertainties about the model parameters and the estimate for the BH jet power. In this section, we clarify the range of the uncertainties by considering four representative effects: the initial BH spin (Section 5.1), the velocity distribution depending on the binary BH formation scenario (Section 5.2), the accretion rate profile changed by the disc wind (Section 5.3) and the feedback on the ISM by the outflow (Section 5.4). These effects on the model parameters and the resulting total power are summarized in Table 2. We enclose the uncertainty of the total power for the BH jets within a factor of 10±3, which is much better than before. 5.1 Initial spin We have discussed only the merged BH population because the merged BHs are undoubtedly spinning fast by bringing the orbital angular momentum at the merger, while the pre-merger BHs might not. However, if the BHs have spins before the mergers, the BHs can launch BZ jets without the mergers. Such spinning BHs could result from the massive stellar collapse. The total number of BHs in our Galaxy is about NBH ∼ 108 (Shapiro & Teukolsky 1983), ∼103 times larger than that of the merged BHs in equation (2).9 In addition the GW recoil is absent without a merger, increasing the total power by a factor of 10 as shown in Fig. 3. Then the total power is larger than the fiducial value by a factor of |${\sim } 10^{4} (a_*^{\rm i}/0.7)^2$| altogether, i.e. Ptot ∼ 1041 erg |$(a_*^{\rm i}/0.7)^2$|⁠, where the |$(a_*^{{\rm i}})^2$| dependence comes from that of the BZ luminosity in equation (14). In other words, the total power dominates the fiducial one if the spin parameter is larger than a* > 7 × 10−3. GW observations show no evidence for large spins. Probably the initial spin would be small for most BHs because the massive star progenitors with solar metallicity lose the angular momentum by stellar wind (Heger et al. 2003; Hirschi, Meynet & Maeder 2005). Because of the same reason, the BH mass is also smaller than the fiducial case (Abbott et al. 2016c), reducing the total power of the BH jets. In low metallicity, the wind is weak and the resulting BH spin may be high (Hirschi, Meynet & Maeder 2005; Yoon & Langer 2005; Kinugawa, Nakano & Nakamura 2016a). A rapid rotation of the progenitors could lead to a chemically homogeneous evolution without a common envelope phase, avoiding a merger before the BH formation (Mandel & de Mink 2016; Marchant et al. 2016). However, the number of such BHs is less than NBH ∼ 108. The event rate of GRBs, which likely produce spinning BHs, is comparable to that of the BH mergers. Although some BHs in X-ray binaries might have high spins, these measurements are subject to systematic errors (Remillard & McClintock 2006). For GW 151226, the effective spin parameter is |$0.21^{+0.20}_{-0.10}$| (The LIGO Scientific Collaboration et al. 2016). So we tentatively take |$a_*^{\rm i} \sim 0.2$| as an upper limit in Table 2. This is the most extreme case because the total power is comparable to that of SN explosions, ESN/100 yr ∼ 3 × 1041 erg s−1. 5.2 Binary BH formation scenario The accretion rate and the resulting jet luminosity sensitively depend on the velocity of the BH in equation (4). We have adopted σv = 40 km s−1 as a fiducial value for the isolated binary formation scenario (Tutukov & Yungelson 1993; Dominik et al. 2015; Belczynski et al. 2016; Lipunov et al. 2017) in equation (20). The GW 150914 masses favour low metallicity below 0.5 Z⊙ (Abbott et al. 2016c). The extreme case is zero metallicity Population III stars (Kinugawa et al. 2014, 2016b; Kinugawa, Nakano & Nakamura 2016a,c). If BHs form in very low metallicity <0.01 Z⊙, the GW events may be dominated by recent BH mergers in dwarf galaxies (Lamberts et al. 2016) because the low metallicity allows a small initial separation of a BH binary. Then the merged, spinning BHs are incorporated into our Galaxy relatively recently, joining in the halo component with a velocity dispersion of σv ∼ 200 km s−1. Another scenario is the dynamical binary formation in a dense stellar cluster (Kulkarni, Hut & McMillan 1993; Sigurdsson & Hernquist 1993; Portegies Zwart & McMillan 2000; Rodriguez et al. 2015; Mapelli 2016; Rodriguez, Chatterjee & Rasio 2016). In a high-density stellar environment, BHs dynamically interact and form binaries. Since the interaction is frequent in the clusters, most of the BH mergers may occur outside the clusters following dynamical ejection. The escape velocity of the clusters is smaller than that of our Galaxy. Thus the merged BHs are floating in our halo with a velocity dispersion of σv ∼ 200 km s−1. Primordial BHs are also a possible candidate (e.g. Nakamura et al. 1997; Ioka et al. 1998; Ioka, Tanaka & Nakamura 1999; Bird et al. 2016; Sasaki et al. 2016), although this scenario requires a fine tuning in the primordial density fluctuation. In this case, the BHs reside in our halo with σv ∼ 200 km s−1. Fig. 5 shows the case of σv ∼ 200 km s−1. Compared with the fiducial case σv ∼ 40 km s−1 (grey dashed line), the luminosity and hence the total power are reduced by a factor of ∼102. This factor is roughly given by the velocity dependence of the accretion rate, ∼(40/200 km s−1)3 ∼ 0.008. The GW recoil effect becomes less significant than the fiducial case because the velocity dispersion is comparable with the recoil velocity. Figure 5. Open in new tabDownload slide Same as the fiducial case in Fig. 2 except for the dispersion of the velocity distribution σv = 200 km s−1. The total luminosity function for the fiducial case σv = 40 km s−1 is also plotted by a grey dashed line. Figure 5. Open in new tabDownload slide Same as the fiducial case in Fig. 2 except for the dispersion of the velocity distribution σv = 200 km s−1. The total luminosity function for the fiducial case σv = 40 km s−1 is also plotted by a grey dashed line. 5.3 Wind It remains highly uncertain how much of the accreting matter at the Bondi radius reaches the BH (Yuan & Narayan 2014). Some supermassive BH systems with jets seem to require the Bondi accretion rates calculated from the observed gas temperature and density to power the observed jets (Allen et al. 2006; Rafferty et al. 2006; Russell et al. 2013). On the other hand, the ADAF model implies positive Bernoulli parameter for the inflow in the self-similar regime, which suggests that hot accretion flows could have outflows (Narayan & Yi 1994, 1995). The mass outflows make the accretion profile decrease inward approximately in a power-law form, |${\dot{M}}(r)={\dot{M}}(r_{\rm disc}) \left({r}/{r_{\rm disc}}\right)^{s}$|⁠, as in the adiabatic inflow–outflow solutions (ADIOS) model (Blandford & Begelman 1999, 2004). The index is limited to 0 ≤ s < 1 by the mass and energy conservation, but has not been determined yet (Yuan & Narayan 2014). The least accretion case corresponds to s ≈ 1. Recent 3D general relativistic magnetohydrodynamics (MHD) simulations suggest that s ≈ 1 continues down to 20rS, below which the mass flux is constant s = 0 (Yuan et al. 2015). This is also implied by an analytical study (Begelman 2012). If we adopt this least accretion case, the accretion rate of the BH is given by \begin{eqnarray} {\dot{M}}_{\rm BH} \approx {\dot{M}} \left(\frac{20 r_{\rm S}}{r_{\rm disc}}\right)^{s} \quad {\rm if}\ r_{\rm disc}>20r_{\rm S}, \end{eqnarray} (31) with s ≈ 1 where the disc radius rdisc is given by equation (8). Correspondingly, the luminosity of the BH jet is reduced by the same factor (20rS/rdisc)s. Fig. 6 shows the luminosity function using the accretion rate of a BH in equation (31). Compared with the fiducial case s = 0 (grey dashed line), the luminosity and the total power are reduced by a factor of 100. This factor is roughly given by the ratio rdisc/20rS in equation (8). The GW recoil effect becomes less significant than the fiducial case because the disc radius rdisc and the accretion rate |${\dot{M}}$| have similar dependencies on the velocity. Figure 6. Open in new tabDownload slide Same as the fiducial case in Fig. 2 except for the accretion rate of the BH |${\dot{M}}_{\rm BH}={\dot{M}} (20 r_{\rm S}/r_{\rm disc})^{s}$| with s = 1, which is reduced by the wind. The total luminosity function for the fiducial case s = 0 is also plotted by a grey dashed line. Figure 6. Open in new tabDownload slide Same as the fiducial case in Fig. 2 except for the accretion rate of the BH |${\dot{M}}_{\rm BH}={\dot{M}} (20 r_{\rm S}/r_{\rm disc})^{s}$| with s = 1, which is reduced by the wind. The total luminosity function for the fiducial case s = 0 is also plotted by a grey dashed line. 5.4 Feedback Feedback from radiation, jets and winds on the surrounding ISM could be crucial for estimating the total power of the BH jets, as frequently argued in the context of supermassive BHs (Yuan & Narayan 2014). In the Galactic BH case, the radiative feedback is weak because the ADAF disc is much fainter than the Eddington luminosity (e.g. Milosavljević et al. 2009). The radiation may ionize the ISM around the Bondi radius, but once ionized, the cross-section for the interaction between the ISM and photons decreases by many orders of magnitude, reducing the radiative feedback. The radiation does not heat the ISM so much, so that the radiative feedback is negligible except for extreme cases. The jet feedback is also not strong because, although the jet dominates the energy output, its penetration ability makes the dissipation scale large beyond the Bondi radius as shown in equation (29). A large amount of ISM is capable of radiating the injected energy. The most influential feedback would be due to the wind from the disc, if it exists. If the wind is efficient with s ≈ 1 in equation (31), even a small efficiency of the wind feedback εw ≳ 10−6(M/10 M⊙)2/3(v/10 km s−1)−4/3 is able to heat the ISM at the Bondi radius to blow away, |$\epsilon _{\rm w} {\dot{M}}_{\rm BH} c^2 > {\dot{M}} V^2$|⁠. A larger efficiency εw ∼ 0.03–0.001 is implied by simulations (Ohsuga & Mineshige 2011; Yuan et al. 2015; Sa̧dowski et al. 2016). On the other hand, the wind will stop if the mass accretion at the Bondi radius is terminated. Therefore we expect that the BH activity is intermittent with some duty cycle |$\mathscr {D}$| if the wind feedback exists. A rough estimate of the duty cycle is as follows. The wind is somewhat collimated initially when it is released from the disc. If it were spherical, the wind would not be launched because the ram pressure of the Bondi accretion on to the disc exceeds that of the wind. The 4π solid angle of the ISM is affected by the wind after the wind is decelerated by the ISM, which will happen outside the Bondi radius because the ram pressure of the accretion is a decreasing function of the radius and the wind goes straight inside the Bondi radius. Thus, the accretion continues at least for the dynamical time at the Bondi radius, \begin{eqnarray} t_{\rm a} = \frac{r_{\rm B}}{V} \sim 40\, {\rm yr} \left(\frac{M}{10 \,\mathrm{M}_{{\odot }}}\right) \left(\frac{V}{10\,{\rm km}\,{\rm s}^{-1}}\right)^{-3}. \end{eqnarray} (32) The injected energy during the active time is about \begin{eqnarray} E_{{\rm a}}&=&\epsilon _{\rm w} {\dot{M}}_{\rm BH} c^2 t_{\rm a} \sim 3 \times 10^{40}\, {\rm erg} \left(\frac{\epsilon _{\rm w}}{3\hbox{\,per\,cent}}\right) \left(\frac{n}{10\,{\rm cm}^{-3}}\right) \nonumber \\ &&\times \left(\frac{M}{10 \,\mathrm{M}_{{\odot }}}\right)^{11/3} \left(\frac{V}{10\,{\rm km}\,{\rm s}^{-1}}\right)^{-8/3}. \end{eqnarray} (33) The injected energy produces a wind remnant and the lifetime of the remnant is about |$t_{\rm merge} \sim t_{\rm PDS} \times 153 (E_{51}^{1/14} n_0^{1/7}/V_6)^{10/7} \sim 2\times 10^{6}$| yr |$E_{51}^{3/14} n_0^{-4/7} (E_{51}^{1/14} n_0^{1/7}/V_6)^{10/7}$| according to the notation in Cioffi, McKee & Bertschinger (1988), which gives \begin{eqnarray} t_{\rm merge} &\sim & 300\, {\rm yr} \left(\frac{E_{\rm a}}{10^{40}\,{\rm erg}}\right)^{31/98} \left(\frac{n}{10\,{\rm cm}^{-3}}\right)^{-18/49} \nonumber \\ && \times \left(\frac{V}{10\,{\rm km}\,{\rm s}^{-1}}\right)^{-10/7}. \end{eqnarray} (34) Therefore the duty cycle is roughly \begin{eqnarray} \mathscr {D} &\sim & \frac{t_{\rm a}}{t_{\rm merge}} \sim 0.1 \left(\frac{V}{10\,{\rm km}\,{\rm s}^{-1}}\right)^{-0.73} \nonumber \\ && \times \left(\frac{n}{10\,{\rm cm}^{-3}}\right)^{-0.051} \left(\frac{M}{10\,\mathrm{M}_{{\odot }}}\right)^{-0.16}. \end{eqnarray} (35) We use |$\mathscr {D} \sim 10^{-1}$| in Table 2. 6 On Fermi GBM events associated with GW 150914 The GBM on board the Fermi satellite reported a 1-s-lasting weak GRB 0.4 s after GW 150914. Assuming the redshift of GW 150914, |$z=0.09^{+0.03}_{-0.04}$|⁠, the luminosity in 1 keV–10 MeV is |$1.8^{+1.5}_{-1.0} \times 10^{49}$| erg s−1 (Connaughton et al. 2016). This was unexpected and prompted many theoretical speculations (Cardoso et al. 2016; Kimura, Takahashi & Toma 2017; Li et al. 2016; Loeb 2016; Perna, Lazzati & Giacomazzo 2016; Veres et al. 2016; Zhang 2016). The anticoincidence shield (ACS) of the spectrometer on board INTEGRAL (SPI) put upper limits on the gamma-ray emission with similar fluxes (Savchenko et al. 2016). The GBM result also depends on the analysis of low count statistics (Greiner et al. 2016). No counterpart is observed for GW 151226 and LVT 151012 (Racusin et al. 2017). Future follow-ups would be finally necessary to confirm or defeat the GBM detection (Morsony, Workman & Ryan 2016; Murase et al. 2016; Yamazaki, Asano & Ohira 2016). If the signal were caused by the merged BH, the BH would be surrounded by matter. The size of the matter distribution is rm ∼ 1 × 108 cm so that the accretion time is \begin{eqnarray} t_{\rm acc}&=&\frac{1}{\alpha \Omega _{\rm K}} \left(\frac{r_{{\rm m}}}{H}\right)^2 \sim 1.4\,{\rm s} \left(\frac{\alpha }{0.1}\right)^{-1} \left(\frac{M}{60 \,\mathrm{M}_{{\odot }}}\right)^{-1/2} \nonumber \\ &&\times \left(\frac{r_{\rm m}}{1 \times 10^{8}\,{\rm cm}}\right)^{3/2} \left(\frac{H/r_{\rm m}}{0.3}\right)^{-2}, \end{eqnarray} (36) where α is the viscosity parameter, H is the disc scale height and |$\Omega _{\rm K}=\sqrt{GM/r_{\rm m}^3}$| is the Kepler rotation frequency. The mass of the matter should be larger than |$M_{{\rm m}} \gtrsim 10^{-5} \theta _{\rm j}^2 \,\mathrm{M}_{{\odot }}$|⁠, where θj is the opening angle of the GRB jet. The accretion from the ISM affects the matter surrounding a BH. In particular, it can evaporate a possible dead disc that was invoked for the GBM event (Kimura et al. 2017; Perna et al. 2016). A dead disc is assumed to be cold and neutral due to the small mass, suppressing the magnetorotational instability and hence the viscosity, and remains unaccreted, keeping matter for producing the gamma-ray event. However, the accretion from the ISM forms a hot disc sandwiching the dead disc and heating its surface.10 The surface temperature develops a gradient, being greater than T ≳ 104 K (corresponding to the sound velocity |$v_{\rm i} \sim \sqrt{k_{\rm B} T/m_{\rm u}} \sim 10$| km s−1) for the ionized atmosphere. The density ni at the base of the ionized atmosphere is determined by the pressure balance |$n_{\rm i} v_{\rm i}^{2} \sim n(r) v_{\rm h}^2$|⁠, where n(r) ∼ n(rB/r)3/2 and |$v_{\rm h} \sim \sqrt{GM/r}$| are the density and the thermal velocity of the hot disc. Given ni and vi, we can estimate the mass evaporation rate (cf. Hollenbach et al. 1994) from the dead disc as \begin{eqnarray} {\dot{M}}_{\rm eva} &\sim & 2\pi r^2 n_{\rm i} v_{\rm i} m_{\rm u} \sim 1\times 10^{15}\,{\rm g}\,{\rm s}^{-1} \left(\frac{r}{10^{12}\,{\rm cm}}\right)^{-5/2} \nonumber \\ && \times \left(\frac{M}{60 \,\mathrm{M}_{{\odot }}}\right)^{4} \left(\frac{n}{1\,{\rm cm}^{-3}}\right)^{5/2} \left(\frac{V}{40\,{\rm km}\,{\rm s}^{-1}}\right)^{-9/2}. \end{eqnarray} (37) Then the evaporation time of the dead disc with mass Mm is \begin{eqnarray} t_{\rm eva} &\sim & 10^{6}\,{\rm yr} \left(\frac{M_{{\rm m}}}{10^{-5} \,\mathrm{M}_{{\odot }}}\right) \left(\frac{r}{10^{12}\,{\rm cm}}\right)^{5/2} \left(\frac{M}{60 \,\mathrm{M}_{{\odot }}}\right)^{-4} \nonumber \\ && \times \left(\frac{n}{1\,{\rm cm}^{-3}}\right)^{-5/2} \left(\frac{V}{40\,{\rm km}\,{\rm s}^{-1}}\right)^{9/2}, \end{eqnarray} (38) which is shorter than ∼1010 yr, the merger time of the BH binary with a separation r ∼ 1012 cm for the fiducial case. One should keep in mind that the above equation is rather sensitive to parameters M, n and v. For example, BH binaries could have a dead disc if they are formed in a low-density environment. However, for typical parameters, the merged BH would not have a dead disc, implying that the GBM event is not related with GW 150914 in the dead disc scenario. We can also make a second argument that a time reversal of this event seems to encounter physical difficulty. Let us go back in time, say tb ∼ 1000 s before the merger. Still the two BHs should be surrounded by the matter. The size of the matter distribution should be larger, rm ∼ 1010 cm (tb/103 s)2/3(α/0.1)2/3(M/60 M⊙)1/3(H/rm/0.3)4/3, otherwise the matter is swallowed by the BHs before the merger. The bounding energy of this matter is only a fraction of the rest mass energy of the matter, \begin{eqnarray} \frac{G M M_{{\rm m}}/r_{{\rm m}}}{M_{\rm m} c^2} &\sim & 10^{-3} \left(\frac{t_{\rm b}}{10^{3}\,{\rm s}}\right)^{-2/3} \left(\frac{\alpha }{0.1}\right)^{-2/3} \nonumber \\ && \times \left(\frac{M}{60 \,\mathrm{M}_{{\odot }}}\right)^{2/3} \left(\frac{H/r_{\rm m}}{0.3}\right)^{-4/3}. \end{eqnarray} (39) This ratio is much smaller than the wind efficiency εw ∼ 0.1 of a super-Eddington accretion disc, so that such matter is easily blown away by the disc wind. As long as a possible dead disc is ionized by the ISM accretion (that occurs unless we consider low n and high v), we have encountered an unlikely set-up. Note that a fraction of the matter Mm should accrete on to the BHs before the merger, otherwise a fine-tuning is needed because the time tb is much larger than the event duration tacc. The accretion is super-Eddington, even if only a fraction of the matter accretes, and should be accompanied by a strong disc wind as suggested by numerical simulations (Ohsuga et al. 2005; Jiang, Stone & Davis 2014; Sa̧dowski et al. 2014). Therefore, it is difficult to keep the matter near the BH before the merger and the BH mergers unlikely accompany observable prompt electromagnetic counterparts. 7 SUMMARY AND DISCUSSIONS We suggest possible connections between the BH mergers observed by GWs and the high-energy sources of TeV–PeV particles in our Galaxy. The GW observations give a lower limit on the number of merged and hence highly spinning BHs as |${\sim } 70000 (\mathscr {R}_{\rm GW}/70\,{\rm Gpc}^{-3}\,{\rm yr}^{-1})$|⁠, and the spinning BHs produce relativistic jets by accreting matter and magnetic fields from the ISM. We calculate the luminosity function, the total power and the maximum acceleration energy of the BH jets, and find that the BH jets are eligible for PeVatrons, sources of CR positrons and TeV unIDs. The BH jets form extended nebulae-like PWNe. If they are observed as TeV unIDs, additional ∼300 nebulae will be discovered by CTA. We quantify the uncertainties of the estimate for the total power of the BH jets within a factor of 10±3, which is much better than before the GW detections, by considering the initial BH spin, the velocity distribution depending on the formation scenario, the accretion profile changed by the wind and the feedback by the outflow (Table 2). The uncertainties will be reduced by the GW observations, in particular, of the BH spins. It is also important to clarify the feedback by the wind from the sub-Eddington accretion disc on the Bondi–Hoyle accretion. Our considerations on the BH accretion and jet imply that the electromagnetic counterparts to BH mergers including the Fermi GBM event after GW 150914 are difficult to detect with the current sensitivity. The accretion from the ISM can evaporate the cold neutral dead disc around the BH. A slight accretion before the merger can also blow away the surrounding matter if any. These should be considered as constraints on dead disc models for prompt electromagnetic counterparts of the BH–BH merger. Although we do not go into detail in this paper, there are several sites of particle acceleration for a BH jet. First, the BH magnetosphere acts as a particle accelerator like pulsars if a gap arises with an electric field along the magnetic field (Hirotani & Pu 2016; Hirotani et al. 2016). The gamma-ray emission associated with leptonic acceleration may be detectable for nearby sources although its luminosity is usually much smaller than the BZ luminosity. Secondly, the internal shocks in the jet are possible like GRBs and AGNs. As long as B ∝ Γ/r during the propagation, the maximum acceleration energy is the same as equation (27). Thirdly, the jet dissipates the magnetic energy when the MHD approximation breaks down. This happens when the plasma density drops below the Goldreich–Julian density (Goldreich & Julian 1969), which is the minimum density required for shielding the electric field. The comoving plasma density is given by |$n^{\prime }_{\rm p} \sim L/4\pi r^2 m_{\rm u} c^3 \Gamma ^2 (1+\sigma )$|⁠, where Lσ/(1 + σ) is the BZ luminosity in equation (14), σ is the ratio of the Poynting to particle energy flux, Γ is the Lorentz factor of the jet, and we should make an appropriate correction if jets are leptonic. The comoving Goldreich–Julian density beyond the light cylinder rℓ = c/ΩH = 2rH/a* is |$n^{\prime }_{\rm GJ} \sim (\Omega _{\rm H}/2\pi q c) (r_{\rm H}/r_{\ell })^3 (r_{\ell }/r \Gamma )$|⁠. By equating |$n^{\prime }_{\rm p}$| with |$n^{\prime }_{\rm GJ}$|⁠, we obtain the radius at which the MHD breaks down, \begin{eqnarray} r_{\rm MHD} &\sim & \sqrt{\frac{\pi \kappa L}{\sigma (1+\sigma ) c}} \frac{q}{m_{\rm u} c^2} \frac{r_{\rm H}}{a_*^2 \Gamma } \sim 2 \times 10^{13}\,{\rm cm} \left[\sigma (1+\sigma )\right]^{-1/2} \nonumber \\ && \times \, a_{*}^{-2} \Gamma ^{-1} \left(\frac{L}{10^{35}\,{\rm erg}\,{\rm s}}\right)^{1/2} \left(\frac{M}{10 \,\mathrm{M}_{{\odot }}}\right). \end{eqnarray} (40) Fourthly, the termination (reverse) shock of the jet at the radius in equation (29) is also a plausible site like a hotspot of AGNs and a pulsar wind nebula for pulsars. The jet could be subject to instability, injecting energy into a cocoon/lobe before reaching the termination shock. The shocks between the cocoon and the ISM are also possible sites of particle acceleration. Note that the BH Cygnus X-1 is surrounded by a ring-like structure in radio, which may be formed by the interaction between a jet/cocoon and the ISM (Gallo et al. 2005). We do not discuss the disc emission in detail. Nearby BH discs with bremsstrahlung, synchrotron and inverse Compton emission could be detected in the future surveys (Matsumoto et al. 2017; Matsumoto & Ioka, in preparation). The accretion discs could also accelerate non-thermal particles and contribute to the observed CRs (Teraki et al., in preparation). An on-axis BH jet may be also observable if the beaming factor is larger than ∼0.01. These are interesting future problems. Acknowledgements We thank Takashi Nakamura, Tsvi Piran and Masaru Shibata for helpful comments. This work is supported by KAKENHI 24103006, 24000004, 26247042, 26287051, 17H01126, 17H06131, 17H06362, 17H06357 (KI), 17K14248 (KK), Grant-in-Aid for JSPS Research Fellow 17J09895 (TM) and by NSF Grant No. PHY-1620777 (KM). The authors also thank the Yukawa Institute for Theoretical Physics at Kyoto University. Discussions during the Long-term and Nishinomiya–Yukawa memorial workshop YITP-T-16-02 on ‘Nuclear Physics, Compact Stars, and Compact Star Mergers 2016’ were useful to complete this work. 1 " O1 was officially 2015 September 18 to 2016 January 12 before the detection of GW 150914. 2 " As we will discuss, only a small fraction of the BH spin energy is released. However, the required energy for TeV–PeV sources is only ∼10−4–10−5 of the SN energy. Thus the problem is what fraction of the spin energy is converted into high-energy particles. 3 " A factor 1/4 comes from the average over the accretion cylinder, |$\int _0^{r_{\rm B}} r {\rm d}r \int _0^{2\pi } {\rm d}\theta (r \cos \theta )^2 / r_{\rm B}^2 \int _0^{r_{\rm B}} r {\rm d}r \int _0^{2\pi } {\rm d}\theta = \frac{1}{4}$|⁠. If a turbulent velocity dominates v, the factor has a fluctuation depending on the velocity direction. 4 " We have confirmed that the results are almost similar even if we use |$L_{{\rm j}} \approx a_*^2 {\dot{M}} c^2$|⁠. 5 " LIGO imposes m1 + m2 < 100 M⊙ instead of Mmax  = 50 M⊙. 6 " For the cases without n2 in Table 1, we use a δ function. 7 " The intrinsic spectral index of CRs s is unobservable and different from the observed one by the diffusion coefficient index, which is usually obtained from observations of the boron-to-carbon ratio (Evoli, Gaggero & Grasso 2015; Genolini et al. 2015; Oliva 2016). 8 " There could be possible hotspots (Abbasi et al. 2014; Aab et al. 2015). 9 " The number would be NBH ∼ 1010 if dark matter were composed of primordial BHs that are relevant to GWs, while it is suggested that the fraction of primordial BHs in dark matter are small ∼10−4 by observations of GWs and cosmic microwave background spectral distortion (Sasaki et al. 2016). 10 " Note that the accretion does not stop even on to a binary because the disc is thick. REFERENCES Aab A. et al. , 2015 , ApJ , 804 , 15 Crossref Search ADS Abbasi R. U. et al. , 2014 , ApJ , 790 , L21 Crossref Search ADS Abbasi R. U. et al. , 2017 , Astropart. Phys. , 86 , 21 Crossref Search ADS Abbott B. P. et al. , 2016a , Phys. Rev. Lett. , 116 , 061102 Crossref Search ADS Abbott B. P. et al. , 2016b , Phys. Rev. Lett. , 116 , 241103 Crossref Search ADS Abbott B. P. et al. , 2016c , ApJ , 818 , L22 Crossref Search ADS Abbott B. P. et al. , 2016d , ApJ , 826 , L13 Crossref Search ADS Abe K. et al. , 2016 , ApJ , 830 , L11 Crossref Search ADS Abramowicz M. A. , Chen X., Kato S., Lasota J.-P., Regev O., 1995 , ApJ , 438 , L37 Crossref Search ADS Abreu P. et al. , 2010 , Astropart. Phys. , 34 , 314 Crossref Search ADS Acharya B. S. et al. , 2013 , Astropart. Phys. , 43 , 3 Crossref Search ADS Ackermann M. et al. , 2016 , ApJ , 823 , L2 Crossref Search ADS Adriani O. et al. , 2009 , Nature , 458 , 607 Crossref Search ADS PubMed Adriani O. et al. , 2016 , ApJ , 829 , L20 Crossref Search ADS Adrián-Martínez S. et al. , 2016 , Phys. Rev. D , 93 , 122010 Crossref Search ADS Agol E. , Kamionkowski M., 2002 , MNRAS , 334 , 553 Crossref Search ADS Aguilar M. et al. , 2013 , Phys. Rev. Lett. , 110 , 141102 Crossref Search ADS PubMed Aharonian F. A. , 2013 , Astropart. Phys. , 43 , 71 Crossref Search ADS Aharonian F. et al. , 2005 , Science , 307 , 1938 Crossref Search ADS PubMed Aharonian F. et al. , 2006 , ApJ , 636 , 777 Crossref Search ADS Aharonian F. et al. , 2008 , A&A , 477 , 353 Crossref Search ADS Allen S. W. , Dunn R. J. H., Fabian A. C., Taylor G. B., Reynolds C. S., 2006 , MNRAS , 372 , 21 Crossref Search ADS Armitage P. J. , Natarajan P., 1999 , ApJ , 523 , L7 Crossref Search ADS Armstrong J. W. , Rickett B. J., Spangler S. R., 1995 , ApJ , 443 , 209 Crossref Search ADS Arons J. , Lea S. M., 1976 , ApJ , 207 , 914 Crossref Search ADS Baker J. G. , Boggs W. D., Centrella J., Kelly B. J., McWilliams S. T., Miller M. C., van Meter J. R., 2007 , ApJ , 668 , 1140 Crossref Search ADS Bamba A. et al. , 2007 , PASJ , 59 , 209 Bamba A. , Yamazaki R., Kohri K., Matsumoto H., Wagner S., Pühlhofer G., Kosack K., 2009 , ApJ , 691 , 1854 Crossref Search ADS Bardeen J. M. , Petterson J. A., 1975 , ApJ , 195 , L65 Crossref Search ADS Barkov M. V. , Khangulyan D. V., Popov S. B., 2012 , MNRAS , 427 , 589 Crossref Search ADS Beckwith K. , Hawley J. F., Krolik J. H., 2009 , ApJ , 707 , 428 Crossref Search ADS Begelman M. C. , 2012 , MNRAS , 420 , 2912 Crossref Search ADS Begelman M. C. , Armitage P. J., 2014 , ApJ , 782 , L18 Crossref Search ADS Begelman M. C. , Cioffi D. F., 1989 , ApJ , 345 , L21 Crossref Search ADS Bekenstein J. D. , 1973 , ApJ , 183 , 657 Crossref Search ADS Belczynski K. , Holz D. E., Bulik T., O’Shaughnessy R., 2016 , Nature , 534 , 512 Crossref Search ADS PubMed Berkhuijsen E. M. , 1999 , in Ostrowski M., Schlickeiser R., eds, Plasma Turbulence and Energetic Particles in Astrophysics . Obserwatorium Astronomiczne, Uniwersytet Jagielloński, Kraków , Poland , p. 61 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Binney J. , Merrifield M., 1998 , Galactic Astronomy . Princeton Univ. Press , Princeton, NJ Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Bird S. , Cholis I., Muñoz J. B., Ali-Haïmoud Y., Kamionkowski M., Kovetz E. D., Raccanelli A., Riess A. G., 2016 , Phys. Rev. Lett. , 116 , 201301 Crossref Search ADS PubMed Bisnovatyi-Kogan G. S. , Ruzmaikin A. A., 1976 , Ap&SS , 42 , 401 Crossref Search ADS Blandford R. D. , 2000 , Phys. Scr., T 85 , 191 Crossref Search ADS Blandford R. D. , Begelman M. C., 1999 , MNRAS , 303 , L1 Crossref Search ADS Blandford R. D. , Begelman M. C., 2004 , MNRAS , 349 , 68 Crossref Search ADS Blandford R. D. , Znajek R. L., 1977 , MNRAS , 179 , 433 Crossref Search ADS Bland-Hawthorn J. , Reynolds R., 2000 , Encyclopedia of Astronomy and Astrophysics, article 2636 . IoP Publishing , Bristol Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Blümer J. , Engel R., Hörandel J. R., 2009 , Progress Part. Nucl. Phys. , 63 , 293 Crossref Search ADS Bondi H. , 1952 , MNRAS , 112 , 195 Crossref Search ADS Bondi H. , Hoyle F., 1944 , MNRAS , 104 , 273 Crossref Search ADS Bonnor W. B. , Rotenberg M. A., 1961 , Proc. R. Soc. Lond. Ser. A , 265 , 109 Crossref Search ADS Brandt W. N. , Podsiadlowski P., Sigurdsson S., 1995 , MNRAS , 277 , L35 Campana S. , Pardi M. C., 1993 , A&A , 277 , 477 Campanelli M. , Lousto C. O., Zlochower Y., Merritt D., 2007 , ApJ , 659 , L5 Crossref Search ADS Cao X. , 2011 , ApJ , 737 , 94 Crossref Search ADS Cardoso V. , Macedo C. F. B., Pani P., Ferrari V., 2016 , J. Cosmol. Astropart. Phys. , 5 , 054 Crossref Search ADS Carr B. J. , 1979 , MNRAS , 189 , 123 Crossref Search ADS Chevalier C. , Ilovaisky S. A., 1998 , A&A , 330 , 201 Chisholm J. R. , Dodelson S., Kolb E. W., 2003 , ApJ , 596 , 437 Crossref Search ADS Cioffi D. F. , McKee C. F., Bertschinger E., 1988 , ApJ , 334 , 252 Crossref Search ADS Clemens D. P. , Sanders D. B., Scoville N. Z., 1988 , ApJ , 327 , 139 Crossref Search ADS Connaughton V. et al. , 2016 , ApJ , 826 , L6 Crossref Search ADS Cordes J. M. , Chernoff D. F., 1998 , ApJ , 505 , 315 Crossref Search ADS de Jager O. C. et al. , 2009 , preprint (arXiv:0906.2644) Dominik M. et al. , 2015 , ApJ , 806 , 263 Crossref Search ADS Draine B. T. , 2011 , Physics of the Interstellar and Intergalactic Medium . Princeton Univ. Press , Princeton, NJ Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Esin A. A. , McClintock J. E., Narayan R., 1997 , ApJ , 489 , 865 Crossref Search ADS Evans P. A. et al. , 2016a , MNRAS , 460 , L40 Evans P. A. et al. , 2016b , MNRAS , 462 , 1591 Crossref Search ADS Evoli C. , Gaggero D., Grasso D., 2015 , J. Cosmol. Astropart. Phys. , 12 , 039 Crossref Search ADS Fender R. P. , Belloni T. M., Gallo E., 2004 , MNRAS , 355 , 1105 Crossref Search ADS Fender R. P. , Maccarone T. J., Heywood I., 2013 , MNRAS , 430 , 1538 Crossref Search ADS Ferrière K. M. , 2001 , Rev. Modern Phys. , 73 , 1031 Crossref Search ADS Fitchett M. J. , 1983 , MNRAS , 203 , 1049 Crossref Search ADS Fujinaga T. et al. , 2011 , PASJ , 63 , S857 Fujita Y. , Inoue S., Nakamura T., Manmoto T., Nakamura K. E., 1998 , ApJ , 495 , L85 Crossref Search ADS Fujita Y. , Kohri K., Yamazaki R., Ioka K., 2009 , Phys. Rev. D , 80 , 063003 Crossref Search ADS Gaisser T. K. , Stanev T., Tilav S., 2013 , Frontiers Phys. , 8 , 748 Crossref Search ADS Gallo E. , Fender R. P., Pooley G. G., 2003 , MNRAS , 344 , 60 Crossref Search ADS Gallo E. , Fender R., Kaiser C., Russell D., Morganti R., Oosterloo T., Heinz S., 2005 , Nature , 436 , 819 Crossref Search ADS PubMed Genolini Y. , Putze A., Salati P., Serpico P. D., 2015 , A&A , 580 , A9 Crossref Search ADS Goldreich P. , Julian W. H., 1969 , ApJ , 157 , 869 Crossref Search ADS Gonzalez J. A. , Sperhake U., Bruegmann B., Hannam M., Husa S., 2007 , Phys. Rev. Lett. , 98 , 091101 Crossref Search ADS PubMed Greiner J. , Burgess J. M., Savchenko V., Yu H.-F., 2016 , ApJ , 827 , L38 Crossref Search ADS Grindlay J. E. , 1978 , ApJ , 221 , 234 Crossref Search ADS Han J. L. , Ferriere K., Manchester R. N., 2004 , ApJ , 610 , 820 Crossref Search ADS Heckler A. F. , Kolb E. W., 1996 , ApJ , 472 , L85 Crossref Search ADS Heger A. , Fryer C. L., Woosley S. E., Langer N., Hartmann D. H., 2003 , ApJ , 591 , 288 Crossref Search ADS HESS Collaboration et al. , 2016 , Nature , 531 , 476 Crossref Search ADS PubMed Heyer M. , Dame T. M., 2015 , ARA&A , 53 , 583 Crossref Search ADS Hillas A. M. , 2005 , J. Phys. G: Nucl. Phys. , 31 , R95 Crossref Search ADS Hirotani K. , Pu H.-Y., 2016 , ApJ , 818 , 50 Crossref Search ADS Hirotani K. , Pu H.-Y., Chun-Che Lin L., Chang H.-K., Inoue M., Kong A. K. H., Matsushita S., Tam P.-H. T., 2016 , ApJ , 833 , 142 Crossref Search ADS Hirschi R. , Meynet G., Maeder A., 2005 , A&A , 443 , 581 Crossref Search ADS Hollenbach D. , Johnstone D., Lizano S., Shu F., 1994 , ApJ , 428 , 654 Crossref Search ADS Hoyle F. , Lyttleton R. A., 1939 , Proc. Cambridge Philos. Soc. , 35 , 405 Crossref Search ADS Ichimaru S. , 1977 , ApJ , 214 , 840 Crossref Search ADS Inutsuka S.-i. , Inoue T., Iwasaki K., Hosokawa T., 2015 , A&A , 580 , A49 Crossref Search ADS Ioka K. , 2010 , Progress Theor. Phys. , 123 , 743 Crossref Search ADS Ioka K. , Mészáros P., 2010 , ApJ , 709 , 1337 Crossref Search ADS Ioka K. , Chiba T., Tanaka T., Nakamura T., 1998 , Phys. Rev. D , 58 , 063003 Crossref Search ADS Ioka K. , Tanaka T., Nakamura T., 1999 , Phys. Rev. D , 60 , 083512 Crossref Search ADS Ipser J. R. , Price R. H., 1977 , ApJ , 216 , 578 Crossref Search ADS Ipser J. R. , Price R. H., 1982 , ApJ , 255 , 654 Crossref Search ADS Jiang Y.-F. , Stone J. M., Davis S. W., 2014 , ApJ , 796 , 106 Crossref Search ADS KamLAND Collaboration et al. , 2016 , ApJ , 829 , L34 Crossref Search ADS Kashiyama K. , Ioka K., Kawanaka N., 2011 , Phys. Rev. D , 83 , 023002 Crossref Search ADS Kasliwal M. M. et al. , 2016 , ApJ , 824 , L24 Crossref Search ADS Kato S. , Fukue J., Mineshige S., 2008 , Black-Hole Accretion Disks — Towards a New Paradigm — . Kyoto Univ. Press , Kyoto, Japan Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Kauffmann J. , Bertoldi F., Bourke T. L., Evans N. J., II, Lee C. W., 2008 , A&A , 487 , 993 Crossref Search ADS Kawanaka N. , Ioka K., Nojiri M. M., 2010 , ApJ , 710 , 958 Crossref Search ADS Kawanaka N. , Ioka K., Ohira Y., Kashiyama K., 2011 , ApJ , 729 , 93 Crossref Search ADS Kimura S. S. , Takahashi S. Z., Toma K., 2017 , MNRAS , 465 , 4406 Crossref Search ADS Kinugawa T. , Inayoshi K., Hotokezaka K., Nakauchi D., Nakamura T., 2014 , MNRAS , 442 , 2963 Crossref Search ADS Kinugawa T. , Nakano H., Nakamura T., 2016a , Progress Theor. Exp. Phys. , 2016 , 103E01 Crossref Search ADS Kinugawa T. , Miyamoto A., Kanda N., Nakamura T., 2016b , MNRAS , 456 , 1093 Crossref Search ADS Kinugawa T. , Nakano H., Nakamura T., 2016c , Progress Theor. Exp. Phys. , 2016 , 031E01 Crossref Search ADS Kobayashi T. , Komori Y., Yoshida K., Nishimura J., 2004 , ApJ , 601 , 340 Crossref Search ADS Kobulnicky H. A. et al. , 2014 , ApJS , 213 , 34 Crossref Search ADS Kohri K. , Ioka K., Fujita Y., Yamazaki R., 2016 , Progress Theor. Exp. Phys. , 2016 , 021E01 Crossref Search ADS Koide S. , Shibata K., Kudoh T., Meier D. L., 2002 , Science , 295 , 1688 Crossref Search ADS PubMed Kotera K. , Silk J., 2016 , ApJ , 823 , L29 Crossref Search ADS Kuijken K. , Gilmore G., 1989a , MNRAS , 239 , 571 Crossref Search ADS Kuijken K. , Gilmore G., 1989b , MNRAS , 239 , 605 Crossref Search ADS Kulkarni S. R. , Hut P., McMillan S., 1993 , Nature , 364 , 421 Crossref Search ADS Kyutoku K. , Seto N., 2016 , MNRAS , 462 , 2177 Crossref Search ADS Lamberts A. , Garrison-Kimmel S., Clausen D., Hopkins P., 2016 , MNRAS , 463 , L31 Li X. , Zhang F.-W., Yuan Q., Jin Z.-P., Fan Y.-Z., Liu S.-M., Wei D.-M., 2016 , ApJ , 827 , L16 Crossref Search ADS Lipunov V. M. , Kornilov V., Gorbovskoy E., Tiurina N., Balanutsa P., Kuznetsov A., 2017 , New Astron. , 51 , 122 Crossref Search ADS Loeb A. , 2016 , ApJ , 819 , L21 Crossref Search ADS Longair M. S. , 2011 , High Energy Astrophysics . Cambridge Univ. Press , Cambridge Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Lubow S. H. , Papaloizou J. C. B., Pringle J. E., 1994 , MNRAS , 267 , 235 Crossref Search ADS Lynden-Bell D. , 1969 , Nature , 223 , 690 Crossref Search ADS McDowell J. , 1985 , MNRAS , 217 , 77 Crossref Search ADS McKinney J. C. , Tchekhovskoy A., Blandford R. D., 2012 , MNRAS , 423 , 3083 Crossref Search ADS McKinney J. C. , Tchekhovskoy A., Blandford R. D., 2013 , Science , 339 , 49 Crossref Search ADS PubMed Mandel I. , de Mink S. E., 2016 , MNRAS , 458 , 2634 Crossref Search ADS Mapelli M. , 2016 , MNRAS , 459 , 3432 Crossref Search ADS Marchant P. , Langer N., Podsiadlowski P., Tauris T. M., Moriya T. J., 2016 , A&A , 588 , A50 Crossref Search ADS Matsumoto H. et al. , 2007 , PASJ , 59 , 199 Matsumoto H. et al. , 2008 , PASJ , 60 , S163 Matsumoto T. , Teraki Y., Ioka K., 2017 , preprint (arXiv:1704.05047) Mészáros P. , 1975 , A&A , 44 , 59 Michel F. C. , 1972 , Ap&SS , 15 , 153 Crossref Search ADS Miller-Jones J. C. A. , 2014 , Publ. Astron. Soc. Aust. , 31 , e016 Crossref Search ADS Milosavljević M. , Bromm V., Couch S. M., Oh S. P., 2009 , ApJ , 698 , 766 Crossref Search ADS Mirabel I. F. , Rodrigues I., 2003 , Science , 300 , 1119 Crossref Search ADS PubMed Mirabel I. F. , Dhawan V., Mignani R. P., Rodrigues I., Guglielmetti F., 2001 , Nature , 413 , 139 Crossref Search ADS PubMed Mirabel I. F. , Mignani R., Rodrigues I., Combi J. A., Rodríguez L. F., Guglielmetti F., 2002 , A&A , 395 , 595 Crossref Search ADS Mizuta A. , Ioka K., 2013 , ApJ , 777 , 162 Crossref Search ADS Morokuma T. et al. , 2016 , PASJ , 68 , L9 Morsony B. J. , Workman J. C., Ryan D. M., 2016 , ApJ , 825 , L24 Crossref Search ADS Murase K. , Kashiyama K., Mészáros P., Shoemaker I., Senno N., 2016 , ApJ , 822 , L9 Crossref Search ADS Nakamura T. , Sasaki M., Tanaka T., Thorne K. S., 1997 , ApJ , 487 , L139 Crossref Search ADS Nakamura T. , Nakano H., Tanaka T., 2016a , Phys. Rev. D , 93 , 044048 Crossref Search ADS Nakamura T. et al. , 2016b , Progress Theor. Exp. Phys. , 2016 , 093E01 Crossref Search ADS Nakano T. , Nishi R., Umebayashi T., 2002 , ApJ , 573 , 199 Crossref Search ADS Narayan R. , Yi I., 1994 , ApJ , 428 , L13 Crossref Search ADS Narayan R. , Yi I., 1995 , ApJ , 452 , 710 Crossref Search ADS Narayan R. , Igumenshchev I. V., Abramowicz M. A., 2003 , PASJ , 55 , L69 Norman C. A. , Melrose D. B., Achterberg A., 1995 , ApJ , 454 , 60 Crossref Search ADS Novikov I. D. , Thorne K. S., 1973 , in Dewitt C., Dewitt B. S., eds, Black Holes (Les Astres Occlus) . Gordon & Breach , New York , p. 343 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Ohsuga K. , Mineshige S., 2011 , ApJ , 736 , 2 Crossref Search ADS Ohsuga K. , Mori M., Nakamoto T., Mineshige S., 2005 , ApJ , 628 , 368 Crossref Search ADS Oliva A. , 2016 , in Proceeding of Science, 34th International Cosmic Ray Conference, PoS(ICRC2015)265 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Palliyaguru N. T. et al. , 2016 , ApJ , 829 , L28 Crossref Search ADS Peres A. , 1962 , Phys. Rev. , 128 , 2471 Crossref Search ADS Perna R. , Lazzati D., Giacomazzo B., 2016 , ApJ , 821 , L18 Crossref Search ADS Pfahl E. , Rappaport S., Podsiadlowski P., Spruit H., 2002 , ApJ , 574 , 364 Crossref Search ADS Portegies Zwart S. F. , McMillan S. L. W., 2000 , ApJ , 528 , L17 Crossref Search ADS Pringle J. E. , Rees M. J., 1972 , A&A , 21 , 1 Racusin J. L. et al. , 2017 , ApJ , 835 , 82 Crossref Search ADS Rafferty D. A. , McNamara B. R., Nulsen P. E. J., Wise M. W., 2006 , ApJ , 652 , 216 Crossref Search ADS Ramirez-Ruiz E. , Celotti A., Rees M. J., 2002 , MNRAS , 337 , 1349 Crossref Search ADS Reid M. J. , McClintock J. E., Narayan R., Gou L., Remillard R. A., Orosz J. A., 2011 , ApJ , 742 , 83 Crossref Search ADS Remillard R. A. , McClintock J. E., 2006 , ARA&A , 44 , 49 Crossref Search ADS Repetto S. , Davies M. B., Sigurdsson S., 2012 , MNRAS , 425 , 2799 Crossref Search ADS Rodriguez C. L. , Morscher M., Pattabiraman B., Chatterjee S., Haster C.-J., Rasio F. A., 2015 , Phys. Rev. Lett. , 115 , 051101 Crossref Search ADS PubMed Rodriguez C. L. , Chatterjee S., Rasio F. A., 2016 , Phys. Rev. D , 93 , 084029 Crossref Search ADS Russell H. R. , McNamara B. R., Edge A. C., Hogan M. T., Main R. A., Vantyghem A. N., 2013 , MNRAS , 432 , 530 Crossref Search ADS Sa̧dowski A. , Narayan R., McKinney J. C., Tchekhovskoy A., 2014 , MNRAS , 439 , 503 Crossref Search ADS Sa̧dowski A. , Lasota J.-P., Abramowicz M. A., Narayan R., 2016 , MNRAS , 456 , 3915 Crossref Search ADS Sakai M. , Yajima Y., Matsumoto H., 2011 , PASJ , 63 , S879 Salpeter E. E. , 1964 , ApJ , 140 , 796 Crossref Search ADS Sana H. et al. , 2012 , Science , 337 , 444 Crossref Search ADS PubMed Sasaki M. , Suyama T., Tanaka T., Yokoyama S., 2016 , Phys. Rev. Lett. , 117 , 061101 Crossref Search ADS PubMed Savchenko V. et al. , 2016 , ApJ , 820 , L36 Crossref Search ADS Scoville N. Z. , Sanders D. B., 1987 , in Hollenbach D. J., Thronson H. A. Jr, eds, Astrophysics and Space Science Library Vol. 134, Interstellar Processes . Reidel , Dordrecht , p. 21 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Serpico P. D. , 2012 , Astropart. Phys. , 39 , 2 Crossref Search ADS Sesana A. , 2016 , Phys. Rev. Lett. , 116 , 231102 Crossref Search ADS PubMed Shakura N. I. , Sunyaev R. A., 1973 , A&A , 24 , 337 Shapiro S. L. , 1973 , ApJ , 180 , 531 Crossref Search ADS Shapiro S. L. , Teukolsky S. A., 1983 , Black Holes, White Dwarfs, and Neutron Stars: The Physics of Compact Objects . Wiley-Interscience , New York Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Shvartsman V. F. , 1971 , SvA , 15 , 377 Sigurdsson S. , Hernquist L., 1993 , Nature , 364 , 423 Crossref Search ADS Sikora M. , Begelman M. C., 2013 , ApJ , 764 , L24 Crossref Search ADS Tavani M. et al. , 2016 , ApJ , 825 , L4 Crossref Search ADS Tchekhovskoy A. , Narayan R., McKinney J. C., 2011 , MNRAS , 418 , L79 Crossref Search ADS The LIGO Scientific Collaboration et al. , 2016 , Phys. Rev. X , 6 , 041015 The Pierre Auger Collaboration et al. , 2016 , Phys. Rev. D , 94 , 122007 Crossref Search ADS Troja E. , Read A. M., Tiengo A., Salvaterra R., 2016 , ApJ , 822 , L8 Crossref Search ADS Tutukov A. V. , Yungelson L. R., 1993 , MNRAS , 260 , 675 Crossref Search ADS Veres P. , Preece R. D., Goldstein A., Mészáros P., Burns E., Connaughton V., 2016 , ApJ , 827 , L34 Crossref Search ADS Waxman E. , 2004 , Pramana , 62 , 483 Crossref Search ADS White N. E. , van Paradijs J., 1996 , ApJ , 473 , L25 Crossref Search ADS Yamazaki R. , Kohri K., Bamba A., Yoshida T., Tsuribe T., Takahara F., 2006 , MNRAS , 371 , 1975 Crossref Search ADS Yamazaki R. , Asano K., Ohira Y., 2016 , Progress Theor. Exp. Phys. , 2016 , 051E01 Crossref Search ADS Yoon S.-C. , Langer N., 2005 , A&A , 443 , 643 Crossref Search ADS Yuan F. , Narayan R., 2014 , ARA&A , 52 , 529 Crossref Search ADS Yuan F. , Gan Z., Narayan R., Sadowski A., Bu D., Bai X.-N., 2015 , ApJ , 804 , 101 Crossref Search ADS Zamaninasab M. , Clausen-Brown E., Savolainen T., Tchekhovskoy A., 2014 , Nature , 510 , 126 Crossref Search ADS PubMed Zel’dovich Y. B. , 1964 , Soviet Phys. Doklady , 9 , 195 Zhang B. , 2016 , ApJ , 827 , L31 Crossref Search ADS Zlochower Y. , Lousto C. O., 2015 , Phys. Rev. D , 92 , 024022 Crossref Search ADS APPENDIX A: COMPARISON WITH PREVIOUS WORKS Studying the accreting BHs in our Galaxy has a long history from 1920s Eddington era. However, the observations of GWs set a lower limit on the number of spinning BHs for the first time. Our paper gives the first considerations on the Galactic BHs after the discovery of GWs. Hoyle & Lyttleton (1939) considered the effect of the ISM accretion on the Sun’s radiation for explaining changes in terrestrial climate. Bondi & Hoyle (1944) investigated the accretion in detail including the effect of perturbations. Bondi (1952) included the pressure effects to complete the Bondi–Hoyle formula. Zel’dovich (1964) and Salpeter (1964) suggested accretion on to a BH as an important source of radiation. Shvartsman (1971) treated both the fluid dynamics and radiative processes by employing non-relativistic approximations. Michel (1972) considered the general relativistic version of the Bondi–Hoyle accretion. Shapiro (1973) made a general relativistic treatment of both the fluid mechanics and radiative processes. Mészáros (1975) included the heating due to dissipation of turbulence and magnetic reconnection, and the radiation from CRs captured with the accretion. The formation of a disc around a massive BH is considered by Lynden-Bell (1969). Pringle & Rees (1972) developed the accretion disc model. Shakura & Sunyaev (1973) completed the standard disc model. Novikov & Thorne (1973) made a general relativistic version of the standard disc model. Ipser & Price (1977, 1982) carried out detailed calculations of radiation from Galactic massive BHs that are spherically accreting from the ISM. Grindlay (1978) discussed X-ray limits on Galactic BHs that are spherically accreting from the ISM. Carr (1979) placed a limit on BHs with the background radiation density. McDowell (1985) suggested that molecular clouds offer high density for spherically accreting BHs to produce detectable fluxes. Campana & Pardi (1993) estimated the number of BHs in the molecular clouds. Heckler & Kolb (1996) proposed a search strategy for isolated stellar mass BHs in the solar neighbourhood using optical surveys. Chisholm, Dodelson & Kolb (2003) gave possible BH candidates using optical and X-ray surveys. Modern estimates for radiation from isolated stellar mass BHs began after Fujita et al. (1998) adopted the ADAF model and Armitage & Natarajan (1999) considered a jet from an isolated accreting BH, although they do not estimate the luminosity function, the total power and the acceleration energy. They did not also consider the MAD state. Agol & Kamionkowski (2002) calculated the luminosity function of isolated accreting BHs, although they assumed a constant efficiency of the disc emission and do not consider a BH jet. Barkov et al. (2012) also considered a jet from an isolated accreting BH, although they do not consider the luminosity function and hence their estimate on the total power is not correct. They did not also consider the MAD state. Fender, Maccarone & Heywood (2013) calculated the X-ray luminosity distribution and suggested a discrepancy between the theoretical expectation and the hard X-ray surveys, although their prescription for radiatively inefficient accretion is very simple. Nakamura, Nakano & Tanaka (2016a) considered the Bondi–Hoyle accretion for the optical counterparts of nearby BH mergers. © 2017 The Authors Published by Oxford University Press on behalf of the Royal Astronomical Society
TI - GW 150914-like black holes as Galactic high-energy sources
JF - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/stx1337
DA - 2017-09-21
UR - https://www.deepdyve.com/lp/oxford-university-press/gw-150914-like-black-holes-as-galactic-high-energy-sources-2NmopZC0YP
SP - 3332
VL - 470
IS - 3
DP - DeepDyve
ER -