# Interaction of Recoiling Supermassive Black Holes with Stars in Galactic Nuclei

## Abstract

Supermassive black hole binaries (SMBHBs) are the products of frequent galaxy mergers. The coalescence of the SMBHBs is a distinct source of gravitational wave (GW) radiation. The detections of the strong GW radiation and their possible electromagnetic counterparts are essential. Numerical relativity suggests that the post-merger supermassive black hole (SMBH) gets a kick velocity up to due to the anisotropic GW radiations. Here we investigate the dynamical co-evolution and interaction of the recoiling SMBHs and their galactic stellar environments with one million direct -body simulations including the stellar tidal disruption by the recoiling SMBHs. Our results show that the accretion of disrupted stars does not significantly affect the SMBH dynamical evolution. We investigate the stellar tidal disruption rates as a function of the dynamical evolution of oscillating SMBHs in the galactic nuclei. Our simulations show that most of stellar tidal disruptions are contributed by the unbound stars and occur when the oscillating SMBHs pass through the galactic center. The averaged disruption rate is , which is about an order of magnitude lower than that by a stationary SMBH at similar galactic nuclei. Our results also show that a bound star cluster is around the oscillating SMBH of about the black hole mass. In addition, we discover a massive cloud of unbound stars following the oscillating SMBH. We also investigate the dependence of the results on the SMBH masses and density slopes of the galactic nuclei.

## 1 Introduction

Supermassive black hole binaries (SMBHBs) are predicted by the hierarchical galaxy formation model in cold dark matter (CDM) cosmology (Begelman et al., 1980; Volonteri et al., 2003). For merging galaxies, their two SMBHs with galactic cores will firstly approach each other by dynamical friction, and then get close enough to form a bound, compact binary system. After that, the SMBHB may stall at a so called “hard binary separation” for a time even longer than the Hubble time (Begelman et al., 1980). However, recent investigations suggest that the hardening rates of SMBHBs can be boosted and they may coalesce within the Hubble time either due to various stellar dynamical processes other than spherical two-body relaxation (Yu, 2002; Chatterjee et al., 2003; Merritt & Poon, 2004; Berczik et al., 2006; Sesana et al., 2008), or gas dynamics (Gould & Rix, 2000; Colpi & Dotti, 2009, and references therein).

Since the evolution of SMBHBs deeply impacts the evolution of host galaxies, it is very important for us to find observational evidences to constrain evolution models of SMBHBs. A statistic way is the calculation of tidal disruption rates. A dormant SMBH could be temporarily activated by tidally disrupting a star passing by and accreting the disrupted stellar debris (Hills, 1975; Rees, 1988; Evans & Kochanek, 1989; Lodato et al., 2009), which may have been observed in several non-active galaxies (Komossa & Bade, 1999; Komossa, 2002; Gezari et al., 2008, 2009). Chen et al. (2008) and Chen et al. (2009) calculated the tidal disruption rate in SMBHB systems at different evolutionary stages, and found that is significantly different from the typical rate in a single SMBH for several orders of magnitude.

The most certain proof for detecting SMBHB individually may come from gravitational wave (GW) radiation observation. Coalescing SMBHBs are important sources of GW radiation (Peters, 1964; Begelman et al., 1980), and can be detected within the next decades by the Laser Interferometer Space Antenna (LISA) and the Pulsar Timing Array (PTA) program (Berentzen et al., 2009; Amaro-Seoane et al., 2010). Because of the poor spatial resolution of both LISA and PTA for locating GW radiation sources, it is of key importance to detect electromagnetic counterparts (EMCs) of GW radiation sources. Besides, identifying SMBHBs by their EMCs is also essential to constrain the poorly understood galaxy-merger history.

Several EMCs have been suggested in the literature to probe SMBHBs, (1) precession of jet orientation and its acceleration in radio galaxies during the in-spiraling of SMBHBs (Begelman et al., 1980; Liu & Chen, 2007), (2) optical periodic outbursts in blazars due to the interaction of SMBHBs and accretion disk (Sillanpaa et al., 1988; Liu et al., 1995, 2006; Liu & Wu, 2002; Valtonen et al., 2008; Haiman et al., 2009), (3) jet reorientation in X-shaped radio galaxies due to the exchange of angular momentum between SMBHBs and accretion disk (Liu, 2004), (4) interruption of tidal disruption flares in SMBHBs systems (Liu et al., 2009). Besides, there are also some EMCs to probe the coalescence of SMBHBs and its remnant, (1) intermittent activity in double-double radio galaxies at binary coalescence (Liu et al., 2003), (2) X-ray, UV, optical, and IR afterglow following binary coalescence (Milosavljević & Phinney, 2005; Shields & Bonning, 2008; Lippai et al., 2008; Schnittman & Krolik, 2008), (3) systematically shifted broad emission lines relative to narrow emission lines (Komossa et al., 2008) and off-center active galactic nuclei (AGNs) (Madau & Quataert, 2004; Loeb, 2007) because of SMBH GW radiation recoil, (4)tidal disruption flares and hypercompact stellar systems for recoiling black holes (Komossa & Merritt, 2008; O’Leary & Loeb, 2009; Merritt et al., 2009; Stone & Loeb, 2011).

The breakthrough on numerical relativity in the past few years reveals the final stage of SMBHB’s coalescence (Pretorius, 2005; Campanelli et al., 2006; Baker et al., 2006a). The coalescence remnant SMBH can be recoiled due to the anisotropically GW emission during the inspiral and final coalescence (Peres, 1962). For nonspinning SMBHs, the recoil velocity is (Baker et al., 2006b; González et al., 2007a; Herrmann et al., 2007), which is just as the same order of magnitude as the stellar dispersion velocity in the galactic center. However, if both of two SMBHs are rapidly spinning, the recoil velocity can be as large as thousands kilometers per second. This recoil velocity depends sensitively on the intersection angle between spin vectors of two SMBHs and their linear momenta. For some extreme cases, with maximally spinning equal-mass SMBHs and antialigned spins in orbital plane, the recoiling velocity can achieve to (Campanelli et al., 2007; González et al., 2007b).

A recoiling SMBH with high velocity has significant displacement relative to the galactic center. For most of massive galaxies with escape velocity (Merritt et al., 2004), the recoiling SMBH has opportunity to escape from its host galaxy. Those recoiling SMBHs which are still bound to host galaxies will oscillate around the galactic centers and change the core density profiles of the host galaxies due to dynamical friction (Boylan-Kolchin et al., 2004; Gualandris & Merritt, 2008). In addition, the recoiling SMBH with a fraction of its accretion disk can produce sets of emission lines separated by relative velocity (Komossa et al., 2008) or off-center active galactic nuclei (AGNs) (Madau & Quataert, 2004; Loeb, 2007).

In gas poor environment, a recoiling SMBH might be observed as a “hypercompact stellar system” (HCSS). That is a kind of star cluster which is bound to the recoiling SMBH, with similar luminosity as a globular cluster or ultracompact dwarf galaxy, and very high velocity dispersion (Merritt et al., 2009; O’Leary & Loeb, 2009). The other signature is an offset tidal flare due to the recoiling SMBH tidally disrupting surrounding stars. Komossa & Merritt (2008) have used both analysis and -body simulations to calculate the tidal disruption rates of the recoiling SMBH for bound and unbound stars. They found that the tidal disruption rates, in most cases, are smaller than in stationary systems. Recently, O’Leary & Loeb (2011) have investigated the bound stellar density profile and tidal disruption evolution around the recoiling SMBH. They found that the tidal disruption rate will monotonically fall as a power-law after the recoil. There is also an analytical work focusing on the tidal disruption rate during a short period just after coalescence (Stone & Loeb, 2011), which predicts multiple tidal disruption flares in few years or decades after the coalescence.

Most of the above works about tidal disruption of the recoiling SMBH have focused on theoretical analysis and neglected the impact of background stars. Since the real evolution process in such kind of system is very complex, a detailed study on the recoiling SMBH co-evolving with the galactic center is essential. Unlike previous work, we evolve the entire system before and after the recoil with unbound stars included. Our work focuses on investigating the co-evolution between the recoiling SMBH and surrounding stars with tidal disruption and accretion processes. To it come true, we use a special high-accuracy, parallel direct -body code accelerated by many-core hardware ( -Grape and -GPU) (Berczik et al., 2005; Harfst et al., 2007; Spurzem et al., 2009; Just et al., 2011), including a simplified tidal disruption scheme (Fiestas et al., 2011). Most of our simulations are calculated on the GPU cluster in National Astronomical Observatories of China (NAOC). With up to million of particles in simulation, we can find out the dynamical evolution of stars near the recoiling SMBHs and estimate their tidal disruption rates, which may give us some hints for observation.

Our results show that the accretion of disrupted stars does not significantly impact the dynamical evolution of recoiling SMBH. However, it changes the tidal disruption rate compared to a stationary SMBH system in a galactic center. Most of tidal disruption events are contributed by unbound stars when the SMBH passing through the galactic center, and the disruption rate in average is , which is roughly an order of magnitude lower than the stationary case, . Besides, the tidal disruption rate for oscillating SMBH far from the galactic center is roughly one or two orders of magnitude lower than stationary SMBH. We also find a bound stellar system around the recoiling SMBH with a mass of black hole mass, which is consistent with Merritt et al. (2009). Except for bound stellar systems, we find that there is a “polarization cloud” composed of stars dynamically perturbed by the oscillating SMBH. This “polarization cloud” is axisymmetrically diffuse in a large spatial area which is beyond the gravitational influence region of the oscillating SMBH. For this reason, most of the stars in the cloud are not permanently bound to the SMBH. This kind of clouds are the echo of dynamical friction. The calculations for parameter dependence show that the variation of initial mass of the recoiling SMBH can speed up or slow down the dynamical evolution, and impact the tidal disruption rate. In addition, the variation of density slope can obviously change the dynamical evolution and tidal disruption rate of the recoiling SMBH.

We give our galactic model and simulation method in Section 2. Some results about bound/unbound stellar systems and tidal disruption rates for the stationary/recoiling SMBH are showed in Section 3. The dependence of our results on other parameters are also investigated. Section 4 gives discussion about the observational implications and a short summary.

## 2 Galactic Model and -Body Simulations

### 2.1 Galactic Model

For simplicity, we adopt a spherical Dehnen model to describe the ellipticals or bulges of galaxies (Dehnen, 1993), which is different from a Sérsic model used in Gualandris & Merritt (2008). In a Dehnen model, the space density profile follows

(1) |

where is the scaling radius, is the total mass of galaxies, and the slope index is a constant between interval [0,3). The mass distribution is proportional to for and for . In this model, it is easy for us to derive the galactic potential, cumulative mass, and half-mass radius respectively (Dehnen, 1993),

(2) |

(3) |

(4) |

where is the gravitational constant. From Equation (2) the escape velocity at radius for is

(5) | |||||

Thus the escape velocity from the galactic center for is

(6) |

If there is a SMBH with mass
at the galactic center, we can estimate its influence radius by using the definition with stellar mass
as in Merritt (2006)^{11}

(7) |

For simplicity, we adopt model units thereafter. In this new unit, the quantities relate to physical quantities with the scale factor of time, velocity, and length, respectively

(8) | |||||

(9) | |||||

(10) |

(11) |

(12) |

(13) |

For a fiducial model adopted in the following calculations with , , and , we can derive and .

### 2.2 Tidal Disruption Scheme in -Body Simulation

Compared to a scattering experiment, direct -body simulation can well represent the dynamical co-evolution of the recoiling SMBH and surrounding stars, especially for the stellar interactions. However, it can not deal with other process well, for example, the tidal disruption by SMBH. To solve this problem, we implement a special tidal disruption scheme in our simulation.

A star with mass and radius will be tidally disrupted if it approaches a BH within the tidal radius (Hills, 1975; Rees, 1988)

(14) |

Here is a dimensionless parameter of order unity. It is roughly for the real conditions in galactic nuclei (Merritt, 2006). Based on Equation (14), to make sure that the tidal radius is larger than the event horizon of the BH, there is a limitation for . For a disrupted star with solar radius and mass , a Schwarzschild BH should have . This limitation mass can be heavier when the BH has rapid spin.

In order to investigate the tidal disruption of the recoiling SMBH with -body simulation, we have to find a proper way to represent the tidal radius in our simulations. As shown in Equation (14), the tidal radius is proportional to the radius of the disrupted star. However, in our -body simulations without stellar evolution scheme, stars are point like particles, whose radii are not defined. We can use the approximate scaling relation between and discussed above to set a proper value for , but this value is so small that there will be only few tidal disruption events recorded in the simulation. That is because of the limitation on particle resolution. Our simulation particle number is not enough even with .

To solve this problem, a simplified strategy is adopted. We assume a larger in our -body simulation to represent the tidal radius. The particle which get close to SMBH within will be tidally disrupted and removed from the system with its mass added onto the SMBH. With this method, we can directly follow the growth of and calculate the tidal disruption rate. However, the large we adopted here will overestimate the tidal disruption rate. The solution is to carry out a series of simulations with decreasing and extrapolate the results to the regime corresponding to real galaxies. Our extrapolation scheme is explained in Section 3.3.1 in detail. We do not consider the growth of for accreting SMBH or the momentum transfer from tidal disrupted stars to the recoiling SMBH, because the total mass of disrupted stars is negligible compared to the BH. A test simulation taking into account momentum transfer shows no significant difference in either the dynamical evolution of recoiling SMBH or the stellar disruption rate.

### 2.3 Numerical Method

To investigate the evolution of the recoiling SMBH with tidal disruption carefully, we make a series of integrations with different parameters which are listed in Table 1. The first column is the sequence number for different models. Columns (2) - (4) give particle number , tidal radius , and initial stellar density slope respectively. Column (5) is the initial mass ratio between and the total mass of system. Column (6) gives initial recoil velocity in the unit of escape velocity at the galactic center.

All of our integrations adopt a parallel direct -body -Grape/ -GPU code with fourth-order Hermite integrator and simplified tidal disruption scheme. They all have time-step accuracy parameter and a softening length . Most of our models adopt equal mass particles, with initial mass of SMBH to represent the ellipticals/bulges and central SMBH respectively. To find out the dependency of the results on particle number, several sets of integrations with and particles have been included. There is also a calculation with to check the dependency on . As mentioned before, for the purposes of extrapolating to the smaller tidal radius in the real galactic center conditions, we have varied the from to , which is roughly corresponding to to for . Since the density slope prior to the recoil is under debate, we choose as our fiducial value, and also investigate other cases with and . The recoil velocity is set to for most of the integrations, and there are calculations without recoiling velocity for comparison. Besides, an integration with and is included for comparing with Komossa & Merritt (2008).

Since Dehnen model provides the stellar distribution in a galaxy without central SMBH, we have to first to set up consistent initial conditions for our system. At first, we put a SMBH into the center of a stellar system with Dehnen profile, and evolve the entire system for 50 simulation time units to relax the core region. After that the system is dynamically and thermally relaxed in relation to the SMBH gravitational potential. At this stage, we do not consider stellar tidal disruption process. However, we want to make sure that stars inside the tidal disruption loss cone of the central SMBH have already been removed by tidal accretion before the recoil of SMBH. Therefore, as a second step, we switch on the tidal disruption scheme in the program and run it again with the SMBH still in equilibrium in the center - this costs a few dynamical timescales at the influence radius. Thus all particles within the tidal disruption loss cone can be removed before we apply a kick to the SMBH. After setting up the initial configuration of the stellar distribution of the systems, we artificially give a recoil velocity to the SMBH and follow its dynamical evolution.

## 3 Results

### 3.1 Dynamical Evolution of the recoiling SMBH

To carefully investigate the dynamical evolution of recoiling SMBH, we have a long time integration with model to , as shown in Figure 1. It is clear that the evolution of the recoiling SMBH can be easily divided into two different phases. The phase I has obviously damping trajectory from to . After that, the SMBH evolves into phase II, where the trajectory damped very slow for a long time. In phase I, the trajectory of recoiling SMBH can be well predicted by Chandrasekhar’s dynamical friction theory (Chandrasekhar, 1943). However, this stage will just keep till the SMBH’s oscillation decayed to the size of core, where a slowly decayed phase II begins. That is because the stellar mass interior to SMBH’s orbit is roughly equal to at the end of phase I, where the assumptions of Chandrasekhar’s dynamical friction theory are invalid. We have tried to fit the trajectory analytically in phase I with the same method as Gualandris & Merritt (2008) used, and got a similar result even though a Dehnen model is adopted here instead of their core-Sérsic model.

If we continue the integration, the orbital oscillation of the SMBH will slowly decay and finally achieve to a thermal equilibrium with surrounding stars, which is the so-called phase III in Gualandris & Merritt (2008). Here we have not continue our integrations into that stage because phase III is actually similar to the case of a stationary central SMBH, which is out of the scope of this paper. We will just focus on phase I and phase II in the following integrations. Thus all of our models integrate to in order to make sure that the evolution can achieve to phase II.

### 3.2 Compact Star Clusters around the Recoiling SMBHs

A recoiling SMBH can carry a retinue of bound stars. Merritt et al. (2009) predict that there will be a “hypercompact stellar system” (HCSS) bound to the recoiling SMBH. We can analytically estimate the bound population of the HCSS with Equations (1a) and (6a) in Merritt et al. (2009). For example, the roughly boundary radius of HCSS should be

(15) |

with our Model , where , we have in our simulation unit. Thus we can define a factor to describe the fraction of bound stars,

(16) |

where is total mass of bound stars. Based on Equations (15) and (13), we have . For the case of , with in Figure 1a of Merritt et al. (2009), we have . That means the ratio for HCSS mass to should be , corresponding to just few particles for simulation. It should be aware that this fraction is just a lower limit because our “two times ” assumption actually overestimate . To confirm the conclusion above, we investigate this problem with -body simulations.

In the simulation, it is very difficult to distinguish whether a star is bound to SMBH or not. That is because most of star particles with negative total energies relative to the SMBH are loosely bound. Only few of them could orbit the SMBH for more than one orbital period. For this reason, we identify a star to be bound to the recoiling SMBH if following two criteria are both satisfied: the particle (1) has negative total energy relative to the SMBH and (2) remains bound for at least one orbital oscillation period of the recoiling SMBH. We find that most of the “bound” particles around SMBH particle will be evaporated before the end of simulation. This effect is more significant for our escaping SMBH calculation with model 10, where the number of bound particles decreases monotonically. For other oscillating SMBH calculations, there are many unbound particles can be captured by the SMBH. That means the interaction between the HCSS and stellar background is very important.

Figure 2 shows the evolution of bound stars around a recoiling SMBH in model . At , when , there are five bound particles, roughly consistent with the estimation from Equation (16). The number of bound particles, , slightly increases during the evolution from phase I to phase II. Because in phase II (1) the low relative velocity between stars and the recoiling SMBH increases the probability of dynamical capture and (2) the stellar density around the SMBH in phase II is higher. Our results suggest that a HCSS may grow via capturing unbound single stars. However, the low particle resolution and large statistical error prevent us from further investigation of the growth of HCSS.

The HCSS may be detected as an off-nuclear compact system of similar size and stellar mass of a globular clusters but of having very high internal velocities (Merritt et al., 2009; O’Leary & Loeb, 2011). However, that is not the whole story. As shown in the Figures 3 and 4, we find that there is a stellar cloud composed of stars which are strongly impacted by the recoiling SMBH at the first apocenter. The total mass of this cloud is comparable to . They form a quasi-axisymmetric stellar cloud and always have a dense region around SMBH. This kind of clouds gives an evidence for the echo of dynamical friction. For this reason, we name that stellar cloud as “polarization cloud”. In our simulation of model 06, we choose a time snapshot when the SMBH arrives at its first apocenter, and pick up those polarization cloud members from the phase space which includes all of star particles. All of polarization cloud members in the phase space are outliers around SMBH particle and can be easily distinguished. Thus we can collect those polarization cloud particles and trace their evolution. However, it is very difficult to find an automatical solution to efficiently distinguish those polarization cloud members for every time snapshot. To deeply investigate this kind of polarization clouds, an efficient method to extract those outlier particles is needed. This problem should be solved in our following work.

In Figure 3, we plot that polarization cloud at time , when the recoiling SMBH arrives at its first apocenter. All the parameters here are obtained from model . In order to investigate the origin and evolution of those cloud members, we look back the data and follow the evolution of the stars from time . Figure 4 gives the evolution of those stars and the recoiling SMBH from until , when our computation stopped. The snapshots with in Figure 4 show the evolutions of the cloud members from beginning toward the first apocenter of the recoiling SMBH. The following three snapshots correspond to the time when the recoiling SMBH comes back from the first apocenter and passes through the density center for the second time. The snapshot at shows the SMBH passing through galactic center for the third time, and the last snapshot is for the end of computation. Our calculation results and Figure 4 imply that the shape of polarization cloud varies with the oscillations of the recoiling SMBH. The comparison between the first and the last two panels vividly shows the changes of the cloud members during the oscillation.

As shown in Figure 2, there are only small fraction of stars really bound to the SMBH, whereas the mass of the polarization cloud is roughly equal to . Therefore most of these cloud members are unbound to the recoiling SMBH. After the recoil of SMBH, those polarization cloud members are accelerated and fall behind the oscillating SMBH. During the oscillation, some of the stars gradually expand to larger area while others continuously follow the SMBH to oscillate around the galactic center, and form an axisymmetrical distribution along the velocity direction of the recoiling SMBH. After multiple interactions with oscillating SMBH, these stars obtain energy from oscillating SMBH and diffuse to a large area. This is a vivid example that an oscillating SMBH transfer its orbital energy to surround stars through dynamical friction.

In order to investigate the observational properties of this polarization cloud, we focus on its most compact stage, when the recoiling SMBH arrive at the apocenter for the first time. The distribution of those stars are shown in Figure 3. Here we scale the -body system to our physical fiducial galaxy model with stellar mass and half mass radius . For such system, the polarization cloud has mass of , and a size of diameter . Further calculation shows that, its half mass radius is . This size and mass are much larger than a globular cluster and comparable to an ultracompact dwarf galaxy. We also calculate the velocity dispersion of this polarization cloud. The line-of-sight velocity dispersion along X, Y and Z axes are , and individually. That is significantly higher than a typical ultracompact dwarf galaxy. Comparing with bound HCSS, this polarization cloud has larger size, heavier mass and similar velocity dispersion. Since this cloud is mixed with HCSS, it may bring troubles to the detection of HCSS. We also estimate the average surface density inside the half mass radius of the polarization cloud. The overdensity of the polarization cloud relative to the stellar background is only 4%. This may be because of that the apocenter of the recoiling SMBH here is not very far from the galactic center, where the surface density of the stellar background is still relatively high in model 06 with . For a SMBH with a higher recoil velocity and a larger apocenter, the situation may be improved. To investigate this problem, more simulations are needed. We will investigate this problem in our future work.

It should be noted that our example above is just for one time snapshot of the evolving SMBH. Every snapshot could have a similar stellar cloud. The composition of these clouds may form an interesting structure which can trace the trajectory of the SMBH. The evolution and observation characters of this kind of structures will be discussed in our future work. For a recoiling SMBH with highly enough velocity to escape from the galaxy, that kind of polarization cloud may not be distinguishable because many stars will be dropped behind the escaping SMBH. Without the multiple interactions between oscillating SMBH and surrounding stars, the polarization cloud can not survive for a long time.

### 3.3 Tidal Disruption

#### Calculate Tidal Disruption Rate Numerically

For a solar type star disrupted by a SMBH with mass , based on Equation (14) and (10), the tidal radius is , corresponding to in simulation units with our fiducial model. As mentioned in Section 2.1, because of the limitation on particle resolution, it is practically impossible with present computing capabilities to calculate the tidal disruption rate directly (with a realistically small tidal radius) through direct -body simulations.

To solve this problem, we perform a parameterized study, which treats both the particle number and the size of the tidal disruption radius as free parameters. All of the simulations have been done with these varying parameters, and physical conclusions are drawn from extrapolating to the real galactic conditions. Fiestas & Spurzem (2010) have shown in Fokker-Planck models that such scaling is reliable and provides physically correct results. Besides, tidal disruption events in our simulations, as well as the case in real galactic center, do not occur frequently, which means that the increase of is discrete. Thus it is difficult to compute the tidal disruption rate through calculating accreted masses in time bins. Instead, we measure the tidal disruption rate in a cumulative, averaged way - fitting the mass growth of the black hole and derive its time derivative.

In order to find out the extrapolation relation, the dependence of tidal disruption rate on particle numbers and tidal radius should be investigated. Figure 5 shows the disruption rates changes with different and , with parameters are , , , and . The left panel shows that the dependence of averaged tidal disruption rates on and . Here we calculate the averaged disruption rates for phase I, phase II and entire process separately with different and . These integrations show that the dependence of tidal disruption rate on particle number is very weak. Further discussions are showed in Section 3.4.

Based on the loss cone feeding theory (Frank & Rees, 1976; Lightman & Shapiro, 1977), contrary to empty loss cone case, the tidal disruption rate for full loss cone does not depend on . It should be notice that the tidal disruption rate here is for mass disruption rate, and the disrupted particles indeed depend on . The slight different average rates between and cases means that the loss cones in our integrations are nearly full and the changes of particle numbers do not bring significant impact.

Opposed to particle numbers, the changes on strongly influence the disruption rates. A detailed study is in the right panel of Figure 5, which shows the evolution of disruption rate - dependence for several special points with integrations for . Five special snapshot have been picked out. Here “1st. APO.”, “1st. BACK”, “1st. DOWN”, “2nd APO.” and “Phase Tran.” refer respectively to the snapshot that the oscillating SMBH the first time arrives at the apocenter, the first time returns to the density center, the first time reaches the opposite apocenter, the second time return to the apocenter and the time of transition from phase I to phase II. Our results show that the disruption rates increase for every value during the evolution of the recoiling SMBH. That is because a decayed orbit of SMBH leads to a growing stellar background around it. Besides, the similar linear disruption rate - relation as left panel has been found. However, the relation here for the first four snapshots are not as good as the average results in the left panel. That may because our particle resolution for small is not good enough in phase I.

Based on our result data, it seems that there is an approximate linear relation between and disruption rates. As a result of full loss cone condition, the disruption rates will be proportional to cross section and density near the recoiling SMBH,

(17) |

where is average relative velocity between the recoiling SMBH and the surrounding stars.

It is difficult for us to analytically derive stellar distribution around a recoiling SMBH. However, the cross section can be wrote as

(18) |

which can write the tidal disruption rate as

(19) |

For that is (even if one takes into account that contains a contribution from the SMBH kick velocity), so we are in the gravitational focusing regime, where

(20) |

Thus we have

(21) |

If and are constant with time, the tidal accretion rate will scale linearly with . This relation is clearly only correct in a time averaged sense, because when the SMBH oscillating in the galaxy, all quantities will be strongly changed within a single orbit (density, stellar velocity dispersion and SMBH velocity, the two latter quantities determining ). However, in a running time averaged over several orbits, our results show that there is not a large variation in these quantities. So we can define a dimensionless tidal accretion parameter by the following relation:

(22) |

Here is the dynamic time scale around the position of the recoiling SMBH. Our simulations imply that across all of our model families does not depend strongly on or , but only on the evolutionary phase. Thus can be used in extrapolations to the real system with very small and large . Based on the analysis above, though it is complicated to get an accurate scaling relation, we still can estimate the tidal disruption rate through linear extrapolation of tidal radius.

#### Tidal Disruption Rate for the Stationary SMBH

For comparison purposes, we calculate the tidal disruption rate of a stationary SMBH without recoil in the galactic center. A series of simulations with different values are carried out. Figure 6 shows the static disruption rate for , corresponding to model , , and respectively. As mentioned in Section 3.3.1, the tidal radius for our fiducial model is , corresponding to in our simulation units. Our simulation results with different particle numbers for stationary case show that the tidal disruption rates weakly depend on the particle numbers. Thus we can loosely accept a full loss cone condition for the stationary case. As discussed in Section 3.3.1, there is a roughly linear correlation between tidal disruption rate and for full loss cone case. Based on Equation (21), we can extrapolate our simulation results to the real galaxy conditions with smaller tidal radii. Figure 6 shows that the tidal disruption rate for is in our simulation units. That can be extrapolated to , and gives us the tidal disruption rate , which corresponds to the order of magnitude for our fiducial model.

#### Tidal Disruption Rate for the recoiling SMBH

In a recoiling SMBH system, there are two populations of stars contributed to the tidal disruption rate. One is the bound stars around the SMBH, and the other is those unbound stars encountered with the recoiling SMBH. If the bound stars dominate, as demonstrated in Komossa & Merritt (2008) for an escaping SMBH, the tidal disruption rate should roughly keep a constant level at the beginning, and finally drop down because of the consumption of bound stars. On the contrary, if the tidal disruption rate is dominated by unbound stars, the results will depend on stellar density, relative velocity of SMBH, and the size of tidal radius, which described by Equation (19).

If we just consider about the contribution of the unbound stars, for an oscillating SMBH, its tidal disruption rate near the galactic center should be higher than the apocenter region because of the different stellar density. However, can also change the result. For instance, we have an estimation in our Model with , , and . Our simulation result shows that the first apocenter of the recoiling SMBH is around . With the assumptions that and , we can estimate the disruption rate ratio from the galactic center to the first apocenter region. By substituted into Equations (1) and (19), that ratio is . When scaling to the real galactic center conditions with , it should be . Moreover, this ratio also depends on . If near apocenters can be larger than , that ratio may be higher. Unfortunately, we can not obtain an accurate value for . For estimation, it can be seen that the tidal disruption rates near the galactic center are several times or even an order of magnitude higher than the apocenter region.

That conclusion is confirmed by our simulation results. Figure 7 gives the distribution of the tidal disruption events count relative to distance in model . It shows that most of the tidal disruption events appear around the density center, which indicates a boosted tidal disruption rate when the oscillating SMBH pass through the galactic center. A SMBH around the apocenter with a nearly zero velocity always corresponds to a low , which should give us a high disruption rate. However, as Equation (19) shown, the low stellar density environment makes the tidal disruption events less than at galactic center. The higher disruption rate near the galactic center means that, for an oscillating SMBH, most of tidal disruption events are contributed by unbound stars. For this reason, unbound stars are very important for calculating tidal disruption rate.

Figure 8 shows the tidal disruption rates for a recoiling SMBH with different , which belong to model and respectively. To increase the accuracy of our fitting results for disruption rate in phase I and phase II, we neglect the transition boundary between the two phases in the fitting. Therefore, we fit disruption rate for phase I and phase II separately. As a result, the fitted results are not continuous between two phases.

The left part of Figure 8 shows that the tidal disruption rates in phase I are linearly increasing. For different values, there is also a similar linear relation to stationary SMBH case, even though there is an increasing tidal disruption rate instead of a constant value. Here the calculation with model () is not included because the particle resolution is too low to estimate disruption rate. With the same extrapolating relation described in Section 3.3.1, we estimate that the averaged tidal disruption rate for phase I should around the order of magnitude , which is about an order of magnitude lower than stationary SMBH in the galactic center.

For the phase II, when the orbit of the recoiling SMBH damped to relative small region, the variations of stellar density and velocity of SMBH are smaller than phase I. Thus the difference of disruption rate between apocenters and galactic center is not so intensive. Since the tidal disruption events during entire phase I are few in our simulation, the accuracy for calculating tidal disruption rate is not very good. However, there are plenty of tidal disruption events in phase II, which provides a better accuracy to calculate the tidal disruption rate in phase II. Actually, as shown in the right part of Figure 8, there are constant disruption rates as similar as the stationary case. That is because (1) the stellar densities around SMBHs in phase II and in the stationary SMBH case are similar and (2) their loss cones are both full.

### 3.4 Dependence on Particle Number

Since particles do not have enough high resolution to represent a real galaxy, it is important for us to find out the dependence of our results on the particle number. Only with this considered, the direct -body simulation results for dynamical properties of the recoiling SMBHs and surrounding stars can be extrapolated to real conditions in galactic nuclei. To bring this forth, a series of integrations with various and same other parameters have been done, which relative to models , and with respectively. Figure 9 shows the particle number dependence of the recoiling SMBH orbital oscillation and mass increasing for model , and respectively. It can be seen that the and cases have similar oscillate amplitudes in phase II, whereas the run can not well represent the damped evolution of the recoiling SMBH. Besides, only the calculations with and particles can give relatively smooth mass increasing curve, and they all achieve to the similar final mass.

Limited by the particle resolution, we can not obtain good enough mass evolution data to fit the disruption rates in phase I both for and integrations with model and . The only thing can be done is to calculate average tidal disruption rates for phase I and II, as discussed in Section 3.3 with Figure 5. However, those problems do not exist for all of integrations and other integrations. For this reason, most of our calculations adopt particles. As shown in Figure 5, the difference of average disruption rates between and is very small, which also indicates that the loss cones are nearly full.

With particles simulation, the mass ratio of SMBH to stars should be . This ratio is smaller than the real galaxy condition which should be . The small mass ratio between SMBH and stars may lead to a higher Brownian velocity comparing with real galaxy condition (Merritt et al., 2007). However, most of our simulations are terminated before phase III, which makes the amplitude of oscillating SMBH is greater than Brownian amplitude. Detailed discussion in Gualandris & Merritt (2008) with series of calculations for different particle numbers shows that the particle mass can not significantly impact the dynamical evolution for both phase I and II.

The limitation of particle number can also influnece the relaxation time . For a homogenous isotropic distribution with equal-mass stars, the two-body relaxation time in our simulation units is (Spitzer, 1987)

(23) |

where is Coulomb logarithm and we can estimate it as

(24) |

When is large enough, for instance , the relaxation time will be roughly proportional to . It means that our -body simulation with particles always give us a shorter relaxation time comparing with real galaxy condition. This fast relaxation in -body simulation may bring troubles to our results. If the core region in our simulation has relaxed before the recoiling SMBH return, our results will departure from the real galaxy condition, because the relaxation time in a real galactic core is always very long.

To solve this problem, we calculate the relaxation time of the central core region for both the simulation model and real galaxy condition. We choose , and in model to calculate the relaxation time . Since we know little about the exactly value of velocity dispersion around , we range the from to 1. For the real galaxy condition, our estimations show that, the is always longer than the half period of the oscillating SMBH. That means the core region with a real galaxy condition can not relax before the SMBH return. The same result applied to our simulations with for . Based on Equation (9), corresponds to , which is smaller than most of the galactic center cases. Thus we can conclude that the core region can not relax before the recoiling SMBH return both for simulation and real galaxy condition.

The limitation of particle number in -body simulation may also impacts the results for the stationary SMBH. In recoil case, the oscillating SMBH always corresponds to the full loss cone status both for the simulation and the real galaxy condition. For a stationary SMBH, things will be different. Both the analytical calculations and the numerical simulations for long term evolution of tidal disruption from stationary SMBH show that the loss cone in real galaxy is empty for most of the cases (Frank & Rees, 1976; Lightman & Shapiro, 1977; Cohn & Kulsrud, 1978; Baumgardt et al., 2004). However, our calculations are just focus on a relative short period, which correspond to for our fiducial model. During this short period, the loss cone could be roughly full. This prediction has been confirmed by our simulation results.

Based on discussions above, there are some conclusions: (1) both and integrations are good enough to represent dynamical evolution of the recoiling SMBH, (2) in order to investigate the tidal disruption and co-evolution of the recoiling SMBH with surrounding stars, it is better to use integrations with , (3) the small difference on average disruption rates between and integrations indicates that the loss cones are nearly full, otherwise the results should change obviously with various particle numbers, (4) the particle number we adopt in simulation is less than the real galaxy condition, which can lead to a shorter relaxation time scale. However, it can not significantly impact our results.

### 3.5 Dependence on other Parameters

Based on the Chandrasekhar’s dynamical friction theory (Chandrasekhar, 1943), the dynamical friction time scale is inversely proportional to the mass of the recoiling SMBH. For this reason, a heavier SMBH should has shorter evolution time scale in phase I. In order to investigate the dependence of evolution time scale on black hole mass , we run a test model with . Comparing with model , we find that model gives a shorter oscillation period in phase I, which is roughly two times smaller than that in model . That is consistent with the prediction of dynamical friction theory. Besides, the tidal disruption rate in model is roughly two times higher than that in model for both of the two phases, which is roughly consistent with Equation (21). For a galaxy with lighter SMBH, which has , there will be a longer dynamical evolution time in phase I and smaller tidal disruption rate comparing with our results above. Limited by the calculation capability, simulations with the lower mass SMBH will cost so much time that we can not investigate them carefully.

We note that in the simulations with very large tidal radius, for example, in model and , rapid increase on lead to the declined disruption rates. This significant growth of is artificial, which makes the declination of disruption rate is unphysical.

When considering different density slopes, the results are quite different. The higher of value is, the faster of the system evolved. Our integration with model shows that the recoiling SMBH will enter phase II around , which is more faster than evolution of SMBH in model with . Thus our models for and actually reach phase III at the end of integration. That can be easily understood because a larger means a denser cusp in the galactic center. For the same reason, a larger also leads to a higher disruption rate.

Figure 10 shows the disruption rates with different values for both the stationary and the recoiling cases. The left panel is for stationary integrations and the right panel is for the recoiling SMBH with time duration . Here we can not compare the same phase for different models because of their different evolving times. The results in the right panel are phase II for model and phase III for model and . It can be seen that, for both panels, there are several times differences between every two adjacent values.

In order to compare with the work in the literature with higher , a simulation with similar initial parameters comparing with model has been done in model , where and . Our result shows that there are just five particles tidally disrupted at the beginning of recoil, and following with zero event recorded for the rest of time. Because of the low particle resolution, we can only give a lower limit to the tidal disruption rate. For a galaxy with mass , and , the kick velocity should be . Thus our lower limit for tidal disruption rate is , which is lower than in Figure of Komossa & Merritt (2008) with similar conditions. For the bound stars, we find that there are only stars still bound to SMBH after the recoil, which is roughly as the same order of magnitude as the result in Merritt et al. (2009). The small amount of bound stars, low stellar density outside and poor particle resolution result in the rare disruption events in computation. To solve this problem, a series of simulations with larger particle numbers are needed.

## 4 Discussions and Conclusions

Supermassive binary black holes in galactic nuclei may under certain conditions get very close and coalesce due to strong emission of gravitational waves. Because of the anisotropy of GW emission, the coalescence remnant SMBH experiences a strong recoil velocity depending on spins of the original pair of black holes. In this paper we investigate the interaction and co-evolution of the recoiling SMBHs with their galactic stellar environments, using very large direct -body simulations with a simplified tidal disruption scheme. It is the first numerical investigation for the evolution of the tidal disruption by an oscillating SMBHs with unbound stars included. Moreover, we discover a polarization cloud of stars surrounding the oscillating SMBHs, which is not mentioned in the literature.

As argured by Gualandris & Merritt (2008), the dynamical evolution of the recoiling SMBH can be divided into three phases. During the phase I, the trajectory of the recoiling SMBH can be well predicted by Chandrasekhar’s dynamical friction theory (Chandrasekhar, 1943). The phase II begins when the SMBH’s oscillation decayed to the size of galactic core, and the motion of the SMBH damps slowly. At phase III, the orbital oscillation of SMBH slowly decay to achieve a thermal equilibrium with surrounding stars. In our simulation, we obtain similar phase I and II as Gualandris & Merritt (2008) have got. Besides, our results show that the mass fraction of bound stars relative to SMBH is roughly at the beginning. This value is consistent with Gualandris & Merritt (2008) and Merritt et al. (2009).

In addition to following the dynamical evolution of the recoiling SMBH, we estimate the stellar tidal disruption rates by the recoiling SMBHs using direct -body simulation together with our simple tidal disruption scheme. Our results show that the stellar tidal disruptions do not significantly change the dynamical evolution of the recoiling SMBH. The stellar disruption rate by the recoiling SMBH in phase II is , which is about the same order of magnitude of the stellar disruption rate by the stationary SMBH at a spherical galactic center. However, in phase I, the stellar tidal disruption rate is about an order of magnitude lower.

In our simulation, a spherical stellar distribution is assumed. This simplified model is not accurate to represent a real galaxy. For most of real bulges/ellipticals, their stellar distributions are non-spherical. As argued by Vicari et al. (2007), the orbital decay time scales of the recoiling SMBHs are prolonged in triaxial galaxies comparing with spherical galaxies. That is because most of the recoiling SMBHs in triaxial galaxies can not return directly to their galactic centers, where the influences of dynamical friction are the strongest. For the same reason, a non-spherical dark matter halo can also significantly increase the decay time of a recoiling SMBH (Guedes et al., 2009). Since most of the tidal disruption events happened during the SMBHs passing through the galactic centers, the triaxiality may significantly reduce the tidal disruption rate from unbound stars. To investigate the impacts of triaxiality, a good solution is to have series of simulations with triaxial stellar distributions. However, it is very complicated to generate a stationary triaxial distribution for -body simulation, and the prolonged orbital decay time needs more calculation resources. Here we just follow previous works to investigate the tidal disruptions from the recoiling SMBHs in spherical stellar distribution. Further works with triaxial stellar distribution will be included in our next paper.

Actually, the triaxiality can not change our results significantly because the influence of triaxiality is moderate in galaxy models with shallow central density profiles (Vicari et al., 2007). Thus our calculations with will not significantly departure from the triaxial model. Besides, in most of our simulations, the maximal displacement of the recoiling SMBH relative to the density center is no more than . That keeps the trajectory of the SMBH inside the bulge, where the impact of dark matter halo are very weak comparing with baryonic matter. For this reason, we neglect the influence of the triaxial dark matter halo which is mentioned by Guedes et al. (2009).

Our simulation for escaping SMBH does not have enough particle resolution to estimate disruption rates. That needs much larger -body particle number in the future simulations. Here our integrations are focus on relatively low recoiling velocity. For the assumption that and , recoil velocity in our fiducial model is . This velocity produces the oscillating SMBH an amplitude less than the effective radius of elliptical/bulge. Our results show that most of tidal disruption events occur when the recoiling SMBH passes through density center, but not exact the original center because of the orbital precession. The tidal disruption rates around each apocenter are much lower. An estimate for average tidal disruption rate at about apocenters in phase I over several periods for model has been done. We calculate the average disruption rate when SMBH have large distant relative to density center. The disruption rate at large distance is several times or even an order of magnitude lower than the average value over whole phase I, which implies that the rate for apocenters is one or two order of magnitude lower than that by stationary SMBHs.

Komossa & Merritt (2008) argued that off-nuclear X-ray tidal flares could be one of the observational signatures for an escaping SMBH. Unfortunately, our simulation of model with particles and the escapable recoil velocity does not have sufficient resolution to estimate the disruption rate. We can just give a lower limit , which does not conflict to their work with similar configurations. Besides, our model with and recoil velocity corresponds to a tidal disruption rate . This value is similar to the result of Model 1 with in Komossa & Merritt (2008). Based on our results, an oscillating SMBH with enough large recoil velocity may also contribute to off-nuclear X-ray flares around its apocenters with the low disruption rates.

An HCSS around the recoiling SMBH in our integrations has similar mass as predicted by Merritt et al. (2009). However, we can not give more information about the star clusters limited by particle resolution. In addition, we find a larger cloud of stars around the recoiling SMBHs, which are not permanently bound to the black hole but keep a kinematic correlation to it for more than an orbital time. It can be described as polarization cloud reacting to the gravitational disturbance by the recoiling SMBH, in other words, the echo of dynamical friction. The size and mass of the polarization clouds are comparable to ultracompact dwarf galaxies, except significantly higher velocity dispersion. The existence of such kind of system may bring difficulties for distinguishing the HCSS. More detailed properties of the polarization cloud will be studied in our following work.

Our simulations for different particle numbers show that both and integrations are sufficient in following the dynamical evolution of the recoiling SMBH, but the former do not have enough resolution for investigating the tidal disruption and co-evolution of SMBH with surrounding stars. Different initial masses of the recoiling SMBHs can change the results moderately. We find that a higher galactic density slope can significantly speed up the dynamical evolution of the system, and the tidal disruption rate is also boosted.

(1) | (2) | (3) | (4) | (5) | (6) |
---|---|---|---|---|---|

01 | 0.5 | 0.001 | 0.7 | ||

02 | 0.5 | 0.001 | 0.7 | ||

03 | 0.5 | 0.001 | 0.7 | ||

04 | 0.5 | 0.001 | 0.7 | ||

05 | 0.5 | 0.001 | 0.0 | ||

06 | 0.5 | 0.001 | 0.7 | ||

07 | 0.5 | 0.002 | 0.7 | ||

08 | 1.0 | 0.001 | 0.0 | ||

09 | 1.0 | 0.001 | 0.7 | ||

10 | 1.0 | 0.001 | 1.1 | ||

11 | 1.5 | 0.001 | 0.0 | ||

12 | 1.5 | 0.001 | 0.7 | ||

13 | 0.5 | 0.001 | 0.0 | ||

14 | 0.5 | 0.001 | 0.7 | ||

15 | 0.5 | 0.001 | 0.0 | ||

16 | 0.5 | 0.001 | 0.7 | ||

17 | 0.5 | 0.001 | 0.0 | ||

18 | 0.5 | 0.001 | 0.7 | ||

19 | 0.5 | 0.001 | 0.0 | ||

20 | 0.5 | 0.001 | 0.7 |

Note. – Col.(1): Model sequence number. Col.(2): Particle numbers adopted in calculations. Col.(3): Tidal disruption radius in the units . Col.(4): Density slope . Col.(5): Initial BH mass in the unit of total mass. Col.(6): Initial recoil velocity in the unit of escape velocity. All the simulations set and .

### Footnotes

- affiliation: Astronomy Department, Peking University, 100871 Beijing, China; fkliu@bac.pku.edu.cn, lis@bac.pku.edu.cn, chenx@bac.pku.edu.cn
- affiliation: Astronomisches Rechen-Institut, Zentrum für Astronomie, Universität Heidelberg, Mönchhofstr. 12-14, D-69120 Heidelberg, Germany
- affiliation: Astronomy Department, Peking University, 100871 Beijing, China; fkliu@bac.pku.edu.cn, lis@bac.pku.edu.cn, chenx@bac.pku.edu.cn
- affiliation: National Astronomical Observatories of China, Chinese Academy of Sciences, 20A Datun Lu, Chaoyang District, 100012, Beijing, China
- affiliation: Astronomisches Rechen-Institut, Zentrum für Astronomie, Universität Heidelberg, Mönchhofstr. 12-14, D-69120 Heidelberg, Germany
- affiliation: Main Astronomical Observatory, National Academy of Sciences of Ukraine, 27 Akademika Zabolotnoho St., 03680 Kyiv, Ukraine
- affiliation: Kavli Institute for Astronomy and Astrophysics, Peking University, 100871 Beijing, China
- affiliation: National Astronomical Observatories of China, Chinese Academy of Sciences, 20A Datun Lu, Chaoyang District, 100012, Beijing, China
- affiliation: Astronomisches Rechen-Institut, Zentrum für Astronomie, Universität Heidelberg, Mönchhofstr. 12-14, D-69120 Heidelberg, Germany
- affiliation: Kavli Institute for Astronomy and Astrophysics, Peking University, 100871 Beijing, China
- This kind of definition is validated only for a singular isothermal sphere nucleus. Here we adopt it also for to estimate at an order of magnitude.

### References

- Amaro-Seoane, P., Sesana, A., Hoffman, L., Benacquista, M., Eichhorn, C., Makino, J., & Spurzem, R. 2010, MNRAS, 402, 2308
- Baker, J. G., Centrella, J., Choi, D.-I., Koppitz, M., & van Meter, J. 2006a, Physical Review Letters, 96, 111102
- Baker, J. G., Centrella, J., Choi, D.-I., Koppitz, M., van Meter, J. R., & Miller, M. C. 2006b, ApJ, 653, L93
- Baumgardt, H., Makino, J., & Ebisuzaki, T. 2004, ApJ, 613, 1133
- Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307
- Berczik, P., Merritt, D., & Spurzem, R. 2005, ApJ, 633, 680
- Berczik, P., Merritt, D., Spurzem, R., & Bischof, H.-P. 2006, ApJ, 642, L21
- Berentzen, I., Preto, M., Berczik, P., Merritt, D., & Spurzem, R. 2009, ApJ, 695, 455
- Boylan-Kolchin, M., Ma, C.-P., & Quataert, E. 2004, ApJ, 613, L37
- Campanelli, M., Lousto, C. O., Marronetti, P., & Zlochower, Y. 2006, Physical Review Letters, 96, 111101
- Campanelli, M., Lousto, C., Zlochower, Y., & Merritt, D. 2007, ApJ, 659, L5
- Chandrasekhar, S. 1943, ApJ, 97, 255
- Chatterjee, P., Hernquist, L., & Loeb, A. 2003, ApJ, 592, 32
- Chen, X., Liu, F. K., & Magorrian, J. 2008, ApJ, 676, 54
- Chen, X., Madau, P., Sesana, A., & Liu, F. K. 2009, ApJ, 697, L149
- Cohn, H., & Kulsrud, R. M. 1978, ApJ, 226, 1087
- Colpi, M., & Dotti, M. 2009, arXiv:0906.4339
- Dehnen, W. 1993, MNRAS, 265, 250
- Evans, C. R., & Kochanek, C. S. 1989, ApJ, 346, L13
- Fiestas, J., Porth, O., Berczik, P., & Spurzem, R. 2011, MNRAS, 1694
- Fiestas, J., & Spurzem, R. 2010, MNRAS, 405, 194
- Frank, J., & Rees, M. J. 1976, MNRAS, 176, 633
- Gezari, S., et al. 2008, ApJ, 676, 944
- Gezari, S., et al. 2009, ApJ, 698, 1367
- González, J. A., Hannam, M., Sperhake, U., Brügmann, B., & Husa, S. 2007b, Physical Review Letters, 98, 231101
- González, J. A., Sperhake, U., Brügmann, B., Hannam, M., & Husa, S. 2007a, Physical Review Letters, 98, 091101
- Gould, A., & Rix, H.-W. 2000, ApJ, 532, L29
- Gualandris, A., & Merritt, D. 2008, ApJ, 678, 780
- Guedes, J., Madau, P., Kuhlen, M., Diemand, J., & Zemp, M. 2009, ApJ, 702, 890
- Haiman, Z., Kocsis, B., & Menou, K. 2009, ApJ, 700, 1952
- Harfst, S., Gualandris, A., Merritt, D., Spurzem, R., Portegies Zwart, S., & Berczik, P. 2007, New Astronomy, 12, 357
- Herrmann, F., Hinder, I., Shoemaker, D., & Laguna, P. 2007, Classical and Quantum Gravity, 24, 33
- Hills, J. G. 1975, Nature, 254, 295
- Just, A., Khan, F. M., Berczik, P., Ernst, A., & Spurzem, R. 2011, MNRAS, 411, 653
- Komossa, S. 2002, Reviews in Modern Astronomy, 15, 27
- Komossa, S., & Bade, N. 1999, A&A, 343, 775
- Komossa, S., & Merritt, D. 2008, ApJ, 683, L21
- Komossa, S., Zhou, H., & Lu, H. 2008, ApJ, 678, L81
- Lightman, A. P., & Shapiro, S. L. 1977, ApJ, 211, 244
- Lippai, Z., Frei, Z., & Haiman, Z. 2008, ApJ, 676, L5
- Liu, F. K. 2004, MNRAS, 347, 1357
- Liu, F. K., & Chen, X. 2007, ApJ, 671, 1272
- Liu, F. K., & Wu, X.-B. 2002, A&A, 388, L48
- Liu, F. K., Li, S., & Chen, X. 2009, ApJ, 706, L133
- Liu, F. K., Wu, X.-B., & Cao, S. L. 2003, MNRAS, 340, 411
- Liu, F. K., Xie, G. Z., & Bai, J. M. 1995, A&A, 295, 1
- Liu, F. K., Zhao, G., & Wu, X.-B. 2006, ApJ, 650, 749
- Lodato, G., King, A. R., & Pringle, J. E. 2009, MNRAS, 392, 332
- Loeb, A. 2007, Physical Review Letters, 99, 041103
- Madau, P., & Quataert, E. 2004, ApJ, 606, L17
- Merritt, D. 2006, Reports on Progress in Physics, 69, 2513
- Merritt, D., & Poon, M. Y. 2004, ApJ, 606, 788
- Merritt, D., Berczik, P., & Laun, F. 2007, AJ, 133, 553
- Merritt, D., Milosavljević, M., Favata, M., Hughes, S. A., & Holz, D. E. 2004, ApJ, 607, L9
- Merritt, D., Schnittman, J. D., & Komossa, S. 2009, ApJ, 699, 1690
- Milosavljević, M., & Phinney, E. S. 2005, ApJ, 622, L93
- O’Leary, R. M., & Loeb, A. 2009, MNRAS, 395, 781
- O’Leary, R. M., & Loeb, A. 2011, arXiv:1102.3695
- Peres, A. 1962, Physical Review, 128, 2471
- Peters, P. C. 1964, Physical Review , 136, 1224
- Pretorius, F. 2005, Physical Review Letters, 95, 121101
- Rees, M. J. 1988, Nature, 333, 523
- Schnittman, J. D., & Krolik, J. H. 2008, ApJ, 684, 835
- Sesana, A., Haardt, F., & Madau, P. 2008, ApJ, 686, 432
- Shields, G. A., & Bonning, E. W. 2008, ApJ, 682, 758
- Sillanpaa, A., Haarala, S., Valtonen, M. J., Sundelius, B., & Byrd, G. G. 1988, ApJ, 325, 628
- Spitzer, L. 1987, Princeton, NJ, Princeton University Press, 1987, 191 p.
- Spurzem, R., et al. 2009, Computer Science - Research and Development, Volume 23, Numbers 3-4, 231
- Stone, N., & Loeb, A. 2011, MNRAS, 412, 75
- Valtonen, M. J., et al. 2008, Nature, 452, 851
- Vicari, A., Capuzzo-Dolcetta, R., & Merritt, D. 2007, ApJ, 662, 797
- Volonteri, M., Haardt, F., & Madau, P. 2003, ApJ, 582, 559
- Yu, Q. 2002, MNRAS, 331, 935