AGN Disks Harden the Mass Distribution of Stellar-mass Binary Black Hole Mergers

TL;DR

This paper demonstrates that mergers of stellar-mass black holes within AGN disks preferentially select heavier black holes, leading to a top-heavy mass distribution and potentially explaining a significant fraction of gravitational-wave detections.

Contribution

The study introduces a Monte Carlo simulation showing how AGN disks harden the black hole mass function, a novel insight into black hole merger environments.

Findings

AGN disks preferentially select heavier black holes.

The black hole mass function becomes more top-heavy due to AGN influence.

Estimated AGN-assisted merger rate is about 4 Gpc$^{-3}$yr$^{-1}$.

Abstract

The growing number of stellar-mass binary black hole mergers discovered by Advanced LIGO and Advanced Virgo are starting to constrain the binaries' origin and environment. However, we still lack sufficiently accurate modeling of binary formation channels to obtain strong constraints, or to identify sub-populations. One promising formation mechanism that could result in different black hole properties is binaries merging within the accretion disks of Active Galactic Nuclei (AGN). Here we show that the black holes' orbital alignment with the AGN disks preferentially selects heavier black holes. We carry out Monte Carlo simulations of orbital alignment with AGN disks, and find that AGNs harden the initial black hole mass function. Assuming an initial power law mass distribution $M_{\rm bh}^{-\beta}$, we find that the power law index changes by $\Delta \beta\sim1.3$, resulting in a more top-heavy population of merging black holes. This change is independent of the mass of, and accretion rate onto, the supermassive black hole in the center of the AGN. Our simulations predict an AGN-assisted merger rate of $\sim4$Gpc$^{-3}$yr$^{-1}$. With its hardened mass spectra, the AGN channel could be responsible for $10-50$% of gravitational-wave detections.

Figures

Equations

$$\sigma T_{\mathrm{eff}}^4$$

$$T^4$$

$$\tau$$

$${\beta }^b{c}_{\mathrm{s}}^2\Sigma$$

$$p_{\mathrm{rad}}$$

$$p_{\mathrm{gas}}$$

$$\beta$$

$$\Sigma$$

$$c_{\mathrm{s}}$$

$$\kappa (\rho ,T)=\{\begin{array}{rlr}& {\kappa }_0{T}^2& T<100K\\ & {\kappa }_1& T<1000K\\ & \frac{{\kappa }_{\mathrm{m}}-{\kappa }_1}{500K}(T-1000K)+{\kappa }_1& T<1500K\\ & {\kappa }_{\mathrm{m}}+({\kappa }_{{\mathrm{H}}^{-}}^{-1}+({\kappa }_{\mathrm{e}}+{\kappa }_{\mathrm{K}}{)}^{-1}{)}^{-1}& \end{array}$$

$${\kappa }_{H^{-}}\approx 2Z^{0.5}{\rho }_{-11}^{0.5}{T}_4^{7.7}\text{mbox}{g}^{-1}\text{mbox}{cm}^2,$$

$$\begin{gathered} {\kappa }_e=0.2(1+X){\left(1+2.7\times 10^{-8}\frac{{\rho }_{-11}}{T_4^2}\right)}^{-1} \\ \times {\left[1+{\left(\frac{T_4}{4.5\times 10^4}\right)}^{0.86}\right]}^{-1}\text{mbox}{g}^{-1}\text{mbox}{cm}^2 \end{gathered}$$

$${\kappa }_{\mathrm{K}}\approx 4(1+X)(Z+0.001){\rho }_{-11}{T}_4^{-3.5}\text{mbox}{g}^{-1}\text{mbox}{cm}^2.$$

$$\rho =\frac{{\Omega }^2}{2\pi G}$$

$$\hat{\mathbf{n}}=(\sin \psi ,-\cos \theta \cos \psi ,\sin \theta \cos \psi )$$

$$\cos \varphi =\frac{\sin \psi \sin \theta }{\sqrt{{\sin }^2\psi {\sin }^2\theta +{\cos }^2\theta }}$$

$$r_{\mathrm{cross}}=\frac{a(1-e^2)}{1\mp e\cos \varphi }$$

$${\hat{\mathbf{r}}}_{\mathbf{cross}}=\frac{\pm \left(\cos \theta \cos \psi ,\sin \psi ,0\right)}{\sqrt{{\sin }^2\psi {\sin }^2\theta +{\cos }^2\theta }}$$

$$E=-\frac{G{M}_{\bullet }{M}_{\mathrm{bh}}}{2a},$$

$$v_{\mathrm{bh}}=\sqrt{G{M}_{\bullet }}\sqrt{\frac{2}{r}-\frac{1}{a}},$$

$${\hat{\mathbf{v}}}_{\mathbf{bh}}=\mp \cos \alpha \hat{r}+\sin \alpha \hat{r}\times \vec{n}$$

$$v_{\mathrm{bh},\mathrm{z}}=v_{\mathrm{bh}}\sin \alpha \sqrt{{\sin }^2\psi {\sin }^2\theta +{\cos }^2\theta }$$

$${\mathbf{v}}_{\mathbf{gas}}={\left(\frac{G{M}_{\bullet }}{r}\right)}^{1/2}\hat{\mathbf{r}}\times \hat{\mathbf{z}}.$$

$$\Delta v=\sqrt{v_{\mathrm{bh}}^2+v_{\mathrm{gas}}^2-2v_{\mathrm{bh}}{v}_{\mathrm{gas}}\sin \alpha \sin \theta \cos \psi }.$$

$$\Delta M_{\mathrm{cross}}=\Delta v{t}_{\mathrm{cross}}{r}_{\mathrm{BHL}}^2\pi \Sigma /(2H),$$

$$\mathbf{\Delta }{\mathbf{v}}_{\mathbf{bh}}$$

$$\mathbf{\Delta J}$$

$${\mathbf{v}}^{\prime }$$

$$a^{\prime }$$

$$e^{\prime }$$

$${\varphi }^{\prime }$$

$${\mathbf{n}}^{\prime }$$

$${\hat{\mathbf{a}}}^{\prime }$$

$${\theta }^{\prime }$$

$${\psi }^{\prime }$$

$$\frac{dn}{da}\propto a^{-3/2-0.5M_{\mathrm{bh}}/M_{\max }},$$

Full text

1 Introduction

Stellar-mass binary black holes (BBHs) can form in multiple distinct astrophysical sites and processes. Formation channels include isolated binary stellar systems in which both stars end their lives as black holes (Belczynski et al., 2002; Sadowski et al., 2008; Postnov & Yungelson, 2014), and dynamical formation in which the black holes become gravitationally bound following a chance encounter (Portegies Zwart & McMillan, 2000; O’Leary et al., 2007; Miller & Lauburg, 2009; Benacquista & Downing, 2013). With the increasing rate of BBH mergers discovered by Advanced LIGO (Aasi et al., 2015) and Advanced Virgo (Acernese et al., 2015a), our prospects to probe these formation mechanisms is rapidly improving (Abbott et al., 2018a).

Active galactic nuclei (AGNs) represent a promising site for the dynamical formation and/or merger of BBHs. Galactic centers harbor a large, dense population of stellar-mass black holes (BHs) due to mass segregation (Bahcall & Wolf, 1977; Kocsis et al., 2006; O’Leary et al., 2009; Hailey et al., 2018). In active galaxies the orbit of some of these BHs will align with the accretion disk of the central supermassive black hole (SMBH) due to gas damping (Syer et al., 1991; Levin, 2007; McKernan et al., 2012). BHs within the disk can migrate due to angular momentum exchange with the disk, towards so-called migration traps (McKernan et al., 2012). As BHs accumulate near these traps, they can form BBH systems (Bellovary et al., 2016; Secunda et al., 2018). Additionally, a significant fraction of the BHs may already enter the disk as binaries (Bartos et al., 2017b), while some BH binaries could be born within the disk itself (Stone et al., 2017). Once a binary enters or is formed within the SMBH accretion disk, it rapidly ($\lesssim 1$ Myr) merges due to gaseous torques (Bartos et al., 2017b; Stone et al., 2017).

An interesting feature of BBH mergers in AGNs is the possibility of detectable electromagnetic emission. As the BBH system is surrounded by a dense gaseous environment, the BHs can form mini-accretion disks and may reach super-Eddington accretion rates (Bartos et al., 2017b; Stone et al., 2017). Emission from the mini-disks must be very bright, and possibly relativistically beamed, in order to be detected, as the resulting radiation needs to outshine the AGN in which the BBH was embedded in. BBH mergers detected by Advanced LIGO (Aasi et al., 2015) and Advanced Virgo (Acernese et al., 2015b) have been regularly followed-up by electromagnetic and neutrino observations (Abbott et al., 2016; Adrián-Martínez et al., 2016; Albert et al., 2017), increasing the chance that any electromagnetic emission from these mergers will be found. Already, there have been claims of the coincident detection of gamma-ray bursts with BBH mergers (Connaughton et al., 2016; Verrecchia et al., 2017). More coincident observations will help clarify whether these associations were astrophysical.

Apart from electromagnetic counterparts, determining the formation channel of BBH mergers relies on the mergers’ reconstructed properties. An important distinguishing feature can be the host galaxy of the binaries. While the localization of BBH mergers with Earth-based interferometers is mostly inadequate to constrain their origin (Chen & Holz, 2016), a set of the best-localized mergers can be sufficient to establish an association with certain rare host galaxy types (Bartos et al., 2017a). This might be particularly helpful with BBH mergers in AGNs, as AGNs represent only a small fraction of galaxies (Hao et al., 2005).

A powerful differentiator between formation scenarios is the mass and spin distribution of the merging BHs (Arca Sedda & Benacquista, 2019). BBH systems from isolated stellar binaries are expected to favor close to equal masses (Dominik et al., 2012; Mandel & de Mink, 2016), and to have BH spins aligned with the orbital axis (O’Shaughnessy et al., 2017; Gerosa et al., 2018; Kushnir et al., 2017; Roulet & Zaldarriaga, 2019; Antonini et al., 2016; Hoang et al., 2018). Dynamically formed binaries typically have unequal BH masses, a top-heavy mass distribution compared to the initial mass function (O’Leary et al., 2016; Rodriguez et al., 2016), and randomly oriented spins compared to their orbital axis. In particular, spin orientation has been explored as a means to differentiate between formation scenarios (Vitale et al., 2017; Farr et al., 2017). The spin components parallel to the orbital axis for BBH mergers detected by Advanced LIGO are mostly consistent with zero, which suggest either no spin alignment, which is the expected outcome of dynamical formation (Farr et al., 2017). It could also be explained by low BH spins at merger (Hotokezaka & Piran, 2017).

The expected properties of BBH mergers in AGNs are currently poorly understood. McKernan et al. (2018) considered the possible range of astrophysical parameters that determine the merger rate, and found that a very broad range of rates, within $\Gamma=10^{-3}-10^{4}\,$Gpc*-3yr-1*, is possible. They also showed via a simplified model that the migration of the black holes within the disk, and gas accretion can modify the spins and masses of merging black holes. In particular, they pointed out that if the post-merger BH stays within the disk then the mass distribution of merging black holes will shift to significantly heavier objects due to repeated mergers.

In this paper we explore the role of the alignment of the orbital plane of BHs around SMBHs with AGN disk planes on the parameters and rate of BBH mergers. These dependencies have not been explored previously even though it has major influence on BBH mergers. We considered a realistic distribution of SMBH masses, accretion rates onto the SMBH, and multiple initial BH mass distributions. We focused on the scenario in which singleton BHs ground down to alignment with the AGN disk and form binaries in migration traps within the disk.

The paper is organized as follows. In Section 2 we present our analytical framework to calculate the alignment rate of black holes as a function of the properties of the BH, the SMBH and the SMBH accretion disk. In Section 3 we present our findings on BHs in the AGN disk and on the properties of binary mergers within the disk. In Section 4 we briefly describe an alternative source of stellar-mass BHs in AGN disks from star formation within the disk. We conclude in Section 5.

2 Method

2.1 AGN disk

We consider a geometrically thin optically thick, radiatively efficient and steady-state accretion disk (Shakura & Sunyaev, 1973), which is expected in AGNs. We assume the disk has viscosity parameter $\alpha=0.1$, and radiation efficiency $\epsilon=L_{\bullet,Edd}/\dot{M}_{\bullet,Edd}c^{2}$, where $L_{\bullet,Edd}$ and $\dot{M}_{\bullet,Edd}$ are the Eddington luminosity and Eddington accretion rate of a SMBH with mass $M_{\bullet}$. We further assume the accretion rate $\dot{m}_{\bullet}=\dot{M}_{\bullet}/\dot{M}_{\bullet,Edd}$ is constant during the lifetime $\tau_{\rm agn}$ of an AGN. These parameters determine the disk surface density $\Sigma(r)$, scale height $H(r)$, isothermal sound speed $c_{\rm s}(r)$, and midplane temperature $T(r)$, which are dependent on the distance $r$ from the SMBH. To derive these quantities we follow the model of Sirko & Goodman (2003), assuming that viscosity in the disk is proportional to the total pressure:

[TABLE]

Where $T_{\rm eff}(r)$ is the black body temperature of the disc at r. $\dot{M}^{\prime}=\dot{M}_{\bullet}\sqrt{1-r_{\rm min}/r}$, $r_{\rm min}$ is the inner radius of the disc, we set $r_{\rm min}=3r_{\rm s}$ here, where $r_{\rm s}=2GM/c^{2}$. $\Omega$ is the Keplerian angular velocity of the The disc, which is $\sqrt{GM/r^{3}}$. $\tau$ is the optical depth at midplane. We assume the mean molecular mass $m=0.62m_{\rm H}$. $\sigma$ is the Stefan-Boltzmann constant and $k$ is Boltzmann constant. In order to solve these equations, we still need to provide the opacity $\kappa$. We approximate the opacity with (c.f. Fig. 1 in Thompson et al. 2005):

[TABLE]

Where $\kappa_{0}=2.4\times 10^{-4}$ g*-1cm2K-2*, $\kappa_{1}=\kappa_{0}(100$K$)^{2}=2.4$ g*-1*cm2. $\kappa_{m}$ is the molecular opacity, which is about $0.1Z$, where $Z$ indicates mass fraction of heavy elements. We adopt $Z=0.01$ in our simulations (Iglesias & Rogers, 1996). $\kappa_{H^{-}}$ is the opacity due to the negative hydrogen ion:

[TABLE]

where $\rho_{-11}\equiv\rho/(10^{-11}\,\mbox{g\,cm}^{-3})$ and $T_{4}=T/(10^{4}\,\mbox{K})$. The electron scattering opacity $\kappa_{e}$ can be written as

[TABLE]

$X$ is the mass fraction of hydrogen. The Kramers opacity $\kappa_{K}$ is due to free-free, bound-free, and bound-bound electronic transitions:

[TABLE]

The equations above are valid where the Toomre Q-parameter $Q=C_{\rm s}\Omega/(\pi G\Sigma)\gtrsim 1$. Beyond that region, out to 0.1 pc, we assume that the disk self-regulates, which we take into account by fixing the Q-parameter at 1 and replace Eq. (1) with (see Eq. 15 of Sirko & Goodman 2003):

[TABLE]

2.2 Orbital variation due to AGN disk crossing

Consider a Cartesian coordinate system with the SMBH in its center, in which the AGN disk lies in the XY plane, and the major axis of the orbit of the stellar mass BH orbiting the SMBH is in the YZ plane (see Fig. 1). Let $\hat{\bf{a}}=(0,\sin{\theta},\cos{\theta})$ be the unit vector pointing in the direction of the major axis of the BH’s orbit in Cartesian coordinates, with $\theta$ being the angle is the orbit’s inclination angle compared to the AGN disk plane. Let $a$ be the semi-major axis, $e$ the eccentricity and $\hat{\bf{n}}$ is the normal vector of the BH’s orbital plane. Let $\psi\equiv\pi/2+\hat{\bf{n}}\hat{\bf{x}}$ the angle between $\hat{\bf{n}}$ and YZ plane, where $\hat{\bf{x}}$ is the unit vector along the x axis111Here we adopt the coordinates $(\theta,\psi,\phi)$ that are the most suitable for the derivation of the results, even though these are not standard Keplerian orbital elements. They can be converted to Keplerian coordinates as follows. Longitude of the ascending mode: $\arccos[\cos\theta\cos\psi/(\sin^{2}\psi\sin^{2}\theta+\cos^{2}\theta)]$. Argument of periapsis: $\phi$. Inclination: $\sin\theta\sin\psi$.. We have

[TABLE]

Let $\phi$ be the angle between the major axis of the BH orbit and the intersection of the orbital plane and the AGN disk. We find that

[TABLE]

The BH orbit crosses the AGN disk plane at two points. The distance of these crossing points from the focal point of the BH orbit are

[TABLE]

where $-$($+$) corresponds to the farther(closer) point. The normal vectors pointing from the focal point to these two crossing points are

[TABLE]

The energy of the BH orbit is

[TABLE]

where $M_{\bullet}$ is the SMBH mass and $M_{\rm bh}$ is the BH mass. Thus, the speed of the BH at the crossing points is

[TABLE]

where $r$ is the BH-SMBH distance.

Let $\alpha$ be the angle between the BH velocity and radius vectors, so the unit vector in the direction of velocity is

[TABLE]

Therefore, the $z$ component of the BH’s velocity is

[TABLE]

The time of crossing is $t_{\rm cross}=2H/v_{\rm z}$, where H is the scale height and we used the fact that $v_{\rm z}$ does not change significantly in a single crossing. The velocity of the gas in the disk is

[TABLE]

The relative speed of the gas in the disk and BH is

[TABLE]

The mass accreted during the crossing is

[TABLE]

where $r_{\rm BHL}=2GM_{\rm bh}/(\Delta v^{2}+c_{s}^{2})$ is the BH’s Bondi-Hoyle-Lyttleton radius and $\Sigma$ is the surface density of the AGN disk. For Eq. 25 we also made use of the fact that $r_{\rm BHL}\ll H$. We define a dimensionless factor $\lambda\equiv\Delta M_{\rm cross}/M_{\rm bh}$. The change of velocity and angular momentum of the BH after a crossing are

[TABLE]

The crossing slightly changes most orbital parameters. Denoting the new parameters with $()^{\prime}$, we obtain

[TABLE]

These updated values above allow for the update of orbital parameters after each crossing. Repeating this procedure, we obtain the orbital evolution.

2.3 Monte Carlo simulations

We carried out Monte Carlo simulations to obtain the distribution of BH masses and mass ratios in BBH mergers within AGN disks, and the timescale of orbital alignment. We assumed that the orbits of singleton BHs align with the AGN disk, in which they form binaries near migration traps. Following the alignment, we consider the random interaction of two BHs, for simplicity without further interactions of the merged objects. This could occur, for example, if the merger produces a sufficiently large natal kick to move the object out of the disk. Additional interactions will result in heavier BHs within mergers. We assume that all BHs whose orbit aligns with the disk merge once, within a negligible time compared to the lifetime of the AGN disk. We neglect the mass increase of the black holes due to accretion from the AGN disk. We used a sample size of $10^{4}$ for a given parameter combination, except for the case of $M_{\bullet}=10^{8}$ M⊙ for which we used a sample size of 2000.

We considered a BH initial mass function $dN/dM_{\rm bh}\propto M_{\rm bh}^{-\beta}$ with BH masses $M_{\rm bh}\in[5\mbox{M}_{\odot},50\mbox{M}_{\odot}]$. The mass distribution of black holes is currently poorly understood. Black hole masses from Galactic X-ray binaries suggests a soft power law spectrum with $\beta\sim 5$ (Kochanek, 2015; Farr et al., 2011), while LIGO/Virgo’s observations indicate an index of $\beta=0-3$ (Abbott et al., 2018b). Here we adopted the index of the stellar initial mass function, $\beta=2.35$, as out fiducial model, and also considered $\beta=2$ and 3 to examine the $\beta$-dependence of our results.

As this initial BH population undergoes mass segregation (Bahcall & Wolf, 1977; Panamarev et al., 2019), its mass distribution will vary with distance from the SMBH. We take into account this mass segregation in the spatial distribution of BHs. Following O’Leary et al. (2009), we adopt an initial distribution as a function of the orbit’s semi-major axis (Gondán et al., 2018)

[TABLE]

with $M_{\rm max}=50$ M*⊙. We see that higher-mass BHs will typically be closer to the SMBH. We considered the maximal semi-major axis of interest to be the radius of influence of the SMBH $R_{\rm inf}=1.2M_{6}^{1/2}$ pc, where $M_{6}=M_{\bullet}/10^{6}$M⊙*. We further adopted a thermal eccentricity distribution of $n(e)=2e$ (see Gillessen et al. 2009). Finally, we assumed that the AGN disk becomes inhomogeneous (clumpy) beyond a radius $R_{\rm disk}=0.1M_{6}^{1/2}$ pc, once the disk’s self-gravity becomes non-negligible (Haiman et al., 2009). The nature of BH-disk interactions in such inhomogeneous (clumpy) disks is unclear. As massive stars begin forming in the clumping disk, they generate pressure support through winds, radiation pressure and supernova explosions, which can keep the disk vertically supported and allow some gas to continue to accrete through into the inner, stable regions of the disk. Such inhomogeneous (clumpy) regions could still produce sufficient drag to align BH orbits with the disk. This may not reduce, on average, the orbital alignment timescale, but may broaden its distribution (for a fixed initial orbit). Nevertheless, due to the uncertainties in this picture, we conservatively assume that no interaction takes place in this inhomogeneous (clumpy) region of the disk.

Another source of inhomogeneities in the accretion disk is the presence of annular gaps associated with massive embedded objects (Kocsis et al., 2011). While stellar-mass objects are unlikely to open these gaps at distances $r\gtrsim 0.01~{}{\rm pc}$ (Stone et al., 2017), it is theoretically possible for embedded black holes to grow, through accretion, up to the isolation mass, which corresponds to an intermediate-mass black hole (Goodman & Tan, 2004; McKernan et al., 2012). The presence of such a gap would dramatically reduce the drag on any inclined orbits that pass through it, but we neglect this effect because the small geometric size of these gaps means that only a minor fraction of inclined orbits will interact with them at any given time.

Given the initial parameters of an orbit, we can calculate the orbit’s alignment time $\tau_{\rm align}$ with the AGN disk, which we define as $1-\hat{\bf n}\cdot\hat{\bf z}<10^{-3}$. We carried out simulations assuming a fiducial AGN lifetime of $\tau_{\rm agn}=10^{7}$ yr (Martini, 2004) and check which BHs aligned with the disk. We note that the lifetime of lower-luminosity AGNs is believed to be longer (Hopkins et al., 2005a, b). To leading order, the alignment timescale for a stellar-mass BH with orbital period $t_{\rm orb}$ is $\tau_{\rm align}\sim t_{\rm orb}\lambda^{-1}/2$ (Syer et al., 1991; Bartos et al., 2017b).

3 Results

3.1 Black hole mass function

We calculated the mass distribution of BHs that align with the AGN disk for a several initial mass function indices $\beta$ and SMBH masses. We found that after alignment the BH mass distribution still follows a power law, but the mass function index in the AGN disk changes from the initial value by $\Delta\beta\approx-1.3$. We found this change to be independent of both the initial $\beta$ index and the SMBH mass. We show the simulated mass distributions for different cases in Fig. 2.

There are three reasons that explain the flattening of the mass function. The first is mass segregation: since heavier BHs are closer to the SMBH, more of them cross the AGN disk plane within $R_{\rm disk}$ where the disk is continuous. The second is alignment rate: heavier BHs align with the AGN disk more quickly, therefore they have larger chance to finish the evolution within the AGN lifetime. Third, heavier BHs accrete more matter while crossing the disk since $\Delta M_{\rm cross}\propto M_{\rm bh}^{2}$, making them slow down more during each crossing.

3.2 Binary mass ratio

Assuming that binaries are formed within the AGN disk from the random combination of captive BHs in the disk, we calculated the mass ratio $q\equiv M_{1}/M_{2}$ for BBH mergers, with BH masses $M_{1}\leq M_{2}$. We found that the distribution of $q$ peaks around $q\approx 0.2$, with almost no mergers having lower $q<0.2$ and with a shallow distribution for $q>0.2$. The lowest allowed mass ratio in our model is $q=0.1$. This means that the typical mass ratio for our AGN-merger model is low compared to mergers from isolated BH binaries (Dominik et al., 2012; Mandel & de Mink, 2016). The mass ratio distribution we obtain is on average lower than what would be expected from the random pairing of BHs drawn from the initial mass function due to the harder mass spectrum in the AGN.

With the above assumptions the mass ratio distribution is largely independent of the initial mass function index $\beta$ and the SMBH mass, making this distribution universal. Fig. 3 shows the obtained $q$ distributions for different SMBH masses and initial mass functions. We also show in Fig. 4 the 2D mass ratio distribution as a function of total binary mass and mass ratio for the representative case of $M_{\bullet}=10^{7}$ M⊙ and $\beta=2.35$.

In reality, binaries formed in the AGN disk will not represent a random pairing of BHs throughout the disk, as we assumed above. For example as heavier BHs are typically closer to the SMBH, the orbital axis of a heavier BH should typically be closer to another heavier BH than a lower-mass BH, making them preferentially more likely to merge. In this case the resulting mass-ratios will be typically closer to 1 than what we reported above.

3.3 Connection between AGN accretion rate and merger rate

In order to calculate the expected merger rate in AGN, we assume that stellar mass in galactic centers are comparable to $2M_{\bullet}$, and that BHs represent 4% of this mass (Miralda-Escudé & Gould, 2000). The obtained BBH merger rate as a function of $\dot{m}_{\bullet}$ is shown in Figure 5 for different SMBH masses. We find that the merger rate is relatively insensitive of the AGN accretion rate: it only changes by about a factor of 10 over three decades of change in the accretion rate from $\dot{m}_{\bullet}=10^{-3}$ to $\dot{m}_{\bullet}=1$.

The reason for the merger rate’s weak dependence on $\dot{m}_{\bullet}$ is that a significant fraction of the BHs in the galactic center—those which are close to alignment at the beginning—take much less time to become fully aligned than $\tau_{\rm agn}$. While the alignment time of these BHs depends on $\dot{m}_{\bullet}$, they will nevertheless merge within $\tau_{\rm agn}$ independently of the accretion rate for the $\dot{m}_{\bullet}$ range we considered. More specifically, we find the approximate proportionality $\tau_{\rm align}\propto\dot{m}_{\bullet}^{0.4}$. This makes the average merger rate during $\tau_{\rm agn}$ only weakly dependent on $\dot{m}_{\bullet}$. For comparison, the average alignment time for BHs that merge within $\tau_{\rm agn}=10^{7}$ yr in an AGN disk with $\dot{m}_{\bullet}=0.1$ and $M_{\bullet}=10^{7}$ M⊙ is 4.4 Myr (see also Syer et al. 1991; Bartos et al. 2017b).

This means that the merger rate does not depend strongly on the AGNs’ actual accretion rate, therefore knowing the accretion rate for different galaxies does not help differentiate between binary BH formation channels.

3.4 Merger rate as a function of SMBH mass

Since $M_{\bullet}$ affects the AGN disk parameters and the size of the population of BHs close to the disk, the AGN alignment rate will depend on $M_{\bullet}$. We characterized this dependence by calculating the BBH merger rate as a function of $M_{\bullet}$ for the range $M_{\bullet}\in[10^{5}\mbox{M}_{\odot},10^{9}\mbox{M}_{\odot}]$.

Our results are shown in Figure 6. We see that the alignment rate continuously, albeit slowly, increases with $M_{\bullet}$.

3.5 The overall merger rate in our local universe

The properties of AGN disks and stellar-mass BH populations in galactic nuclei are not strongly constrained observationally, and the astrophysical parameters determining the merger rate can vary in a broad range (McKernan et al., 2018). Therefore, here we aim to give an estimate of the expected BBH merger rate in accretion disks based on our results, using our fiducial parameters.

As we find that low AGN accretion rate does not diminish BH alignment with the disk, we considered the broader group of Seyfert galaxies as opposed to focusing only on highly accreting quasars as was done in previous work (Bartos et al., 2017b). We adopted a number density of $n_{\rm Seyfert}=0.018$ Mpc*-3* for Seyfert galaxies (Hao et al., 2005). This number density includes galaxies as faint as $M_{\rm r}=-14$. As observations indicate that the SMBH accretion rate for Seyfert galaxies varies within $\dot{m}_{\bullet}=10^{-3}-1$ (Woo & Urry, 2002), and the merger rate only changes by a factor of few between the two extreme values, for simplicity we adopted a fiducial accretion rate of $\dot{m}_{\bullet}=0.1$ for all Seyfert galaxies. Further, we considered a fiducial Seyfert SMBH mass of $M_{\bullet}=10^{6}$ M*⊙*. These choices characterize the Seyfert population and allow us not to rely on the uncertain distribution of these parameters.

For these choices, we find a fiducial overall BBH merger rate of $\sim 4$ Gpc*-3yr-1*. With the limited lifetime of AGN disks we assumed here, this merger rate does not appreciably deplete the BH population in galactic nuclei.

The hardening of the BH mass spectrum by AGNs also means that the LIGO-Virgo detection rate of AGN-assisted mergers will be higher than for a similar merger rate with BHs drawn from the initial mass spectrum.

4 Star formation in AGN disks

At radii $\gtrsim 0.1M_{6}^{1/2}$ pc from the SMBH, self-gravity becomes important in AGN disks, resulting in fragmentation (Goodman, 2003). Some of these fragments will form massive stars, that will in turn produce BH within the AGN disk. Some of these BHs may be born into binaries, which can be hardened to the point of merger by hydrodynamic drag with the surrounding AGN gas (Stone et al., 2017). The disk of massive young O and B stars in the central parsec of the Milky Way provides local evidence for a recent episode of disk-mode star formation (Levin & Beloborodov, 2003).

The mass function of stars born from fragments in AGN disks is observed, at least in our own Galactic Center, to have a zero-age main sequence spectral slope of $\sim 1.7\pm 0.2$, i.e. harder than the Salpeter index of 2.35 (Lu et al., 2013). This top-heavy initial mass function is clearly conducive to compact remnant formation, and all else equal, will produce a shallower BH mass distribution.

However, the mass spectrum of BHs produced in AGN disks will depend both on this harder stellar mass spectrum and stellar metallicity. If higher metallicity stars are produced in AGN disks, this can make the BH mass spectrum significantly steeper (Belczynski et al., 2010).

One important distinction between the in situ formation channel, and the calculations of §3, is the BBH merger rate dependence on $\dot{m}_{\bullet}$. Rates of BBH production due to fragmentation of Toomre unstable disks scale linearly with the $\dot{m}_{\bullet}$ at high accretion rates, and have an even steeper dependence at significantly sub-Eddington accretion rates (Thompson et al., 2005). If BBHs are predominantly formed in situ, their merger rate $\Gamma$ will thus depend on a steep power of the AGN accretion rate, and the volume-averaged merger rate will be dominated by bright quasars rather than the low-luminosity Seyferts we have found to be of importance for realigning pre-existing BHs into the disk.

5 Conclusion

We carried out Monte Carlo simulations of the interaction of stellar-mass BHs and AGN disks to study the effect of these interactions on the properties of BBH mergers. We explored the effect of the initial BH mass function, the mass of the supermassive black hole in the center of the AGN, and the AGN accretion rate. Specifically, we focus on the scenario where singleton stellar-mass BHs from a pre-existing, spherical and isotropic nuclear star clusters are ground down into alignment with the AGN disk, and then efficiently assemble into binaries in disk migration traps (Bellovary et al., 2016). Our findings are summarized in the following:

• •

The BH mass function hardens due to the mass-dependent impact of the AGN disk. Assuming a power-law initial mass function $M_{\rm bh}^{-\beta}$, the index of the BH mass distribution hardens by $\Delta\beta\approx 1.3$. This change is independent of $\beta$ at least within the $\beta\in[2,3]$ range we tested. The change is also independent of $M_{\bullet}$. This means that independently of the merger distribution in different galaxies, the overall effect will be similar.

• •

The BBH merger rate in AGN shows no dependence on the SMBH mass.

• •

For our fiducial parameters we find that the BBH merger rate in AGN is $\Gamma\sim 4$ Gpc*-3yr-1*, i.e $4-40\%$ of the total BBH merger rate (Abbott et al., 2018a, b). Given the harder mass spectrum, this corresponds to a fraction of $\sim 10-50$% of BBH detections by LIGO/Virgo, assuming that the rest of the detections come from a population of black holes drawn from an initial mass spectrum with $\beta=2.35$.

• •

BHs formed within AGN disks may add a comparable number of mergers, further increasing the AGN contribution to the total BBH merger rate (Stone et al., 2017). This formation channel is strongly correlated with $\dot{m}_{\bullet}$, making it differentiable from BH mergers from disk-alignment.

• •

The merger rate of BBHs produced from BHs aligned by gas drag into AGN disks is only weakly dependent on the SMBH accretion rate. This means that AGNs with lower accretion rates will contribute a substantial part of the overall BBH merger rate from AGNs. This enhances the possibility of finding electromagnetic emission from some BBH mergers as lower-accreting AGNs represent weaker background radiation.

• •

Because the merger rate we find has only a weak dependence on $\dot{m}_{\bullet}$ and $M_{\bullet}$, the volumetric rate of AGN-aligned BBH mergers will be dominated by the more abundant low-mass Seyferts. The relatively high volumetric density of these galaxies will necessitate more BBH detections for host galaxy identification techniques than expected if other factors (e.g. in situ BH formation) we have neglected bias the BBH merger rate in AGNs towards a sub-population of high-$\dot{m}_{\bullet}$ or high-$M_{\bullet}$ galaxies (Bartos et al., 2017a).

Future work to refine this analysis will additionally needs to take into account that some BHs are already in binaries before moving into the AGN disk, changing the mass and mass ratio distributions. Spin and the speed of migration within the disk can also affect the final distributions, as can repeated mergers within migration traps (Secunda et al., 2018). It will also be useful to better understand how AGN-assisted mergers can be differentiated from other dynamical formation channels that also result in hardened mass spectra (O’Leary et al., 2016), for example by using their localization (Bartos et al., 2017a; Corley et al., 2019).

The authors thank Brian Metzger and Saavik Ford for useful discussions. The authors are thankful for the generous support of the University of Florida and Columbia University in the City of New York. The Columbia Experimental Gravity group is grateful for the generous support of the National Science Foundation under grant PHY-1708028. This work received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No 638435 (GalNUC) and was supported by the Hungarian National Research, Development, and Innovation Office grant NKFIH KH-125675 (to BK).