Binary super-massive black hole environments diminish the gravitational-wave signal in the pulsar timing band
We assess the effects of super-massive black hole (SMBH) environments on the gravitational-wave (GW) signal from binary SMBHs. To date, searches with pulsar timing arrays for GWs from binary SMBHs, in the frequency band nHz, include the assumptions that all binaries are circular and evolve only through GW emission. However, dynamical studies have shown that the only way that binary SMBH orbits can decay to separations where GW emission dominates the evolution is through interactions with their environments. We augment an existing galaxy and SMBH formation and evolution model with calculations of binary SMBH evolution in stellar environments, accounting for non-zero binary eccentricities. We find that coupling between binaries and their environments causes the expected GW spectral energy distribution to be reduced with respect to the standard assumption of circular, GW-driven binaries, for frequencies up to nHz. Larger eccentricities at binary formation further reduce the signal in this regime. We also find that GW bursts from individual eccentric binary SMBHs are unlikely to be detectable with current pulsar timing arrays. The uncertainties in these predictions are large, owing to observational uncertainty in SMBH-galaxy scaling relations and the galaxy stellar mass function, uncertainty in the nature of binary-environment coupling, and uncertainty in the numbers of the most massive binary SMBHs. We conclude, however, that low-frequency GWs from binary SMBHs may be more difficult to detect with pulsar timing arrays than currently thought.
keywords:black hole physics — galaxies: evolution — gravitational waves — methods: data analysis
The merger of a pair of galaxies hosting central super-massive black holes (SMBHs) is expected to result in the formation of a binary SMBH (Begelman, Blandford & Rees 1980). The central SMBHs sink in the merger remnant potential well through the action of dynamical friction, and form a bound binary when the mass within the orbit of the lighter SMBH is dominated by the heavier SMBH. As stars within the binary orbit are quickly ejected, the binary will decay further only if another mechanism to extract binding energy and angular momentum exists. Proposed mechanisms include slingshot scattering of stars on radial, low angular momentum orbits intersecting the binary (Frank & Rees, 1976; Quinlan, 1996; Yu, 2002), and friction against a spherical Bondi gas accretion flow (Escala et al., 2004) or a circum-nuclear gas disk (e.g., Roedig et al., 2011). If the orbital decay process can drive the binary to a small separation, gravitational-wave (GW) emission will eventually cause the binary to coalesce (e.g., Peters & Mathews, 1963; Baker et al., 2006).
In a cosmological context, merging dark matter halos follow parabolic trajectories (e.g., van den Bosch et al., 1999; Wetzel, 2011), implying large initial eccentricities (typically , Hashimoto et al., 2003) for the orbits of SMBHs sinking towards galaxy merger remnant centres. Steep stellar density gradients in merging galaxies may reduce this eccentricity; indeed, some models suggest that binary SMBHs are likely to be close to circular upon formation (Casertano, Phinney & Villumsen 1987; Hashimoto, Funato & Makino 2007; Polnarev & Rees, 1994). Slingshot interactions between binaries and individual stars again grow the eccentricities (e.g., Sesana, Haardt & Madau 2006; Quinlan, 1996; Berentzen et al., 2009; Khan et al., 2012), because binaries spend more time, and hence lose more energy, at larger separations. Roedig et al. (2011) found that binary SMBHs embedded in massive self-gravitating gas disks will have large eccentricities, between 0.6 and 0.8, at the onset of GW-dominated evolution. In the GW-dominated regime, however, binaries are expected to quickly circularise (e.g., Peters & Mathews, 1963; Baker et al., 2006).
There is no direct observational evidence for the existence of binary SMBHs (Dotti, Sesana & Decarli 2012). However, the GW emission from binaries prior to coalescence is an unambiguous signature of their existence. Observing GWs from binary SMBHs will enable binary SMBH physics, as well as models for the formation and evolution of the cosmological SMBH population, to be observationally tested. Here, we focus on the possibility of detecting GWs from binary SMBHs in the early parts of their GW-dominated evolutionary stages with radio pulsar timing arrays (PTAs; Foster & Backer, 1990; Manchester et al., 2013). PTAs are currently sensitive to GWs in the frequency band nHz, which is complementary to other GW detection experiments.
PTAs target both a stochastic, isotropic background of GWs from binary SMBHs (e.g., Yardley et al., 2011; van Haasteren et al., 2011) and GWs from individual binary systems (Yardley et al. 2010; Ellis, Siemens & Creighton 2012). The summed GW signal from all binary SMBHs in the Universe is expected to approximate an isotropic background, although individual binaries are potentially detectable at all frequencies within the PTA band (Ravi et al., 2012). Recent PTA results suggest that a large fraction of existing models for the GW background from binary SMBHs is inconsistent with observations (Shannon et al., 2013). However, most current predictions for the spectral shape (Phinney, 2001), statistical nature (Sesana, Vecchio & Volonteri 2009; Ravi et al., 2012) and strength (Sesana, 2013b) of the GW signal from binary SMBHs assume that all binaries are in circular orbits, and losing energy and angular momentum only to GWs. These assumptions correspond to the well-known power law GW background characteristic strain spectrum from binary SMBHs that is proportional to , where is the GW frequency.
Here, we present an examination of the properties of the GW signal from binary SMBHs given a realistic model for binary orbital evolution. We use a semi-analytic model for galaxy and SMBH formation and evolution (Guo et al., 2011, hereafter G11) implemented in the Millennium simulation (Springel et al., 2005) to specify the coalescence rate of binary SMBHs, and augment this with a framework (Sesana, 2010) for the evolution of binary SMBHs in stellar environments. We neglect gas-driven binary evolution. This is because massive galaxies at low redshifts, which are expected to dominate the GW signal from binary SMBHs, will typically be late-type and gas-poor (e.g., Yu et al. 2011; McWilliams, Ostriker & Pretorius 2012).
Two key phenomena in binary SMBH evolution affect the summed GW signal relative to the case of circular binaries evolving under
GW emission alone
Interactions between binary SMBHs and their environments will accelerate orbital decay compared to purely GW-driven binaries, reducing the time each binary spends radiating GWs. This may reduce the energy density in GWs at the lower end of the PTA frequency band (e.g., Wyithe & Loeb, 2003; Sesana, 2013a).
While circular binaries emit GWs at the second harmonics of their orbital frequencies, eccentric binaries emit GWs at multiple harmonics (Peters & Mathews, 1963). Given a population of binary SMBHs, this is expected to transfer GW energy density from lower frequencies in the PTA frequency band to higher frequencies (Enoki & Nagashima, 2007; Sesana, 2013a).
We consider the effects of both these phenomena on the GW signal from binary SMBHs relative to the circular, GW-driven case. We also examine the possibility of detecting bursts of GWs from individual eccentric, massive binaries. In §2, we outline the binary population model. We give our predictions for the summed GW signal in §3, along with a discussion of uncertainties in our model. We consider the possibility of detectable GW bursts in §4. Finally, we summarise our results in §5 and present our conclusions §6. Summaries of key PTA implications can be found at the ends of §3, §4 and §5. We adopt a concordance cosmology consistent with the Millennium simulation (Springel et al., 2005), with , , , and km s Mpc.
2 Description of modelling methods
2.1 The gravitational wave background from a cosmological source population
Consider a population of GW sources with comoving volume density at redshift , each radiating a GW luminosity per unit rest-frame frequency, , of . The specific intensity of GWs at the Earth from sources between redshifts and is
Here, , where is the observed GW frequency, and the comoving volume element is
where the Hubble parameter, , is given by , is the vacuum speed of light, and is the luminosity distance at redshift . Now, following Phinney (2001) and re-arranging, the energy density in GWs at the Earth per logarithmic frequency unit is, in any homogeneous and isotropic universe,
where , (with as the gravitational constant) is the critical mass density of the Universe, and is the proper time. The redshift is related to as
Although our Equation (4) is directly comparable to Equation (5) of Phinney (2001), we note that ours and Phinney’s expressions are only mathematically identical if we explicitly assume that each source radiates for an infinitesimal (proper) time.
2.2 The binary SMBH population at formation
To begin, we consider binary SMBHs with component masses , orbital semi-major axes and eccentricities , embedded in isotropic, unbound cuspy stellar distributions with velocity dispersions . Quinlan (1996) found that binary hardening caused by slingshot interactions with individual stars becomes effective at binary component separations of
We adopt a simple, one-parameter distribution for the eccentricities of binary SMBHs at formation,
based on the postulate that the semi-major and semi-minor axes of the orbits ( and respectively) are each log-normally distributed.
This is justified because (a) while many stellar encounters influence the values of and , the effects of these encounters
on the parameter values are heterogeneous, and (b) both and are strictly positive.
We hence model the ratio using a probability density function given by
Here, is the free parameter; larger values of correspond to typically larger binary eccentricities, and corresponds to a population of circular binaries. We do not consider any variation of with binary component masses or redshift, because there are no strong motivations for such variations. The eccentricity of a binary at is given by .
Let be a vector of parameters of binaries at formation. We denote the distribution of binaries in these parameters as . In this notation, the multivariate density function for a parameter vector with components indexed by an integer is given by
Binaries at formation have semi-major axes .
We use the results of the semi-analytic model of G11 to specify . As outlined in Ravi et al. (2012) and Shannon et al. (2013), the G11 results can be used to predict the coalescence rate of binary SMBHs. For this work, we only use coalescences with both and greater than , and only draw from the implementation of the G11 model in the Millennium simulation. We scale all SMBH masses by a factor of 1.9 (Shannon et al., 2013) to account for recent SMBH and galaxy bulge measurements.
We count the coalescing pairs of SMBHs in bins of , and (with widths , , and respectively) within the entire Millennium simulation box. Binaries are also randomly assigned values of using Equation (7). In this work, we consider four different initial binary eccentricity distributions defined by . Denoting the binary counts for different values of , , and by the discrete distribution , we have
where is the comoving volume of the Millennium simulation box (Springel et al., 2005). We average the distribution over merger order (Ravi et al., 2012) and 1000 realisations of the initial -distribution. We do not fit an analytic function to , as was done by Ravi et al. (2012) and Shannon et al. (2013). We discuss the possible consequences of this for our results in §3.3.3.
We relate the dark matter halo virial velocities, , of galaxies in the G11 model to spheroid stellar velocity dispersions (Baes et al., 2003; Marulli et al., 2008). This assumption is discussed further in §3.3.2. For each bin of , and , we find the average velocity dispersions of recently-merged galaxies in the G11 model hosting an SMBH of mass . We use these values to specify for each bin of the discrete distribution .
2.3 Evolution of binary SMBH orbits to the GW regime
We assume that all SMBH binary orbits decay through interactions with fixed, isotropic, unbound cuspy stellar backgrounds, and through GW emission. The former scenario has been extensively studied numerically by Quinlan (1996), Sesana et al. (2006) and Sesana (2010). We assume a power-law stellar density distribution within the binary gravitational influence radius for all galaxies prior to mergers. For the majority of this paper, we additionally assume a stellar density profile power-law index of corresponding to a mild stellar cusp (see Equation 1 of Sesana, 2010). We consider variations in these assumptions further in §3.3.2.
We evolve the binary eccentricities, , and semi-major axes, , through scattering by unbound stars and loss of energy and angular momentum to GWs using expressions for and from Equations (15) and (16) of Sesana (2010). The effects of the ejection of stars that are bound to the SMBHs (Sesana, Haardt & Madau 2008) are significant only for binary separations greater than , and we hence neglect this phenomenon.
We use the fits of Sesana et al. (2006) for the rates of evolution of binary semi-major axes and eccentricities based on numerical scattering experiments (the ’’ and ’’ coefficients respectively from Tables 1 and 3 of Sesana et al. 2006). We log-interpolate the published values at binary component mass ratios of interest. As Sesana et al. (2006) only provide rates of semi-major axis evolution for circular binaries, we assume here that the rate of semi-major axis evolution at a given semi-major axis is independent of eccentricity. This approximation leads to the semi-major axis evolution rate being underestimated by at most 20% for the most eccentric binaries (see Figure 3 of Sesana et al., 2006). We also only use the seven values for our initial binary eccentricities (i.e., ) considered by Sesana et al. (2006); see their Table 3. These are 0, 0.15, 0.3, 0.45, 0.6, 0.75 and 0.9.
By numerically integrating the expressions for and for each combination of and , we first calculate the binary eccentricities, , at a rest-frame orbital frequency of Hz. Binaries with this orbital frequency emit negligible GW power in the PTA frequency band. The orbital frequency of a binary is given by
Letting , we hence form the distribution function of binaries with orbital frequencies of Hz, , from the distribution of binaries at formation. If a binary at formation has Hz, we do not evolve the binary backwards in time to an orbital frequency of Hz.
To then specify the population of GW-emitting binary SMBHs, we need to calculate the numbers of binaries with different orbital frequencies. The GW luminosity, , per unit frequency, , of a binary SMBH depends on the masses and , the eccentricity , and the orbital frequency (Peters & Mathews, 1963). The functional form of is given in, for example, Equation (2.6) of Enoki & Nagashima (2007). We now define a new parameter vector .
The distribution can be used to specify the distribution function using a continuity equation similar to Equation (35) of Phinney (2001):
where is the number of coalescences of binary SMBHs with parameters per unit proper time . The derivative is equivalent to , where is given in Equation (15) of Sesana (2010). The solution is
We also associate each value of with a unique value of by further integrating the expression for from Sesana (2010).
Then, from Equation (3), we have
Recall that . For consistency with other works, we calculate the characteristic strain spectrum, defined as
We perform the integral in Equation (13) over between Hz. The upper orbital frequency limit corresponds to GW emission that is outside the PTA frequency band, even for binaries at high redshifts. For eccentric binaries, we consider radiation up to the 100th harmonic of (Peters & Mathews, 1963). We assume that binaries reach their last stable orbits at separations of three Schwarzschild radii of the more massive SMBH (Hughes, 2002), and neglect GW emission at smaller separations.
3 Predictions for the characteristic strain spectrum
As stated above, we consider four different initial eccentricity distributions: . Recall that the case corresponds to all binaries being circular. Initially circular binaries are not expected to become eccentric (e.g., Sesana, 2010). The probability mass functions of the binary eccentricities, , at in the three cases with are shown in Figure 1. For comparison, we also show in the bottom panel of Figure 1 a ‘thermal’ probability mass function for , derived from the probability density function for . This would be expected if binary systems followed a purely Maxwell-Boltzmann distribution of energies (e.g., Ambartsumian, 1937), as is roughly the case for galactic stellar binaries (Duquennoy & Mayor, 1991).
In Figure 2, we plot the characteristic strain spectra for each initial eccentricity distribution. Also depicted is the prediction in the circular (i.e., ), GW-driven case (i.e., for including only GW-driven orbital decay for all ). This latter prediction corresponds to the standard power-law. In order to help highlight the physical effects at work, Figure 3 shows the characteristic strain spectra for each assumed contributed by binaries with combined masses in the ranges and respectively.
The model we utilise for interactions between binaries and their stellar environments results in an attenuation of in the PTA frequency band compared to the power-law obtained in the circular, GW-driven case. For , the signal is attenuated at frequencies Hz. At these frequencies, stellar interactions are the dominant binary orbital decay process, increasing in Equation (12) and reducing the number of binaries observed per unit orbital frequency. For increasing , the signal is further attenuated at low frequencies, although a slight ( dex), increasing excess is present at frequencies between Hz and Hz. This is caused by two effects: eccentric binaries evolve faster than circular binaries, and eccentric binaries radiate GWs at higher harmonics of their orbital frequencies than circular binaries.
The ‘substructure’, or two bumps, in the characteristic strain spectra is a direct consequence of the mass-distribution of the binaries in our model. If were smooth and analytic, the characteristic strain spectra would have only one clear peak. Here, however, we evaluate this distribution from the G11 semi-analytic model outputs (see Equation (8)), which results in the distribution being incomplete at the high-mass end. These gaps in the distribution lead to the two apparent peaks in the characteristic strain spectra.
As is evident in Figure 3, the first peaks of the spectra in Figure 2 are dominated by the highest-mass binaries, whereas the second peaks are dominated by lower-mass binaries. This is because the evolution of the highest-mass binaries begins to be GW-driven at lower frequencies than for less massive binaries. There are expected to be very few binaries in the (combined) mass range ; only with Hz are expected to be present in the observable Universe according to the G11 model. In contrast, binaries are expected the the range . The effects of sparsity in are discussed further in §3.3.3.
3.2 Comparison with previous work
Our results for the GW characteristic strain spectra from an eccentric binary SMBH population are broadly consistent with similar previous studies (Enoki & Nagashima, 2007; Sesana, 2013a). Both these works find spectra which depart from the standard power law of the circular, GW-driven case at frequencies Hz. The results of Sesana (2013a) for binaries with eccentricities at formation of are in fact very similar to ours (see their Figure 2), with slight substructure evident along with the slight excess for Hz.
Our results for , however, differ somewhat from those of Sesana (2013a). Whereas the maximum separation between the zero-eccentricity and high-eccentricity curves (the red solid and dashed curves in Figure 2 of Sesana 2013a) is approximately 0.5 dex, the maximum difference between our curves for and in our Figure 2 is 0.35 dex. We also find similarly-shaped spectra for all , whereas Sesana (2013a) has a clear single peak in their zero-eccentricity curve.
The differences between our results and those of Sesana (2013a) for are caused by the nature of the respective binary SMBH mass distributions used. As discussed above, if the mass-distribution of binary SMBHs () were smooth and analytic, which is the case in Sesana (2013a), only a single peak is expected. The reason for the similarity between our results and those of Sesana (2013a) for non-zero eccentricities may be because of some discreteness in the eccentric binary SMBH distribution used by Sesana (2013a), as evidenced by the jagged nature of their strain spectrum at low frequencies.
The characteristic strain spectrum we predict in the circular, GW-driven case is dex lower than that predicted by Shannon et al. (2013). This difference is because we do not fit an analytic function to the discrete binary distribution . We discuss this point further in §3.3.3.
3.3 Uncertainties in the model predictions
In this section, we describe the key uncertainties in our prediction of , which are summarised in Figure 4. We consider in turn the accuracy of the model predictions for SMBH demographics and coalescence rates and for the rate of evolution of binary systems, and the effects of incomplete high-mass binary SMBH distributions.
SMBH demographics and coalescence rates
The merger rate of massive galaxies predicted by galaxy formation models (Bertone, De Lucia & Thomas 2007) implemented in the Millennium simulation (Springel et al., 2005) has been shown to be consistent with observational estimates at redshifts (Bertone & Conselice, 2009). Marulli et al. (2008) found that the model matches the observed quasar bolometric luminosity function at redshifts for a variety of assumed quasar lightcurves. This, together with the reproduction of the local SMBH-galaxy scaling relations, suggests that the rate of formation of massive binary SMBHs at low redshifts is satisfactorily reproduced by the G11 semi-analytic model, which is used as the basis for this paper. Furthermore, the characteristic strain spectrum expected in the case for binaries with combined masses at redshifts has a maximum disparity with the unrestricted spectrum of 0.02 dex. Hence, our model robustly predicts the contribution to the GW signal from massive, low-redshift binaries, which are likely to dominate the total GW signal (see also Wyithe & Loeb, 2003; McWilliams et al., 2012; Sesana, 2013b).
However, there remain a range of theoretical uncertainties. For example, the G11 model treatment of SMBHs does not include physically-motivated prescriptions for SMBH formation (e.g., Haiman, 2013), SMBH ejection caused by gravitational recoil following the coalescence of binary systems (e.g., Kulier et al., 2013), and does not account for any mass accreted onto SMBHs in merging galaxies prior to coalescence (e.g., Van Wassenhove et al., 2012).
There are also specific observational uncertainties in tuning the semi-analytic model. The current sample of SMBH and host galaxy bulge mass measurements, which is used to tune the quasar-mode SMBH accretion efficiency, allows for a confidence interval of dex in the SMBH masses (Shannon et al., 2013). Similarly, the galaxy stellar mass function predicted by the G11 model is matched to Sloan Digital Sky Survey observations in the nearby Universe (e.g., Li & White, 2009). These observations have a 0.2 dex systematic uncertainty, with negligible contribution from cosmic variance (Li & White, 2009), which corresponds (to first order) to a dex uncertainty in the galaxy merger rate.
The uncertainty in SMBH masses corresponds to a dex uncertainty in , while the uncertainty in the merger rate translates directly to the range of predictions for allowed by the observed galaxy stellar mass function. Combining both ranges results in a 0.4 dex () uncertainty in , which corresponds to a 0.2 dex uncertainty in .
The binary evolution model
In this paper, we assume that all galaxies hosting SMBHs have spherically-symmetric central stellar density profiles that are power-law functions of radius, , following Sesana (2010). That is, the stellar density, , is proportional to , where we have hitherto assumed . These profiles are equivalent to the central (asymptotic) behaviour of the Dehnen (1993) stellar potential and density models, which correspond well to high-resolution observations of the centres of nearby galaxy bulges (Faber et al., 1997). Our assumption of a universal is, however, not in agreement with observations, which typically show , with corresponding to the most extreme ‘core’ galaxies and corresponding to the most extreme ‘power-law’ galaxies (Dehnen, 1993; Faber et al., 1997). Furthermore, ‘core’ galaxies are generally more massive, early-type systems with more massive SMBHs, and ‘power-law’ galaxies are generally less massive, late-type systems with less massive SMBHs (e.g., Faber et al., 1997; McConnell & Ma, 2013). While we do not attempt to correlate with galaxy properties from the G11 model, we show in Figure 4 characteristic strain spectra in the case for and . The logarithmic differences between the spectra for these -values and the spectrum for may be applied only approximately to the spectra for other -values, because varying varies both the rate of semi-major axis decay and the rate of eccentricity evolution for binaries.
The model that we use (Sesana et al., 2006; Sesana, 2010) for binaries evolving through separations less than due to interactions with fixed stellar backgrounds is qualitatively similar to the results of recent numerical simulations of dry (i.e., free of dynamically significant gas) galaxy merger remnants (Khan et al., 2012). However, as we show in Appendix A, it is apparent that the model we use includes stronger stellar-driven orbital decay than the simulations of Khan et al. (2012). This is despite our assumption (see §2.2) that the rate of semi-major axis evolution is independent of binary eccentricity. This is not surprising, because the assumption of a fixed stellar background is qualitatively equivalent to the assumption of a full stellar loss-cone (cf. Quinlan & Hernquist, 1997; Sesana, 2010). Hence, the model we use maximises binary orbital decay rates, in particular for spherically symmetric stellar distributions.
We are likely therefore to be overestimating the effects of stellar interactions on the binary SMBH population. The numerical simulations of Khan et al. (2012) suggest that the orbital frequencies at which binary SMBH evolution begins to be predominantly GW-driven are up to 0.45 dex less than the corresponding frequencies that result from the model we use. This implies that the frequency below which the characteristic strain spectrum turns over from the power law may be up to 0.45 dex lower than we predict.
While the relation that we assume is established in the local Universe (Baes et al., 2003), it has not been studied at higher redshifts. Given the expected decrease in the stellar mass in a halo of a given mass with increasing redshift (Moster et al., 2010), it is possible that we are overestimating the velocity dispersions of the stellar cores of merger remnants beyond the local Universe. This would imply that higher-redshift binaries decay more slowly than in our model, again increasing the low-frequency parts of the presented characteristic strain spectra. Further work is required to quantify the magnitude of this increase.
Finally, while the assumed functional form of the -distribution (Equation (7)) is physically motivated, there may be some correlation between the orbital eccentricities of binaries with separations , and their masses and redshifts. Additionally, a variety of studies find physical reasons for binaries to be quite circular upon formation (Casertano et al., 1987; Polnarev & Rees, 1994; Hashimoto et al., 2003), which suggests that low- values may be preferred. Current numerical simulations (e.g., Khan et al., 2012) have not been run with a sufficient range of initial conditions to provide conclusive results on this point.
Accounting for discreteness in the binary SMBH distributions from the G11 model
Given the distribution (Equation (8), §2.2) of binary SMBHs, we have examined the expected value of the GW characteristic strain spectrum. However, the binary SMBH distribution that we use is not exactly the distribution expected from the G11 semi-analytic model implemented within the Millennium simulation. The Millennium simulation provides a single realisation of the dark matter halo merger history within a large comoving volume, and the G11 prescriptions specify properties of the galaxies, and SMBHs, associated with the halos. To form the discrete binary distribution , we count the numbers of binary SMBHs forming in the entire Millennium volume between redshift snapshots in bins of and , assigning values of to each binary according to our -distribution. However, despite the large volume, is poorly populated for high and at every redshift. To estimate the expected nature of this distribution, statistical modelling is required. This was done by Shannon et al. (2013), who considered the circular, GW-driven case, and found that the modelling resulted in a characteristic strain spectrum increased by 0.15 dex. However, the distribution has two more dimensions (an extra mass dimension and the eccentricity dimension) than that considered by Shannon et al. (2013), which significantly complicates the modelling. Instead, we simply consider it possible that the strain spectrum we have derived may be up to 0.15 dex larger.
A qualitatively similar effect was pointed out by (Sesana, Vecchio & Colacino 2008), who compared characteristic strain spectra generated from realisations of the binary SMBH population of the Universe to the spectrum expected on average, in the circular, GW-driven case. Whereas the average spectrum was a power-law proportional to , individual realisations had a lower amplitude at higher frequencies. This was because the numbers of binaries radiating GWs at a given frequency (per unit frequency) decreases with increasing frequency, implying that, for example, there is a frequency above which the expected number of sources is less than unity. However, the correct model for the average characteristic strain spectrum still had for all , despite all realisations of the spectrum being below this power-law at high frequencies. This situation is analogous to our suggestion of an increase in the characteristic strain spectrum if the average behaviour of were correctly modelled.
We also do not attempt here to describe the statistical nature of the GW signal, as was done by Ravi et al. (2012) in the circular, GW-driven case. Ravi et al. (2012) modelled a GW signal that was mildly non-Gaussian, with individual sources dominant at all GW frequencies of interest to PTAs. Shannon et al. (2013) further showed that assuming non-Gaussian statistics for the GW signal caused constraints on to degrade by . This reflects the fact that realisations of at a particular frequency would have a larger variance in the non-Gaussian case than in the Gaussian case.
As discussed in §3.1, environment-driven binary evolution causes the highest-mass binaries to dominate at low frequencies to a greater extent than in the purely GW-driven case. This, coupled with the sparsity of these binaries in our calculations, causes the low-frequency substructure in the characteristic strain spectra for all in Figure 2. Our results, however, suggest a more general conclusion: that, at low frequencies, environment-driven binary evolution causes the variance in realisations of to be significantly increased relative to the assumption of only GW-driven evolution. Including this increased variance in at low frequencies in the calculation of PTA upper bounds on (e.g., Shannon et al., 2013) would cause these constraints to be further degraded relative to constraints based on the work of Ravi et al. (2012).
Synthesis of uncertainties in
We refer the reader to Figure 4, where we show an approximate confidence interval on the characteristic strain spectrum according to the model we describe. This interval represents our uncertainty in the expected value of the signal, not the realisation-to-realisation uncertainty. The interval encompasses the maximum possible ranges of and (see §3.3.2), and also includes observational uncertainties in the SMBH-bulge mass relation and in the galaxy stellar mass function (see §3.3.1). We also include our assertion that the characteristic strain spectrum could be up to 0.15 dex larger than what we calculate if the binary SMBH distribution were correctly specified (see §3.3.3).
It is clear that that there is relatively more uncertainty in our prediction at frequencies Hz, where environmental interactions and binary eccentricities may affect the signal. We have also not included our uncertainty in the specific model for environment-driven binary SMBH evolution. As discussed in §3.3.2, the model we use may represent the maximum level of binary-environment coupling; other models may result in the characteristic strain spectrum being boosted at low frequencies relative to our prediction. For example, the model of Khan et al. (2012) suggests that the effects of environmental interactions may only be relevant for Hz (also see Appendix A). We have also weighted each -value equally, whereas it is possible that low- values are preferred over high- values.
In Figure 4, we also indicate the best upper bound on a stochastic, Gaussian GW background from binary SMBHs, published recently by the Parkes Pulsar Timing Array (PPTA; Shannon et al., 2013). This upper bound corresponds to with 95% confidence. While PTA bounds are traditionally shown as wedges (e.g., Sesana et al., 2008, Figure 13) on characteristic strain spectrum plots, Shannon et al. (2013) argued that their limit was applicable only at a single GW frequency. We hence display this limit as a single dot.
Our prediction for the characteristic strain spectrum at a frequency of of (with approximately 68% confidence) is broadly consistent with previous results (e.g., Wyithe & Loeb, 2003; Sesana et al., 2008; Sesana, 2013b; Shannon et al., 2013) that considered the circular, GW-driven case. Indeed, for Hz, the predicted characteristic strain spectrum closely resembles the power law expected in the circular, GW-driven case, with the exception that for larger slightly more signal is present. The departure from the power law at Hz is caused by binary SMBHs radiating at these frequencies reaching their last stable orbits and not being included in our calculations (see, e.g., Wyithe & Loeb, 2003).
3.4 Summary of PTA implications
Our results suggest a challenging future for attempts at detecting the GW background from binary SMBHs with PTAs. The frequency of optimal sensitivity for PTAs generally corresponds to the inverse of the characteristic observation time (e.g., Shannon et al., 2013). Typical observation times of yr (Manchester et al., 2013) imply that the properties of the GW signal at frequencies in the range Hz to Hz are of primary importance for PTA work. The model we use in this paper implies that the GW characteristic strain spectrum may be reduced throughout this frequency range relative to the circular, GW-driven case (i.e., with respect to a power law). For example, Shannon et al. (2013) presented a single-frequency constraint on that is inconsistent with a variety of astrophysical predictions assuming circular, GW-driven binaries. However, the constraint is at a frequency where the characteristic strain spectrum we predict (see Figure 4) is reduced by at least 0.08 dex relative to the circular, GW-driven case. More generally, the gains in sensitivity to a GW background with observing time, estimated assuming (e.g., Siemens et al., 2013), are likely to be overestimated. Furthermore, it is possible that at low frequencies the increased contribution to the total GW signal of rare, massive binary SMBHs relative to the circular, GW case will cause the signal at these frequencies to be more non-Gaussian than suggested by Ravi et al. (2012).
However, our results require significant refinement. It is clear from Figure 4 that our prediction for the characteristic strain spectrum at low frequencies is quite uncertain. Besides this uncertainty, the model we use in this paper for the coupling between binary SMBHs and stellar environments (Sesana et al., 2006; Sesana, 2010) may in fact maximise the strength of this coupling (see §3.3.2 and Appendix A). Also, it is possible that lower-eccentricity scenarios may be preferred over the higher-eccentricity scenarios. Both the above possibilities would result in the low-frequency parts of the characteristic strain spectrum being increased relative to our predictions. We strongly urge further work on modelling the evolution of binary SMBH orbits in a variety of realistic galaxy merger scenarios. This is of significant importance for predicting the strength of the GW signal from binary SMBHs in the PTA frequency band.
4 Predictions for GW bursts
4.1 The distribution of GW bursts
The prospect of detecting GW bursts with PTAs has been pursued recently by a number of authors (e.g., Finn & Lommen, 2010; Pitkin, 2012). GW burst detection algorithms generally contain few assumptions about the source properties, except that they search for a strong signal confined to a short time-period. Here, we focus on the properties of GW bursts from eccentric binary SMBHs, and use our distribution of binary SMBHs, , to predict the distribution of burst events.
It is necessary to form a definition of a GW burst from an eccentric binary SMBH in terms pulsar timing data products. Pulsar timing is the practice of measuring the times of arrival (ToAs) of pulses from millisecond radio pulsars and fitting a physical model to these measurements. In Appendix B, we describe how GWs from eccentric binaries affect pulsar timing measurements by inducing variations to ToAs. We present an expression for the rms deviation, , of the ToA variations as a function of time caused by a given binary SMBH in Equation B11, averaged over all orientation parameters.
In Figure 5, we show the orbital phase, , the expected energy flux in GWs at the Earth and as functions of time for a binary SMBH with eccentricity 0.8 and orbital period at the Earth of 3.1 yr at the starting time, component masses and , and redshift 0.1. We also show the induced ToA variations corresponding to the GW metric perturbation at the Earth, , for arbitrary orientation parameters (binary inclination rad, line of nodes orientation rad). The orbital evolution of the binary was calculated using the work of Peters & Mathews (1963), and the energy flux at the Earth is averaged over binary inclination. The time-intervals considered to be GW bursts are highlighted in all panels of Figure 5. These ‘bursts’ correspond to the times of the largest change in the shortest amount of time in the ToA variations (see the top panel of Figure 5), and can be identified using . We define the true burst amplitude, , to be twice the peak value of , because that represents the expected peak-to-peak variation for a burst. The burst duration, , is the time-interval between peaks in , represented by the widths of the shaded intervals in Figure 5. It is interesting that the GW bursts in the ToAs correspond to the motion of the binary through apastron, rather than periastron. This is because the GW-induced ToA variations are given by the time-integral of the GW amplitude as a function of time, as outlined in Appendix B.
The qualitative properties of in Figure 5 apply to binaries with any component masses, orbital period and eccentricity. That is, there are two peaks per rotation period, separated by less in orbital phase for more eccentric binaries, and separated by half an orbital phase for circular binaries. For a binary specified by and , we integrate the equations for the evolution of the orbit (Peters & Mathews, 1963) from zero orbital phase to numerically calculate as the mean of the first two peaks in , and as the time-interval between peaks.
We use the distribution of binary SMBHs, , to calculate the distribution of GW bursts. As described above, this distribution is specified as the number of binaries in discrete bins of width , , , and , where the eccentricity bin-widths depend on the other parameters. Scaling this distribution by the comoving volume shell between redshifts and specifies the number of observable binary SMBHs. For parameters at the midpoints of each bin, we calculate and . We approximate the burst rate from binaries in a bin as the number of binaries divided by their period observed at the Earth, and record the expected number of bursts in a 10 yr time-span.
Using the distributions of binary SMBHs for and , we calculated the numbers of GW bursts for different values of the expected maximum level of ToA variations, , and the duration, . We depict the distributions of GW bursts in Figure 6 as the number of bursts per 10 yr observation time, , per dex, in bins of 0.075 dex in and 0.05 dex in . We only considered bursts with ns and 0.1 yr 10 yr. An rms ToA variation of 40 ns corresponds to the best timing precisions currently achieved for millisecond radio pulsars (Osłowski et al., 2013; Hobbs, 2013).
In total, we predict 0.06 bursts per 10 yr in the case with these strengths and durations, as compared with 0.12 bursts per 10 yr in the case. This difference in the total number of bursts is because of the smaller number of binary SMBHs that we expect to observe if the population is generally more eccentric. However, we note that bursts from low-eccentricity binaries, which will dominate the burst population in the case, may be less detectable than bursts from high-eccentricity binaries. There are proportionally more short-duration bursts in the case than in the case, because larger binary eccentricities result in shorter bursts.
In both cases, the burst distribution is quite heavily skewed towards long bursts, with approximately a factor of 100 more bursts expected with yr durations than with yr durations. There are also fewer bursts with durations longer than yr in both cases. This typical burst duration corresponds to binaries with separations where GW-driven evolution is equivalent to evolution driven by stellar environments.
The typical combined masses of the binary SMBHs that produce GW bursts are . In Figure 7, we show the distributions of the combined masses of all binary SMBHs producing the bursts in the distributions in Figure 6. The distributions are similar in shape, although the distribution for includes relatively more high-mass binaries than the distribution for . This is because, in the case, lower-mass binaries are less likely to be able to produce strong GW bursts because they are likely to be more eccentric. More eccentric binaries of a given mass and orbital period produce typically weaker bursts (see Equation (B11) in Appendix B).
In summary, we find:
For bursts with durations between 0.1 yr and 10 yr, and with expected maximum ToA variations of ns, we predict between 0.06 and 0.12 bursts per 10 yr observation, with lower burst rates corresponding to higher-eccentricity binary SMBH populations.
Higher-eccentricity binary populations result in relatively more shorter duration bursts than lower-eccentricity populations.
However, the burst rate decreases by a factor of 10 per dex below a duration of yr. This also appears to be the most likely duration, with few bursts longer than 8 yr expected.
The burst rate decreases by a factor of 10 per dex increase in amplitude.
Various uncertainties discussed in §3.3 also apply to these calculations. The uncertainty in the galaxy merger rate will also directly correspond to the uncertainty in the GW burst rate (i.e., 0.3 dex). Given that the high-end power-law logarithmic slope of the SMBH mass function in the G11 model is , the 0.2 dex uncertainty in the SMBH masses will, to first order, correspond to an uncertainty of 0.4 dex in the merger rate. Therefore, the uncertainty in the burst rate from the model is approximately 0.5 dex.
To our knowledge, only one study has attempted to predict the properties of GW bursts from a population of eccentric binary SMBHs (Liu et al., 2012). While our modelling methods and definition of GW bursts differ substantially from this work, we agree with these authors that it is unlikely that current PTAs will be able to detect GW bursts from binary SMBHs. The rarity of short GW bursts (lasting around 1 yr) from binary SMBHs suggests that high-cadence PTA observations targeting such bursts are not well motivated.
5 Summary of results
We have used a semi-analytic model for galaxy and SMBH formation and evolution (Guo et al., 2011) implemented in the Millennium simulation (Springel et al., 2005), augmented with a model for the evolution of binary SMBHs within fixed stellar backgrounds (Sesana et al., 2006; Sesana, 2010), to predict the properties of low-frequency GWs from binary SMBHs. We specify the form of a phenomenological distribution of initial binary eccentricities, and consider a selection of cases with differing levels of typical binary eccentricity.
Our quantitative results are uncertain due to a variety of factors. The range of initial binary eccentricity distributions that we consider corresponds to a 0.4 dex variation in the characteristic strain spectrum at low frequencies. Moreover, uncertainties in the tuning of the G11 model provide another 0.2 dex of uncertainty in the spectrum at all frequencies. There is also uncertainty in our estimate of the binary SMBH distribution predicted by the G11 model, in particular for the most massive binaries. Finally, while the G11 model is likely to provide a satisfactory representation of the merger rate of massive, low-redshift galaxies, the binary evolution model that we use may overestimate binary hardening caused by stellar interactions.
Our specific findings are as follows:
The expected characteristic strain spectrum of the GW background from binary SMBHs will turn over from the standard power law at a frequency up to Hz. The turn-over frequency depends on the efficiency of stellar interactions in extracting energy and angular momentum from binary SMBHs, as well as the typical binary eccentricities at formation.
The nature of the spectrum at frequencies below the turn-over frequency is extremely uncertain, and depends on the numbers of massive () binaries and on binary eccentricities. The most massive binary SMBHs predominantly produce the lowest-frequency parts of the spectrum, and their numbers depend strongly on the strength of their coupling to their environments. The spectrum will be attenuated if binaries with typically larger eccentricities are present.
The most massive eccentric binaries will very rarely produce GW bursts detectable in pulsar timing data. A larger-eccentricity binary population will produce fewer bursts that are typically shorter and weaker. Our results suggest that GW bursts from binary SMBHs do not provide viable targets for PTA observations.
We emphasise a set of key implications of our work for PTAs:
Given the expected low-frequency turn-over in the GW characteristic strain spectrum, along with the large uncertainty in the signal at these frequencies, the increase with time of PTA sensitivities to a GW background from binary SMBHs will not be as strong as currently thought.
Short-duration, strong GW bursts from eccentric binary SMBHs are unlikely to occur during typical PTA dataset lengths.
PTA data analysts cannot assume when searching for a GW background from binary SMBHs. Indeed, model-independent searches cannot assume any particular spectral shape.
Model-dependent searches and constraints need to carefully account for the uncertainty in model predictions.
In this paper, we predict both the GW background characteristic strain spectrum and the distribution of strong GW bursts from eccentric binaries. At a GW frequency of (1 yr), we predict a characteristic strain of with approximately 68% confidence.
Accelerated binary evolution driven by three-body stellar interactions causes the characteristic strain spectrum to be diminished with respect to a power-law at Hz. At these low frequencies, the signal is further attenuated if binary SMBHs are typically more eccentric at formation. The low-frequency signal may be dominated by a few binaries with combined masses () greater than , to a larger extent than predicted in the circular, GW-driven case (Ravi et al., 2012). Numerous uncertainties, however, affect our results. These include observational uncertainties in parameters of our model, and theoretical uncertainties in the efficiency of coupling between binary SMBHs and their environments.
We also expect between 0.06 and 0.12 GW bursts that produce 40 ns amplitude ToA variations over a 10 yr observation time. Larger typical binary eccentricities at formation will result in fewer events than if binaries are less eccentric at formation. These bursts are caused by binary SMBHs with combined masses of , and typically last yr. Shorter, stronger bursts are significantly less likely, as are longer bursts.
Upcoming radio telescopes with extremely large collecting areas, such as the Five hundred metre Aperture Spherical Telescope (FAST, Li, Nan & Pan 2013) and the Square Kilometre Array (SKA, Cordes et al., 2004) are likely to significantly expand the sample of pulsars with sufficient timing precision for GW detection as compared to current instruments. PTAs formed with FAST and the SKA will hence be sensitive to a stochastic GW signal at much higher frequencies than current PTAs, which is desirable given the results we present.
The mechanism by which binary SMBHs are driven to the GW-dominated regime must involve some form of binary-environment coupling. Hence, independent of the exact model, there will always be some low-frequency attenuation of the GW signal relative to the circular, GW-driven binary case. Our results indicate that this attenuation occurs within the PTA frequency band. However, the strength of the binary-environment coupling is quite uncertain, and we urge future work on this topic. Finally, as also emphasised in previous works (Enoki & Nagashima, 2007; Sesana, 2013a), constraining or measuring the spectrum of the GW background at a number of frequencies would provide an excellent test of models for the binary SMBH population of the Universe.
The authors thank Alberto Sesana, Sarah Burke-Spolaor, Yuri Levin, Simon Mutch and Jonathan Khoo for useful discussions. V.R. is a recipient of a John Stocker Postgraduate Scholarship from the Science and Industry Endowment Fund, and J.S.B.W. acknowledges an Australian Research Council Laureate Fellowship. The Millennium and Millennium-II Simulation databases used in this paper and the web application providing online access to them were constructed as part of the activities of the German Astrophysical Virtual Observatory. This work was performed on the swinSTAR supercomputer at the Swinburne University of Technology.
Appendix A Testing our implementation of a binary SMBH evolution model
In the main text, we use the results of Sesana et al. (2006) and Sesana (2010) (hereafter collectively referred to as S06) to model the evolution of binary SMBHs in fixed stellar backgrounds for separations less than (see Equation (6)). S06 numerically solved three-body scattering problems for binary SMBHs interacting with stars on radial, intersecting orbits drawn from a spherically-symmetric, fixed distribution, and provided fitting formulae for the binary hardening and eccentricity growth rates as functions of binary properties. We use these fitting formulae to evolve binary SMBH orbits as described in §2.3.
Here, we compare this method of evolving binary SMBHs with recent -body simulations of binary SMBH evolution in merging galaxies of various mass ratios and stellar density distributions (Khan et al., 2012, hereafter K12). K12 simulated the mergers of spherical galaxies with various mass ratios, power-law stellar cusp density profiles with various indices, with typical approach trajectories from cosmological simulations. The SMBHs were traced until separations close to, and in some cases beyond, where the GW-driven orbital decay dominated the orbital decay caused by three-body stellar interactions. By extrapolating the binary orbits assuming constant stellar-driven hardening rates and eccentricities, K12 estimated binary semi-major axes, , below which GW-driven evolution dominated. The accuracy of these extrapolations was confirmed using a selection of simulations including post-Newtonian corrections to the binary SMBH orbits.
Here, we take the final eccentricities and semi-major axes of the binaries in each of the scenarios considered by K12 for which was estimated, and evolve the binaries using the binary evolution model of S06 to estimate an equivalent quantity to , . We list the ratios in Table A1 both without (column 4) and with the binary eccentricity held fixed (column 5) for each relevant scenario of K12. The cusp density profile indices, , and the galaxy and SMBH mass ratios, , are given in columns 2 and 3 respectively.
We find in all cases. This implies that the S06 model that we use in our work has stronger stellar-driven binary evolution than the K12 model. We also have smaller ratios when we hold the binary eccentricities fixed. This is because the binary eccentricities invariably grow when allowed to evolve, and lower eccentricities imply smaller GW-driven hardening rates. While an intuitive explanation of the difference between the S06 and K12 models is difficult to attain, the K12 work involves a more sophisticated, and possibly more realistic, treatment of the distribution of stellar orbits in the cores of merged galaxies than the S06 work. Given (Equation (9)), a difference of a factor of two in the semi-major axes at which binary SMBH evolution begins to be GW-dominated corresponds to a difference of a factor of in orbital frequency.
Appendix B Effects on pulsar timing measurements of GW bursts from a binary SMBH
Here, we provide a mathematical description of GW bursts from eccentric binary SMBHs. The spatial metric perturbation tensor, or GW strain, corresponding to a binary SMBH was defined by Wahlquist (1987) to lowest order in the slow-motion, far-field regime using the quadrupole formula. This tensor, , can be written as (Wahlquist, 1987)
where and are the ‘plus’ and ‘cross’ polarisation tensors respectively, and and are the time-varying polarisation amplitudes, which depend on the orbital phase (which is a function of time), the value of at the line of nodes (), the value of at periastron (), the orientation of the line of nodes (), the binary inclination (), and a factor
where is the comoving coordinate distance to redshift and is the binary chirp mass.
GWs at the Earth and at a pulsar cause a fractional shift in the observed pulsar rotation frequency. Here, we neglect the effects of GWs at the pulsar, because, as outlined in, e.g., Finn & Lommen (2010), GW bursts will generally affect pulsar timing data at vastly different times for different pulsars. This means that an approach that seeks to detect GW bursts by observing correlated effects in multiple pulsar datasets will only need to consider the effects of GWs at the Earth. For a pulsar with location defined by the unit direction tensor , the observed fractional pulsar rotation frequency shift is given by (Wahlquist, 1987; Hobbs et al., 2009)
where is the cosine of the angle between the pulsar and GW source directions, and we follow the Einstein summation convention over the tensor indices.
Fractional shifts in will cause cumulative variations in ToAs. That is,
where is the ToA variation at a time . In order to calculate the average duration and strength of a GW burst from a binary SMBH as manifested in , we need to calculate the variance
where the angle brackets signify averaging over the subscripted quantities, and and are the right ascension and declination of the pulsar assuming that the GW source is located along the z-axis. To simplify this, we set , and express as
The linear independence of the polarisation tensors implies that
We find that
Then, we have
To summarise, Equation (B11) gives the variance of the ToA variations caused by GWs from an individual eccentric binary SMBH at the Earth, averaged over all orientation parameters.
- pagerange: 1–15
- pubyear: 2014
- We refer to this as the “circular, GW-driven case” throughout the paper.
- See Gaddum (1945) for a discussion of the ubiquity of log-normal distributions in nature.
- Ambartsumian V. A., 1937, Astronomicheskii Zhurnal, 14, 207
- Baes M., Buyle P., Hau G. K. T., Dejonghe H., 2003, MNRAS, 341, L44
- Baker J. G., Centrella J., Choi D.-I., Koppitz M., van Meter J. R., Miller M. C., 2006, ApJ, 653, L93
- Begelman M. C., Blandford R. D., Rees M. J., 1980, Nat, 287, 307
- Berentzen I., Preto M., Berczik P., Merritt D., Spurzem R., 2009, ApJ, 695, 455
- Bertone S., Conselice C. J., 2009, MNRAS, 396, 2345
- Bertone S., De Lucia G., Thomas P. A., 2007, MNRAS, 379, 1143
- Callegari S., Mayer L., Kazantzidis S., Colpi M., Governato F., Quinn T., Wadsley J., 2009, ApJ, 696, L89
- Casertano S., Phinney E. S., Villumsen J. V., 1987, in de Zeeuw P. T., Tremaine S. D., eds, Structure and Dynamics of Elliptical Galaxies Vol. 127 of IAU Symposium, Dynamical Friction and Orbit Circularization. p. 475
- Cordes J. M., Kramer M., Lazio T. J. W., Stappers B. W., Backer D. C., Johnston S., 2004, New Astronomy Review, 48, 1413
- Dehnen W., 1993, MNRAS, 265, 250
- Dotti M., Sesana A., Decarli R., 2012, Advances in Astronomy, 2012
- Duquennoy A., Mayor M., 1991, A&A, 248, 485
- Ellis J. A., Siemens X., Creighton J. D. E., 2012, ApJ, 756, 175
- Enoki M., Nagashima M., 2007, Progress of Theoretical Physics, 117, 241
- Escala A., Larson R. B., Coppi P. S., Mardones D., 2004, ApJ, 607, 765
- Faber S. M., Tremaine S., Ajhar E. A., Byun Y.-I., Dressler A., Gebhardt K., Grillmair C., Kormendy J., Lauer T. R., Richstone D., 1997, AJ, 114, 1771
- Finn L. S., Lommen A. N., 2010, ApJ, 718, 1400
- Foster R. S., Backer D. C., 1990, ApJ, 361, 300
- Frank J., Rees M. J., 1976, MNRAS, 176, 633
- Gaddum J. H., 1945, Nat, 156, 3964
- Guo Q., White S., Boylan-Kolchin M., De Lucia G., Kauffmann G., Lemson G., Li C., Springel V., Weinmann S., 2011, MNRAS, 413, 101
- Haiman Z., 2013, in Wiklind T., Mobasher B., Bromm V., eds, Astrophysics and Space Science Library Vol. 396 of Astrophysics and Space Science Library, The Formation of the First Massive Black Holes. p. 293
- Hashimoto Y., Funato Y., Makino J., 2003, ApJ, 582, 196
- Hobbs G., 2013, Classical and Quantum Gravity, 30, 224007
- Hobbs G., Jenet F., Lee K. J., Verbiest J. P. W., Yardley D., Manchester R., Lommen A., Coles W., Edwards R., Shettigara C., 2009, MNRAS, 394, 1945
- Hughes S. A., 2002, MNRAS, 331, 805
- Khan F. M., Preto M., Berczik P., Berentzen I., Just A., Spurzem R., 2012, ApJ, 749, 147
- Kulier A., Ostriker J. P., Natarajan P., Lackner C. N., Cen R., 2013, ArXiv e-prints
- Li C., White S. D. M., 2009, MNRAS, 398, 2177
- Li D., Nan R., Pan Z., 2013, in IAU Symposium Vol. 291 of IAU Symposium, The Five-hundred-meter Aperture Spherical radio Telescope project and its early science opportunities. pp 325–330
- Liu J., Zhang Y., Zhang H., Sun Y., Wang N., 2012, A&A, 540, A67
- Manchester R. N., Hobbs G., Bailes M., Coles W. A., van Straten W., Keith M. J. e. a., 2013, PASA, 30, 17
- Marulli F., Bonoli S., Branchini E., Moscardini L., Springel V., 2008, MNRAS, 385, 1846
- McConnell N. J., Ma C.-P., 2013, ApJ, 764, 184
- McWilliams S. T., Ostriker J. P., Pretorius F., 2012, ArXiv e-prints
- Moster B. P., Somerville R. S., Maulbetsch C., van den Bosch F. C., Macciò A. V., Naab T., Oser L., 2010, ApJ, 710, 903
- Osłowski S., van Straten W., Demorest P., Bailes M., 2013, MNRAS, 430, 416
- Peters P. C., Mathews J., 1963, Physical Review, 131, 435
- Phinney E. S., 2001, ArXiv Astrophysics e-prints
- Pitkin M., 2012, MNRAS, 425, 2688
- Polnarev A. G., Rees M. J., 1994, A&A, 283, 301
- Quinlan G. D., 1996, New Astronomy, 1, 35
- Quinlan G. D., Hernquist L., 1997, New Astronomy, 2, 533
- Ravi V., Wyithe J. S. B., Hobbs G., Shannon R. M., Manchester R. N., Yardley D. R. B., Keith M. J., 2012, ApJ, 761, 84
- Roedig C., Dotti M., Sesana A., Cuadra J., Colpi M., 2011, MNRAS, 415, 3033
- Sesana A., 2010, ApJ, 719, 851
- Sesana A., 2013a, ArXiv e-prints
- Sesana A., 2013b, MNRAS, 433, L1
- Sesana A., Haardt F., Madau P., 2006, ApJ, 651, 392
- Sesana A., Haardt F., Madau P., 2008, ApJ, 686, 432
- Sesana A., Vecchio A., Colacino C. N., 2008, MNRAS, 390, 192
- Sesana A., Vecchio A., Volonteri M., 2009, MNRAS, 394, 2255
- Shannon R. M., Ravi V., Coles W. A., Hobbs G., Keith M. J., Manchester R. N., Wyithe J. S. B., Bailes M., Bhat N. D. R., Burke-Spolaor S., Khoo J., Levin Y., Oslowski S., Sarkissian J. M., van Straten W., Verbiest J. P. W., Want J.-B., 2013, Science, 342, 334
- Siemens X., Ellis J., Jenet F., Romano J. D., 2013, ArXiv e-prints
- Springel V., White S. D. M., Jenkins A., Frenk C. S., Yoshida N., Gao L., Navarro J., Thacker R., Croton D., Helly J., Peacock J. A., Cole S., Thomas P., Couchman H., Evrard A., Colberg J., Pearce F., 2005, Nat, 435, 629
- van den Bosch F. C., Lewis G. F., Lake G., Stadel J., 1999, ApJ, 515, 50
- van Haasteren R., Levin Y., Janssen G. H., Lazaridis K., Kramer M., Stappers B. W. e. a., 2011, MNRAS, 414, 3117
- Van Wassenhove S., Volonteri M., Mayer L., Dotti M., Bellovary J., Callegari S., 2012, ApJ, 748, L7
- Wahlquist H., 1987, General Relativity and Gravitation, 19, 1101
- Wetzel A. R., 2011, MNRAS, 412, 49
- Wyithe J. S. B., Loeb A., 2003, ApJ, 590, 691
- Yardley D. R. B., Coles W. A., Hobbs G. B., Verbiest J. P. W., Manchester R. N., van Straten W., Jenet F. A., Bailes M., Bhat N. D. R., Burke-Spolaor S., Champion D. J., Hotan A. W., Oslowski S., Reynolds J. E., Sarkissian J. M., 2011, MNRAS, 414, 1777
- Yardley D. R. B., Hobbs G. B., Jenet F. A., Verbiest J. P. W., Wen Z. L., Manchester R. N., Coles W. A., van Straten W., Bailes M., Bhat N. D. R., Burke-Spolaor S., Champion D. J., Hotan A. W., Sarkissian J. M., 2010, MNRAS, 407, 669
- Yu Q., 2002, MNRAS, 331, 935
- Yu Q., Lu Y., Mohayaee R., Colin J., 2011, ApJ, 738, 92