Gravitational wave emission from binary supermassive black holes
Massive black hole binaries (MBHBs) are unavoidable outcomes of the hierarchical structure formation process, and according to the theory of general relativity are expected to be the loudest gravitational wave (GW) sources in the Universe. In this article I provide a broad overview of MBHBs as GW sources. After reviewing the basics of GW emission from binary systems and of MBHB formation, evolution and dynamics, I describe in some details the connection between binary properties and the emitted gravitational waveform. Direct GW observations will provide an unprecedented wealth of information about the physical nature and the astrophysical properties of these extreme objects, allowing to reconstruct their cosmic history, dynamics and coupling with their dense stellar and gas environment. In this context I describe ongoing and future efforts to make a direct detection with space based interferometry and pulsar timing arrays, highlighting the invaluable scientific payouts of such enterprises.
pacs:04.70.-s – 98.65.Fz – 04.30.-w – 04.30.Db – 04.30.Tv – 04.80.Nn
Today, massive black holes (MBHs) are ubiquitous in the nuclei of nearby galaxies , and we see them shining as quasars along the whole cosmic history up to redshift . In the last decade, MBHs were recognized as fundamental building blocks in hierarchical models of galaxy formation and evolution, but their origin remains largely unknown. In fact, our current knowledge of the MBH population is limited to a small fraction of objects: either those that are active, or those in our neighborhood, where stellar- and gas-dynamical measurements are possible. According to the current paradigm, structure formation proceeds in a hierarchical fashion , in which massive galaxies grow by accreting gas through the filaments of the cosmic web and by merging with other galaxies. As a consequence, the MBHs we see in today’s galaxies are expected to be the natural end-product of a complex evolutionary path, in which black holes (BHs) seeded in proto-galaxies at high redshift grow through cosmic history via a sequence of MBH-MBH mergers and accretion episodes [4, 5]. In this framework, a large number of MBH binaries naturally form following the frequent galaxy mergers.
According to Einstein’s theory of General Relativity, accelerating masses cause modifications of the spacetime that propagate at the speed of light, better known as gravitational waves (GWs). However, in Einstein’s equations, the matter-metric coupling constant is of the order of (where is the gravitational constant and is the speed of light), which is of the order ! As a matter of fact the spacetime is extraordinarily stiff, therefore only massive, compact astrophysical object can produce a sizable strain that would be observable with advanced technology . Two MBHs orbiting each other in a bound binary (MBHB) system carry a huge time varying quadrupole momentum and are therefore expected to be the loudest gravitational wave (GW) sources in the Universe . The frequency spectrum of the emitted radiation covers several order of magnitudes, from the sub-nano-Hz up to the milli-Hz. The Hz window is going to be probed by spaceborne interferometers like the recently proposed European eLISA [8, 9, 10]. At Hz, joint precision timing of several ultrastable millisecond pulsars (i.e. a pulsar timing array, PTA) provides a unique opportunity to get the very first low-frequency detection. The European Pulsar Timing Array (EPTA) , the Parkes Pulsar Timing Array (PPTA)  and the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) , joining together in the International Pulsar Timing Array (IPTA) , are constantly improving their sensitivities, getting closer to their ambitious target.
This focus issue contribution aims at covering all the relevant aspects of MBHBs intended as GW sources, in the spirit of providing a broad overview. Given the extent of the topic, we will just skim through its several facets, providing the appropriate references for in-depth reading. In Section 2 we introduce the concept of GWs at a very basic level, defining the relevant astrophysical scales implied by MBHBs. A general overview of MBH formation and evolution (both their masses and spin), together with a brief description of MBHB dynamics is provided in Section 3. We return on GWs in Section 4, where we describe in deeper detail the expected signal from MBHBs and its dependencies on the relevant parameter of the system. There, we also introduce some basics of parameter estimation theory, showing how the rich astrophysical information enclosed in individual signals can be recovered. Section 5 is then devoted to the scientific payouts of GW detection both in the milli-Hz and in the nano-Hz regime. We wrap-up in Section 6 with some brief conclusion remarks.
2 Gravitational waves: basics
The existence of GWs was one of the first predictions of Einstein’s General Relativity (GR), since they arise as natural solutions of the linearized Einstein equations in vacuum. Expanding the metric tensor as
where represent the Minkowski flat metric and , and switching to the appropriate Lorentz gauge, the perturbation satisfies
where is the d’Alambertian operator and the source term is the stress-energy tensor. Equation 2 represents a set of wave equations and therefore admits wave solutions. These solutions are ripples in the fabric of the spacetime propagating at the speed of light: GWs. GWs are transverse, i.e. they act in a plane perpendicular to the wave propagation, and (at least in GR) have two distinct polarizations, usually referred to as and (see cartoon in  and Section 4.2). By expanding the mass distribution of the source into multipoles, conservation laws enforce GWs coming from the mass monopole and mass dipole to be identically zero, so that the first contribution to GW generation comes from the mass quadrupole moment . The GW amplitude is therefore proportional to the second time derivative (acceleration) of . Moreover, energy conservation enforces the amplitude to decay as the inverse of the distance to the source, . A straightforward dimensional analysis shows that the amplitude of a GW is of the order of 
In order to generate GWs we therefore need accelerating masses with a time varying mass-quadrupole moment. The prefactor implies that these waves are tiny, so that the only detectable effect is produced by massive compact astrophysical objects.
We now specialize to the case of a binary astrophysical source. From now on, we shall use Geometric units ; in such units s and pcs. For sources at cosmological distances (i.e., with non negligible redshift ), all the following equations are valid if all the masses are replaced by their redshifted counterparts (e.g., ), the distance is taken to be the luminosity distance (i.e., ) and is kept to be the observed GW frequency (the frequency in the source emission restframe is ). Keeping this in mind, consider a binary system of masses , mass ratio and total mass , in circular orbit at a Keplerian frequency at a distance to the observer. A detailed calculation (in the quadrupole approximation) shows that the system emits a monochromatic wave at a frequency , with inclination–polarization averaged GW strain given by :
where we introduced the chirp mass . For a pair of Schwarzschild BHs, the maximum frequency of the wave is emitted at the innermost stable circular orbit (ISCO) and can be written as:
where . GWs carry away energy from the system, with total luminosity given by:
Equating the energy loss to the shrinking of the binary semimajor axis ,
and converting into frequency using Kepler’s law yields
The integral of equation (8) from to defines the remaining lifetime of a binary emitting at a frequency before its final coalescence.
Putative eccentricity plays an important role in the evolution of the binary and in the emitted GWs. In this case the luminosity, at a fixed MBHB semimajor axis, is boosted to , where is given by equation (6) and
Accordingly, the evolution of the binary orbit is a factor faster than in the circular case. The energy is radiated in form of a rather complicated GW spectrum, covering the spectral range , where is an integer index (see Section 4.2 for more details). In particular, the emission is stronger close to the binary periastron (there, the acceleration is larger, and so is the derivative of the quadrupole moment of the source), which leads to efficient circularization according to
Examples of GW driven circularization are illustrated in the MBHB evolutionary paths shown in panels ’a3’ and ’b3’ of figure 2.
Normalizing equation (4) to typical astrophysical MBHB values gives a strain of
If we now estimate the frequency change in an observation time as , we find
where now yr. Equations (11), (12) and (13) define the properties of typical MBHB signals. A light binary spans a frequency range of Hz, in a Gyr before the coalescence. For a source at 1Gpc, the strain is at Hz, and , implying a chirping signal rapidly sweeping through the frequency band during a putative observation. On the other hand, a massive binary covers a frequency range of Hz, in a Gyr before the coalescence. For a source at 1Gpc, the strain is at Hz, and Hz, implying a non evolving, monochromatic source. It is likely that many such sources accumulates at these low frequency, resulting in an incoherent superposition of monochromatic waves. MBHBs are therefore loud primary GW targets in the nano-Hz–milli-Hz regime, and they can manifest themselves both as rapidly chirping signals (at milli-Hz frequencies) as well as incoherent superposition of monochromatic waves (at nano-Hz frequencies).
3 Massive black hole binaries
The mechanism responsible for the formation of the first seed BHs is not well understood. These primitive objects started to form at the onset of the cosmic dawn, around , according to current cosmological models . At an epoch of , the earliest stars formed in small, metal-poor protogalactic halos may have had masses exceeding 100 , ending their lives as comparable stellar mass BHs, providing the seeds that would later grow into MBHs . However, as larger, more massive and metal enriched galactic discs progressively formed, other paths for BH seed formation became viable (see  for a review). Global gravitational instabilities in gaseous discs may have led to the formation of quasi-stars of that later collapsed into seed BHs . Alternative scenarios are the collapse of massive stars formed in run-away stellar collisions in young, dense star clusters  or the collapse of unstable self-gravitating gas clouds in the nuclei of gas-rich galaxy mergers at later epochs . Thus, the initial mass of the seeds remains one of the largest uncertainties in the present theory of MBH formation.
However, once formed, these seed BHs inevitably took part in the hierarchical structure formation process, growing along the cosmic history through a sequence of mergers and accretion episodes [4, 5]. Figure 1 shows examples of expected MBHB merger rates as a function of redshift for a sample of selected seed BH models. The uncertainty is large, with numbers ranging from ten to several hundreds events per year. Multiplying by the Hubble time and dividing by the number of galaxies within our Hubble horizon (), figures imply that each galaxy experienced few to few hundred mergers in its past life, placing mergers among the crucial mechanisms in galaxy evolution.
3.1 Mass and spin evolution
Astrophysical BHs are extremely simple objects, described by two quantities only, namely their mass, , and angular momentum, 111We use boldface to express vectors and standard font to express their magnitude.. The magnitude of the latter can be expressed by the dimensionless parameter 222Not to be confused with the binary semimajor axis, it will be clear case by case when we refer to one or the other.. By definition . Along the cosmic history, MBH mass and spin inevitably evolve according to three principal evolution mechanisms: (i) merger with other MBHs, (ii) episodic accretion of compact objects, disrupted stars, or gas clouds, and (iii) prolonged accretion of large supplies of gas via accretion disks. As shown by , coalescences of MBHs with random spin directions result in a broad remnant spin distribution; in particular highly spinning MBHs tend to spin-down. Despite the important of MBH-MBH mergers, The dominant role in the mass and spin evolution of MBHs can be attribute to accretion. Continuous Eddington limited accretion implies an exponential mass growth , where Gyr and is the mass-radiation conversion efficiency ( for ). If this happens in a coherent fashion through, e.g., a thin disk  , an initially Schwarzschild BH becomes maximally spinning after accreting an amount of mass of the order of . However high spins imply , considerably slowing down the mass growth, making it impossible to produce a MBH of at (i.e. in yrs). The problem is avoided if mass is accreted in a series of small incoherent packets (chaotic accretion ). In this case, depending on the angular momentum of the accreted material, the MBH is spun up or down, performing a random walk in spin magnitude that keep it close to zero. However this is true only if the angular momentum direction of the packets is nearly isotropically distributed on the sphere. Real galaxies usually show large coherent gas structures, and a significant amount of rotation (see, e.g., [29, 30]). If the spin vectors of the accreting packets have, on average, a preferential direction (i.e., they angular momenta do not sum up to zero), then the spin evolution is more complicated, and high spin values might still be preferred , as shown in the left panel of figure 2. Rapid mass growth is difficult to reconcile with measurements of high spins (although the latter involve galaxies in the local Universe ), and the requirement of high spins to power energetic relativistic jets in many theoretical models (e.g., ), and a complete joint understanding of the MBH mass and spin evolution is still missing.
3.2 Massive black hole binary dynamics
MBHs become loud sources of GWs when they are in bound, sub–pc binaries. Forming after galaxy mergers, those binaries sit at the center of the stellar bulge of the remnant, possibly surrounded by massive gas inflows triggered by dynamical instabilities related to the strong variations of the gravitational potential during the merger episode . The interaction with the environment imprints distinctive signatures in the binary orbital elements and in the individual spins of the holes. We will see later how this information can be recovered by GW observation, which will therefore allow to directly probe the complex physics underlying the evolution of these spectacular objects.
3.2.1 Stellar driven binaries.
Ignoring technical details related to the ’loss cone evolution’ (see D. Merritt contribution to this issue), a background of stars scattering off the binary drives its semimajor axis evolution according to the equation 
where is the density of the background stars, is the stellar velocity dispersion and is a numerical coefficient of order 15. The eccentricity evolution in stellar environments has been tackled by several authors by means of full N-body simulations. In general, equal mass, circular binaries tend to stay circular or experience a mild eccentricity increase , while binaries that form already eccentric, or with (regardless of their initial eccentricity) tend to grow more eccentric [38, 39], in reasonable agreement with the prediction of scattering experiments [36, 40]. The interaction with stars is unlikely to significantly affect the individual spins of the holes. Therefore, star driven binaries are expected to grow to quite high eccentricities, while the spins of the individual holes are likely randomly oriented.
3.2.2 Gas driven binaries.
In the case of circumbinary disks, the detailed evolution of the system depends on the complicated and uncertain dissipative physics of the disk itself. The simple case of a coplanar prograde circumbinary disk admits a selfconsistent, non stationary solution that was derived by . In this case, the binary semimajor axis evolution can be approximated as [41, 42]
Here, is the mass accretion rate at the outer edge of the disk, is the semimajor axis at which the mass of the unperturbed disk equals the mass of the secondary MBH, and is the reduced mass of the binary. In the circumbinary disk scenario, eccentricity excitation has been seen in several simulations [43, 44]. In particular, the existence of a limiting eccentricity has been found in , in the case of massive selfgravitating disks. If the accretion flow is coherent, and and the spins of the two MBHs are misaligned, the Bardeen-Petterson effect  will act to align to in a very short timescale (yr [47, 48]). Therefore, in gaseous rich environments, mildly eccentric binaries might be the norm, and the MBH individual spins tend to align with the orbital angular momentum.
Compared to the GW driven case, , equations (14) and (15) have a very different (milder and positive) dependence. Therefore, equating equations (14) and (15) to gives the transition frequency between the external environment driven and the GW driven regimes. For typical astrophysical systems one gets:
We therefore see that very massive nano-Hz MBHBs might still be influenced by their environment, and therefore have high eccentricities (, see panel ’b3’ in figure 2). Even though GW emission efficiently circularizes binaries (see Section 2), systems in the milli-Hz range can still retain substantial residual eccentricities (, see panel ’a3’ in figure 2).
4 MBHB waveforms
Having introduced the basics of GW emission from a binary system in Section 2, we turn now in some more detail to the gravitational waveform modeling. In particular we show how eccentricity and spins affect the detectable GW signal and we describe the basic theory of information recovery, that enables us to dig out the parameters of the source from the detected waveform.
4.1 The stages of the binary coalescence
The evolution of MBHBs is customarily divided into three phases: inspiral, merger, and ring-down . The inspiral is a relatively slow, adiabatic process. Different techniques have been employed to describe this stage, ranging from classic Post Newtonian (PN) expansions of the energy–balance equation , to non–adiabatic resummed methods in which the equations of motion are derived from an effective one body (EOB) relativistic Hamiltonian . A detailed description of such methods is beyond the scope of this article, and an excellent overview can be found in . The inspiral is followed by the dynamical coalescence, in which the MBHs plunge and merge together, forming a highly distorted, perturbed remnant. At this stage, all analytical approximations break down, and the system can only be described solving by directly the Einstein equation using numerical simulations [54, 55, 56]. The distorted remnant settles into a stationary Kerr BH as it rings down, by emitting gravitational radiation. This latter stage can be, again, modeled analytically using BH perturbation theory . An example of the full waveform with the identification of the various stages is given in figure 3.
In recent years there has been a major effort in constructing accurate waveforms inclusive of all three phases. ”Complete” waveforms can be designed by stitching together analytical PN waveforms for the early inspiral with a (semi)phenomenologically described merger and ring-down phase calibrated against available numerical data (known as PhenomB-PhenomC waveforms, ). Alternatively, complete waveforms can be constructed within the EOB formalism by adding free parameters to be calibrated against NR simulations and by attaching a series of damped sinusoidals describing the ringdown (known as EOBNR waveforms [59, 60]). A detailed overview is given in . What is relevant to our discussion is that the full evolution of MBHBs can be tackled with a combination of analytical and numerical methods, and accurate waveforms encoding all the parameters of the system can be computed. In the following we concentrate on the inspiral signal only, which is the richest in terms of encoded information.
4.2 The adiabatic inspiral: impact of eccentricity and spin
For circular binaries, the evolution of the adiabatic inspiral is completely determined by the energy-balance equation that relates the derivative of the energy function to the gravitational flux radiated away 333For eccentric binaries, an angular momentum balance equation must also be imposed to compute the evolution in eccentricity.
from which we may derive the binary acceleration and phase evolution as:
For orbital velocities , and can be expanded in powers of to a given order in . This results in a corresponding expansion for the binary acceleration of the form :
where , , and are the Newtonian, (post)1-Newtonian, and (post)2-Newtonian contributions to the equations of motion, is the contribution due to the radiation-reaction force, and and are the spin-orbit and spin-spin coupling contributions. This is reflected in the radiated wave, that can be similarly expanded in the form:
where is the standard quadrupole moment of the source, define the spatial components of the perturbation tensor, and the subscripts have the same meaning as in equation (19). Choosing the appropriate orthonormal radiation frame, can be written as two independent polarizations only; and . To the leading quadrupole order, in the circular case one obtains the familiar form
where (usually referred as inclination) is the angle defined by the line of sight with respect to the orbital angular momentum vector, , and as defined in Section 2.
The eccentricity enters directly in the computation of since it affects the velocity of the MBHs along the orbit. In fact affects the computation of at all orders, starting from the simple quadrupole term, by ”splitting” each polarization amplitude and into harmonics according to (see, e.g., equations (5-6) in  and references therein):
The amplitude coefficients , , and are linear combinations of the Bessel functions of the first kind , and , and is an additional precession term to the phase given by . For , and we recover the circular limit given by equation (21). As GW emission tends to decrease eccentricity, this is likely to be mostly important at large separations (i.e., , see panels ’a3’ and ’b3’ in figure 2). An example of an eccentric waveform is shown in the upper panel of figure 4.
Turning now to spins, equation (19) shows that spins enter as higher order corrections in the computation of the acceleration and, consequently, of the waveform. Therefore they become important only when is substantial, significantly affecting the signal only for . Such corrections generate additional small terms in the phase evolution, but most importantly, cause the orbital angular momentum to precess around the total conserved angular momentum (in which and the spins of the two MBHs) according to a precession equation, that to the leading order reads 
(in this case is the binary separation). The typical precession timescale is given by
where we used , and we assumed equal mass binaries in the last equality. It is clear that at, say, , is of the order of several tens of years, making precession effects negligible. An example of waveform modulated by plane precession in the late inspiral is shown in the lower panel of figure 4.
The most general detectable signal from a spinning eccentric binary is a function of 17 parameters (some describing the intrinsic properties of the binary and some others related to the relative binary-detector position and orientation): two combination of the redshifted masses, and 444As discussed in Section 2, for sources at cosmological distances, the GW depends on the redshifted masses, the intrinsic ones can be extracted by measuring and then measuring according to a cosmological model, or by obtaining an independent measurement of through, e.g., the identification of an electromagnetic counterpart to the GW signal.; six parameters defining the individual spin vectors, two magnitudes and and four angles; two parameters related to the eccentricity of the orbit, initial eccentricity and an additional angle defining the line of nodes; source inclination with respect to the line of sight, ; polarization angle, ; sky location, two angle usually labeled and 555Of particular interest is the source sky localization errorbox , defined, following , in terms of and as: (according to the notation used in Section 4).; luminosity distance ; initial orbital phase ; and, depending on the type of signal, initial frequency or time to coalescence (the two are related, in the quadrupole approximation, by equation (8)). We saw examples of how some of these parameters are imprinted in the waveform, in the following subsection we turn to the problem of how accurately they can be extracted given some signal observation.
4.3 Parameter extraction
We briefly review the basic theory regarding the estimate of the statistical errors that affect the measurements of the source parameters. For a comprehensive discussion of this topic we refer the reader to . The data collected in a detector is given by the superposition of the noise and a signal determined by a vector of parameters :
In the following, we make the usual (simplifying) assumption that is a zero-mean Gaussian and stationary random process characterized by the one-sided power spectral density . The inference process in which we are interested is how well one can infer the actual value of the unknown parameter vector , based on the data , and any prior information on available before the experiment. Within the Bayesian framework, one is therefore interested in deriving the posterior probability density function (PDF) of the unknown parameter vector given the data set and the prior information. Bayes’ theorem yields
where is the likelihood function, is the prior probability density of , and is the marginal likelihood or evidence. Under the assumption of Gaussian noise
where the inner product between two functions is defined as the integral
applied to the Fourier Transform of the functions, e.g.,
In the neighborhood of the maximum-likelihood estimate value , the likelihood function can be approximated as a multi-variate Gaussian distribution,
where and the matrix is the Fisher information matrix (FIM). Here label the components of (i.e., the parameters defining the shape of the signal), and we have used Einstein’s summation convention (and we do not distinguish between covariant and contravariant indices). The FIM is simply related to the derivatives of the GW signal with respect to the unknown parameters integrated over the observation:
In terms of the inner product the maximal signal-to-noise ratio (SNR) at which a signal can be observed is obtained by matched filtering of the data against a template equal to the waveform signal, . The optimal matched filtering SNR achievable in this way is . In the limit of large SNR, tends to , and the inverse of the FIM provides a lower limit to the error covariance of unbiased estimators of , the so-called Cramer-Rao bound. The variance-covariance matrix is simply the inverse of the FIM, and its elements are
where () are the correlation coefficients. We can therefore interpret as a way to quantify the expected uncertainties on the measurements of the source parameters. We refer the reader to  and references therein for an in-depth discussion of the interpretation of the inverse of the FIM in the context of assessing the prospect of the estimation of the source parameters for GW observations. When combining different pieces of independent (i.e., having uncorrelated noise) information (for example, by observing several pulsars, or by combining the inspiral and the ringdown portion of a signal), the FIM that characterizes the joint observations in equation (30) is simply given by the sum of the matrices of all the individual pieces
and all the theory follows unchanged.
5 The gravitational wave landscape: observations and scientific payouts
As shown in Section 2, GW emission from MBHBs covers several decades in frequency, ranging from sub-nano-Hz to milli-Hz. As shown in figure 5, this range is (or it will be) covered by multiple probes. The ground-based network of advanced interferometric detectors (three LIGO detectors, VIRGO , and the Kamioka Gravitational wave Detector, KAGRA ) and possibly the third-generation Einstein Telescope (ET, ) will observe inspiralling binaries up to around few. The milli-Hz regime will be the hunting territory of spaced based detectors such as eLISA, whereas PTAs are already probing the nano-Hz portion of the frequency band. In particular, space based interferometers and PTAs will provide a complementary, complete census of the MBHB population throughout the Universe. In this section we focus on the prospects of GW observation in these two bands and on the related scientific payouts.
5.1 The milli-Hz regime: science with space based interferometry
Space based interferometry will open a revolutionary new window on the Universe. In the following we refer to the eLISA design presented in  to describe the extraordinary scientific payouts of milli-Hz GW observations. Figure 6 highlights the exquisite capabilities of eLISA in covering almost all the mass-redshift parameter space relevant to MBH astrophysics. GW observations will catch sources with at early cosmological times, prior to reionization. A binary with can be detected out to with a SNR , making an extensive census of the MBH population in the Universe possible.
As detailed in Section 4, detected waveforms carry information on all the relevant source parameters, including redshifted masses and spins of the individual BHs prior to coalescence, the distance to the source and its sky location. The left panel of figure 7 shows error distributions in the source parameter estimation, for events collected in a meta-catalog of sources, based on state of the art MBH evolution models (see  for details). Here circular precessing spinning binary were considered (i.e. waveforms determined by 15 parameters), and ”hybrid” waveforms of the PhenomC family were used to evaluate uncertainties based on the FIM approximation, as outlined in the previous section. Individual redshifted masses can be measured with unprecedented precision, i.e. with an error of , on both components. The spin of the primary hole can be measured with an exquisite accuracy, to a 0.01-0.1 absolute uncertainty. This precision mirrors the big imprint left by the primary MBH spin in the waveform. The measurement is more problematic for that can be either determined to an accuracy of 0.1, or remain completely undetermined, depending on the source mass ratio and spin amplitude. The source luminosity distance error has a wide spread, usually ranging from being undetermined (but see  for possible shortcomings of the FIM approximation in these cases) to a stunning few percent accuracy (note that this is a direct measurement of ) GW detectors are full sky monitors, and the localisation of the source in the sky is also encoded in the waveform pattern. Sky location accuracy is typically estimated in the range 10-1000 square degrees666These numbers assume a ’single Michelson’ (four laser links) configuration for eLISA, the full triangular ’two Michelsons’ (six laser links) configuration results in a significant improvement of the estimation of all parameters, in particular luminosity distance and sky location (by 1-2 orders of magnitude)..
While measurements of individual systems are extremely interesting and very useful for making, e.g., strong-field tests of GR, it is the properties of the whole set of MBHB mergers that are observed which will carry the most information for astrophysics. GW observations of multiple MBHBs may be used together to learn about their formation and evolution through cosmic history, as demonstrated by [72, 73]. We briefly provide here an illustrative example from . As argued above, in the general picture of MBH cosmic evolution, the population is shaped by the seeding process and the accretion history.  therefore consider a set of 4 models with distinctive properties: (i) small seeds and extended (coherent) accretion (SE), (ii) light seeds and chaotic accretion (SC); (iii) large seeds and extended accretion (LE), (iv) large seeds and chaotic accretion (LC). Each model predicts a theoretical distribution of coalescing MBHBs. A given dataset of observed events can be compared to a given model by computing the likelihood that the observed dataset is a realization of model . When testing a dataset against a pair of models and , one assigns probability to model , and probability to model . The probabilities and are a measure of the relative confidence one has in model and , given an observation . Setting a confidence threshold of one can count what fraction of the 1000 realizations of model yield a confidence when compared to an alternative model . Results are shown in the left-hand panel of table 1 for all pairs of models, assuming one year observation and circular non-spinning waveforms (i.e., for an extremely conservative waveform model). The vast majority of the pair comparisons yield a confidence in the true model for almost all the realizations, with the exception of comparisons LE to LC and SE to SC, i.e., comparisons among models differing by accretion mode only. This is because the accretion mode (efficient versus chaotic) particularly affects the spin distribution of the coalescing systems, which is not considered in the circular non-spinning waveform model. It is sufficient to add a measurement of the remnant spin parameter to make those pairs easily distinguishable (Right-hand panel of table 1). The right panel of figure 7 shows the evolution of the fraction of correctly identified models as a function of observation time (no spin information included). Small versus large seed scenarios (SE vs LE and SE vs LC) can be easily discriminated after only 1 year of observation.
5.2 The nano-Hz regime: science with pulsar timing arrays
PTAs are sensitive at much lower frequencies (Hz), where the expected signal is given by a superposition of a large number of massive (), relatively nearby () sources overlapping in frequency. As argued in Section 3, at such low frequencies the properties of the MBHBs are likely to be severely affected by their coupling with their stellar and gaseous environment. In particular binaries can be highly eccentric, which might suppress the low frequency portion of the spectrum, crucial to PTA detection . Here we consider circular GW driven binaries for simplicity. The overall expected characteristic strain of the GW signal can be written as 
where , is the comoving number of binaries emitting in a given logarithmic frequency interval with chirp mass and redshift in the range and , respectively; and is the inclination–polarization averaged strain given by equation (4).
The GW spectrum has a characteristic shape , where is the signal normalization at , which depends on the details of the MBH binary population only. A Montecarlo realization of the signal is shown in figure 8 for a selected MBHB population model. In the timing residual of each individual pulsar, the signal appears as a structured red noise (left panel), but a representation of the characteristic strain in the Fourier domain reveals the complexity of its nature (right panel). Although several millions of sources contribute to it, the bulk of the strain comes from few hundred sources only. Therefore, the signal is far from being a Gaussian isotropic background ; a handful of sources dominates the strain budget, and some of them might be individually identified. Detection techniques have been developed for stochastic signals [81, 82, 83, 84] and individual sources [85, 86, 87], and more sophisticated schemes accounting for signal anisotropy have recently been proposed [88, 89, 90]. In terms of level of the stochastic signal, recent works [80, 78, 76] set a plausible range , the upper limit being already in tension with current PTA measurement [77, 91]. This is shown in the left panel of figure 9, where observations are compared to theoretically predictions. Here, the difference between the top-left and the top-right panel is given by the recent upgrades in the MBH mass-host relation [92, 93] to include the overmassive black holes measured in brightest cluster galaxies (BCGs) , that boosts the range of expected signal by a factor of two. In the lower panels instead, we consider two subset of the models featuring these upgraded relations: (i) those in which accretion does not occur prior to binary coalescence, and (ii) those in which accretion precedes the formation of the binary, and is more prominent on the secondary MBH . In the latter case, binaries observed by PTA are much more massive (and with a larger ) implying a much larger (by almost a factor of three) signal.
An handful of sources might be bright enough to be individually resolved, and in this case some of their parameters can be determined according to the scheme described in Section 4. A pioneering investigation was performed by  assuming circular, non spinning monochromatic systems. In this case the waveform is function of 7 parameters only: the amplitude 777For monochromatic signals the two masses and the luminosity distance degenerate into a single amplitude parameter., sky location , polarization , inclination , frequency and phase , defining the parameter vector . Results about typical parameter estimation accuracy are shown in the right panel of figure 9. For SNR The source amplitude is determined to a accuracy, whereas are only determined within a fraction of a radian. For bright enough sources (SNR) sky location within few tens to few deg is possible (see also [85, 87]), and even sub deg determination, under some specific conditions . Even though this is a large chunk of the sky, these systems are extremely massive and at relatively low redshift (), making any putative electromagnetic signature of their presence (e.g., emission periodicity related to the binary orbital period, peculiar emission spectra, peculiar K line profiles, etc.) detectable [76, 97].
We provided a general overview of massive black hole binaries as gravitational wave sources. MBHs are today ubiquitous in massive galaxies, and power luminous quasars up to . Although they are believed to play a central role in the process of structure formation, their origin and early growth is largely unknown. According to our current understanding, MBHBs must form in large numbers along the cosmic history, providing the loudest sources of GWs in the Universe in a wide range of frequencies spanning from the sub-nano-Hz up to the milli-Hz. GWs carry precise information about the parameters of the emitting systems. We showed how those parameters are imprinted in the phase (and amplitude) modulation of the wave, and can therefore be efficiently extracted and determined to high accuracy with ongoing and future GW probes. From those we will learn about MBH formation and evolution through cosmic history, about the nature of the first BH seeds, their subsequent accretion history, and, more generally, about the early hierarchical structure formation at high redshift. We will also learn about the complex interplay of physical processes, including stellar and gas dynamics and GW emission, that leads to the dynamical formation and evolution of MBHBs. Direct GW detection will open a new era in MBH and MBHB astrophysics.
I would like to thank S. Babak, M. Dotti, F. Ohme and A. Petiteau for providing useful material. This work is supported by the DFG grant SFB/TR 7 Gravitational Wave Astronomy and by DLR (Deutsches Zentrum fur Luft- und Raumfahrt).
-  J. Magorrian, S. Tremaine, D. Richstone, R. Bender, G. Bower, A. Dressler, S. M. Faber, K. Gebhardt, R. Green, C. Grillmair, J. Kormendy, and T. Lauer. The Demography of Massive Dark Objects in Galaxy Centers. The Astronomical Journal, 115:2285–2305, June 1998.
-  D. J. Mortlock, S. J. Warren, B. P. Venemans, M. Patel, P. C. Hewett, R. G. McMahon, C. Simpson, T. Theuns, E. A. Gonzáles-Solares, A. Adamson, S. Dye, N. C. Hambly, P. Hirst, M. J. Irwin, E. Kuiper, A. Lawrence, and H. J. A. Röttgering. A luminous quasar at a redshift of z = 7.085. Nature, 474:616–619, June 2011.
-  S. D. M. White and M. J. Rees. Core condensation in heavy halos - A two-stage theory for galaxy formation and clustering. Mon. Not. Roy. Astron. Soc. , 183:341–358, May 1978.
-  G. Kauffmann and M. Haehnelt. A unified model for the evolution of galaxies and quasars. Mon. Not. Roy. Astron. Soc. , 311:576–588, January 2000.
-  M. Volonteri, F. Haardt, and P. Madau. The Assembly and Merging History of Supermassive Black Holes in Hierarchical Models of Galaxy Formation. Astrophys. J., 582:559–573, January 2003.
-  K. S. Thorne. Gravitational Waves. In E. W. Kolb and R. D. Peccei, editors, Particle and Nuclear Astrophysics and Cosmology in the Next Millenium, page 160, 1995.
-  S. A. Hughes. Untangling the merger history of massive black holes with LISA. Mon. Not. Roy. Astron. Soc. , 331:805–816, April 2002.
-  P. Amaro-Seoane et al. eLISA: Astrophysics and cosmology in the millihertz regime. GW Notes, Vol. 6, p. 4-110, 6:4–110, May 2013.
-  P. Amaro-Seoane et al. Low-frequency gravitational-wave science with eLISA/NGO. Classical and Quantum Gravity, 29(12):124016, June 2012.
-  T. e. Consortium. The Gravitational Universe. ArXiv e-prints, May 2013.
-  R. D. Ferdman et al. The European Pulsar Timing Array: current efforts and a LEAP toward the future. Classical and Quantum Gravity, 27(8):084014, April 2010.
-  R. N. Manchester et al. The Parkes Pulsar Timing Array Project. ArXiv e-prints 1210.6130, October 2012.
-  F. A. Jenet et al. The North American Nanohertz Observatory for Gravitational Waves. ArXiv e-prints 0909.1058, September 2009.
-  G. Hobbs et al. The International Pulsar Timing Array project: using pulsars as a gravitational wave detector. Classical and Quantum Gravity, 27(8):084013–+, April 2010.
-  S. A. Hughes. Listening to the universe with gravitational-wave astronomy. Annals of Physics, 303:142–178, January 2003.
-  K. S. Thorne. Gravitational radiation., pages 330–458. 1987.
-  M. Tegmark, J. Silk, M. J. Rees, A. Blanchard, T. Abel, and F. Palla. How Small Were the First Cosmological Objects? Astrophys. J., 474:1, January 1997.
-  V. Bromm, P. S. Coppi, and R. B. Larson. Forming the First Stars in the Universe: The Fragmentation of Primordial Gas. Astrophys. J. Lett., 527:L5–L8, December 1999.
-  P. Madau and M. J. Rees. Massive Black Holes as Population III Remnants. Astrophys. J. Lett., 551:L27–L30, April 2001.
-  M. Volonteri. Formation of supermassive black holes. The Astronomy and Astrophysics Review, 18:279–315, July 2010.
-  M. C. Begelman, M. Volonteri, and M. J. Rees. Formation of supermassive black holes by direct collapse in pre-galactic haloes. Mon. Not. Roy. Astron. Soc. , 370:289–298, July 2006.
-  B. Devecchi and M. Volonteri. Formation of the First Nuclear Clusters and Massive Black Holes at High Redshift. Astrophys. J., 694:302–313, March 2009.
-  L. Mayer, S. Kazantzidis, A. Escala, and S. Callegari. Direct formation of supermassive black holes via multi-scale gas inflows in galaxy mergers. Nature, 466:1082–1084, August 2010.
-  A. Sesana, M. Volonteri, and F. Haardt. The imprint of massive black hole formation models on the LISA data stream. Mon. Not. Roy. Astron. Soc. , 377:1711–1716, June 2007.
-  S. A. Hughes and R. D. Blandford. Black Hole Mass and Spin Coevolution by Mergers. Astrophys. J. Lett., 585:L101–L104, March 2003.
-  N. I. Shakura and R. A. Sunyaev. Black holes in binary systems. Observational appearance. Astronomy & Astrophysics, 24:337–355, 1973.
-  K. S. Thorne. Disk-Accretion onto a Black Hole. II. Evolution of the Hole. Astrophys. J., 191:507–520, July 1974.
-  A. R. King, S. H. Lubow, G. I. Ogilvie, and J. E. Pringle. Aligning spinning black holes and accretion discs. Mon. Not. Roy. Astron. Soc. , 363:49–56, October 2005.
-  S. A. Kassin, B. J. Weiner, S. M. Faber, J. P. Gardner, C. N. A. Willmer, A. L. Coil, M. C. Cooper, J. Devriendt, A. A. Dutton, P. Guhathakurta, D. C. Koo, A. J. Metevier, K. G. Noeske, and J. R. Primack. The Epoch of Disk Settling: z ~ 1 to Now. Astrophys. J., 758:106, October 2012.
-  M. H. Fabricius, R. P. Saglia, D. B. Fisher, N. Drory, R. Bender, and U. Hopp. Kinematic Signatures of Bulges Correlate with Bulge Morphologies and Sérsic Index. Astrophys. J., 754:67, July 2012.
-  M. Dotti, M. Colpi, S. Pallini, A. Perego, and M. Volonteri. On the Orientation and Magnitude of the Black Hole Spin in Galactic Nuclei. Astrophys. J., 762:68, January 2013.
-  C. S. Reynolds. Measuring Black Hole Spin using X-ray Reflection Spectroscopy. ArXiv e-prints, February 2013.
-  R. D. Blandford and R. L. Znajek. Electromagnetic extraction of energy from Kerr black holes. Mon. Not. Roy. Astron. Soc. , 179:433–456, May 1977.
-  J. C. Mihos and L. Hernquist. Gasdynamics and Starbursts in Major Mergers. Astrophys. J., 464:641, June 1996.
-  A. Sesana. Self Consistent Model for the Evolution of Eccentric Massive Black Hole Binaries in Stellar Environments: Implications for Gravitational Wave Observations. Astrophys. J., 719:851–864, August 2010.
-  G. D. Quinlan. The dynamical evolution of massive black hole binaries I. Hardening in a fixed stellar background. New Astronomy, 1:35–56, July 1996.
-  D. Merritt, S. Mikkola, and A. Szell. Long-Term Evolution of Massive Black Hole Binaries. III. Binary Evolution in Collisional Nuclei. Astrophys. J., 671:53–72, December 2007.
-  T. Matsubayashi, J. Makino, and T. Ebisuzaki. Orbital Evolution of an IMBH in the Galactic Nucleus with a Massive Central Black Hole. Astrophys. J., 656:879–896, February 2007.
-  M. Preto, I. Berentzen, P. Berczik, and R. Spurzem. Fast Coalescence of Massive Black Hole Binaries from Mergers of Galactic Nuclei: Implications for Low-frequency Gravitational-wave Astrophysics. Astrophys. J. Lett., 732:L26, May 2011.
-  A. Sesana, F. Haardt, and P. Madau. Interaction of Massive Black Hole Binaries with Their Stellar Environment. I. Ejection of Hypervelocity Stars. Astrophys. J., 651:392–400, November 2006.
-  P. B. Ivanov, J. C. B. Papaloizou, and A. G. Polnarev. The evolution of a supermassive binary caused by an accretion disc. Mon. Not. Roy. Astron. Soc. , 307:79–90, July 1999.
-  Z. Haiman, B. Kocsis, and K. Menou. The Population of Viscosity- and Gravitational Wave-driven Supermassive Black Hole Binaries Among Luminous Active Galactic Nuclei. Astrophys. J., 700:1952–1969, August 2009.
-  P. J. Armitage and P. Natarajan. Eccentricity of Supermassive Black Hole Binaries Coalescing from Gas-rich Mergers. Astrophys. J., 634:921–927, December 2005.
-  J. Cuadra, P. J. Armitage, R. D. Alexander, and M. C. Begelman. Massive black hole binary mergers within subparsec scale gas discs. Mon. Not. Roy. Astron. Soc. , 393:1423–1432, March 2009.
-  C. Roedig, M. Dotti, A. Sesana, J. Cuadra, and M. Colpi. Limiting eccentricity of subparsec massive black hole binaries surrounded by self-gravitating gas discs. Mon. Not. Roy. Astron. Soc. , 415:3033–3041, August 2011.
-  J. M. Bardeen and J. A. Petterson. The Lense-Thirring Effect and Accretion Disks around Kerr Black Holes. Astrophys. J. Lett., 195:L65, January 1975.
-  P. A. G. Scheuer and R. Feiler. The realignment of a black hole misaligned with its accretion disc. Mon. Not. Roy. Astron. Soc. , 282:291, September 1996.
-  A. Perego, M. Dotti, M. Colpi, and M. Volonteri. Mass and spin co-evolution during the alignment of a black hole in a warped accretion disc. Mon. Not. Roy. Astron. Soc. , 399:2249–2263, November 2009.
-  F. Ohme. Analytical meets numerical relativity: status of complete gravitational waveform models for binary black holes. Classical and Quantum Gravity, 29(12):124002, June 2012.
-  É. É. Flanagan and S. A. Hughes. Measuring gravitational waves from binary black hole coalescences. I. Signal to noise for inspiral, merger, and ringdown. Phys. Rev. D , 57:4535–4565, April 1998.
-  L. Blanchet. Gravitational Radiation from Post-Newtonian Sources and Inspiralling Compact Binaries. Living Reviews in Relativity, 9:4, June 2006.
-  A. Buonanno and T. Damour. Effective one-body approach to general relativistic two-body dynamics. Phys. Rev. D , 59(8):084006, April 1999.
-  A. Buonanno, Y. Chen, and M. Vallisneri. Detecting gravitational waves from precessing binaries of spinning compact objects: Adiabatic limit. Phys. Rev. D , 67(10):104025, May 2003.
-  F. Pretorius. Evolution of Binary Black-Hole Spacetimes. Physical Review Letters, 95(12):121101, September 2005.
-  M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlochower. Accurate Evolutions of Orbiting Black-Hole Binaries without Excision. Physical Review Letters, 96(11):111101, March 2006.
-  J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and J. van Meter. Gravitational-Wave Extraction from an Inspiraling Configuration of Merging Black Holes. Physical Review Letters, 96(11):111102, March 2006.
-  E. Berti, V. Cardoso, and C. M. Will. Gravitational-wave spectroscopy of massive black holes with the space interferometer LISA. Phys. Rev. D , 73(6):064030, March 2006.
-  L. Santamaría, F. Ohme, P. Ajith, B. Brügmann, N. Dorband, M. Hannam, S. Husa, P. Mösta, D. Pollney, C. Reisswig, E. L. Robinson, J. Seiler, and B. Krishnan. Matching post-Newtonian and numerical relativity waveforms: Systematic errors and a new phenomenological model for nonprecessing black hole binaries. Phys. Rev. D , 82(6):064016, September 2010.
-  T. Damour, A. Nagar, M. Hannam, S. Husa, and B. Brügmann. Accurate effective-one-body waveforms of inspiralling and coalescing black-hole binaries. Phys. Rev. D , 78(4):044039, August 2008.
-  A. Buonanno, Y. Pan, H. P. Pfeiffer, M. A. Scheel, L. T. Buchman, and L. E. Kidder. Effective-one-body waveforms calibrated to numerical relativity simulations: Coalescence of nonspinning, equal-mass black holes. Phys. Rev. D , 79(12):124028, June 2009.
-  L. E. Kidder. Coalescing binary systems of compact objects to (post)-Newtonian order. V. Spin effects. Phys. Rev. D , 52:821–847, July 1995.
-  B. Willems, A. Vecchio, and V. Kalogera. Probing White Dwarf Interiors with LISA: Periastron Precession in Eccentric Double White Dwarfs. Physical Review Letters, 100(4):041102, February 2008.
-  A. Vecchio. LISA observations of rapidly spinning massive black hole binary systems. Phys. Rev. D, 70(4):042001, August 2004.
-  C. Cutler. Angular resolution of the LISA gravitational wave detector. Phys. Rev. D , 57:7089–7102, June 1998.
-  E. T. Jaynes and G. L. Bretthorst. Probability Theory. April 2003.
-  M. Vallisneri. Use and abuse of the Fisher information matrix in the assessment of gravitational-wave parameter-estimation prospects. Phys. Rev. D , 77(4):042001, February 2008.
-  J. Aasi, J. Abadie, B. P. Abbott, R. Abbott, T. D. Abbott, M. Abernathy, T. Accadia, F. Acernese, C. Adams, T. Adams, and et al. Search for gravitational waves from binary black hole inspiral, merger, and ringdown in LIGO-Virgo data from 2009-2010. Phys. Rev. D , 87(2):022002, January 2013.
-  K. Somiya. Detector configuration of KAGRA-the Japanese cryogenic gravitational-wave detector. Classical and Quantum Gravity, 29(12):124007, June 2012.
-  M. Punturo et al. The Einstein Telescope: a third-generation gravitational wave observatory. Classical and Quantum Gravity, 27(19):194002, October 2010.
-  M. Volonteri and P. Natarajan. Journey to the M- relation: the fate of low-mass black holes in the Universe. Mon. Not. Roy. Astron. Soc. , 400:1911–1918, December 2009.
-  A. Sesana. Detecting Massive Black Hole Binaries and Unveiling their Cosmic History with Gravitational Wave Observations. In G. Auger, P. Binétruy, and E. Plagnol, editors, 9th LISA Symposium, volume 467 of Astronomical Society of the Pacific Conference Series, page 103, January 2013.
-  J. E. Plowman, R. W. Hellings, and S. Tsuruta. Constraining the black hole mass spectrum with gravitational wave observations - II. Direct comparison of detailed models. Mon. Not. Roy. Astron. Soc. , 415:333–352, July 2011.
-  A. Sesana, J. Gair, E. Berti, and M. Volonteri. Reconstructing the massive black hole cosmic history through gravitational waves. Phys. Rev. D , 83(4):044036, February 2011.
-  A. Sesana. Insights on the astrophysics of supermassive black hole binaries from pulsar timing observations. ArXiv e-prints, July 2013.
-  A. Sesana, A. Vecchio, and C. N. Colacino. The stochastic gravitational-wave background from massive black hole binary systems: implications for observations with Pulsar Timing Arrays. Mon. Not. R. Astron. Soc., 390:192–209, October 2008.
-  A. Sesana, C. Roedig, M. T. Reynolds, and M. Dotti. Multimessenger astronomy with pulsar timing and X-ray observations of massive black hole binaries. Mon. Not. R. Astron. Soc., 420:860–877, February 2012.
-  R. van Haasteren et al. Placing limits on the stochastic gravitational-wave background using European Pulsar Timing Array data. Mon. Not. R. Astron. Soc., 414:3117–3128, July 2011.
-  S. T. McWilliams, J. P. Ostriker, and F. Pretorius. Gravitational waves and stalled satellites from massive galaxy mergers at z 1. ArXiv e-prints, November 2012.
-  A. Sesana and A. Vecchio. Measuring the parameters of massive black hole binary systems with pulsar timing array observations of gravitational waves. Phys. Rev. D, 81(10):104008–+, May 2010.
-  V. Ravi, J. S. B. Wyithe, G. Hobbs, R. M. Shannon, R. N. Manchester, D. R. B. Yardley, and M. J. Keith. Does a ”Stochastic” Background of Gravitational Waves Exist in the Pulsar Timing Band? Astrophys. J., 761:84, December 2012.
-  R. W. Hellings and G. S. Downs. Upper limits on the isotropic gravitational radiation background from pulsar timing analysis. Astrophys. J. Lett., 265:L39–L42, February 1983.
-  F. A. Jenet, G. B. Hobbs, K. J. Lee, and R. N. Manchester. Detecting the Stochastic Gravitational Wave Background Using Pulsar Timing. Astrophys. J. Letter, 625:L123–L126, June 2005.
-  R. van Haasteren, Y. Levin, P. McDonald, and T. Lu. On measuring the gravitational-wave background using Pulsar Timing Arrays. Mon. Not. R. Astron. Soc., 395:1005–1014, May 2009.
-  M. Anholm, S. Ballmer, J. D. E. Creighton, L. R. Price, and X. Siemens. Optimal strategies for gravitational wave stochastic background searches in pulsar timing data. Phys. Rev. D, 79(8):084030–+, April 2009.
-  J. A. Ellis, X. Siemens, and J. D. E. Creighton. Optimal Strategies for Continuous Gravitational Wave Detection in Pulsar Timing Arrays. Astrophys. J., 756:175, September 2012.
-  S. Babak and A. Sesana. Resolving multiple supermassive black hole binaries with pulsar timing arrays. Phys. Rev. D , 85(4):044034, February 2012.
-  A. Petiteau, S. Babak, A. Sesana, and M. de Araújo. Resolving multiple supermassive black hole binaries with pulsar timing arrays. II. Genetic algorithm implementation. Phys. Rev. D , 87(6):064036, March 2013.
-  N. J. Cornish and A. Sesana. Pulsar Timing Array Analysis for Black Hole Backgrounds. ArXiv e-prints, May 2013.
-  C. M. F. Mingarelli, T. Sidery, I. Mandel, and A. Vecchio. Characterising gravitational wave stochastic background anisotropy with Pulsar Timing Arrays. ArXiv e-prints, June 2013.
-  S. R. Taylor and J. R. Gair. Searching For Anisotropic Gravitational-wave Backgrounds Using Pulsar Timing Arrays. ArXiv e-prints, June 2013.
-  P. B. Demorest, R. D. Ferdman, M. E. Gonzalez, D. Nice, S. Ransom, I. H. Stairs, Z. Arzoumanian, A. Brazier, S. Burke-Spolaor, S. J. Chamberlin, J. M. Cordes, J. Ellis, L. S. Finn, P. Freire, S. Giampanis, F. Jenet, V. M. Kaspi, J. Lazio, A. N. Lommen, M. McLaughlin, N. Palliyaguru, D. Perrodin, R. M. Shannon, X. Siemens, D. Stinebring, J. Swiggum, and W. W. Zhu. Limits on the Stochastic Gravitational Wave Background from the North American Nanohertz Observatory for Gravitational Waves. Astrophys. J., 762:94, January 2013.
-  N. J. McConnell and C.-P. Ma. Revisiting the Scaling Relations of Black Hole Masses and Host Galaxy Properties. Astrophys. J., 764:184, February 2013.
-  A. W. Graham and N. Scott. The M -L Relation at High and Low Masses, the Quadratic Growth of Black Holes, and Intermediate-mass Black Hole Candidates. Astrophys. J., 764:151, February 2013.
-  J. Hlavacek-Larrondo, A. C. Fabian, A. C. Edge, and M. T. Hogan. On the hunt for ultramassive black holes in brightest cluster galaxies. Mon. Not. R. Astron. Soc., 424:224–231, July 2012.
-  S. Callegari, S. Kazantzidis, L. Mayer, M. Colpi, J. M. Bellovary, T. Quinn, and J. Wadsley. Growing Massive Black Hole Pairs in Minor Mergers of Disk Galaxies. Astrophys. J., 729:85, March 2011.
-  K. J. Lee, N. Wex, M. Kramer, B. W. Stappers, C. G. Bassa, G. H. Janssen, R. Karuppusamy, and R. Smits. Gravitational wave astronomy of single sources with a pulsar timing array. Mon. Not. R. Astron. Soc., 414:3251–3264, July 2011.
-  T. Tanaka, K. Menou, and Z. Haiman. Electromagnetic counterparts of supermassive black hole binaries resolved by pulsar timing arrays. ArXiv e-prints, July 2011.