Gravitational waves from resolvable massive black hole binary systems and observations with Pulsar Timing Arrays
aa r X i v : . [ a s t r o - ph ] M a r Mon. Not. R. Astron. Soc. , 000–000 (0000) Printed 22 October 2018 (MN L A TEX style file v2.2)
Gravitational waves from resolvable massive black holebinary systems and observations with Pulsar TimingArrays
A. Sesana , A. Vecchio and M. Volonteri Center for Gravitational Wave Physics, The Pennsylvania State University, University Park, PA 16802, USA School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham, B15 2TT, UK Dept. of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA
Received —
ABSTRACT
Massive black holes are key components of the assembly and evolution of cosmicstructures and a number of surveys are currently on-going or planned to probe thedemographics of these objects and to gain insight into the relevant physical processes.Pulsar Timing Arrays (PTAs) currently provide the only means to observe gravita-tional radiation from massive black hole binary systems with masses > ∼ M ⊙ . Thewhole cosmic population produces a stochastic background that could be detectablewith upcoming Pulsar Timing Arrays. Sources sufficiently close and/or massive gener-ate gravitational radiation that significantly exceeds the level of the background andcould be individually resolved. We consider a wide range of massive black hole binaryassembly scenarios, we investigate the distribution of the main physical parametersof the sources, such as masses and redshift, and explore the consequences for PulsarTiming Arrays observations. Depending on the specific massive black hole populationmodel, we estimate that on average at least one resolvable source produces timingresiduals in the range ∼ −
50 ns. Pulsar Timing Arrays, and in particular the futureSquare Kilometre Array (SKA), can plausibly detect these unique systems, althoughthe events are likely to be rare. These observations would naturally complement on thehigh-mass end of the massive black hole distribution function future surveys carriedout by the Laser Interferometer Space Antenna (
LISA ). Key words: black hole physics, gravitational waves – cosmology: theory – pulsars:general
Massive black hole (MBH) binary systems with masses inthe range ∼ − M ⊙ are amongst the primary candi-date sources of gravitational waves (GWs) at ∼ nHz - mHzfrequencies (see, e.g., Haehnelt 1994; Jaffe & Backer 2003;Wyithe & Loeb 2003, Sesana et al. 2004, Sesana et al. 2005).The frequency band ∼ − Hz − Laser Interferometer Space Antenna ( LISA , Bender et al.1998), a space-borne gravitational wave laser interferome-ter being developed by ESA and NASA. The observationalwindow 10 − Hz − − Hz is already accessible with Pul-sar Timing Arrays (PTAs; e.g. the Parkes radio-telescope,Manchester 2008). PTAs exploit the effect of GWs on thepropagation of radio signals from a pulsar to the Earth (e.g.Sazhin 1978, Detweiler 1979, Bertotti et al. 1983), produc-ing a characteristic signature in the time of arrival (TOA)of radio pulses. The timing residuals of the fit of the actualTOA of the pulses and the TOA according to a given model, carry the physical information about unmodelled effects, in-cluding GWs (e.g. Helling & Downs 1983, Jenet et al. 2005).The complete Parkes PTA (Manchester 2008), the EuropeanPulsar Timing Array (Janssen et al. 2008), and NanoGrav are expected to improve considerably on the capabilitiesof these surveys and the planned Square Kilometer Array(SKA; ) will produce a major leap insensitivity.Popular scenarios of MBH formation and evolution (e.g.Volonteri, Haardt & Madau 2003; Wyithe & Loeb 2003,Koushiappas & Zentner 2006, Malbon et al. 2007, Yoo etal. 2007) predict the existence of a large number of mas-sive black hole binaries (MBHB) emitting in the frequencyrange between ∼ − Hz and 10 − Hz. PTAs can gaindirect access to this population, and address a number of http://arecibo.cac.cornell.edu/arecibo-staging/nanograv/c (cid:13) A. Sesana et al. unanswered questions in astrophysics (such as the assem-bly of galaxies and dynamical processes in galactic nuclei),by detecting gravitational radiation of two forms: (i) thestochastic GW background produced by the incoherent su-perposition of radiation from the whole cosmic population ofMBHBs and (ii) GWs from individual sources that are suf-ficiently bright (and therefore massive and/or close) so thatthe gravitational signal stands above the root-mean-square(rms) value of the background. Both classes of signals areof great interest, and the focused effort on PTAs could leadto the discovery of systems difficult to detect with othertechniques.The possible level of the GW background, and the con-sequences for observations have been explored by severalauthors (see e.g.
Rajagopal & Romani 1995; Phinney 2001,Jaffe & Backer 2003; Jenet et al. 2005; Jenet et al. 2006;Sesana et al. 2008). Recently, Sesana Vecchio & Colacino(2008, hereinafter PaperI) studied in details the propertiesof such a signal and the astrophysical information encodedinto it, for a comprehensive range of MBHB formation mod-els. As shown in PaperI, there is over a factor of 10 uncer-tainty in the characteristic amplitude of the MBHB gener-ated background in the PTA frequency window. However,the most optimistic estimates yield an amplitude just a fac-tor ≈ × − Hz − × − Hzfor essentially every assembly scenario that is considered atpresent.The background generated by the cosmic population ofMBHBs is present across the whole observational window ofPTAs (cf. PaperI). The Monte Carlo simulations reported inPaperI show clearly the presence of distinctive strong peakswell above the average level of the stochastic contribution(cf. Figure 1 and 4 in PaperI). This is to be expected, asindividual sources can generate gravitational radiation suf-ficiently strong to stand above the rms value of the stochas-tic background. These sources are of great interest becausethey can be individually resolved and likely involve the mostmassive MBHBs in the Universe. Their observation can of-fer further insight into the high-mass end of the MBH(B)population, galaxy mergers in the low-redshift Universe anddynamical processes that determine the formation of MBHpairs and the evolution to form close binaries with orbitalperiods of the order of years.Some exploratory studies have been carried out aboutdetecting individual signals from MBHBs in PTA data(Jenet et al. 2004, 2005). In this paper we study system-atically for a comprehensive range of assembly scenarios theproperties, in particular the distribution of masses and red-shift, of the sources that give rise to detectable individualevents; we compute the induced timing residuals and the ex-pected number of sources at a given timing residual level. Tothis aim, the modelling of the high-mass end of the MBHBpopulation at relatively low redshift is of crucial impor-tance. We generate a statistically significant sample of merg-ing massive galaxies from the on-line Millennium database ( ) and populate them withcentral MBHs according to different prescriptions (Tremaineet al. 2002, Mclure et al. 2006, Lauer et al. 2007, Tundo et al.2007). The Millennium simulation (Springel et al. 2005) cov-ers a comoving volume of (500 /h ) Mpc ( h = H / − Mpc − is the normalized Hubble parameter), ensur-ing a number of massive nearby binaries adequate to con-struct the necessary distribution. For each model we com-pute the stochastic background, the expected distribution ofbright individual sources and the value of the characteristictiming residual δt gw , see Equation (20), for an observationtime T . The signal-to-noise ratio at which a source can thenbe observed scales as SNR ≈ δt gw /δt rms where δt rms is theroot-mean-square level of the timing residuals noise, bothcoming from the receiver and the GW stochastic backgroundcontribution. In the following we summarise our main re-sults:(i) The number of detectable individual sources for dif-ferent thresholds of the effective induced timing residuals δt gw is shown in Table 1. Depending on the specific MBHpopulation model, we estimate that on average at least oneresolvable source produces timing residuals in the range ∼ − M > × M ⊙ . Here M = M / M / / ( M + M ) / is the chirp mass of the binaryand M > M are the two black hole masses. Most ofthe resolvable sources are located at relatively high redshift(0 . < z < . z ≪ /T , where T is the observational time. Using thisdefinition we find that at frequencies less than 10 − Hz thereare typically a few resolvable sources, considering T = 5 yrs,with residuals in the range ∼ − µ Hz. As the level ofthe background decreases for increasing frequencies, faintersources become visible individually.(iv) As a sanity check, we have compared the MBHB pop-ulations and stochastic background levels obtained usingdata from the Millennium simulation (adopted in this pa-per) with those derived by means of merger tree realisationsbased on the Extended Press & Schechter (EPS) formalism(considered in PaperI) and have found good agreement. Thisprovides an additional validation of the results of this paperand PaperI. Moreover it supports that EPS merger trees, ifhandled sensibly, can offer a valuable tool for the study ofMBH evolution even at low redshift.The paper is organised as follows. In Section 2 we de-scribe MBHB population models, in particular the range ofscenarios considered in this paper. A short review of the tim-ing residuals produced by GWs generated by an individualbinary (in circular orbit) in the data collected by PTAs isprovided in Section 3. Section 4 contains the key results of c (cid:13) , 000–000 assive black holes and pulsar timing the paper: the expected timing residuals from the estimatedpopulation of MBHBs, including detection rates for currentand future PTAs. We also provide a comparison between thestochastic background computed according to the prescrip-tions considered here and the results of PaperI. Summaryand conclusions are given in Section 5. In this section we introduce the population synthesis modeladopted to estimate the number, and astrophysical param-eters of MBHBs that emit GWs in the frequency regionprobed by PTAs. The two fundamental ingredients to com-pute the merger rate of MBHBs are (i) the merger historyof galaxy haloes, and (ii) the MBH population associatedto those haloes. We discuss them in turn. In building andevolving the population of sources we follow exactly thesame method as in PaperI, to which we refer the reader forfurther details, with one important difference: the galactichaloes merger rates are derived using the data provided bythe Millennium simulation, and not EPS-based models. Wewill justify this choice in the next sub-Section, but note thatthe two methods yield (within the statistical error) the sameresults. This is a result that is important in itself and hasfar reaching consequences (outside the specific issues relatedto PTAs).
In PaperI we used models based on the Extended Press &Schechter (EPS) formalism (Press & Schechter 1974, Lacey& Cole 1993, Sheth & Tormen 1999), that trace the hier-archical assembly of dark matter haloes through a Monte-Carlo approach. Although EPS–based models tend to over-predict the bright end of the quasar luminosity function at z < d N/d M dzdt (seePaperI), that is the coalescence rate (the number of coales-cences N per time interval dt ) in the chirp mass and redshiftinterval [ M , M + d M ] and [ z, z + dz ], respectively. The re-sulting distribution is reasonably smooth over most of the( M , z ) plane, but small number statistics becomes impor-tant at z < . M > M ⊙ , which is an important region of the parameter space when one deals with individ-ual sources.To avoid this problem, in this paper we generate distri-butions of coalescing MBHBs using the galaxy haloes mergerrates derived from the on-line Millennium run database.The Millennium simulation (Springel et al. 2005) covers avolume of (500 /h ) Mpc and is the ideal tool to con-struct a statistically representative distribution of massivelow/medium–redshift sources. In fact, the typical ensembleof events available to construct the mass function of coa-lescing binaries is ∼
100 times larger than in a typical EPS-based merger tree realization. As a first step we compile cat-alogues of galaxy mergers from the semi-analytical model ofBertone et al. (2007) applied to the Millennium run.
We need to associate to each merging galaxy in our cataloguea central MBH, according to some sensible prescription. TheBertone et al. 2007 catalogue contains many properties ofthe merging galaxies, including the bulge mass M bulge , andthe bulge rest frame magnitude M V both of the progenitorsand of the merger remnant. This is all we need in orderto populate a galaxy with a central MBH. The process istwofold.(i) We populate the coalescing galaxies with MBHs ac-cording to four different MBH-host prescriptions: • M BH − M bulge in the version given by Tundo et al.(2007, “Tu” models, see Table 1): M BH M ⊙ = 8 .
31 + 1 .
12 log (cid:18) M bulge M ⊙ (cid:19) , (1)with a dispersion ∆ = 0 .
33 dex. • M BH − M bulge , with a redshift dependence in the ver-sion given by Mclure et al. (2006, “Mc” models, see Table1): M BH M bulge = 2 .
07 log(1 + z ) − . , (2)with a redshift dependent dispersion ∆ = 0 . z + 0 . • M BH − M V as given by Lauer et al. (2007, “La” mod-els, see Table 1): M BH M ⊙ = 8 . − . (cid:16) M V + 222 . (cid:17) , (3)with dispersion ∆ = 0 .
35 dex. • M BH − σ as given by Tremaine et al. (2002, “Tr”models, see Table 1): M BH M ⊙ = 8 .
13 + 4 .
02 log σ , (4)where σ is the velocity dispersion in units of 200 kms − , and the assumed dispersion of the relation is ∆ = 0 . σ is obtained applying the Faber-Jackson (Faber &Jackson 1976) relation in the form reported by Lauer etal. 2007 to the values of M V obtained by the catalogue.To each merging system we assign MBH masses accordingto equations (1)-(4) so that we have the masses of the twoMBH progenitors, M and M . For each prescription we also c (cid:13) , 000–000 A. Sesana et al. M BH − M bulge M BH − M bulge M BH − M V M BH − σ Tundo et (2007) Mclure et al. (2006) Lauer et al. (2007) Tremaine et al. (2002)Single BH accretion Tu-SA Mc-SA La-SA Tr-SA6 . .
7) 6 . .
5) 8 . .
0) 6 . . . .
1) 1 . .
8) 1 . .
2) 0 . . . .
4) 0 .
02 (0 .
1) 0 . .
6) 0 .
01 (0 . .
04 (0 .
2) 0 .
002 (0 .
04) 0 . .
2) 0 .
002 (0 . . .
9) 7 . .
7) 9 . .
2) 7 . . . .
4) 1 . .
1) 2 . .
5) 1 . . . .
7) 0 . .
4) 0 . .
8) 0 .
06 (0 . . .
4) 0 .
03 (0 .
2) 0 . .
5) 0 .
007 (0 . . .
5) 6 . .
4) 6 . .
7) 6 . . . .
0) 0 . .
6) 1 . .
1) 0 . . . .
3) 0 .
07 (0 .
1) 0 . .
3) 0 .
003 (0 . .
02 (0 .
1) 0 .
001 (0 .
03) 0 .
02 (0 . − ( − ) Table 1.
The table summarises the 12 models of assembly of massive black hole binary populations considered in the paper – ”Tu”,”Mc”, ”La” and ”Tr” identify the MBH-host relation; “SA”, “DA”, and “NA” label the accretion mode; full details on the models aregiven in Section 2 – and the total number of individually resolvable systems N ( δt gw ) for selected values of the characteristic timingresiduals (for each model, from top to bottom δt gw = 1 , ,
50 and 100 ns considering an integration time of 5 yrs). The values in thetable are the sample mean and standard deviation within brackets computed over the 1000 Monte-Carlo realisations for each model. calculate the mass of the MBH remnant, M r , using the sameequations. In all cases (1)-(4), the remnant mass is M r >M + M , consistent with the fact that MBHs are expectedto grow predominantly by accretion. We also emphasize thatthe observed scatter is included in each relation, accordingto the observational evidence that similar bulges may hostsignificantly different MBHs.(ii) For each MBH-host relation we consider three differ-ent accretion scenarios: • The masses of the coalescing MBHs are M and M .That is, either no accretion occurs, and the merger rem-nant, M + M < M r , sits below the predicted mass, oraccretion is triggered after the MBHB coalescence. We la-bel this accretion mode as “NA” (no accretion, see Table1).Post-coalescence accretion is expected for gas-richmergers, where MBH pairing and coalescence is believedto occur on very short timescales (Mayer et al. 2007). Ifwe are to assume that during a galaxy merger the MBHremnant is always brought on the correct correlation withits host, by the combination of merging and accretion, twoadditional routes are possible. • Accretion is triggered before the MBHB coalescenceand only the more massive MBH ( M ) accretes mass; inthis case the masses of the coalescing MBHs are αM and M , where α = M r − M M − . (5)We label this accretion mode as “SA” (single BH accre-tion, see Table 1). • Accretion is triggered before the MBHB coalescenceand both MBHs are allowed to accrete the same fractionalamount of mass; in this case the masses of the coalescingMBHs are βM and βM , where β = M r M + M − . (6)We label this accretion mode as “DA” (double BH accre-tion, see Table 1).The ”SA” and the ”DA” modes are to be expected in gas-poor mergers, especially in non-equal-mass mergers, wherethe 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 MBHBsis expected to be much longer than in a gaseous environment(e.g., Sesana et al. 2007, Dotti et al. 2006).The MBHB models that we consider here relay on twoassumptions: all bulges host a MBH, and a MBHB alwayscoalesces following the hosts’ merger. Regarding the firstassumption, dynamical processes such as gravitational re-coil and triple MBH interactions may deplete bulges fromtheir central MBH. However, if one compares the mass func-tion of coalescing binaries obtained here to the results ofEPS merger tree models, where both gravitational recoil andtriple interactions are consistently taken into account, onefinds that the two distributions (shown in the upper–leftpanel of figure 1) are in excellent agreement, within the sta-tistical uncertainties. This is because triple interactions arelikely to eject the lighter MBH from the host, leaving behinda massive binary in the vast majority of the cases (Volon-teri Haardt & Madau 2003, Hoffman & Loeb 2007). Grav-itational recoil, on the other hand, may be effective in ex-pelling light MBHs from protogalaxies at high redshift, buthas probably a negligible impact on the population of MBHsin the mass range of interest for PTA observations (Volonteri2007). The assumption that MBHBs coalesce within an Hub-ble time following the hosts’ merger is justified by severalrecent studies of MBHB dynamics. The stellar distributioninteracting with the binary may be efficiently repopulatedas a consequence of non-axisymmetric or triaxial galaxy po- c (cid:13) , 000–000 assive black holes and pulsar timing tentials (Merritt & Poon 2004, Berzcik et al. 2006), or bymassive perturbers (Perets & Alexander 2008); moreover, ifone considers post-Newtonian corrections to the binary evo-lution and the effects of eccentricity, one finds that the coa-lescence timescale is significantly reduced (Berentzen et al.2008). If the binary evolution is gas–driven, typical harden-ing timescales are expected to be shorter than 10 yr (Escalaet al. 2005, Dotti et al. 2006).In summary, we build a total of 12 models (4 MBH-hostprescriptions × M BH − σ relation (”Tr”) with single BH accretion prescrip-tion (”SA”) will be referred to as Tr-SA, see Table 1 for asummary. Assigning a MBH to each galaxy, we obtain a list of coa-lescences (labelled by MBH masses and redshift); the sameoutput quantity given by the EPS-based merger trees, usedas a starting point in PaperI. Each event in the list can benow properly weighted over the observable volume shell ateach redshift to obtain the distribution d N/d M dzdt alongthe lines described in PaperI.A further technical detail has to be considered tosmooth the coalescence distributions. The Millennium sim-ulation provides better statistics for constructing the massfunction of merging objects, but the redshift sampling israther poor. In fact, the Millennium database consists of 63snapshots of the whole simulation. The most recent ones aretaken at z = 0 , . , . , .
064 and 0.089; we need, at leastat 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 ( M , z ) dependence in the dif-ferential mass function at z < . d N/d M dzdt in the following form: d Nd M dzdt = Ψ( z ) × Z . dz d Nd M dzdt = Ψ( z ) × F ( M , t ) (7)By doing so we are redistributing according to a givenfunction Ψ( z ) the average mass function obtained at z < .
3. This is justified since, as expected, the mass func-tion does not show any significant evolution for redshiftsbelow 0 .
3. The dependence on z should be of the formΨ( z ) ∝ n G × dV C /dz , where n G is the galaxy/MBH num-ber density and dV C /dz is the differential comoving volumeshell. At such small redshifts the impact of merger activityon galaxy/MBH number density is negligible (of order of0.1/Gyr, e.g., Wake et al. 2008; White et al. 2007; Masjediet al. 2006; Bell et al. 2006); therefore, we assume n G tobe constant. On the other hand the Universe can be consid-ered Euclidean, so that the differential volume shell is justproportional to z . We then obtain the coalescing MBHBdistribution in the form d Nd M dzdt = C × z × F ( M , t ) , (8)where C is a normalization factor set by the condition Z d M Z . dz d Nd M dzdt = Z d M Z . dz C × z ×F ( M , t ) . (9) Figure 1.
Mass function of coalescing MBHBs according toMBH–host relations reported by several authors. The thick his-tograms refer to the model labeled at the top of the panel anddescribed in the text. Error bars are calculated assuming a Pois-son error in the number count of events from the coalescence cat-alogues that contribute to the chirp mass interval. In the top-leftpanel the thin histogram is the coalescing MBHB mass functionpredicted by the VHMhopk model studied in PaperI.
The M distributions of coalescing binaries are shown infigure 1 for all the MBH-host prescriptions assuming accre-tion on M only (models Tr-SA, Tu-SA, Mc-SA, La-SA ).The top-left panel also shows the distribution obtained by areference EPS merger tree model (VHMhopk; Volonteri, Sal-vaterra & Haardt 2006) used in PaperI. The agreement withthe M BH − σ prescription (”Tr”) is good for M > M ⊙ ,although the statistical errors are large due to low numberstatistics. The discrepancy for M < M ⊙ is due both tothe resolution limit of the Millennium simulation and to thefact that we relate MBHs to bulges. However, the fact thatour sample may be incomplete for M < M ⊙ has little (ifany) impact on the results of this paper: the MBHB popu-lation observable with PTAs is by far dominated by sourceswith M > M ⊙ , as we have shown in PaperI, and evenmore so when we consider systems that can be individu-ally resolved. We will further discuss this point in Section4. Note that, as discussed by, e.g., Lauer et al. 2007 andTundo et al. 2007, the high mass end of the population de-rived using the M BH − σ relation (”Tr” models) drops verysteeply. The drop is much faster than in all the other cases,that is for distributions inferred from the M BH − M bulge orthe M BH − M V relations (”Tu”, ”Mc”, ”La” models). Thisis because σ seems to converge to a plateau in the limit ofvery massive galaxies (Lauer et al. 2007). c (cid:13) , 000–000 A. Sesana et al.
The search for GWs using timing data exploits the effect ofgravitational radiation on the propagation of the radio wavesfrom one (or more) pulsar(s). The characteristic signatureof GWs on the time of arrival (TOA) of radio pulses (e.g.Sazhin 1978, Detweiler 1979, Bertotti et al. 1983) is a linearcombination of the two independent GW polarisations. Inpractice, the analysis consists in computing the differencebetween the expected and actual TOA of pulses; the timingresiduals contain information on all the effects that have notbeen included in the fit, including GWs. In this section wesummarise the observed signal produced by GWs, followingclosely Jenet et al. (2004), to which we refer the reader forfurther details.The observed timing residual generated by a GW sourcedescribed by the independent polarisation amplitudes h + , × is r ( t ) = 12 (1 + cos µ ) [ r + ( t ) cos(2 ψ ) + r × ( t ) sin(2 ψ ) , ] (10)where t is the time at the receiver, µ is the opening anglebetween the GW source and the pulsar relative to Earth and ψ is the GW polarisation angle. The two functions r + , × ( t )are defined as r + , × ( t ) = r ( e )+ , × ( t ) − r ( p )+ , × ( t ) , (11) r ( e )+ , × ( t ) = Z t dt ′ h ( e )+ , × ( t ′ ) , (12) r ( p )+ , × ( t ) = Z t dt ′ h ( p )+ , × h t ′ − dc (1 − cos µ ) i . (13)Note that r ( e )+ , × ( t ) and r ( p )+ , × ( t ) have the same functionalform, and result from the integration of the time evolution ofthe polarisation amplitudes at different times, with a delay ≃ . × ( d/ − cos µ ) yr, where d is the distanceof the pulsar from the Earth. For GWs propagating exactlyalong the Earth-pulsar direction, there is no effect on theTOAs ( r ( t ) = 0 for cos µ = ± circular orbit. We modelgravitational radiation at the leading quadrupole Newto-nian order, that is fully justified by the fact that binariesin the mass and frequency range considered here are farfrom the final merger; in fact, the time to coalescence is ≃ M / M ⊙ ) − / ( f/ × − Hz) − / yr. The timingresiduals (10) can be written as (Jenet et al, 2004): r ( t ) = r ( e ) ( t ) − r ( p ) ( t ) (14)where r ( e ) ( t ) = α ( t ) (cid:2) a + (1 + cos ι ) cos Φ( t ) + 2 a × cos ι sin Φ( t ) (cid:3) (15)In the previous expressions ι is the source inclination angle,Φ( t ) the GW phase related to the frequency f ( t ) (twice theorbital frequency) byΦ( t ) = 2 π Z t f ( t ′ ) dt ′ , (16)and α [ f ( t )] = M / D [ πf ( t )] − / ≃ . (cid:18) M M ⊙ (cid:19) / (cid:18) D
100 Mpc (cid:19) − × (cid:16) f × − Hz (cid:17) − / ns (17)is an overall scale factor that sets the size of the residuals; D is the luminosity distance to the GW source. The expressionfor r ( p ) ( t ) can be simply obtained from Equations (15), (16)and (17) by shifting the time t → t − d (1 − cos µ ) /c , seeequation (12) and (13). The two functions a + and a × arethe ”antenna beam patterns” that depend on the sourcelocation in the sky and the polarisation of the wave.MBHBs observable with PTAs produce a quasi-monochromatic signal – the frequency change is ≃ × − ( M / M ⊙ ) / ( f/ × − Hz) / nHz yr − – ofknown form (though unknown parameters). The optimaldata analysis approach to search for these signals is by meansof the well known technique of matched-filtering. The dataset can be represented as: δt ( t ) = r ( t ) + δt n ( t ) (18)where r ( t ) is the contribution from GWs, and δt n ( t ) ac-counts for the fluctuations due to noise; the latter contribu-tion is the superposition of the intrinsic noise in the measure-ments and the GW stochastic background from the wholepopulation of MBHBs. The angle-averaged optimal signal-to-noise ratio at which a signal from a MBHB radiating at(GW) frequency ≈ f can be detected using a single pulsaris h ρ i = (cid:20) δt gw ( f ) δt rms ( f ) (cid:21) . (19)In the previous expression δt rms ( f ) is the root-mean-squarevalue of the noise level δt n at frequency f , h . i represents theaverage over the source position in the sky and orientation ofthe orbital plane, and δt gw ( f ) is the characteristic amplitudeof the timing residual over the observation time T definedas: δt gw ( f ) = 815 α ( f ) p fT , (20)where the numerical pre-factor comes from the angle averageof the amplitude of the signal: (cid:10) a (1 + cos ι ) + 4 a × cos ι (cid:11) = 815 . (21)Equation 19 is appropriate to describe observations us-ing a single pulsar. In reality one can take advantage of theseveral pulsars that are continuously monitored to increasethe signal-to-noise ratio, and therefore the confidence of de-tection: adding coherently the residuals from several pul-sars – currently the Parkes PTA contains about 20 pulsars,and more are expected to available with future instruments– yields an increase in signal-to-noise ratio proportional tothe square-root of the number of pulsars used in the observa-tions. We will use the characteristic amplitude of the resid-uals δt gw to quantify the strength of a GW signal in PTAobservations; δt gw can be used to compute in a straightfor-ward way the signal-to-noise ratio, as a function of the noiselevel and number of pulsars in the array (all of which arequantities that do not depend on the astrophysical model),and therefore asses the probability of detection of sources inthe context of a given MBHB assembly scenario. c (cid:13) , 000–000 assive black holes and pulsar timing For each of the twelve models considered in Section 2 – themodels are the result of four MBH-host galaxy prescriptionsand three different accretion scenarios – that encompass avery broad range of MBHB’s assembly scenarios, we com-pute the number of sources that are potentially resolvableindividually and several statistical properties of the popu-lation, such as the redshift, chip-mass and frequency distri-butions, by means of Monte-Carlo realisations of the wholepopulation of MBHBs according to a given model. Beforedescribing the results, we provide details about each step inthe computation of the relevant quantities.The distribution given by equation (8) is straightfor-wardly converted into d N/dzd M d ln f r , i.e. the comovingnumber of binaries emitting in a given logarithmic (rest-frame) frequency interval with chirp mass and redshift inthe range [ M , M + d M ] and [ z, z + dz ], along the lines de-scribed in Section 3 of PaperI. As in PaperI, we then assumethat in the frequency range of interest ( f > × − Hz)the binary evolution is driven by GW emission only. Thisis a reasonable assumption, since the coalescence timescalefor these systems is typically ∼ < yr; any other putativemechanism (i.e. star ejection, gas torques) of angular mo-mentum removal must have an enormous efficiency to com-pete with radiation reaction on such short timescales. Foreach of the twelve models we estimate the GW stochasticbackground (and the corresponding rms value of the tim-ing residuals as a function of frequency) generated by thesources following the scheme described in Section 4 of Pa-perI. Finally we generate distribution of bright, individuallyresolvable sources by running 1000 Monte-Carlo realisationsof the whole population of MBHBs and by selecting onlythose sources whose characteristic timing residuals, equation( 20), exceed the stochastic background level. We note thatthe result of which and how many sources raise above the av-erage noise level depends on the duration of the observation T , for two reasons: (i) T affects the size of the observationalwindow in frequency space, in particular the minimum fre-quency 1 /T that can be reached, and (ii) the backgroundlevel decreases as the observation time increases (as the sizeof the frequency resolution bin ∆ f = 1 /T decreases), en-hancing the number of individually resolvable sources. Forthe results that are described here and summarised in table1, we set T = 5 yrs; increasing the data span to 10 years,the background level would be slightly lower, and the ad-ditional resolved sources would be barely brighter than thebackground. The statistics of bright, well resolvable sourcesis basically unaffected. We can then cast the results in termsof the cumulative number of resolvable sources as a functionof the timing residuals, according to: N ( δt gw ) = Z ∞ δt gw dNδt ′ gw δt ′ gw , (22)where the integral is restricted to the sources that produceresiduals above the rms level of the stochastic background;if we do not consider this additional constraint, then, forany given δt gw , N ( δt gw ) simply returns the total number ofsources in the Monte-Carlo realisation above that particularresidual threshold.Each Monte-Carlo realisation clearly yields a differentvalue for N ( δt gw ) (or its distribution according to a given Figure 2.
Summary of the properties of the population of mas-sive black hole binary systems – according to model Tu-SA –that generate gravitational waves in the frequency window cov-ered by Pulsar Timing Arrays.
Top-left panel : characteristic am-plitude of the timing residuals δt gw (equation( 20)) as a functionof frequency; the asterisks are the residuals generated by individ-ual sources and the solid line is the estimated level of the GWstochastic background. Top-right panel : δt gw as a function of thenumber N ( δt gw ) of total (dotted line) and individually resolvable(solid line) sources, see equation 22. Bottom-left panel : Distri-bution of the number of total (dotted lines) and resolvable (solidlines) sources per logarithmic frequency interval dN ( δt gw ) /d log f as a function of the GW frequency for different values of δt gw :from top to bottom 1, 10 and 100 ns, respectively. Bottom-rightpanel : distribution of the total (dotted lines) and individuallyresolvable (solid lines) number of sources per logarithmic chirpmass interval dN ( δt gw ) /d log M as a function of chirp mass fordifferent values of δt gw : from top to bottom 1, 10 and 100 ns,respectively. An observation time of 5 years is assumed. parameter); the values quoted in the next section and sum-marised in table 1 refer to the sample mean computed overthe set of Monte-Carlo realisations and the sample standarddeviation. Figure 4 quantifies the typical 1- σ error in ourestimate of the number of sources. The large number of Monte-Carlo realisations allows usto study the details of the properties of the individualsources in a statistical sense. We concentrate in particularon the physical properties of the population, such as the ex-pected number of sources per logarithmic frequency interval dN ( δt gw ) /d log f , chirp mass range dN ( δt gw ) /d log M andredshift dN ( δt gw ) /dz , and the observable parameters, suchas the timing residuals produced by each system and theoverall expected number of resolvable MBHBs at a givenlevel of timing residual noise. A summary of the typicalrange of information that can be extracted from the sim-ulations is shown in figure 2 for the specific model Tu-SA. c (cid:13) , 000–000 A. Sesana et al.
Figure 3.
Left panel : The effect of the MBH-host galaxy relation, assuming that accretion always takes place for a single black holebefore merger (”SA” models), on the number of observable systems. The plot shows the number of total (thin lines) and resolvable(thick lines) sources N ( δt gw ) as a function of δt gw , see equation 22. Four different MBH merger scenarios are considered: Tu-SA (solidline), Tr-SA (long–dashed lines), Mc-SA (short–dashed lines) and La-SA (dotted lines), see Section 2 for a description of the models. Right panel : The effect of the MBH accretion model on the number of individually observable systems. The plot shows the number ofresolvable sources only, N ( δt gw ) as a function of δt gw . As reference for the MBH-host galaxy relation, models ”La” (thick lines) and ”Tr”(thin lines) are considered. The line style is as follow: model La-SA and Tr-SA (solid lines), La-DA and Tr-DA (short–dashed lines), andLa-NA and Tr-NA (long–dashed lines). The duration of the observation is set to T = 5 yr The top-left panel shows the induced residuals of each sin-gle source compared to the level produced by the stochas-tic background from the whole population; the plot clearlyshows the importance of taking into account the additional”noise contribution” from the brightness of the GW sky inconsidering the detectability of resolvable systems. Thereare many sources inducing residuals above, say the 5 nslevel, however most of them contribute to the build-up ofthe background, and are not individually resolvable. Theexpected number of bright resolvable MBHBs at frequencies < − Hz, and at a timing level > δt gw from 1000 Monte-Carlo realizations of the emitting population. In this partic-ular case a sensitivity of ≈
10 ns is required to resolve anindividual source; for a timing precision of 100 ns there isa 5% chances to observe a particularly bright source. Notethat at the 1 ns level, there are ∼
100 MBHBs contribut-ing to the signal, however 90% of them contribute to thebackground and only about 10 sources are individually re-solvable. In the two bottom panels of figure 2 we plot thefrequency and chirp mass distributions of resolvable sourcesfor selected values of the minimum detectable residual am-plitude δt gw . Not surprisingly, the chirp mass of observablesystems decreases for smaller values of the considered resid-ual threshold, however even for an rms level of 1 ns all thesystems are characterised by M > M ⊙ . The frequencydistribution shows instead a peak corresponding to the fre- quency at which the background level equals the selectedvalue of δt gw . At higher frequencies the number of sourcesdrastically drops because the number of emitting binaries isa quite steep function of frequency, N ( f ) ∝ f − / ; the de-crease at lower frequencies is because most of the emittersactually contribute to the background and are not individ-ually resolved (as clearly shown by the dotted lines). Thequalitative behaviour of the results obtained using differentastrophysical models is very similar to the one describedin figure 2, with differences that affect only the numericalvalues of the different quantities.The central question that we want to address in thispaper is what is the expected number of individually resolv-able sources that produce an effective timing residual abovea given value, as a function of different models of MBHBformation and evolution. We summarise these results in fig-ures 3 and 4, and in Table 1, where we show the mean totalnumber of individual sources that exceed a given level oftiming residual (as defined by equation (22)), as a functionof the timing residual. The qualitative behaviour of the re-sults is similar for all the scenarios, but the actual numbersvary significantly. In fact, both MBH-host relations and ac-cretion prescriptions have a strong impact on the statisticsof the bright sources and, consequently, on their detection.We analyse the effect of the black holes populations and ofthe accretion in turn. The left panel of figure 3 shows re-sults where the accretion prescription is the same (”SA”),but the underlying MBH-host relation changes; on the otherhand, in the right panel, we select two MBH-host relations, c (cid:13) , 000–000 assive black holes and pulsar timing Figure 4.
The number of individually resolvable sources for se-lected MBHBs assembly models. The plots show N ( δt gw ) as afunction of δt gw , see equations 22 and 20, for a typical model(Tu-SA, top panel), and, in the lower panel, the two models thatyield the largest (La-DA: top curve) and smallest (Tr-NA: bottomcurve) number of sources. The solid lines represent the mean valueof N ( δt gw ) and the shaded area the 1 − σ region as computed from1000 Monte-Carlo realisations of each MBHB population, see alsoTable 1. Figure 5.
Top panel
Cumulative mean number of resolvablesources N ( δt gw ) as a function of the characteristic timing resid-ual δt gw for different mass cuts: solid line: M > M ⊙ ; dashedline: 10 M ⊙ < M < M ⊙ ; dotted line 10 M ⊙ < M < M ⊙ . Bottom panel : same as top panel but for different red-shift intervals: solid line z < .
1; dashed line 0 . < z <
1; dottedline z > Figure 6.
Distribution of the number of resolvable sources asa function of the gravitational wave frequency. The plots show dN ( δt gw ) /d log f for different values of the characteristic ampli-tude of the timing residuals (from right to left δt gw = 1 , , Figure 7.
Chirp mass distribution of the number of resolvablesources. The plots show dN ( δt gw ) /d log M for the same charac-teristic amplitude of the timing residuals and models as in fig-ures 6 and 8.c (cid:13) , 000–000 A. Sesana et al.
Figure 8.
Redshift distribution of the resolvable sources. Theplots show dN ( δt gw ) /d log z for the same characteristic amplitudeof the timing residuals and models as in figures 6 and 7. but change the accretion history. The left panel shows thatassuming a sensitivity threshold of 30 ns, one expects to ob-serve of the order of one source in the La-SA model, whilethere is only a probability ≈
5% for the Tr-SA model. Ifwe maintain the same MBH-host relation and we considerdifferent accretion scenarios for e.g. the Lauer et al. popula-tion, the mean number of expected sources varies by a fac-tor of ≈ . . ≈ ≈ − ≈
10 uncertainty in thebackground level estimated in PaperI.The typical spread around the mean values obtained inthe Monte-Calro realisations is shown in figure 4 for selectedmodels. When N ( δt gw ) ≫ σ range is roughly thePoisson error around the mean (reflecting the uncorrelatednature of sources in each Monte Carlo realisation). When N ( δt gw ) < δt gw value if the actualMBHB population in the Universe follows the prescriptiongiven by the considered model; in this case a non-detectionis trivially consistent with the model predictions. Table 1provides a summary of the results for the 12 models.It is also interesting to investigate the mass-redshift dis-tribution of the detectable sources. Figure 5 shows the ex-pected number of detectable individual sources (as a func-tion of the residuals amplitude) for different redshift andmass ranges. Obviously, higher timing residuals correspondto higher M , since the strength of the signal is proportionalto M / , and the most likely sources to be detected fallin the range 10 M ⊙ < ∼ M < ∼ M ⊙ (MBHBs with M > Figure 9.
The characteristic amplitude of the GW stochasticbackground from the population of MBHB systems. In each panelthe thin lines identify the estimated level of the stochastic back-ground assuming “SA” (solid line), “DA” (dashed line) and “NA”(dotted line) accretion modes. The total GW amplitude from asingle Monte–Carlo realisation of the signal corresponding to the“SA” accretion mode is also shown as thick solid line. M ⊙ produce indeed larger timing residuals, but are alsomuch rarer). Interestingly, the vast majority of detectablesources are at redshift 0 . < ∼ z < ∼
1, which shows that PTAscould probe the medium-redshift Universe, and are unlikelyto discover nearby sources. The reason is simply that, atleast at small redshift, the Universe volume increases as z .A summary of the properties of individual re-solvable sources, dN ( δt gw ) /d log f , dN ( δt gw ) /d log M and dN ( δt gw ) /d log z , is given in figures 6, 7 and 8, respectively,for all the four MBHB population models considered here,with accretion limited to a single black hole prior to mergerconsidering different residual thresholds δt gw = 1 , , , ≈ × M ⊙ for allmodels assuming δt gw = 1 ns. Increasing δt gw , the distribu-tion peak shifts towards higher M and the mean numberof events is strongly model dependent for δt gw >
10 ns, cf.the values in Table 1. The redshift distributions (figure 8)consistently show a broad peak in the range 0 . < z < δt gw increases, because in this model similar galaxies arepopulated by more massive black holes if found at higherredshifts, see equation (2). c (cid:13) , 000–000 assive black holes and pulsar timing Figure 10.
The characteristic amplitude of the GW stochasticbackground compared to the estimate given in PaperI. The thickdashed line is the stochastic background predicted by the Tr-SAmodel; the dotted lines bound the background levels computed forall the models investigated in this paper. For comparison, the solidthick line is the background predicted by the VHMhopk modeland the shaded area the range of uncertainty of the strength ofthe signal, as reported in PaperI (see Sections 4 and 5 of PaperIfor details).
As a sanity check, we compare the stochastic backgroundsderived according to the MBH populations inferred usingthe Millennium Simulation to the predictions of EPS–basedmodels reported in PaperI (the reader is referred to PaperIfor the technical details). In figure 9 we show Monte Carlogenerated signals for each model; in each panel we plot thestochastic levels according to the three different accretionmodes discussed in Section 2. Both the accretion prescrip-tion and the adopted MBH-host correlation influence thelevel of the background. If accretion occurs onto the rem-nant (i.e. after coalescence, ”NA” models), the predictedcharacteristic amplitude of the GW background can be upto a factor of 3 lower with respect to models in which boththe MBHs accrete before the final coalescence (”DA” mod-els); on the other hand, M BH − σ (”Tr”) models predictlower backgrounds compared to M BH − M bulge (”Tu”, ”Mc”)and M BH − M V (”La”) models. A comparison of these re-sults with those presented in PaperI is given in figure 10. At10 − Hz the models studied here cover a characteristic am-plitude range consistent with the uncertainty estimated inPaperI. The Tr-SA model predicts a stochastic backgroundthat agrees with the typical EPS–model within 30% for f < − Hz. The high frequency end is instead steeper; thiseffect is caused by incompleteness in the low mass end of theMBH population. As shown in figure 1, the mass function ofcoalescing MBHB derived from the Millennium simulationis not consistent with the one obtained using the EPS for-malism for M < M ⊙ . The weight of a single dark matter particle in the Millennium simulation is 8 . /h × M ⊙ ,allowing the reconstruction of haloes with minimum mass ofthe order of ≈ × M ⊙ . Assuming a barionic fractionof 0.1, the simulation is then incomplete for barionic struc-tures less massive than ≈ × M ⊙ . We checked this byplotting the mass function of barionic structures and find-ing a sudden drop below 10 M ⊙ . It is then inevitable thatin the results derived from the Millennium Simulation mostof the MBHs with mass below a few × M ⊙ are missing.Since many of these MBHs are expected to merge with moremassive ones during cosmic history, the (spurious) lack ofMBHs in this mass rage explains the flattening of the massfunction d ˙ N M /d log M shown in the top-left panel of figure1. All the backgrounds are rather similar at f > − Hzbecause all the MBH prescriptions adopted lead to similarMBH mass functions at M BH < M ⊙ (this fact is inde-pendent of the incompleteness issue). This means that theslope of the background on the right of the knee has a welldefined dependence upon the adopted MBH population: themore pronounced is the high mass tail of the MBH massfunction, the steeper is the high frequency end of the GWbackground. In turn, models constructed using the Millen-nium simulation confirm the findings of PaperI. We have investigated the ability of Pulsar Timing Arrays(PTAs) to resolve individual massive black hole binary sys-tems by detecting gravitational radiation produced duringthe in-spiral phase through its effect on the residuals of thetime of arrival of signals from radio pulsars. We have con-sidered a broad range of assembly scenarios, using the dataof the Millennium simulation to evaluate the galactic haloesmerger rates, and a total of twelve different models thatcontrol the relations between the mass of the central blackholes and the galactic haloes, and the evolution of the blackhole masses through accretion. These models therefore coverqualitatively (and to large extent quantitatively) the wholespectrum of MBH assembly scenarios currently considered.Regardless of the model, we estimate that at least one re-solvable source is expected to produce timing residuals inthe range ∼ −
50 ns, and therefore future PTAs, and inparticular the Square Kilometre Array may be able to ob-serve these systems. A whistle-stop summary of the mod-els and results is contained in Table 1. The total numberof visible events clearly depend on the sensitivity of PTAsand on the astrophysical scenario. As expected, the bright-est sources (for PTAs) are very massive binaries with chirpmass M > × M ⊙ . However, (initially) quite surprisinglymost of the resolvable sources are located at relatively highredshift ( z > . c (cid:13) , 000–000 A. Sesana et al. of this paper and PaperI, it shows that we can now havein hand self-consistent predictions for stochastic and deter-ministic signals from the cosmic population of MBHBs, andsuggests that EPS merger trees could provide a valuable ap-proach to the studies of MBH evolution at low-to-mediumredshift.As a final word of caution, we would like to stress thatthe results of this paper clearly suffer from considerable un-certainties determined by the still poor quantitative infor-mation about several parameters that control the models.The spread of the predictions of the expected events is there-fore likely dominated by the lack of knowledge of the modelparameters, rather than the differences between the assem-bly scenarios.
REFERENCES
Bell E. F., Phleps S., Somerville R. S., Wolf C., Borch A. &Meisenheimer K., 2006, ApJ, 652, 270Bender, P. et al., 1998,
LISA , Laser Interferometer Space Antennafor gravitational wave measurements: ESA Assessment StudyReportBerczik P., Merritt D., Spurzem R. & Bischof H. P., 2006, ApJ,642, 21Berentzen I., Preto M., Berczik P., Merritt D. & Spurzem R.,2008, arXiv0812.2756Bertone S., De Lucia G. & Thomas P. A., 2007, MNRAS, 379,1143Bertotti B., Carr B. J. & Rees M. J., 1983, MNRAS, 203, 945Cox, T. J., Jonsson, P., Somerville, R. S., Primack, J. R., & Dekel,A. 2008, MNRAS, 384, 386De Propris R., Conselice C. J., Driver S. P., Liske J., Patton D.,Graham A. & Allen P., 2007, ApJ, 666, 212Detweiler S., 1979, ApJ, 234, 1100Dotti, M., Colpi, M., & Haardt, F. 2006, MNRAS, 367, 103Escala A., Larson R. B., Coppi P. S. & Mardones D., 2005, ApJ,630, 152Faber S. M. & Jackson R. E., 1976, ApJ, 204, 668Haehnelt M. G., 1994, MNRAS, 269, 199Hellings R. W. & Downs G. S., 1983, ApJ, 265, 39Hoffman L. & Loeb A., 2007, MNRAS, 377, 957Jaffe A. H. & Backer D. C., 2003, ApJ, 583, 616Janssen G. H., Stappers B. W., Kramer M., Purver M., Jessner A.,Cognard I., 2008, in ”40 YEARS OF PULSARS: MillisecondPulsars, Magnetars and More”, AIP Conference Proceedings,983, 633Jenet, F. A., Lommen, A., Larson, S. L., & Wen, L. 2004, ApJ,606, 799Jenet F. A., Hobbs G. B., Lee K. J. & Manchester R. N., 2005,ApJ, 625, 123Jenet F. A. et al., 2006, ApJ, 653, 1571Koushiappas S. M. & Zentner, A. R. 2006, ApJ, 639, 7Lacey C. & Cole S., 1993, MNRAS, 262, 627Lauer T. R., Tremaine S., Richstone D. & Faber S. M., 2007, ApJ,670, 249Lin L. et al., 2004, ApJ, 617, 9Malbon, R. K., Baugh, C. M., Frenk, C. S., & Lacey, C. G. 2007,MNRAS, 382, 1394Manchester R. N., 2008, in ”40 YEARS OF PULSARS: Millisec-ond Pulsars, Magnetars and More”, AIP Conference Proceed-ings, 983, 584Marulli F., Crociani D., Volonteri M., Branchini E. & MoscardiniL., 2006, MNRAS, 368, 1269Merritt D. & Poon M. Y., 2004, ApJ, 606, 788Masjedi M. et al., 2006, ApJ, 644, 54 Mayer, L., Kazantzidis, S., Madau, P., Colpi, M., Quinn, T., &Wadsley, J. 2007, Science, 316, 1874McLure R. J., Jarvis M. J., Targett T. A., Dunlop J. S. & BestP. N., 2006, MNRAS, 368, 1395Phinney E. S., 2001, arXiv:astro-ph/0108028Press W. H. & Schechter P., 1974, ApJ, 187, 425Rajagopal M. & Romani R. W., 1995, ApJ, 446, 543Sazhin M. V., 1978, Soviet Astron., 22, 36Sesana A., Haardt F., Madau P. & Volonteri M., 2004, ApJ, 611,623Sesana A., Haardt F., Madau P. & Volonteri M., 2005, ApJ, 623,23Sesana, A., Haardt, F., & Madau, P. 2007, ApJ, 660, 546Sesana A., Vecchio A. & Colacino C. N., 2008, MNRAS, 390, 192Sheth R. K. & Tormen G., 1999, MNARS, 308, 119Springel V. et al., 2005, Nature, 435, 629Tremaine S. et al., 2002, ApJ, 574, 740Tundo E., Bernardi M., Hyde J. B., Sheth R. K. & Pizzella A.,2007, ApJ, 663, 57Volonteri M., Haardt F. & Madau P., 2003, ApJ, 582, 599Volonteri M., Salvaterra R. & Haardt F., 2006, MNRAS, 373, 121Volonteri M., 2007, ApJ, 663, 5Wake D. A. et al. 2008, MNRAS, 387, 1045White M., Zheng Z., Brown M. J. I., Dey A., Jannuzi B. T., 2007,ApJ, 655, L6Wyithe J.S.B. & Loeb A., 2003, ApJ, 590, 691Yoo J., Miralda-Escud´e J., Weinberg D. H., Zheng Z., Morgan,C. W. 2007, ApJ, 667, 813Yu, Q. 2002, MNRAS, 331, 935c (cid:13)000