Where Are the Electromagnetic-wave Counterparts of Stellar-mass Binary Black Hole Mergers?

Since the first gravitational-wave (GW) event detected in 2015, there have been 11 published events (LIGO/Virgo science runs O1–O2; Abbott et al. 2019), and this number is increasing in the current O3 run.³ The observation of the electromagnetic-wave (EM) counterpart of GW170817 offered promising prospects for multimessenger in fundamental physics (e.g., Baker et al. 2017; Nishizawa & Kobayashi 2018), astrophysics (e.g., Abbott et al. 2018; Wu et al. 2019), and cosmology (Vitale & Chen 2018). Among the targets of LIGO/Virgo and other ground-based GW detectors, double neutron star (DNS) mergers, and more rarely neutron star-black hole (NBH) mergers, were the only kinds of sources expected to be accompanied by EM counterparts. This means that the majority of targets, i.e., mergers of stellar-mass binaries black holes (BBHs) were thought to have no hope of being detected via EM. As a result, follow-up searches for EM counterparts are mainly focused on a GW detection that shows a high possibility of being a DNS or a NBH. On the other hand, if some of the BBH merger can also be detected with EM, multimessenger astronomy can be performed distances one order of magnitude farther than that of DNSs. Sample sizes would also increase significantly.

There are three well-known channels for forming coalescing BBHs: (1) isolated evolution of massive binaries via the common envelope phase (Ivanova et al. 2013; Postnov & Yungelson 2014); (2) isolated evolution of massive binaries via rotational chemical mixing (de Mink & Mandel 2016; Mandel & de Mink 2016; Marchant et al. 2016); (3) dynamical capture in dense clusters (see Chatterjee et al. 2017 and references therein). The abovementioned channels produce BBH mergers in "clean" environments, where no electromagnetic radiation is expected (Liu & Zhang 2009; Zhang et al. 2016). We refer to those channels as "clean channels" (see Mandel & Farmer 2018 for a pleasant review). Recently, a new channel of producing BBH mergers has received attention, i.e., BBHs embedded in the disk surrounding a supermassive black hole (SMBH), or an active galactic nucleus (AGN). We refer to this channel as "the AGN-disk channel" (Bartos et al. 2017; Stone et al. 2017; Leigh et al. 2018; Yi et al. 2018). In this circumstance, some authors predicted that there would be EM radiation accompanying GW. Farris et al. (2010) found that the shock due to the orbital motion of the binary could provide X-ray radiation. Such radiation may not be luminous enough to be detected from a distance larger than several hundred Mpc. Bartos et al. (2017) suggested there is thermal emission from the transient accretion disk of the BBH and/or Doppler boosted emission from the relativistic outflow. The fluxes of those high-energy (HE) emissions would be well under the limit of current detectors, unless the accretion rate is highly super-Eddington (we would also like to mention some less-standard scenarios, where gamma-ray radiations are predicted as EM counterparts of BBHs; see Veres et al. 2019). Differing from previous studies that focused on HE radiation, in this Letter we propose that these BBH mergers will produce flares in the radio bands (several hundred MHz to GHz), like those in fast radio bursts (FRBs) but with a longer duration. The spectral flux densities of such radio flares are above the detecting limit of radio telescopes with large collecting areas. In Section 2, we will present the emission mechanism and make predictions about the EM counterparts.

It is beneficial to ask the following question: how can we distinguish between the AGN-disk channel BBH-originating GW events and ones from clean channels, solely with GW observation? Here we suggest that their total masses and effective spin parameters (M_tot and χ_eff) can serve as indicators of their origins. The simulation of Postnov & Kuranov (2019) showed that BBH mergers from isolated massive binary evolution have χ_eff clustering around zero. For the dynamical channel, a similar distribution is expected. Meanwhile, the AGN-disk channel BBH has the opportunity to accrete mass and angular momenta from surrounding materials. As a result, we expect them to possess larger M_tot and χ_eff. Indeed, GW170729 has the largest M_tot, whose χ_eff also stands out of the cluster around zero, where others lie. Hence, we suspect that GW170729 could be the first example of AGN-disk channel BBH mergers.

In Section 3, we give the probability density distribution of M_tot and χ_eff of BBH mergers that originated from the AGN-disk channel. An analytic formalism and Monte-Carlo simulation are used in this study. In Section 4, we discuss the observational aspects of searching for EM counterparts of GWs from BBHs and their host galaxies. We also discuss caveats of our results. We use geometric units (G = c = 1), unless otherwise indicated.

The scenario that we are studying in this Letter is the same as that in previous literature: a BBH, either formed in or migrated to and trapped in the disk around a SMBH, inspirals and coalesces with the assistance of accretion from the surrounding matter. As shown in previous numerical simulations, several features are expected (see Figure 1 for an illustration): a region with less dense gas in the AGN disk is cleaned by the BBH. It can be a gap-ring in the orbit where the BBH lies, or a cavity around the BBH (Baruteau et al. 2011). Gas inflows through overdense streams from the AGN disk, forming a circumbinary disk (CD). Some of the gas in the CD is transferred onto individual BHs via their mini-disks. Meanwhile, other gases in the CD are barricaded from being accreted due to the slingshot effect of the binary orbit, which gives rise to an inner gap in the CD.

Zoom In Zoom Out Reset image size

After the coalescence of the BBH, remnant materials in mini-disks and those accumulated in the CD will be quickly accreted onto the newly formed, highly spinning Kerr BH. The accretion can be transient and temporarily super-Eddington, analogous to that in tidal disruption events (Guillochon & Ramirez-Ruiz 2013; Kara et al. 2016; Wu et al. 2018) or microquasars (Soria et al. 2014). As proposed in Yi et al. (2019), a transient high-mass transfer onto a stellar-mass BH from its stellar companion could trigger a clumpy jet, from which FRBs could arise. A similar scenario also applies here: The transient accretion of the remnant gas onto the newly formed Kerr BH can launch an inhomogeneous Blandford–Znajek jet. Gaseous clumpy ejecta in the jet have diverse velocities and are ejected at different instances. Therefore, they might collide among each other at some distance l_col. l_col is related to the time interval between the ejected instances δt and the average Lorentz factor of the plasma bulk motion γ:

$$\begin{eqnarray}&&{l}_{\mathrm{col}}\sim c{\gamma }^{2}\delta t.\end{eqnarray} \tag{ 1 }$$

We denote the height of the transient accretion disk at the innermost radius as h. The typical separation between clumps being accreted is of the same order of magnitude of h. Therefore, δt ∼ h/v, where v is the freefalling velocity at the innermost radius.

The collision between clumps in the jet could trigger plasma oscillation, which will form temporarily charge-separated bunches. The plasma frequency of the oscillation is

$$\begin{eqnarray}&&{\nu }_{\mathrm{pls}}(l)=2\sqrt{\displaystyle \frac{\gamma {e}^{2}{n}_{{\rm{e}}}(l)}{\pi {m}_{{\rm{e}}}}},\end{eqnarray} \tag{ 2 }$$

where n_e(l) is the number density of electrons, which is simply related to the mass density of the fully ionized plasma ρ(l). Because the dimension of the clump expands along the jet with increasing cone radius, the density of the plasma decreases accordingly as

$$\begin{eqnarray}&&\rho (l)={\rho }_{0}{\left(\displaystyle \frac{h}{s(l)}\right)}^{3},\end{eqnarray} \tag{ 3 }$$

where ρ₀ is the mass density of the gas at the innermost radius, which can be related to the accretion rate with

$$\begin{eqnarray}&&\dot{M}=2\pi {r}_{\mathrm{in}}h{\rho }_{0}v.\end{eqnarray} \tag{ 4 }$$

Placing the above equations into Equation (2), we obtain the plasma frequency of the ejecta as a function of the distance to the BH:

$$\begin{eqnarray}&&{\nu }_{\mathrm{pls}}(l)=3.45\sqrt{\displaystyle \frac{{f}_{\mathrm{Edd}}}{{m}_{\bullet }\eta \tilde{h}}}{\gamma }_{100}^{-2.5}{\left(\displaystyle \frac{{l}_{\mathrm{col}}}{{{\rm{\Theta }}}_{0.1}l}\right)}^{1.5}\,\mathrm{GHz}.\end{eqnarray} \tag{ 5 }$$

In the above equation, f_Edd and η are the Eddington ratio and the radiation efficiency of accretion, respectively.

When such a bunch slides along the curling magnetic field lines in the jet, curvature radiation will be emitted. The frequency of the curvature radiation in the rest frame is

$$\begin{eqnarray}&&{\nu }_{\mathrm{cur}}(l)=2{\gamma }_{\mathrm{los}}{\gamma }_{\parallel }^{3}c/s(l),\end{eqnarray} \tag{ 6 }$$

where γ_los is the Lorentz factor of the bulk motion along the line of sight, and γ_∥ is the Lorentz factor corresponding to the sliding of the plasma along the magnetic field lines. For simplicity, we assume γ_los ∼ γ_∥ ∼ γ; s(l) is the local radius of the curvature of the magnetic field lines. Because the magnetic field lines are highly spiral in the jet, we can relate s(l) approximately with the opening angle of the cone of jet as s(l) ∼ Θl.

As a result, the curvature radiation frequency as a function of distance from the BH is

$$\begin{eqnarray}&&{\nu }_{\mathrm{cur}}(l)=3.8{\tilde{h}}^{-1}\displaystyle \frac{{\gamma }_{100}^{2}}{{m}_{\bullet }}\displaystyle \frac{{l}_{\mathrm{col}}}{{{\rm{\Theta }}}_{0.1}l}\,\mathrm{GHz},\end{eqnarray} \tag{ 7 }$$

where $\tilde{h}$ is the height of the transient accretion disk scaled with the innermost radius r_in, γ₁₀₀ is the Lorentz factor of the bulk motion of the bunch scaled with 100, m_• is the mass of the merged BH in unit of M_⊙, and Θ_0.1 is the opening angle of the jet cone divided by 0.1.

When the condition that ν_cur = ν_pls occurs at some distance where l > l_col, the plasma instability will grow, and the kinetic energy of the bunch will be coherently emitted in a μs timescale. It corresponds to a spike in the radio light curve. Each collision among the jet clumps has the opportunity to give such a spike, thus we expect a multicomponent radio flare, which is composed by many spikes, like the sub-structures seen in FRBs. The clumpy jet is thought to arise due to the sudden increase of the accretion rate at the moment of merger. The follow-up accretion is expected to be less temporally varying, and the accretion rate is as high as it is in the onset.

The duration and the power of the flare are proportional to the mass of the BH. The former corresponds to the freefalling timescale near the innermost stable orbit (ISCO):

$$\begin{eqnarray}&&{\tau }_{\mathrm{dur}}\sim {r}_{\mathrm{ISCO}}/c.\end{eqnarray} \tag{ 8 }$$

For a BH with 60 M_⊙, τ_dur ≈ 1.8 ms, the apparent luminosity is

$$\begin{eqnarray}&&L=\displaystyle \frac{8\eta {m}_{\bullet }{f}_{\mathrm{Edd}}}{2\pi (1-\cos {\rm{\Theta }})}\times {10}^{38}\,\mathrm{erg}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 9 }$$

With a typical mass for coalescent BHs and a super-Eddington accretion rate, L can be 10⁴¹–10⁴² erg s⁻¹. Therefore, we expect the EM counterparts to be longer and brighter FRBs than usual. Such FRBs, if localized with high accuracy, are expected to be found at the core regions of active galaxies, and to be simultaneous with a GW chirp. We encourage the reader to refer to Yi et al. (2019) for a more detailed derivation and discussion of the physics of this model.

The spin evolution of a BH under accretion along the equatorial plane was first investigated by Bardeen (1970). His result did not take into account the torque exerted by photons emitted from the accretion disk. This term will require a correction to Bardeen's equation, yielding the extreme value of a* = 0.998 instead of unity (Thorne 1974). We neglect this small correction here. When the angular momentum of the accretion disk is misaligned with the spin of the BH, the evolution of the BH angular momentum is (Perego et al. 2009)

$$\begin{eqnarray}&&\displaystyle \frac{d{{\boldsymbol{J}}}_{\mathrm{BH}}}{{dt}}=\dot{M}{L}_{\mathrm{ISCO}}\hat{{\boldsymbol{l}}}+4\pi {\int }_{\mathrm{disk}}\displaystyle \frac{{\boldsymbol{L}}(R)\times {{\boldsymbol{J}}}_{\mathrm{BH}}}{{R}^{2}}{dR},\end{eqnarray} \tag{ 10 }$$

where L_ISCO is the specific angular momentum brought onto the black hole from ISCO; $\hat{{\boldsymbol{l}}}$ is a unit vector, which is guaranteed to be parallel with J_BH due to the Bardeen–Petterson effect (Bardeen & Petterson 1975); and L(R) is the angular momentum per unit area in the accretion disk at distance R from the BH. The first term on the right is the change rate of the spin modulus, and the second term is the change rate of the spin direction.

Note that

$$\begin{eqnarray}&&{a}_{* }=\pm \displaystyle \frac{{J}_{\mathrm{BH}}}{{M}^{2}},\end{eqnarray} \tag{ 11 }$$

where J_BH is the modulus of the spin angular momentum. When the spin is counter-rotating with the binary orbit, the signature of a_* is minus. The evolution of a_* just follows Bardeen's equation:

$$\begin{eqnarray}&&\displaystyle \frac{{{da}}_{* }}{d\mathrm{ln}M}=\displaystyle \frac{1}{M}\displaystyle \frac{{L}_{\mathrm{ISCO}}}{{E}_{\mathrm{ISCO}}}-2{a}_{* },\end{eqnarray} \tag{ 12 }$$

where E_ISCO is the specific energy that was brought onto the black hole from ISCO. With the explicit expression of L_ISCO and E_ISCO (Bardeen et al. 1972), the above equation can be integrated to give a_* as function of x and initial a_*,0, where $x\equiv {M}_{{\rm{f}}}/{M}_{0}$ is the ratio between the mass after accretion and the initial mass of the BH. The explicit expression of a_*(a_*,0, x) is lengthy, thus we just use the numerical integration of Equation (12) for our purposes.

The second term on the right side of Equation (10) governs the alignment of BH spin and angular momentum of the disk. The angle between the two-vector θ declines to zero exponentially on an alignment timescale:

$$\begin{eqnarray}&&\theta (t)={\theta }_{0}\exp \left(-\displaystyle \frac{t}{{\tau }_{\mathrm{align}}}\right),\end{eqnarray} \tag{ 13 }$$

where θ₀ is the initial θ at t = 0, and τ_align is the alignment timescale, which is

$$\begin{eqnarray}&&{\tau }_{\mathrm{align}}=1.13\times {10}^{5}{\alpha }_{0.1}^{58/35}{f}_{{\nu }_{2}}^{-5/7}{M}_{6}^{-2/35}{\left(\displaystyle \frac{{f}_{\mathrm{Edd}}}{{\eta }_{0.1}}\right)}^{-\tfrac{32}{35}}{a}_{* }^{5/7}\,\mathrm{yr}.\end{eqnarray} \tag{ 14 }$$

Because we denote the parallel spin with positive a_* and anti-parallel spin with negative a_*, the range of θ is from 0° to 90°.

In the above equation, α_0.1 is the α parameter in a standard thin disk divided by 0.1 (Shakura & Sunyaev 1973); ${f}_{{\nu }_{2}}$ is a nonlinear-effect-related coefficient defined from a simulation; M₆ is the BH mass in units of 10⁶ M_⊙ and ${\eta }_{0.1}\equiv \eta /0.1$. Those parameters are all of order unity, therefore we set them as one in the following Monte-Carlo simulation.

Equation (9) of Yi et al. (2018) gives the mass of an individual BH as a function of its initial mass, orbital separation, and properties of the AGN disk:

$$\begin{eqnarray}&&{m}_{{\rm{f}}}={\left(\displaystyle \frac{{a}_{\mathrm{12,0}}^{4}{f}_{\mathrm{Edd}}\gamma {m}_{0}^{4\gamma }}{3.94\times {10}^{-6}}\right)}^{1/(3+4\gamma )},\end{eqnarray} \tag{ 15 }$$

where m_f and m₀ are the masses of the individual BH in units of M_⊙, at coalescence and initially respectively; a_12,0 is the initial orbital separation of the BBH in units of 10¹² cm; f_Edd is the Eddington ratio of accretion onto individual BHs during inspiral; γ is a parameter used in the numerical simulation of Tang et al. (2017) to describe the mass sink rate. In the above equation, we fix η_0.1 = 1, thus absorbing any change of η into the definition of f_Edd. In Yi et al. (2018), the authors studied equal-mass BBHs. Here we inherent this simplification. The total mass is just twice the individual mass at coalescence to be calculated with Equation (15): M_tot = 2m_f M_⊙.

The effective spin χ_eff is defined as

$$\begin{eqnarray}&&{\chi }_{\mathrm{eff}}=\displaystyle \frac{\left({M}_{1}{a}_{* ,1}\cos {\theta }_{1}+{M}_{2}{a}_{* ,2}\cos {\theta }_{2}\right)}{{M}_{\mathrm{tot}}},\end{eqnarray} \tag{ 16 }$$

where M_i and a_*,i (i = 1, 2) are the mass and spin parameter of individual BHs, respectively, and θ_i is the angle between the spin and the orbital angular momenta.

Given a distribution of initial mass, spin parameters, misalignment angles, Eddington accretion rate of individual BHs, and the initial separation, we can calculate a distribution of the χ_eff and M_tot at merger using Equations (12)–(16). From a distribution of these parameters, the probability density of χ_eff and M_tot at merger can be obtained using the Monte-Carlo method.

We assume that the initial mass function of BHs in the AGN-disk channel is the same as that in clean channels. We use a parameterization proposed by Talbot & Thrane (2018):

$$\begin{eqnarray}&&p({M}_{0})=(1-\lambda ){p}_{\mathrm{pow}}({M}_{0})+\lambda {p}_{\mathrm{Gaussian}}({M}_{0}),\end{eqnarray} \tag{ 17 }$$

where p_pow is a power law with the index −α, a low mass cutoff at M_gap, and a high-mass cutoff at M_cap; p_Gaussian is a Gaussian peak. This peak distribution attributes BHs formed via pulsational pair-instability supernovae (PPSNe). We denote the mean and variance of the peak as m_pp and σ_pp. λ is the portion of PPSNe BHs, which can be estimated through

$$\begin{eqnarray}&&\lambda ={\left(\displaystyle \frac{{m}_{\mathrm{gap}}}{{m}_{\mathrm{cap}}}\right)}^{\alpha -1}.\end{eqnarray} \tag{ 18 }$$

We take the fiducial parameters: M_gap = 5 M_⊙, M_cap = 40 M_⊙, α = 2, m_pp = M_cap, σ_pp = 1 M_⊙. The probability density function (PDF) of the initial mass of individual BHs is plotted as a blue line in the upper panel of Figure 2. The distribution of other parameters is as follows: γ is calculated from τ_sink as (Tang et al. 2017)

$$\begin{eqnarray}&&\gamma =-0.496+1.68{\tau }_{\mathrm{sink}},\end{eqnarray} \tag{ 19 }$$

and τ_sink is sampled uniformly from 3 to 5 (Yi et al. 2018); f_Edd is sampled uniformly in log space from -3 to 1. For the distribution of the binary separation a_12,0, we use the same value used in Yi et al. (2018), i.e., Öpik's law with an upper limit from ionization of the binary, with the upper limit of a_12,0 = 5. The PDF of individual BH masses at merger is plotted as an orange line in the upper panel of Figure 2. The initial a_* of individual BHs is assumed to follow a Gaussian distribution centering at zero and with a width 0.1, which is plotted together with the pdf of a_* at merger in the lower panel of Figure 2. The initial misalignment angle is assumed to be isotropic. Nearly all BHs have their spin axis aligned with the orbital axis at merger, as shown in Figure 3. The time of growth is calculated with

$$\begin{eqnarray}&&{\tau }_{\mathrm{growth}}=\mathrm{log}(x)/{f}_{\mathrm{Edd}}\times {10}^{8}\,\mathrm{yr},\end{eqnarray} \tag{ 20 }$$

We exclude those samples whose τ_growth are larger than the typical lifetimes of AGNs, τ_AGN. Here we use a fiducial value of τ_AGN = 10⁸ yr.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also want to include the effect that the detectable volume has when increasing with the chirp mass of BBH: $V\propto {{ \mathcal M }}^{5/2}$ (in the local universe, where ${ \mathcal M }$ is the chirp mass of the BBH and it is proportional to M_tot). As a result, the detected distribution will be biased toward higher total mass:

$$\begin{eqnarray}&&{p}_{\det }({M}_{\mathrm{tot}},{\chi }_{\mathrm{eff}})\propto p({M}_{\mathrm{tot}},{\chi }_{\mathrm{eff}}){M}_{\mathrm{tot}}^{5/2}.\end{eqnarray} \tag{ 21 }$$

The detector-biased mass PDF of BHs is plotted as a green line in the upper panel of Figure 2.

In Figure 4, the density plot represents the probability density of detecting an AGN-disk channel BH with M_tot and χ_eff with GW observation. We also plot the observed data with their uncertainty in the same figure.

Zoom In Zoom Out Reset image size

It is shown that GW170729 falls far from the cluster of other events in the parameter space, and toward the peak of the predicted distribution of the AGN-disk channel events. This hints that GW170729 is from a different channel. Stone et al. (2017) found that the BBH merger rate density in the AGN-disk channel is ∼3 yr⁻¹ Gpc⁻¹, with large uncertainty. This is implied to occur among all BBH merger GW events; the ratio of AGN-disk events can be a few percent to tens of percent (the BBH merger rate calculated from the LIGO/Virgo O1–O2 runs is 9.7–101 yr⁻¹ Gpc⁻¹ according to Abbott et al. 2019). We suggest a search in the archival data for an FRB around the time of GW170729, and we expect to see some more AGN-disk channel candidates after the O3 run.

4.1. Observation of EM Counterparts and Host Galaxies

4.2. Caveats

S.X.Y. acknowledge support from the Netherlands Organisation for Scientific Research (NWO). K.S.C. is supported by a GRF grant under 17310916.