Astrophysical constraints on massive black hole binary evolution from pulsar timing arrays

Abstract

We consider the information that can be derived about massive black hole binary (MBHB) populations and their formation history solely from current and possible future pulsar timing array (PTA) results. We use models of the stochastic gravitational-wave background from circular MBHBs with chirp mass in the range 10⁶–10¹¹ M_⊙ evolving solely due to radiation reaction. Our parametrized models for the black hole merger history make only weak assumptions about the properties of the black holes merging over cosmic time. We show that current PTA results place an upper limit on the black hole merger density which does not depend on the choice of a particular merger history model; however, they provide no information about the redshift or mass distribution. We show that even in the case of a detection resulting from a factor of 10 increase in amplitude sensitivity, PTAs will only put weak constraints on the source merger density as a function of mass, and will not provide any additional information on the redshift distribution. Without additional assumptions or information from other observations, a detection cannot meaningfully bound the massive black hole merger rate above zero for any particular mass.

black hole physics, gravitational waves, methods: data analysis, pulsars: general, Galaxy: evolution

INTRODUCTION

Massive black holes (MBHs) reside at the centre of most galaxies (see e.g. Kormendy & Ho 2013, and references therein), and are believed to have a central role in their evolution (see e.g. Volonteri 2012, and references therein for a recent review). Mapping the population of MBHs, studying their properties, demographics, and their connection to the broader formation of structure is one of the open problems of modern astrophysics. This is however difficult to tackle, due to the large range of scales and the wide variety of physical processes involved. The MBH evolutionary path remains a highly debated subject with many competing hypotheses still in play. Currently, favoured hierarchical structure formation scenarios imply frequent galaxy mergers (White & Rees 1978). As a result, MBH binaries (MBHBs) should be quite common in the Universe (Begelman, Blandford & Rees 1980; Volonteri, Haardt & Madau 2003). To date there is no confirmed observed MBHB, although a number of candidates exist (see e.g. Dotti, Sesana & Decarli 2012, and references therein) and tantalizing claims have been recently made (Graham et al. 2015a,b; Liu et al. 2015).

A means to survey MBHBs is through observation of gravitational waves (GWs) that these systems generate as they inspiral towards their final merger. The accurate timing of an array of highly stable millisecond pulsars – a pulsar timing array (PTA; Foster & Backer 1990) – provides a direct observational means to probe the cosmic population of MBHBs on orbital time-scales of order of several years. Astrophysical modelling suggests that the radiation emitted by an ensemble of MBHBs produces a GW stochastic background (GWSB) in the frequency range ∼10⁻⁹–10⁻⁷ Hz, where PTAs operate (Sesana, Vecchio & Colacino 2008; Sesana, Vecchio & Volonteri 2009; Ravi et al. 2012; Sesana 2013b). Such a background affects the time of arrival of radio pulses in a characteristic fashion (Sazhin 1978; Detweiler 1979; Hellings & Downs 1983), which can be used to discriminate the signal from a plethora of other undesired effects (Lentati et al. 2015).

Over the last decade pulsar timing has placed progressively tighter constraints on gravitational radiation in this frequency regime (see e.g. Jenet et al. 2006). More recently, the three international consortia consisting of the Parkes Pulsar Timing Array (PPTA; Shannon et al. 2013, 2015), NANOGrav (Demorest et al. 2013; Arzoumanian et al. 2015), and the European Pulsar Timing Array (EPTA; van Haasteren et al. 2011, 2012; Lentati et al. 2015), which in collaboration form the International Pulsar Timing Array (Hobbs et al. 2010), have used data from observations of unprecedented sensitivity to place constraints that are starting to probe astrophysically interesting regions of the parameter space (Sesana 2013b).

In this Letter, we consider a GWSB produced by MBHBs in circular orbits losing energy and angular momentum purely through GW emission. We use an analytical merger rate model which makes minimal assumptions about the cosmological history of MBHB evolution and can capture the key characteristics of simulation results to investigate the astrophysical implications of current (Shannon et al. 2013; Arzoumanian et al. 2015; Lentati et al. 2015), and future plausible (Siemens et al. 2013; Ravi et al. 2015) PTA results (either an upper limit or a detection). Because our model is fully general – not committing to any particular cosmological MBHB merger history – we can identify and separate features of the merger history that are constrained by PTA data alone, from those that can only be constrained by adopting a particular merger history (as e.g. done in Shannon et al. 2013; Arzoumanian et al. 2015) – in other words, by applying a particular cosmological prior. Because our model is capable of reproducing the MBHB cosmic population found in cosmological simulations for certain choices of parameters, our results will be consistent with (but much broader than) those that would be obtained under a choice of specific classes of MBHB merger history models.

In Section 2, we describe our method and merger rate model. In Section 3, we present our results for several upper limits and a possible future dectection, and discuss implications for the MBHB population. We present our conclusions in Section 4.

MODEL AND METHOD

Astrophysical model

For the standard scenario of circular binaries driven by radiation reaction only, the characteristic strain of the GWSB, h at frequency f is given by (Phinney 2001):

\begin{eqnarray} h^2(f) &=& \frac{4G^{5/3}}{3\pi ^{1/3}c^2} f^{-4/3} \int {\rm d} \log _{10} \mathcal {M}\int {\rm d}z (1+z)^{-1/3} \mathcal {M}^{5/3} \nonumber \\ &&\times \frac{{\rm d}^3N}{{\rm d}V_{\mathrm{c}}\;{\rm d}z\;{\rm d}\log _{10}\mathcal {M}}, \end{eqnarray}

(1)

where z is the redshift and |$\mathcal {M}$| is the chirp mass related to the binary component masses (m₁, m₂) by |$\mathcal {M}=(m_1m_2)^{3/5}/(m_1+m_2)^{1/5}$|⁠. The integral sums over the sources in z and |$\mathcal {M}$| weighted by the distribution of the source population, |${\rm d}^3N/{\rm d}V_{\mathrm{c}}\;{\rm d}z\;{\rm d}\log _{10}\mathcal {M}$|⁠, the number of binary mergers per comoving volume, redshift, and (rest-frame) chirp mass interval. We choose a simple model for this, described by

\begin{eqnarray} \frac{{\rm d}^3N}{{\rm d}V_{\mathrm{c}}\;{\rm d}z\;{\rm d}\log _{10}\mathcal {M}} &=& {\dot{n}}_0\, \left[\left( \frac{\mathcal {M}}{10^7\,\mathrm{M}_{{\odot }}}\right)^{-\alpha } \exp ^{-(\mathcal {M}/{\cal M}_{*})}\right]\nonumber \\ &&\times \left[(1+z)^{\beta } \exp ^{-(z/z_0)}\right] \frac{{\rm d}t_R}{{\rm d}z}\,, \end{eqnarray}

(2)

where t_R is the time in the source rest frame (here we use H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_M = 0.3, Ω_Λ = 0.7, and Ω_k = 0). Following general astrophysical assumptions, we consider a scenario where the GW background is produced by MBHBs in the redshift and chirp mass range of 0 ≤ z ≤ 5 and |$10^6\le \mathcal {M}/\mathrm{M}_{{\odot }}\le 10^{11}$|⁠. These ranges set the integration limits of equation (1).

The model is described by five parameters. The parameter |${\dot{n}}_0$| is the normalized merger rate per unit rest-frame time, comoving volume, and logarithmic (rest-frame) chirp mass interval. The parameters β and z₀ describe the distribution of the sources in redshift. The parameter β controls the low-redshift power-law slope and the parameter z₀ the high-redshift cut-off for the distribution; the peak of the merger rate d²N/dt_R dV_c, corresponds to a redshift (z₀β − 1). The parameters α and |${\cal M}_{*}$| provide a similar description of the chirp mass distribution. The model was chosen to capture the expected qualitative features of the cosmic MBH merger rate without restricting to any particular merger history; for example, it can reproduce rates extracted from merger tree models (Volonteri et al. 2003; Sesana et al. 2008), and large-scale cosmological simulations of structure formation (Springel et al. 2005; Sesana et al. 2009).

The characteristic amplitude has a simple power-law scaling, and we can rewrite equation (1) as h(f) = A_1yr(f/f_1yr)^{− 2/3}, where A_1yr is the characteristic amplitude at the reference frequency f_1yr = 1 yr⁻¹, which is customarily used when quoting limits in the PTA literature. A single number, the amplitude A_1yr, carries the whole information about the merging history of MBHBs (within the model considered in this Letter), that one wishes to reconstruct from the observations.

Method

The objective is to put constraints on the population parameters, which we denote by θ, given the results of PTA analyses. In our case θ is a five-dimensional parameter space, |$\mathbf {\theta } = \lbrace {\dot{n}}_0, \beta , z_0, \alpha , {\cal M}_{*}\rbrace$|⁠. We want to compute the posterior density function (PDF) of the parameters given PTA observations denoted by d. The population parameters fully specify the GW signal h(f; θ) (equation 1), which in turn specifies the statistical properties of the GW-induced deviations to pulse arrival times, the PTA observable. Given data from pulsar timing and our model for the merger rate (equation 2), we use Bayes’ theorem to find the posterior distribution of the model parameters, p(θ|d) = p(θ)p(d|A_1yr(θ))/p(d), where p(d|A_1yr(θ)) is the PTA likelihood for a given stochastic background, h(f; θ), p(θ) is the prior on the model parameters and p(d) is the evidence. In standard analysis of the PTA data, constraints are put on the GW characteristic amplitude at periods of 1 yr, A_1yr, which in turn is a function of the parameters of the underlying population, specified by h(f; θ). The PTA analysis uses a likelihood function, p(d|A_1yr(θ)), which we approximate as described below. Our method does not rely on this approximation; we use it only for analytical convenience in this Letter. If a given PTA analysis provides a posterior distribution for A_1yr, then a straightforward reweighting can produce the corresponding likelihood required for our analysis (if flat priors on A_1yr are used in the analysis, then the reweighting is trivial because the posterior and the likelihood are proportional to each other).

In this Letter, we consider the two cases in which the PTA analysis provides either an upper limit or a detection. For the upper-limit scenario, we model p(d|A_1yr) using a Fermi-like distribution,

\begin{equation} p_{\mathrm{ul}}(d|A_{\mathrm{1yr}}) \propto \left\{\exp[(A_{\mathrm{1yr}}-A_{\mathrm{ul}})/ \sigma_{\mathrm{ul}}]+1\right\}^{-1}, \end{equation}

(3)

where A_ul is the upper-limit value returned by the actual analysis. The sharpness of the tail-off, σ_ul can be adjusted to give an upper limit with a chosen confidence, which we set at 95 per cent. We model a detection scenario using a Gaussian in the logarithm of A_1yr,

\begin{eqnarray} p_{\mathrm{det}}\left(d|A_{\mathrm{1yr}} \right) \propto \exp\{-[\log_{10}(A_{\mathrm{1yr}})-\log_{10}(A_{\mathrm{det}})]^2/2\sigma_{\mathrm{det}}^2\}, \end{eqnarray}

(4)

at a chosen level of detection A_det. We choose the width of the detection to be σ_det = 0.2. We compute the marginalized distribution on the model parameters θ using two independent sampling techniques, to verify the results of our analysis: a nested sampling approach (Veitch & Vecchio 2010) and emcee, an ensemble Markov Chain Monte Carlo sampler (Foreman-Mackey 2013).

Our priors on the model parameters are set as follows. We use a prior on |${\dot{n}}_0$| that is flat in |$\log _{10} {\dot{n}}_0$| down to a lower limit, which we set to |${\dot{n}}_0 = 10^{-20}\, \mathrm{Mpc}^{-3}\, \mathrm{Gyr}^{-1}$|⁠, after which it is flat in |${\dot{n}}_0$| to zero. The prior upper bound is set to 10³ Mpc⁻³ Gyr⁻¹. This value is set by the ultraconservative assumption that all the matter in the Universe is formed by MBHs. Our prior allows for the number of mergers to span many orders of magnitude (flat in log ) but avoids divergence as |${\dot{n}}_0 \rightarrow 0$|⁠. It also allows for the absence of MBHB binaries merging within an Hubble time. The priors for the other parameters are uniform within ranges that incorporate values that give a good fit to semi-analytical merger tree models (see e.g. Sesana et al. 2008, 2009; Sesana 2013b): α ∈ [−3.0, 3.0], β ∈ [−2.0, 7.0], z₀ ∈ [0.2, 5.0], and |$\log _{10}{\cal M}_{*}/\,\mathrm{M}_{\odot } \in [6.0, 9.0]$|⁠. While our prior allows for parameter values that can reproduce the merger rates of detailed models, it is uninformative in that we do not assume that the merger rate distribution must take values from those models. Our priors reflect large theoretical uncertainties about MBHB formation and evolution scenarios, and the lack of any confirmed MBHB candidate (see however Graham et al. 2015a,b; Liu et al. 2015).

Our method is summarized as follows: (i) produce a likelihood for A_1yr (or use smoothed posterior samples from PTA results reweighted if necessary), (ii) choose a MBHB merger rate model, (iii) produce PDFs for the model parameters from which we infer MBHBs population properties.

RESULTS

Current upper limits on the GWSB obtained recently are A_1yr^{(95 per cent)} = 1 × 10^{− 15}, 1.5 × 10⁻¹⁵, 3 × 10⁻¹⁵ for the PPTA (Shannon et al. 2015), NANOGrav (Arzoumanian et al. 2015), and the EPTA (Lentati et al. 2015), respectively. The sensitivity gain provided by the addition of new pulsars to the PTAs and more recent data sets may allow in the short- to mid-term to reach a sensitivity below A_1yr = 1.0 × 10⁻¹⁵, and in the more distant future A_1yr ∼ 10⁻¹⁶ (Siemens et al. 2013; Ravi et al. 2015). As a consequence, here we consider three PTA analysis outcomes: (i) an upper limit at 95 per cent confidence of 1 × 10⁻¹⁵, which represents the present state of play and either (ii) an upper limit (at 95 per cent confidence) of 1 × 10⁻¹⁶ or (iii) a detection at the same level, that is A_det = 1 × 10⁻¹⁶, and σ_det = 0.2 in equation (4), which describes possible results coming from the expected improvements of the PTA sensitivity in the next 5–10 yr.

The main results of our analysis are summarized in Fig. 1, which shows the inferred posterior distribution of the merger history of MBHBs in terms of the MBHBs comoving volume merger density per redshift and (logarithm of) chirp mass intervals. Fig. 2 provides PDFs on selected parameters based on current PTA limits, and Fig. 3 provides a similar summary for a future limit or a detection at the level described above.

Figure 1.

Posteriors for the merger rate density. The top row shows the merger rate density in chirp mass (integrated over redshift), |${\rm d}^2N/{\rm d}V_{\mathrm{c}}{\rm d}\log _{10}{\cal M}$| and the bottom row in redshift (integrated over chirp mass), d²N/dV_cdz, for two 95 per cent confidence upper limits at 1 × 10⁻¹⁵ (left) and 1 × 10⁻¹⁶ (centre) and a detection at 1 × 10⁻¹⁶ (right), as described in the text. We consider contributions to the GW background from MBHs in the chirp mass range |$10^{6}\le \mathcal {M}/ \mathrm{M}_{\odot }\le 10^{11}$| and redshift range 0 ≤ z ≤ 5. The solid black line gives the posterior median; dark grey, mid-grey, and light-grey bands show the central 68 per cent, 95 per cent, and 99 per cent credible interval, respectively. The dashed lines show draws from the posterior. For comparison, the overlaid dark areas represent the 99.7 per cent confidence regions predicted by the MBH assembly models of Sesana (2013b). For these models, we show only the chirp mass in the range ≈10⁶–10⁹ M_⊙ as outside this interval the lower percentile is zero. For the redshift range, these models only consider MBHB mergers for z ≲ 1.3.

Figure 2.

Marginalized posterior distributions for selected astrophysical parameters for the case of 95 per cent upper limit of 1 × 10⁻¹⁵, which corresponds to the current status of the observations. The marginalized PDF on the merger rate parameter, |$\dot{n}_0$| is shown in the left-hand panel, and the marginalized PDF on |$({\cal M}_{*}, \dot{n}_0)$| in the right-hand panel, where the contours mark the 67 per cent (solid) and 95 per cent (dashed) confidence regions. In the left-hand panel, the dashed lines mark the 95 per cent confidence width (⁠|$-20.8 \le \log_{10} [\dot{n}_0/\mathrm{Mpc}^{-3}\, \mathrm{Gyr}^{-1}]\le -2.3$|⁠) while the dotted line marks the 95 per cent upper limit (⁠|$\log_{10} [\dot{n}_0/\mathrm{Mpc}^{-3}\, \mathrm{Gyr}^{-1}]\le -3.3$|⁠). The left-hand side of the distribution in |$\log_{10} (\dot{n}_0/\mathrm{Mpc}^{-3}\, \mathrm{Gyr}^{-1})$| follows our prior, while the right-hand side is determined by the PTA upper limit.

Figure 3.

Posterior distribution for the upper limit (left) and detection (right) at 1 × 10⁻¹⁶. The top panels show the one-dimensional posterior distribution for the merger rate parameter, |$\dot{n}_0$|⁠. The dashed lines mark the 95 per cent confidence width (upper limit: |$-20.9 \le \log_{10} [\dot{n}_0/\mathrm{Mpc}^{-3}\, \mathrm{Gyr}^{-1}]\le -4.2$|⁠; detection: |$-11.7 \le \log_{10} [\dot{n}_0/\mathrm{Mpc}^{-3}\, \mathrm{Gyr}^{-1}]\le -1.3$|⁠) and the dotted line the 95 per cent upper limit (upper limit: |$\log_{10} [\dot{n}_0/\mathrm{Mpc}^{-3}\, \mathrm{Gyr}^{-1}] = -5.0$|⁠; detection: |$\log_{10} [\dot{n}_0/\mathrm{Mpc}^{-3}\, \mathrm{Gyr}^{-1}] = -1.9$|⁠). The central and bottom panels show the two-dimensional posterior distributions for |${\dot{n}}_0$| with the mass parameters α and |${\cal M}_{*}$|⁠. The solid and dashed contours mark the 67 per cent and 95 per cent confidence regions, respectively.

We consider first the implications of current limits. The PDFs on the parameters |$\dot{n}_0$| and |${\cal M}_{*}$| of the model are shown in Fig. 2; we do not provide the equivalent plots for α, β, and z₀ as they are equivalent to the prior. Fig. 2 clearly shows that the present PTA limits enable us to reduce the allowed normalization of the MBHB merger rate density to |$\dot{n}_0 \lesssim 5\times 10^{-3}\, \mathrm{Mpc}^{-3}\, \mathrm{Gyr}^{-1}$| with 95 per cent confidence, but yield no additional constraints on the other parameters of the model.

Our model parameters describe the shape of the merger rate distribution in redshift and chirp mass. The PDFs of those parameters induce a posterior density on |${\rm d}^2N/{\rm d}V_{\mathrm{c}}{\rm d}\log_{10}{\cal M}$| and d²N/dV_cdz, integrating over redshift and chirp mass, respectively, as shown in Fig. 1. We see that current observations limit the maximum merger density as a function of mass, but place no constraints on the shape of the distribution. The corresponding number of sources per frequency bin that contribute to the signal is |${\rm d}^2N/{\rm d}f {\rm d}\log_{10}\mathcal{M} \propto f^{-11/3}$|⁠, and we find that for masses above a few × 10⁹ M_⊙, our upper limit on |${\rm d}^2\,N / {\rm d}V_{\mathrm{c}}{\rm d}\log _{10}{\cal M}$| implies that at a frequency around 1.8 nHz there is fewer than one source per frequency bin (taken to be Δf = 1/T, with T = 17.66 yr, the timespan of current EPTA data sets; Lentati et al. 2015). This means that at those large masses, the assumption that the observed GW signal is stochastic is violated, and our analysis cannot be used to constrain the exact shape of the mass function here (in this case a different PTA search approach would be necessary, see e.g. Arzoumanian et al. 2014; Babak et al. 2015; Taylor et al. 2015). While current PTA observations provide feeble constraints on the shape of the mass distribution, they yield no information about the redshift distribution. The bottom-left panel of Fig. 1 shows no structure in d²N/dV_cdz.

It is useful to compare these results to limits on the MBHB merger rates implied by binary candidates reported in the literature and to specific theoretical models. Let us first consider what is known observationally today. A few MBHB candidates have been reported recently. Graham et al. (2015b) suggested the possible observation of a MBHB at redshift z = 0.2784 with (rest frame) total mass log (M/M_⊙) ∼ 8.5 and period of ∼1884 d. Liu et al. (2015) reported the observation of a potential MBHB at z = 2.060 with a shorter period of 542 d and primary MBHB mass log (M/M_⊙) ∼ 9.97. Using the redshift to calculate the enclosed volume and the binary parameters for the time to merger, we can estimate the predicted rate from each of these observations. Assuming that these two systems are indeed MBHBs, and that their constituents are of comparable mass, they imply merger rates of ≈3 × 10⁻⁷ and ≈0.1 Mpc⁻³ Gyr⁻¹ (the latter number takes into account that the source has been found in an analysis of only one of the Pan-STARRS1 Medium Deep Survey fields of 8 deg²). In turn they yield a merger density |${\rm d}^2N/{\rm d}V_{\mathrm{c}}{\rm d}\log _{10}{\cal M}$| ≈10⁻⁶ Mpc⁻³ at |${\cal M}\approx 3\times 10^{8}\,\mathrm{M}_{{\odot }}$| and ≈1 Mpc⁻³ at |${\cal M}\approx 10^{10}\,\mathrm{M}_{{\odot }}$|⁠, respectively. The upper-left panel of Fig. 1 clearly shows that the rate density inferred from Graham et al. (2015b) is consistent with current upper limits, while that inferred from Liu et al. (2015) is several orders of magnitude above the 99 per cent credible interval implied by current PTA results. It is therefore unlikely that this source is an MBHB with the claimed parameters. Other proposed MBHBs in the literature (Valtonen, Ciprini & Lehto 2012; Kun et al. 2014) imply merger density estimates of 5 × 10⁻⁵ Mpc⁻³ at |$\mathcal {M}\approx 10^9 \,\mathrm{M}_{\odot }$| and 3 × 10⁻⁶ Mpc⁻³ at |$\mathcal {M}\approx 3.5 \times 10^{8} \, \,\mathrm{M}_{\odot }$|⁠, which are consistent with the current upper limits.

On the theoretical side, current limits are consistent with the assumption that most Milky Way-like galaxies contain an MBH in the mass range considered here that undergoes ∼1 major merger in an Hubble time. The density of Milky Way-like galaxies is 10⁻² Mpc⁻³, which yields an estimate of |${\rm d}^2N/{\rm d}V_{\mathrm{c}}{\rm d}\log _{10}{\cal M} \sim 10^{-3}\,\mathrm{Mpc}^{-3}$|⁠, which is consistent with our results at |$\mathcal {M}\sim 10^6\,\mathrm{M}_{{\odot }}$|⁠, appropriate for a typical MBHB forming in the merger of Milky Way-like galaxies. We also compare the limits on the |${\rm d}^2N/{\rm d}V_{\mathrm{c}}{\rm d}\log_{10}{\cal M}$| and d²N/dV_cdz with specific distributions obtained from predictions of astrophysical models for the cosmic assembly of MBHs. We consider the models presented in Sesana (2013b), extended to include the most recent MBH-galaxy scaling relations (Kormendy & Ho 2013). These models produce a central 99 per cent interval of A_1yr ∈ [2 × 10⁻¹⁶, 4 × 10⁻¹⁵]. The 99.7 per cent confidence region in the merger density from those models is marked by a dark-shaded area in each panel of Fig. 1. Two conclusions can be drawn: (i) present MBHB population models are consistent with current PTA limits and (ii) those models are drawn from a very restricted prior range of the parameters that control the evolution of MBHBs, driven by specific assumptions on their assembly history. For example, in those models there is a one-to-one correspondence between galaxy and MBH mergers. Our results are consistent with the conclusions drawn by Shannon et al. (2013) about the implications of the PPTA limit for the MBHB merger history. However, since they consider specific models that lie close to the upper end of the 99 per cent credible range allowed by current limits, they emphasize the fact that PTA limits might soon be in tension with those specific classes of models.

We turn now to consider what we could infer about the MBHB merger history in the future as PTA sensitivity increases. For definiteness, we consider both an upper limit and a detection at the level of A_1yr = 10⁻¹⁶. Selected marginalized PDFs on the model parameters are shown in Fig. 3, where we see a slight correlation in the 2D marginalized PDF of |$({\cal M}_{*}, \dot{n}_0)$|⁠, as expected. This is simply explained by considering the Schechter-like mass profile of equation (2): as the characteristic mass |${\cal M}_{*}$| decreases, and therefore the exponential cut-off of MBHB progressively depletes the high-mass portion of the population, a given value of the GW characteristic amplitude allows for a larger overall normalization, |$\dot{n}_0$|⁠. The posterior MBHB merger densities per logarithm of chirp mass and redshift are shown in Fig. 1. For the case of an upper limit, the results are qualitatively similar to the case of the present PTA upper limit, simply scaled accordingly. In particular, despite the much tighter limit on the overall merger rate we are still unable to place any meaningful constraint on the redshift distribution of merging MBHBs. The overall merger density as a function of redshift shifts by two orders of magnitude and the same is true for the merger density as a function of mass. Note that a non-detection at this level might pose a serious challenge to currently favoured theoretical MBH assembly models with simple black hole dynamics, as shown in the upper-centre panel of Fig. 1.

In the case of detection, the posterior on the shapes of the merger rate distribution in redshift and chirp mass are plotted on the right-hand panels of Fig. 1. We still obtain essentially no bounds on the shape of the merger rate density in redshift. We also obtain no meaningful lower bound on the merger rate density for chirp masses. That is, there is no chirp mass at which we can bound the merger density above a rate physically indistinguishable from zero; we know that some MBHBs merge, but we cannot determine which ones. Additional information, such as theoretical assumptions, electromagnetic observations constraining the mass spectrum of merging black holes (like those discussed earlier in this section), or GW observations that measure the binary mass spectrum directly (such as those of an eLISA-like instrument Amaro-Seoane et al. 2013), is required to place any constraints on the masses of the merging systems. For example, if we accept the priors provided by Sesana (2013b), the mass function of merging MBHBs can be determined more precisely, as shown by the overlap between our posterior and the dark band in the upper–right panel of Fig. 1.

CONCLUSIONS

We have considered the implications of current PTA limits on the GWSB to constrain the merger history of MBHBs. Using a general model for the mass and redshift evolution of MBHBs in circular orbit driven by radiation reaction, we find that existing PTA results alone place essentially no constraints on the merger history of MBHBs. We also find that even with an increase in amplitude sensitivity of an order of magnitude, and assuming that a detection is made, no bounds can be put on the functional form of the merger rate density in redshift and chirp mass unless additional information coming through a different set of observations is available.

Finally, we want to caution the reader that the results presented here apply only within the model assumptions that have been made. We have considered a generic (and well justified) functional form for the MBHB merger rate density, but if one chooses a significantly different form (and associated priors for the parameters), results could be different (even radically). Moreover, it has been suggested that physical effects other than radiation reaction, such as gas and/or interactions with stars (e.g. Kocsis & Sesana 2011; Sesana 2013a; Sampson, Cornish & McWilliams 2015), could affect the evolution of MBHBs. These effects are not included in our model, and their impact on astrophysical inference needs to be evaluated in the future.