Constraining stellar mass black hole mergers in AGN disks detectable with LIGO
Black hole mergers detectable with LIGO can occur in active galactic nucleus (AGN) disks. Here we parameterize the merger rates, the mass spectrum and the spin spectrum of black holes (BH) in AGN disks. The predicted merger rate spans , so upper limits from LIGO () already constrain it. The predicted mass spectrum has the form of a broken power-law consisting of a pre-existing BH powerlaw mass spectrum and a harder powerlaw mass spectrum resulting from mergers. The predicted spin spectrum is multi-peaked with the evolution of retrograde spin BH in the gas disk playing a key role. We outline the large uncertainties in each of these LIGO observables for this channel and we discuss ways in which they can be constrained in the future.
The gravitational wave (GW) events detected by the Advanced Laser Interferometer Gravitational-Wave Observatory (LIGO) correspond to the merger of stellar mass black holes (BH) considerably more massive than those observed in our own Galaxy. The upper end of the range of BH merger rates derived from LIGO observations of 212 Gpc yr (Abbott et al., 2016b) requires consideration of locations where BH mergers can occur faster than expected from GW emission alone. Among the first few LIGO detections are possible low value spin or misaligned spins, which may be problematic for models of binary evolution (O’Shaugnessey, Gerosa & Wysocki, 2017). While BHs with larger than expected masses can occur naturally in the field (Belczynski et al., 2010; deMink & Mandel, 2016), they are more likely to form in regions with concentrations of BHs, such as galactic nuclear star clusters (Hopman & Alexander, 2006; O’Leary et al., 2009; Antonini & Rasio, 2016; Rodriguez et al., 2016). Massive gas disks in active galactic nuclei (AGN) provide natural locations for gas accretion and repeated mergers because the gas disk can drive migration of BH towards migration traps, reduce the inclination of intersecting orbits, enable binary formation, and harden existing binaries. Together, these effects can result in rapid increase in the mass of embedded BHs, potentially to observed values (e.g. McKernan et al., 2012, 2014; Bellovary et al., 2016; Bartos et al., 2017; Stone et al., 2017).
In this paper we parameterize the expected merger rate, and the mass and spin distributions from this channel for comparison with the LIGO observations, and we discuss how observations and simulations can constrain these predictions.
2 Model Outline
Galactic nuclei likely contain some of the densest concentrations of BHs in the Universe (e.g. Morris, 1993; Miralda-Escudé & Gould, 2000, and references therein), so it is natural to look for BH mergers in galactic nuclei (O’Leary et al., 2009; McKernan et al., 2012; Antonini et al., 2014). While BH binary mergers can occur at modestly enhanced rates (compared to the field) in nuclear star clusters just from dynamical binary hardening (Antonini & Rasio, 2016; Rodriguez et al., 2016), or capture from single-single (O’Leary et al., 2009) and binary-single encounters (Samsing et al., 2014), a dense nuclear disk of gas can greatly accelerate the rate of BHB formation and merger (McKernan et al., 2012, 2014).
The simplest picture of this LIGO channel begins with a spherical distribution of BH, stars and other stellar remnants orbiting in the central of a galactic nuclei around a supermassive black hole (SMBH). Next, around the SMBH, we add a massive gas disk, which can be geometrically thin or thick. A fraction of the initial number of BH in the nucleus , will have orbits coincident with the disk and approximately half of these orbits should be retrograde compared to the disk gas. Yet another fraction of the population intersect the disk on their orbits and are ground down into the plane of the disk within the AGN disk lifetime (). Thus an overall fraction of nuclear BH end up embedded in the disk, and quickly have their orbits damped and circularized by gas drag (e.g. McKernan et al., 2012). The net torques from disk gas causes BH to migrate within the disk and encounter each other at low relative velocities (McKernan et al., 2012; Bellovary et al., 2016). BH binaries that form in the disk are expected to merge efficiently due to gas torques (e.g. Haiman et al., 2009; Kocsis et al., 2011; McKernan et al., 2011; Stahler, 2010; Baruteau et al., 2011; McKernan et al., 2012). BH mergers may preferentially occur in convergence zones containing migration traps Bellovary et al. (2016) which occur in semi-realistic models of AGN disks (Sirko & Goodman, 2003; Thompson et al., 2005). Multiple objects trapped in such orbits collide efficiently rather than being ejected (Horn et al. (2012); Secunda, Bellovary et al. (2018) in prep.). In this paper, we examine what constraints can be put on the merger rate and the BH spin and mass distributions for this AGN channel.
3 Rate of black hole binary mergers in AGN disks
We parameterize the rate of BH-BH mergers in AGN disks simply as:
where is the average number density of galactic nuclei in the Universe, is the fraction of galactic nuclei that have active AGNs which last for time , is the fraction of nuclear BH that end up in the disk, is the fraction of BH in BH-BH binaries in the disk, and represents the fractional change in number of BH in the central region () over a full AGN duty cycle 111If then is approximately conserved between AGN episodes. If (grows) shrinks between AGN phases due to the net effect of mergers, infall of new BH, stellar evolution etc.. can be parameterized as:
However, if we want to constrain the constributions of this channel to LIGO observations, it is much more useful to show the allowed range of and the range of each of the contributing factors from eqn. (1), which we list in Table 1.
The lower limit corresponds to galaxies with stellar mass greater than or equal to that of the Milky Way (Baldry et al., 2012) as measured from Schechter function fits to galaxy luminosity functions (e.g. Cole et al., 2001). The upper limit corresponds to dwarf galaxies with stellar mass (Baldry et al., 2012), which includes all locally observed SMBH () inferred from studies of galaxies and dwarf galaxies (Reines & Volonteri, 2015).
Also in Table 1, corresponds to the number of BH allowed of Sgr A* according to the distribution of the S-star orbits (Antonini et al., 2014), whereas pc seems to be the maximal density allowed by simulations (Antonini et al., 2014).
The lower limit to assumes only quasar disks are efficient BH merger sites and assumes all LINER galactic nuclei (Ho, 2008) consist of advection dominated accretion flows (ADAFs) with high accretion rate (Paczynski & Witta, 1980; Narayan & Yi, 1995; Lasota et al., 2016), capable of driving BH mergers.
The binary fraction of BH has been estimated to be as high as (Antonini et al., 2014)), but dynamically hot environments such as star clusters, could actually yield very low binary fractions over time in the absence of gas (Miller & Davies, 2012; Leigh et al., 2016) due to the large number of ’ionizing’ interactions, so we choose in Table 1.
Reasonable estimates of span 0.1-100Myr (Haehnelt & Rees, 1993; King & Nixon, 2015; Schawinski et al., 2015). will be highest if AGN episodes are short-lived but frequently repeated and efficient at BH mergerse. These circumstances ensure that there are multiple opportunities for BH in a galactic nucleus to encounter each other at low relative velocity and merge in a disk.
Note. – Range of parameters in Eqn. (1) and range of merger rate (see text). from Baldry et al. (2012). from Miralda-Escudé & Gould (2000); Antonini et al. (2014). for Seyfert AGN (Ho, 2008). with all LINERs and other low luminosity AGNs. . comes from h/R, the disk aspect ratio. h/R 0.01–0.1 (Sirko & Goodman, 2003). – (Thompson et al., 2005). –0.7 in super-Eddington ADAFs (Lasota et al., 2016). depends on , and .
From Table 1, the allowed range from Eqn. (1) is –. The upper bound to the LIGO BH binary merger rate of already rules out upper limits to most parameters in Table 1 222The LIGO rate upper bound places a lower limit on , since a small value of suggests most BH in AGN are consumed in mergers and would imply a much greater than observed and allows actual astrophysical limits to be placed on models of AGN disks by LIGO BH merger detections. Future observational constraints and simulation results will, however, be required to figure out which upper limits are ruled out by LIGO. For example, the upper limit to could be reduced by contrasting activity rates as a function of galactic mass in a complete sample. The inferred can be constrained via population studies of the X-ray emission from binaries around Sgr A* and in M31, as well as via dynamics studies of the number density of BH allowed from the orbital parameters of stars in galactic nuclei. The upper limit on can be reduced if we can observationally distinguish between high- and low-accretion rate LINERs. Simulations that include a spherical component of individual stars and BH as well as migrating objects in the disk are required to properly constrain . Encounters between objects from the spherical dynamical component and the disk dynamical component will occur at relatively high velocity and can therefore ionize sufficiently soft, large radius, binaries. Thus, in order for to be moderately large in this channel, we require to be large, since otherwise the rate of ionizing encounters can ionize binaries (Leigh et al., 2017). So limits on from semi-analytic approaches or simulations (Kennedy et al., 2016) can also help constrain .
Uncertainties in are dominated mainly by lack of knowledge of the distribution and number of BH in galactic nuclei, how efficiently gas disks can grind down orbits, and whether geometrically thick disks can efficiently merge BHs. Understanding multiple-object migration and the role of retrograde orbiters is another key area for future work.
4 Constraining BH masses
By merging BHs in AGN disks, we expect ’overweight’ BH to result (McKernan et al., 2012). To investigate the range of BH masses involved in mergers in this channel, we use a toy model calculation of the evolution of a population of BH embedded and migrating in an AGN disk. We made many simplifying assumptions: there are no BH binaries to begin with (), BH remain in the disk after merger, tertiary encounters are neglected, no BHs merge with the SMBH, no new BH are added to the population () and we ignore mass growth due to gas accretion. We began with a uniform distribution of BH drawn from a Kroupa (2002) initial mass function , with distributed over three mass bins (5,10,15) and chose normalization .
where is a numerical factor of order . So the toy model population outlined above will evolve over time. If BH are uniformly distributed across a disk of radius , BH orbits are separated by on average. This separation could be closed in Myr from eqn. (3). Our initial distribution of singleton BH separated by on average will therefore evolve from towards within Myr due to migration. The probability of encounter between BH of masses in time is
When a pair of BHs approaches within their binary Hill radius , where is the binary mass ratio and is the radius of the binary center of mass, gas drag can cause them to merge rapidly. Baruteau et al. (2011) showed that binary semi-major axis halves due to gas drag in only 200(1000) orbits about the binary center of mass for a retrograde (prograde) binary compared to gas velocity. Using this result, a BH binary with at has a characteristic timescale for binary hardening of 0.4 kyr (8 kyr) in the retrograde(prograde) case. Only 20–25 such halvings (corresponding to –0.2 Myr, naively assuming a constant gas hardening rate) would shrink sufficiently that GW emission takes over and the merger happens promptly. The gas hardening rate may be even faster than this estimate since more gas enters the binary’s Hill sphere as it shrinks (Baruteau et al., 2011), which may pump binary eccentricity. However, gas torques may decrease in efficiency once the binary has hardened sufficiently that the binary velocity is substantially supersonic compared to most gas within the Hill radius (Sánchez-Salcedo & Chametla, 2014). For our toy model, we therefore assume Myr is the minimum gas hardening timescale to merger, but we note that the actual gas hardening timescale could take up to an order of magnitude longer.
In our toy model, if the typical time for a BH to encounter another BH in the disk is Myr, then adding an additional Myr for a gas-hardening timescale, yields a characteristic time to merger of Myr in our model. So, we expect that around half the initial population of our toy model will have encountered each other and merged in this time. In calculating the evolution of our toy model, we chose Myr to correspond to a time when of the initial population of lowest mass BHs () have encountered each other and merged. All other encounters are normalized to this encounter rate. For simplicity, we assume all binaries formed in merge within that time, and we neglect the mass-energy loss from the mergers. After , all BH that merged are removed from their original mass bins, and the newly merged object is added to the appropriate mass bin.
Figure 1 demonstrates the simplistic evolution expected as the initial BH distribution (black line) evolves to the red curve in time step Myr, where of the lowest mass BHs in the initial (black) distribution have merged. The red curve evolves to the blue curve after an additional Myr, when of the lowest mass BH on the red curve are expected to merge. The BH mass distribution in our toy model flattens from to as low-mass BH are consumed.
Now assume that BH from the non-disk spherical population, interact with the disk and their orbits are ground down into the disk, i.e. . The addition of some of the (initially) spherical BH population into the disk will support the BH mass distribution in the disk at the low mass end. So an initial power law distribution of BH mass will evolve towards a broken-power law distribution of the form
where , , where is the fraction of BH initially in the disk and on average is the fraction of BH ground down into the disk over and lies near the upper end of the inital mass range ( in our toy model).
In order to include gas accretion in this toy model, we assumed a gas accretion rate for BH on [retrograde, prograde] orbits of , where
is the Eddington mass accretion rate with the proton mass and the accretion luminosity efficiency. Over an AGN disk lifetime of Myr, we can neglect gas accretion onto BH on retrograde orbits.
Note. – Parameter ranges predicted for BH binaries in this channel, assuming initial BH mass range 5–15 M and uniform distribution of BH (see text).
In Table 2 we list parameter ranges for BH masses on the basis of the probabilistic toy model outlined above for three different assumptions: 1) (roughly the blue curve in Fig. 1), corresponding to a short lived disk with . 2) , corresponding either to a long lived disk (Myr) or efficient gas hardening with a low rate of orbit grind down (). 3) for , corresponding either to efficient orbit grind down (), or efficient stellar formation and evolution in the disk with a new top-heavy IMF. In Table 2 we list the binary mass ratio range for each set of assumptions. The lower limit to is trivially the lowest possible mass binary drawn from the initial mass distribution, with no growth from gas accretion and the upper limit to is simply the highest mass binary in the distribution. Also listed in Table 2 are the range of mass ratios () of the binaries in the three different scenarios, with the lower limit given by the range of BH masses allowed in the three different distributions and is the trivial upper limit.
If the fraction of BH ground down into the disk , the fraction of BH coincident with the disk, which will be true for relatively long-lived, thin () disks, the BH mass spectrum evolves from an initial power-law distribution to a broken power-law as in Eqn. (5) with . The uncertainty in mass estimates for this channel is driven mainly by the initial mass distribution of BH in the central region, as well as the ratio of , which in turn depends on disk density and .
5 Range of BH spins
As black holes in the AGN disk accrete gas and merge with each other, their initial spin distribution will change with time. Assuming a uniform distribution of spins () and angular momenta () for BH in galactic nuclei, there will be four distinct populations of BHs in AGN disks as follows:
Prograde spin, on prograde orbits, denoted by ().
Prograde spin, on retrograde orbits ().
Retrograde spin, on prograde orbits ().
Retrograde spin, on retrograde orbits ().
We expect the fraction of BH co-orbital with the AGN disk should have an initial uniform distribution across all four BH populations.
The four BH populations will evolve differently due to gas accretion. The () population rapidly accretes gas, spins up, and aligns spins with the disk gas once the BH has accreted a few of its own mass(Bogdanovic et al., 2007), i.e. in . An initially uniform spin distribution evolves towards at an average rate where is the average gas accretion rate as a fraction of the Eddington rate (which takes Myr to double mass). By contrast, the () population faces a strong headwind, so it accretes very weakly from the gas. An initially uniform distribution of spins in this population will remain uniform over . The () population spins down towards after an increase of mass by a factor (Bardeen, 1970) and will then join the () population. The () population spins down more slowly due to the headwind and so an initial uniform distribution of spins remains uniform over .
BH mergers will further complicate the spin evolution of the four BH populations. The four populations interact due to migration and form binaries if captured within the binary Hill sphere. Binary orbital angular momentum () is the dominant contributor to the spin of the merged BH binary so equal mass BH mergers yield merger products with (Hofmann et al., 2016). Binaries can form with prograde or retrograde orbital angular momentum compared to the disk gas (denoted by ). If a binary forms with retrograde orbital angular momentum (), the merger is faster than in the prograde case (Baruteau et al., 2011), and the merger product will have (i.e. retrograde spin compared to disk gas). Thus the fastest growing of the four populations of BH in the disk due to mergers will actually be . This population evolves towards low spin () due to gas accretion, at an average rate . Among the initial fraction of co-orbital BHs, we expect equal numbers of prograde to retrograde orbits. However, since prograde orbits are ground down faster (smaller headwind, greater Bondi radius), we expect .
Applying all of this to our toy model above allows us to construct the spin distribution in Fig. 2. An initial uniform spin distribution (black line) evolves towards the solid red curve after Myr. The corresponding mass distribution is the red curve in Fig. 1. The red solid curve in Fig. 2 shows a prominent peak at due to a faster merger rate of retrograde binaries and a smaller peak at due to mergers of prograde binaries. Both peaks are smeared out towards the right by gas accretion during and will consist of BH masses from the initial mass distribution. Some pile-up is happening at due to gas accretion onto the already near maximal spinners of the () population. The red dashed curve shows what happens if we assume gas accretion can occur at super-Eddington rates onto BH in the disk ( the Eddington rate). In particular the more massive merged population at gets quickly smeared out and driven towards low spin. Thus, from Fig. 2 if LIGO constrain the spins of most merger precursor BHs to be small, the AGN channel requires super-Eddington accretion onto initially retrograde spin BH to grow this population.
Only the () population will align or anti-align relatively quickly with the AGN disk gas. Assuming the () population are all aligned or anti-aligned with the disk gas, by drawing randomly from a uniform distribution across (), there is a chance that both BH have (anti-)aligned spins and represents our lower limit for the fraction of BH (anti-) aligned with disk gas. If , then effectively the two populations () will dominate so , which is our approximate upper limit for the fraction of BH (anti-) aligned with disk gas. Our estimates of suggest that a larger population of mergers will be requied to test this channel in population spin studies than estimated by Fishbach et al. (2017); Gerosa & Berti (2017). Anti-aligned binaries in the AGN disk allow LIGO a unique chance to test the spin precession instability (Gerosa et al., 2015).
Once a BH binary merges, the resulting merger product can experience a gravitational radiation recoil kick of –400 km s, depending on relative spins and mass ratios (e.g. Merritt et al., 2004; Campanelli et al., 2007). The result of kicks from mergers between aligned and anti-aligned objects is to incline the merger product’s orbit relative to the AGN disk by where is the orbital velocity of the binary center of mass. Since km/s in most of the disk, the orbital inclination perturbation is at most a few degrees and the merger product could be ground back down into the disk in time . Mergers of BH with spins out of alignment with the plane of the disk and each other can produce the largest magnitude kicks (up to several thousand kilometers per second) (e.g. Schnittman & Buonnano, 2007; Lousto et al., 2012). Such mergers will be rare, but will produce large kicks ( in the mass ratio , Campanelli et al. (2007)), escape the disk at angle and may not be ground back down within .
Table 3 summarizes the ranges allowed for spins in this LIGO channel. The typical spin distribution depends on the relative fractions of the four populations of BH in the disk () and their evolution as changes, driven in turn by disk aspect ratio () and the disk gas density and . We expect an initial population uniform across (), but () will grow with the fraction of BH ground-down into the disk. Peaks will arise in the spin distribution at due to mergers and gas accretion will drive and independent of mergers. Gas accretion at super-Eddington rates plus faster mergers by retrograde binaries may be required to generate a population of overweight, low spin BH in the AGN disk.
Note. – Parameter ranges allowed for BH spins in this channel (see text).
6 Observational Constraints: GW
Binary black hole mergers in an AGN disk imply unique, testable predictions that would not be expected from other BH merger channels, including: 1. A spin distribution (see §5) that includes aligned/anti-aligned spin binaries and 2. a population of overweight BH or IMBH orbiting SMBHs, generating GWs detectable with the Laser Interferometer Space Antenna (LISA) (McKernan et al., 2014). A circularized IMBH-SMBH binary at a migration trap () around a SMBH with will be detectable with LISA at modest signal-to-noise ratio in a year’s observation (McKernan et al., 2014). If AGN disks are efficient at gas-driven mergers of BH, we expect that every AGN must contain one or more IMBH-SMBH binaries, implying an approximate rate comparable to that in Portegies Zwart et al. (2006).
7 Observational Constraints:EM
The brightest AGN are too bright compared to any short-term EM signal that might result from a BH merger in a gas disk. Low luminosity AGN might permit short timescale EM events from BH mergers to be visible. As IMBHs grow in migration traps, gaps and cavities in the accretion flow can form and oscillations on the dynamical timescale of the accreting IMBH can be detected in optical, UV, and X-ray spectral signatures (e. g. McKernan et al., 2013, 2014; McKernan & Ford, 2015). Temporal and energetic asymmetries in the X-ray signatures are best detected using micro-calorimeters, such as the one that will fly on the X-ray Astronomy Recovery Mission succeeding Hitomi. Perturbations of the innermost disk will occur as migrators in the disk plunge into the SMBH and temporarily dominate the local co-rotating mass, detectable in large UV-optical quasar surveys (Drake et al., 2009) as well as the X-ray band. Large optical surveys of quasar disks can also limit total supernova rates due to migrating/accreting/colliding stars (Graham et al., 2017), in turn placing limits on the disk populations of stars and stellar remnants. Estimates of the rates of transits by bloated stars, best detected in the X-ray band (McKernan & Yaqoob, 1998), can put limits on the population on spherical orbits around and passing through AGN disks.
As the AGN phase ends, remaining BH will interact dynamically, so the distribution of orbital parameters of the BHs and stars entrained in the disk will relax. Alexander et al. (2007) show that if very massive stars () exist in our own Galactic nucleus, they can pump the eccentricity distribution of massive stars to even within 5 Myrs. However, such stars are short-lived and observed stellar eccentricities reach (Paumard et al., 2006). On the other hand, a population of overweight BHs caused by merger in an AGN disk can rapidly pump stellar orbital eccentricites post-AGN and inflate the thickness () of stellar disks in galactic nuclei. Thus, if this BH merger channel is efficient, thin disks of stars will not be observed in post-AGN galactic nuclei.
Neutron stars (NS) should also exist in AGN disks, and can migrate. So there should be a correlation between NS-NS and NS-BH mergers in AGN disks and the rate of BH-BH mergers expected from this channel. No correlation has been observed so far between short gamma-ray bursts in the local universe and AGNs (Berger, 2014), but so far, only a handful of short gamma-ray bursts have sufficiently accurate positions in the sky to rule out an association with AGN in these cases. The efficiency of this LIGO channel could be further constrained by ongoing studies of the correlation of short gamma-ray bursts with AGN. Future simulations could usefully focus on the expected distribution of NS in mass segregating clusters in galactic nuclei, and ultimately on determining the expected NS merger rate in AGN disks.
We parameterize the rate of black hole mergers within AGN disks and the mass and spin distributions that result. The strongest observational constraints can be placed on this channel by: 1. ruling out a population of maximal spin BH via LIGO, 2. ruling out a correlation betwen short gamma-ray bursts and AGN, 3. constraining the rate of obscured supernovae in AGN disks via studies of large samples of AGN, 4. ruling out a population of high accretion rate ADAFs in galactic nuclei and 5. observing very thin disks of stars in nearby Galactic nuclei. Future simulations should focus on 1. the ratio of NS/BH in nuclear star clusters undergoing mass segregation, 2. encounters between prograde and retrograde orbiters in AGN disks and 3. interactions and binary formation between BHs with pro- and retro-grade spins and orbits at migration traps in a range of AGN disk models. If AGN are efficient at merging BH, LISA will detect a large population of IMBH in disks around SMBH in the nearby Universe.
Thanks to Maya Fishbach, Davide Gerosa, Matthew Graham, Daniel Holz, Dan Stern and Nick Stone for useful conversations. BM & KESF are supported by NSF PAARE AST-1153335 and NSF PHY11-25915. BM & KESF thank CalTech/JPL and NASA GSFC for support during sabbatical. M-MML is partly supported by NSF AST11-09395.
- Abbott et al. (2016a) Abbott B.P. et al., 2016, PhRvL, 116,1102
- Abbott et al. (2016b) Abbott B.P. et al., 2016, ApJL, 833,L1
- Antonini et al. (2014) Antonini F., 2014, ApJ, 794, 106
- Antonini & Rasio (2016) Antonini F. & Rasio F., 2016, ApJ (submitted), arXiv:1606.04889
- Alexander et al. (2007) Alexander R.D., Begelman M.C. & Armitage P.J., 2007, ApJ, 654, 907
- Baldry et al. (2012) Baldry I.K. et al., 2012, MNRAS, 421, 621
- Bardeen (1970) Bardeen J.M., 1970, Nature, 226,64
- Bartos et al. (2017) Bartos I., Kocsis B., Haiman Z. & Márka S., 2017, ApJ, 835,165
- Baruteau et al. (2011) Baruteau C., Cuadra J. & Lin D.N.C., 2011, ApJ, 726, 28
- Bellovary et al. (2016) Bellovary J., Mac Low M.-M., McKernan B. & Ford K.E.S., 2016, ApJ, 819, L17
- Belczynski et al. (2010) Belczynski K., Bulik T., Fryer C.L., Ruiter A., Valsecchi F., Vink J.S. & Hurley J.R., 2010, ApJ, 714, 1217
- Berger (2014) Berger E., 2014, ARA&A, 52, 43
- Bogdanovic et al. (2007) Bogdanovic T., Reynolds C.S. & Miller M.C., 2007, ApJ, 661, L147
- Campanelli et al. (2007) Campanelli M., Lousto C.O., Zlochower Y. & Merritt D., 2007, PhRvL, 98, 231102
- Cole et al. (2001) Cole S. et al., 2001, MNRAS, 326, 255
- Drake et al. (2009) Drake A.J. et al., 2009, ApJ, 696, 870
- deMink & Mandel (2016) deMink S. & Mandel I., 2016, MNRAS, 460, 3545
- Fishbach et al. (2017) Fishbach M., Holz D.E. & Farr B., 2017, ApJ, 840, L24
- Gerosa et al. (2015) Gerosa D. et al., 2015, Phys Rev Lett, 115, 141102
- Gerosa & Berti (2017) Gerosa D. & Berti E., 2017, Phys Rev D (submitted), arXiv:1703.06223
- Graham et al. (2017) Graham M. et al., 2017, MNRAS, submitted
- Haehnelt & Rees (1993) Haehnelt M.G. & Rees M., 1993, MNRAS, 263, 168
- Haiman et al. (2009) Haiman Z., Kocsis B. & Menou K., 2009, ApJ, 700, 1952
- Ho (2008) Ho L.C., 2008, ARA&A, 46, 475
- Hofmann et al. (2016) Hofmann F., Barausse E. & Rezzolla L., 2016, arXiv:1605.01938
- Hopman & Alexander (2006) Hopman C. & Alexander T., 2006, ApJ, 645, L133
- Horn et al. (2012) Horn B., Lyra W., Mac Low M.-M. & Sándor Z., 2012, ApJ, 750, 34
- Kennedy et al. (2016) Kennedy G. et al., 2016, MNRAS, 460, 240
- King & Nixon (2015) King A. & Nixon C.J., 2015, MNRAS, 453, L46
- Kocsis et al. (2011) Kocsis B., Yunes N. & Loeb A., 2011, PRD, 84, 024032
- Kroupa (2002) Kroupa P. 2002, Science, 295, 82
- Lasota et al. (2016) Lasota J.-P. et al., 2016, A&A, 587, 13
- Leigh et al. (2016) Leigh N. W. C., Antonini F., Stone N. C., Shara M. M., Merritt D. 2016, MNRAS, 463, 1605
- Leigh et al. (2017) Leigh N. W. C., Geller, A. M., McKernan, B., Ford, K. E. S., Mac Low, M.-M., Bellovary, J., Haiman, Z., Lyra, W., Samsing, J., O’Dowd, M., Kocsis, B., Endlich, S. MNRAS, submitted (ArXiv:TBD)
- Lousto et al. (2012) Lousto C.O., Zlochower Y., Dotti M. & Volonteri M., 2012, PRD, 85, 084015
- McKernan & Yaqoob (1998) McKernan B. & Yaqoob T., 1998, ApJ, 501, L29
- McKernan et al. (2011) McKernan B. et al, 2011, MNRAS, 417, L103
- McKernan et al. (2012) McKernan B., Ford K.E.S., Lyra W. & Perets H.B., 2012, MNRAS, 425, 460
- McKernan et al. (2013) McKernan B., Ford K.E.S., Kocsis B. & Haiman Z., 2013, MNRAS, 432, 1468
- McKernan et al. (2014) McKernan B., Ford K.E.S., Kocsis B., Lyra W. & Winter L.M., 2014, MNRAS, 441, 900
- McKernan & Ford (2015) McKernan B. & Ford K.E.S., 2015, MNRAS, 452, L1
- Miller & Davies (2012) Miller M. C., Davies M. B. 2012, ApJ, 755, 81
- Merritt et al. (2004) Merritt D., Milosavljević M., Favata M., Hughes S.A. & Holz D.E., 2004, ApJ, 607, L9
- Miralda-Escudé & Gould (2000) Miralda-Escudé J. & Gould A., 2000, ApJ, 545, 847
- Morris (1993) Morris M., 1993, ApJ, 408, 496
- Narayan & Yi (1995) Narayan R. & Yi I., 1995, ApJ, 444, 231
- O’Leary et al. (2009) O’Leary R.M., Kocsis B. & Loeb A., 2009, MNRAS, 395, 2127
- O’Shaugnessey, Gerosa & Wysocki (2017) O’Shaughnessey R., Gerosa D. & Wysocki D., 2017, Phys. Rev. Lett, accepted, arXiv:1704.03879
- Paardekooper et al. (2010) Paardekooper S.-J., Baruteau C., Crida A. & Kley W., 2010, MNRAS, 401, 1950
- Paczynski & Witta (1980) Paczynski B. & Witta P.J., 1980, A&A, 88, 23
- Paumard et al. (2006) Paumard T. et al., 2006, ApJ, 643, 1011
- Portegies Zwart et al. (2006) Portegies Zwart S.F. et al., 2006, ApJ, 641, 319
- Reines & Volonteri (2015) Reines A.E. & Volonteri M., 2015, ApJ, 813, 82
- Rodriguez et al. (2016) Rodriquez C. et al. 2016, arXiv. etc.
- Samsing et al. (2014) Samsing J., MacLeod M. & Ramirez-Ruiz E. 2014, ApJ, 784, 71
- Sánchez-Salcedo & Chametla (2014) Sánchez-Salcedo F. J. & Chametla R.O. 2014, ApJ, 794, 167
- Schawinski et al. (2015) Schawinski K., Koss M., Berney S. & Sartori L.F. 2015, ApJ, 451, 2517
- Schnittman & Buonnano (2007) Schnittman J.D. & Buonanno A. 2007, ApJ, 662, L63
- Sirko & Goodman (2003) Sirko E. & Goodman J. 2003, MNRAS, 341, 501
- Stahler (2010) Stahler S.W., 2010, MNRAS, 402, 1758
- Stone et al. (2017) Stone N.C. Metzger B.D. & Haiman Z., 2017, MNRAS, 464, 946
- Thompson et al. (2005) Thompson T.A., Quataert E. & Murray N. 2005, ApJ, 630, 167