DETECTION OF IMBHs WITH GROUND-BASED GRAVITATIONAL WAVE OBSERVATORIES: A BIOGRAPHY OF A BINARY OF BLACK HOLES, FROM BIRTH TO DEATH

By following the stellar dynamics at the center of our Galaxy, we now have the most well-established evidence for the existence of a supermassive black hole (SMBH). The close examination of the Keplerian orbits of the so-called S stars (also called S0 stars, where the letter "S" stands simply for source) has revealed the nature of the central dark object located at the Galactic center. By following S2 (S02), the mass of Sgr A* was estimated to be about 3.7 × 10⁶ M_☉ within a volume with radius no larger than 6.25 light hours (Schödel et al. 2003; Ghez et al. 2003). More recent data based on 16 years of observations set the mass of the central SMBH to ∼4 × 10⁶ M_☉ (Eisenhauer et al. 2005; Ghez et al. 2005, 2008; Gillessen et al. 2009).

Massive black holes in a lower range of masses may exist in smaller stellar systems such as globular clusters. These are called intermediate-mass black holes (IMBHs) because their masses range between M ∼ 10² and M ∼ 10⁴ M_☉, assuming that they follow the observed correlations between SMBHs and their host stellar environments. Nevertheless, the existence of IMBHs has never been confirmed, though we have some evidence that could favor their existence (see Miller & Colbert 2004; Miller 2009, and references therein).

If we wanted to apply the same technique to detect IMBHs in globular clusters as we use for SMBHs in galactic centers, we would need ultra-precise astronomy, since the sphere of influence of an IMBH is ∼few arcseconds. For instance, for a 10⁴ M_☉ IMBH, the influence radius is of ∼5'' assuming a central velocity dispersion of σ = 20 km s⁻¹ and a distance of ∼5 kpc. The number of stars enclosed in that volume is only a few. Currently, with adaptive optics, one can aspire—optimistically—to have a couple of measurements of velocities if the target is about ∼5 kpc away on a timescale of about 10 yr. The measures depend on a number of factors, such as the required availability of a bright reference star in order to have a good astrometric reference system. Also, the sensitivity limits correspond to a K-band magnitude of ∼15, (B MS stars at 8 kpc, like, e.g., S2 in our Galactic center).

This means that, in order to detect an IMBH or, at least, a massive dark object in a globular cluster center by following the stellar dynamics around it, one has to have recourse to the Very Large Telescope Interferometer (VLTI) and to one of the next-generation instruments, the VLTI Spectro Imager (VSI) or GRAVITY (Gillessen et al. 2006; Eisenhauer et al. 2008). In this case, we can hope to improve the astrometric accuracy by a factor of ∼10. Only in that scenario would we be in the position of following closely the kinematics around a potential IMBH so as to determine its mass.

Gravitational-wave (GW) astronomy could contribute to IMBH detection. In recent years, the field has reached a milestone with the construction of an international network of GW interferometers that have achieved or are close to achieving their design sensitivity. Moreover, the first-generation ground-based detectors Laser Interferometer Gravitational-Wave Observatory (LIGO) and Virgo will undergo major technical upgrades in the next five years that will increase the volume of the observable universe by a factor of 1000.³

The data that will be taken by the advanced interferometers are expected to transform the field from GW detection to GW astrophysics. The availability of accurate waveform models for the full binary of IMBH (hereafter BBH) coalescence in order to construct templates for match-filtering is crucial in the GW searches for compact binaries. The construction of these kind of templates has recently been made possible, thanks to the combination of post-Newtonian (PN) calculations of the BBH inspiral and numerical relativity (NR) simulations of the merger and ringdown. Two approaches to this problem are the effective-one-body techniques (Buonanno & Damour 1999; Buonanno et al. 2009) and the phenomenological matching of PN and NR results (Ajith et al. 2008, 2009; Santamaria et al. 2010). In this paper, we use the results of the latter.

The structure of this paper is as follows. In Section 2, we expand the astrophysical context of this problem and give a description of the different efforts made to address the evolution of a BBH in a stellar cluster from its birth to its final coalescence. In Section 3, we introduce the techniques used in the data analysis of the waveform modeling of BBH coalescences, present our hybrid waveform model, and use the model to compute and discuss expected signal-to-noise ratios (S/Ns) in present and future GW detectors. The use of the new waveform model allows us to give an improved estimate for the number of events one can expect for several physical configurations in Advanced LIGO and the Einstein Telescope (ET), which we present in Section 4. We conclude with a summary of our results and future prospects of our work in Section 5.

The aim of this section is not to give a detailed explanation of the processes of formation of IMBHs and BBHs, but a general introduction of the two different scenarios that play a role in the formation of BBHs.

2.1. Birth

2.2. Growing Up (Shrinking Down): The Role of Triaxiality on Centrophilic Orbits

In the case of the double-cluster channel, the cluster, which is in rotation, results from the merger of the two initial clusters and may develop a triaxiality sufficient to produce sufficient centrophilic orbits. These "boxy" orbits, as seen by Berczik et al. (2006), are typical of systems that do not possess symmetry around any of their axes. In contrast to loop orbits, a characteristic of spherically symmetric or axisymmetric systems, "boxy" orbits bring stars arbitrarily close to the center of the system, since it oscillates independently along the three different axes. Therefore, such stars, due to the fact of being potential sling shots, can feed the process of shrinkage of the BBH semi-major axes by removing energy and angular momentum from it after a strong interaction. In the strong triaxial systems of Berczik et al. (2006), the rotation caused in the process of merger creates an unstable structure in the form of a bar. Within the bar the angular momentum will not be conserved and thus the BBH loss cone is full due to stars on centrophilic orbits, independent of the number of stars N_⋆. On the other hand, for BBH systems in clusters, the role of Brownian motion has an impact in the loss cone (see Amaro-Seoane et al. 2010). In the models of Amaro-Seoane & Freitag (2006), the initial conditions are a realistic parabolic merger of two stellar clusters. The resulting merged cluster does not show the strong axisymmetry of Berczik et al. (2006). In the simulations we address for the results of this work, the BBH (of IMBHs) is not stalling, in spite of the reduced number of centrophilic orbits due to the architecture of the stellar system.

In Figure 1, we show the role of the cluster symmetry explicitly by depicting the evolution of the triaxiality of the cluster formed as a result of the merger of the two clusters for our fiducial model in the case of the double-cluster channel (which is the reference model of Amaro-Seoane & Freitag 2006). After a merger which is the result of a parabolic orbit, the final system is oblate rather than prolate; i.e., a ∼ b>c, where a, b, and c are the cluster axes. At the outskirts, the resulting merged cluster is flatter and at the center the BBH makes it rather spherical. Amaro-Seoane et al. (2009) addressed the single-cluster channel scenario after the formation of the IMBHs and used additional simulations to further evolve the BBH. They used scattering experiments of three bodies including relativistic precession to the first PN order, as well as radiation reaction caused by GW, so that they did not have to integrate every single star in the cluster to understand the posterior evolution of the BBH. In their work, between the strong encounters, a and e of the BBH were evolved by resorting to the quadrupolar formulas of Peters (1964).

Zoom In Zoom Out Reset image size

The BBH will be completely circularized when it reaches the frequencies probed by Advanced LIGO and the ET, because the emission of GWs takes over the dynamics of the system.

2.3. Death

While the emission of GWs is present all the time from the very first moment in which the BBH is formed, the amplitude and frequency of the waves is initially so low that no present or planned detector would be able to register any information from the system. Only when the semi-major axis shrinks sufficiently does the frequency increase enough to "enter" the LISA band, which we assume starts at 10⁻⁴ Hz. The BBH then crosses the entire detector window during its inspiral phase, as we can see in Figure 2. We depict the evolution of a BBH of mass 439.2 + 439.2 M_☉. The reason for this particular choice of masses is to give the reader a point of reference to understand the whole picture. Recently, Amaro-Seoane et al. (2010) included the effect of the rotation of the host cluster and addressed the dynamical evolution of the global system with a BBH of that mass. The authors have shown that LISA will see the system of Figure 2 with a median S/N of a few tens. The fact that the system merges outside its band prevents LISA from observing the loudest part of the BBH coalescence. In order to follow the system at this early stage of its evolution in the LISA band, a simple PN approach suffices for modeling the GW radiation. We are far enough from the highly relativistic regime and only the inspiral phase of the BBH coalescence is visible to the space antenna, with a rather low S/N compared to the posterior phases in the evolution of the BBH.

Zoom In Zoom Out Reset image size

As the binary system depicted in Figure 2 leaves the LISA band and enters the strong field regime, higher order PN corrections and eventually input from NR simulations need to be considered in order to model the GW waveform. Three reference frequencies in the evolution of a compact BBH that approaches its merger are the innermost stable circular orbit (f_ISCO) of a test particle orbiting a Schwarzschild black hole, the light-ring frequency (f_LR) corresponding to the smallest unstable orbit of a photon orbiting a Kerr black hole and the fundamental ringdown frequency (f_FRD) of the decay of the quasi-normal modes computed by Berti et al. (2005).

For the binary system shown in Figure 2, the values of these three frequencies are $f_{\rm ISCO}|_{878.4 M_{\odot }} \simeq 5 {\rm Hz}$, $f_{\rm LR}|_{878.4 M_{\odot }} \simeq 14.2\, {\rm Hz}$, and $f_{\rm FRD}|_{878.4 M_{\odot }} \simeq 21.4 \, {\rm Hz}$. Should such a binary exist at a distance of 100 Mpc (the choice for this number is based solely on the fact that is easily scalable), and if it were to be detected with Advanced LIGO, it should produce a sky-averaged S/N of ∼450, assuming a low frequency cutoff of 10 Hz. To that total S/N, the contribution of parts of the inspiral happening before the system reaches the characteristic frequencies f_ISCO, f_LR, and f_FRD would be 0%, 37%, and 95%, respectively. Figure 3 illustrates the same percentages for binaries with total masses between 100 and 2000 M_☉. It is immediately noticed that for the binaries of IMBHs of interest to this study, most of the S/N that these binaries will produce in Advanced LIGO comes from the last stages of the BBH coalescence. In Figure 4 we depict the effect that the IMBH binary systems of different total masses and orientations would have in the time-domain output of the Advanced LIGO detector if such sources existed at a distance of 1 Gpc.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We can estimate the time that the binary system takes to evolve from $f=0.01\,\rm {Hz}$, a frequency where the BBH can be seen by LISA, to the lower cutoff frequency of 10 Hz of Advanced LIGO or of 1 Hz of the ET. A lower order approximation based on the Newtonian quadrupole formula (Peters 1964) leads to the following expression for the evolution of the frequency in terms of the chirp mass ${\cal {M}}= (m_1 m_2)^{3/5}M_{\rm tot}^{-1/5}$ and frequency of the system

$$\begin{equation} \frac{df}{dt}=\frac{96}{5}\pi ^{8/3}{\cal {M}}^{5/3}f^{11/3}. \end{equation} \tag{ 1 }$$

We find a delay of only 0.2 yr (80 days) for a BBH with total mass M = 878.4 M_☉ to go from 0.01 Hz to the beginning of the Advanced LIGO band and almost similar numbers to the beginning of the ET band (the evolution of the system is extraordinarily quick in the late inspiral phase, which explains the fast evolution from 1 to 10 Hz). In view of these figures, LISA could be used as an "alarm" to prepare ground-based detectors to register the final coalescence in detail, the death of the BBH as such, by adjusting their sweet spots (the most sensitive part of the detector) to the particular BBH. The high accuracy of parameter estimation during the inspiral phase of which LISA is capable could be combined with the information obtained from the large-S/N triggers that the BBH merger and ringdown will produce in Advanced LIGO or ET and thereby achieve a more complete characterization of the system.

Accurate theoretical modeling of the gravitational radiation h(t) emitted by a BBH is key to improving BBH detectability and parameter estimation. While the PN theory is valid to model the early inspiral phase of the BBH evolution in the LISA band, an exact description of the merger and ringdown stages is only possible via NR calculations. Several approaches have been proposed to match PN and NR in order to obtain a full waveform for the entire coalescence, most notably for non-spinning configurations in Buonanno & Damour (1999, 2000), Buonanno et al. (2009), Ajith et al. (2007, 2008), and recently introducing spins for non-precessing configurations in Pan et al. (2010) and Ajith et al. (2009).

For the purpose of the S/N and horizon distance calculations shown in this paper, we have chosen a new procedure for the construction of hybrid PN–NR waveforms in the frequency domain proposed by Santamaria et al. (2010). The construction procedure developed for matching PN and NR data is sketched in Figure 5; further details on the fitting procedure can be found in Santamaria et al. (2010). In essence, the model consists of a phenomenological fit to hybrid PN+NR waveforms for spinning, non-precessing BBH systems with comparable masses (0.15 ≲ η < 0.25). This is compatible with the simulations of Gürkan et al. (2006), which find pairs of VMSs with mass ratios close to 1. The phenomenological waveforms are parameterized by three physical parameters: the total mass of the binary M, the symmetric mass ratio η, and the spin parameter χ ≡ (1 + δ)χ₁/2 + (1 − δ)χ₂/2, where δ ≡ (m₁ − m₂)/M, χ_i = S_i/m²_i and S_i is the spin angular momentum of the ith black hole. Only the dominant mode ℓ = 2,  m = 2 of the gravitational radiation enters in the model; the effects of the higher modes, which become increasingly significant as η decreases, are neglected. For comparable-mass scenarios, this restriction does not substantially affect our results. The waveform model used in our calculations is valid for binary systems with aligned spins only, which represents a first step toward the incorporation of the spins of the BHs in full inspiral-merger-ringdown (IMR) models. A discussion including precessing systems with arbitrary spins to this date still requires work regarding availability of NR simulations and development appropriate techniques to match precessing merger data with PN. Work in this regard is ongoing.

Zoom In Zoom Out Reset image size

For the results presented in this article we choose to focus on three configurations: (1) equal-mass, non-spinning, (2) equal-mass, equal-value spins aligned with the direction of the total angular momentum and magnitude χ = 0.75, and (3) mass ratio 1:3 (η = 0.19), non-spinning. Provided with the full IMR waveforms given by our model, we are interested in assessing the detectability of these systems by current and future GW observatories. The S/N of a model waveform with respect to the output stream of the detector is the quantity typically quoted to signify the detectability of a signal. The S/N ρ produced by a GW signal h(t) in a detector can be computed as (see, e.g., Thorne 1987; Finn 1992)

$$\begin{equation} \rho ^{2} \equiv 4\int _0^\infty \frac{\, \tilde{h}(f)\, \tilde{h}^\ast (f)}{S_n(f)}\, d f = \int _0^\infty \frac{\vert 2\, \tilde{h}(f)\, \sqrt{f} \vert ^{2}}{S_n(f)}\, d\ln {f}, \end{equation} \tag{ 2 }$$

where $\tilde{h}(f)$ is the Fourier transform of the strain h(t). In the last equation, S_n(f) represents the one-sided noise spectral density of the detector.

Figure 6 provides a graphical representation of the detectability of several BBHs by current and future generations of ground-based GW detectors. Displayed are the design sensitivity of current initial LIGO and Virgo (sensitivities that have been met or approximately met during the S5/VSR1 data taking), the proposed noise curves of Advanced Virgo and two possible configurations of Advanced LIGO (broadband or "base" and optimized for 30–30 M_☉ BBHs), and the designed noise budget for the ET in its broadband, "base," and "xylophone" (Hild et al. 2010) configurations. The hybrid waveforms constructed according to the model of Figure 5 for the three chosen configurations has been conveniently scaled to represent BBH systems with total mass 200, 500, and 1000 M_☉.

Zoom In Zoom Out Reset image size

As the right-hand side of Equation (2) suggests, plotting the quantity $2\, \vert \tilde{h}(f)\vert \, \sqrt{f}$ versus $\sqrt{S_h(f)}$ allows for direct visual comparison of the importance of each of the stages of the BBH coalescence. The three frequencies f_ISCO, f_LR, and the Lorentzian ringdown frequency f_LRD = 1.2f_FRD are marked on top of the configuration (1) waveforms with solid circles, squares, and open squares, respectively. One can immediately appreciate that systems with total mass above 500 M_☉ fall almost completely below the 40 Hz "seismic wall" of the initial LIGO detectors; however they will become very interesting sources for the second generation of GW interferometers and the proposed ET. Indeed, as we show in Section 2, they will also be seen by the future space-borne LISA. Additionally, it is easy to appreciate why the hang-up configuration (3) will produce larger S/Ns and therefore will be seen to further distances, for it merges at higher frequencies with respect to a non-spinning configuration with the same η. In contrast, for equal spin values, the S/N decreases with smaller symmetric mass ratios (compare configurations (1) and (3)).

In Figure 7, we compute the S/N expected for these sources in each of the above-mentioned detectors as a function of the redshifted total mass of the system M_z = (1 + z) M_BBH, for optimally oriented and located sources at a distance of 6.68 Gpc (z = 1). We have cross-checked our S/N results with those computed by Boyle et al. (2009) for initial LIGO at 100 Mpc, obtaining S/N values within 1% of those quoted by them. A direct comparison with Gair et al. (2009) shows a disagreement of ∼30% in the computed S/Ns for Advanced LIGO and the ET, which might be attributed to the different waveforms and cosmological model used.

Zoom In Zoom Out Reset image size

Unsurprisingly, the S/Ns calculated for the third generation of ground-based detectors surpass the expectations for initial and Advanced Virgo and LIGO at all masses. At z = 1, S/Ns of the order of 0.1 are expected for current LIGO and Virgo interferometers for binaries with total mass up to a few hundred M_☉. These values would scale to S/Ns of ∼1 at the closer distances typically surveyed by the initial interferometers, which are of the order of tens to few hundreds of Mpc. Signals with single-interferometer S/Ns below the commonly used threshold of S/N = 5.5 are most likely to be missed; therefore, IMBH binary coalescences are of limited interest for the first-generation detectors—and even for their current, enhanced counterparts—which are in turn most sensitive to neutron star binaries and stellar-mass back hole binaries.

Advanced LIGO and Virgo will be able to measure averaged S/Ns of the order of 1–20 at z = 1, with a maximal response to BBH systems with total mass in the range of 300–800 M_☉. S/N values of tens are well above the threshold established for detection. Therefore, should binaries of IMBHs exist in our neighboring universe, the second-generation of GW interferometers should be able to detect them, for they will be loud enough to stand above the detector noise. As we will see later, the non-negligible IMBH binary coalescence rate for Advanced LIGO indicates that these sources should be taken into account in future match-filtered searchers for inspiral binaries. Due to the moderate S/N of the signals, a potential detection might be sufficient to confirm the existence of binaries of IMBHs but not enough to determine all its parameters, such as mass ratio and spins, with sufficient accuracy; For the advanced detectors, Cutler & Flanagan (1994) estimate an accuracy of ∼1% in the reduced mass μ = m₁m₂/M using PN templates with negligible spins at S/N ∼ 10. When they take the spins into account, Δμ increases by ∼50. It is expected that the use of full IMR templates will improve these figures, but the final accuracy will always be limited by the slightly above-threshold S/Ns expected in Advanced LIGO and Virgo.

The ET will instead measure S/N values within the 10² range at z ∼ 1, and it is expected to be sensitive to binaries with total masses of the order of 10⁴ M_☉, a significantly larger range than that surveyed by Advanced LIGO and Virgo. It is noticeable how the ET xylophone configuration increases the detectability of binaries with masses above 1000 M_☉ with respect to the broadband ET configuration. This is due to its improved sensitivity precisely at frequencies in the range of 1–30 Hz, which is where systems of mass above thousands of solar masses accumulate most of their S/N (see Figure 6). The large-S/N events that binaries of IMBHs would produce in the ET are of great importance for astrophysics; the reason being that as the S/N increases, the accuracy in parameter estimation also does. In the limit of large S/N, the inverse of the Fisher information matrix Γ_a b = (h_, a|h_, b) is an estimator of the errors in the recovered parameters. At the relative large S/Ns potentially produced by IMBH binaries in the ET, the possibility of extracting their mass ratios and spins to high accuracies would be revolutionary for characterizing the IMBHs populations.

As for LISA, BBHs with masses of hundreds of solar masses will be seen with a moderate S/N (see Amaro-Seoane et al. 2010, for a detailed study of the parameter extraction)—it is only at masses above tens of thousands of solar masses that LISA will start taking over the ground-based observatories. Although the space antenna will be most sensitive to BBH binaries with masses in the range of 10⁶⁻⁷ M_☉, the possibility that it can act as a complementary observatory for the ET for IMBH binaries is very promising. Parameter accuracy studies for IMBHs in LISA are already available using the inspiral part of the coalescence (including also non-negligible eccentricities; see Amaro-Seoane et al. 2010), and indicate that masses and sky positions will be recovered with a high accuracy level. In order to complete the characterization of IMBHs with the information given by the second and third generations of ground-based detectors, a comprehensive study of parameter recovery taking the BBH coalescence into account is very much desirable.

Miller (2002) estimated for the first time the event rate for intermediate-mass mergers of IMBHs (i.e., typically stellar black holes merging with IMBHs) in clusters by calculating the luminosity distance for the inspiral, merger, and ringdown (Flanagan & Hughes 1998) out to which these three stages can be detected with an S/N larger than 10. In his approach, the maximum distance for the detector was 3 Gpc (z ∼ 0.53), with no cosmological corrections. The event rate was calculated as R = ∫(4/3) D(M)³ ν(M) n_ng f(M) dM. In this equation, n_ng is the number density of globular clusters, which was taken to be n_ng ∼ 8h³/Mpc³, as in the work of Portegies Zwart & McMillan (2000). The rate of coalescence of stellar-mass compact objects with the IMBH is ν(M) and f(M) = dN/dM is the mass distribution of massive enough black holes in clusters. Obviously, ∫f(M) dM = f_tot < 1. Miller (2002) uses the estimation of Flanagan & Hughes (1998) for it and finds that a few per year should be detectable during the last phase of their inspiral. Two years later, Will (2004) revisited the problem using matched filtering for the parameter estimation, an updated curve for the sensitivity of the detector and restricted PN waveforms to calculate an analytical expression for the luminosity distance D_L; his estimation is a detection rate for binaries of about 1 Myr⁻¹ in a mass-range of (10:100) M_☉.

Subsequently, Fregeau et al. (2006) calculated the number of events that initial and Advanced LIGO (and LISA) could see from the single-cluster channel. In their estimation, they assume that the VMSs formed in the runaway scenario do not merge into one, but evolve separately and eventually each form an individual IMBH, following the numerical results of the Monte Carlo experiments of Gürkan et al. (2006). They derive a generalized form of the event rate which can be summarized as follows:

$$\begin{eqnarray} R &=& \frac{d N_{\rm event}}{dt_0} = \int _0^{z_{\rm max}}\frac{d^2 M_{\rm SF}}{dV_c dt_e} \, g_{\rm cl}\, g \,\frac{dt_e}{dt_0}\, \frac{dV_c}{dz} \nonumber\\ &&\times\int _{M_{\rm cl,\,min}}^{M_{\rm cl,\,max}}\frac{dN^2_{\rm cl}}{dM_{\rm SF,\,cl} dM_{\rm cl}}\,dM_{\rm cl} \,dz. \end{eqnarray} \tag{ 3 }$$

In this expression, d²M_SF/dV_cdt_e is the star formation rate (SFR) per unit of comoving volume per unit of local time, g_cl is the fraction of mass that goes into the massive clusters of interest, g is the fraction of massive clusters which form IMBHs, dt_e/dt₀ = (1 + z)⁻¹ is the relation between local and observed time, dV_c/dz is the change of comoving volume with redshift, dN²_cl/dM_SF cl dM_cl is the distribution function of clusters over individual cluster mass M_cl and total star-forming mass in clusters M_SF cl, and z_max is the maximum redshift that the detector is capable of observing.

In order to compute the integral above, an estimation of the observable volume of each individual detector is required. A commonly used measure of the reach of a detector is the horizon distance D_h, defined as the distance at which a detector measures an S/N = 8 for an optimally oriented and optimally located binary, i.e., an overhead, face-on orbit. Suboptimally located and oriented sources are detected with S/N = 8 at closer distances. However, in order to compute merger rates, an average distance over all possible orientations is more meaningful. The orientation-averaged distance represents the cube root of the total volume to which a detector is sensitive, assuming uniformly distributed sources, and is 2.26 times smaller than the horizon distance D_h (Finn & Chernoff 1993). Moreover, at the distances that Advanced LIGO and the ET are expected to survey, a certain cosmological model needs to be assumed. We adopt the standard ΛCDM universe with parameters given by the first five years of the WMAP sky survey (Hinshaw et al. 2009). These are Ω_Λ = 0.73, Ω_b = 0.046, Ω_c = 0.23, H₀ = 70.5 km s⁻¹ Mpc⁻¹, and t₀ = 13.72 Gyr. Using the full IMR waveforms described in Section 3 and the corresponding redshift function z(d) for the ΛCDM model, we compute the orientation-averaged distance D_h/2.26 for non-spinning systems with symmetric mass ratio η = 0.25, 0.19 and for an equal-mass system with spins aligned in the direction of the angular moment and total spin χ = 0.75 (hang-up configuration). Computing D_h(M_z) is essentially equivalent to inverting Equation (2) for ρ(d, M(d)) for ρ = 8 using the relation z(d) given by the cosmological model. The results for D_h(M_z) for our three configurations can be seen in Figure 8 for Advanced LIGO and ET.

Zoom In Zoom Out Reset image size

The limit for the integral in z is thus given by the observable volume of the detector of interest, quantified using the orientation-averaged distance. The maximum redshift z_max is in turn a function of the mass of the binary system, its configuration—mass ratio, spin—and the particular waveform model used in the calculation of the horizon distance. We present our values for z_max in Figure 9, which essentially shows the same data as in Figure 8 but expressed in terms of redshift versus total mass of the binary. As we see in Figure 8, the maximum values for the orientation-averaged distance for ET are as large as z ∼ 10 for configuration (2) at M_z ∼ 300 M_☉. This implies that the ET will be able to probe the different proposed scenarios to produce the first generation of black hole seeds, as pointed out by Sesana et al. (2009). However, at these large cosmological distances, the stellar formation rate is unknown and the validity of the rate integral cannot be stated. We therefore set a maximum cutoff value of z_max = 5 in the computation of Equation (5). The final value of z_max(M_BBH) that we have used in the computation of the rates is shown in Figure 9 for our three particular physical configurations.

Zoom In Zoom Out Reset image size

Regarding the term in the dM_cl integral in Equation (5), three different parameterizations of the SFR are available in the literature, i.e., R_SF1,2,3 as summarized by Equations (4), (5), and (6) of Porciani & Madau (2001). The three models are similar for distances up to z ∼ 2, where SFR peaks, differing from there on. For the results shown here, we have compared the three of them. As for the distribution of cluster masses, the factor can be approximated as

$$\begin{equation} \frac{dN^2_{\rm cl}}{dM_{\rm SF,\,cl} dM_{\rm cl}} = \frac{f(M_{\rm cl})}{\int M_{\rm cl}\,f(M_{\rm cl})\,dM_{\rm cl}}. \end{equation} \tag{ 4 }$$

In order to compute the integral in the denominator, we take dN_cl/dM_cl ∝ 1/M²_cl, following the power law form observed by Zhang & Fall (1999) for young star clusters in the Antennæ. The validity of assuming the same law for the larger volume of the universe surveyed by Advanced LIGO or the ET is, however, a generalization not based on direct observations. Thus, we should take this premise with care. By assuming an efficiency factor of f_GC ∼ 2 × 10⁻³, based on the results of Gürkan et al. (2006), we can set the values for the limits in M_cl = f_GC/M_BBH according to the masses of the IMBH binaries of interest. In our case, taking the standard definition of IMBH into account, this is M_cl, max/M_cl, min = M_BBH, max/M_BBH, min = 2 × 10⁴M_☉/2 × 10² M_☉. The integral can now be expressed as

$$\begin{eqnarray} R &=& \frac{f_{\rm GC}}{\ln (M_{\rm cl,\,max}/M_{\rm cl,\,min})} \, g \, g_{\rm cl} \int _{M_{\rm BBH,\,min}}^{M_{\rm BBH,\,max}}\frac{dM_{\rm BBH}}{M^2_{\rm BBH}}\nonumber\\ &&\times\int _0^{z_{\rm max}(M_{\rm BBH})} \mbox{SFR}_i(z)\, F(z)\,\frac{1}{(1+z)}\, \frac{dV_c}{dz} dz, \end{eqnarray} \tag{ 5 }$$

where M_{BBH, max(min)} is the range of total mass of the BBH that we are considering, SFR_i(z),  i = 1, 2, 3 is any of the three considered stellar formation rates of Porciani & Madau (2001) and

$$\begin{equation} F(z) = \frac{\sqrt{\Omega _M(1+z)^3 + \Omega _k(1+z)^2 + \Omega _\Lambda }}{(1+z)^{3/2}} \end{equation} \tag{ 6 }$$

is the factor that relates the stellar formation rate function in different cosmologies with respect to the Einstein-de Sitter universe.

Regarding the values of the terms g and g_cl, it is unfortunately very little that we know about the initial cluster conditions required to form an IMBH binary. These factors therefore have large uncertainties. Following the existing literature, we leave the fraction of massive clusters that form IMBHs, g, as a parameter and set it to 0.1 as an example. Nevertheless, as proven in the simulations of Freitag et al. (2006b), it could be as large as 0.5. As for the fraction of mass going into massive clusters, g_cl, previous works have also tentatively set it to 0.1. Nevertheless, the results of McLaughlin (1999) seem to indicate that this value might be too optimistic. McLaughlin (1999) estimates an empirical cluster formation efficiency by mass of ≡ M^init_GCs/M^init_GAS, where M^init_GCs is the initial globular cluster population and M^init_GAS is the initial reservoir available in the protogalaxy. He then infers a universal value of ≃ 0.25% after evaluating a sample of 97 giant ellipticals, brightest cluster galaxies, and faint dwarfs. He finds identical results for the Population II spheroid of the Milky Way. Therefore, we have chosen to set our g_cl to 1/400, noting that, in any case, the final rates scale trivially with both g and g_cl.

We have evaluated the integral in Equation (5) for the three stellar formation rates we consider and the three binary configurations (1), (2), and (3) for Advanced LIGO and the ET. Our results are summarized in Table 1. We find that the event rate does not depend strongly on the assumed stellar formation rate; the differences are of the order of ∼5% for Advanced LIGO and ∼30% for the ET. In the case of Advanced LIGO, the insensitivity to the stellar formation model is easily explained by the fact that the three models of Porciani & Madau (2001) are very similar until z ∼ 2; our integration in z stops at z ∼ 1 (see the red curves in Figure 9), therefore ignoring contributions at higher z. For the ET, the differences among SFR_i(z),  i = 1, 2, 3 between z ∼ 2 and 5 are very much attenuated by the rapid decrease in comoving volume at high redshift. In order to assess the effect of stopping the integration at z = 5, we have also computed the rates for the ET without imposing this condition. The differences amount to ∼10%, therefore we conclude that the rates estimations are not very sensitive to very high-z SFRs. Nonetheless, we should stress again that the lack of stellar formation data at these large cosmological distances makes assumptions at z>5 largely speculative and thus we restrict ourselves to presenting the rates calculated with z_max clamped at a value of 5. Given the inherent uncertainties present in our approach, which can amount to a few orders of magnitude, the factors associated with the choice of SFR do not represent a major source of error. We therefore quote the results found for SFR₂(z) only, corresponding to the SFR that increases up to z ∼ 2 and keeps constant afterward.

Table 1. Event Rates Formed in the Single-cluster Channel for IMBH Binaries Potentially Observable by Advanced LIGO and the ET per year

Detector	Configuration	$R\left[\left(\displaystyle {\frac{g}{0.1}}\right)\, \left(\displaystyle {\frac{g_{\rm cl}}{1/400}}\right)\,{\rm yr}^{-1}\right]$
Advanced LIGO	η = 0.25,  χ = 0	0.85
	η = 0.25,  χ = 0.75	2.36
	η = 0.19,  χ = 0	0.31
ET	η = 0.25,  χ = 0	21.5
	η = 0.25,  χ = 0.75	24.0
	η = 0.19,  χ = 0	16.8

Note. We take g = 0.1, g_cl = 1/400 as standard scaling values.

Download table as:  ASCIITypeset image

The event rates do, however, strongly depend on the spins and mass ratio of the binary. As expected, "loud" configurations like the hang-up case increase the event rate by a factor of ∼3 in the case of Advanced LIGO. The effect of the spin in the rates is more notable for larger values of the total spin of the binary; the correct determination of the spin distribution of IMBH binary systems, which is at present far from being completely understood, would be extremely useful to further quantify the impact of the spin on the total event rates. Smaller mass ratios decrease the event rates, the difference between η = 0.25 and η = 0.19 also being approximately a factor of 3. The differences are not so extreme in the case of the ET, due to the fact that we are cutting off z_max at a value of 5 and, thus, neglecting contributions at higher redshift which might amount to another ∼10%, as discussed above. Even so, it is evident that the expected rates increase with the total spin χ and the symmetric mass ratio η. The exact computation of the total rates for all possible physical configurations of the binary would necessitate further integration on the mass ratios and spins of the system. At present, it is not clear what that distribution might be; therefore, we simply summarize our results for all configurations under consideration in Table 1. We find event rates of the order of ∼1(g/0.1)(g_cl/1/400) yr⁻¹ for Advanced LIGO and of the order of ∼20(g/0.1)(g_cl/1/400) yr⁻¹ for the ET. These rates assume formation of IMBH binaries in the single-cluster channel.

For Advanced LIGO and a non-spinning configuration, Fregeau et al. (2006) found a expected detection rate of 10(g/0.1)(g_cl/0.1) yr⁻¹ (0.25 yr⁻¹ in our units), adopting a reach of 2 Gpc (z_max ∼ 0.37). This estimation, calculated assuming a ringdown-only search, underestimates the volume of the universe that Advanced LIGO would observe using full IMR templates, as we show in Figure 8. Our results for non-spinning systems are compatible with this observation; we quote rates larger than those of Fregeau et al. (2006) by a factor of 3 but otherwise of a compatible order of magnitude. In the case of the ET, Gair et al. (2009) quote a rate of ∼500(g/0.1)(g_cl/0.1) yr⁻¹ (∼13 yr⁻¹ in our units) for non-spinning configurations, a value of the same order of magnitude as the one found by us. The fact that Gair et al. (2009) use a different family of waveforms to model the coalescence and a fitted formula for the averaged distance based on EOBNR waveforms, together with the slightly different values in the integration limits explain the discrepancies in the exact figures. Therefore, our new results for the rates of IMBH binary coalescence for non-spinning systems agree reasonably well with previous works that used similar methods but different waveform and detector models. In addition, we have now quantified the effect that the spins of the BHs will have on the expected rates and have found it to be of a factor ∼3 for total spins of the binary as high as χ = 0.75

Thus far we have concentrated on the single-cluster channel scenario. Amaro-Seoane & Freitag (2006) give a prescription to calculate an estimate of the event rates for the double-cluster channel. This was based on the fact that the only difference, in terms of the event calculation, between both astrophysical scenarios involves, first, that in the double-cluster channel there is one single IMBH in one cluster and, second, that these two clusters must collide so that the IMBHs form a BBH when they sink to the center due to dynamical friction. As explained in Section 4 of Amaro-Seoane & Freitag (2006), the connection between the event rate estimation of the two channels is

$$\begin{equation} \Gamma ^{\rm \,doub} = P_{\rm merg}\,g\,\Gamma ^{\rm \,sing}, \end{equation} \tag{ 7 }$$

where Γ^doub is the event rate of the double-cluster channel, Γ^sing of the single channel and P_merg is the probability for two clusters to collide in the scenario of Amaro-Seoane & Freitag (2006). They find that the unknown parameter P_merg could have values in the range P_merg ∈ [0.1, 1]. The total event rate assuming that both formation channels are possible would be

$$\begin{equation} \Gamma ^{\rm \,tot} = \Gamma ^{\rm \,sing} \left(\frac{g}{0.1}\right) \left(\frac{g_{\rm cl}}{1/400}\right) \left(1+ P_{\rm merg}\,\left(\frac{g}{0.1}\right)\right). \end{equation} \tag{ 8 }$$

Assigning parameter P_merg its pessimistic and optimistic limits, we can compute the lower limit and upper limit of the event rates for Advanced LIGO and the ET, assuming contributions of the two channels. Taking into account the double-cluster channel increases the rates' upper limit by a factor of 2. For instance, for the equal-mass, non-spinning case, the values are

$$\begin{eqnarray} \Gamma _{\rm Adv.\,LIGO}^{\rm total} & \in [(0)\,0.94,\,1.7]\,\left(\frac{g}{0.1}\right) \left(\frac{g_{\rm cl}}{1/400}\right)\,{\rm yr}^{-1}\\ \Gamma _{\rm ET}^{\rm total} & \in [(0)\,23.7,\,43.0]\,\left(\frac{g}{0.1}\right) \left(\frac{g_{\rm cl}}{1/400}\right)\,{\rm yr}^{-1}. \end{eqnarray} \tag{ 9 }$$

Note that the results are quoted in terms of g and g_cl; the already-mentioned uncertainty in their values could increase these rates by at least 1 order of magnitude. These event rates are encouraging to address the problem of detection and characterization of systems of IMBH binaries with ground-based GW observatories. On the other hand, one should bear in mind that the existence of IMBHs has not yet been altogether corroborated, so that the pessimistic estimate is still somewhat "optimistic." This is why we have added a (0) in the previous rates as the absolute lower limit.

Even though we do not have any evidence of IMBHs so far, a number of theoretical works have addressed their formation in dense stellar clusters. If we were to follow the same techniques that have led us to discover the SMBH in our own Galaxy, we would need the VLT interferometer and next-generation instruments, such as the VSI or GRAVITY, which should be operative in the next ∼10 yr. An alternative or even complementary way of discovering IMBHs is via their emission of GWs when they are in a BBH system.

The identification and characterization of these systems rely on accurate waveform modeling of their GW emission, which has been made possible by the success of NR in simulating the last orbits of the BBH coalescence and the coupling of these results to analytical PN calculations of the inspiral phase. We use a PN–NR hybrid waveform model of the BBH coalescence based on a construction procedure in the frequency domain (see Santamaria et al. 2010, for details).

Using this hybrid waveform, we have estimated the S/N corresponding to the current and Advanced LIGO and Virgo detectors, the proposed ET and the space-based LISA at a distance of z = 1, i.e., 6.68 Gpc. The results indicate that IMBH binaries will produce S/Ns sufficient for detection in Advanced LIGO and Virgo and notably larger S/Ns in the ET, thus making them interesting sources to follow up. Eventual observations of IMBH binaries with future ground-based detectors could be complementary to those of LISA, which is expected to detect these systems with moderate S/Ns and to be more sensitive to SMBH binaries. More remarkably, in principle, if LISA and the ET are operative at the same time, they could complement each other and be used to track a particular event.

Furthermore, we have revisited the event rate of BBHs for various detectors and find encouraging results, within the inherent uncertainties of the approach. Our estimations are consistent with previous works, and additionally we have quantified the effect of the total spin of the binary in the expected event rates. We have estimated the distance to which Advanced LIGO and the ET will be able to see binaries of IMBHs. This quantity depends strongly on the mass ratio and spins of the binary. For Advanced LIGO, equal-mass, non-spinning configurations of observed total mass ∼200–700 M_☉ can be seen up to z ∼ 0.8. If the spins are aligned with the total angular momentum and significant (χ_1,2 ∼ 0.75), the reach increases to z ∼ 1.5 for observed total masses of ∼400 M_☉. The ET will be able to explore even more remote distances, reaching to z ∼ 5 and further. Our present knowledge of stellar formation at such large redshifts is incomplete; therefore, we have computed the event rates for the ET integrating only until z = 5. We have compared three star formation models and three different configurations of the binary. The effect of the particular formation model is not very significant; remarkably, the physical configuration does however influence the final rates. We provide the results for three physical configurations, taking into account both single- and double-cluster channels in the binary formation. For a fully correct calculation of the event rate integrated over all possible configurations, more detailed knowledge of the distribution of spins and mass ratios of IMBH binaries formed in globular clusters is required.

Advanced ground-based detectors are designed to be able to operate in different modes so that their sensitivity can be tuned to various kinds of astrophysical objects. Considering the importance of an eventual detection of a BBH, the design of an optimized Advanced LIGO configuration for systems with M ∼ 10²⁻⁴ M_☉ would be desirable to increase the possibility of observing such a system. In case an IMBH binary coalescence were detected, the recovery and study of the physical parameters of the system could serve to test general relativity and prove or reject alternative theories, such as scalar–tensor type or massive graviton theories. The direct identification of an IMBH with GWs will be a revolutionary event not only due to the uncertainty that surrounds their existence or their potential role in testing general relativity, but also because the information encoded in the detection will provide us with a detailed description of the environment of the BBH/IMBH itself.

We are indebted to C. Cutler, B. Sathyaprakash, and B. Krishnan for discussions about the work. We thank C. Reisswig and D. Pollney for providing us with their Cauchy-extrapolated data and F. Ohme for helping us to implement the PN formalism used for the waveform modeling. We are grateful to E. Robinson and D. Vanecek for comments on the manuscript. P.A.S.'s work was partially supported by the DLR (Deutsches Zentrum für Luft- und Raumfahrt). L.S. has been partially supported by DAAD grant A/06/12630, and the author thanks the IMPRS for Gravitational Wave Astronomy for support.