Abstract
Massive black holes are key components of the assembly and evolution of cosmic structures, and a number of surveys are currently on going or planned to probe the demographics of these objects and to gain insight into the relevant physical processes. Pulsar Timing Arrays (PTAs) currently provide the only means to observe gravitational radiation from massive black hole binary systems with masses ≳10^{7} M_{⊙}. The whole cosmic population produces a stochastic background that could be detectable with upcoming PTAs. Sources sufficiently close and/or massive generate gravitational radiation that significantly exceeds the level of the background and could be individually resolved. We consider a wide range of massive black hole binary assembly scenarios, investigate the distribution of the main physical parameters of the sources, such as masses and redshift, and explore the consequences for PTAs observations. Depending on the specific massive black hole population model, we estimate that on average at least one resolvable source produces timing residuals in the range ∼5–50 ns. PTAs, and in particular the future Square Kilometre Array, can plausibly detect these unique systems, although the events are likely to be rare. These observations would naturally complement on the highmass end of the massive black hole distribution function future surveys carried out by the Laser Interferometer Space Antenna.
1 INTRODUCTION
Massive black hole (MBH) binary systems with masses in the range ∼10^{4}–10^{10} M_{⊙} are amongst the primary candidate sources of gravitational waves (GWs) at ∼ nHz–mHz frequencies (see, e.g., Haehnelt 1994; Jaffe & Backer 2003; Wyithe & Loeb 2003; Sesana et al. 2004, 2005). The frequency band ∼10^{−5}–1 Hz will be probed by the Laser Interferometer Space Antenna (Bender et al. 1998), a spaceborne GW laser interferometer being developed by European Space Agency and NASA. The observational window 10^{−9}–10^{−6} Hz is already accessible with Pulsar Timing Arrays (PTAs; e.g. the Parkes radiotelescope; Manchester 2008). PTAs exploit the effect of GWs on the propagation of radio signals from a pulsar to the Earth (e.g. Sazhin 1978; Detweiler 1979; Bertotti, Carr & Rees 1983), producing a characteristic signature in the time of arrival (TOA) of radio pulses. The timing residuals of the fit of the actual TOA of the pulses and the TOA according to a given model carry the physical information about unmodelled effects, including GWs (e.g. Hellings & Downs 1983; Jenet et al. 2005). The complete Parkes PTA (Manchester 2008), the European PTA (Janssen et al. 2008) and NanoGrav1 are expected to improve considerably on the capabilities of these surveys, and the planned Square Kilometre Array (SKA; http://www.skatelescope.org) will produce a major leap in sensitivity.
Popular scenarios of MBH formation and evolution (e.g. Volonteri, Haardt & Madau 2003; Wyithe & Loeb 2003; Koushiappas & Zentner 2006; Malbon et al. 2007; Yoo et al. 2007) predict the existence of a large number of massive black hole binaries (MBHBs) emitting in the frequency range between ∼10^{−9} and 10^{−6} Hz. PTAs can gain direct access to this population, and address a number of unanswered questions in astrophysics (such as the assembly of galaxies and dynamical processes in galactic nuclei), by detecting gravitational radiation of two forms: (i) the stochastic GW background produced by the incoherent superposition of radiation from the whole cosmic population of MBHBs and (ii) GWs from individual sources that are sufficiently bright (and therefore massive and/or close) so that the gravitational signal stands above the rms value of the background. Both classes of signals are of great interest, and the focused effort on PTAs could lead to the discovery of systems difficult to detect with other techniques.
The possible level of the GW background and the consequences for observations have been explored by several authors (see e.g. Rajagopal & Romani 1995; Phinney 2001; Jaffe & Backer 2003; Jenet et al. 2005, 2006; Sesana, Vecchio & Colacino 2008). Recently, Sesana et al. (2008, hereafter Paper I) studied in detail the properties of such a signal and the astrophysical information encoded into it, for a comprehensive range of MBHB formation models. As shown in Paper I, there is over a factor of 10 uncertainty in the characteristic amplitude of the MBHBgenerated background in the PTA frequency window. However, the most optimistic estimates yield an amplitude just a factor of ≈3 below the upper bound placed using current data (Jenet et al. 2006), and nearterm future observations could either detect such a stochastic signal or start ruling out selected MBHB population scenarios. Based on our current astrophysical understanding of the formation and evolution of MBHBs and the estimates of the sensitivity of SKA, one could argue that this instrument guarantees the detection of this signal in the frequency range 3 × 10^{−9}–5 × 10^{−8} Hz for essentially every assembly scenario that is considered at present.
The background generated by the cosmic population of MBHBs is present across the whole observational window of PTAs (cf. Paper I). The Monte Carlo simulations reported in Paper I clearly show the presence of distinctive strong peaks well above the average level of the stochastic contribution (cf. figs 1 and 4 in Paper I). This is to be expected as individual sources can generate gravitational radiation sufficiently strong to stand above the rms value of the stochastic background. These sources are of great interest because they can be individually resolved and likely involve the most massive MBHBs in the Universe. Their observation can offer further insight into the highmass end of the MBH(B) population, galaxy mergers in the lowredshift Universe and dynamical processes that determine the formation of MBH pairs and the evolution to form close binaries with orbital periods of the order of years.
Some exploratory studies have been carried out about detecting individual signals from MBHBs in PTA data (Jenet et al. 2004, 2006). In this paper, we systematically study for a comprehensive range of assembly scenarios the properties, in particular the distribution of masses and redshift, of the sources that give rise to detectable individual events; we compute the induced timing residuals and the expected number of sources at a given timing residual level. To this aim, the modelling of the highmass end of the MBHB population at relatively low redshift is of crucial importance. We generate a statistically significant sample of merging massive galaxies from the online Millennium data base (http://www.gvo.org/Millennium) and populate them with central MBHs according to different prescriptions (Tremaine et al. 2002; McLure et al. 2006; Lauer et al. 2007; Tundo et al. 2007). The Millennium simulation (Springel et al. 2005) covers a comoving volume of (500/h_{100})^{3} Mpc^{3} (h_{100}=H_{0}/100 km s^{−1} Mpc^{−1} is the normalized Hubble parameter), ensuring a number of massive nearby binaries adequate to construct the necessary distribution. For each model, we compute the stochastic background, the expected distribution of bright individual sources and the value of the characteristic timing residual δt_{gw} (see equation 20) for an observation time T. The signaltonoise ratio (SNR) at which a source can then be observed scales as SNR ≈δt_{gw}/δt_{rms}, where δt_{rms} is the rms level of the timing residuals noise, both coming from the receiver and the GW stochastic background contribution. In the following we summarize our main results.

The number of detectable individual sources for different thresholds of the effective induced timing residuals δt_{gw} is shown in Table 1. Depending on the specific MBH population model, we estimate that on average at least one resolvable source produces timing residuals in the range ∼5–50 ns. Future PTAs, and in particular SKA, can plausibly detect these unique systems; the detection is however by no means guaranteed, events will be rare and just above the detection threshold.

As expected, the brightest signals come from very massive systems with . Here, is the chirp mass of the binary and M_{1} > M_{2} are the two black hole masses. Most of the resolvable sources are located at relatively high redshift (0.2 < z < 1.5), and not at z≪ 1 as one would naively expect, giving the opportunity to probe the universe at cosmological distances.

The number of resolvable MBHBs depends on the actual level of the stochastic background generated by the whole population; here, we have used the standard simplified assumption that the background level is determined by having more than one source per frequency resolution bin of width 1/T, where T is the observational time. Using this definition, we find that at frequencies less than 10^{−7} Hz there are typically a few resolvable sources, considering T= 5 yr, with residuals in the range ∼1 nHz−1 μHz. As the level of the background decreases for increasing frequencies, fainter sources become visible individually.

As a sanity check, we have compared the MBHB populations and stochastic background levels obtained using data from the Millennium simulation (adopted in this paper) with those derived by means of merger tree realizations based on the Extended Press & Schechter (EPS) formalism (considered in Paper I) and have found good agreement. This provides an additional validation of the results of this paper and Paper I. Moreover, it supports that EPS merger trees, if handled sensibly, can offer a valuable tool for the study of MBH evolution even at low redshift.
M_{BH}−M_{bulge}Tundo et al. (2007)  M_{BH}−M_{bulge}McLure et al. (2006)  M_{BH}−M_{V}Lauer et al. (2007)  M_{BH}−σ Tremaine et al. (2002)  
Single BH accretion  TuSA  McSA  LaSA  TrSA 
6.9(2.7)  6.6(2.5)  8.1(3.0)  6.2(2.5)  
1.5(1.1)  1.0(0.8)  1.7(1.2)  0.8(0.8)  
0.2(0.4)  0.02(0.1)  0.5(0.6)  0.01(0.1)  
0.04(0.2)  0.002(0.04)  0.1(0.2)  0.002(0.04)  
Double BH accretion  TuDA  McDA  LaDA  TrDA 
8.3(2.9)  7.3(2.7)  9.6(3.2)  7.0(2.8)  
2.2(1.4)  1.6(1.1)  2.6(1.5)  1.2(1.0)  
0.6(0.7)  0.2(0.4)  0.8(0.8)  0.06(0.2)  
0.2(0.4)  0.03(0.2)  0.3(0.5)  0.007(0.1)  
No accretion (before merger)  TuNA  McNA  LaNA  TrNA 
6.4(2.5)  6.0(2.4)  6.8(2.7)  6.0(2.5)  
1.3(1.0)  0.5(0.6)  1.5(1.1)  0.5(0.6)  
0.1(0.3)  0.07(0.1)  0.1(0.3)  0.003(0.05)  
0.02(0.1)  0.001(0.03)  0.02(0.1)  −(−) 
M_{BH}−M_{bulge}Tundo et al. (2007)  M_{BH}−M_{bulge}McLure et al. (2006)  M_{BH}−M_{V}Lauer et al. (2007)  M_{BH}−σ Tremaine et al. (2002)  
Single BH accretion  TuSA  McSA  LaSA  TrSA 
6.9(2.7)  6.6(2.5)  8.1(3.0)  6.2(2.5)  
1.5(1.1)  1.0(0.8)  1.7(1.2)  0.8(0.8)  
0.2(0.4)  0.02(0.1)  0.5(0.6)  0.01(0.1)  
0.04(0.2)  0.002(0.04)  0.1(0.2)  0.002(0.04)  
Double BH accretion  TuDA  McDA  LaDA  TrDA 
8.3(2.9)  7.3(2.7)  9.6(3.2)  7.0(2.8)  
2.2(1.4)  1.6(1.1)  2.6(1.5)  1.2(1.0)  
0.6(0.7)  0.2(0.4)  0.8(0.8)  0.06(0.2)  
0.2(0.4)  0.03(0.2)  0.3(0.5)  0.007(0.1)  
No accretion (before merger)  TuNA  McNA  LaNA  TrNA 
6.4(2.5)  6.0(2.4)  6.8(2.7)  6.0(2.5)  
1.3(1.0)  0.5(0.6)  1.5(1.1)  0.5(0.6)  
0.1(0.3)  0.07(0.1)  0.1(0.3)  0.003(0.05)  
0.02(0.1)  0.001(0.03)  0.02(0.1)  −(−) 
This paper is organized as follows. In Section 2, we describe MBHB population models, in particular the range of scenarios considered in this paper. A short review of the timing residuals produced by GWs generated by an individual binary (in circular orbit) in the data collected by PTAs is provided in Section 3. Section 4 contains the key results of this paper: the expected timing residuals from the estimated population of MBHBs, including detection rates for current and future PTAs. We also provide a comparison between the stochastic background computed according to the prescriptions considered here and the results of Paper I. Summary and conclusions are given in Section 5.
2 THE MASSIVE BLACK HOLE BINARY POPULATION
In this section, we introduce the population synthesis model adopted to estimate the number, and astrophysical parameters of MBHBs that emit GWs in the frequency region probed by PTAs. The two fundamental ingredients to compute the merger rate of MBHBs are (i) the merger history of galaxy haloes and (ii) the MBH population associated with those haloes. We discuss them in turn. In building and evolving the population of sources we follow exactly the same method as in Paper I, to which we refer the reader for further details, with one important difference: the galactic halo merger rates are derived using the data provided by the Millennium simulation, and not EPSbased models. We will justify this choice in the next section, but note that the two methods yield (within the statistical error) the same results. This is a result that is important in itself and has farreaching consequences (outside the specific issues related to PTAs).
2.1 From merger trees to the Millennium data base
In Paper I we used models based on the EPS formalism (Press & Schechter 1974; Lacey & Cole 1993; Sheth & Tormen 1999) that trace the hierarchical assembly of dark matter haloes through a Monte Carlo approach. Although EPSbased models tend to overpredict the bright end of the quasar luminosity function at z < 1 (e.g. Marulli et al. 2006), we showed that EPS halo merger rates at low redshift are consistent with observations of close galaxy pairs (Lin et al. 2004; Bell et al. 2006; De Propris et al. 2007). In this paper we focus on MBHBs, whose GWs induce timing residuals above the stochastic signal from the whole population, and are therefore detectable as individual sources. The population of low/mediumredshift and highmass sources will particularly impact on the results. At low redshift, the EPSbased merger tree outputs need to be handled with care. In models such as those considered in Paper I, each realization of the Universe is obtained by reconstructing the merger history of about 200 dark matter haloes [see Volonteri et al. (2003) for details]. The outputs of the models are a list of coalescences labelled by MBH masses (for a given recipe that associates a MBH mass and with a given dark matter halo) and redshift. These events are then properly weighted over the observable volume shell at each redshift to obtain the distribution (see Paper I) that is the coalescence rate (the number of coalescences N per time interval d t) in the chirp mass and redshift interval and [z, z+ d z], respectively. The resulting distribution is reasonably smooth over most of the plane, but small number statistics become important at z < 0.5 and , which is an important region of the parameter space when one deals with individual sources.
To avoid this problem, in this paper we generate distributions of coalescing MBHBs using the galaxy halo merger rates derived from the online Millennium run data base. The Millennium simulation (Springel et al. 2005) covers a volume of (500/h_{100})^{3} Mpc^{3} and is the ideal tool to construct a statistically representative distribution of massive low/mediumredshift sources. In fact, the typical ensemble of events available to construct the mass function of coalescing binaries is ∼100 times larger than in a typical EPSbased merger tree realization. As a first step, we compile catalogues of galaxy mergers from the semianalytical model of Bertone, De Lucia & Thomas (2007) applied to the Millennium run.
2.2 Populating galaxies with massive black holes
We need to associate with each merging galaxy in our catalogue a central MBH, according to some sensible prescription. The Bertone et al. (2007) catalogue contains many properties of the merging galaxies, including the bulge mass M_{bulge}, and the bulge restframe magnitude M_{V} both of the progenitors and of the merger remnant. This is all we need in order to populate a galaxy with a central MBH. The process is twofold.

We populate the coalescing galaxies with MBHs according to four different MBH–host prescriptions.

M _{BH}−M_{bulge} in the version given by Tundo et al. (2007; ‘Tu’ models, see Table 1):
with a dispersion Δ= 0.33 dex. 
M _{BH}−M_{bulge}, with a redshift dependence in the version given by McLure et al. (2006; ‘Mc’ models, see Table 1):
with a redshiftdependent dispersion Δ= 0.125z+ 0.25 dex (see fig. 3 of McLure et al. 2006). 
M _{BH}−M_{V} as given by Lauer et al. (2007; ‘La’ models, see Table 1):
with dispersion Δ= 0.35 dex. 
M _{BH}−σ as given by Tremaine et al. (2002; ‘Tr’ models, see Table 1):
where σ_{200} is the velocity dispersion in units of 200 km s^{−1} and the assumed dispersion of the relation is Δ= 0.3 dex. σ is obtained by applying the Faber–Jackson (Faber & Jackson 1976) relation in the form reported by Lauer et al. (2007) to the values of M_{V} obtained by the catalogue.To each merging system, we assign MBH masses according to equations (1)–(4) so that we have the masses of the two MBH progenitors, M_{1} and M_{2}. For each prescription, we also calculate the mass of the MBH remnant, M_{r}, using the same equations. In all cases (1)–(4), the remnant mass is M_{r} > M_{1}+M_{2}, consistent with the fact that MBHs are expected to grow predominantly by accretion. We also emphasize that the observed scatter is included in each relation, according to the observational evidence that similar bulges may host significantly different MBHs.


For each MBH–host relation, we consider the following three different accretion scenarios.

The masses of the coalescing MBHs are M_{1} and M_{2}. That is, either no accretion occurs, and the merger remnant, M_{1}+M_{2} < M_{r}, sits below the predicted mass, or accretion is triggered after the MBHB coalescence. We label this accretion mode as ‘NA’ (no accretion; see Table 1).
Postcoalescence accretion is expected for gasrich mergers, where MBH pairing and coalescence is believed to occur on very short timescales (Mayer et al. 2007). If we are to assume that during a galaxy merger the MBH remnant is always brought on the correct correlation with its host, by the combination of merging and accretion, two additional routes are possible.

Accretion is triggered before the MBHB coalescence and only the more massive MBH (M_{1}) accretes mass; in this case, the masses of the coalescing MBHs are αM_{1} and M_{2}, where
We label this accretion mode as ‘SA’ (single BH accretion; see Table 1). 
Accretion is triggered before the MBHB coalescence, and both MBHs are allowed to accrete the same fractional amount of mass; in this case, the masses of the coalescing MBHs are βM_{1} and βM_{2}, where

We label this accretion mode as ‘DA’ (double BH accretion; see Table 1).
The ‘SA’ and the ‘DA’ modes are to be expected in gaspoor mergers, especially in nonequalmass mergers, where the dynamical evolution of the binary is much slower (e.g. Yu 2002) than the infall timescale of the gas (e.g. Cox et al. 2008). In a stellar environment, the orbital decay of MBHBs is expected to be much longer than in a gaseous environment (e.g. Dotti, Colpi & Haardt 2006; Sesana, Haardt & Madau 2007).
The MBHB models that we consider here relay on two assumptions: all bulges host a MBH, and a MBHB always coalesces following the hosts' merger. Regarding the first assumption, dynamical processes, such as gravitational recoil and triple MBH interactions, may deplete bulges from their central MBH. However, if one compares the mass function of coalescing binaries obtained here to the results of EPS merger tree models, where both gravitational recoil and triple interactions are consistently taken into account, one finds that the two distributions (shown in the upperleft panel of Fig. 1) are in excellent agreement, within the statistical uncertainties. This is because triple interactions are likely to eject the lighter MBH from the host, leaving behind a massive binary in the vast majority of the cases (Volonteri et al. 2003; Hoffman & Loeb 2007). Gravitational recoil, on the other hand, may be effective in expelling light MBHs from protogalaxies at high redshift, but has probably a negligible impact on the population of MBHs in the mass range of interest for PTA observations (Volonteri 2007). The assumption that MBHBs coalesce within an Hubble time following the hosts' merger is justified by several recent studies of MBHB dynamics. The stellar distribution interacting with the binary may be efficiently repopulated as a consequence of nonaxisymmetric or triaxial galaxy potentials (Merritt & Poon 2004; Berczik et al. 2006), or by massive perturbers (Perets & Alexander 2008); moreover, if one considers postNewtonian corrections to the binary evolution and the effects of eccentricity, one finds that the coalescence timescale is significantly reduced (Berentzen et al. 2008). If the binary evolution is gasdriven, typical hardening timescales are expected to be shorter than 10^{8} years (Escala et al. 2005; Dotti et al. 2006).
In summary, we build a total of 12 models (four MBH–host prescriptions × three accretion modes). Hereafter, we will refer to each model using the labels associated with the prescriptions employed in it. As an example, the model based on the M_{BH}−σ relation (‘Tr’) with single BH accretion prescription (‘SA’) will be referred to as TrSA (see Table 1 for a summary).
2.3 Computing the coalescence distributions
Assigning a MBH to each galaxy, we obtain a list of coalescences (labelled by MBH masses and redshift); the same output quantity is given by the EPSbased merger trees, used as a starting point in Paper I. Each event in the list can be now properly weighted over the observable volume shell at each redshift to obtain the distribution along the lines described in Paper I.
A further technical detail has to be considered to smooth the coalescence distributions. The Millennium simulation provides better statistics for constructing the mass function of merging objects, but the redshift sampling is rather poor. In fact, the Millennium data base consists of 63 snapshots of the whole simulation. The most recent ones are taken at z= 0, 0.02, 0.04, 0.064 and 0.089; we need, at least at low redshift, to spread the events over a continuum in z, to obtain a sensible distance distribution of the GW sources. To this aim, we decouple the dependence in the differential mass function at z < 0.3 and rewrite in the following form:
By doing so, we are redistributing according to a given function Ψ(z) the average mass function obtained at z < 0.3. This is justified since, as expected, the mass function does not show any significant evolution for redshifts below 0.3. The dependence on z should be of the form Ψ(z) ∝n^{2}_{G}× dV_{C}/dz, where n_{G} is the galaxy/MBH number density and dV_{C}/dz is the differential comoving volume shell. At such small redshifts, the impact of merger activity on galaxy/MBH number density is negligible (of the order of 0.1 Gyr^{−1}; e.g. Bell et al. 2006; Masjedi et al. 2006; White et al. 2007; Wake et al. 2008); therefore, we assume n_{G} to be constant. On the other hand, the Universe can be considered Euclidean, so that the differential volume shell is just proportional to z^{2}. We then obtain the coalescing MBHB distribution in the form where C is a normalization factor set by the conditionThe distributions of coalescing binaries are shown in Fig. 1 for all the MBH–host prescriptions assuming accretion on M_{1} only (models TrSA, TuSA, McSA, LaSA). The topleft panel also shows the distribution obtained by a reference EPS merger tree model (VHMhopk; Volonteri, Salvaterra & Haardt 2006) used in Paper I. The agreement with the M_{BH}−σ prescription (‘Tr’) is good for , although the statistical errors are large due to low number statistics. The discrepancy for is due to both the resolution limit of the Millennium simulation and the fact that we relate MBHs to bulges. However, the fact that our sample may be incomplete for has little (if any) impact on the results of this paper: the MBHB population observable with PTAs is by far dominated by sources with , as we have shown in Paper I, and even more so when we consider systems that can be individually resolved. We will further discuss this point in Section 4. Note that, as discussed by, e.g., Lauer et al. (2007) and Tundo et al. (2007), the highmass end of the population derived using the M_{BH}−σ relation (‘Tr’ models) drops very steeply. The drop is much faster than in all the other cases that is for distributions inferred from the M_{BH}−M_{bulge} or the M_{BH}−M_{V} relations (‘Tu’, ‘Mc’, ‘La’ models). This is because σ seems to converge to a plateau in the limit of very massive galaxies (Lauer et al. 2007).
3 TIMING RESIDUALS FROM RESOLVED MASSIVE BLACK HOLE BINARIES
The search for GWs using timing data exploits the effect of gravitational radiation on the propagation of the radio waves from one (or more) pulsar(s). The characteristic signature of GWs on the TOA of radio pulses (e.g. Sazhin 1978; Detweiler 1979; Bertotti et al. 1983) is a linear combination of the two independent GW polarizations. In practice, the analysis consists in computing the difference between the expected and actual TOA of pulses; the timing residuals contain information on all the effects that have not been included in the fit, including GWs. In this section, we summarize the observed signal produced by GWs, following closely Jenet et al. (2004), to which we refer the reader for further details.
The observed timing residual generated by a GW source described by the independent polarization amplitudes h_{+,×} is
where t is the time at the receiver, μ is the opening angle between the GW source and the pulsar relative to the Earth and ψ is the GW polarization angle. The two functions, r_{+,×}(t), are defined as Note that r^{(e)}_{+,×}(t) and r^{(p)}_{+,×} (t) have the same functional form, and result from the integration of the time evolution of the polarization amplitudes at different times, with a delay ≃3.3 × 10^{3}(d/1 kpc) (1 − cos μ) yr, where d is the distance of the pulsar from the Earth. For GWs propagating exactly along the Earth–pulsar direction, there is no effect on the TOAs [r(t) = 0 for cos μ=±1].From now on we will concentrate on the timing residuals produced by binary systems in circular orbit. We model gravitational radiation at the leading quadrupole Newtonian order that is fully justified by the fact that binaries in the mass and frequency range considered here are far from the final merger; in fact, the time to coalescence is . The timing residuals (10) can be written as (Jenet et al. 2004)
where In the previous expressions, ι is the source inclination angle, Φ(t) is the GW phase related to the frequency f(t) (twice the orbital frequency) by and is an overall scalefactor that sets the size of the residuals and D is the luminosity distance to the GW source. The expression for r^{(p)}(t) can be simply obtained from equations (15)–(17) by shifting the time t→ t − d (1 − cos μ)/c (see equations 12 and 13). The two functions, a_{+} and a_{×}, are the ‘antenna beam patterns’ that depend on the source location in the sky and the polarization of the wave.MBHBs observable with PTAs produce a quasimonochromatic signal – the frequency change is – of known form (though unknown parameters). The optimal data analysis approach to search for these signals is by means of the wellknown technique of matchedfiltering. The data set can be represented as
where r(t) is the contribution from GWs and δt_{n}(t) accounts for the fluctuations due to noise; the latter contribution is the superposition of the intrinsic noise in the measurements and the GW stochastic background from the whole population of MBHBs. The angleaveraged optimal SNR at which a signal from a MBHB radiating at (GW) frequency ≈f can be detected using a single pulsar is In the previous expression, δt_{rms}(f) is the rms value of the noise level δt_{n} at frequency f; 〈. 〉 represents the average over the source position in the sky and orientation of the orbital plane and δt_{gw}(f) is the characteristic amplitude of the timing residual over the observation time T defined as where the numerical prefactor comes from the angle average of the amplitude of the signalEquation (19) is appropriate to describe observations using a single pulsar. In reality one can take advantage of the several pulsars that are continuously monitored to increase the SNR, and therefore the confidence of detection: coherently adding the residuals from several pulsars – currently the Parkes PTA contains about 20 pulsars, and more are expected to be available with future instruments – yields an increase in SNR proportional to the square root of the number of pulsars used in the observations. We will use the characteristic amplitude of the residuals δt_{gw} to quantify the strength of a GW signal in PTA observations; δt_{gw} can be used to compute the SNR in a straightforward way, as a function of the noise level and number of pulsars in the array (all of which are quantities that do not depend on the astrophysical model), and therefore assesses the probability of detection of sources in the context of a given MBHB assembly scenario.
4 RESULTS
For each of the 12 models considered in Section 2– the models are the result of four MBH–host galaxy prescriptions and three different accretion scenarios – that encompass a very broad range of MBHB's assembly scenarios, we compute the number of sources that are potentially resolvable individually and several statistical properties of the population, such as the redshift, chipmass and frequency distributions, by means of Monte Carlo realizations of the whole population of MBHBs according to a given model. Before describing the results, we provide details about each step in the computation of the relevant quantities.
The distribution given by equation (8) is straightforwardly converted into , i.e. the comoving number of binaries emitting in a given logarithmic (restframe) frequency interval with chirp mass and redshift in the range and [z, z+ dz], along the lines described in section 3 of Paper I. As in Paper I, we then assume that in the frequency range of interest (f > 3 × 10^{−9} Hz) the binary evolution is driven by GW emission only. This is a reasonable assumption, since the coalescence timescale for these systems is typically ≲10^{6} years; any other putative mechanism (i.e. star ejection, gas torques) of angular momentum removal must have an enormous efficiency to compete with radiation reaction on such short timescales. For each of the 12 models, we estimate the GW stochastic background (and the corresponding rms value of the timing residuals as a function of frequency) generated by the sources following the scheme described in section 4 of Paper I. Finally we generate distribution of bright, individually resolvable sources by running 1000 Monte Carlo realizations of the whole population of MBHBs and by selecting only those sources whose characteristic timing residuals, equation (20), exceed the stochastic background level. We note that the result of which and how many sources raise above the average noise level depends on the duration of the observation T for two reasons: (i) T affects the size of the observational window in frequency space, in particular the minimum frequency 1/T that can be reached, and (ii) the background level decreases as the observation time increases (as the size of the frequency resolution bin δf= 1/T decreases), enhancing the number of individually resolvable sources. For the results that are described here and summarized in Table 1, we set T= 5 yr; increasing the data span to 10 years, the background level would be slightly lower and the additional resolved sources would be barely brighter than the background. The statistics of bright, wellresolvable sources are basically unaffected. We can then cast the results in terms of the cumulative number of resolvable sources as a function of the timing residuals, according to
where the integral is restricted to the sources that produce residuals above the rms level of the stochastic background; if we do not consider this additional constraint, then, for any given δt_{gw}, N(δt_{gw}) simply returns the total number of sources in the Monte Carlo realization above that particular residual threshold.Each Monte Carlo realization clearly yields a different value for N(δt_{gw}) (or its distribution according to a given parameter); the values quoted in the next section and summarized in Table 1 refer to the sample mean computed over the set of Monte Carlo realizations and the sample standard deviation. Fig. 4 quantifies the typical 1σ error in our estimate of the number of sources.
4.1 Single resolvable binaries
The large number of Monte Carlo realizations allows us to study the details of the properties of the individual sources in a statistical sense. We concentrate in particular on the physical properties of the population, such as the expected number of sources per logarithmic frequency interval dN(δt_{gw})/dlog f, chirp mass range and redshift d N(δt_{gw})/d z, and the observable parameters, such as the timing residuals produced by each system and the overall expected number of resolvable MBHBs at a given level of timing residual noise. A summary of the typical range of information that can be extracted from the simulations is shown in Fig. 2 for the specific model TuSA. The topleft panel shows the induced residuals of each single source compared to the level produced by the stochastic background from the whole population; the plot clearly shows the importance of taking into account the additional ‘noise contribution’ from the brightness of the GW sky in considering the detectability of resolvable systems. There are many sources inducing residuals above, say the 5 ns level, however most of them contribute to the buildup of the background and are not individually resolvable. The expected number of bright resolvable MBHBs at frequencies <10^{−7} Hz, or at a timing level >1 ns, is typically around 10, with fainter sources resolvable at higher frequencies. The topright panel shows the mean number of individual sources detectable as a function of δt_{gw} from 1000 Monte Carlo realizations of the emitting population. In this particular case, a sensitivity of ≈10 ns is required to resolve an individual source; for a timing precision of 100 ns, there is a 5 per cent chance to observe a particularly bright source. Note that at the 1 ns level, there are ∼100 MBHBs contributing to the signal, however 90 per cent of them contribute to the background and only about 10 sources are individually resolvable. In the two bottom panels of Fig. 2, we plot the frequency and chirp mass distributions of resolvable sources for selected values of the minimum detectable residual amplitude δt_{gw}. Not surprisingly, the chirp mass of observable systems decreases for smaller values of the considered residual threshold, however even for a rms level of 1 ns all the systems are characterized by . The frequency distribution shows instead a peak corresponding to the frequency at which the background level equals the selected value of δt_{gw}. At higher frequencies, the number of sources drastically drops because the number of emitting binaries is a quite steep function of frequency, N(f) ∝f^{−8/3}; the decrease at lower frequencies is because most of the emitters actually contribute to the background and are not individually resolved (as clearly shown by the dotted lines). The qualitative behaviour of the results obtained using different astrophysical models is very similar to the one described in Fig. 2, with differences that affect only the numerical values of the different quantities.
The central question that we want to address in this paper is what is the expected number of individually resolvable sources that produce an effective timing residual above a given value, as a function of different models of MBHB formation and evolution. We summarize these results in Figs 3 and 4; and in Table 1, where we show the mean total number of individual sources that exceed a given level of timing residual (as defined by equation 22), as a function of the timing residual. The qualitative behaviour of the results is similar for all the scenarios, but the actual numbers vary significantly. In fact, both MBH–host relations and accretion prescriptions have a strong impact on the statistics of the bright sources and, consequently, on their detection. We analyse the effect of the black hole populations and of the accretion in turn. The lefthand panel of Fig. 3 shows results where the accretion prescription is the same (‘SA’), but the underlying MBH–host relation changes; on the other hand, in the righthand panel, we select two MBH–host relations, but change the accretion history. The lefthand panel shows that assuming a sensitivity threshold of 30 ns, one expects to observe of the order of 1 source in the LaSA model, while there is only a probability of ≈5 per cent for the TrSA model. If we maintain the same MBH–host relation and consider different accretion scenarios for, e.g., the Lauer et al. population, the mean number of expected sources varies by a factor of ≈5 between 0.3 (LaNA) and 1.5 (LaDA) (see righthand panel). Fig. 3 shows that in the most pessimistic case – TrNA model that has the sharpest cut off of the MBHB mass function bright end – a precision of ≈5 ns should guarantee a positive detection; in the optimistic LaDA case, a precision of ≈50 ns should be sufficient. In turn, the timing precision required for positive detection is in the range 5−50 ns that is basically consistent with a factor of ≈10 uncertainty in the background level estimated in Paper I.
The typical spread around the mean values obtained in the Monte Calro realizations is shown in Fig. 4 for selected models. When N(δt_{gw}) ≫ 1, the 1σ range is roughly the Poisson error around the mean (reflecting the uncorrelated nature of sources in each Monte Carlo realization). When N(δt_{gw}) < 1, it can be interpreted as the probability to find a single source above the considered δt_{gw} value if the actual MBHB population in the Universe follows the prescription given by the considered model; in this case, a nondetection is trivially consistent with the model predictions. Table 1 provides a summary of the results for the 12 models.
It is also interesting to investigate the massredshift distribution of the detectable sources. Fig. 5 shows the expected number of detectable individual sources (as a function of the residuals amplitude) for different redshift and mass ranges. Obviously, higher timing residuals correspond to higher , since the strength of the signal is proportional to , and the most likely sources to be detected fall in the range (MBHBs with produce indeed larger timing residuals, but are also much rarer). Interestingly, the vast majority of detectable sources are at redshift 0.1 ≲z≲ 1, which shows that PTAs could probe the mediumredshift Universe, and are unlikely to discover nearby sources. The reason is simply that, at least at small redshift, the Universe volume increases as z^{3}.
A summary of the properties of individual resolvable sources, and dN(δt_{gw})/d log z, is given in Figs 6, 7 and 8, respectively, for all the four MBHB population models considered here, with accretion limited to a single black hole prior to merger considering different residual thresholds δt_{gw}= 1, 10, 50, 100 ns. All the models show the same qualitative features, as we have highlighted before. The frequency distribution shown in Fig. 6 was discussed above and is the same for all the models. The distribution of the detectable sources as a function of chirp mass (Fig. 7) peaks at ≈ 3 × 10^{8} M_{⊙} for all models assuming δt_{gw}= 1 ns. Increasing δt_{gw}, the distribution peak shifts towards higher and the mean number of events is strongly modeldependent for δt_{gw} > 10 ns (cf. the values in Table 1). The redshift distributions (Fig. 8) consistently show a broad peak in the range 0.2 < z < 1, due to the volume effect previously discussed. Note that in the ‘Mc’ model the peak shifts towards higher redshifts as δt_{gw} increases, because in this model similar galaxies are populated by more MBHs if found at higher redshifts (see equation 2).
4.2 Stochastic background
As a sanity check, we compare the stochastic backgrounds derived according to the MBH populations inferred using the Millennium simulation to the predictions of EPSbased models reported in Paper I (the reader is referred to Paper I for the technical details). In Fig. 9, we show Monte Carlo generated signals for each model; in each panel, we plot the stochastic levels according to the three different accretion modes discussed in Section 2. Both the accretion prescription and the adopted MBH–host correlation influence the level of the background. If accretion occurs on to the remnant (i.e. after coalescence, ‘NA’ models), the predicted characteristic amplitude of the GW background can be up to a factor of 3 lower with respect to models in which both the MBHs accrete before the final coalescence (‘DA’ models); on the other hand, M_{BH}−σ (‘Tr’) models predict lower backgrounds compared to M_{BH}−M_{bulge} (‘Tu’, ‘Mc’) and M_{BH}−M_{V} (‘La’) models. A comparison of these results with those presented in Paper I is given in Fig. 10. At 10^{−8} Hz, the models studied here cover a characteristic amplitude range consistent with the uncertainty estimated in Paper I. The TrSA model predicts a stochastic background that agrees with the typical EPS model within 30 per cent for f < 10^{−7} Hz. The highfrequency end is instead steeper; this effect is caused by incompleteness in the lowmass end of the MBH population. As shown in Fig. 1, the mass function of coalescing MBHB derived from the Millennium simulation is not consistent with the one obtained using the EPS formalism for . The weight of a single dark matter particle in the Millennium simulation is 8.66/h_{100}× 10^{8} M_{⊙}, allowing the reconstruction of haloes with minimum mass of the order of ≈5 × 10^{10} M_{⊙}. Assuming a barionic fraction of 0.1, the simulation is then incomplete for barionic structures less massive than ≈5 × 10^{9} M_{⊙}. We checked this by plotting the mass function of barionic structures and finding a sudden drop below 10^{9} M_{⊙}. It is then inevitable that in the results derived from the Millennium simulation most of the MBHs with mass below a few ×10^{6} M_{⊙} are missing. Since many of these MBHs are expected to merge with more massive ones during cosmic history, the (spurious) lack of MBHs in this mass rage explains the flattening of the mass function shown in the topleft panel of Fig. 1. All the backgrounds are rather similar at f > 10^{−7} Hz because all the MBH prescriptions adopted lead to similar MBH mass functions at M_{BH} < 10^{8} M_{⊙} (this fact is independent of the incompleteness issue). This means that the slope of the background on the right of the knee has a welldefined dependence upon the adopted MBH population: the more pronounced is the highmass tail of the MBH mass function, the steeper is the highfrequency end of the GW background. In turn, models constructed using the Millennium simulation confirm the findings of Paper I.
5 SUMMARY AND CONCLUSIONS
We have investigated the ability of PTAs to resolve individual MBH binary systems by detecting gravitational radiation produced during the inspiral phase through its effect on the residuals of the TOA of signals from radio pulsars. We have considered a broad range of assembly scenarios, using the data of the Millennium simulation to evaluate the galactic haloes merger rates, and a total of 12 different models that control the relations between the mass of the central black holes and the galactic haloes, and the evolution of the black hole masses through accretion. These models therefore cover qualitatively (and to large extent quantitatively) the whole spectrum of MBH assembly scenarios currently considered. Regardless of the model, we estimate that at least one resolvable source is expected to produce timing residuals in the range ∼5–50 ns, and therefore future PTAs, and in particular the SKA may be able to observe these systems. A whistlestop summary of the models and results is contained in Table 1. The total number of visible events clearly depends on the sensitivity of PTAs and on the astrophysical scenario. As expected, the brightest sources (for PTAs) are very massive binaries with chirp mass . However (initially) quite surprisingly most of the resolvable sources are located at relatively high redshift (z > 0.2). In conjunction with the observation of the stochastic GW background from the whole population of MBHBs, the identification of individual MBHBs could provide new constraints on the populations of these objects and the relevant physical processes.
As a byproduct of the analysis, we have also estimated the level of the GW stochastic background produced by the different models, finding good agreement with the estimates derived using merger tree realizations based on the EPS formalism considered in Paper I. Such agreement provides a further validation of the results of this paper and Paper I, it shows that we can now have in hand selfconsistent predictions for stochastic and deterministic signals from the cosmic population of MBHBs, and suggests that EPS merger trees could provide a valuable approach to the studies of MBH evolution at lowtomedium redshift.
As a final word of caution, we would like to stress that the results of this paper clearly suffer from considerable uncertainties determined by the still poor quantitative information about several parameters that control the models. The spread of the predictions of the expected events is therefore likely dominated by the lack of knowledge of the model parameters, rather than the differences between the assembly scenarios.
AS gratefully acknowledges the support of National Science Foundation Grant no. PHY 0653462 and no. PHY 0555615, and NASA Grant no. NNG 05GF71G, awarded to the Pennsylvania State University.