Issue 
A&A
Volume 638, June 2020



Article Number  A119  
Number of page(s)  8  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/202037681  
Published online  23 June 2020 
Binary black hole mergers in AGN accretion discs: gravitational wave rate density estimates
^{1}
Department of Physics, ETH Zürich, Zürich, Switzerland
email: matthias.groebner@student.ethz.ch
^{2}
PhysikInstitut, Universität Zürich, Zürich, Switzerland
email: wako.ishibashi@physik.uzh.ch
Received:
7
February
2020
Accepted:
4
May
2020
The majority of gravitational wave (GW) events detected so far by LIGO/Virgo originate from binary black hole (BBH) mergers. Among the different binary evolution paths, the merger of BBHs in accretion discs of active galactic nuclei (AGNs) is a possible source of GW detections. We consider an idealised analytical model of the orbital evolution of BBHs embedded in an AGN accretion disc. In this framework, the disc–binary interaction increases the orbital eccentricity and decreases the orbital separation, driving the BBH into a regime where GW emission eventually leads to coalescence. We compute the resulting GW merger rate density from this channel based on a weighted average of the merger timescales of a population of BBHs radially distributed within the AGN accretion disc. The predicted merger rates broadly lie in the range ℛ ∼ (0.002−18) Gpc^{−3} yr^{−1}. We analyse the dependence of the merger rate density on both the accretion disc and binary orbital parameters, emphasising the important role of the orbital eccentricity. We discuss the astrophysical implications of this particular BBHinAGN formation channel in the broader context of binary evolution scenarios.
Key words: gravitational waves / black hole physics / accretion, accretion disks / methods: analytical
© 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 stellarmass 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^{−3} yr^{−1} (Abbott et al. 2019a), with an updated estimate of (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 nearzero eccentricity close to coalescence, while nonnegligible eccentricities are expected for binaries formed through dynamical interactions in dense star clusters (Breivik et al. 2016; Samsing & RamirezRuiz 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 stellarmass black holes has been recently uncovered in the central parsec of our own Galaxy (Hailey et al. 2018), suggesting that ∼10^{4} BHs may be accumulated in the Galactic centre (Generozov et al. 2018). Galactic nuclei are also known to harbour supermassive black holes (SMBHs) with typical masses of M_{SMBH} ∼ (10^{6}−10^{9}) 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^{−3}−10^{4}) Gpc^{−3} yr^{−1}. A recent study by Tagawa et al. (2019) combined Nbody simulations with semianalytical 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^{−3} yr^{−1}.
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 discdriven evolution equations to the corresponding equations of GWemission. 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 discdriven and GWdriven 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 discdriven and GWdriven 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 semimajor axis a, its orbital eccentricity e, and the total binary mass M_{b} = m_{1} + m_{2}. The total energy of the BBH is given by , where is the reduced mass and q = m_{2}/m_{1} is the mass ratio. The orbital angular momentum of the binary is given by , where 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 , where (a_{*}, e_{*}) are the semimajor 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 . The disc is characterised by its scaleheight H and its gas surface density , 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):
The accretion disc is modelled as a geometrically thin and optically thick disc, characterised by the αviscosity prescription ν = αc_{s}H, 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
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 . 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 .
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, .
The assumption that nonaxisymmetric 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. Nonlinear 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 semimajor 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 threedimensional magnetohydrodynamic (3D MHD) simulations by Shi et al. (2012), even nearcircular 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 semimajor 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 . 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 discbinary interaction give
We note that 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 gravitationalwavedriven evolution equations of the orbital parameters (Peters 1964):
To investigate the dynamics of the interplay between the GW emission and the discbinary interaction, we combine the discdriven (Eqs. (3) and (4)) with the corresponding GWdriven (Eqs. (5) and (6)) equations. The coupled “disc+GW” evolution is thus described by
This coupled system of differential equations can be integrated numerically; this yields the temporal evolution of the binary semimajor 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 semimajor axis as a function of the orbital eccentricity for different values of the initial eccentricity (e_{0} = 0.01−0.9). The black dots indicate the critical radius where GWemission takes over the discdriven 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_{0} = 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^{7} 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 semimajor axis, the disc–binary interaction dominates the gravitationalwave 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 discdriven regime the orbital eccentricity grows, and in the GWdriven 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 semimajor 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 (e_{0} = 0.01−0.9). 

Open with DEXTER 
At a critical semimajor axis a_{c}, where ȧ_{disc}(a_{c}) = ȧ_{GW}(a_{c}), the binary dynamics evolve from being predominantly discdriven to predominantly GWdriven; the black circles in Fig. 1 indicate these transition points. We observe that the critical semimajor axis a_{c} increases for increasing initial eccentricities. Thus, for smaller e_{0}, the disc interaction drives the binary to lower separations as compared to orbits with higher e_{0}. However, the higher absolute value of the orbital eccentricity at a_{c} for orbits with higher e_{0} 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} = 0.01, blue curve in Fig. 1), there is only a small subsequent eccentricity growth in the discdriven regime.
The combined “disc+GW”driven evolution reduces the purely gravitational wavedriven merger time by several orders of magnitude. Figure 2 shows, as a function of the initial eccentricity of the orbit, the purely gravitational wavedriven merger times (blue) and those where disc interactions are also taken into account (green). In the case of purely GWdriven evolution, binaries in this particular configuration would not merge within a Hubble time (t_{H} ∼ 10^{10} yr), unless the initial eccentricity is higher than e_{0} ≳ 0.9. In contrast, in the coupled “disc+GW” evolution, binaries can efficiently merge on timescales of τ ∼ (10^{6}−10^{8}) yr for all values of the initial eccentricity. For instance, at e_{0} = 0.01, the merger time is τ ∼ 10^{13} years for the GWonly evolution, while it is τ ∼ 2 × 10^{8} years for the combined “disc+GW” evolution. Hence, if located in the field, such a loweccentricity 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} = 0.99, the merger time is τ ∼ 2 × 10^{7} years for the purely GWdriven evolution, compared to τ ∼ 2 × 10^{5} years for the coupled disc+GW evolution.
Fig. 2. Merger timescales as a function of the initial orbital eccentricity: merger times resulting from GWonly decay (blue curve), and both GW+disc evolution (green curve). 

Open with DEXTER 
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)
where n_{AGN} is the number density of AGNs (per Gpc^{−3}), N_{BH} is the number of stellarmass BHs in the galactic nucleus (within the central ∼pc^{3}), 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_{0}, a_{0}, 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 powerlaw function. Here, we adopt the initial stellar mass function of Kroupa (2002)
with the normalisation parameters k_{0}, k_{1} = 0.08k_{0} and k_{2} = 0.04k_{0}. 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 and the total stellar mass is . For the stellar density profile ρ_{*}(r), we adopt the “Nuker law” parametrization (Lauer et al. 1995)
and simplify it by using its asymptotic slopes,
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}/σ^{2}.
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
where we define
The total mass of all BHs is
The BH mass distribution is denoted as f_{BH}(m) and is normalised such that ∫dmf_{BH}(m) = 1. We adopt a powerlaw mass function for the BH distribution:
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
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
The number of BHs with masses between m and m + dm in a radial shell from radius r to r + dr is accordingly
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_{j}f_{b}f_{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 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
where ϑ_{0} is a constant and q_{min} ≤ q ≤ q_{max}.
Orbital separation distribution: The distribution of the initial orbital separation a_{0} is assumed to be logarithmically flat,
between the limits a_{0, min} and a_{0, max}. This distribution is biased towards low semimajor axis.
Orbital eccentricity distribution: We choose the distribution of the initial orbital eccentricity e_{0} to follow a thermal distribution,
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^{2}; 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 or a_{0}) 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 . 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
where the sum goes over all 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
We operate with the binning ℐ_{r} = ℐ_{e0} = 20 and , 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
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_{AGN}N_{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 lowluminosity 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^{4} Gpc^{−3} [4 × 10^{4} Gpc^{−3}, 3 × 10^{6} Gpc^{−3}].
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^{6}. 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^{3}, 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^{4}−10^{6}) 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^{7} 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,min}m_{BH,max}] = [5, 50] M_{⊙},

[q_{min}q_{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}.
Merger rate densities and merger timescales for variations in the binary orbit parameters and AGN disc parameters.
Overall, our predicted BBHinAGN merger rates broadly lie in the range ℛ ∼ (0.002−18) Gpc^{−3} yr^{−1}, with a characteristic value of ℛ ∼ 0.2 Gpc^{−3} yr^{−1} for the fiducial parameter choice described in the previous section. This is to be compared with the corresponding rate of ℛ ∼ 5.5 × 10^{−10} Gpc^{−3} yr^{−1} expected for BBH mergers driven solely by GW emission. Thus the merger rate seems to be astrophysically irrelevant in the case of purely GWdriven inspiral, while taking into account a disc–binary interaction can yield astrophysically significant rates. This again highlights the importance of including the discdriven 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 quasicircular orbits (δfunction at e_{0} = 0.001), the resulting merger rate density is ℛ ∼ 3 × 10^{−6} Gpc^{−3} yr^{−1}. This astrophysically low value may also be related to the fact that for nearzero initial eccentricities, there is negligible eccentricity growth in the discdriven 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 higheccentricity binaries (the mean eccentricity in a thermal distribution is twothirds, while it is onethird in a uniform distribution).
Similarly, a logarithmically flat distribution in the initial orbital separation is biased toward lower semimajor axis compared to a uniform distribution. As the orbital decay is more efficient at large separations in the discdriven regime, the resulting merger rate is higher in the latter case (cf. Table 1). On the other hand, variations in the powerlaw 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 discdriven regime, both the orbital decay and eccentricity growth scale as ∝αf_{g}h^{2}. Therefore, binaries embedded in a gasrich 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 GWdriven regime the binary orbital decay is more efficient at small semimajor axis. According to our model, the orbital decay in the discdriven regime is instead more efficient at large separations. There are two interesting physical trends by which the discbinary 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 discdriven regime. Consequently the binary enters the GWdriven 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 GWdriven 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^{−3} yr^{−1}. The current LIGOVirgo estimate for the BBH merger rate density lies in the range ℛ = (9.7−101) Gpc^{−3} yr^{−1}. 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 gasassisted orbital evolution could be easily one or more orders of magnitudes higher.
For comparison, a merger rate of ℛ ∼ 3 Gpc^{−3} yr^{−1} is expected for BBHs embedded in a selfgravitating disc of AGNs (Stone et al. 2017); while a similar rate of ℛ ∼ 1.2 Gpc^{−3} yr^{−1} 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^{−3}−10^{4}) Gpc^{−3} yr^{−1} for BBH mergers in AGN discs (McKernan et al. 2018). Recently, Tagawa et al. (2019) combined Nbody simulations with a semianalytical 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^{−3}yr^{−1}. 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 razorthin) allowing larger BBH merger rates. As already noted, the merger rate is also a slowly increasing function of the central supermassive black hole mass M_{SMBH}. Summarising, in our simplified model, BBH mergers are most effective when the binaries are embedded in a highviscosity, gasrich, 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 discdriven 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 nonzero 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 KozaiLidov 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 BBHinAGN 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 superEddington 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 highenergy emission, such as Xrays or gammarays, could potentially be detected by current or upcoming spaceborne 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 corotation 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 nonlinear 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). Threedimensional 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 timeaveraged 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^{−2}, 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 nonnegligible 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 GWdriven regime. Indeed, BBHs are ultimately driven to merger by a combination of at least three physical processes operating on different physical scales: threebody 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. FK19114. MH acknowledges support from Swiss National Science Foundation (SNSF) grant Nr. IZCOZ0177057.
References
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 061102 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019a, Phys. Rev. X, 9, 031040 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019b, ApJ, 882, L24 [NASA ADS] [CrossRef] [Google Scholar]
 Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651 [NASA ADS] [CrossRef] [Google Scholar]
 Artymowicz, P., & Lubow, S. H. 1996, ApJ, 467, L77 [NASA ADS] [CrossRef] [Google Scholar]
 Artymowicz, P., Clarke, C. J., Lubow, S. H., & Pringle, J. E. 1991, ApJ, 370, L35 [NASA ADS] [CrossRef] [Google Scholar]
 Bartos, I., Kocsis, B., Haiman, Z., & Márka, S. 2017, ApJ, 835, 165 [CrossRef] [Google Scholar]
 Baruteau, C., Cuadra, J., & Lin, D. N. C. 2011, ApJ, 726, 28 [CrossRef] [Google Scholar]
 Breivik, K., Rodriguez, C. L., Larson, S. L., Kalogera, V., & Rasio, F. A. 2016, ApJ, 830, L18 [CrossRef] [Google Scholar]
 Cuadra, J., Armitage, P. J., Alexander, R. D., & Begelman, M. C. 2009, MNRAS, 393, 1423 [NASA ADS] [CrossRef] [Google Scholar]
 Dermine, T., Izzard, R., Jorissen, A., & Winckel, H. 2013, A&A, 551, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Duffell, P. C., D’Orazio, D., Derdzinski, A., et al. 2019, ArXiv eprints [arXiv:1911.05506] [Google Scholar]
 Ford, K. E. S., & McKernan, B. 2019, MNRAS, 490, L42 [CrossRef] [Google Scholar]
 Fragione, G., Grishin, E., Leigh, N. W. C., Perets, H. B., & Perna, R. 2019, MNRAS, 488, 47 [CrossRef] [Google Scholar]
 Generozov, A., Stone, N. C., Metzger, B. D., & Ostriker, J. P. 2018, MNRAS, 478, 4030 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425 [Google Scholar]
 Hailey, C. J., Mori, K., Bauer, F. E., et al. 2018, Nature, 556, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Hayasaki, K. 2009, PASJ, 61, 65 [NASA ADS] [Google Scholar]
 Hayasaki, K., Mineshige, S., & Sudou, H. 2007, PASJ, 59, 427 [NASA ADS] [Google Scholar]
 Hoang, B.M., Naoz, S., Kocsis, B., Rasio, F. A., & Dosopoulou, F. 2018, ApJ, 856, 140 [NASA ADS] [CrossRef] [Google Scholar]
 Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511 [NASA ADS] [CrossRef] [Google Scholar]
 Kroupa, P. 2002, Science, 295, 82 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Lauer, T. R., Ajhar, E. A., Byun, Y.I., et al. 1995, ApJ, 110, 2622 [NASA ADS] [CrossRef] [Google Scholar]
 Lubow, S. H., & Artymowicz, P. 2000, in Protostars and Planets IV, eds. V. Mannings, A. P. Boss, & S. S. Russell, 731 [Google Scholar]
 MacFadyen, A. I., & Milosavljević, M. 2008, ApJ, 672, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Mandel, I., & Farmer, A. 2018, ArXiv eprints [arXiv:1806.05820] [Google Scholar]
 Mapelli, M. 2018, ArXiv eprints [arXiv:1809.09130] [Google Scholar]
 McKernan, B., Ford, K. E. S., Bellovary, J., et al. 2018, ApJ, 866, 66 [CrossRef] [Google Scholar]
 Moody, M. S. L., Shi, J.M., & Stone, J. M. 2019, ApJ, 875, 66 [NASA ADS] [CrossRef] [Google Scholar]
 Muñoz, D. J., Miranda, R., & Lai, D. 2019, ApJ, 871, 84 [NASA ADS] [CrossRef] [Google Scholar]
 Muñoz, D. J., Lai, D., Kratter, K., & Miranda, R. 2020, ApJ, 889, 114 [CrossRef] [Google Scholar]
 Murase, K., Kashiyama, K., Mészáros, P., Shoemaker, I., & Senno, N. 2016, ApJ, 822, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Peters, P. C. 1964, Phys. Rev., 136, 1224 [NASA ADS] [CrossRef] [Google Scholar]
 Pringle, J. E. 1981, ARA&A, 19, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Rafikov, R. R. 2013, ApJ, 774, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Ragusa, E., Lodato, G., & Price, D. J. 2016, MNRAS, 460, 1243 [CrossRef] [Google Scholar]
 Rasskazov, A., & Kocsis, B. 2019, ApJ, 881, 20 [CrossRef] [Google Scholar]
 Samsing, J., & RamirezRuiz, E. 2017, ApJ, 840, L14 [CrossRef] [Google Scholar]
 Schödel, R., GallegoCano, E., Dong, H., et al. 2018, A&A, 609, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Secunda, A., Bellovary, J., Mac Low, M.M., et al. 2019, ApJ, 878, 85 [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33 [NASA ADS] [Google Scholar]
 Shi, J.M., & Krolik, J. H. 2015, ApJ, 807, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Shi, J.M., Krolik, J. H., Lubow, S. H., & Hawley, J. F. 2012, ApJ, 749, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, N. C., Metzger, B. D., & Haiman, Z. 2017, MNRAS, 464, 946 [NASA ADS] [CrossRef] [Google Scholar]
 Tagawa, H., Haiman, Z., & Kocsis, B. 2019, ApJ, submitted [arXiv:1912.08218] [Google Scholar]
 Terquem, C., & Papaloizou, J. C. B. 2017, MNRAS, 464, 2429 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Merger rate densities and merger timescales for variations in the binary orbit parameters and AGN disc parameters.
All Figures
Fig. 1. Orbital semimajor 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 (e_{0} = 0.01−0.9). 

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

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.