© ESO 2020

1. Introduction

Following the first observation of gravitational waves (GW) in 2015 (Abbott et al. 2016), GW events are now regularly detected by LIGO/Virgo. The majority of the detections originate from the merger of two stellar-mass black holes or binary black hole (BBH) mergers. Ten confirmed BBH mergers were observed in the first and second observing runs of Advanced LIGO and Advanced Virgo, leading to an estimated BBH merger rate density in the range ℛ = (9.7−101) Gpc⁻³ yr⁻¹ (Abbott et al. 2019a), with an updated estimate of $\mathcal{R}=53.2_{-28.8}^{+58.5}{\mathrm{Gpc}}^{-3}{\mathrm{yr}}^{-1}$ (at 90% credibility) (Abbott et al. 2019b).

Two main BBH formation scenarios are currently debated in the literature: isolated evolution in galactic fields and dynamical formation in dense stellar environments (e.g. see the recent reviews by Mandel & Farmer 2018; Mapelli 2018, and references therein). The predicted merger rates are broadly similar for the two formation channels, but the expected spin and eccentricity distributions are different and may help discriminate between the two paths. For instance, isolated binaries formed in galactic fields are likely to have near-zero eccentricity close to coalescence, while non-negligible eccentricities are expected for binaries formed through dynamical interactions in dense star clusters (Breivik et al. 2016; Samsing & Ramirez-Ruiz 2017).

In galactic nuclei a large concentration of stellar remnants can be found as a result of dynamical friction and mass segregation. Indeed, a density cusp of stellar-mass black holes has been recently uncovered in the central parsec of our own Galaxy (Hailey et al. 2018), suggesting that ∼10⁴ BHs may be accumulated in the Galactic centre (Generozov et al. 2018). Galactic nuclei are also known to harbour super-massive black holes (SMBHs) with typical masses of M_SMBH ∼ (10⁶−10⁹) M_⊙ (e.g. Kormendy & Ho 2013, and references therein). Accretion of matter from the surrounding environment leads to the formation of an accretion disc around the central SMBH, powering the active galactic nucleus (AGN). Some of the stars orbiting in the nuclear star cluster can then be ground down into the AGN accretion disc, with a significant fraction ending up in binary systems.

In this context, an alternative formation path for binary black holes was recently discussed in the literature: BBH formation within AGN accretion discs (Stone et al. 2017; Bartos et al. 2017; McKernan et al. 2018). This novel formation channel has exceptional physical implications: BBH mergers can be considerably accelerated in the gaseous environment of AGN accretion discs and this is the only scenario in which we can plausibly expect electromagnetic counterparts from BBH mergers. The gaseous torques operating in the AGN accretion disc drive the binary into the regime where GW inspiral can take over, eventually leading to coalescence. McKernan et al. (2018) previously assumed that the BBHs merge within the AGN disc lifetime and estimated by this means a merger rate density of ℛ = (10⁻³−10⁴) Gpc⁻³ yr⁻¹. A recent study by Tagawa et al. (2019) combined N-body simulations with semi-analytical tools to investigate the formation, disruption, and evolution of binaries in AGN discs. Therein, these latter authors found a BBH merger rate density of ℛ = (0.02−60) Gpc⁻³ yr⁻¹.

In this paper we study the evolution of BBHs embedded in AGN accretion discs. A simple analytical model of the disc–binary interaction is established using existing numerical studies as an inspiration. We then couple the disc-driven evolution equations to the corresponding equations of GW-emission. We compute the merger rate density that is expected from this coupled channel based on a weighted average of the merger timescales of a population of BBHs radially distributed in the AGN accretion disc. We analyse the physical dependence of the merger rate density on both the binary orbit and the accretion disc parameters, emphasising the important role of the eccentricity. The paper is structured as follows. In Sect. 2, we summarise the main equations governing the BBH evolution in the disc-driven and GW-driven regimes. In Sect. 3, we present a detailed computation of the BBH merger rate density estimates for this coupled evolution channel. We analyse the dependence of the resulting merger rate on the underlying physical parameters (Sect. 4). We discuss comparisons with other work and the astrophysical implications in the broader context of BBH formation scenarios in Sect. 5.

2. Binary black hole evolution in AGN accretion discs

We briefly introduce an idealised model for the evolution of BBHs embedded in an AGN accretion disc surrounding a central SMBH. We recall the main equations governing the orbital evolution of BBHs throughout the disc-driven and GW-driven regimes. We then emphasise the important role of the eccentricity on the orbital evolution in this channel.

2.1. The coupled “disc+GW”-driven evolution

We consider a BBH system on an elliptic orbit characterised by its semi-major axis a, its orbital eccentricity e, and the total binary mass M_b = m₁ + m₂. The total energy of the BBH is given by $E_b=-\frac{G{M}_b\mu }{2a}$, where $\mu =\frac{m_1{m}_2}{M_b}=\frac{q}{{(1+q)}^2}{M}_b$ is the reduced mass and q = m₂/m₁ is the mass ratio. The orbital angular momentum of the binary is given by $L_b=\mu a^2{\mathrm{\Omega }}_b\sqrt{1-e^2}$, where ${\mathrm{\Omega }}_b=\sqrt{\frac{G{M}_b}{a^3}}$ is the angular frequency of the BBH.

The binary evolves in the gaseous environment of the AGN accretion disc surrounding a central object of mass M_SMBH. We require a stable binary system that is bound to the central mass; that is, when considering the binary system as a whole, orbiting around M_SMBH, the orbital separation a must be smaller than the Roche lobe limit of the central AGN. The stability criterion is given by $\frac{a_{*}}{a}\gtrsim \frac{1+e}{1-e_{*}}{\left(\frac{3M_{\mathrm{SMBH}}}{M_b}\right)}^{1/3}$, where (a_*, e_*) are the semi-major axis and eccentricity of the outer orbit involving the SMBH (Hoang et al. 2018; Fragione et al. 2019). For simplicity, we assume e_* = 0 throughout. The AGN accretion disc is assumed to follow Keplerian rotation, that is, the orbital frequency of the gas at the radius r is $\mathrm{\Omega }=\sqrt{\frac{G{M}_{\mathrm{SMBH}}}{r^3}}$. The disc is characterised by its scale-height H and its gas surface density ${\mathrm{\Sigma }}_g=\frac{f_g{\sigma }^2}{\pi Gr}$, where f_g is the gas mass fraction, and where σ stands for the stellar velocity dispersion. We relate the central mass M_SMBH to the stellar dispersion σ through the empirical “M − σ relation” (Kormendy & Ho 2013):

$$\begin{array}{c}\hfill \frac{M_{\mathrm{SMBH}}}{{10}^9{M}_{\odot }}=(0.{31}_{-0.033}^{+0.037}){\left(\frac{\sigma }{200\mathrm{k}m{s}^{-1}}\right)}^{4.38\pm 0.29}.\end{array}$$(1)

The accretion disc is modelled as a geometrically thin and optically thick disc, characterised by the α-viscosity prescription ν = αc_sH, where 0 <  α <  1 is a dimensionless parameter and c_s = HΩ is the local sound speed (Shakura & Sunyaev 1973). We recall that in such a geometrically thin disc the aspect ratio is always h = H/r ≪ 1. From standard accretion disc theory (Pringle 1981), the viscous torque acting to transfer angular momentum outwards is given by

$$\begin{array}{c}\hfill T_{\mathrm{visc}}(r)=-3\pi \alpha c_s^2(r){\mathrm{\Sigma }}_g(r)r^2.\end{array}$$(2)

Gravitational torques due to the disc–binary interaction may clear an inner cavity in the gas distribution of the AGN disc. This cavity inside of which the BBH resides is surrounded by a circumbinary accretion disc. We assume the tidal torques of the binary to suppress mass accretion into the central cavity, a finding initially noted by Artymowicz & Lubow (1994) as a result of their smoothed particle hydrodynamic (SPH) simulations. We also assume that neither the primary nor the secondary mass have their own minidiscs. We discuss possible effects of the mass inflow on the orbital evolution in Sect. 5.2. Henceforth, we work in the reference frame of the binary system. In this model, no matter is accreting onto the binary and so the binary is injecting angular momentum into the circumbinary disc, and accordingly we require ${\dot{L}}_b<0$. We note that according to this convention, in the binary reference frame, the sign of the viscous torque must be opposite to that of Eq. (2). By conservation of angular momentum, the binary orbital angular momentum is absorbed by the circumbinary disc and transferred through its inner edge r_in, thus $-{\dot{L}}_b\sim T_{\mathrm{visc}}(r_{\mathrm{in}})$.

To relate the viscous angular momentum flux in the inner region of the circumbinary disc to the orbital evolution of the binary, we impose a further assumption to the particular evolution model: we assume that the complex disc–binary interaction can be solely characterised through the angular momentum of the binary. Concretely, we assume that the binary–disc interaction may well be approximated as an adiabatic process. In addition, we assume that torques act on average axisymmetrically onto the disc. This treatment is also employed for similar evolution models outlined in Rafikov (2013) and Hayasaki (2009). Based on these two idealised assumptions, the binary energy dissipation can be related to its change in angular momentum through the orbital frequency, that is, ${\dot{E}}_b={\mathrm{\Omega }}_b{\dot{L}}_b$.

The assumption that non-axisymmetric potential perturbations of the binary system are small around the average binary potential requires that the cavity shape stay circular throughout the evolution. A circular cavity shape may be an acceptable assumption for circular orbits and for orbits with low eccentricities. Non-linear and backcoupling effects are very likely to occur for higher eccentricities and those complicate an analytic treatment considerably. Nevertheless, a circular cavity morphology is adopted in a toy model of Hayasaki (2009), where the orbital eccentricity of binaries is allowed to approach e = 1. Furthermore, gap sizes are commonly reported in terms of the binary orbital semi-major axis up to a multiplicative factor (even for eccentricities up to 0.7, see e.g. Artymowicz & Lubow (1994)). On those grounds, we adopt a circular cavity shape throughout. In more realistic situations, the circumbinary disc can be distorted and become eccentric, as a result of its interaction with the central binary (MacFadyen & Milosavljević 2008; Shi et al. 2012). Eccentric binaries can directly drive the disc eccentricity growth, for example via bar potential (Lubow et al. 2000). As observed in three-dimensional magnetohydrodynamic (3D MHD) simulations by Shi et al. (2012), even near-circular binaries may induce disc eccentricity growth through the impact of gas streams hitting the inner edge of the disc.

Numerical simulations suggest that the inner edge of the circumbinary disc is usually located around twice the semi-major axis, with the cavity size increasing with increasing eccentricity (Artymowicz & Lubow 1994, Hayasaki et al. 2007). In accordance with the above simulations we assume r_in  =  2a(1  +  e), and as a consequence the viscous torque is $T_{\text{visc}}(r_{\text{in}})=12\pi \alpha c_s^2{\mathrm{\Sigma }}_g{a}^2{(1+e)}^2$. We further assume that the values of the circumbinary accretion disc parameters can be set equal to the corresponding local values of the background AGN disc (cf. Baruteau et al. 2011).

The standard evolution equations of an elliptic Keplerian orbit applied to the thus outlined disc-binary interaction give

$$\begin{array}{cc}& \frac{\dot{a}}{a}=-\frac{24\pi \alpha c_s^2{\mathrm{\Sigma }}_g{(1+e)}^2}{\mu {\mathrm{\Omega }}_b},\hfill \end{array}$$(3)

$$\begin{array}{cc}& \frac{\dot{e}e}{1-e^2}=\frac{12\pi \alpha c_s^2{\mathrm{\Sigma }}_g{(1+e)}^2}{\mu {\mathrm{\Omega }}_b}(\frac{1}{\sqrt{1-e^2}}-1).\hfill \end{array}$$(4)

We note that ${\dot{L}}_b<0$ implies Ė_{b < 0}, and therefore in this channel a binary system always evolves toward lower energy. Therefore, even though we have ė ≥ 0, the orbit cannot become unbound. The trend of binary eccentricity growth (ė > 0) was first found by Artymowicz et al. (1991), who used SPH simulations to analyse the effects of the circumbinary disc on the evolution of the binary orbital elements. We also note that it would not be meaningful to assert that the disc dynamics alone can provide a full description of the orbital evolution. In fact, a binary evolution premised purely on Eqs. (3) and (4) seems unphysical, for e → 1 and a → 0 indicates that the orbit enters a regime where gravitational wave emission cannot be neglected.

We recall the gravitational-wave-driven evolution equations of the orbital parameters (Peters 1964):

$$\begin{array}{cc}& \dot{a}=-\frac{64}{5}\frac{G^3}{c^5}\frac{\mu M_b^2}{a^3}\frac{1}{{(1-e^2)}^{7/2}}(1+\frac{73}{24}{e}^2+\frac{37}{96}{e}^4),\hfill \end{array}$$(5)

$$\begin{array}{cc}& \dot{e}=-\frac{304}{15}\frac{G^3}{c^5}\frac{\mu M_b^2}{a^4}\frac{e}{{(1-e^2)}^{5/2}}(1+\frac{121}{304}{e}^2).\hfill \end{array}$$(6)

To investigate the dynamics of the interplay between the GW emission and the disc-binary interaction, we combine the disc-driven (Eqs. (3) and (4)) with the corresponding GW-driven (Eqs. (5) and (6)) equations. The coupled “disc+GW” evolution is thus described by

$$\begin{array}{cc}& \dot{a}={\dot{a}}_{\mathrm{disc}}+{\dot{a}}_{\mathrm{GW}},\hfill \end{array}$$(7)

$$\begin{array}{cc}& \dot{e}={\dot{e}}_{\mathrm{disc}}+{\dot{e}}_{\mathrm{GW}}.\hfill \end{array}$$(8)

This coupled system of differential equations can be integrated numerically; this yields the temporal evolution of the binary semi-major axis and orbital eccentricity. Of particular interest is the merger timescale, which is taken as the point where the numerical solution a(t), or equivalently e(t), approaches the abscissa axis. In the following, we consider the coupled “disc+GW”-driven evolution, focusing on the important role of the orbital eccentricity.

2.2. The role of eccentricity

Figure 1 shows the evolution of the semi-major axis as a function of the orbital eccentricity for different values of the initial eccentricity (e₀ = 0.01−0.9). The black dots indicate the critical radius where GW-emission takes over the disc-driven interaction. We take as fiducial values for the parameters of the accretion disc: viscosity parameter α = 0.1, gas fraction f_g = 0.1, and aspect ratio h = H/r = 0.01. The binary parameters are set to: a₀ = 1 AU, M_b = 50 M_⊙, and q = 1. In this fiducial model, the binary is located at a distance of r = 0.1 pc in an AGN disc surrounding a SMBH of M_SMBH = 10⁷ M_⊙. We note that the stability condition (cf. Sect. 2.1) is satisfied for the fiducial model parameters, and so the binary remains bound to the central object. From Fig. 1 we see that for most values of the semi-major axis, the disc–binary interaction dominates the gravitational-wave emission. Only at the latest stages of inspiral when a is sufficiently low does gravitational emission become important. As for the eccentricity evolution we observe two distinct trends: in the disc-driven regime the orbital eccentricity grows, and in the GW-driven regime the eccentricity decays. The disc–binary interaction shrinks the orbital separation, while at the same time it increases the orbital eccentricity, such as to rapidly bring the BBH into the regime where GW emission takes over.

Fig. 1. Orbital semi-major axis as a function of the orbital eccentricity for a binary system based on the coupled evolution equations (Eqs. (7) and (8)). In this plot we show a(e) for different initial eccentricities (e0 = 0.01−0.9).

At a critical semi-major axis a_c, where ȧ_disc(a_c) = ȧ_GW(a_c), the binary dynamics evolve from being predominantly disc-driven to predominantly GW-driven; the black circles in Fig. 1 indicate these transition points. We observe that the critical semi-major axis a_c increases for increasing initial eccentricities. Thus, for smaller e₀, the disc interaction drives the binary to lower separations as compared to orbits with higher e₀. However, the higher absolute value of the orbital eccentricity at a_c for orbits with higher e₀ results eventually in a merging that is faster than for orbits with lower initial eccentricities. We also note that for a small initial eccentricity (e₀ = 0.01, blue curve in Fig. 1), there is only a small subsequent eccentricity growth in the disc-driven regime.

The combined “disc+GW”-driven evolution reduces the purely gravitational wave-driven merger time by several orders of magnitude. Figure 2 shows, as a function of the initial eccentricity of the orbit, the purely gravitational wave-driven merger times (blue) and those where disc interactions are also taken into account (green). In the case of purely GW-driven evolution, binaries in this particular configuration would not merge within a Hubble time (t_H  ∼  10¹⁰ yr), unless the initial eccentricity is higher than e₀  ≳  0.9. In contrast, in the coupled “disc+GW” evolution, binaries can efficiently merge on timescales of τ  ∼  (10⁶−10⁸) yr for all values of the initial eccentricity. For instance, at e₀  =  0.01, the merger time is τ  ∼   10¹³ years for the GW-only evolution, while it is τ  ∼  2  ×   10⁸ years for the combined “disc+GW” evolution. Hence, if located in the field, such a low-eccentricity binary would not merge within the age of the universe. A gaseous environment would provide a mechanism for such binaries to merge within a plausible timescale. From Fig. 2 we also observe that in both channels the merger timescale decreases as the initial eccentricity increases: at e₀ = 0.99, the merger time is τ ∼ 2 × 10⁷ years for the purely GW-driven evolution, compared to τ ∼ 2 × 10⁵ years for the coupled disc+GW evolution.

Fig. 2. Merger timescales as a function of the initial orbital eccentricity: merger times resulting from GW-only decay (blue curve), and both GW+disc evolution (green curve).

3. Binary black hole merger rate density estimates

In this section, we estimate the BBH merger rate densities that can be expected from the particular evolution channel described in Sect. 2. A natural and common way to parametrise the rate density of BBH mergers in AGN discs is (cf. McKernan et al. 2018)

$$\begin{array}{c}\hfill \mathcal{R}=\frac{n_{\mathrm{AGN}}{N}_{\mathrm{BH}}{f}_d{f}_b}{{\tau }_{\mathrm{avg}}},\end{array}$$(9)

where n_AGN is the number density of AGNs (per Gpc⁻³), N_BH is the number of stellar-mass BHs in the galactic nucleus (within the central ∼pc³), f_d is the fraction of BHs embedded in the AGN accretion disc, f_b is the fraction of such BHs residing in binaries, and τ_avg is the average BBH merger time.

In McKernan et al. (2018) no attempt was made to constrain the actual merger timescale, and binaries were assumed to merge within the AGN disc lifetime (τ_AGN ∼ 1−10 Myr). Based on the model of the previous section, we are able to explicitly compute the merger timescales resulting from the coupled “disc+GW”-driven evolution, taking into account various parameters such as e₀, a₀, f_g, M_SMBH, and so on. We modify the prescription of Eq. (9) and estimate the BBH merger rate density from a weighted average of the merger times of a population of BBHs radially distributed within the AGN accretion disc. Below, we provide a description of the computation steps involved in estimating the BBH merger rate density ℛ.

3.1. Mass distributions

We combine different mass distributions and mass density profiles to derive an expression for the number of black holes with a given mass range that reside within a radial shell of an AGN disc.

3.1.1. Stellar distribution

We start with an initial stellar population. The initial stellar mass distribution seems to follow a multiple power-law function. Here, we adopt the initial stellar mass function of Kroupa (2002)

$$\begin{array}{c}\hfill f_{\mathrm{IMF}}(m)=\{\begin{array}{cc}{k}_0{\left(\frac{m}{M_{\odot }}\right)}^{-0.3}\hfill & \mathrm{if}m<0.08M_{\odot }\hfill \\ k_1{\left(\frac{m}{M_{\odot }}\right)}^{-1.3}\hfill & \mathrm{if}0.08M_{\odot }<m<0.5M_{\odot }\hfill \\ k_2{\left(\frac{m}{M_{\odot }}\right)}^{-2.3}\hfill & \mathrm{if}m>0.5M_{\odot }\hfill \end{array},\end{array}$$(10)

with the normalisation parameters k₀, k₁ = 0.08k₀ and k₂ = 0.04k₀. We assume that the lowest mass of the stars m_s, min is smaller than 0.08 M_⊙ and that the highest mass m_s, max is above 0.5 M_⊙. The total number of stars is $N_{*}={\int }_{m_{\text{s},\text{min}}}^{m_{\text{s},\text{max}}}dm{f}_{\text{IMF}}(m)$ and the total stellar mass is $M_{*}={\int }_{m_{\text{s},\text{min}}}^{m_{\text{s},\text{max}}}dmm{f}_{\text{IMF}}$. For the stellar density profile ρ_*(r), we adopt the “Nuker law” parametrization (Lauer et al. 1995)

$$\begin{array}{c}\hfill {\rho }_{\ast }(r)\propto {\left(\frac{r_b}{r}\right)}^{\gamma }{[1+{\left(\frac{r}{r_b}\right)}^{\alpha }]}^{\frac{\gamma -\beta }{\alpha }},\end{array}$$(11)

and simplify it by using its asymptotic slopes,

$$\begin{array}{c}\hfill {\rho }_{\ast }(r)=\{\begin{array}{cc}{\rho }_0{\left(\frac{r}{r_b}\right)}^{-\gamma }\hfill & \mathrm{if}r\le r_b\hfill \\ {\rho }_0{\left(\frac{r}{r_b}\right)}^{-\beta }\hfill & \mathrm{if}r>r_b\hfill \end{array}.\end{array}$$(12)

The break radius r_b is approximately on the order of the radius of influence of the SMBH (Schödel et al. 2018), characterised as r_b = GM_SMBH/σ².

3.1.2. Black hole distribution

If we assume that every star heavier than m_s, cr ∼ 5 M_⊙ evolves into a BH, the total number of BHs is

$$\begin{array}{c}\hfill N_{\mathrm{BHs}}={\displaystyle {\int }_{m_{\mathrm{s},\mathrm{cr}}}^{m_{\mathrm{s},\max }}\mathrm{d}m{f}_{\mathrm{IMF}}=k{M}_{\ast },}\end{array}$$(13)

where we define

$$\begin{array}{c}\hfill k=\frac{1}{M_{\odot }}\frac{-\frac{0.04}{1.3}[{\left(\frac{m_{\mathrm{s},\max }}{M_{\odot }}\right)}^{-1.3}-{\left(\frac{m_{\mathrm{s},\mathrm{cr}}}{M_{\odot }}\right)}^{-1.3}]}{[0.22-\frac{1}{1.7}{\left(\frac{m_{\mathrm{s},\min }}{M_{\odot }}\right)}^{1.7}-\frac{0.04}{0.3}{\left(\frac{m_{\mathrm{s},\max }}{M_{\odot }}\right)}^{-0.3}]}.\end{array}$$(14)

The total mass of all BHs is

$$\begin{array}{c}\hfill M_{\mathrm{BH}}={\displaystyle {\int }_{m_{\mathrm{s},\mathrm{cr}}}^{m_{\mathrm{s},\max }}\mathrm{d}mm{f}_{\mathrm{IMF}}.}\end{array}$$(15)

The BH mass distribution is denoted as f_BH(m) and is normalised such that ∫dmf_BH(m) = 1. We adopt a power-law mass function for the BH distribution:

$$\begin{array}{c}\hfill f_{\mathrm{BH}}(m)={\xi }_0{m}^{-\kappa }.\end{array}$$(16)

If N(r) denotes the total number of BHs inside some radius r, we take the number dN of BHs lying inside r with masses between m and m + dm to be

$$\begin{array}{cc}& \mathrm{d}N=N(r)f_{\mathrm{BH}}(m)\mathrm{d}m\hfill \\ \hfill & =k{M}_{\ast }(r)f_{\mathrm{BH}}(m)\mathrm{d}m,\hfill \end{array}$$(17)

where M_*(r) is the total stellar mass inside radius r (M_*(∞)≡M_*). We note that for the last equality we assume that the mass distribution of the BHs is independent of the radius at which they reside. We can write dN in the following form

$$\begin{array}{c}\hfill \mathrm{d}N=\mathrm{d}mk{f}_{\mathrm{BH}}(m){\displaystyle {\int }_0^r\mathrm{d}r{\rho }_{\ast }(r)4\pi r^2.}\end{array}$$(18)

The number of BHs with masses between m and m + dm in a radial shell from radius r to r + dr is accordingly

$$\begin{array}{c}\hfill \{\begin{array}{cc}\mathrm{d}mk{f}_{\mathrm{BH}}(m)\frac{4\pi {\rho }_0}{3-\gamma }{r}_b^{\gamma }[{(r+dr)}^{3-\gamma }-r^{3-\gamma }]\hfill & \mathrm{if}r\le r_b\hfill \\ \mathrm{d}mk{f}_{\mathrm{BH}}(m)\frac{4\pi {\rho }_0}{3-\beta }{r}_b^3[{\left(\frac{r+dr}{r_b}\right)}^{3-\beta }-{\left(\frac{r}{r_b}\right)}^{3-\beta }]\hfill & \mathrm{if}r>r_b\hfill \end{array}.\end{array}$$(19)

3.2. Rate per galactic nucleus

To estimate the merger rate per galactic nucleus, we add up the rate contributions coming from different radial shells. We assume that the AGN disc around the SMBH extends from r_min to r_max. We partition the disc into ℐ_r intervals. We use N_j to denote the number of BHs within a radial shell [r_min + jΔr, r_min + (j + 1)Δr], where Δr = (r_max − r_min)/|ℐ_r|. The mean radius of a radial shell is r_j = r_min + (j + 0.5)Δr. The number of BBHs within each radial shell that are embedded in the disc is given by N_BBH, j = N_jf_bf_d/2.

3.2.1. Binary black hole distributions in a radial shell

Within each radial shell labelled “j” we distribute the number of BBHs N_BBH, j for the binary parameter $(m_1,m_2,e_0,a_0)\cong (m_1\equiv \stackrel{\sim }{M},q,e_0,a_0)$ according to some canonical binary parameter distribution functions. We assume for simplicity that the binary parameter distributions are independent from each other. We normalise the distributions in this section with respect to N_BBH, j.

Primary mass distribution: We adopt the initial mass function (16) for the primary mass distribution.

Secondary mass distribution: The secondary mass distribution is determined by the distribution of the mass ratio q, assumed to follow a uniform distribution

$$\begin{array}{c}\hfill f(q)={\vartheta }_0,\end{array}$$(20)

where ϑ₀ is a constant and q_min ≤ q ≤ q_max.

Orbital separation distribution: The distribution of the initial orbital separation a₀ is assumed to be logarithmically flat,

$$\begin{array}{c}\hfill f(a_0)={\chi }_0/a_0,\end{array}$$(21)

between the limits a_0, min and a_0, max. This distribution is biased towards low semi-major axis.

Orbital eccentricity distribution: We choose the distribution of the initial orbital eccentricity e₀ to follow a thermal distribution,

$$\begin{array}{c}\hfill f(e_0)={\varsigma }_0{e}_0,\end{array}$$(22)

between the limits e_0, min and e_0, max. The number of BBHs below a given eccentricity e scale in a thermal distribution as e²; a thermal distribution thus favours high eccentricities.

3.2.2. Calculation of rate within a radial shell

In order to calculate the rate of coalescences in a radial shell j, we partition each binary parameter into bins and distribute N_BBH, j among the different bins. The merger rate is then determined as follows:

The quantities ℐ_x (where x is one of $\stackrel{\sim }{M},q,e_0,$ or a₀) will denote the number of bins for the corresponding binary distributions x. We use N_x, i to denote the number of BBHs residing within the bin i ∈ ℐ_x of the binary distribution x, and use x_i to denote the average of the variable x within bin i. We need to calculate the merger time τ_merger, λ for all possible combinations of $\lambda \in {\mathcal{I}}_{\stackrel{\sim }{M}}\times {\mathcal{I}}_q\times {\mathcal{I}}_{a_0}\times {\mathcal{I}}_{e_0}$. We take the average quantities x_i of the orbital parameter x in bin i as the input values needed for the numerical calculation of the merger time. The total merger rate resulting from a radial shell j is therefore

$$\begin{array}{c}\hfill {\mathcal{R}}_j={\displaystyle \sum _{\lambda }\frac{N_{\lambda }}{N_{\mathrm{BBH},j}^3}\frac{1}{{\tau }_{\mathrm{merger},\lambda }},}\end{array}$$(23)

where the sum goes over all $\lambda \in {\mathcal{I}}_{\stackrel{\sim }{M}}\times {\mathcal{I}}_q\times {\mathcal{I}}_{a_0}\times {\mathcal{I}}_{e_0}$ and where N_λ is the number of BBHs residing with the binning configuration λ.

3.2.3. Merger rate for a galactic nucleus

The merger rate for a galactic nucleus is obtained by summing up the rate (23) for each radial shell j

$$\begin{array}{c}\hfill {\mathcal{R}}_{\mathrm{gal}}={\displaystyle \sum _{j\in {\mathcal{I}}_r}{\mathcal{R}}_j.}\end{array}$$(24)

We operate with the binning ℐ_r = ℐ_e0 = 20 and ${\mathcal{I}}_{\stackrel{\sim }{M}}={\mathcal{I}}_q={\mathcal{I}}_{a_0}=10$, as it is both stable with respect to finer binnings and computationally economical.

3.3. Merger rate density

The merger rate density ℛ is given by

$$\begin{array}{c}\hfill \mathcal{R}={\mathcal{R}}_{\mathrm{gal}}{n}_{\mathrm{AGN}}.\end{array}$$(25)

We note that the number N_* of stars within r_max, the fraction f_b of BH in BBHs, the fraction f_d of BHs embedded in the disc, and the number density n_AGN of AGNs affect the merger rate density ℛ only as scaling factors, and we report rate densities in a form where we set these to a fixed value.

The number density of AGNs is given by n_AGN = f_AGNN_GN, where f_AGN is the fraction of galactic nuclei that are active, and N_GN is the average number density of galactic nuclei. The lower limit on n_AGN would correspond to luminous quasars (f_AGN ∼ 0.01), while the upper limit corresponds to the more common low-luminosity AGNs (f_AGN ∼ 0.3). A significant fraction of BHs in galactic nuclei reside in binary systems, with a typical binary fraction in the range f_b ∼ (0.01−0.2) (McKernan et al. 2018), and possibly up to f_b ∼ (0.6−0.8) for binaries located in the innermost regions of AGN gaseous discs (Secunda et al. 2019). The fraction of BBHs that end up embedded in the AGN disc scales with the disc aspect ratio and is roughly given by f_d ≥ H/r, covering the range f_d ∼ (0.01−0.7) for geometrically thin/thick discs (Ford & McKernan 2019). We therefore adopt the following values for the population parameters (indicating in parenthesis the lower bound and upper bound):

• f_b = 0.1 [0.01, 0.2],

• f_d = 0.01 [0.01, 0.7],

• n_AGN = 4 ×10⁴ Gpc⁻³ [4 × 10⁴ Gpc⁻³, 3 × 10⁶ Gpc⁻³].

In addition we fix the maximal radius as r_max = 1 pc. The inner edge of the accretion disc can in principle extend down to the innermost stable circular orbit, but we do not expect binaries to reside there. We therefore set fiducially r_min to 100 times the Schwarzschild radius of the SMBH. We estimate the number N_* of stars residing within 1 pc as 1 × 10⁶. Furthermore, we use the fiducial values m_s, min = 0.01 M_⊙, m_s, max = 200 M_⊙, and m_s, cr = 5 M_⊙. The number of BHs within 1 pc, obtained from formula (13) is then about 7.5 × 10³, which is close to what is observed in our Galactic centre (Generozov et al. 2018) and to the lower end of the range N_BH ∼ (10⁴−10⁶) estimated in McKernan et al. (2018).

Additionally, throughout we use the fiducial values β = 3.2 and γ = 1.5 for the mass distribution parameters. Also, unless subject to variation, we adopt a viscosity parameter α = 0.1, a gas fraction f_g = 0.1, and a disc aspect ratio H/r = 0.01. We assume M_SMBH = 1 × 10⁷ M_⊙ for our fiducial model, and adopt for the BH mass distribution a Salpeter mass function, that is, we set κ = 2.35. Unless a binary parameter is subject to variation, we use the respective distributions of Sect. 3.2.1 and the following ranges of the orbital parameters

• [m_BH,minm_BH,max] = [5, 50] M_⊙,

• [q_minq_max] = [0.1, 1],

• [a_{0, min}a_{0, max}] = [1, 100] AU,

• [e_{0, min}e_{0, max}] = [0.001, 0.99].

4. Results: GW event rates

In Table 1, we display the resulting BBH merger rate densities ℛ alongside the average and median merger timescales (τ_avg, τ_med) for variations in the binary orbit and AGN disc parameters. We note that the average and median merger timescales generally differ by about two orders of magnitude. This attests to the inadequacy of resorting to average merger timescales for rate estimates. Resorting to an average timescale when determining rates leads to underestimation of BBH mergers with merger times smaller than τ_avg and overestimation of BBHs with merger times higher than τ_avg.

Overall, our predicted BBH-in-AGN merger rates broadly lie in the range ℛ ∼ (0.002−18) Gpc⁻³ yr⁻¹, with a characteristic value of ℛ ∼ 0.2 Gpc⁻³ yr⁻¹ for the fiducial parameter choice described in the previous section. This is to be compared with the corresponding rate of ℛ ∼ 5.5 × 10⁻¹⁰ Gpc⁻³ yr⁻¹ expected for BBH mergers driven solely by GW emission. Thus the merger rate seems to be astrophysically irrelevant in the case of purely GW-driven inspiral, while taking into account a disc–binary interaction can yield astrophysically significant rates. This again highlights the importance of including the disc-driven regime in the coupled “disc+GW” evolution scenario.

In addition, the eccentricity distribution plays an important role in determining the BBH merger rates (see Table 1). For instance, assuming that all binaries are initially on quasi-circular orbits (δ-function at e₀ = 0.001), the resulting merger rate density is ℛ ∼ 3 × 10⁻⁶ Gpc⁻³ yr⁻¹. This astrophysically low value may also be related to the fact that for near-zero initial eccentricities, there is negligible eccentricity growth in the disc-driven regime (cf. Fig. 1). Considering a more realistic distribution of initial eccentricities, such as the thermal or uniform distributions, leads to BBH merger rates that are ≳4 orders of magnitude higher. We assume a thermal distribution for our fiducial model. Assuming a uniform distribution instead results in a merger rate that is a factor of approximately two lower than the fiducial choice. This can be attributed to the fact that a thermal distribution tends to favour high-eccentricity binaries (the mean eccentricity in a thermal distribution is two-thirds, while it is one-third in a uniform distribution).

Similarly, a logarithmically flat distribution in the initial orbital separation is biased toward lower semi-major axis compared to a uniform distribution. As the orbital decay is more efficient at large separations in the disc-driven regime, the resulting merger rate is higher in the latter case (cf. Table 1). On the other hand, variations in the power-law slope of the BH mass distribution seem to have only a modest effect on the merger rates (within a plausible range of κ).

With regards to the AGN parameters, we see that the merger rate is an increasing function of the SMBH mass: ℛ is higher for binaries orbiting more massive central objects. Likewise, we observe that increasing the viscosity parameter α, gas fraction f_g, and disc aspect ratio h lead to higher BBH merger rates. Among the accretion disc parameters, the aspect ratio seems to be the most important factor determining the global ℛ range. In physical terms, the predicted trends can be interpreted as follows: in the disc-driven regime, both the orbital decay and eccentricity growth scale as ∝αf_gh². Therefore, binaries embedded in a gas-rich viscous accretion disc with a large aspect ratio are more likely to merge rapidly, which is reflected in higher GW event rates.

5. Discussion

5.1. The coupled “disc+GW” evolution of BBHs in AGN discs

We consider the evolution of BBHs embedded in AGN accretion discs governed by the coupled “disc+GW”-driven evolution equations (Eqs. (7) and (8)). It is well known that in the GW-driven regime the binary orbital decay is more efficient at small semi-major axis. According to our model, the orbital decay in the disc-driven regime is instead more efficient at large separations. There are two interesting physical trends by which the disc-binary interaction facilitates the merging of the BBH: orbital decay and eccentricity growth. On the one hand, gaseous torques act to shrink the binary separation to small enough radii where GW inspiral can take over. On the other hand, there can be considerable eccentricity growth in the disc-driven regime. Consequently the binary enters the GW-driven regime with high eccentricity, where the merging is fostered due to the steep dependence of the orbital decay on the eccentricity (Eq. (5)). Both effects combine to accelerate the BBH merger in the coupled “disc+GW” evolution scenario, as compared to a purely GW-driven evolution.

As a result, binaries that otherwise would not merge within the age of the universe (if located in galactic fields), can be efficiently driven to undergo a merger in the gaseous environment of AGN accretion discs. This observation is not new and has been suggested by various authors (Cuadra et al. 2009; McKernan et al. 2018). A novelty of this work lies in the quantification of this statement by explicitly computing the resulting BBH merger rates based on the weighted averages of binary populations embedded within the AGN accretion disc. From our calculations we expect typical GW event rates in the range ℛ ∼ (0.002−18) Gpc⁻³ yr⁻¹. The current LIGO-Virgo estimate for the BBH merger rate density lies in the range ℛ = (9.7−101) Gpc⁻³ yr⁻¹. Our predicted merger rates are clearly lower than the observationally inferred ones. Importantly, not all black hole binaries originate from this particular AGN formation channel. Nevertheless, this peculiar path may significantly contribute to the total merger rate. In fact, our rate estimates should be viewed as lower limits, since we adopt the lower bounds for most of the physical parameters (e.g. n_AGN, f_d, see Sect. 3). Therefore the actual BBH merger rates resulting from gas-assisted orbital evolution could be easily one or more orders of magnitudes higher.

For comparison, a merger rate of ℛ ∼ 3 Gpc⁻³ yr⁻¹ is expected for BBHs embedded in a self-gravitating disc of AGNs (Stone et al. 2017); while a similar rate of ℛ ∼ 1.2 Gpc⁻³ yr⁻¹ is predicted for BBHs trapped in the innermost regions of the AGN disc (Bartos et al. 2017). Nevertheless, large uncertainties are always involved, and taking into account the entirety of the allowed range of physical parameters leads to a span of ℛ = (10⁻³−10⁴) Gpc⁻³ yr⁻¹ for BBH mergers in AGN discs (McKernan et al. 2018). Recently, Tagawa et al. (2019) combined N-body simulations with a semi-analytical model to investigate how binaries form and merge in AGN discs. The setup of these latter authors incorporates the formation and disruption of binaries, gas dynamical friction, torques from circumbinary discs, migration in the AGN disc, and several different types of stellar interaction. By this means they estimate a merger fraction per BH in an AGN lifetime, and thereby find a volumetric BBH merger rate of ℛ = (0.02−60) Gpc⁻³yr⁻¹. In contrast to this latter study, we do not assume zero orbital eccentricity in our evolution channel, although indeed our treatment rests on idealised assumptions (see Sects. 2.1 and 5.2). From our toy model, we explicitly compute the weighted average of the merger timescales of BBH populations radially distributed in the AGN accretion disc, and calculate the resulting GW event rates. Overall, our predicted merger rates are broadly comparable with similar scenarios of BBHs in AGN discs, although our fiducial rate lies on the lower side (but we also note that a favourable combination of the accretion disc parameters could easily increase the fiducial rate, bringing it more in line with other estimates).

One of the purposes of this study is to explore how the rate density ℛ varies as a function of the underlying physical parameters, that is, the binary and accretion disc configurations. For what concerns the accretion disc parameters, we observe that higher viscosities, larger gas fractions, and larger aspect ratios can lead to enhanced merger rates. According to our parameter space study, the key factor is the disc aspect ratio h = H/r, with geometrically thicker discs (h ∼ 0.1, that is not razor-thin) allowing larger BBH merger rates. As already noted, the merger rate is also a slowly increasing function of the central super-massive black hole mass M_SMBH. Summarising, in our simplified model, BBH mergers are most effective when the binaries are embedded in a high-viscosity, gas-rich, and relatively thick accretion disc surrounding a massive central object.

A distinctive feature of the BBH evolution in AGN accretion discs in this channel is the eccentricity growth in the disc-driven regime. The orbital eccentricity plays a major role in this scenario: it is responsible for considerably reducing the merger timescales, thereby increasing the corresponding merger rate densities. Importantly, including eccentric binaries leads to GW event rates enhanced by several orders of magnitude with respect to circular binaries (Sect. 4). We note that a finite non-zero value of the initial eccentricity is also required in order to get efficient eccentricity growth in the disc regime (Fig. 1). Such initial seed eccentricities could be acquired through different stellar dynamical processes operating in galactic nuclei, such as the Kozai-Lidov mechanism, and other triple systems formed via binary–binary or binary–single interactions (e.g. Rasskazov & Kocsis 2019, and references therein).

Another unique signature of a BBH-in-AGN channel is the potential association with electromagnetic (EM) counterparts (Stone et al. 2017; Bartos et al. 2017). In general, BBH mergers are not expected to produce EM counterparts, unlike mergers involving neutron stars. However, in the dense gaseous environment of AGN accretion discs, binary black holes may accrete significant amounts of gas at high rates, possibly exceeding the Eddington limit. Such super-Eddington accretion episodes may lead to transient EM counterparts, through the development of ultrafast outflows and/or relativistic jets, radiating over a broad range of wavelengths (e.g. Murase et al. 2016). The associated high-energy emission, such as X-rays or gamma-rays, could potentially be detected by current or upcoming space-borne instruments (Bartos et al. 2017).

5.2. Caveats and outlook

Finally, because many simplifying assumptions are adopted in this work, our results should be considered as tentative. As already mentioned in Sect. 2.1, our disc–binary interaction is mainly based on angular momentum conservation. A more refined and correct treatment of the problem requires taking into account resonant interaction mechanisms, such as co-rotation and Lindblad resonances, between the binary and its surrounding disc. The outer Lindblad resonances usually tend to increase the binary eccentricity, whereas the inner Lindblad and corotational resonances tend to damp the eccentricity growth (Artymowicz et al. 1991). At high eccentricities, several competing resonances can be simultaneously present. Around e ≳ 0.5, the opposite effects of such resonances may partially cancel out, considerably reducing the eccentricity growth rate (Lubow et al. 2000). Moreover, the eccentricity pumping tends to decline, and eventually saturates for higher eccentricities (Dermine et al. 2013). Due to the potential overestimation of the eccentricity growth rate, our reported merger timescales and GW rate densities may be too high overall. Further work is required to see how such non-linear effects can be properly incorporated into our toy model (one possibility is to use the resonance formalism of Goldreich & Tremaine (1980), but this formalism is inapplicable unless e ≪ 1).

Furthermore, we assume that the binary torque is efficient at suppressing the mass inflow from the disc. That is, we completely neglect gas inflows into the cavity within which the binary resides, and the subsequent accretion onto the individual black holes. Based on this, in our model the binary is only losing angular momentum to the disc and is thus invariably migrating inwards. More generally, the binary loses angular momentum to the disc via gravitational torques, while it gains angular momentum through accretion. The binary orbital decay rate is then roughly determined by the balance between these two opposite processes. Several numerical studies suggest that the trend of inward migration holds even when including the effects of mass transfer, because the angular momentum gain by accretion is usually not enough to offset the angular momentum loss to the disc (MacFadyen & Milosavljević 2008; Cuadra et al. 2009; Shi et al. 2012). Three-dimensional MHD simulations of circumbinary discs indicate the development of two narrow streams flowing from the inner edge of the disc towards the binary, supplying substantial mass accretion. The time-averaged accretion rate onto the binary is found to be comparable to that on a single central mass (Shi & Krolik 2015). The large accretion rates imply a smaller orbital shrinkage rate, although the net result is ȧ/a < 0 (Shi et al. 2012).

In contrast to the above conclusions, recent 2D and 3D viscous hydrodynamical simulations of circumbinary accretion indicate that accreting binaries consistently gain angular momentum from the disc (Moody et al. 2019; Muñoz et al. 2019, 2020). This implies that the binaries expand as they accrete, with positive ⟨ȧ⟩ > 0 (for moderate mass ratios). On the other hand, Ragusa et al. (2016) argue that in the case of thin discs with small aspect ratios (h ≲ 10⁻², typical of AGN discs) the accretion rate is suppressed in their 3D SPH simulations (note that most numerical studies rather assume h ∼ 0.1). A similar conclusion, namely that the accretion rate onto the binary is significantly reduced for thin discs, is also reached in 2D hydrodynamic simulations (Terquem & Papaloizou 2017). Such a dependence on the disc aspect ratio was previously noted by Artymowicz & Lubow (1996), who suggested that the development of efficient accreting gas streams requires warm and thick discs (with aspect ratio h ≳ 0.05). The overall binary evolution is also influenced by the uncertain evolution of the eccentricity, which may assist the orbital decay (Duffell et al. 2019), but further investigations are required to obtain a more complete picture.

In addition, we do not consider the stellar processes occurring on larger scales. These are also likely at the origin of the non-negligible initial eccentricities, which are relevant for the subsequent binary evolution. The inclusion of such additional mechanisms should lead to more physically motivated initial conditions. Ideally, one would follow the evolution of the BBHs all the way from galactic scales down into the GW-driven regime. Indeed, BBHs are ultimately driven to merger by a combination of at least three physical processes operating on different physical scales: three-body stellar scatterings at large radii, gaseous torques in AGN accretion discs at intermediate radii, and finally GW emission at small radii.

Acknowledgments

WI and PJ acknowledge support from the University of Zurich. ST is supported by Forschungskredit Nr. FK-19-114. MH acknowledges support from Swiss National Science Foundation (SNSF) grant Nr. IZCOZ0-177057.

References

All Tables

All Figures

	Fig. 1. Orbital semi-major axis as a function of the orbital eccentricity for a binary system based on the coupled evolution equations (Eqs. (7) and (8)). In this plot we show a(e) for different initial eccentricities (e0 = 0.01−0.9).
In the text

	Fig. 2. Merger timescales as a function of the initial orbital eccentricity: merger times resulting from GW-only decay (blue curve), and both GW+disc evolution (green curve).
In the text