High velocity stars from GC-SMBHB interaction

# High velocity stars from the interaction of a globular cluster and a massive black hole binary

## Abstract

High velocity stars are stars moving at velocities so high to require an acceleration mechanism involving binary systems or the presence of a massive central black hole. In the frame of a galaxy hosting a supermassive black hole binary (of total mass M), we investigated a mechanism for the production of high velocity stars due to the close interaction between a massive and orbitally decayed globular cluster and the super massive black hole binary. Some stars of the cluster acquire high velocities by conversion of gravitational energy into kinetic energy deriving from their interaction with the black hole binary. After the interaction, few stars reach a velocity sufficient to overcome the galactic gravitational well, while some of them are just stripped from the globular cluster and start orbiting around the galactic centre.

###### keywords:
galaxies: haloes – galaxies: nuclei – galaxies: star clusters: general; stars: kinematics and dynamics

## 1 Introduction

Massive Black Holes (BHs) are likely in the majority of galactic nuclei over the whole Hubble sequence and are recognized as fundamental building blocks in models of galaxy formation and evolution (Kormendy & Ho, 2013; Dubois et al., 2014; Sesana et al., 2014). According to the present cosmological model CDM, structures formation involves hierarchical mergers of galaxies, which follow the hierarchical growth of their parent dark matter halos (Mayer at al., 2007; Roškar at al., 2015). Galaxies may experience multiple mergers during their lifetime and, if more than one of these interacting galaxies host a massive BH, the formation of a Black Hole Binary (BHB) is a natural consequence of the hierarchical paradigm (Begelman, Blandford & Rees, 1980; Volonteri, Haardt & Madau, 2003). How long BHBs survive and whether they merge are key points for several questions in high-energy and extragalactic astronomy, as the detection of gravitational waves (Sesana et al., 2005; Sesana, 2013; Rosado & Sesana, 2014).

BHBs lifetime is governed by the loss of energy and angular momentum, by means of which two massive BHs become gravitationally bound and finally merge. BHBs are thought to undergo four dynamical stages before their coalescence (Begelman, Blandford & Rees, 1980; Yu, 2002). The first one is the dynamical friction stage, during which each BH inspirals independently toward the centre of the common gravitational potential on the Chandrasekhar time-scale (Chandrasekhar, 1943)

 tdf∼4×106logN(σc200 km s−1)(rc100 pc)2(108M⊙mi)yr , (1)

where , , are the number of stars, the one-dimensional velocity dispersion and the radius of the core, respectively, while () is the BH mass. During the second stage, usually referred to as non-hard binary stage, BHs velocities increase, while their orbital period shortens, because more and more stars in the galactic nuclear core are scattered off through gravitational slingshots, while the dynamical friction becomes less efficient. The third stage, labelled hard binary stage, begins when the orbital separation of BHs is about their influence radius (Quinlan, 1996)

 ah=2.8(m2108M⊙)(200 km s−1σc)2pc , (2)

where is the mass of the lighter BH. Hard BHBs lose energy mainly by the three-body slingshot effect with stars passing in their vicinity, which are expelled after one or more encounters. The span of this phase depends on the loss-cone refilling. If it is dominated by two-body relaxation in a spherical and symmetric system, energy will not be lost efficiently and the binary may stall at parsec scale (the so-called ”final parsec problem”). However, more realistic calculations show that this stall does not sussist when, for example, asymmetric potentials or gas dynamics are taken into account (Merritt & Milosavljević, 2005; Berczik et al., 2006; Merritt, Mikkola & Szell, 2007; Khan et al., 2013). Finally, the last stage is characterised by the energy loss due to gravitational radiation. A BHB will coalesce within the time (Peters, 1964)

 tcoal∼5.8 m21×106m2M(a0.01 pc)4(108M⊙m1)3, (3)

where is the mass of the heavier BH and is the total mass.

On the point of view of observations, Super Massive BHBs (SMBHBs) are quite hard to detect, in particular in dry environments. However, the presence of surrounding gas enhances the possibilities of finding and resolving BHBs in galactic centres thanks to their energetic interactions (Komossa et al., 2003; Rodriguez et al., 2006; Fabbiano et al., 2011; Deane et al., 2014). Moreover, in these gas-rich systems, the BHB evolution is likely to be driven primarily by gas (Goicovic et al., 2015). On the other hand, in gas-poor environments, when tidal disruption events enlighten the cores of galaxies hosting quiescent BHBs, observers are able to detect and study them (Chen et al., 2009; Liu, Li & Chen, 2009; Chen et al., 2011; Wegg & Nate Bode, 2011; Liu & Chen, 2013; Liu, Li & Komossa, 2014).

Many galaxies show nucleated central regions, the so-called Nuclear Star Clusters (NSCs), which are among the densest stellar populations observed in the Universe (Matthews & Gallagher, 1997; Carollo et al., 1997, 1998; Côté et al., 2006; Turner et al., 2012). Two different processes are thought to give birth to NSCs. The first mechanism involves radial gas inflow into the galactic centre and predicts that NSCs consist mostly of stars formed locally (Loose, Kruegel & Tutukov, 1982; Milosavljević, 2004; Bekki, 2007; Schinnerer et al., 2006, 2008). The second mechanism suggests that massive stellar clusters, such as Globular Clusters (GCs), spiral into the galactic centre and merge to form a dense nucleus (Tremaine et al., 1975; Capuzzo-Dolcetta, 1993; Capuzzo-Dolcetta & Miocchi, 2008; Antonini et al., 2012; Antonini, 2013). Observations show that both processes occur in nature and contribute to the formation of NSCs. In the latter scenario, GCs are expected to be totally disrupted by BH tidal forces, but some stars may be accelerated to high velocities and ejected in jets from the inner galactic regions (Capuzzo-Dolcetta & Fragione, 2015; Arca-Sedda, Capuzzo-Dolcetta, Spera, 2015).

High velocity stars have been observed in the Galactic halo. Such high velocities may be gained thanks to several physical mechanisms, as three-body interactions among binary systems or with the massive black hole in the Galactic centre, or kicks due to supernovae explosions. High velocity stars can be divided in two different categories, i.e. runaway stars (RS) and hypervelocity stars (HVSs).

RSs, historically defined in the context of O and B stars (Humason & Zwicky, 1947), are Galactic halo stars with peculiar motions higher than km s. Such young massive stars are not expected to be present far from star-forming regions and are thought to have travelled far from their birthplace. Two mechanisms are expected to produce RSs: supernova ejections and dynamical ejections (Silva & Napiwotzki, 2011). In the first scenario (Blaauw, 1961; Portegies Zwart, 2000; Scheck et al., 2006; Przybilla et al., 2008), a RS has origin in a binary system, where it receives a velocity kick when its companion explodes as a supernova. In the dynamical ejection scenario (Poveda, Ruiz & Allen, 1967), stars are accelerated thanks to a three- or four-body interaction (Leonard & Duncan, 1990; Gvaramadze, 2009; Gvaramadze, Gualandris & Portegies Zwart, 2009; Gvaramadze & Gualandris, 2011; Perets & Subr, 2012). Observations show that both the ejection mechanisms operate in nature (Hoogerwerf et al., 2001), but, in any case, RSs velocities are under the Galaxy escape velocity.

HVSs are stars escaping the host Galaxy. Hills (1988) was the first to predict theoretically their existence, while Brown et al. (2005) serendipitously discovered the first HVS in the outer halo. Hills’ mechanism involves the tidal breakup of a binary passing close to a massive BH and has been investigated in literature by several authors (Yu & Tremaine, 2003; Gualandris, Portegies Zwart & Sipior, 2005; Bromley et al., 2006; Sari, Kobayashi & Rossi, 2009; Kobayashi et al., 2012; Rossi, Kobayashi & Sari, 2014). Moreover, Hills’ scenario predicts the existence of a population of stars orbiting in the inner Galactic regions around the central BH (Gould & Quillen, 2003; Ginsburg & Loeb, 2006; Perets, Hopman & Alexander, 2007). Also other mechanisms have been proposed to explain the production of HVSs (Tutukov & Federova, 2009; Brown, 2015), as the interaction of a SMBHB with a single star (Gualandris, Portegies Zwart & Sipior, 2005; Baumgardt, Gualandris & Portegies Zwart, 2006; Sesana, Haardt & Madau, 2006), the arrival from another nearby galaxy (Gualandris & Portegies Zwart, 2007; Bonanos et al., 2008; Sherwin, Loeb & O’Leary, 2008; Perets, 2009; Brown et al., 2010) and supernova explosion in a close binary (Zubovas, Wynn & Gualandris, 2013).

Since high velocity stars production mechanisms involve different astrophysical phenomena, it would be possible to infer information about different branches of physics, as the physics of regions near massive BHs (Gould & Quillen, 2003; Sesana, Haardt & Madau, 2007; O’Leary & Loeb, 2008) and of supernovae (Portegies Zwart, 2000; Zubovas, Wynn & Gualandris, 2013). Moreover, the study of the proper motions of such fast moving stars can improve the knowledge of the Galaxy gravitational potential shape and of its Dark Matter component (Gnedin et al., 2005; Yu & Madau, 2007).

Observations of high velocity and hypervelocity objects have usually been limited to high-mass, early-type stars due to obvious observational bias (Brown, Geller & Kenyon, 2009, 2010, 2014). Nowadays, observers have started investigating low-mass high velocity stars (Palladino et al., 2014; Zhong et al., 2014; Li et al., 2015; Vickers, Smith & Grebel, 2015), some of which are low-mass HVSs candidates. Moreover, the European ESA satellite Gaia (http://www.cosmos.esa.int/web/gaia), along with Gaia-ESO (http://www.gaia-eso.eu/), is expected to measure proper motions with an unprecedented precision, providing for a larger and less biased sample, and to find HVSs in a sample of stars.

The aim of this paper is to investigate the production mechanism of high velocity stars, which involves a GC and a SMBHB. Actually, when the orbit of a GC delivers it close to a BHB in the centre of its host galaxy due to dynamical friction, some stars are stripped from the cluster and may be ejected with high velocities. For these test cases, we assumed a total mass M with the scope of better underlying the physical mechanism and of making comparison with our previous results in the case of a single BH (Capuzzo-Dolcetta & Fragione, 2015).

The paper is organised as follows. In Section 2 we outline and describe our approach to the study of the consequences of the GC-BHB interaction. In Section 3 the results of our scattering experiments are presented and discussed. Finally, in Section 4 we draw the conclusions.

## 2 Method

Our scattering experiments refer to the interaction of three different components: a SMBHB, a GC and a star. In our simulations the SMBHB centre of mass sits initially in the origin of the reference frame, while the GC follows an elliptical orbit at a relatively close distance around it. The assumption of close distance to the BH is motivated by the fact that the globular cluster is supposedly orbitally decayed by dynamical friction braking suffered in the inner dense galactic regions.

The SMBHB initial circular orbit has radius . Therefore, for the heavier black hole the initial conditions are given by

 r1,c=m2Mah   ;   v1,c=√GMahm2, (4)

while for the lighter black hole by

 r2,c=m1Mah   ;   v2,c=√GMahm1, (5)

where and are the radii of the circular orbits of and , respectively, around their centre of mass.

The mechanical energy (per unit mass) of the GC on a circular orbit of radius , around the BHB center of mass (neglecting the stellar background), is

 Ec≡12v2c−GMrc=−12GMrc, (6)

where is the circular velocity. Thereafter, taking into account that the angular momentum (per unit mass) of the GC for the circular orbit is , the pericenter () and apocenter () distances of the GC on orbits of same energy (), but different angular momentum , are given by

 r±=rc⎛⎝1±√1−(LLc)2⎞⎠, (7)

where the sign gives the pericenter and the sign gives the apocenter. Therefore, by varying the ratio , we can compare the circular orbit with a set of orbits at same energy, but different eccentricity

 e=r+−r−r−+r+=√1−α. (8)

Furthermore, in our simulations, the GC is assumed to have a Plummer (1911) mass profile

 M(r)=MGCr3(r2+b2)3/2, (9)

where is the total mass of the GC and its core radius. The test star is on circular orbit inside the GC sphere of influence.

The cartesian reference frame has been chosen as that with the -axis along the line connecting the GC with the SMBHB center of mass and -axis orthogonal, so that the frame is equiverse to the GC orbital revolution.

To summarize, the values of the relevant initial parameters have been set as follows (see also Table 1):

• the total mass of the SMBHB is M;

• the binary ratio assumes the values of , , , , and ;

• the radius of the SMBHB initial circular orbit is set equal to ;

• the initial phase of the BHB is randomly generated;

• the GC mass, , is fixed to M;

• the GC core radius is set to pc;

• the radius of the GC reference circular orbit is pc;

• the GC orbital eccentricity ranges from () to () and is parametrized varying at steps of ;

• the GC orbits are coplanar with the SMBHB one;

• the test star mass, , is set equal to M;

• the test star circular orbit radii are equal to , , , , where is the radius of the Hill’s sphere;

• the star initial position on the circular orbits is randomly generated;

• the angles , , , which determine the orientation of the star circular orbit in the GC reference frame, are randomly generated.

The choice of the range of , and consequently of , toward large values of , is justified by the fact that the efficiency of the energy transfer on the test star orbiting the GC tends to vanish at eccentricities less than .

Given the above set of initial parameters, we integrated the system of the differential equations of the 4-bodies (SMBHB, GC and star) motion

 ¨ri=−G∑j≠imj(ri−rj)∣∣ri−rj∣∣3, (10)

for , using the fully regularized algorithm of Mikkola & Aarseth (2001). The enormous range of variation of the involved masses requires a regularized algorithm, since any not-regularized direct summation code would fail in dealing with the close star-SMBHB interaction and would carry to a large energy error during the close encounter. Thanks to a transformed leapfrog method, combined with the Bulirsch-Stoer extrapolation method, the Mikkola’s ARW code overcomes this problem and leads to extremely accurate integrations of the bodies trajectories (Bulirsch & Stoer, 1966; Mikkola & Tanikawa, 1999a, b; Mikkola & Merritt, 2006, 2008; Hellström & Mikkola, 2010). Thanks to the regularized algorithm, the fractional energy error is kept below over the whole integration time.

## 3 Results

In our scattering experiments the test star orbiting the GC has three possible fates after the interaction with the SMBHB. Actually, the star can either remain bound to the GC, but on an orbit significantly perturbed respect to the original one, or can be captured by the SMBHB and starts orbiting around the galactic centre on precessing loops, or can be lost by the GC-SMBHB system.

The distinction among these three different situations is made by computing the mechanical energy of the star respect to the SMBHB and the GC after the scattering. If its energy respect to the GC remains negative, the star remains bound to the GC, while if this energy becomes positive, while the star energy respect to the BHB is negative, the star becomes bound to the binary. Finally, if both these energies are positive, the star is able to leave the SMBHB-GC system.

In this paper, we focus our attention on the stars ejected at high velocities after the close interaction with a SMBHB, studying the effects of the mass ratio and of the SMBHB orbital eccentricity.

### 3.1 The role of the mass ratio

From our scattering experiments, it is possible to derive the velocity profile of the ejected stars. In the case of single massive BH, the velocity distribution is narrow and peaked at small velocities, depending on the GC core radius. Actually, when the GC has a Plummer profile (Eq. 9), for the set of parameters chosen in this work, the initial gravitational energy of the star is (in our simulations pc). Therefore, the distribution is peaked at a lower velocity respect to the case of point-mass GC (Capuzzo-Dolcetta & Fragione, 2015). The ejection velocity roughly depends on the initial amount of star gravitational energy and on the pericentre of the GC orbit. Moreover, the dispersion of the nearly Gaussian distribution is determined, besides by the different (randomly generated) angles and initial conditions, by the relative inclination between the angular momentum of the cluster and that of the test star. Actually, if during the close interaction, i.e. near the pericentre of the GC orbit respect to the BHB centre of mass, the ejection velocity of the star will be higher than the cases in which . Fig. 1 shows the resulting nearly Gaussian velocity distribution, for all the values of , peaked at km s in the case of single BH (black line).

The situation is different when a scattering with a BHB is taken into account. Actually, the ejection velocity depends not only on , and , but also on , and so . The generic star of an infalling GC is able to exchange energy with the binary through the gravitational slingshot, which enhances its ejection velocity. On the other hand, the energy transfer makes the SMBHB shrink and reduces the apocentre of its orbit. Fig. 1 shows the resulting velocity distribution for BHB for different values of . If the mass ratio is high (top panel), a considerable tail in the distribution up to km s is produced. On the contrary, if the mass ratio is low (bottom panel), the differences between the distributions are not so pronounced. In particular, if , the velocity distributions for single and binary BH are essentially similar.

Velocity distributions show that for high values of , there is a considerable fraction of the ejected stars, which acquire significant high velocities. Actually. for , the huge tail of the distribution extends up to km s. Therefore, a considerable fraction of ejected stars is unbound not only respect to the SMBHB-GC system, but also respect to the host galaxy, becoming HVSs. The escape velocity (at pc as justified in Capuzzo-Dolcetta & Fragione (2015)) for an elliptical galaxy, as NGC 3377 (Marconi & Hunt, 2003), of total mass M, is km s (Hernquist, 1990), while for a spiral galaxy, of total mass M (Wang, Sulkanen & Lovelace, 1992; Rownd, Dickey & Helou, 1994; Kornreich & Lovelace, 2008; Lingam, 2014), is km s (Miyamoto & Nagai, 1975; Binney, 1981; Hernquist, 1990; Fujita, 2009). Fig. 1 shows that a not-negligible fraction of ejected stars is beyond the local escape velocity. The fate of these stars is to escape the SMBHB-GC system, to travel across the halo and to leave the host galaxy.

The introduction of a secondary BH, comparable in mass with the primary, leads to a peculiar distribution velocity for the stars ejected at high velocities. In principle, these distributions can be used, when data for proper motions and radial velocities will be available for galaxies whose central object(s) total mass is M, to distinguish if the central object is a single BH or a BHB.

The effect of the presence of a secondary BH is not limited to the velocity distribution of ejected stars. Actually, also the branching ratios of ejected stars, i.e. the probabilities that the system BHB-GC loses stars, depend on the mass ratio . Fig. 2 shows the branching ratios for the case of single and binary BH. In the case of single BH (), the branching ratio is for the set of parameter studied in this work. For low values of the mass ratio, i.e. , the branching ratio remains nearly constant, while for higher values, it is an increasing function of . Therefore, only if the mass ratio is sufficiently high, the probability increases respect to the case of single BH.

In conclusion, the effect of the binariety, of course if , is dual. On the one hand it enhances the probability of stars ejections, on the other it produces a considerable tail in the velocity distribution, leading to the production of HVSs.

### 3.2 The role of the binary eccentricity

In order to explore the role of the BHB orbital eccentricity , we performed the same set of simulations, in the case of mass ratio , for , , elliptical orbits. The angular momenta are chosen such that, as in the case of GC, the energy of elliptical orbits is equal to the energy of the reference circular orbit () of radius pc. In this case, the initial conditions, assuming that the black holes start their orbital motion at the apocentre, are

 r1=(1+η)r1,c   ;   v1=F(η)v1,c=√1−η1+ηv1,c, (11)

for the heavier black hole , while for the lighter black hole

 r2=(1+η)r2,c   ;   v2=F(η)v2,c=√1−η1+ηv2,c. (12)

On the contrary, if the black holes start their orbital motion at the pericentre, and . Finally, for this set of simulations, the angle between the semi-major axis of the SMBHB and the x-axis of the cartesian reference frame is randomly generated.

Fig. 3 shows the resulting velocity distributions for a BHB with mass ratio and different orbital eccentricities. It is clear that the tail in the velocity distribution is produced independently on , with velocities up to km s. However, different eccentricities do not change the shape of the distribution function respect to the case of the circular orbit. Therefore, these results suggest that, although the different shape of the binary orbit, the energy exchange between the generic star of the GC and the BHB is independent on the eccentricity, but depends only on the initial total energy of the binary.

For what regards the branching ratios, Tab. 2 shows the resulting branching ratios of our scattering experiments as function of . From Tab. 2, it is clear that the branching ratio remains nearly constant for different eccentricities. Therefore, not only the shape of the velocity distribution is preserved, but also the probability of stars ejection.

In conclusion, the BHB orbital eccentricity does not effect neither the shape of the velocity distribution nor the branching ratio of ejected stars. Therefore, our simulations suggest that the features of the ejected stars depend only on the binary mass ratio , but not on the binary eccentricity.

## 4 Conclusions

The existence of high velocity stars (Silva & Napiwotzki, 2011) has been explained thanks to a dynamical ejection mechanism (Poveda et al., 1967; Gvaramadze et al., 2009) or a supernova ejection mechanism (Blaauw, 1961; Portegies Zwart, 2000), which are able to accelerated stars up to several hundreds km s. On the other hand, HVSs (Hills, 1988) require the presence of a massive BH, or massive BHB, due to their extreme velocities up to thousands km s (Yu & Tremaine, 2003; Brown, 2015).

In this paper, we extended the recently study Capuzzo-Dolcetta & Fragione (2015), referred to the case of a single BH. Here we generalised the study to the ejection of high velocity stars due to the close passage of a massive GC near a massive BHB. In the frame of the CDM cosmological model, SMBHBs are a natural consequence of the hierarchical paradigm (Begelman, Blandford & Rees, 1980; Volonteri, Haardt & Madau, 2003). Actually, galaxies may experience multiple mergers during their lifetime and, if more than one of them host a massive BH, the formation of a Black Hole Binary is a natural result.

In our study, we assumed a total mass M of the SMBHB to have a direct comparison with previous scattering experiments with a single BH (Capuzzo-Dolcetta & Fragione, 2015). The underlying mechanism is a four-body interaction, where the bodies are the SMBHB, the GC (of mass M in our study) and a generic test star ( M) belonging to the cluster. We performed a series of high precision scattering experiments in order to investigate the ejection probability, after a close interaction with a SMBHB, and the distribution velocity of the ejected stars.

In Capuzzo-Dolcetta & Fragione (2015) results, the velocity distribution was nearly Gaussian, narrow and peaked at low values in dependence on the GC core radius. The ejection velocity is roughly determined by the initial amount of star gravitational energy and by the pericentre of the GC orbit, while the dispersion of the velocity distribution by the different initial conditions and by the relative inclination angle between the angular momentum of the cluster and that of the test star. Moreover, when the interaction occurs with a SMBHB, GC stars are able to exchange energy through gravitational slingshots, which enhances its ejection velocity and produces a tail in the distribution, which depends on the mass ratio of the BHB. Actually, the higher is the mass ratio, the more extended the distribution toward high velocities will be. At the same time, also the branching ratio for ejected stars depends on the mass ratio . While for low values of , the branching ratio remains nearly at the value of the single BH case (), for , it is an increasing function of the mass ratio.

Finally, we performed the same set of simulations in the case of mass ratio , but assuming different eccentricities . The shape of the velocity distribution and the values of the branching ratios show that eccentricity has not a substantial effect on the results. Therefore, our simulations suggest that the features of the ejected stars depend only on the binary mass ratio.

In conclusion, the effect of the binariety (if ) is dual, since it enhances the probability of stars ejections and produces an extended tail in the velocity distribution, corresponding to the production of HVSs. Finally, we suggest that the infall of a GC on a SMBHB could enhance the energy loss by the BHB and conduce them to the final stage, where gravitational waves emission is the main mechanism to cause the energy loss of the binary.

## Acknowledgments

We thank S. Mikkola for making available to us his ARW code and for useful discussions about his use.

### References

1. Antonini F., 2013, ApJ, 763.1, 62
2. Antonini F., Capuzzo-Dolcetta R., Mastrobuono-Battisti A., Merritt D., 2012, ApJ, 750, 111
3. Arca-Sedda M., Capuzzo-Dolcetta R., Spera M., 2015, preprint (arXiv:1510.01137v1)
4. Baumgardt H., Gualandris A., Portegies Zwart S.F., 2006, J. Phys.: Conf. Ser., 54.1, 301
5. Begelman M.C., Blandford R.D., Rees M.J., 1980, Nature, 287, 307
6. Bekki K. 2007, Publ. Astron. Soc. Austr., 24, 77B
7. Berczik P., Merrit D., Spurzem R., Bischof H.-P., 2006, ApJ Lett., 642.1, L21
8. Binney J., 1981, MNRAS, 196, 455
9. Blaauw A., 1961, Bull. Astron. Inst. Netherlands, 15, 265
10. Bonanos A.Z., Lopez-Morales M., Hunter I., Ryans R.S.I., 2008, ApJ, 675, L77
11. Bromley B.C., Kenyon S.J., Geller M.J., Barcikowski E., Brown W.R., Kurtz M.J., 2006, ApJ, 653, 1194
12. Brown W.R., 2015, Annu. Rev. Astron. Astrophys., 53, 15
13. Brown W.R., Geller M.J., Kenyon S.J., 2009, ApJ, 690.2, 1639
14. Brown W.R., Geller M.J., Kenyon S.J., 2014, ApJ, 787, 89
15. Brown W.R., Geller M.J., Kenyon S.J., Kurtz M.J., 2005, ApJ Lett., 622, L33
16. Brown W.R. et al., 2010, ApJ Lett., 719, L23
17. Bulirsch R., Stoer J., 1966, Numer. Math., 8, 1
18. Capuzzo-Dolcetta R., 1993, ApJ, 415, 616
19. Capuzzo-Dolcetta R., Fragione G., 2015, MNRAS, 2015, 454.3, 2677-2690
20. Capuzzo-Dolcetta R., Miocchi P., MNRAS Lett., 388, L69
21. Carollo C. M., Stiavelli M., de Zeeuw P.T., Mack J., 1997, AJ, 114, 2366
22. Carollo C.M., Stiavelli M., Mack J. 1998, AJ, 116, 68
23. Chandrasekhar S., AJ, 1943, 97, 255
24. Chen X., Madau P., Sesana A., Liu F.K., 2009, ApJ, 697, L149
25. Chen X., Sesana A., Madau P., Liu F.K., 2011, ApJ, 729, 13
26. Côté P. et al., 2006, ApJ Suppl. Ser., 165, 57
27. Deane R.P. et al., 2014, Nature, 511, 57
28. Dubois Y., Volonteri M., Silk J., 2014, MNRAS, 440, 1590
29. Fabbiano G., Wang J., Elvis M., Risaliti G., 2011, Nature, 477, 431
30. Fujita Y., 2009, ApJ, 691, 1050
31. Ginsburg I., Loeb A., 2006, MNRAS, 368, 221
32. Gnedin O.Y., Gould A., Miralda-EscudÃ© J., Zentner A.R., 2005, ApJ, 634, 344.
33. Goicovic F.G. et al., 2015, preprint (arXiv:1507.05596)
34. Gould A., Quillen A.C., 2003, ApJ, 592, 935
35. Gualandris A., Portegies Zwart S.F., 2007, MNRAS Lett., 376.1, L29
36. Gualandris A., Portegies Zwart S.F., Sipior M.S., 2005, MNRAS, 363, 223
37. Gvaramadze V.V., 2009, MNRAS, 395, L85
38. Gvaramadze V.V., Gualandris A., 2011, MNRAS, 410, 304
39. Gvaramadze V.V., Gualandris A., Portegies Zwart S.F., 2009, MNRAS, 396, 570
40. Hellström C., Mikkola S., 2010, Celest. Mech. Dyn. Astron., 106, 143
41. Hernquist L., 1990, ApJ, 356, 359
42. Hills J.G., 1988, Nature, 331, 687
43. Hoogerwerf R., de Bruijne J.H.J., de Zeeuw P.T., 2001, A & A, 365, 49
44. Humason M.L., Zwicky F., 1947, ApJ, 105, 85
45. Khan F.M., Holley-Bockelmann K., Berczik P., Just A., 2013, ApJ, 773.2, 100
46. Kobayashi S., Hainick Y., Sari R., Rossi E.M., 2012, ApJ, 670, 747
47. Komossa S. et al., 2003, ApJ, 582, L15
48. Kormendy J., Ho L.C., 2013, Annu. Rev. Astron. Astrophys., 51, 511
49. Kornreich D.A., Lovelace R.V.E., 2008, ApJ, 681, 104
50. Leonard P.J.T., Duncan M.J., 1990, AJ, 99, 608
51. Li Y. et al., 2015, preprint (arXiv:1506.01818v2)
52. Lingam M., 2014, Astrophys. Space Sci., 354, 561
53. Liu F.K., Chen X., 2013, ApJ, 767, 18
54. Liu F.K., Li S., Chen X., 2009, ApJ, 706, L133
55. Liu F.K., Li S., Komossa K., 2014, ApJ, 786, 103
56. Loose H.H., Kruegel E., Tutukov A., 1982, A& A, 105, 342
57. Marconi A., Hunt L.K., 2003, ApJ Lett. 589, L21
58. Matthews L.D., Gallagher J.S., 1997, AJ, 114, 1899
59. Mayer L. at al., 2007, Science, 316, 1874
60. Merritt D., Mikkola S., Szell A., 2007, ApJ, 671, 53
61. Merritt D., Milosavljević M., 2005, Living Rev. Relativ., 8, 8
62. Mikkola S., Aarseth S., 2001, Celest. Mech. Dyn. Astron., 84, 343
63. Mikkola S., Merritt D., 2006, MNRAS, 372, 219
64. Mikkola S., Merritt D., 2008, AJ, 135, 2398
65. Mikkola S., Tanikawa K., 1999a, Celest. Mech. Dyn. Astron., 74, 287
66. Mikkola S., Tanikawa K., 1999b, MNRAS, 310, 745
67. Milosavljević M., 2004, ApJ Lett., 605, 13 MilosavljeviÂ´
68. Miyamoto M., Nagai R., 1975, Publ. Astron. Soc. Jpn., 27, 533
69. O’Leary R. M., Loeb A., 2008, MNRAS, 383, 86
70. Palladino L.E. et al., 2014, ApJ, 780, 7
71. Perets H.B., 2009, ApJ, 698, 1330
72. Perets H.B., Hopman C., Alexander T., 2007, ApJ, 656, 709
73. Perets H.B., Subr L., 2012, ApJ, 751, 133
74. Peters P.C., 1964, Phys. Rev. B, 136, 1224
75. Plummer H.C., 1911, MNRAS, 71, 140
76. Portegies Zwart S.F., 2000, ApJ, 544, 437
77. Poveda A., Ruiz J., Allen C., 1967, Bol. Obser. Tonantzintla y Tacubaya, 4, 86
78. Przybilla N., Nieva M.F., Heber U., Butler K., 2008, ApJ Lett., 684, L103
79. Quinlan G.D., 1996, New Astron., 1, 35
80. Rodriguez C. et al., 2006, ApJ, 646, 49
81. Rosado P.A., Sesana A., 2014, MNRAS, 439, 3986
82. Roškar R. at al., 2015, MNRAS, 449, 494
83. Rossi E.M., Kobayashi S., Sari R., 2014, ApJ, 795.2, 125
84. Rownd B.K., Dickey J.M., Helou G. 1994, AJ, 108, 1638
85. Sari R., Kobayashi S., Rossi E.M., 2010, ApJ, 708, 605
86. Scheck L., Kifonidis K., Janka H.-T., Müller E., 2006, A & A, 457, 963
87. Schinnerer E. et al., 2006, ApJ, 649, 181
88. Schinnerer E. et al., 2008, ApJ Lett., 684, 21
89. Sesana A., 2013, Class. Quant. Grav., 30.24, 244009
90. Sesana A. et al., 2005, ApJ, 623.1, 23
91. Sesana A., Barausse E., Dotti M., Rossi E.M., 2014, ApJ, 794.2, 104
92. Sesana A., Haardt F., Madau P., 2006, ApJ, 651.1, 392
93. Sesana A., Haardt F., Madau P., 2007, MNRAS Lett., 379, L45
94. Sherwin B., Loeb A., O’Leary R., 2008, MNRAS, 386, 1179
95. Silva M.D.V., Napiwotzki R., 2011, MNRAS, 411, 2596
96. Tremaine S.D., Ostriker J.P., Spitzer L.J., 1975, ApJ, 196, 407
97. Turner M.L. et al., 2012, ApJ Suppl. Ser., 203, 5
98. Tutukov A.V., Fedorova A.V., 2009, Astron. Rep., 53.9, 839
99. Vickers J.J., Smith M.C., Grebel E.K., 2015, AJ, 150.3, 77
100. Volonteri M., Haardt F., Madau P., 2003, ApJ, 582, 559
101. Wang J.C.L., Sulkanen M.E., Lovelace R.V.E., 1992, ApJ, 390, 46
102. Wegg C. & Nate Bode J., 2011, ApJ, 738, L8
103. Yu Q., 2002, MNRAS, 331, 935
104. Yu Q., Madau P., 2007, MNRAS, 379, 1293
105. Yu Q., Tremaine S., 2003, ApJ, 599, 1129
106. Zhong J. et al., 2014, ApJ, 789, L2
107. Zubovas K., Wynn G.A., Gualandris A., 2013, ApJ, 771.2, 118
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters