Binary supermassive black hole environments diminish the gravitationalwave signal in the pulsar timing band
Abstract
We assess the effects of supermassive black hole (SMBH) environments on the gravitationalwave (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 nonzero 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, GWdriven 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 SMBHgalaxy scaling relations and the galaxy stellar mass function, uncertainty in the nature of binaryenvironment coupling, and uncertainty in the numbers of the most massive binary SMBHs. We conclude, however, that lowfrequency 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 analysis1 Introduction
The merger of a pair of galaxies hosting central supermassive 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 circumnuclear gas disk (e.g., Roedig et al., 2011). If the orbital decay process can drive the binary to a small separation, gravitationalwave (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 selfgravitating gas disks will have large eccentricities, between 0.6 and 0.8, at the onset of GWdominated evolution. In the GWdominated 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 GWdominated 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 wellknown 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 semianalytic 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 gasdriven binary evolution. This is because massive galaxies at low redshifts, which are expected to dominate the GW signal from binary SMBHs, will typically be latetype and gaspoor (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 GWdriven 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, GWdriven 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 restframe frequency, , of . The specific intensity of GWs at the Earth from sources between redshifts and is
(1) 
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 rearranging, the energy density in GWs at the Earth per logarithmic frequency unit is, in any homogeneous and isotropic universe,
(2)  
(3)  
(4) 
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
(5) 
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 semimajor 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
(6) 
We assume that dynamical friction is effective in driving binaries to mean separations (e.g., Callegari et al., 2009; Khan et al., 2012), and consider binaries at this stage to be newly formed.
We adopt a simple, oneparameter distribution for the eccentricities of binary SMBHs at formation,
based on the postulate that the semimajor and semiminor axes of the orbits ( and respectively) are each lognormally 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
(7) 
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 semimajor axes .
We use the results of the semianalytic 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
(8) 
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 recentlymerged 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 powerlaw 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 powerlaw 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 semimajor 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 semimajor axes and eccentricities based on numerical scattering experiments (the ’’ and ’’ coefficients respectively from Tables 1 and 3 of Sesana et al. 2006). We loginterpolate the published values at binary component mass ratios of interest. As Sesana et al. (2006) only provide rates of semimajor axis evolution for circular binaries, we assume here that the rate of semimajor axis evolution at a given semimajor axis is independent of eccentricity. This approximation leads to the semimajor 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 restframe 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
(9) 
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 GWemitting 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):
(10) 
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
(11)  
(12) 
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
(13) 
Recall that . For consistency with other works, we calculate the characteristic strain spectrum, defined as
(14) 
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
3.1 Results
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 MaxwellBoltzmann 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., ), GWdriven case (i.e., for including only GWdriven orbital decay for all ). This latter prediction corresponds to the standard powerlaw. 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 powerlaw obtained in the circular, GWdriven 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 massdistribution 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 semianalytic model outputs (see Equation (8)), which results in the distribution being incomplete at the highmass 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 highestmass binaries, whereas the second peaks are dominated by lowermass binaries. This is because the evolution of the highestmass binaries begins to be GWdriven 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, GWdriven 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 zeroeccentricity and higheccentricity 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 similarlyshaped spectra for all , whereas Sesana (2013a) has a clear single peak in their zeroeccentricity 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 massdistribution 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 nonzero 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, GWdriven 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 highmass 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 SMBHgalaxy scaling relations, suggests that the rate of formation of massive binary SMBHs at low redshifts is satisfactorily reproduced by the G11 semianalytic 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, lowredshift 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 physicallymotivated 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 semianalytic model. The current sample of SMBH and host galaxy bulge mass measurements, which is used to tune the quasarmode 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 sphericallysymmetric central stellar density profiles that are powerlaw 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 highresolution 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 ‘powerlaw’ galaxies (Dehnen, 1993; Faber et al., 1997). Furthermore, ‘core’ galaxies are generally more massive, earlytype systems with more massive SMBHs, and ‘powerlaw’ galaxies are generally less massive, latetype 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 semimajor 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 stellardriven orbital decay than the simulations of Khan et al. (2012). This is despite our assumption (see §2.2) that the rate of semimajor 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 losscone (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 GWdriven 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 higherredshift binaries decay more slowly than in our model, again increasing the lowfrequency 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 semianalytic 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, GWdriven 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, GWdriven case. Whereas the average spectrum was a powerlaw 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 powerlaw 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, GWdriven case. Ravi et al. (2012) modelled a GW signal that was mildly nonGaussian, with individual sources dominant at all GW frequencies of interest to PTAs. Shannon et al. (2013) further showed that assuming nonGaussian 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 nonGaussian case than in the Gaussian case.
As discussed in §3.1, environmentdriven binary evolution causes the highestmass binaries to dominate at low frequencies to a greater extent than in the purely GWdriven case. This, coupled with the sparsity of these binaries in our calculations, causes the lowfrequency substructure in the characteristic strain spectra for all in Figure 2. Our results, however, suggest a more general conclusion: that, at low frequencies, environmentdriven binary evolution causes the variance in realisations of to be significantly increased relative to the assumption of only GWdriven 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 realisationtorealisation uncertainty. The interval encompasses the maximum possible ranges of and (see §3.3.2), and also includes observational uncertainties in the SMBHbulge 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 environmentdriven binary SMBH evolution. As discussed in §3.3.2, the model we use may represent the maximum level of binaryenvironment 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, GWdriven case. Indeed, for Hz, the predicted characteristic strain spectrum closely resembles the power law expected in the circular, GWdriven 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, GWdriven case (i.e., with respect to a power law). For example, Shannon et al. (2013) presented a singlefrequency constraint on that is inconsistent with a variety of astrophysical predictions assuming circular, GWdriven 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, GWdriven 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 nonGaussian 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 lowereccentricity scenarios may be preferred over the highereccentricity scenarios. Both the above possibilities would result in the lowfrequency 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 timeperiod. 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 timeintervals 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 peaktopeak variation for a burst. The burst duration, , is the timeinterval 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 GWinduced ToA variations are given by the timeintegral 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 timeinterval 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 binwidths 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 timespan.
4.2 Results
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 loweccentricity binaries, which will dominate the burst population in the case, may be less detectable than bursts from higheccentricity binaries. There are proportionally more shortduration 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 GWdriven 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 highmass binaries than the distribution for . This is because, in the case, lowermass 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 highereccentricity binary SMBH populations.

Highereccentricity binary populations result in relatively more shorter duration bursts than lowereccentricity 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 highend powerlaw 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 highcadence PTA observations targeting such bursts are not well motivated.
5 Summary of results
We have used a semianalytic 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 lowfrequency 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, lowredshift 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 turnover 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 turnover frequency is extremely uncertain, and depends on the numbers of massive () binaries and on binary eccentricities. The most massive binary SMBHs predominantly produce the lowestfrequency 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 largereccentricity 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 lowfrequency turnover 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.

Shortduration, 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, modelindependent searches cannot assume any particular spectral shape.

Modeldependent searches and constraints need to carefully account for the uncertainty in model predictions.
6 Conclusions
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 threebody stellar interactions causes the characteristic strain spectrum to be diminished with respect to a powerlaw at Hz. At these low frequencies, the signal is further attenuated if binary SMBHs are typically more eccentric at formation. The lowfrequency signal may be dominated by a few binaries with combined masses () greater than , to a larger extent than predicted in the circular, GWdriven 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 GWdominated regime must involve some form of binaryenvironment coupling. Hence, independent of the exact model, there will always be some lowfrequency attenuation of the GW signal relative to the circular, GWdriven binary case. Our results indicate that this attenuation occurs within the PTA frequency band. However, the strength of the binaryenvironment 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.
Acknowledgements
The authors thank Alberto Sesana, Sarah BurkeSpolaor, 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 MillenniumII 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 threebody scattering problems for binary SMBHs interacting with stars on radial, intersecting orbits drawn from a sphericallysymmetric, 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, powerlaw 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 GWdriven orbital decay dominated the orbital decay caused by threebody stellar interactions. By extrapolating the binary orbits assuming constant stellardriven hardening rates and eccentricities, K12 estimated binary semimajor axes, , below which GWdriven evolution dominated. The accuracy of these extrapolations was confirmed using a selection of simulations including postNewtonian corrections to the binary SMBH orbits.
Here, we take the final eccentricities and semimajor 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.
Model  (fixed )  

A1  0.5  0.1  0.51  0.47 
A2  0.5  0.25  0.72  0.53 
A3  0.5  0.5  0.62  0.57 
A4  0.5  1.0  0.68  0.65 
B1  1.0  0.1  0.58  0.36 
B2  1.0  0.1  0.74  0.38 
B3  1.0  0.1  0.87  0.41 
B4  1.0  0.1  0.93  0.47 
D1  1.75  0.1  0.72  0.41 
D2  1.75  0.1  0.68  0.44 
D3  1.75  0.1  0.62  0.48 
D4  1.75  0.1  0.63  0.51 
We find in all cases. This implies that the S06 model that we use in our work has stronger stellardriven 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 GWdriven 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 semimajor axes at which binary SMBH evolution begins to be GWdominated 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 slowmotion, farfield regime using the quadrupole formula. This tensor, , can be written as (Wahlquist, 1987)
(15) 
where and are the ‘plus’ and ‘cross’ polarisation tensors respectively, and and are the timevarying 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
(16) 
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)
(17) 
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,
(18) 
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
(19) 
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 zaxis. To simplify this, we set , and express as
(20) 
where
(21) 
and
(22) 
The linear independence of the polarisation tensors implies that
(23) 
We find that
(24) 
Then, we have
(25) 
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.
Footnotes
 pagerange: 1–15
 pubyear: 2014
 We refer to this as the “circular, GWdriven case” throughout the paper.
 See Gaddum (1945) for a discussion of the ubiquity of lognormal distributions in nature.
References
 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., BoylanKolchin 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 eprints
 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 Fivehundredmeter 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 eprints
 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 eprints
 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 eprints
 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., BurkeSpolaor 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 eprints
 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., BurkeSpolaor 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., BurkeSpolaor 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