The nightmare scenario: measuring the stochastic gravitationalwave background from stalling massive blackhole binaries with pulsartiming arrays
Abstract
Massive blackhole binaries, formed when galaxies merge, are among the primary sources of gravitational waves targeted by ongoing Pulsar Timing Array (PTA) experiments and the upcoming spacebased LISA interferometer. However, their formation and merger rates are still highly uncertain. Recent upper limits on the stochastic gravitationalwave background obtained by PTAs are starting to be in marginal tension with theoretical models for the pairing and orbital evolution of these systems. This tension can be resolved by assuming that these binaries are more eccentric or interact more strongly with the environment (gas and stars) than expected, or by accounting for possible selection biases in the construction of the theoretical models. However, another (pessimistic) possibility is that these binaries do not merge at all, but stall at large ( pc) separations. We explore this extreme scenario by using a semianalytic galaxy formation model including massive black holes (isolated and in binaries), and show that future generations of PTAs will detect the stochastic gravitationalwave background from the massive blackhole binary population within years of observations, even in the “nightmare scenario” in which all binaries stall at the hardening radius. Moreover, we argue that this scenario is too pessimistic, because our model predicts the existence of a subpopulation of binaries with small mass ratios () that should merge within a Hubble time simply as a result of gravitationalwave emission. This subpopulation will be observable with large signaltonoise ratios by future PTAs thanks to nextgeneration radio telescopes such as SKA or FAST, and possibly by LISA.
keywords:
binaries, black holes, gravitational waves, galaxies: evolution1 Introduction
Massive black holes (MBHs) with masses in the range are ubiquitous in the nuclei of nearby and distant galaxies (Kormendy & Richstone, 1995). In the accepted framework of hierarchical structure formation, massive galaxies are formed by continuous accretion of dark matter and gas from cosmic filaments, and by (minor and major) galaxy mergers. The latter process is expected to lead to the formation of a population of MBH binaries in the nuclei of postmerger galaxies (Begelman, Blandford & Rees, 1980). If these binaries are at sufficiently close separations, they efficiently emit gravitational waves (GWs), which may be observable by ongoing Pulsar Timing Array (PTA) experiments (Hellings & Downs, 1983) for total binary masses of  and separations of hundreds to thousands of gravitational radii. These experiments include the European Pulsar Timing Array (EPTA; Desvignes et al. 2016), the Parkes Pulsar Timing Array (PPTA; Reardon et al. 2016) and the North American Nanohertz Observatory for Gravitational Waves (NANOGrav; The NANOGrav Collaboration et al. 2015), joining together in the International Pulsar Timing Array (IPTA; Verbiest et al. 2016). Moreover, MBH binaries with total masses  will also be observable in the late inspiral, merger and ringdown phases by the upcoming spaceborne Laser Interferometer Space Antenna (LISA; Audley et al. 2017; Klein et al. 2016). More precisely, both LISA and PTAs will not only target the GWs from individual resolved sources, but also the stochastic background resulting from the superposition of the GWs produced by all the unresolved sources that exist in the universe. Still, although LISA is expected to detect e.g. a significant stochastic background of Galactic whitedwarf binaries (Audley et al., 2017), unresolved MBH binaries are expected to be relatively rare in the LISA frequency range (i.e. LISA will detect the majority of the MBH mergers in the Universe, or even all of them depending on the astrophysical model, c.f. Klein et al. 2016). Conversely, PTAs are expected to first detect the stochastic GW background from MBH binaries, though they are sensitive also to signals from individual loud sources (Rosado, Sesana & Gair, 2015).
Indeed, PTA upper limits on the stochastic GW background from MBH binaries have steadily improved over the past few years, and they have recently started being in marginal tension with the predictions of theoretical models (Lentati et al., 2015; Arzoumanian et al., 2016; Shannon et al., 2015). This is not surprising in itself, because when two galaxies merge, the MBHs are expected to be deposited in the outskirts of the newly formed galaxy, at separations that could be as large as kpc. Early on, dynamical friction from the stellar and gas background is probably very efficient at driving the MBHs towards the galactic centre (since most MBHs will be still surrounded by a relic stellar cluster inherited from their previous host galaxy). However, when the MBHs form a bound binary, dynamical friction becomes inefficient, because the binary’s orbital velocity exceeds the typical velocity of the stars. The subsequent shrinking of the binary is then assured by threebody interactions with stars. Indeed, stars with angular momentum in an appropriate region of parameter space (the “loss cone”) will interact strongly with the binary and typically extract energy from it. As a result, the binary will shrink (Quinlan, 1996), while the stars will be scattered away and possibly ejected from the galaxy as hypervelocity stars (Sesana, Haardt & Madau, 2006). After a phase of fast shrinking, the binary starts hardening at a rate when its separation reaches the hardening radius ( pc for total masses of , c.f. Equation 9 below). This will eventually deplete the “loss cone” in the phase space of the surrounding stars, which will result in threebody interactions also becoming inefficient. While the loss cone will be replenished naturally by the scattering between stars on the relaxation timescale, this exceeds the Hubble time for galaxies hosting MBHs with masses above . Since GW emission does not become efficient enough to drive the binary to merger within a Hubble time until the separation becomes of the order of pc (for total masses of , c.f. Equation 10 below), binaries with large masses () may therefore stall at separations pc. This is known as the “lastparsec problem” (c.f. Milosavljević & Merritt, 2001; Merritt & Milosavljević, 2005; Preto et al., 2011; Colpi, 2014).
A way around this problem is provided by processes that could help replenish the loss cone, e.g. galaxy rotation (HolleyBockelmann & Khan, 2015) or a triaxiality in the galactic gravitational potential (Yu, 2002; Khan, Just & Merritt, 2011; Vasiliev, 2014; Vasiliev, Antonini & Merritt, 2014, 2015, 2014; Sesana & Khan, 2015), induced for instance by mergers. This replenishment would make stellar scattering efficient again at driving the orbital evolution down to the separation needed for the binary to merge as a result of GW emission. Other possibilities to overcome the lastparsec problem are the presence of a gaseous disc, which would result in planetarylike migration of the MBHs towards the centre (Haiman, Kocsis & Menou 2009; see however Lodato et al. 2009 for complications arising in this scenario); or the interaction with a third incoming MBH coming from a subsequent galaxy merger, which would trigger the coalescence of the binary via the combined action of KozaiLidov resonances and GW emission (Hoffman & Loeb, 2007; Bonetti et al., 2016).
Nevertheless, while there is a consensus that the lastparsec problem will be somehow solved, the exact mechanism by which this would happen is still debated. As a matter of fact, the aforementioned PTA limits on the stochastic GW background are starting to probe (and even in some cases to be in marginal tension with) models for MBH binary formation and evolution that assume that all binaries merge efficiently under the effect of GW emission alone (Lentati et al., 2015; Arzoumanian et al., 2016; Shannon et al., 2015). Several ways to explain the PTA limits have been proposed. MBH binaries are normally assumed to be almost circular when they enter the PTA band, but they could have a significant nonzero eccentricity (Taylor, Simon & Sampson, 2017), which could be left over for instance from triple MBH KozaiLidov oscillations and which would move at least part of the radiated power outside the PTA band. Binaries may also interact more strongly with the environment (gas and stars) than expected (Kocsis & Sesana, 2011; Sampson, Cornish & McWilliams, 2015; Ravi et al., 2014; Taylor et al., 2016; Kelley et al., 2017), and these interactions may still be important at the separations that are most relevant for PTAs ( pc). As a result, a binary’s orbital energy would be lost at least partly to the environment (as heating of the gas or increase in the stars’ velocities) as the system inspirals, rather than be emitted in GWs alone. Yet another possibility is that the theoretical predictions for the number of MBH binaries are off, since they are produced with models calibrated to the MBH scaling relations, which may be biasedhigh due to selection effects (Sesana et al. 2016; see also Shankar et al. 2016; Shankar, Bernardi & Sheth 2017; Barausse et al. 2017). However, the simplest and most pessimistic possibility is that the lastparsec problem may not be solved after all, and that MBHs may stall at separations or even larger.
In this paper we explore this “nightmarescenario” by using a comprehensive semianalytic galaxyformation model (Barausse, 2012), which includes MBHs (in isolation and in binaries) and their coevolution with their host galaxies. We show that while the stochastic GW background predicted within this pessimistic scenario is way outside the reach of current experiments, it will still be detectable within years of observations by future PTA experiments, thanks to nextgeneration radio telescopes with large collecting areas such as the Square Kilometre Array radiotelescope (SKA; Smits et al., 2009) or the Five hundred meter Aperture Spherical Telescope (FAST; Nan et al., 2011). Moreover, we show that even if we insist that all binaries should stall at the hardening radius , our semianalytic model predicts the existence of a nonnegligible subpopulation of binaries with small mass ratios (). These binaries would coalesce in less than a Hubble time if initially placed at a separation . This would significantly increase the expected stochastic background signal, which should be observable with very high signaltonoise ratios by SKA or FASTbased PTAs. We also show that the formation of this subpopulation of binaries is not an artifact of the simplified prescriptions used in the semianalytic model to account for the orbital evolution of merging galaxies and MBH binaries. Moreover, we show that this subpopulation may also give rise to a few events detectable by the LISA mission, in the form of intermediate massratio inspirals (IMRIs; AmaroSeoane et al., 2007; Miller, 2009; Mandel & Gair, 2009), if MBHs form from “light” seeds (), e.g. the remnants of popIII stars.
The structure of this paper is as follows. In Section 2 we derive a general expression for the stochastic GW background from a population of stalling binaries. In Section 3 we present our semianalytic galaxy formation model and show that it predicts the existence of a subpopulation of binaries that merge within a Hubble time from their hardening radius. We also outline in more detail the model for stalling MBH binaries that we utilise in this paper. The stochastic GW background from stalling and merging MBH binaries and its detection prospects are presented in Section 4. We lay out our conclusions in Section 5. All cosmological parameters are taken from Planck Collaboration et al. (2016).
2 Gravitationalwave background
The stochastic background of GWs with energy density can be characterised by the dimensionless parameter
(1) 
where is the critical mass density of the Universe, km/s/Mpc is the Hubble constant and is the frequency measured in the detector frame, which is related to the frequency in the source frame
by , where is the redshift.
Let us consider a population of binary systems with comoving number density , where each binary is
characterised by the masses of the components , the chirp mass ,
and the separation (we will assume circular orbits throughout this paper). The total energy density resulting from the emission of GWs by these sources
is (Phinney, 2001; Sesana, Vecchio & Colacino, 2008; Rosado, 2011):
(2) 
If a binary overcomes the lastparsec problem and the merger takes place at redshift on timescales much shorter than the Hubble time, the observed spectrum is related to the emitted spectrum via
(3) 
with
(4) 
If all binaries merge efficiently, Equations 2, 3 and 4 therefore give the power law
(5) 
In practice, this power law can be suppressed at low frequencies, depending on whether environmental effects (i.e. interactions with stars and gas) are still important at the separations pc that are most important for PTAs. On the other hand, at high frequencies the signal is dominated by highmass binaries ( for Hz, c.f. Sesana, Vecchio & Colacino 2008, hereafter SVC08). Since these systems are intrinsically rare, a given realisation of the Universe may contain (on average) less than one such source close enough to contribute significantly to an observed frequency bin. SVC08 showed that above Hz, this small number statistics effect causes simulated realisations of to display excess power (relative to Equation 5) in (few) frequency bins, and lower power in the remaining ones. As a result, the slope of at Hz is flatter than predicted by Equation 5, and may even become zero or change sign at high frequencies (SVC08; c.f. also black lines in Figure 1 later on).
To model binaries that do not merge in a Hubble time but stall at a separation , we adopt two complementary approaches (in practice, as we will see, the results are close). In the first, we simply use Equations 2, 3 and 4, but cut off the spectrum given by Equation 4 outside the small frequency interval . Here, , where is the orbital frequency , while the (small) frequency shift is computed by evolving the binary under GW emission, from its formation redshift (at which the separation is ) to the present time.
This approach relies on Equations 3 and 4, which assume that GW emission happens on timescales much shorter than the Hubble time. Since this is only approximately valid for stalling binaries, we also do the calculation by assuming that the binary emits at the stalling frequency at all times after formation. This assumption, while approximate in a different way (as in reality the frequency will slowly shift due to GW emission), allows us to account for the changing redshift of the universe, i.e. each stalling binary contributes to a range of detectorframe frequencies . In more detail, since for an unevolving stalling binary and thus , the detectorframe spectrum at frequency can be expressed as
(6) 
where the emitted power is given by the quadrupole formula
(7) 
while
with and (where and are the density parameters of matter and cosmological constant). Note that as a result of the Heaviside function in the definition above, if , i.e. stalling binaries only emit at redshifted frequencies . The energy density due to a population of stalling binaries is then simply given by integrating over all sources, i.e.
(8) 
3 The model
3.1 Semianalytic galaxy formation model
We follow the mergers of MBHs by a stateoftheart semianalytic galaxy formation model introduced in Barausse (2012), with later updates to improve the spin of evolution of MBHs (Sesana et al., 2014) and to include nuclear star clusters in the centre of galaxies (Antonini, Barausse & Silk, 2015a, b). This model accounts for the cosmological evolution and merger history of galaxies inside darkmatter halos, which are produced with PressSchechter algorithms calibrated to match the results of Nbody simulations (Press & Schechter, 1974; Parkinson, Cole & Helly, 2008). Galaxies form from the cooling of a “hot” unprocessed intergalactic medium shock heated to the halo’s virial temperature, or by cold accretion flows in lowmass and highredshift systems (Dekel & Birnboim, 2006; Cattaneo et al., 2006). Once cooled or accreted to the halo’s centre, the gas settles on a disclike geometry by conservation of angular momentum, and eventually undergoes star formation. Bulges form as a result of either major galactic mergers or bar instabilities, which both destroy the stellar and gaseous discs and typically trigger bursts of star formation as the gas is funneled towards the central ( kpc) region of the galaxy. Whenever star formation takes place (in the bulge or in the disc) we account for the feedback of supernova explosions, which remove gas and tend to quench star formation, preferentially in lowmass systems.
The model also accounts for the presence of MBHs, which grow from highredshift seeds
with mass of either (“light seeds”, arising e.g. from the remnants of popIII stars; Madau & Rees, 2001)
or (“heavy seeds”, resulting for instance from protogalactic disc instabilities; Volonteri, Lodato & Natarajan, 2008).
These seed black holes are assumed to only form at , with halo occupation fractions that
depend on their exact formation mechanism (c.f. Klein et al., 2016, for details). Note however that at
the predictions of our model are essentially independent of the seed model, at least in large systems, because accretion
and mergers wash out the effect of the initial conditions as time progresses. In this paper we will see
an example of this fact, already noted e.g. in Barausse (2012); Sesana et al. (2014); Antonini, Barausse & Silk (2015a, b).
The model assumes that MBHs accrete gas from a nuclear reservoir of cold gas, whose feeding correlates linearly with bulge star formation (Granato et al., 2004; Lapi et al., 2014). Since the latter takes place in our model following major mergers and disc instabilities, our MBHs undergo long periods of quiescent activity occasionally interrupted by active quasar periods. The feedback of MBH activity on the surrounding gas (“AGN feedback”) is accounted for in both phases (radiomode and quasar feedback), and quenches star formation in preferably highmass systems.
Galaxy and blackhole mergers are modelled by starting from the halo merger history.
Whenever two halos coalesce in the extended PressSchechter merger tree, we assume that the smaller one (together with the galaxy it hosts) survives as a subhalo/satellite galaxy inside the more massive host halo. We then account for the slow infall of this satellite to the centre of the host halo by dynamical friction, by using the expressions of BoylanKolchin, Ma & Quataert (2008). Note that
these expressions are calibrated against Nbody simulations and account for the dynamical friction due to both Dark Matter and baryons (c.f. discussion in Section 2.3 of BoylanKolchin, Ma & Quataert, 2008).
As already mentioned in the Introduction, once a bound MBH binary is formed, dynamical friction becomes inefficient at driving the binary’s evolution any further, but threebody interactions of the binary with stars become important. These interactions tend to transfer energy from the binary, whose orbit shrinks, to the stars, which may even get ejected from the galaxy as hypervelocity stars (Sesana, Haardt & Madau, 2006). In more detail, after an initial fast shrinking of the orbit, the binary will harden at a constant rate after reaching the hardening radius (Merritt, 2006):
(9) 
where is the mass ratio and is the stellar velocity dispersion. However, unless mechanisms such as galaxy rotation or mergerinduced triaxiality in the galactic gravitational potential help refill the loss cone, or unless other processes that tend to shrink the binary (e.g. KozaiLidov resonances due to triple MBH interactions, or gasinduced planetarylike migration) are at play, the binary may stall at separations (“lastparsec problem”).
3.2 Stalling vs merging binaries
Unlike previous studies that were conducted with the same semianalytic model, where MBHs were either assumed to merge at the same time as the host galaxies (Barausse, 2012; Sesana et al., 2014), or after a suitable “delay” time (accounting for the hardening due to threebody interactions with stars, gasdriven planetarylike migration, and interactions with a third “intruder” MBH; Antonini, Barausse & Silk, 2015a, b; Klein et al., 2016), we hereby assume that the lastparsec problem is not solved, and consider three models that bracket the range of possible stalling scenarios for MBH binaries.
In model we assume that all MBH binaries stall exactly at the separation from which GWs would drive them to coalescence in a Hubble time Gyr (assuming circular orbits):
(10) 
This is intentionally an artificial and pessimistic model, since there is nothing special, physically, about the separation . Nonetheless, it will allow us to prove an often underappreciated point, i.e. the fact that even if all MBH binaries stall, they may still produce a stochastic GW background detectable by PTAs.
One may, however, argue that model is actually the most optimistic among all the models where MBH binaries stall, as the stalling radius may be much larger than . (Clearly, the stalling radius may not be smaller than otherwise binaries would not stall, but rather coalesce in less than a Hubble time under the effect of GW emission alone). We therefore consider an even more pessimistic model , where all binaries stall at . The hardening radius comes about because in this model we are implicitly assuming that the stalling of MBH binaries is due to losscone depletion. Moreover, note that Equations (9) and (10) imply that becomes smaller than for small mass ratios (i.e. these binaries would merge in less than a Hubble and not stall; see also Fig. 1 in Sesana 2010). Therefore, in order to be on the conservative side, we take the stalling radius to be the larger between and .
To assess what happens when this last assumption is not made, we also consider a model , where we place initially all binaries at the hardening radius when two galaxies coalesce, and evolve from there under GW emission alone.
A few comments are in order here. First, in both models and , binaries essentially always emit GWs with frequency twice the orbital frequency at the stalling radius (in the source frame). The signal in model will instead be dominated by the binaries with , for which and which therefore merge efficiently.
Second, we note that the stalling radius due to losscone depletion might actually be even smaller than (by a factor ) for comparable mass ratios (Merritt & Milosavljević, 2005). (Note however that our , given by Equation 9, agrees to within 20% with the stalling radius given by Equation 12 in Merritt 2006, at all mass ratios.) We do not account for this possible effect (which would anyway tend to increase our signal) in neither model nor , again in order to be on the conservative side.
Third, we note that our semianalytic galaxyformation model, despite accounting for the dynamical friction
on the satellite halo/galaxy from the dark matter and the baryon distributions as well as for tidal effects on the satellite, still predicts that a nonnegligible
number of unequalmass galaxy mergers should take place in a Hubble time.
As already mentioned, these systems will in turn form MBH binaries with small mass ratios
, which indeed constitute % (%) of all the binaries with total mass in the lightseed (heavyseed) case
However, the stellar bulge of the satellite galaxy only starts getting tidally disrupted when its tidal radius becomes comparable to its halflight radius . Since the tidal radius is related to the distance between the satellite and the centre of the host halo by (Henriques & Thomas, 2010)
(11) 
we obtain that the satellite’s bulge starts getting tidally disrupted at a separation
(12) 
From that separation onwards, the satellite evolution is driven by the dynamical friction of the “naked” satellite MBH against the stellar background of the host. The time needed for the satellite MBH to fall to the centre is therefore (Binney & Tremaine, 1987)
(13) 
where the subscripts “s” and “h” denote the satellite and the host, respectively. As in McWilliams, Ostriker & Pretorius (2014), we use the results of Oser et al. (2012); Nipoti et al. (2009) to assume
(14) 
Note that the redshift dependence is valid for (Oser et al., 2012). We then use (for both the host and the satellite) the correlation between blackhole and bulge stellar mass of Kormendy & Ho (2013), lowering the normalisation by a factor –3 to account for the selection bias (on the resolvability of the MBH sphere of influence) pointed out in Shankar et al. (2016), to obtain the final result
(15) 
Choosing and , becomes comparable to or larger than the lookback time only for for (the uncorrected relation of Kormendy & Ho 2013), or for for .
However, if we take into account that the selection bias highlighted by Shankar et al. (2016) may not only change the normalisation but might also make the blackhole – stellar mass relation steeper, the dynamical friction time may be even longer. If we adopt the intrinsic scaling relation of Equation 6 of Shankar et al. (2016), we find
(16) 
Still, for and , even this expression gives a dynamical friction time lower than the lookback time already at . To be on the conservative side, when considering the model where we place all binaries at the hardening radius (model ), we use Equation 16 to discard all systems for which is longer than the lookback time (c.f. discussion in Section 4). Note that since the redshift dependence is valid for , we actually replace in Equation 16, in order to avoid artificially short dynamical friction times at high redshift.
Overall, this discussion shows that it makes sense to assume that MBH binaries with efficiently reach the hardening radius after their host galaxies merge. Nonetheless, even in the absence of such systems, i.e. if all MBH binaries were to stall and not coalesce, we would fall back onto our “most pessimistic” model . We will show in the next section that even this model would still be detectable by future PTAs.
4 Results
We now address the prospects of using future PTA experiments to detect the GW stochastic backgrounds from the models discussed above, namely model , in which ; model , in which ; and model , where all binaries are assumed to form at a separation , and are let evolve under GW emission alone from there (i.e. most of the binaries will not merge by , unless they have low mass ratios , c.f. discussion above). These models are represented in Figure 1 by respectively purple, red and green bands, in the lightseed (left panel) and heavyseed (right panel) model for MBHs.
To compute the GW background for models and , we use two different approximations, as explained in Section 2. In the first, we use Equations 2, 3 and 4, but we cut off the spectrum given by Equation 4 outside the frequency interval . As explained previously, this method accounts for the orbital evolution of the binary, which sweeps a finite (albeit small) interval in sourceframe frequency from its formation to the present time, but neglects the change in cosmological redshift during the lifetime of source. In the second approximation, we use Equations 6–8, which account for the varying cosmological redshift during the lifetime of source, but neglect the orbital evolution of the binary, which is assumed to stall at a fixed separation (i.e. emit at fixed GW frequency in the source frame) after formation. The difference between these two approximations, which can be thought of as an uncertainty in our predictions, is illustrated by the width of the purple and red bands. Note that since the evolution under the influence of GW emission is not significant at the separations considered in these two models, the two approximations yield very similar results.
As for model , to be conservative we neglect the contribution of the binaries for which (which do not merge in a Hubble time and therefore give a negligible contribution from the signal from this model), and account only for those with . For this subset of binaries, we compute the background by using Equations 2, 3 and 4, but we only consider the spectrum given by Equation 4 between the initial frequency of the binary, corresponding to the hardening radius, and the final frequency that it has at (be that finite, if the binary has not yet merged by , or formally infinite, if the binary has merged by ). We have checked that neglecting these cutoffs does not change our results significantly. Note that in this model, unlike in models and , neglecting the change in cosmological redshift during the evolution of the binary is a very good approximation, since the bulk of the signal comes from binaries that merge. Note also that the finite width of the green band in Figure 1 accounts for the effect of including or not including binaries for which is longer than the lookback time at formation (c.f. discussion at the end of Section 3). As can be seen, excluding those systems only makes a small difference.
As far as the spectral shapes of the predictions in Figure 1 are concerned, observe that models and have somewhat similar behaviour (with energy density decreasing with frequency), though the normalisation is different because of the different stalling radii adopted. As for model , the signal is dominated by the subset of merging binaries, hence the spectral dependence is similar to that of a scenario in which all binaries merge efficiently (i.e. successfully overcome the lastparsec problem), shown by the blue line (and given analytically by Equation 5).
We also consider the possibility of a highfrequency turnover/flattening of the GW background, as a result of the small number of sources that may contribute at high frequencies. To this purpose, we follow SVC08 and generate several MonteCarlo realisations of the signal as it would be detected, in each frequency bin of width , by a PTA experiment of duration .
On the one hand, we find that in models and , all realisations of the signal are essentially identical to the predictions obtained by using the expressions in Section 2 (which neglect this finitestatistics effect). This is because unlike in the case discussed by SVC08, in these models there are always many sources contributing to the highfrequency bins of the spectrum. The reason is twofold. First, in models and the bulk of the signal comes from lowermass binaries than in the scenario considered in SVC08, which assumes that all binaries merge efficiently. Since the MBH mass function decreases with the MBH mass, lowermass binary systems are more numerous. Second, in the case of binaries merging efficiently, a given MBH binary sweeps the entire frequency range of PTAs very quickly, while for a stalling binary the change in frequency (in the detector frame) is much slower and to be ascribed almost entirely to the change in cosmological redshift over the system’s lifetime. As a result, simply by using the continuity equation as in SVC08, the expected number of binaries per frequency bin is much lower for merging binaries than for stalling ones.
On the other hand, for model the highfrequency part of the signal’s spectrum shows features qualitatively similar to those found in SVC08, with a general flattening or even turnover of the power law, and few pronounced “spikes” in few highfrequency bins (several realisations of the signal for an experiment duration yr are shown in black in Figure 1; note that our bins have width , therefore the lowest plotted frequency is the midpoint of the first bin, ). The resemblance to the results of SVC08 is not surprising, because in model the bulk of the signal comes from merging binaries, like in the case of SVC08. (Note however that the number of merging sources in model is lower than in SVC08, which is reflected by the different normalisation of the green bands and black lines with respect to the blue line.)
Current upper limits from ongoing PTA observations are indicated by black stars in Figure 1. We show results from PPTA (P15, Shannon et al. 2015; and P13, Shannon et al. 2013), EPTA (E15; Lentati et al., 2015) and NANOGrav (N16; Arzoumanian et al., 2016). Note that the hypothesis of efficient, circular mergers (blue line) is already excluded in our model by the PPTA limits. However, all the other scenarios (models , and ) considered in this paper are still below the observed upper limits, though they may be tested with more sensitive experiments.
In order to estimate the detection prospects with future PTAs, we consider an SKAlike experiment monitoring pulsars with ns timing accuracy for yr, and calculate a powerlaw integrated sensitivity curve (by using the procedure in Thrane & Romano, 2013, thick black line). By construction, any powerlaw spectrum tangent to this sensitivity curve has a signaltonoise ratio (SNR) of , and a powerlaw spectrum that crosses it would have . Therefore, as can be seen from Figure 1, models and would be easily detectable by such an SKAbased PTA, and even the most pessimistic scenario (model ) may be marginally detectable, as we discuss in detail below.
We note also that the results from the heavy and lightseed models are very similar to one another in the range of frequencies relevant for present and future PTAs. Indeed, the main difference between the two models is the modest bump at nHz that can be seen in the heavyseed model. This feature is due to binaries with chirp mass , i.e. stalling binaries of MBHs that have not evolved much from their initial seed masses. (“Seed binaries” are also present in the lightseed model, but because of their lower masses they radiate at higher frequencies). In the following we show only the results for the lightseed model.
The powerlaw integrated sensitivity curve shown in Figure 1 is computed in the weaksignal limit, i.e. under the assumption that the GW background is subdominant with respect to the intrinsic whitenoise component (Chamberlin et al., 2015). The whitenoise power spectrum is where is the pulsar timing accuracy and is the cadence of the pulsar measurement. For an SKAbased PTAs, reasonable typical values may be ns and yr, so the background signal produced by our models would be in the intermediate regime: It dominates the noise at low frequencies, but is subdominant at high frequencies. In this regime, we therefore have to use the general expression for the SNR (Anholm et al., 2009; Siemens et al., 2013; Chamberlin et al., 2015):
(17) 
where is the total observation time, is the Hellings and Downs coefficient for pulsars and (Hellings & Downs, 1983) and is the power spectrum of the signal. Note that the lower limit of the integral is set by the total observation time .
In Figure 2, we show the SNR by assuming an SKAbased PTA with ns and yr, as a function of number of pulsars (which we distribute isotropically in the sky) and observation time. The upper and middle panels correspond to models and , respectively. In both cases, the SNR is very sensitive to the total observation time because of the steep frequency dependence of the signal (as seen in Figure 1). As a result of this frequency dependence, only a small frequency range around contributes to the integral in Equation 17. We find that the signal can be detected with an SKAbased PTA experiment in both cases: for binaries stalling at (model ) or (model ) an SNR of () can be obtained with () pulsars, and or years of observations respectively for models and .
Detection prospects are even better for our model (bottom panel in Figure 2): can be obtained after only years of monitoring pulsars.
Current PTAs (PPTA, EPTA and NANOGrav) have worse timing accuracies than those assumed above, but have already been gathering data for several years and have built long timing baselines. Assuming a timing accuracy of only ns, current experiments will need to monitor pulsars for a duration of years in order to detect the signal from model C with an SNR of . Taking into account the data already gathered by the different PTAs, this detection might thus be not too far in the future. Detecting the signal from stalling binaries will be more challenging: in the case of models A and B, a detection with will require monitoring pulsars (all with timing accuracy of ns) for a duration of and years, respectively.
Let us now explore the range of masses of the MBH binaries that contribute to the background signal. In Figure 3, we show the contributions to the energy density for different ranges of the binary’s chirp mass. The vertical line denotes the limiting frequency corresponding to years of observation. For binaries that stall at (model ) or (model ), which are represented respectively in the upper and middle panels, the signal at is dominated by systems in the chirpmass range , with a smaller contribution from the range. This is exactly the range of masses targeted by LISA. In other words, our results suggest that these MBH binaries will be observed either by LISA if they merge efficiently (i.e. within a Hubble time; c.f. Klein et al. 2016), or by SKAbased PTAs if they stall. In the case of model (bottom panel), the mass distribution is instead quite different, with the signal being dominated by systems in the range at all frequencies. This is very similar to the masses of the binaries that contribute to the PTA GW background under the hypothesis that the finalparsec problem is efficiently solved (blue lines in Figure 1), c.f. SVC08. Clearly, this resemblance comes about because, as already mentioned, the signal in model is dominated by a subpopulation of binaries for which the finalparsec problem is not relevant, because they have , and thus coalesce in less than a Hubble time under the effect of GW emission alone.
In Figure 4, we also show the GW energy density distribution at yr in different chirpmass and redshift ranges. Note that in our model (bottom panel) the distribution is sharply peaked at , whereas in the case of stalling binaries (i.e. models and ; upper and middle panels) the distribution is significantly broader.
Finally, we have investigated whether the subpopulation of merging binaries with predicted by model contains systems that would be observable by LISA with significant detection rates. Indeed, binaries with may in principle be detectable by LISA as IMRIs. We find that the lightseed model predicts that about 1 such event would be detectable every 2 years, if we adopt the latest LISA sensitivity curve described in Audley et al. (2017). (By comparison, the LISA mission will last at least 4 years, with a possible extension to up to 10 years.) The detectable events typically have total (sourceframe) masses between a few and a few , mass ratios between a few and a few , redshift distribution peaked around and extending up to , and typical SNR . Therefore, they are formed by a MBH with mass corresponding the LISA frequency range, and a second black hole of mass , which has not accreted much during its previous history and whose mass is therefore close to the seed mass . Conversely, in our heavyseed scenario the seed mass is , so systems with mass ratio and including a seed black hole would have a total mass too large to emit a strong signal in the LISA band. In fact, we have verified that in the heavyseed scenario we only obtain IMRIs per year that are detectable by LISA.
5 Summary
Recent upper limits from PTA experiments (and especially PPTA, Shannon et al. 2015) are starting to be in tension with some of the current theoretical estimates for the merger rate of MBH binaries. This finding can be interpreted as due to the influence of the environment (gas and/or stars) on MBH binaries while they are in the PTA band (Kocsis & Sesana, 2011; Ravi et al., 2014; Taylor et al., 2016; Sampson, Cornish & McWilliams, 2015; Kelley et al., 2017); to a possible residual eccentricity (Taylor et al., 2016), resulting for instance from triple MBH interactions; or to a wrong normalisation of the theoretical predictions due to selection biases in the observations against which they are calibrated (Sesana et al., 2016). However, a much more worrisome possibility (not only for PTA experiments but also for LISA) is that MBH binaries may be unable to evolve past the hardening radius pc (lastparsec problem) and stall there. In this paper, we have used a stateoftheart semianalytic galaxyformation model including MBHs (isolated and in binaries) to study the stochastic GW background from populations of stalling binaries and its detection prospects with future PTAs. We have presented two major findings:

Even in the least favorable scenarios, the GW background produced by stalling MBH binaries might be observable with the next generation of PTAs (see also Taylor et al., 2016). Specifically, if MBH binaries stall at the separation from which they would still need a Hubble time to merge, the resulting background can be observed with a SNR of () with an SKAlike experiment that monitors () pulsars after years of observations. Even in the most pessimistic case, in which the binaries stall at , the same SNR of () can be achieved with () pulsars after years (Figure 2). This signal is dominated by binaries in the chirp mass range . Therefore, according to our results, binaries in this mass range will be detected either by LISA if they merge, or by an SKAbased PTA experiment if they stall. Observations with the timing accuracies of current PTAs will require a timing baseline of years with pulsars.

Our model predicts the existence of a subpopulation of MBH binaries with low mass ratios and hardening radii sufficiently small to allow these binaries to merge within a Hubble time under the effect of GW emission alone (). This subpopulation of binaries produces a strong GW background signal that will be easily observable by the next generation of PTAs, requiring only years of observations with pulsars at SKA sensitivity to obtain an SNR of . The timing accuracies achievable with current PTAs will likely require about years of observing time for the same number of pulsars, but in view of the data already gathered by these experiments, the time to detection might actually be shorter. We have also shown that the formation of these binaries is not an artifact of the simplified prescriptions used in the semianalytic model to account for the orbital evolution of merging galaxies and MBH binaries. Moreover, this subpopulation may yield a few detectable events for the LISA mission, if MBHs form from the remnants of popIII stars at high redshifts.
Acknowledgements
We thank Alberto Sesana, Joe Silk and Marta Volonteri for many invaluable insights and discussions about the issues presented in this paper. This work has been financially supported by the Programme National Hautes Energies (PNHE) funded by CNRS/INSUIN2P3, CEA and CNES, France. The work of ID has been done within the Labex ILP (reference ANR10LABX63) part of the Idex SUPER, and received financial state aid managed by the Agence Nationale de la Recherche, as part of the programme Investissements d’avenir under the reference ANR11IDEX000402. EB acknowledges support from the H2020MSCARISE2015 Grant No. StronGrHEP690904 and from the APACHE grant (ANR16CE310001) of the French Agence Nationale de la Recherche. This work has made use of the Horizon Cluster, hosted by the Institut d’Astrophysique de Paris. We thank Stephane Rouberol for running smoothly this cluster for us.
Footnotes
 pagerange: The nightmare scenario: measuring the stochastic gravitationalwave background from stalling massive blackhole binaries with pulsartiming arrays–References
 pubyear: 2017
 Note that another quantity widely used to characterise a GW stochastic background is the characteristic strain , which is related to by .
 Note however that the merger rate for MBHs with mass between and – i.e. the ones that are targeted by the LISA GW detector – depend more strongly on the seed model, c.f. Klein et al. (2016).
 A subtle point in the implementation of the dynamical friction timescale in a DarkMatter merger tree is given by the treatment of the coalescence of halos (each of which will in general contain its own collection of subhalos). As in Barausse (2012), if the coalescence has (DarkMatter) mass ratio , we reinitialize the dynamical friction times of all the subhalos. Otherwise, we reinitialize the dynamical friction times of the subhalos of the “satellite halo”, but keep those of the subhalos of the “host” unchanged. This corresponds to an intuitive scenario where the incoming halo perturbs and randomizes the orbits of the host’s subhalos, provided that it has a sufficiently large mass. We refer to Barausse (2012) for a more exhaustive description of this and other details of the implementation.
 Note (M. Volonteri, private communication) that MBH binaries with in this mass range were also found in a different semianalytic model, based on Volonteri, Haardt & Madau (2003).
References
 AmaroSeoane P., Gair J. R., Freitag M., Miller M. C., Mandel I., Cutler C. J., Babak S., 2007, Classical and Quantum Gravity, 24, R113
 Anholm M., Ballmer S., Creighton J. D. E., Price L. R., Siemens X., 2009, \prd, 79, 084030
 Antonini F., Barausse E., Silk J., 2015a, \apj, 812, 72
 Antonini F., Barausse E., Silk J., 2015b, \apjl, 806, L8
 Arzoumanian Z. et al., 2016, \apj, 821, 13
 Audley H. et al., 2017, ArXiv:1702.00786
 Barausse E., 2012, \mnras, 423, 2533
 Barausse E., Shankar F., Bernardi M., Dubois Y., Sheth R. K., 2017, \mnras, 468, 4782
 Begelman M. C., Blandford R. D., Rees M. J., 1980, \nat, 287, 307
 Binney J., Tremaine S., 1987, Galactic dynamics
 Bonetti M., Haardt F., Sesana A., Barausse E., 2016, \mnras, 461, 4419
 BoylanKolchin M., Ma C.P., Quataert E., 2008, \mnras, 383, 93
 Cattaneo A., Dekel A., Devriendt J., Guiderdoni B., Blaizot J., 2006, \mnras, 370, 1651
 Chamberlin S. J., Creighton J. D. E., Siemens X., Demorest P., Ellis J., Price L. R., Romano J. D., 2015, \prd, 91, 044048
 Colpi M., 2014, \ssr, 183, 189
 Dekel A., Birnboim Y., 2006, \mnras, 368, 2
 Desvignes G. et al., 2016, \mnras, 458, 3341
 Granato G. L., De Zotti G., Silva L., Bressan A., Danese L., 2004, \apj, 600, 580
 Haiman Z., Kocsis B., Menou K., 2009, \apj, 700, 1952
 Hellings R. W., Downs G. S., 1983, \apjl, 265, L39
 Henriques B. M. B., Thomas P. A., 2010, \mnras, 403, 768
 Hoffman L., Loeb A., 2007, \mnras, 377, 957
 HolleyBockelmann K., Khan F. M., 2015, \apj, 810, 139
 Kelley L. Z., Blecha L., Hernquist L., Sesana A., 2017, arXiv:1702.02180
 Khan F. M., Just A., Merritt D., 2011, \apj, 732, 89
 Klein A. et al., 2016, \prd, 93, 024003
 Kocsis B., Sesana A., 2011, \mnras, 411, 1467
 Kormendy J., Ho L. C., 2013, \araa, 51, 511
 Kormendy J., Richstone D., 1995, \araa, 33, 581
 Lapi A., Raimundo S., Aversa R., Cai Z.Y., Negrello M., Celotti A., De Zotti G., Danese L., 2014, \apj, 782, 69
 Lentati L. et al., 2015, \mnras, 453, 2576
 Lodato G., Nayakshin S., King A. R., Pringle J. E., 2009, \mnras, 398, 1392
 Madau P., Rees M. J., 2001, \apjl, 551, L27
 Mandel I., Gair J. R., 2009, Classical and Quantum Gravity, 26, 094036
 McWilliams S. T., Ostriker J. P., Pretorius F., 2014, \apj, 789, 156
 Merritt D., 2006, \apj, 648, 976
 Merritt D., Milosavljević M., 2005, Living Reviews in Relativity, 8
 Miller M. C., 2009, Classical and Quantum Gravity, 26, 094031
 Milosavljević M., Merritt D., 2001, \apj, 563, 34
 Nan R. et al., 2011, International Journal of Modern Physics D, 20, 989
 Nipoti C., Treu T., Auger M. W., Bolton A. S., 2009, \apjl, 706, L86
 Oser L., Naab T., Ostriker J. P., Johansson P. H., 2012, \apj, 744, 63
 Parkinson H., Cole S., Helly J., 2008, \mnras, 383, 557
 Phinney E. S., 2001, ArXiv:astroph/0108028
 Planck Collaboration et al., 2016, \aap, 594, A13
 Press W. H., Schechter P., 1974, \apj, 187, 425
 Preto M., Berentzen I., Berczik P., Spurzem R., 2011, \apjl, 732, L26
 Quinlan G. D., 1996, \na, 1, 35
 Ravi V., Wyithe J. S. B., Shannon R. M., Hobbs G., Manchester R. N., 2014, \mnras, 442, 56
 Reardon D. J. et al., 2016, \mnras, 455, 1751
 Rosado P. A., 2011, \prd, 84, 084004
 Rosado P. A., Sesana A., Gair J., 2015, \mnras, 451, 2417
 Sampson L., Cornish N. J., McWilliams S. T., 2015, \prd, 91, 084055
 Sesana A., 2010, \apj, 719, 851
 Sesana A., Barausse E., Dotti M., Rossi E. M., 2014, \apj, 794, 104
 Sesana A., Haardt F., Madau P., 2006, \apj, 651, 392
 Sesana A., Khan F. M., 2015, \mnras, 454, L66
 Sesana A., Shankar F., Bernardi M., Sheth R. K., 2016, \mnras, 463, L6
 Sesana A., Vecchio A., Colacino C. N., 2008, \mnras, 390, 192
 Shankar F., Bernardi M., Sheth R. K., 2017, \mnras
 Shankar F. et al., 2016, \mnras, 460, 3119
 Shannon R. M. et al., 2013, Science, 342, 334
 Shannon R. M. et al., 2015, Science, 349, 1522
 Siemens X., Ellis J., Jenet F., Romano J. D., 2013, Classical and Quantum Gravity, 30, 224015
 Smits R., Kramer M., Stappers B., Lorimer D. R., Cordes J., Faulkner A., 2009, \aap, 493, 1161
 Taffoni G., Mayer L., Colpi M., Governato F., 2003, \mnras, 341, 434
 Taylor S. R., Simon J., Sampson L., 2017, Physical Review Letters, 118, 181102
 Taylor S. R., Vallisneri M., Ellis J. A., Mingarelli C. M. F., Lazio T. J. W., van Haasteren R., 2016, \apjl, 819, L6
 The NANOGrav Collaboration et al., 2015, \apj, 813, 65
 Thrane E., Romano J. D., 2013, \prd, 88, 124032
 Vasiliev E., 2014, Classical and Quantum Gravity, 31, 244002
 Vasiliev E., Antonini F., Merritt D., 2014, Astrophys. J., 785, 163
 Vasiliev E., Antonini F., Merritt D., 2015, \apj, 810, 49
 Verbiest J. P. W. et al., 2016, \mnras, 458, 1267
 Volonteri M., Haardt F., Madau P., 2003, \apj, 582, 559
 Volonteri M., Lodato G., Natarajan P., 2008, \mnras, 383, 1079
 Yu Q., 2002, \mnras, 331, 935