Cosmic Evolution of Stellar-mass Black Hole Merger Rate in Active Galactic Nuclei

Observations by the LIGO and Virgo gravitational-wave detectors (Acernese et al. 2014; Aasi et al. 2015) show a high rate of stellar-mass black hole mergers of ∼10–100 Gpc⁻³ yr⁻¹ (Abbott et al. 2019a). Despite the rapidly growing number of detections, however, the origin of these black hole mergers is currently not known. Possible formation mechanisms include isolated stellar binary evolution (Dominik et al. 2012; Kinugawa et al. 2014; Giacobbo et al. 2018; Bavera et al. 2019) and chance encounters in dense stellar clusters such as galactic nuclei or globular clusters (Portegies Zwart & McMillan 2000; O'Leary et al. 2006; Samsing et al. 2014; Rodriguez et al. 2015; O'Leary et al. 2016; Zhang et al. 2019).

Active galactic nuclei (AGNs) represent a unique environment in which the interaction of a dense cusp of stellar-mass black holes near the galactic center (Hailey et al. 2018) and the dense AGN accretion disk result in dramatically altered merger rates and properties (McKernan et al. 2014; Bartos et al. 2017b; Stone et al. 2017). As black holes orbiting the central, supermassive black hole cross the AGN disk, they experience friction that can align their orbit with the disk (Bartos et al. 2017b). Additional black holes can be born within the AGN disk due to gravitational fragmentation (Stone et al. 2017). Once in the disk, black holes migrate inward where their density further increases, enabling a higher rate of interaction (Tagawa et al. 2019).

AGN disks act as black hole assembly lines that collect and concentrate black holes in small volumes, enhancing their merger rate. The resulting rate can be a significant fraction of the total merger rate observed by LIGO and Virgo with estimates ranging within ∼10⁻³–10² Gpc⁻³ yr⁻¹ (McKernan et al. 2018, 2020; Tagawa et al. 2019; Yang et al. 2019b).

The properties of black holes within AGN-assisted mergers are expected to be different from other formation channels, enabling observational probes. Heavier black holes will be overrepresented in mergers within AGN disks compared to the black hole initial mass function (IMF) as they can more efficiently align their orbit with the disk (Yang et al. 2019b). In addition, as multiple black holes are driven toward the same small volume within the disk, consecutive mergers of the same black hole, or so-called hierarchical mergers, are common (Yang et al. 2019a). Such hierarchical mergers will result in characteristic high black hole spins, which will typically be aligned or antialigned with the binary orbit (Tagawa et al. 2019; Yang et al. 2019a). Two particular black hole mergers recorded by LIGO-Virgo so far, GW170729 (Abbott et al. 2019a) and GW170817A (Zackay et al. 2019), have the characteristically high mass and aligned spin expected from hierarchical mergers in AGN disks, albeit they are also consistent with other hierarchical formation channels (Yang et al. 2019a; Gayathri et al. 2020). Additional heavy black hole candidates have been reported that could have also originated from this channel and will require further investigation (Udall et al. 2019; Nitz et al. 2020).

Beyond the properties of the black holes themselves, AGNs provide other means to probe this population. Mergers in AGN disks are only expected in galaxies with active nuclei, which can be used to statistically differentiate them from other formation channels (Bartos et al. 2017a). Additionally, merger in a gas-rich environment may produce detectable electromagnetic emission (Bartos et al. 2017b; McKernan et al. 2019; Yi & Cheng 2019).

Here we investigated a distinct property of a binary merger population: its rate evolution with redshift. Redshift dependence can be used to differentiate between different formation channels and to better understand the underlying mechanisms that result in binary formation and merger. While the binary merger rate's redshift dependence has been previously explored for different binary formation scenarios (e.g., Fishbach et al. 2018), our analysis is the first such investigation for the AGN channel.

The paper is organized as follows. In Section 2 we examine the expected redshift evolution for different formation channels. In Section 3 we discuss the conversion of rate densities to detection rates. In Section 4 we present our results. We conclude in Section 5.

In this section we compute the expected cosmic evolution of the black hole merger rate density for the AGN channel, and review the expected evolution for field binaries and globular clusters from the literature.

2.1. The AGN Channel

The black hole merger rate in AGNs is proportional to the AGN population density n_AGN(z), where z is redshift. It can be evaluated through the AGN luminosity function (LF) ϕ_L(L, z), which is a function of redshift and AGN luminosity L. The bolometric AGN LF can be fitted as (Shen et al. 2020):

$$\begin{eqnarray}&&{\phi }_{{\rm{L}}}(L,z)=\displaystyle \frac{{\phi }_{* }(z)}{{\left[\tfrac{L}{{L}_{* }(z)}\right]}^{{\gamma }_{1}(z)}+{\left[\tfrac{L}{{L}_{* }(z)}\right]}^{{\gamma }_{2}(z)}}{\mathrm{Mpc}}^{-3},\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&{\gamma }_{1}(z)={a}_{0}{T}_{0}(1+z)+{a}_{1}{T}_{1}(1+z)+{a}_{2}{T}_{2}(1+z)\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\gamma }_{2}(z)=\displaystyle \frac{2{b}_{0}}{{\left(\tfrac{1+z}{3}\right)}^{{b}_{1}}+{\left(\tfrac{1+z}{3}\right)}^{{b}_{2}}}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\mathrm{log}{\unicode{x00141}}_{* }(z)=\displaystyle \frac{2{c}_{0}}{{\left(\tfrac{1+z}{3}\right)}^{{c}_{1}}+{\left(\tfrac{1+z}{3}\right)}^{{c}_{2}}}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\mathrm{log}{\phi }_{* }(z)={d}_{0}{T}_{0}(1+z)+{d}_{1}{T}_{1}(1+z),\end{eqnarray} \tag{ 5 }$$

where T_n is the nth order Chebyshev polynomial. The best fit for the 11 parameters in this model are $\{{a}_{0},{a}_{1},{a}_{2};{b}_{0},{b}_{1},{b}_{2};{c}_{0},{c}_{1},{c}_{2};{d}_{0},{d}_{1}\}$ = {0.8396, −0.2519, 0.0198; 2.5432, −1.0528, 1.1284; 13.0124, −0.5777, 0.4545; −3.5148, −0.4045}.

The direct integration of ϕ_L(L, z) will yield the AGN density n_AGN(z). However, the lower end of the LF is subject to large uncertainty, thus we introduced a cutoff L_min when integrating. On the other hand, the mass of SMBHs can be correlated with the AGN luminosity via:

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{\bullet }}{{M}_{\odot }}=3.17\times {10}^{-5}\displaystyle \frac{1-\epsilon }{\dot{m}}\displaystyle \frac{L}{{L}_{\odot }},\end{eqnarray} \tag{ 6 }$$

where is the radiation efficiency of the SMBH, $\dot{m}=\dot{{M}_{\bullet }}/{\dot{M}}_{\mathrm{Edd}}$, $\dot{{M}_{\bullet }}$ is the accretion rate of the SMBH, and ${\dot{M}}_{\mathrm{Edd}}={L}_{\mathrm{Edd}}/\epsilon {c}^{2}$ is the Eddington rate.

Equation (6) could be rewritten to give a relation between the normalized accretion rate $\dot{m}$ and Eddington ratio $\lambda =L/{L}_{\mathrm{Edd}}$:

$$\begin{eqnarray}&&\dot{m}=(1-\epsilon )\lambda .\end{eqnarray} \tag{ 7 }$$

The Eddington ratio λ is found to take the form (Tucci & Volonteri 2017):

$$\begin{eqnarray}&&P(\lambda | L,z)={f}_{\mathrm{uno}}{P}_{1}(\lambda | z)+{f}_{\mathrm{obs}}{P}_{2}(\lambda | z),\end{eqnarray} \tag{ 8 }$$

where ${f}_{\mathrm{uno}}=1-{f}_{\mathrm{obs}}$ is the fraction of unobscured (type-1) AGN and f_obs is the fraction of obscured (type-2) AGN. P₁ and P₂ are the Eddington ratio distributions of type-1 and type-2 AGNs, respectively.

${P}_{1}(\lambda | z)$ follows a log-normal distribution:

$$\begin{eqnarray}&&{P}_{1}(\lambda | z)=\displaystyle \frac{1}{2\pi \sigma (z)\lambda }{e}^{-{[\mathrm{ln}\lambda -\mathrm{ln}{\lambda }_{{\rm{c}}}(z)]}^{2}/2{\sigma }^{2}(z)},\end{eqnarray} \tag{ 9 }$$

with $\mathrm{log}{\lambda }_{{\rm{c}}}(z)=\max (-1.9+0.45z,\mathrm{log}(0.03))$ and $\sigma (z)=\max (1.03-0.15z,0.6)$.

${P}_{2}(\lambda | z)$ follows a gamma distribution with a cutoff at low-Eddington luminosities:

$$\begin{eqnarray}&&{P}_{2}(\lambda | z)={N}_{2}(z){\lambda }^{\alpha (z)}{e}^{-\lambda /{\lambda }_{0}},\end{eqnarray} \tag{ 10 }$$

where λ₀ = 1.5 and N₂(z) is the normalization factor. The slope of the power-law part, α(z), takes the form:

$$\begin{eqnarray}\alpha (z)=\left\{\begin{array}{ll}-0.6 & z\lt 0.6\\ -0.6/(0.4+z) & z\geqslant 0.6.\end{array}\right.\end{eqnarray} \tag{ 11 }$$

We assume that the Eddington ratio distribution has a cutoff (λ_l) at low-Eddington luminosities, which is fixed to be 10⁻⁴ in our study. Our results below are not sensitive to this choice, which we confirmed for the ${\lambda }_{{\rm{l}}}={10}^{-4}\mbox{--}{10}^{-2}$ range.

The fraction f_obs can be parameterized as a function of X-ray luminosity L_X and redshift (Ueda et al. 2014):

$$\begin{eqnarray}&&{f}_{\mathrm{obs}}=\displaystyle \frac{(1+{f}_{\mathrm{CTK}})\psi ({L}_{{\rm{X}}},z)}{1+{f}_{\mathrm{CTK}}\psi ({L}_{{\rm{X}}},z)},\end{eqnarray} \tag{ 12 }$$

where f_CTK is the relative number density of compton thick (CTK, $\mathrm{log}{N}_{{\rm{H}}}\gt 24$) AGNs to that of compton thin (CTN, $\mathrm{log}{N}_{{\rm{H}}}=20\mbox{--}24$) AGNs, N_H is the neutral hydrogen column density in unit of cm⁻². We assume f_CTK = 1 is this work. ψ is the fraction of obscured AGNs ($\mathrm{log}{N}_{{\rm{H}}}=20\mbox{--}22$) in total CTN AGNs and can be expressed as:

$$\begin{eqnarray}\begin{array}{rcl}\psi ({L}_{{\rm{X}}},z) & = & \min ({\psi }_{\max },\max ({\psi }_{43.75}(z)\\ & & -\beta (\mathrm{log}{L}_{{\rm{X}}}-43.75),{\psi }_{\min })),\end{array}\end{eqnarray} \tag{ 13 }$$

where we adopt ${\psi }_{\max }=0.84$, ${\psi }_{\min }=0.2$, and $\beta =0.24$. ψ_43.75(z) can be written as:

$$\begin{eqnarray}{\psi }_{43.75}(z)=\left\{\begin{array}{ll}0.43{\left(1+z\right)}^{0.48} & z\lt 2\\ 0.43{\left(1+2\right)}^{0.48} & z\geqslant 2.\end{array}\right.\end{eqnarray} \tag{ 14 }$$

Since ${f}_{\mathrm{obs}}$ is dependent on X-ray luminosity, we need to convert the bolometric luminosity to the X-ray luminosity using a bolometric correction (Marconi et al. 2004):

$$\begin{eqnarray}&&\mathrm{log}(L/{L}_{{\rm{X}}})=1.54+0.24\xi +0.012{\xi }^{2}-0.0015{\xi }^{3},\end{eqnarray} \tag{ 15 }$$

with $\xi =\mathrm{log}L/{L}_{\odot }-12$. lower limit of the Eddington ratio. The distribution function of λ is independent of the redshift; therefore, the average Eddington ratio is a constant for z ≲ 1, as shown by other studies (Georgakakis et al. 2017).

We additionally need the BH merger rate (Γ) of a single AGN. Yang et al. (2019b) showed that several factors affect the BH merger rate, of which $\dot{m}$ is the dominant one. We assumed that the mean number (N_disk) of stellar black holes in the AGN disk is a univariate function of $\dot{m}$ and obtained a power-law fit (${M}_{\bullet }={10}^{6}{M}_{\odot },\epsilon =0.1$):

$$\begin{eqnarray}&&{N}_{\mathrm{disk}}(\dot{m})=5.5\,{\dot{m}}^{1/3}.\end{eqnarray} \tag{ 16 }$$

We assume that the black holes in AGN disks will hierarchically merge in the migration traps (Yang et al. 2019a) and the number of stellar black holes in AGN disks follows a Poisson distribution with a mean value of N_disk. Consequently, the average BH merger rate is:

$$\begin{eqnarray}&&{\rm{\Gamma }}(\dot{m})=({N}_{\mathrm{disk}}(\dot{m})-1+{e}^{-{N}_{\mathrm{disk}}(\dot{m})})/{\tau }_{\mathrm{AGN}},\end{eqnarray} \tag{ 17 }$$

where ${\tau }_{\mathrm{AGN}}={10}^{7}\,\mathrm{yr}$ is the AGN lifetime.

We note that this dependence is different for several subdominant merger processes associated with AGNs (see Tagawa et al. 2019 for a comparison of different contributing processes); however, we found that our results below only weakly depend on the specific form of ${\rm{\Gamma }}(\dot{m})$. Therefore, other AGN-related processes will not alter our results below.

Combining these above factors, we arrived at an expression for the redshift-dependence of the cosmic BH merger rate in AGN:

$$\begin{eqnarray}&&{R}_{\mathrm{AGN}}(z)={\int }_{L\in {I}_{{\rm{L}}}}{\phi }_{{\rm{L}}}d\mathrm{log}L{\int }_{{\lambda }_{{\rm{l}}}}^{1}{\rm{\Gamma }}(\dot{m})P(\lambda | L,z)d\lambda .\end{eqnarray} \tag{ 18 }$$

Here, the integral domain of the luminosity is ${I}_{{\rm{L}}}=[{L}_{\min },3.15\times {10}^{14}{L}_{\odot }]$.

We show R_AGN(z) in Figure 1 for multiple choices of L_min. We see that the distribution only weakly depends on z, with a maximum around z = 0.8 which is about a factor of two greater than the minimum at z = 0. We also see that the choice of Ł_min does not meaningfully affect the normalized redshift distribution of ${R}_{\mathrm{AGN}}(z)$, although it does change the magnitude of merger rate density. Therefore, in the following, we adopted ${L}_{\min }={10}^{41}\,\mathrm{erg}\ {{\rm{s}}}^{-1}$ as our fiducial model.

Zoom In Zoom Out Reset image size

2.2. Field Binaries

We also compared our model with other formation channels. The first considered mechanism is the field binary channel. We assumed that the formation rate density of field binaries in comoving volume follows the low-metallicity star formation rate (SFR; Fishbach et al. 2018):

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{\mathrm{FB}}(z) & \propto & {\psi }_{\mathrm{MD}}(z){f}_{{\rm{Z}}}(z)\\ & \propto & \displaystyle \frac{{\left(1+z\right)}^{2.7}}{1+{\left(\tfrac{1+z}{2.9}\right)}^{5.6}}\gamma (0.84,{\left(\displaystyle \frac{Z}{{Z}_{\odot }}\right)}^{2}{10}^{0.3z}),\end{array}\end{eqnarray} \tag{ 19 }$$

where ψ_MD(z) is the Madau–Dickinson SFR (Madau & Dickinson 2014) and f_Z(z) is the fraction of star formation occurring at metallicity smaller than Z, γ(s, x) is the lower incomplete gamma function.

Since the binaries do not merge immediately after their formation, we adopted a time-delay model to estimate the merger rate density:

$$\begin{eqnarray}&&{R}_{\mathrm{FB}}={\int }_{{z}_{{\rm{m}}}}^{\infty }{\rho }_{\mathrm{FB}}({z}_{{\rm{f}}})p(t({z}_{{\rm{f}}})-t({z}_{{\rm{m}}}))\displaystyle \frac{{dt}}{{{dz}}_{{\rm{f}}}}{{dz}}_{{\rm{f}}},\end{eqnarray} \tag{ 20 }$$

where p(t) is the distribution of time delay and t(z) is the cosmological look back time. We assumed that the time delay had a flat distribution in log space between 50 Myr and 15 Gyr. The mass distribution was assumed to be (Fishbach et al. 2018):

$$\begin{eqnarray}&&{f}_{\mathrm{FB}}({m}_{1},{m}_{2})=\displaystyle \frac{1-\alpha }{{M}_{\max }^{1-\alpha }-{\left(5{M}_{\odot }\right)}^{1-\alpha }}\displaystyle \frac{{m}_{1}^{-\alpha }}{{m}_{1}-5{M}_{\odot }},\end{eqnarray} \tag{ 21 }$$

where $5{M}_{\odot }\lt {m}_{2}\lt {m}_{1}\lt {M}_{\max }$.

2.3. Globular Clusters

In order to compare the expected redshift evolutions to black hole merger observations via gravitational waves, we need to convert the expected merger rate distribution R_merger to detection rate distribution R_det using the sensitive distance range of LIGO-Virgo.

$$\begin{eqnarray}&&{R}_{\det }(z)=\displaystyle \frac{{R}_{\mathrm{merger}}}{1+z}\displaystyle \frac{{{dV}}_{{\rm{c}}}}{{dz}}\int {P}_{\det }({ \mathcal M })f({m}_{1},{m}_{2}){{dm}}_{1}{{dm}}_{2}.\end{eqnarray} \tag{ 22 }$$

Here, V_c is the comoving volume and $f({m}_{1},{m}_{2})$ is the mass function of binary black holes in AGN disks. ${P}_{\det }({ \mathcal M })$ is the probability of an event at redshift z with detector-frame chirp mass

$$\begin{eqnarray}&&{ \mathcal M }=(1+z)\displaystyle \frac{{\left({m}_{1}{m}_{2}\right)}^{3/5}}{{\left({m}_{1}+{m}_{2}\right)}^{1/5}}\end{eqnarray} \tag{ 23 }$$

being detected by advanced LIGO. In our calculations, we assumed that an event is detectable when its signal-to-noise ratio (S/N) is greater than 8. The S/N of a black hole merger is (O'Shaughnessy et al. 2010):

$$\begin{eqnarray}&&\rho =8\omega \displaystyle \frac{370\,{\rm{Mpc}}}{{D}_{{\rm{L}}}(z)}{\left(\displaystyle \frac{{ \mathcal M }}{{M}_{\odot }}\right)}^{5/6},\end{eqnarray} \tag{ 24 }$$

where D_L(z) is the luminosity distance corresponding to redshift z, and ω ∈ [0, 1] is a geometrical factor determined by the position and orientation of the binary system.⁶ We adopt the noise power spectral density of LIGO at its design sensitivity (Ajith 2011; Martynov et al. 2016) in the evaluation of the S/N above.

Yang et al. (2019b) found that the AGN disk can significantly change the IMF of merging black holes. The hierarchical black hole mergers will then alter the mass distribution of binary black holes (Yang et al. 2019a) since the mass of one component (the remnant of the previous merger) of the binary will increase as the hierarchical merging process continues. We adopted the weighted average binary mass distribution in their work(${N}_{\mathrm{disk}}=2.5$):

$$\begin{eqnarray}&&\bar{f}({m}_{1},{m}_{2})=\sum _{{\rm{n}}=1}^{\infty }{P}_{{\rm{n}}}{f}_{{\rm{n}}}({m}_{1},{m}_{2}),\end{eqnarray} \tag{ 25 }$$

where P_n and f_n are the fraction and mass distribution of nth generation, respectively.

We calculated the expected redshift distributions of detected events for the models described in Section 2 using the conversion described in Section 3. We then compared these distributions to the reconstructed cosmic evolution for black hole mergers observed through gravitational waves by LIGO-Virgo during the O1 and O2 observing periods (Model B in Abbott et al. 2019b).

Our results are shown in Figure 2. We see that the expected rate evolution for LIGO-Virgo observations is currently uncertain and is essentially consistent with all three formation channel models considered here. Looking at the LIGO-Virgo distribution using its expected value, we see that the observed distribution peaks at a higher redshift than the field-binary and globular cluster channels, but at a lower redshift than our AGN model. Taking this expected distribution at face value, the observed distribution is consistent with the linear combination of an AGN and a field-binary population, with 40% and 60% contributions, respectively.

Zoom In Zoom Out Reset image size

A more precise determination of the binaries' origin can be done through the more accurate measurement of their redshift evolution. Current uncertainty on the measured evolution includes all three models considered here. However, looking at Figure 2, we see that the AGN channel falls on the lowest end of this uncertainty region, while field binaries and globular clusters are at the very high end. Therefore, even a modest improvement on the precision of cosmic evolution should be able to rule out either AGNs or the other two channels as the sole origin of the observed binaries. Additionally, for differentiating between formation channels, other information will be useful, such as black hole masses, spins, and host galaxies. We therefore anticipate that with the rapid increase in the number of detections, differentiation between formation channels on the population level will be possible soon.

We computed the expected redshift distribution of the merger rate of stellar-mass black hole mergers in AGN disks. We found that the distribution is close to being uniform out to z ≈ 1, which is distinct from our expectations for field binaries and for some other dynamical merger scenarios, such as in globular clusters. This distinct evolution, together with other differences, can help differentiate between the possible origins of binary mergers and help probe their environment.

The authors are thankful to the University of Florida and Columbia University in the City of New York for their generous support. The Columbia Experimental Gravity group is grateful for the generous support of the National Science Foundation under grant PHY-1708028. This project was supported by funds from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No. 638435 (GalNUC) and by the Hungarian National Research, Development, and Innovation Office grant NKFIH KH-125675 (to B.K. and H.T.). Z.H. acknowledges support from NASA grant NNX15AB19G and NSF grant 1715661.