Post-Newtonian Evolution of Massive Black Hole Triplets in Galactic Nuclei: IV. Implications for LISA

TL;DR

This paper models the evolution of massive black hole triplets in galactic nuclei to predict their detection rates and characteristics for LISA, highlighting the importance of triple interactions and eccentricities in gravitational wave signals.

Contribution

It introduces a semi-analytic model including a novel triple interaction channel, providing new insights into the coalescence rates and eccentricities of MBHBs relevant for LISA.

Findings

LISA may detect 25-75 MBHB events per year.

Triple interactions could account for up to 30% of detections.

High eccentricities (up to 0.9) are possible at LISA entry, affecting waveform templates.

Abstract

Coalescing massive black hole binaries (MBHBs) of $10^{4-7} \rm M_{\odot}$, forming in the aftermath of galaxy mergers, are primary targets of the space mission LISA, the {\it Laser Interferometer Space Antenna}. An assessment of LISA detection prospects requires an estimate of the abundance and properties of MBHBs that form and evolve during the assembly of cosmic structures. To this aim, we employ a semi-analytic model to follow the co-evolution of MBHBs within their host galaxies. We identify three major evolutionary channels driving the binaries to coalescence: two standard paths along which the binary evolution is driven by interactions with the stellar and/or gaseous environment, and a novel channel where MBHB coalescence occurs during the interaction with a third black hole. For each channel, we follow the orbital evolution of MBHBs with physically motivated models that include a self-consistent treatment of the orbital eccentricity. We find that LISA will detect between $\approx 25$ and $\approx 75$ events per year depending on the seed model. We show that triple-induced coalescences can range from a few detected events up to $\sim 30\%$ of the total detected mergers. Moreover, even if the standard gas/stars-driven evolutionary channels should fail and MBHBs were to stall, triple interactions would still occur as a result of the hierarchical nature of galaxy formation, resulting in about $\approx 10$ to $\approx 20$ LISA detections per year. Remarkably, triple interactions among the black holes can produce coalescing binaries with large eccentricities ($\gtrsim 0.9$) upon entrance into the LISA band. This eccentricity will remain significant ($\sim 0.1$) also at merger, requiring suitable templates for parameter estimation.

Figures

Tables

Table 1. Table 1: Merger rate and population composition of the MBHBs with source-frame chirp mass in the range 10 2 − 10 8 ​ M ⊙ superscript 10 2 superscript 10 8 subscript M direct-product 10^{2}-10^{8}~{}~{}\rm M_{\large\odot} . Coalescences are grouped according to the evolutionary channel followed, i.e. gas/stellar-driven systems (binaries) or triplet-driven systems (triplets/quadruplets). The numbers in parenthesis refer to the fraction of mergers involving a quadruple system, which for counting purposes, are considered as a sub-class of triplets. Finally, the Triplet-ej entry denotes binaries that started their evolution as part of a triplet, but at the time of coalescence, the third perturbing object has been ejected (see Bonetti et al., 2018a , for a detailed discussion) . In the following, given the rather small impact of this evolutionary channel on the mass range relevant for LISA, we neglect this sub-sample of sources.

LS
$ℳ$ [$M_{\odot }$]	Model-delayed					Model-stalled
	Rate [yr-1]	Binaries		Triplets		Rate [yr-1]	Binaries		Triplets
		gas	star	Triplet (quad)	Triplet-ej		gas	star	Triplet (quad)	Triplet-ej
${10}^2-{10}^3$	161.8	33.6%	65.0%	1.4%(0.03%)	0.01%	3.856	–	–	99.2%(3.3%)	0.7%
${10}^3-{10}^4$	49.00	54.6%	43.4%	2.0%(0.03%)	0.1%	10.27	–	–	97.2%(11.4%)	2.8%
${10}^4-{10}^5$	8.559	47.3%	48.0%	4.6%(0.05%)	0.1%	5.637	–	–	95.2%(27.9%)	4.8%
${10}^5-{10}^6$	3.430	25.9%	68.4%	5.5%(0.7%)	0.2%	1.436	–	–	95.5%(54.8%)	4.5%
${10}^6-{10}^7$	0.468	15.3%	77.6%	6.4%(0.5%)	0.7%	0.167	–	–	93.7%(81.5%)	6.3%
${10}^7-{10}^8$	0.141	10.5%	80.6%	8.3%(1.0%)	0.7%	0.041	–	–	92.6%(83.2%)	7.4%
HS
$ℳ$ [$M_{\odot }$]	Model-delayed					Model-stalled
	Rate [yr-1]	Binaries		Triplets		Rate [yr-1]	Binaries		Triplets
		gas	star	Triplet (quad)	Triplet-ej		gas	star	Triplet (quad)	Triplet-ej
${10}^2-{10}^3$	0.000	0.0 %	0.0 %	0.0 %(0.0%)	0.0%	0.000	–	–	0.0 %(0.0%)	0.0 %
${10}^3-{10}^4$	0.911	14.3%	57.7%	27.7%(0.0%)	0.0%	0.440	–	–	100.0%(0.0%)	0.0 %
${10}^4-{10}^5$	8.233	19.1%	34.3%	45.6%(1.5%)	1.0%	3.975	–	–	97.6 %(3.8%)	2.1 %
${10}^5-{10}^6$	13.11	20.9%	39.8%	37.4%(7.4%)	1.9%	6.166	–	–	96.6 %(22.7%)	3.4 %
${10}^6-{10}^7$	0.610	23.2%	67.0%	9.0 %(2.2%)	0.8%	0.301	–	–	96.3 %(57.4%)	3.7 %
${10}^7-{10}^8$	0.100	10.7%	78.1%	10.2%(1.3%)	1.0%	0.056	–	–	87.3 %(61.5%)	12.7%

Table 2. Table 2: Number of detections ( N det subscript 𝑁 det N_{\rm det} ) from the 3 different evolutionary paths considered. Binaries that merge during the LISA lifetime (left columns) represent the vast majority of events, though we find a few persistent sources that do not coalesce within the 4-year mission (right columns). For stellar and gas-driven MBHBs, comparisons between γ = 1 𝛾 1 \gamma=1 and γ = 1.75 𝛾 1.75 \gamma=1.75 and between DMM and HKM prescriptions (see Section  3.1 ) are also reported.

Model						$N_{\det }\left(>4\mathrm{yr}\right)$
	star		gas		triple	star		gas		triple
	$\gamma =1$	$\gamma =1.75$	DMM	HKM		$\gamma =1$	$\gamma =1.75$	DMM	HKM
LS-stalled	-	-	-	-	67.1	-	-	-	-	2.2
LS-delayed	135.3	160.5	154.1	153.2	6.8	1.8	2.0	0.5	0.4	0.6
HS-stalled	-	-	-	-	44.7	-	-	-	-	3.8
HS-delayed	36.7	36.9	18.3	17.8	37.3	0.6	0.4	0.2	0.1	3.3

Table 3. Table 3: Number of detections (in 4 years) of triplet-driven MBHBs with e LISA > 0.1 subscript 𝑒 LISA 0.1 e_{\rm LISA}>0.1 ( N det ​ ( e > 0.1 ) subscript 𝑁 det 𝑒 0.1 N_{\rm det}(e>0.1) ) compared to those belonging to the triplet channel only ( N det , triple subscript 𝑁 det triple N_{\rm det,triple} ), and to the total number of detections ( N det , total subscript 𝑁 det total N_{\rm det,total} ), for all four considered models.

Model	$N_{\det }\left(e_{\mathrm{LISA}}>0.1\right)$	$N_{\det ,\mathrm{triple}}$	$N_{\det ,\mathrm{total}}$
LS-stalled	31.9	69.3	69.3
LS-delayed	4.2	7.4	297.7
HS-stalled	21.8	48.5	48.5
HS-delayed	19.3	37.3	96.4

Equations

$$\frac{{\mathrm{d}}^{3}n}{\mathrm{d}z\mathrm{d}{\log }_{10}{m}_1\mathrm{d}{\log }_{10}q}.$$

$$\frac{{\mathrm{d}}^{4}N}{\mathrm{d}z\mathrm{d}{\log }_{10}{m}_1\mathrm{d}{\log }_{10}q\mathrm{d}t}=$$

$$\frac{{\mathrm{d}}^{3}n}{\mathrm{d}z\mathrm{d}{\log }_{10}{m}_1\mathrm{d}{\log }_{10}q}\times \frac{4\pi c{d}_L^2}{(1+z{)}^2},$$

$$\frac{{\mathrm{d}}^{5}N}{\mathrm{d}z\mathrm{d}{\log }_{10}{m}_1\mathrm{d}{\log }_{10}q\mathrm{d}{\log }_{10}{f}_{\mathrm{orb}}\mathrm{d}{\log }_{10}e}=$$

$$\frac{{\mathrm{d}}^{4}N}{\mathrm{d}z\mathrm{d}{\log }_{10}{m}_1\mathrm{d}{\log }_{10}q\mathrm{d}t}\times \frac{1}{\Sigma v},$$

$$\displaystyle\frac{~{}{\rm d}f_{\rm orb}}{~{}{\rm d}t}\Big{|}_{\rm gw}=\frac{96G^{5/3}}{5c^{5}}(2\pi)^{8/3}\mathcal{M}^{5/3}\,f_{\rm orb}^{11/3}\,{\cal F}(e),$$

$$\displaystyle\frac{~{}{\rm d}e}{~{}{\rm d}t}\Big{|}_{\rm gw}=-\frac{G^{5/3}}{15c^{5}}(2\pi)^{8/3}\mathcal{M}^{5/3}\,f_{\rm orb}^{8/3}\,{\cal G}(e),$$

$$\mathcal{M}=(m_1{m}_2{)}^{3/5}(m_1+m_2{)}^{-1/5}$$

$$\mathcal{F}(e)$$

$$\mathcal{G}(e)$$

$$\displaystyle\frac{~{}{\rm d}f_{\rm orb}}{~{}{\rm d}t}\Big{|}_{\rm star}$$

$$\displaystyle\frac{~{}{\rm d}e}{~{}{\rm d}t}\Big{|}_{\rm star}$$

$${\log }_{10}\left(\frac{R_{\mathrm{eff}}}{\mathrm{kpc}}\right)=\{\begin{array}{ll}-5.54+0.56{\log }_{10}\left(\frac{M_{\star }}{{\mathrm{M}}_{\odot }}\right)& \text{if }{M}_{\star }>10^{10.3}{\mathrm{M}}_{\odot },\\ & \\ -1.21+0.14{\log }_{10}\left(\frac{M_{\star }}{{\mathrm{M}}_{\odot }}\right)& \text{if }{M}_{\star }\le 10^{10.3}{\mathrm{M}}_{\odot }.\end{array}$$

$$\rho (r)=\frac{(3-\gamma )M_{\star }}{4\pi }\frac{r_0}{r^{\gamma }(r+r_0{)}^{4-\gamma }},$$

$$r_0\approx 1.3R_{\mathrm{eff}}[2^{1/(3-\gamma )}-1].$$

$$\sigma \simeq 0.4\sqrt{\frac{G{M}_{\star }}{r_0}},$$

$$\displaystyle\frac{~{}{\rm d}f_{\rm orb}}{~{}{\rm d}t}\Big{|}_{\rm gas}=\dfrac{3}{2}\mathcal{D}$$

$$\times {\left[\frac{4q}{(1+q{)}^2}\right]}^{-3/8}{\alpha }_{\nu }^{1/2}{\delta }^{-7/8}{f}_{\mathrm{orb}}^{19/12},$$

$$\displaystyle\frac{~{}{\rm d}f_{\rm orb}}{~{}{\rm d}t}\Big{|}_{\rm gas}=\dfrac{3\dot{M}}{\mu}\sqrt{\dfrac{\delta}{1-e^{2}}}\,\,f_{\rm orb},$$

$${\delta }^3=\frac{(1+e_{\lim }{)}^3}{1-e_{\lim }},$$

$$S_n(f)=\frac{10}{3L^2}\left(P_{\mathrm{OMS}}+\frac{4P_{\mathrm{acc}}}{(2\pi f{)}^4}\right)\left(1+\frac{6}{10}{\left(\frac{f}{f_{\star }}\right)}^2\right)+S_c(f),$$

$$P_{\mathrm{OMS}}=2.25\times 10^{-22}\left[1+{\left(\frac{2\mathrm{mHz}}{f}\right)}^4\right]{\mathrm{m}}^2\mathrm{H}{\mathrm{z}}^{-1},$$

$$\displaystyle P_{\rm acc}=9\times 10^{-30}\Biggl{[}1+$$

$$\times \left[1+{\left(\frac{f}{8mHz}\right)}^4\right]{\mathrm{m}}^2\mathrm{H}{\mathrm{z}}^3,$$

$$S_c(f)=A{f}^{-7/3}{e}^{-f^{\alpha }+\beta f\sin \kappa f}\left[1+\tanh (\gamma (f_k-f))\right]\mathrm{H}{\mathrm{z}}^{-1}.$$

$$A$$

$$\alpha$$

$$\beta$$

$$\kappa$$

$$\gamma$$

$$f_k$$

$$\dot{E_n}=\frac{32G^{7/3}}{5c^5}(2\pi \mathcal{M}{f}_{\mathrm{orb}}{)}^{10/3}{g}_n(e),$$

$$g_n(e)$$

$$\displaystyle+2eJ_{n+1}(ne)-J_{n+2}(ne)\bigg{)}^{2}$$

$$\displaystyle+(1-e^{2})\Big{(}J_{n-2}(ne)-2J_{n}(ne)+J_{n+2}(ne)\Big{)}^{2}$$

$$\displaystyle+\frac{4}{3n^{2}}J_{n}^{2}(ne)\Bigg{]}.$$

$$h_{c,n}=\frac{1}{\pi d_L}\sqrt{\frac{2G\dot{E_n}}{c^3{\dot{f}}_n}}$$

$$(S/N{)}^2=\sum _n\int \frac{h_{c,n}^2}{f_n{S}_n(f_n)}\mathrm{d}\ln f_n,$$

Full text

††pubyear: 2019††pagerange: Post-Newtonian Evolution of Massive Black Hole Triplets in Galactic Nuclei: IV. Implications for LISA–LABEL:lastpage

1 Introduction

The existence of massive black hole binaries (MBHBs) is a key prediction of our current $\Lambda$CDM cosmological model (White & Rees, 1978; Kauffmann & Haehnelt, 2000; Volonteri et al., 2003). The ubiquitous presence of MBHs at the centre of bright nearby galaxies (see, e.g. Kormendy & Richstone, 1995; Magorrian et al., 1998; Ferrarese & Merritt, 2000; Gebhardt et al., 2000) and the existence of luminous quasars at redshift $z$ as large as 7 (Mortlock et al., 2011; Bañados et al., 2018) hint at evolutionary pathways in which galactic haloes, growing through repeated mergers and accretion of cold gas along dark matter filaments, drag their central MBHs to form close pairs on kpc scales. In the interaction with the stellar and gaseous background of the new galaxy, gravitational and hydrodynamical torques control and guide the sinking of the two MBHs toward the galactic core where coalescence takes place (see e.g. Colpi, 2014, for a review).

Coalescing MBHBs forming in the aftermath of galaxy collisions are among the loudest sources of gravitational waves (GWs) in the Universe, and key targets of the ESA led space mission LISA, the Laser Interferometer Space Antenna (Amaro-Seoane et al., 2017). For the first time, LISA will survey the GW frequency window from about 100 $\mu$Hz to 100 mHz, where the late inspiral, merger and ring-down of MBHBs will be best observed for massive black holes (MBHs) with masses between $10^{4}~{}\rm M_{\large\odot}$ and $10^{7}~{}\rm M_{\large\odot}.$ While the galaxy mass function suggests that these low mass MBHs may be the most common, they are the least known in terms of basic demographics, birth, growth, dynamics and connection to their galaxy hosts (Kormendy & Ho, 2013). Nevertheless, LISA’s exploration volume is immense, and will allow us to detect GW signals from MBHBs coalescing at redshifts as large as 20, when the Universe was only 180 Myrs old, down to the climax of star formation at $z\sim 2$ (Madau & Dickinson, 2014). By observing deep into the high-redshift Universe LISA will tell us when the first MBHBs formed and how they evolved across cosmic epochs, providing the first census of this new population of GW sources (Haehnelt, 1994a; Jaffe & Backer, 2003; Wyithe & Loeb, 2003; Enoki et al., 2004; Sesana et al., 2004, 2005; Jenet et al., 2005; Rhook & Wyithe, 2005; Barausse, 2012; Klein et al., 2016).

However, in order to fully plan (and exploit) LISA observations, theoretical models are needed to characterise this new population of sources and make predictions for the rate of MBHB coalescences as a function of redshift. This rate depends sensitively on the MBHB mass and mass ratio, which are related to those of the interacting galaxies, and on the mechanisms leading to coalescence. Tracing the dynamics of MBHBs during a galactic merger is extremely challenging, as the process is multi-scale and depends on a variety of physical processes. Indeed, LISA MBHBs encompass a wide interval of masses, from $10^{4}~{}\rm M_{\large\odot}$ to a few $10^{7}~{}\rm M_{\large\odot}$, mass ratios $q$ from $\sim 0.1$ to $\sim 1$, and coalescence redshifts ranging from $z\sim 20$ to the present, implying very different environments and initial conditions for each binary. Moreover, the dynamical range is huge and MBHs have to travel long distances before coalescence, starting from the typical scale of a galactic merger (i.e. tens to hundred kpc, depending on the encounter geometry), through mpc separations (where GWs start driving the inspiral) and down to $\mu$pc scale where coalescence takes place. Inside the newly formed galaxy, gas and stars concur to extract MBH orbital energy and angular momentum, thus affecting the evolution, but these interactions are complex because of feedback processes (from AGNs and supernovae) and the multi-phase nature of the gas component.

Furthermore, bottlenecks are possible along the way to coalescence. MBHs can wander in the outskirts of the new galaxy if the mass ratio of the progenitors is small (below about $0.1$), as tidal stripping and ram pressure can remove stars and gas in the MBH surroundings, reducing the effectiveness of dynamical friction at dragging the MBHs to the nuclear region (Taffoni et al., 2003; Callegari et al., 2011). This can produce a population of wandering black holes on scales in excess of 10-100 pc (Tremmel et al., 2017; Dosopoulou & Antonini, 2017). In gas-rich environments and again on the large scales, formation of massive gas clumps (before they fragment to form stars) can lead to MBH stochastic dynamics for comparable masses, and to scattering of the MBHs off the disc plane under rather generic conditions. This may widen the distribution of MBHB coalescence times (measured as “delays” from the galaxy merger), which can range from $\approx 100$ Myrs to several Gyrs. Individual mergers, which can be studied using full gas-free N-body simulations (Berczik et al., 2006; Khan et al., 2011; Sesana et al., 2011b; Khan et al., 2012a, b, 2013; Vasiliev, 2014; Vasiliev et al., 2015; Vasiliev, 2016; Gualandris et al., 2017), high resolution SPH simulations of galaxies in isolation (Mayer et al., 2007; Dotti et al., 2007; Callegari et al., 2008; Lodato et al., 2009; Callegari et al., 2011; Fiacconi et al., 2013; Capelo et al., 2015; Roškar et al., 2015; del Valle et al., 2015; Mayer et al., 2016; Goicovic et al., 2017; Lupi et al., 2015; Tamburello et al., 2017) or in fully cosmological contexts (Vogelsberger et al., 2014; Khan et al., 2016; Schaye et al., 2015; Tremmel et al., 2017; Tremmel et al., 2018), can indeed reveal the complexity of the black hole interaction with the (large scale) environment.

On smaller scales $\sim$pc, once the MBHs have formed a binary, stellar hardening can become a major process guiding the journey to coalescence. Stars are scattered away individually from the binary, which leads to loss-cone depletion in phase-space, and eventually binary stalling (Yu, 2002). Recently, however, it has been demonstrated that regimes of full loss-cone can be attained by enhanced stellar diffusion in phase space due to the presence of triaxiality in the galaxy’s gravitational potential (Yu, 2002; Khan et al., 2011; Khan et al., 2013; Vasiliev, 2014; Gualandris et al., 2017) or to galactic rotation (Holley-Bockelmann & Khan, 2015). Alternatively, if a significant amount of gas is present in the form of a (massive) circumbinary disc, the evolution of the MBHB towards the GW-dominated regime may be faster, because MBHBs are expected to undergo planetary-like migration (Haiman et al., 2009; Dotti et al., 2015). A third process can drive a MBHB down to coalescence when the galaxy host is not isolated. Particularly at high redshifts, when galaxy mergers are commonplace, a third MBH can approach the binary as a result of dynamical friction, forming an interacting ‘triplet’. If the triplet is hierarchical, the binary may undergo the Kozai-Lidov oscillations (Kozai, 1962; Lidov, 1962), which tend to secularly increase the eccentricity of the inner binary, eventually driving it to coalescence. In addition, whenever the triplet becomes unstable, chaotic 3-body interactions may further enhance the coalescence rate (Blaes et al., 2002; Hoffman & Loeb, 2007; Amaro-Seoane et al., 2010; Kulkarni & Loeb, 2012; Bonetti et al., 2016, 2018b; Ryu et al., 2018).

If single detailed simulations can resolve the great complexity of the MBH interactions, semi-analytical models (SAMs) provide the large picture of MBHB coalescence (Kauffmann & Haehnelt, 2000; Volonteri et al., 2003; Sesana et al., 2011a; Barausse, 2012; Valiante et al., 2016; Ricarte & Natarajan, 2018a; Dayal et al., 2019). SAMs are optimally suited to perform statistical studies of the formation and evolution of binary and multiple MBH systems. While less realistic than full-fledged simulations, SAMs are vastly superior from the point of view of computational performance, thus allowing a systematic exploration of many different scenarios for the formation and co-evolution of MBHs within their host galaxies. In particular, as shown by Klein et al. (2016), the MBHB merger rate can reduce significantly when astrophyiscally motivated “delays” times are accounted for in a merger scheme, also leading to a shift to lower $z$ of the maximum rate.

In this paper we combine the SAM developed by Barausse (2012) (and later extensions) with a self-consistent treatment of the dynamical interactions driving MBHBs to coalescence. In particular, we highlight the effect and evolution of the orbital eccentricity during the standard stellar and gas-hardening phases, and also in the context of a novel coalescing channel involving the formation of MBH triplets (Blaes et al., 2002; Hoffman & Loeb, 2007; Kulkarni & Loeb, 2012; Bonetti et al., 2018b). This channel arises naturally when MBHBs have a “long” hardening time. Under these circumstances, a subsequent galaxy merger may deliver a third MBH that is close enough to the MBHB to promote coalescence. In this context Bonetti et al. (2018b, a) explored a large sample of 3-body numerical simulations, with the goal of implementing this triplet channel in SAMs.

Our aim is to characterise MBHB coalescences and forecast LISA detection prospects with physically motivated models for the time delays between galaxy and MBHB mergers.111The method was previously applied in the context of the nHz pulsar timing array experiments, which target GW signals and backgrounds emitted by $10^{8-9}~{}\rm M_{\large\odot}$ MBHBs at low redshifts (Bonetti et al., 2018a). We consider two conceptually different dynamical scenarios: a fiducial model in which MBHBs merge on timescales of millions-to-billions of years, depending on the properties of the host galaxy (and in particular on its gas and star content); and an extreme model in which all MBHBs stall at about their hardening radius, so that their coalescence is prompted only by the interaction with a third/multiple MBHs. Clearly, the latter is the most pessimistic scenario, yielding to the lowest LISA rates.

The paper is organised as follows: in Section 2 we briefly review the main features of the SAM adopted to track structure evolution. In Section 3 we describe the physical processes that drive MBHBs to coalescence, including stellar hardening, gas-driven migration and triplet dynamics, and how we generate a mock catalogue of LISA events. In Section 4 we detail the formalism employed to compute waveforms and signal-to-noise ratios (hereafter S/N), and in Section 5 we present our main results. Finally, Section 6 is devoted to a discussion of the results and to our conclusions.

2 Semi-analytic galaxy formation model

SAMs of galaxy and MBH formation and evolution represent a very flexible tools to explore MBH merger rates at a fairly modest computational cost when compared to the more sophisticated, but much dearer, hydrodynamical simulations (e.g. Springel, 2010; Guillet & Teyssier, 2011; Vogelsberger et al., 2014; Schaye et al., 2015; Rodriguez-Gomez et al., 2015; Kelley et al., 2017; Wadsley et al., 2017; Tremmel et al., 2018). For this reason, they have been extensively used, among other applications, to compute projected merger rates for MBHs (Sesana et al., 2004, 2005; Sesana et al., 2007; Sesana et al., 2011a; Tamanini et al., 2016; Klein et al., 2016) and (more recently) extreme mass ratio inspirals (Babak et al., 2017) in the LISA band, as well as stochastic background amplitudes and resolved event rates for pulsar timing arrays (Sesana et al., 2009; Dvorkin & Barausse, 2017; Ryu et al., 2018; Bonetti et al., 2018a).

In this work, we utilise the SAM of Barausse (2012), with the improvements described in Sesana et al. (2014), Antonini et al. (2015a) and Antonini et al. (2015b). This model, already employed in Klein et al. (2016); Dvorkin & Barausse (2017); Babak et al. (2017); Bonetti et al. (2018a) to explore the detection prospects of the future LISA mission and current pulsar timing array campaigns, evolves the baryonic galactic structures (including MBHs) along dark matter merger trees produced with an extended Press-Schechter formalism, modified to reproduce the results of N-body simulations (Press & Schechter, 1974; Parkinson et al., 2008). We have checked that our results are robust against the time step of the merger tree (see e.g. Shimizu et al., 2002). Besides the MBHs, the code also follows the evolution of several baryonic components. These include for instance a hot chemically primordial intergalactic medium, and a cold chemically enriched interstellar medium, which eventually forms stars. Both the cold gas and stellar components exist as discs and/or spheroids. On smaller scales, the model includes a pc-scale nuclear gas reservoir  (Barausse, 2012; Antonini et al., 2015a, b),222If the nuclear gas reservoir is more massive than the hosted MBH, we define a galaxy as gas-rich. whose growth rate is linearly correlated with star formation in the spheroid (Granato et al., 2004; Haiman et al., 2004; Lapi et al., 2014), which in turn takes place in our model after major galactic mergers and bar instabilities in the galactic disc. This reservoir is available to accrete onto the central MBH on the gas viscous timescale (Sesana et al., 2014) and/or to form a nuclear star cluster (Antonini et al., 2015a, a). Besides forming from this in-situ star formation in the nuclear gas, the nuclear star cluster also forms from migration of globular clusters to the galactic centre (Antonini et al., 2015a, b). Supernova and AGN feedback are included in the model. The latter, together with the aforementioned prescriptions for MBH accretion, ensures that MBHs grow to reproduce the local observed correlations with host mass and velocity dispersions at $z=0$ (Barausse, 2012; Sesana et al., 2014; Barausse et al., 2017).

The features that most impact the MBH merger rates for LISA are the initial mass function of the MBH seeds at high redshift, and the mechanisms (with their related timescale) which drive MBHBs to coalesce after a galaxy merger. Regarding the former, we adopt two competing scenarios for the seeds. In a “light-seed” (LS) scenario, we make the hypothesis that MBH seeds originate from the remnants of pop III stars forming in high-redshift, low-metallicity galaxies (e.g. Madau & Rees, 2001). To model this process, we assume that only the rarest, most massive haloes (collapsing from the $3.5\sigma$ peaks of the primordial density field at redshift between 15 and 20) contain a MBH seed (Madau & Rees, 2001). The mass of this seed is assumed to be $\sim 2/3$ of the initial pop III mass (to account for mass losses during stellar collapse), which is in turn chosen randomly from a log-normal distribution centred on 300 $M_{\odot}$ and with r.m.s. of 0.2 dex, with a mass gap between 140 and 260 $~{}\rm M_{\large\odot}$, to account for the formation of pair-instability supernovae (Heger & Woosley, 2002).

In the “heavy-seed” (HS) scenario, we model MBH seeds as forming from the collapse of protogalactic discs as a result of bar instabilities. Following Volonteri et al. (2008) we describe the process by a single free parameter, i.e. the critical Toomre parameter $2\lesssim Q_{c}\lesssim 3$, marking the onset of instability. We set $Q_{c}=2.5$, and we only seed haloes at $15<z<20$ as in Volonteri et al. (2008).

As for the timescale on which MBHBs coalesce after a galaxy merger333In our SAM, galaxy mergers take place a dynamical friction time (Boylan-Kolchin et al., 2008) after their host haloes merge. We also account for environmental effects on satellite haloes (tidal stripping and evaporation) following Taffoni et al. (2003). (i.e. the “delays” between galaxy and MBHB coalescences), we adopt, as in Bonetti et al. (2018a), two different prescriptions. In a more realistic scenario (labelled as Model-delayed) we assume that MBHBs form promptly after a galaxy merger, and that their evolution is driven by a combination of three possible processes: stellar hardening (on timescales $\lesssim$ few Gyr), gas-induced migration (on timescales of $\sim 10^{7}-10^{8}$ yr if enough gas is available), and triple MBH interactions (on timescales of $\sim 10^{8}-10^{9}$ yr), if a third MBH is brought in by a subsequent galaxy merger.

In a more pessimistic scenario (Model-stalled), we instead make the assumption that all MBHBs stall, such that coalescences can only occur because of the formation of a MBH triplet, which can promote binary mergers via either Kozai-Lidov resonances (Kozai, 1962; Lidov, 1962) or chaotic three-body interactions (Bonetti et al., 2016, 2018b). The evolution of MBH triplets is modelled by fitting a set of simulations (Bonetti et al., 2018b) produced with the three-body post-Newtonian code presented in Bonetti et al. (2016).

In summary, we consider four different evolutionary models for MBHs in galactic haloes, according to the seeding and delay prescriptions. For each model, the SAM provides us with a set of events labelled by coalescence redshift, mass of the primary MBH $m_{1}$, mass ratio $q$ (i.e. $q=m_{2}/m1\leq 1$, being $m_{2}$ the less massive binary component), and a flag identifying the evolutionary path followed, i.e. stellar-, gas-, triplet-driven dynamics. This information is then used to construct the comoving number density of coalescing MBHBs $n$, per unit redshift and (logarithmic) unit of $m_{1}$ and $q$:

[TABLE]

Equation (1) above can be easily converted into the number $N$ of coalescences per unit observer time as Haehnelt (1994b)

[TABLE]

where $d_{L}$ denotes the luminosity distance.

3 Synthetic catalogues of coalescing MBHBs

Our aim is to produce a simulated catalogue of events observable by LISA. To that purpose, we need to convert the intrinsic merger rate (equation 2) into an instantaneous number of systems across the sky per (logarithmic) units of eccentricity $e$ and observed orbital frequency $f_{\rm orb}$. It is useful to relate the evolution of binaries in the $(\log_{10}f_{\rm orb},\log_{10}e)$ plane as the flow of a fluid in a two dimensional space. If the flow is stationary, the continuity equation implies that the fluid’s density along a flow tube is inversely proportional to the velocity $v$ and to the flow tube’s cross section $\Sigma$. We can thus write

[TABLE]

where ${1}/(\Sigma v)=\tau\times{\rm d}p/({\rm d}\log_{10}f_{\rm orb}{\rm d}\log_{10}e)$ is the product of the average binary lifetime $\tau$ and the probability density $~{}{\rm d}^{2}p/({\rm d}\log_{10}f_{\rm orb}{\rm d}\log_{10}e)$ that a binary is observed at a given $f_{\rm orb}$ and $e$. The practical computation of this quantity is described in Sec. 3.2. The normalisation ${{\rm d}^{4}N}/({\rm d}z~{}{\rm d}\log_{10}m_{1}{\rm d}\log_{10}q\,{\rm d}t)$ is provided by our SAM. Note that equation (3) generalises the treatment of Sesana et al. (2008) to binaries with eccentricity.

The SAM informs us of the type of evolutionary path (stellar-, gas- or triplet-driven) followed by binaries, without any information about the time evolution of the orbital elements. Such details are however necessary in order to compute $1/(\Sigma v)$ in the specific environment of the galaxy merger remnant. As already discussed, we consider two “standard” evolutionary channels, i.e. stellar-driven dynamics and gas-driven dynamics (c.f. also Antonini et al., 2015b; Klein et al., 2016), and a third, less standard dynamical process involving the formation and evolution of triple systems of MBHs (Bonetti et al., 2016, 2018b, 2018a). All such paths are, in principle, followed by a final, short phase dominated by GW energy losses. In the following, we will discuss in detail our implementation of MBHB dynamics.

3.1 Dynamics of MBHBs

3.1.1 GW-driven dynamics

Irrespective of the actual processes involved in the early stages of the orbital evolution of a MBHB, the last stages will always be dominated by GW radiation losses. For the first two evolutionary channels mentioned above (i.e. stellar and gas driven dynamics), we model the GW phase adopting the orbital averaged equations derived by Peters (1964), which describes the time evolution of the orbital frequency $f_{\rm orb}$ and eccentricity $e$:

[TABLE]

where $\mathcal{M}$ is the binary source-frame chirp mass, defined as

[TABLE]

while ${\cal F}(e)$ and ${\cal G}(e)$ are algebraic functions of the eccentricity:

[TABLE]

In the case of triple MBH interactions, we instead adopt the more refined treatment of Bonetti et al. (2016), briefly described in the next Section 3.1.4.

3.1.2 Stellar-driven dynamics

In the aftermath of a merger involving two gas-poor galaxies, the main process driving the early evolution of a MBHB is the scattering of background stars (Mikkola & Valtonen, 1992; Quinlan, 1996; Sesana et al., 2006). Stars intersecting the MBHB orbit can extract a non-negligible fraction of the binary energy and angular momentum because of the slingshot mechanism. The secular, cumulative effect of a large number of 3-body encounters is effective at hardening the MBHB. Eventually, when the orbital semi-major axis shrinks enough, GW emission takes over and the binary quickly evolves to coalescence. The whole process can be described by the relatively simple, albeit approximate, coupled differential equations:

[TABLE]

where $M=m_{1}+m_{2}$ is the binary total mass, $\sigma$ is the velocity dispersion of the host galaxy, $\rho_{i}$ is the stellar density at the binary influence radius (i.e. the radius containing twice the MBHB mass in stars), whereas $H$ and $K$ represent numerical factor calibrated against scattering experiments (see Quinlan, 1996; Sesana et al., 2006, for complete details). The efficiency of the stellar hardening saturates at the hardening radius, $a_{h}\sim Gm_{2}/4\sigma^{2}$ (where $m_{2}$ is the mass of the secondary, Quinlan, 1996), and beyond that point the MBHB hardens at a constant rate.

In order to make use of equations (9-10), we need to link the properties of the stellar distribution to the mass of the hosted MBHB, i.e. we need to evaluate $\rho_{i}$ and $\sigma$ for a given $M$. These quantities can be inferred once the stellar mass $M_{\star}$ has been specified. For a given $M$, we numerically sample $M_{\star}$ from the stellar mass distribution provided by the SAM,444This procedure is needed since the evaluation of $1/(\Sigma v)$ in equation 3 is associated with discrete values of $m_{1}$ and $q$ (as detailed in Section 3.2). and we then evaluate an effective radius for the stellar spheroid with the following scaling relation (see Shen et al., 2003):

[TABLE]

Next, we model the stellar density profile of the spheroid by assuming a Dehnen profile (being $\gamma$ the inner logarithmic slope; see Dehnen, 1993), i.e.

[TABLE]

where the scale radius $r_{0}$ is inferred from the galaxy effective radius $R_{\rm eff}$ as in Dehnen (1993),

[TABLE]

The density at the influence radius $\rho_{i}$ is then readily obtained from equation (12) by setting $r=r_{\rm inf}$. Finally, the velocity dispersion $\sigma$ is evaluated (under the assumption that the virial theorem holds) from

[TABLE]

where the numerical pre-factor takes into account possible anisotropies in the stellar phase space (see e.g. Baes & Dejonghe, 2002).

Throughout the paper we consider two fairly different density profiles, characterised by different values of the inner logarithmic slope in equation 12: $\gamma=1$, which gives relatively low central stellar densities, resulting in a rather weak coupling between the MBHB and its stellar environment; and $\gamma=1.75$, which instead gives much higher central densities, resulting in much stronger coupling. We initialise the MBHB at a separation such that the enclosed mass in stars matches the mass of the secondary MBH, and at this separation we draw a random eccentricity from a thermal distribution ($dp/de=2e$). We finally evolve the binary until coalescence, using equations (4-5) and equations (9-10).555We assume that the two MBHs merge when they reach a separation lower than six gravitational radii, i.e. $a<6GM/c^{2}$.

3.1.3 Gas-driven dynamics

In a gas-rich galaxy merger, we model the early evolution of the MBHB adopting a simplified description of type II migration in circumbinary discs (Artymowicz & Lubow, 1994; Gould & Rix, 2000; Haiman et al., 2009; Roedig et al., 2011; del Valle & Escala, 2012; Kocsis et al., 2012; Dotti et al., 2015). In order to gain insight into the orbital evolution of the binary, we consider two different models, i.e. the fairly sophisticated description of type II migration proposed by Haiman et al. (2009) (HKM hereinafter), and the simpler model considered by Dotti et al. (2015) (DMM herein after).

In the HKM model, the binary is embedded in a thin circumbinary disc, with the plane of the disc aligned with the binary orbit. The orbital decay is dominated by viscous angular momentum exchange with the disc. In order to describe the evolution quantitatively, HKM assumed that the circumbinary gas forms a standard Shakura & Sunyaev accretion disc (Shakura & Sunyaev, 1973). HKM found that just before the transition to GW-driven evolution, the viscous orbital decay is in the “secondary-dominated” Type-II migration regime (i.e. the mass of the secondary is larger than the enclosed disc mass). This is slower than disc-dominated Type-I migration. For the HKM model, the orbital frequency evolves as

[TABLE]

where $\mathcal{D}=4.5\times 10^{-43}{\rm m^{7/8}kg^{-1/8}s^{-1}}$ is a constant that depends on the details of the model. The numerical value reported here is specific of a $\beta-$disc with fiducial values of the involved parameters We further use $\alpha_{\nu}=0.3$ for the disc viscosity parameter and $\epsilon=0.1$ for the radiative efficiency (see Haiman et al., 2009, for full details). The accretion rate in Eddington units $\dot{m}\equiv\epsilon\dot{M}c^{2}/L_{\rm Edd}$ is modelled as in Appendix A of Tamanini et al. (2016), where we curbed the accretion to Eddington. Finally, $\delta=(1+q)(1+e)$ is the size of the gap (in units of the binary semi-major axis, see e.g. Artymowicz & Lubow, 1994).

DMM give an order of magnitude estimate of the binary coalescence time-scale, proposing a simple model for the interaction between a MBHB and a corotating circumbinary accretion disc. Under this configuration, the gravitational torque exerted by the MBHB stops the gas inflow at a separation of $\delta\sim 2$ (see, e.g. Artymowicz & Lubow, 1994). At this radius the gravitational torque between the binary and the disc is perfectly balanced by the torques that determine the radial gas inflow on large scales. The analogous to equation (3.1.3) for the simpler DMM model is

[TABLE]

where $\dot{M}$ represents here the accretion rate within the disc, and $\mu$ the reduced mass of the binary. In general, the evolution implied by this last equation is somewhat faster compared to equation (3.1.3). An equal mass binary with $M=2\times 10^{6}$ $M_{\odot}$ reaches the GW dominated regime in $\simeq 1.4\times 10^{8}$ yrs, compared to $\simeq 1.8\times 10^{8}$ yrs in the HKM case.

Note that the two considered models describe the orbital evolution as a change in the orbital frequency only, with the eccentricity unaffected. Indeed, as suggested by e.g. Roedig et al. (2011), during type II migration the binary eccentricity converges to a specific value, $e_{\rm lim}\sim 0.7$, with only a mild dependence on the mass ratio. The value of the limiting eccentricity is obtained in Roedig et al. (2011) by balancing the disc angular frequency with that of the secondary at apocentre. This leads to an expression that relates $e_{\rm lim}$ to the size of the gap (in units of the binary semi-major axis) $\delta$:

[TABLE]

which, once solved for $e_{\rm lim}$, yields the limiting eccentricity. We then initialise the initial MBHB eccentricity by sampling from a Gaussian distribution with mean equal $e_{\rm lim}$ and standard deviation 0.05. We note, however, that the particular adopted ansatz does not have a strong influence on the properties of the observed LISA events, as all gas-driven mergers will present a negligible residual eccentricity (see next Section LABEL:sec:results).

3.1.4 Triple MBH dynamics

In the two previous subsections we have described the binary orbital evolution in the cases of dry and wet galaxy mergers, respectively. In gas-rich galaxies, the actual efficiency of gas-MBHB interaction in realistic physical conditions is poorly known, and stalling of the binary is still a possibility (see, e.g. Lodato et al., 2009). In dry galaxy mergers, we have seen that the efficiency of the stellar hardening saturates at the hardening radius and, unless new stars are continuously supplied to the otherwise depleting loss-cone, the whole process is fated to stop. Because typically $a_{h}\sim 0.1$ pc for $\sim 10^{6}~{}\rm M_{\large\odot}$ black holes, it is not guaranteed that the MBHB can eventually close the gap down to separations $a_{\rm gw}\sim 10^{-3}$ pc, i.e. the separations at which GW emission alone can drive the two MBHs to coalesce within a Hubble time.

A possible mechanism that could solve this problem is provided by triple MBH interactions. Triple MBH systems can form when a MBHB stalled at separations $\lesssim a_{\rm h}$ interacts with a third MBH carried by a new galaxy merger (see, e.g. Mikkola & Valtonen, 1990; Heinämäki, 2001; Blaes et al., 2002; Hoffman & Loeb, 2007; Kulkarni & Loeb, 2012). The net outcome of a series of complex three-body interactions can effectively increase the binary eccentricity, reaching values close to unity and leading to copious GW losses, hence driving MBHBs to coalescence. A thorough assessment of the process requires a detailed description of the three-body dynamics, which we have carried out in a series of papers (Bonetti et al., 2016, 2018b, 2018a). Here we summarise the main features and results of that work.

We numerically integrate the orbits of MBH triplets formed by a stalled MBHB at the centre of a stellar spherical potential, and by a third MBH approaching the system from larger distances. The employed numerical scheme directly integrates the three-body (Hamiltonian) equations of motion up to 2.5PN order, introducing velocity-dependant forces to account for the dynamical friction on the intruder during its initial orbital decay toward the galactic centre, and for the stellar hardening (Quinlan, 1996) of the outer binary as described in Section 3.1.2.

We found that over a large sample of $\sim$ 15000 simulated systems, a fraction of about 30% experiences the merger of any one pair of MBHs in the triplet. The typical coalescence timescale is $\sim$ a few $\times 10^{8}$ yrs, a timescale dominated by the orbital decay of the third MBH. Finally, extremely important for the study we are presenting here is the finding that coalescing binaries are generally driven to eccentricities, in excess of 0.9 and up to 0.9999 in some cases. Moreover, even at coalescence, binaries can retain a non-zero eccentricity, up to 0.1. Binaries driven by triple MBH interactions can therefore have large eccentricities when entering the LISA band.

3.2 Catalogue compilation

The quantity ${1}/(\Sigma v)=\tau\times{\rm d}^{2}p/({\rm d}\log_{10}f_{\rm orb}{\rm d}\log_{10}e)$ appearing in equation (3) is evaluated numerically following the evolution of a sample of binaries over the 2-dimensional space $(\log_{10}f_{\rm orb},\log_{10}e)$, covered by a grid in the range $-14\leq\log_{10}f_{\rm orb}\leq-1$ and $-5\leq\log_{10}(1-e)\leq 0$. Binaries are selected according to a combination of primary mass and mass ratio. Typically, for each pair $(m_{1},q)$, we simulate from few tens to few thousands systems, according to the evolutionary channel followed as described in Section 3.1. In practice, for a specific binary characterised by $(m_{1},q)$, we evaluate the time that the system spends in each bin of the 2-dimensional grid of frequency and eccentricity. We then sum over all binaries, and divide by the total number of simulated systems and by the bin area to compute ${1}/(\Sigma v)$ (see Fig. 9 in Bonetti et al., 2018b, for a concrete example.).

Our binaries are simulated for 45 combinations of primary mass and mass ratio, i.e. $\log_{10}m_{1}=2,3,4,5,6,7,8,9,10$ and $\log_{10}q=-1.5,-1.0,-0.5,-0.25,0.0$, resulting in a coarse coverage of the parameter space. Given the strong dependence of the GW losses on the binary parameters, we would need a much finer sampling of the $(m_{1},q)$ plane. In order to interpolate in the $(m_{1},q)$ plane, we exploit the fact that the evolution of GW dominated binaries, which is described by equations (4-5), depends in a simple way on $m_{1}$ and $q$. Indeed, one can rewrite equations (4-5) in terms of the logarithm of the dimensionless frequency $\tilde{f}_{\rm orb}=f_{\rm orb}GM/c^{3}$ and in terms of the dimensionless time $\tilde{t}=t\mathcal{M}^{5/3}M^{-8/3}c^{3}/G=t\eta/(GM/c^{3})$, where $\eta$ is the symmetric mass ratio $\eta=m_{1}\,m_{2}/M^{2}$. In the variables ($\log_{10}\tilde{f}_{\rm orb},e,\tilde{t}$) the dependence on $m_{1}$ and $q$ disappears (as expected as general relativity has no preferred scale). This allows us to interpolate our results for ${1}/(\Sigma v)$ in the $(m_{1},q)$ space employing a bilinear interpolation among the four closest pairs of $(m_{1},q)$.

It is important to stress again that the above procedure is valid only when the binary evolution in the $(\log_{10}f_{\rm orb},\log_{10}e)$ plane is driven by GW losses. As discussed in the previous section, before the GW-dominated regime several astrophysical processes other than GW emission drive the binary evolution. In these earlier evolutionary phases our interpolation scheme is indeed not correct. However, GWs always dominate the dynamics of MBHBs in the LISA frequency band, hence our procedure is justified.

In order to create a catalogue of LISA events, we now need to sample the distribution given by equation (3). Using a standard rejection technique, we randomly generate a number of MBHBs equal to the normalization of our sampled distribution, and average over 100 realizations. In order to speed up the computation, we limit the sampling to a given minimum value of the (observed) orbital frequency. Equation (3) effectively represents a snapshot of all binaries in our synthetic Universe with observed orbital frequency larger than the chosen minimum threshold (i.e. 10-100 nHz). Then, the binaries are evolved under the effect of GW emission for the nominal LISA mission lifetime of 4 years (observer time). For each system we finally compute the S/N, as detailed in Section 4.2.

As an example, we report in Fig. 1 the sampling for the triplet channel of Model-stalled with the LS prescription. In each panel, the upper blue histogram represents the marginalisation of the 5-dimensional distribution (equation 3) over the remaining 4 parameters. The lower blue histogram is the sub-sample of MBHBs that are bound to coalesce in 4 years. The red line in the first three leftmost panels shows the same distribution but obtained directly by multiplying the merger rate (equation 2) for 4 years; as expected, this distribution matches the randomly generated sub-sample.

It is interesting to note that the circularity distribution (where the circularity is defined as $1-e$, rightmost panel of the figure) extends up to very low values, i.e. to eccentricities very close to unity. The implication of highly eccentric binaries for the LISA event rate will be discussed in the next sections.

4 Gravitational wave emission and detection

The synthetic catalogue of binaries observable by LISA described in previous sections is then used to characterise the shape and strength of the emitted GW signal. After a brief review of the LISA sensitivity, we describe the formalism employed to compute waveforms in the various stages of the MBHB orbital evolution.

4.1 LISA sensitivity

Throughout the paper, we adopt a six-link LISA configuration (i.e. one consisting of two independent detectors). For each independent detector, we adopt the sky-averaged, single detector LISA sensitivity curve (updated to February 2018), as detailed in the “LISA Strain Curves” document (referenced as LISA-LCST-SGS-TN-001666See also https://atrium.in2p3.fr/nuxeo/nxpath/default/Atrium/sections/Public/LISA@view_documents?tabIds=%3A). Besides the instrumental noise we also take into account the effect of a large population of unresolved galactic compact binaries (mostly white-dwarf binaries) as given by Cornish & Robson (2018). The stochastic background from that population produces a “confusion noise” that degrades the instrumental sensitivity at frequencies below $\sim 1$ mHz. Combining the instrumental and confusion noise contributions, the total LISA sensitivity as a function of the GW frequency $f$ is

[TABLE]

where $L=2.5$ Gm is the LISA arm-length, $f_{\star}=19.09$ mHz, with

[TABLE]

[TABLE]

and

[TABLE]

The last term, $S_{c}(f)$, is the aforementioned confusion noise. In equation (21), for a 4-year mission, the parameters $A,\alpha,\beta,\kappa,\gamma$ and $f_{k}$ take the following values:

[TABLE]

4.2 Waveforms and S/N calculation

As already discussed, in the LISA frequency domain ($10^{-4}\lesssim f\lesssim 1$ Hz) the dynamics is driven by GW losses. In this section we briefly describe the formalism adopted to compute the S/N in the inspiral, merger and ringdown phases. Because, as we will see in Section 5, triple MBH interactions can lead binaries to develop large eccentricities in the LISA band, we will focus in particular on including orbital eccentricity in the waveforms. Inspiral-merger-ringdown waveforms including eccentricity typically assume small to moderate eccentricities (Huerta et al., 2017, 2018), which is not necessarily adequate for our purposes (c.f. Fig. 9). For this reason, we have decided to rather opt for a simplified model accounting for (arbitrarily large) eccentricity in the inspiral phase at leading Newtonian order, while neglecting higher order PN corrections as well as MBH spins. In more detail, we use the Newtonian fluxes from Peters (1964) to evolve the semi-major axis and eccentricity of a given MBHB, while the S/N (of the inspiral phase) is computed as in Barack & Cutler (2004) by converting the Newtonian fluxes (in the various harmonics) into characteristic GW amplitudes. To further account for the effect of the merger-ringdown phase on the S/N, we also augment the inspiral S/N by a merger-ringdown contribution, which we compute from the post-inspiral portion of PhenomC waveforms. In the following, we will describe our procedure in more detail.

4.2.1 Inspiral

In the inspiral phase, we evolve MBHBs in the GW regime with the standard orbit-averaged dynamics given by equations (4-5) (Peters, 1964). The computation of the S/N requires knowledge of the emitted waveform. Gravitational radiation emitted by an eccentric binary requires a more complicated treatment compared to the standard case of a circular orbit. While for circular binaries most of the GW power is contained in the dominant quadrupolar mode, whose frequency is twice the orbital frequency, several harmonics are excited with comparable amplitudes in the eccentric case, i.e. the GW spectrum contains several dominant frequencies $f_{n}=nf_{\rm orb}$, where $n$ is the harmonic number. Therefore the total emitted power in GW is spread on a broad spectrum of frequencies, with the fraction of power per harmonic given by Peters & Mathews (1963)

[TABLE]

where the dimensionless function $g_{n}(e)$ (see Fig. 2) is related to the power carried by each harmonic, and is given by

[TABLE]

Here $J_{n}$ represents the $n$-th order Bessel function of the first kind. Note that in the case of a circular binary $g_{n}(0)=\delta_{2n}$, where $\delta_{mn}$ is the standard Kronecker delta.

Finally, in order to evaluate the S/N, we employ the formalism developed by Barack & Cutler (2004). The characteristic strain of each harmonic is given by

[TABLE]

where $\dot{f}_{n}=n\dot{f}_{\rm orb}$ is the time derivative of the $n$-th harmonic. The total S/N is then computed as

[TABLE]

where $S_{n}(f_{n})$ is the sky-inclination-polarization averaged power spectral density of LISA, which according to its definition in Section 4.1 already accounts for the fact that the LISA constellation is comprised of two independent interferometers.

4.2.2 Merger and Ringdown

The PhenomC model provides a useful approximation to the frequency domain waveforms from circular black hole binaries with non-precessing spins, by “hybridising” the inspiral post-Newtonian amplitude and phase with numerical relativity results for the coalescence phase (Santamaría et al., 2010). However, since they assume circular binaries, these waveforms cannot be used to describe the inspiral phase of the systems in our catalogues, for which eccentricities can be very large in the early inspiral (c.f. Fig. 9). However, because the eccentricity steadily decreases along the MBH inspiral under the effect of GW emission Peters (1964), all of our sources have an eccentricity $e\lesssim 0.1$ at plunge (c.f. Fig. 8 and also Fig. 10 of Bonetti et al., 2018b). We can therefore use PhenomC waveforms to approximately compute the S/N of the plunge-merger-ringdown phase and sum this contribution to the inspiral S/N computed as in Section 4.2.1. We assume that this phase corresponds to frequencies $>2\,f_{\rm ISCO}$, where $f_{\rm ISCO}$ is the orbital frequency at the last stable circular orbit (ISCO) of a Schwarzschild space-time, i.e.

[TABLE]

Examples of the resulting piecewise waveforms in the frequency domain are plotted in Fig. 3 against the LISA sensitivity. While this procedure is admittedly crude, the waveform is only used to compute indicative source S/N and not for precise parameter estimation. Results are therefore robust against the details of the model. Indeed, we have verified that neglecting the plunge-merger-ringdown S/N does not change the projected LISA detection rates, i.e. typically the additional contribution merely increases the S/N of sources that would be detected anyway, thus potentially affecting parameter estimation – especially for the mass and spin of the remnant MBH (Klein et al., 2016) – but not the detection rate. The reason for this behaviour lies in the shape of LISA sensitivity, which is maximal for MBHB in the range $10^{4}-10^{6}~{}\rm M_{\large\odot}$ (i.e. the typical mass of LISA events), while is relatively poor for higher masses, where the signal coming from the merger and ringdown is much more relevant than that of the inspiral.

5 Results

5.1 Intrinsic MBHB merger rates

We start by discussing the intrinsic cosmic MBHB merger rates (i.e. irrespectively from LISA detection capabilities) predicted by the four MBH evolutionary models under investigation. Those are tabulated, divided in source-frame chirp mass decades, in Table 1. We limit the range to $\mathcal{M}<10^{8}~{}\rm M_{\large\odot}$ because LISA has poor sensitivity to more massive binaries, and at higher masses the rates drop to $\ll 1$ yr*-1* anyway. All results presented in the current section, unless indicated otherwise, refer to models in which we employed the DMM prescription for the gas dynamics, and $\gamma=1$ for the slope of the Denhen’s potential. Results obtained with the alternative prescriptions (i.e. HKM and $\gamma=1.75$) are quantitatively similar.

Differences between the HS and LS models have already been discussed in literature (see, e.g. Sesana et al., 2007; Sesana et al., 2011a; Klein et al., 2016) and our results are in line with previous findings. Focusing on Model-delayed, the LS scenario results in much higher merger rates ($\approx 220$ yr*-1* vs $\approx 23$ yr*-1*) with lighter chirp masses in the source frame (${\cal M}<10^{3}~{}\rm M_{\large\odot}$ vs ${\cal M}\sim 10^{5}~{}\rm M_{\large\odot}$) compared to the HS case. This is a natural consequence of the seeding mechanism, which features lighter and more abundant seeds in the LS model. The comprehensive dynamical treatment described in Section 3.1 allows us to break down the population in the different channels. The two seeding models are very different also in this respect. In the LS scenario, triplets account only for few % of the mergers, their occurrence increasing monotonically with mass. Conversely, in the HS scenario, triplets account for about 40% of the mergers in the mass range most relevant for LISA ($10^{4}~{}\rm M_{\large\odot}\lesssim{\cal M}\lesssim 10^{6}~{}\rm M_{\large\odot}$), and have less impact at higher masses. Comparing gas and stellar channels, the latter appears to be more common, reflecting a preference for gas-poor mergers in our models. For both LS and HS, this is particularly true at high masses, where dry mergers account for $\sim 80$% of the total. At lower masses, the two channels share more similar merger fractions.

Comparison between Model-delayed and Model-stalled reveals that binary stalling has a much stronger impact on the LS model, for which the merger rate drops by one order of magnitude to $\approx 20$ yr*-1*. The impact is minor on the HS model, where the rate is only halved, as expected from the fact that a large fraction of mergers was already induced by triplets in the HS Model-delayed. Also worthy of notice is the statistics of quadruplets. In Model-delayed, irrespective of the seeding scenario, quadruplets are typically more than ten times rarer than triplets; in Model-stalled, however, quadruplets can dominate the rate, especially at the high mass end. This is a natural consequence of the underlying physical implementation. The vast majority of quadruplets are indeed formed by two binaries, each coming from a different halo/galaxy. In Model-stalled a binary can merge only by encountering another systems. As cosmic time proceeds, a progressively larger fraction of MBHs is inevitably in binaries, so that mergers frequently occurs between two MBHBs, resulting in a quadruple system. This occurrence is much rarer in Model-delayed, where standard merger channels are at work. An important point to make is that even in a pessimistic (from the GW detection perspective) scenario in which all binaries stall, triplets (and quadruplets) provide an efficient channel to keep the expected MBHB cosmic merger rate at a minimum level of few yr*-1*.

Further information can be extracted from Fig. 4, where we report, for all considered models, the differential number of mergers as a function of redshift and source-frame chirp mass, expected in 4 years (the planned LISA mission lifetime). By comparing left to right panels, clearly visible is the reduction and the shift toward lower redshift of the distribution peak. This shift naturally arises from the fact that in Model-stalled two subsequent galaxy mergers are required to produce a MBHB merger. This effect is more prominent for the LS case (upper panels). The superimposed “waterfall plot" of S/N contours highlights the window that LISA will open on the cosmic population of MBHBs. Notably, LISA can essentially detect every single MBHB merger in the Universe produced by our HS models (lower panels). In the LS case (upper panels), however, low mass/high redshift systems are beyond LISA capabilities, which is also apparent from the marginalised distributions over masses and redshift, shown respectively by the top and right sub-panel of each triangle plot.

For the sake of clarity, the same marginalised distributions are separately reported in Fig. 5, together with the subdivision into different channels (i.e. star, gas, triplets) for Model-delayed, where all of them are operational. As already mentioned, in the LS scenario (upper panels) most mergers are stellar-driven and gas-driven (with a slightly predominant role played by the former), while the triplet channel only plays a minor role. The situation is very different in the HS scenario (lower panels), where instead most coalescing MBHBs are either stellar-driven or triple-driven.

5.2 LISA detection rates

We now focus on the statistics of systems detectable by LISA. Most of them will be a sub-sample of the coalescences discussed in the previous section. However, few detections will involve purely inspiralling binaries. As customary, we consider a nominal detection threshold of ${\rm S/N}=8$. We evolve the whole sample of MBHBs (equation 3) for the 4-year LISA mission lifetime, and we evaluate the S/N of each event employing the formalism described in Section 4. Results are shown in Fig. 6, where we show all sources with ${\rm S/N}>0.1$, separating MBHBs that merge within 4 years (solid lines) from the rest (dashed lines). Clearly visible is the difference between the LS (upper panel) and the HS (lower panel) cases. For LS scenarios, in fact, a large fraction of merging binaries does not reach the ${\rm S/N}=8$ threshold, while for heavy seeds essentially all MBHBs that merge during the mission are detected, as also highlighted by Fig. 4 and discussed in the previous section.

A steady MBHB merger rate has to be supplied by a continuous population of systems evenly distributed in time to coalescence. It is therefore expected that loud binaries merging in the LISA band are only the tip of a much fainter and much larger population of inspiralling systems, far from coalescence. The S/N distribution of this larger population, shown by the dashed lines in Fig. 6, decays roughly as a power-law, with a tail extending to ${\rm S/N}>8$. This has a couple of potentially interesting consequences. First, there exists a sub-population of detectable sources with coalescence timescale longer than the nominal LISA lifetime. Although rare (a few at most, in our models), these sources are persistent in the detector data stream for the whole duration of the mission. Second, the large population of inspiralling binaries with ${\rm S/N}<8$ will produce an unresolved GW background.

The number of LISA detections, assuming a 4-year mission, is listed in Table 2, divided into the different evolutionary channels, for all models. We also show the detection rates according to the specific prescriptions adopted for the stellar and gas dynamics (see Section 3.1). If binary stalling is a common occurrence in Nature (as in Model-stalled), we predict that LISA will still detect several dozens of merging MBHBs, plus few persistent sources, regardless of the nature of the seeding process. In Model-delayed, about $\sim 100$ and $\sim 300$ detections are expected in the HS and LS scenarios, respectively. Reflecting the properties of the whole population, in the LS scenario stars and gas contribute evenly to the number of detected systems, with triplets relegated to a mere 3% of detections. For the HS scenario, detected systems are more evenly distributed among the different channels (although gas is subdominant). Table 2 shows also how the detection statistics is robust against the different prescriptions considered for stellar and gas-driven dynamics. As already mentioned, we do not find significant differences with the prescription adopted.

Analogue to Fig. 5, in Fig. 7 we report, for the LS (upper panels) and HS (lower panels) scenarios, the differential number of detections as a function of redshift (left panel) and source-frame chirp mass (right panel), but considering only detectable sources (i.e. ${\rm S/N}>8$). For Model-delayed we show the usual breakup into the three sub-populations of our samples. For the HS model the number of detections essentially coincides with the number of mergers, as LISA, as already discussed, can detect almost all MBHBs in the Universe in this case. Conversely, in the LS scenario, the LISA detection efficiency is considerably worse, and the number of detections is significantly lower than the total number of mergers. This is due to the low mass, high redshift sources that contribute to the intrinsic merger rate, but which are not loud enough to be observed with LISA.

It is clear that the properties of the detected merging MBHB population can inform us about the MBH seeding mechanism that takes place in Nature (see, e.g. Sesana et al., 2011a; Klein et al., 2016). In fact, the production of a large number of mergers of low mass binaries at high redshift is a common signature of LS models, a feature that can easily be distilled from LISA data (see top vs. bottom panels in Fig. 7). On the other hand, the impact of a specific dynamical evolution mechanism (e.g. gas vs. stellar-driven) on the merger distribution is more subtle, without clear-cut distinctive features. The challenge of GW astronomy is to reconstruct the MBHB cosmic history from a limited number of detections. Countless ingredients contribute to the outcome (e.g. Berti & Volonteri, 2008; Sesana et al., 2011a; Ricarte & Natarajan, 2018b), thus pinning down the correct physics underlying the observed population will require a deep understanding of all the processes at play and a significant reverse-engineering effort. In this context, minor shifts in the redshift distribution predicted by a specific model (see blue vs. red lines in Fig. 7) are insufficient for astrophysical inference. Detection of electro-magnetic counterparts will be decisive to assess the importance of the gaseous environment for merging MBHBs. Similarly, distinctive properties of individual merging systems can provide useful information about their nature. We will now see that triplet induced mergers retain a clear signature of their dynamical origin: a non-zero (and possibly very large) eccentricity in the LISA band.

5.3 Signatures of triplet-induced mergers in LISA data

Traditionally, in the study of GW emission from MBHBs and of their observability with LISA, eccentricity has been neglected (e.g. Sesana et al., 2004; Klein et al., 2016). The rationale behind this is that GW emission causes efficient orbit circularisation (Peters, 1964) and unless a binary forms with high eccentricity, any trace of it is washed out by the time the source reaches the LISA band. By modelling the dynamical interplay between MBHBs and their stellar and gaseous environments it was soon realised that binaries can indeed be quite eccentric (e.g. Quinlan, 1996; Cuadra et al., 2009; Khan et al., 2012a), and that some residual eccentricity could still be present when a source enters the LISA band (Armitage & Natarajan, 2005; Sesana, 2010; Roedig et al., 2011; Roedig & Sesana, 2012). As we now discuss, triplets can push eccentricity evolution to the extreme, with dramatic consequences for GW parameter estimation and possibly detection.

We therefore proceed to compare the MBHB eccentricity distribution resulting from the 3 considered evolutionary channels (i.e. star, gas and triplets). The properties of triple-induced mergers are the result of secular dynamics and chaotic interactions that primarily depend on the masses of the involved systems. Conversely, the binary initial conditions and the properties of the surrounding environment play an important role in the standard merger channels, as described in Section 3.1.

Since eccentricity changes as the binary shrinks, it is customary to record its value at the last stable orbit, which is a well defined reference. The distribution of $e_{\rm ISCO}$ for all MBHBs that coalesce within 4 years in our models is reported in Fig. 8. Here, we assumed model DMM for the gas channel and $\gamma=1$ for the stellar channel. As expected, the final value of the eccentricity right before merger is well below unity in all cases. It is immediately clear, however, that triplets have a distinctive signature, with a relatively flat distribution extending up to $e_{\rm ISCO}\approx 0.1$, more than an order of magnitude larger than what is found in standard stellar- and gas-driven binary evolution.

Although, in general, MBHBs have $e_{\rm ISCO}\lesssim 0.1$, it should be noted that LISA will observe several thousands of orbital cycles well before final coalescence. The correct identification of a signal early on in the inspiral phase is of critical importance for source pre-localization and electromagnetic follow up planning. It is therefore important to estimate the eccentricity distribution when the signal enters the LISA band. To this end, during the evaluation of the S/N of each source in our sample, we record the value of the eccentricity at selected S/N thresholds. In particular, we record the properties of the sources when the nominal detection threshold of ${\rm S/N}=8$ is met. This is defined as the point of entrance of the signal in the LISA band, and the corresponding eccentricity is denoted as $e_{\rm LISA}$. This procedure yields the distributions shown in Fig. 9. We compare results from Model-stalled and Model-delayed, where for the latter we also visualise the sub-population composition. For both seeding prescriptions (LS in the upper panels and HS in the lower panels), Model-delayed results in a bimodal eccentricity distribution, with a peak at $\log_{10}e\sim-2.5$, mostly formed by stellar- and gas-driven binaries (cf. thin dashed and dotted red lines) and a plateau stretching up to $e_{\rm LISA}=1$ (cf. thick blue line and thin solid red line). This latter component is primarily generated by triplet-driven mergers, which can still retain a remarkably high eccentricity when entering the LISA band. The alternative modelling for stellar- and gas-driven binaries (i.e. $\gamma=1.75$ and HKM prescriptions) provides quantitatively very similar distributions, both in terms of peak position ($\log_{10}e\simeq-2.5$) and cut-off eccentricity ($\log_{10}e\simeq-1$), suggesting that a sizeable residual eccentricity can represent a robust signature of triplet-induced coalescences.

Finally, in Fig. 10 we report, for the triplet channel only, the differential number of events with ${\rm S/N}\geq 8$ as a function of source-frame chirp mass and $e_{\rm LISA}$. The linear scale in eccentricity highlights the bimodality of the distribution, separating those systems that have (low $e_{\rm LISA}$) and those that have not yet (high $e_{\rm LISA}$) undergone circularisation. No significant dependence on the source mass is seen but at ${\cal M}<10^{4}~{}\rm M_{\large\odot}$, where low eccentricity systems are predominant. This is due to the fact that low mass binaries (especially at high redshift) spend their early inspiral phase in the LISA sensitivity bucket, around 1 mHz. If the eccentricity at this point is too high, the GW power is shifted at higher frequencies, where the LISA sensitivity has already significantly deteriorated, ultimately resulting in a non-detection (see examples in Chen & Amaro-Seoane, 2017). The statistics of detected eccentric systems is reported in Table 3. Integrating the distribution for $e_{\rm LISA}>0.1$, we expect about 23 (36) and 20 (5) systems for Model-stalled and Model-delayed, respectively, assuming the HS (LS) scenario. This implies that, in our models, very eccentric binaries can account for between 1.5% and 50% of the total merging MBHBs detectable by LISA.

6 Discussion and Conclusions

In this paper we forecasted the expected detection statistics of MBHBs in the context of future space based gravitational wave observatories such as LISA. In particular, we assessed the impact of massive black hole triplets (and quadruplets) on the expected rate of cosmic MBH mergers.

To this end, we employed the machinery developed in the first three papers of the series (Bonetti et al., 2016, 2018b, 2018a), which combines a state of the art semi-analytic model for the joint evolution of galaxies and MBHs (Barausse, 2012) to a large library of dedicated 3-body simulations tailored to study the dynamics of MBH triplets in galactic nuclei (Bonetti et al., 2018b). In our framework, MBHBs have three possible paths to coalescence: hardening against the dense stellar background (stellar channel); efficient extraction of energy and angular momentum by a circumbinary disc (gaseous channel); interaction with a third MBH brought in by a subsequent merger (triplet channel). The efficiency of the stellar and gaseous channels is set by the central stellar density and gas content of the host remnant right after a galaxy merger, which are consistently followed along each branch of the tree. The detailed MBHB dynamical evolution depends, however, on a number of parameters that cannot be determined by the semi-analytic model, such as the $\alpha$ viscosity of the circumbinary disc or the specific slope of the power-law describing the stellar density profile, as well as on the initial binary eccentricity. We thus explored a number of scenarios employing reasonable choices of those parameters.

We considered two different models for the MBH seeding mechanisms – low mass seeds as the evolutionary endpoint of Pop III stars (LS), and high mass seeds ($\sim 10^{5}~{}\rm M_{\large\odot}$) from the collapse of massive protogalactic discs (HS) – and for each of them we explored two models for the dynamical evolution of the MBHBs – efficient gas and star driven evolution (Model-delayed) vs. inefficient gas and star driven evolution (Model-stalled)– for a total of 4 main models. In Model-delayed we considered several prescriptions for the stellar and gaseous channels, as described in Section 3.1. For each separate channel we constructed the numerical multi-dimensional differential distribution of MBHBs as a function of mass, mass ratio, redshift, (observed) orbital frequency and eccentricity. We then sampled this distribution to obtain Montecarlo realizations of the observable cosmic MBHB population. This procedure allowed us to trace all the GW emitters in the Universe, from loud merging binaries to faint low frequency systems in their early inspiral phase. We used the latest LISA sensitivity curve to study the detectability of this cosmic population and the properties of the detected systems.

For the cosmic MBHB merger and LISA detection rates, our findings are generally consistent with the literature (e.g. Sesana et al., 2007; Sesana et al., 2011a; Klein et al., 2016). Considering the fiducial Model-delayed, Pop III light seeds (LS model) result in high merger rates (hundreds per year) mostly involving low mass systems (${\cal M}<10^{3}~{}\rm M_{\large\odot}$) at high redshift (with a broad peak at $z>10$, extending to $z\approx 20$). Conversely, heavy seeds (HS model) produce an order of magnitude fewer mergers ($\approx 20$ per year), involving much more massive systems at lower redshift (see Fig. 4). As a consequence, LISA can detect almost all MBHB mergers in the Universe in this latter case, but misses more than 50% of the systems in the former. The original and more relevant results of this work concern the impact of triplet dynamics on MBHB mergers, GW emission and detection. In this respect, our main findings can be summarised as follows:

for standard assumptions (Model-delayed) all three dynamical channels (star, gas, triplets) contribute to the cosmic merger rate of MBHBs. Stars and gas are responsible for a comparable share of mergers both in the LS and in the HS scenarios, whereas the contribution of triplets is more prominent in the HS case (accounting for about 40% of the total merger rate); 2. 2.

MBH triplets represent an important “safety net” for GW observations at mHz frequencies. Model-stalled demonstrates that should all binary shrinking mechanisms fail, triple (and quadruple) interactions would still bring a sizeable number of MBHB to final coalescence (about $18\,{\rm yr}^{-1}$ and $10\,{\rm yr}^{-1}$ in models LS and HS respectively); 3. 3.

even in Model-delayed, where the fiducial recipes for gas and star driven hardening are employed, triplets are responsible for a non negligible fraction of MBHB coalescences. This is only about 3% (rate of $\approx 2\,{\rm yr}^{-1}$) for the LS model and about 40% (rate of $\approx 10\,{\rm yr}^{-1}$) for the HS model. The different contribution of triplets can be explained by keeping in mind that their formation requires the existence of long-lived MBHBs, and that the greater the binary mass, the longer its lifetime. In the LS scenario, most BH mergers involve low mass objects that still approximately have their original mass. If a higher BH mass is found, then this means that the MBH has experienced several mergers or, most likely, that the additional mass has been accreted. An abundant reservoir of gas usually determines a strong coupling with the environment and, hence, a faster evolution, therefore suppressing the formation of triplets. On the contrary, in the HS scenario, MBHBs start with a higher mass ($\sim 10^{5}~{}\rm M_{\large\odot}$). This seems to be the sweet spot for triplet-induced mergers, since these binaries are moderately long-lived (as the nuclear gas reservoir is typically negligible compared to these more massive binaries). At the same time, there exists a large population of MBHs with roughly the same mass and thus capable of perturbing the “stalled” binaries. This fact also explains why at higher MBH masses the contribution of triplet-induced mergers decreases. Indeed, as the MBHB mass grows, it is less likely that a possible perturber is massive enough to trigger the coalescence; 4. 4.

triplet-induced MBHB inspirals retain remarkably high eccentricities in the LISA band. In all scenarios considered, about 50% of these systems have an eccentricity $e_{\rm LISA}>0.1$ when entering the LISA band, and eccentricities as large as 0.9 are the norm rather than the exception. Conversely, such high eccentricities are rare among systems evolving via standard shrinking channels (gas and stellar-driven).

Item (ii), in particular, is of paramount importance for future GW observations with LISA. Most of the LISA science case is indeed based on MBHB merger rates derived from semi-analytic models (see e.g. the seminal study of Volonteri et al., 2003) with a number of simplifying assumptions, especially concerning the treatment of the MBHB dynamics . On the numerical side, zoom-in high resolution hydrodynamical simulations have been instrumental in understanding the large scale dynamics leading to the pairing of the binaries, but the parsec and sub-parsec evolution still relies on sub grid recipes (e.g, Bellovary et al., 2013; Tremmel et al., 2018). Large cosmological simulations, on the other hand, simply do not have yet the resolution to catch the bulk of the expected LISA sources at masses below $10^{6}~{}\rm M_{\large\odot}$, besides suffering of even more severe shortcomings in the treatment of the MBHB dynamics (see e.g. Vogelsberger et al., 2014; Schaye et al., 2015, and references therein). Sub-grid recipes are based on analytic modelling and numerical simulations of MBHBs in stellar and gaseous discs. Although the literature on the subject is vast (Berczik et al., 2006; Preto et al., 2011; Khan et al., 2011; Khan et al., 2012a; Vasiliev et al., 2015; Sesana & Khan, 2015; Bortolas et al., 2016; Milosavljević & Phinney, 2005; Cuadra et al., 2009; Roedig et al., 2011; Roedig & Sesana, 2012; Shi et al., 2012; Farris et al., 2014; Tang et al., 2018) and the main physical mechanisms at play are generally understood, the evolution outcomes depend on poorly known initial conditions (stellar density, relaxation mechanisms, properties of the gas feeding) and, in the case of gas, on (micro)physics that is hard to control (e.g. cooling, viscosity, AGN and supernova feedback). Star driven evolution can take up to Gyrs (Vasiliev et al., 2015; Sesana & Khan, 2015; Gualandris et al., 2017), whereas the efficiency of coalescence mediated by gas is hard to predict (Lodato et al., 2009; Goicovic et al., 2018). A legitimate question is therefore what happens if MBHB hardening mechanisms fail altogether. Our results show that triple MBH interactions are an efficient pathway to bring a sizeable number of MBHBs to coalescence in a cosmological context, thus minimising the risk posed by binary stalling to low frequency GW astronomy.

Item (iv) has important implications both for GW astrophysics and signal detection. The expected occurrence of binaries entering in the LISA band with eccentricities as high as 0.9 calls for the development of faithful eccentric waveforms. For such high eccentricities, detection with circular templates may still be feasible, but will be sub-optimal, causing a significant loss in accumulated S/N. This may ultimately delay detection and may severely affect the estimation of source parameters, jeopardising the chances of success of electromagnetic follow-ups. The development of eccentric waveform is ongoing (Key & Cornish, 2011; Nishizawa et al., 2016; Huerta et al., 2017; Klein et al., 2018), but more work is needed in this direction. In terms of GW astrophysics, the detection of several highly eccentric binaries might be the smoking gun of a substantial population of MBH triplets. In general, standard hardening channels do not produce eccentricities larger than 0.1 in the LISA band (Roedig & Sesana, 2012), although for specific conditions this might be possible. For example, MBHBs hardening in counter rotating stellar systems as well as in counter rotating discs may have extremely high eccentricities (Nixon et al., 2011; Sesana et al., 2011b; Holley-Bockelmann & Khan, 2015; Schnittman & Krolik, 2015). It is however unclear how relevant these physical configurations are, especially for LISA MBHBs, since the cumulative effect of star-binary interactions/viscous torques tend to align the the orbital MBHB angular momentum to that of the stellar distribution/disk (Gualandris et al., 2012; Dotti et al., 2013; Roedig & Sesana, 2014; Rasskazov & Merritt, 2017).

Besides MBHBs merging within the LISA mission lifetime, we have also simulated the much larger population of inspiralling systems still far from coalescence. We found that detection of few persistent MBHBs with ${\rm S/N}>8$ is expected. Moreover, the collective signal of the vast population of sub-threshold systems might give rise to an interesting unresolved background (Sesana et al., 2007). We plan to study in detail the statistical properties and emitted GW signal of this persistent MBHB population in a forthcoming publication.

Acknowledgements

MB, MC and FH acknowledge partial financial support from the INFN TEONGRAV specific initiative. MB acknowledges CINECA, under the TEONGRAV initiative, for the availability of high performance computing resources and support. This work was supported by the H2020-MSCA-RISE-2015 Grant No. StronGrHEP-690904. This project has received funding (to E. Barausse) from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. GRAMS-815673; project title “GRavity from Astrophysical to Microscopic Scales”). AS is supported by a University Research Fellowship of the Royal Society. The authors would also like to acknowledge networking support by the COST Action CA16104. 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.