<figure><img src="image_1.png" /></figure>

# Massive Black Hole Merger Rates: The Effect of Kiloparsec Separation Wandering and Supernova Feedback

Enrico Barausse<sup>1,2,3 Ⓓ</sup>, Irina Dvorkin<sup>3</sup>, Michael Tremmel<sup>4 Ⓓ</sup>, Marta Volonteri<sup>3 Ⓓ</sup>, and Matteo Bonetti<sup>5 Ⓓ</sup>

<sup>1</sup>SIŞSA, Via Bonomea 265, 34136 Trieste, Italy and INFN Sezione di Trieste, Italy; barausse@sissa.it

<sup>2</sup>IFPU—Institute for Fundamental Physics of the Universe, Via Beirut 2, I-34014 Trieste, Italy

<sup>3</sup>Institut d'Astrophysique de Paris, CNRS & Sorbonne Université, UMR 7095, 98 bis bd Arago, F-75014 Paris, France

<sup>4</sup>Yale Center for Astronomy and Astrophysics, Physics Department, P.O. Box 208120, New Haven, CT 06520, USA

<sup>5</sup>Dipartimento di Fisica “G. Occhialini,” Università degli Studi di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy

Received 2020 June 4; revised 2020 August 12; accepted 2020 September 20; published 2020 November 18

## Abstract

We revisit the predictions for the merger rate of massive black hole binaries detectable by the Laser Interferometer Space Antenna (LISA) and their background signal for pulsar-timing arrays. We focus on the effect of the delays between the merger of galaxies and the final coalescence of black hole binaries, and on supernova feedback on the black hole growth. By utilizing a semianalytic galaxy formation model, not only do we account for the driving of the evolution of binaries at separations ≲1 pc (gas-driven migration, stellar hardening, and triple/quadruple massive black hole systems), but we also improve on previous studies by accounting for the time spent by black hole pairs from kiloparsec down to parsec separation. We also include the effect of supernova feedback, which may eject gas from the nuclear region of low-mass galaxies, thus hampering the growth of black holes via accretion and suppressing their orbital migration in circumbinary disks. Despite including these novel physical effects, we predict that the LISA detection rate should still be ≳2yr⁻¹, irrespective of the model for the black hole seeds at high redshifts. Scenarios where black holes form from ~100M☉ seeds are more significantly impacted by supernova feedback. We also find that for detectable events, the merging black holes typically have mass ratios between ~0.1 and 1. Predictions for the stochastic background in the band of pulsar-timing array experiments are instead rather robust and show only a mild dependence on the model.

*Unified Astronomy Thesaurus concepts: Astrophysical black holes (98); Gravitational waves (678); Galaxy formation (595); Supermassive black holes (1663)*

## 1. Introduction

The origins of massive black holes (MBHs) and the nature of their coevolution with their host galaxies remain fundamental questions in astrophysics, which current and future gravitational wave detectors may help decipher. MBHs are ubiquitous in massive galaxies (Gehren et al. 1984; Kormendy & Richstone 1995) as well as in a fraction of low-mass dwarf galaxies (Reines et al. 2011, 2013; Baldassare et al. 2020). Feedback from active galactic nuclei (AGNs), powered by growing MBHs, is commonly thought to regulate and quench star formation in massive galaxies (Croton et al. 2006) and possibly in dwarf galaxies as well (Dickey et al. 2019; Sharma et al. 2020). Despite their importance to galaxy formation theory, the mechanisms that drive and regulate MBH growth in galaxies, as well as the physical processes surrounding AGN feedback, are still not well understood. Scaling relations between MBH mass and galaxy properties are indicative of coeval growth (Kormendy & Ho 2013; McConnell & Ma 2013; Schramm & Silverman 2013), although as the census of MBHs in the local universe improves, such relationships are found to be more complicated than previously thought (Volonteri & Reines 2016; Shankar et al. 2016; Barausse et al. 2017; Greene et al. 2019).

An important limitation of our understanding of MBHs and their effect on galaxies is related to the difficulty of studying their early evolution at high redshift, as only the most luminous, rapidly growing MBHs are accessible by observations. Gravitational waves represent an intriguing window into the history of MBHs, because their propagation through the universe, unlike that of electromagnetic radiation, is essentially

unobstructed. The future Laser Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2017) will be able to detect gravitational waves emitted by merging MBHs with masses $10^4–10^7 M☉$ out to redshifts greater than 20, while providing accurate estimates on MBH mass, spin, and MBH binary orbital parameters (Klein et al. 2016; Colpi et al. 2019). Such detections promise new insight into the MBH population at early times and may place unique constraints on models of MBH formation and growth (Sesana et al. 2007; Berti & Volonteri 2008; Volonteri & Natarajan 2009; Plowman et al. 2011; Sesana et al. 2011; Klein et al. 2016; Ricarte & Natarajan 2018a; Bonetti et al. 2019).

Currently, pulsar-timing arrays are searching for (unresolved) gravitational wave signals from lower redshift ($z < 2$), higher mass ($>10^8 M☉$) MBH binaries using Milky Way pulsars (Lentati et al. 2015; Shannon et al. 2015; Arzoumanian et al. 2016; Verbiest et al. 2016; Arzoumanian et al. 2018; Perera et al. 2019). While the sensitivity of pulsar-timing arrays is still increasing with time, the absence of a detection so far and the resulting upper limits on the background of unresolved gravitational waves at nanohertz frequencies are already placing significant constraints on models of MBH mergers (Rajagopal & Romani 1995; Jaffe & Backer 2003; Wyithe & Loeb 2003; Sesana et al. 2008, 2009; Ravi et al. 2012; Sesana 2013; McWilliams et al. 2014; Kulier et al. 2015; Ravi et al. 2015; Sesana et al. 2016; Kelley et al. 2017; Bonetti et al. 2019), even ruling out the most extreme ones (McWilliams et al. 2014). However, it is critical to note that much of the astrophysics that can be inferred from gravitational wave detectors relies heavily on the models used to predict the MBH
binary and merger populations and their various underpinning assumptions.

Semianalytic models of galaxy and MBH formation and evolution are powerful tools for deriving astrophysics from gravitational wave detections, as well as informing the experimental setups themselves (e.g., Sesana et al. 2004; Volonteri et al. 2008; Volonteri & Natarajan 2009; Barausse 2012; Klein et al. 2016; Ricarte & Natarajan 2018b; Bonetti et al. 2019). These models relate the hierarchical formation of dark matter halos to the evolution of galaxies and their MBHs. Semianalytic models have been successful in reproducing the observed evolution of galaxy morphology, color, star formation rate, and luminosity functions and, because they are relatively inexpensive, have been used to explore the wide parameter space of galaxy formation physics (e.g., Somerville & Primack 1999; Croton et al. 2006; Somerville et al. 2008). Critically, modern semianalytic galaxy formation models include prescriptions for MBH growth and feedback, similarly constrained by both galaxy observations and quasar luminosity functions (Somerville et al. 2008; Barausse 2012; Sesana et al. 2014; Ricarte & Natarajan 2018b). By tracking the properties of merging galaxies and the formation of MBH binaries over cosmic time, these simulations can predict how gravitational wave signals relevant to both LISA and pulsar-timing arrays are affected by various physical processes, such as feedback and MBH dynamical evolution on various scales.

The process of forming an MBH binary begins when two galaxies embedded in their dark matter halos merge, and the MBHs are still separated by tens to hundreds of kiloparsecs. The MBH orbital evolution then proceeds down to separations of a few hundred parsecs (Tremmel et al. 2018, 2018). However, in order for two MBHs to merge, they must evolve down to separations of ~0.001–0.01 pc, where emission of gravitational radiation can bring the MBHs to coalescence. When the MBHs form a binary, on ~1–10 pc scales, the dynamical evolution is facilitated by complicated interactions with the binary’s stellar and gaseous environment (Quinlan 1996; MacFadyen & Milosavljević 2008; Cuadra et al. 2009; Lodato et al. 2009; Roedig et al. 2011; Nixon et al. 2011; Sesana & Khan 2015; Vasiliev et al. 2015; Duffell et al. 2020; Muñoz et al. 2019), as well as with other MBHs through three-/four-body interactions (Bonetti et al. 2018b, 2019). The exact evolutionary channel of an MBH binary is partly determined by the morphology and kinematic structure of their host galaxies (Khan et al. 2013; Holley-Bockelmann & Khan 2015). Each of these dynamical processes effectively results in a delay between the merger of two dark matter halos/galaxies and the merger of two MBHs (and the resulting gravitational wave emission). For the larger-scale evolution, semianalytic models are informed by results from cosmological N-body and hydrodynamic simulations, while the evolution of MBH binaries is informed by detailed binary evolution simulations.

Indeed, because of their importance to galaxy formation, MBHs have become an integral part of most large-scale cosmological hydrodynamic simulations, which have been able to reproduce MBH and galaxy scaling relations, as well as implement feedback from MBHs that can successfully regulate and quench star formation in massive galaxies (Di Matteo et al. 2005; Dubois et al. 2012; Vogelsberger et al. 2014; Schaye et al. 2015; Volonteri et al. 2016; Tremmel et al. 2017; Pontzen et al. 2017; Nelson et al. 2019; Ricarte et al. 2019). However, the poor resolution of these simulations requires relatively

simplistic models for MBH formation, such that the earliest phases of growth and mergers often remain unresolved. Because of the complicated physics involved in fully hydrodynamic simulations, they are also generally limited to relatively small volumes with poor statistics. Most importantly, the dynamical evolution of MBHs is often completely ignored even on resolved ~kiloparsec scales, though there are important exceptions to this (Dubois et al. 2012; Tremmel et al. 2017; Hirschmann et al. 2014; Pfister et al. 2019). While some cosmological simulations have been used, with significant postprocessing, to account for unresolved MBH binary evolution and make unique predictions relevant for gravitational wave astrophysics (e.g., Blecha et al. 2016; Kelley et al. 2017; Katz et al. 2019), they can also be useful tools for improving semianalytic simulations, which, because of their low computational cost, provide better statistics and the ability to test the effects of different model assumptions.

Recently, cosmological hydrodynamic simulations have seen important improvements to how MBHs are modeled, which has led to new results regarding their growth and dynamical evolution. In the ROMULUS simulations (Tremmel et al. 2017), the dynamical evolution of MBHs is self-consistently followed by using a new, subgrid model accounting for unresolved dynamical friction (Tremmel et al. 2015). Using these simulations, Tremmel et al. (2018) showed that the formation of MBH pairs with separations below a kiloparsec (the precursors to MBH binaries) often occurs after several gigayears of orbital evolution of MBH pairs at kiloparsec-scale separations, with many MBH pairs failing to ever form a binary within a Hubble time (Tremmel et al. 2018). Furthermore, the inspiral timescale may also be influenced by (sub)kiloparsec-scale structures, such as galactic bars, spirals, and massive star-forming clumps. These might render the orbital decay significantly stochastic, thus enriching the dynamics of MBHs on kiloparsec scales (even leading to the possible ejection of infalling MBHs; see, e.g., the recent hydrodynamic simulations of Tamburello et al. 2017; Bortolas et al. 2020).

Detailed, high-resolution cosmological simulations of MBH formation at high redshift have also shown that supernova (SN) feedback can stunt early MBH growth in low-mass galaxies (Dubois et al. 2015; Habouzit et al. 2017). In this paper, we combine these new results from large-scale simulations with updated models for MBH binary evolution (e.g., Sesana et al. 2014; Antonini et al. 2015a; Bonetti et al. 2018b, 2019) into the semianalytic model for MBH and galaxy evolution of Barausse (2012). The goal of this paper is to understand how the combination of different assumptions for MBH formation, early growth, and dynamical evolution affects the gravitational wave signal expected for LISA and pulsar-timing arrays.

In Section 2, we summarize the semianalytic model and its physical ingredients, as well as the different model variations that we run. Our results are presented in Sections 3–6, including MBH merger rates and predicted mass, redshift, and mass ratio distributions for events detectable by LISA, quasar luminosity functions, as well as predictions for pulsar-timing arrays. In Section 7, we summarize our results and draw our conclusions.

## 2. The Semianalytic Model

We describe the synergistic coevolution of MBH and their host galaxies using the semianalytic model of Barausse (2012), with later updates to specific aspects of it described in Sesana
et al. (2014), Antonini et al. (2015a), and Bonetti et al. (2018b, 2019). In the following, we limit ourselves to a brief review of the model, referring the reader to the aforementioned works for more details, and highlighting the changes with respect to them. Besides slight changes in the gas cooling, star formation, and AGN feedback prescriptions, which we describe in this section and against which the results are robust, the improvements to the model on which this paper hinges are instead described in Sections 2.1–2.4.

The model is built on top of a dark matter merger tree constructed with an extended Press–Schechter algorithm (Press & Schechter 1974; Parkinson et al. 2008) suitably modified to reproduce the results of N-body simulations (Parkinson et al. 2008). Baryonic structures are then evolved along this dark matter merger tree using semianalytic prescriptions. These structures include a chemically pristine intergalactic medium, which accretes onto dark matter halos and gets shock-heated to their virial temperature in high-mass systems at low redshift, or flows into halos on their dynamical time along cold filaments (Cattaneo et al. 2006; Dekel & Birnboim 2006; Dekel et al. 2009). Unlike in Barausse (2012) and other works based on it, we allow here for a smooth transition between the shock-heating and cold-filament regimes by using the results of Correa et al. (2018), based on hydrodynamic cosmological simulations from the EAGLE project, and particularly their Equations (11)–(15) for the fraction of hot mode gas accretion.

Cooling of the intergalactic medium and/or its inflow along cold streams then gives origin to a cold gas phase (an “interstellar medium”) that eventually undergoes star formation. In more detail, the model tracks the evolution of gaseous/stellar disks and spheroids, with major galaxy mergers and bar instabilities disrupting disks and turning them into spheroids. The star formation is described in spheroids via the prescriptions of Barausse (2012) and in disks via those of Dutton & van den Bosch (2009), who assume that star formation takes place in dense molecular clouds.<fnref n="6" /> Star formation and SN feedback also drive the chemical evolution of the interstellar medium.

On smaller scales, the model accounts for the formation of nuclear star clusters (from in situ star formation and/or from migration of globular clusters to galactic nuclei; see Antonini et al. 2015b, 2015a) and the presence of MBHs. MBHs form from high-redshift seeds and then grow by accretion of nuclear gas and coalescences with other MBHs brought by galaxy mergers. The mass growth rate of the reservoir of nuclear gas available within the MBH influence radius for accretion and for in situ nuclear star cluster formation is assumed to be linearly correlated with star formation in the spheroid (Granato et al. 2004; Lapi et al. 2014; Ricarte et al. 2019), i.e., we assume $\dot{M}_{\text{nucl}} = A\dot{M}_{\text{sf,sph}}$, with the model’s calibration against local and high-redshift observables (see Barausse 2012; Sesana et al. 2014; Antonini et al. 2015b, 2015a; and Section 5 for a list) yielding $A = 5 \times 10^{-2}$. Accretion of this nuclear gas onto the MBH is assumed to take place on the viscous timescale evaluated at the influence radius (Sesana et al. 2014), but is topped off at $\dot{M}_{\text{bh,max}} = A_{\text{Edd}}\dot{M}_{\text{Edd}}$, with $\dot{M}_{\text{Edd}}$ the Eddington rate, and $A_{\text{Edd}} = 1$ or 2 depending on the seeding model (see below). The MBH evolution also exerts a feedback on the growth of structures (AGN feedback). While Barausse (2012), and later work based on it, only accounted for AGN feedback by radio jets,

in this paper we also consider the effect of radiative feedback (Bieri et al. 2017). We model this effect by injecting 5% of the AGN luminosity into the surrounding gas and then computing the feedback onto the bulge and diffuse gas distribution via Equations (42) and (43) of Barausse (2012) (where we therefore replace $L_{\text{jet,QSO}}$ with $L_{\text{jet,QSO}} + 0.05L_{\text{bh,QSO}}$).

The most crucial aspects of our model for the prediction of LISA event rates and the amplitude of the pulsar-timing array stochastic backgrounds are the seeding mechanism of MBHs and the timescales on which MBHs coalesce after their host galaxies merge (Sesana et al. 2007; Klein et al. 2016; Bonetti et al. 2019). Several physical models for the mass function of MBH seeds at high redshift have been put forward; see, e.g., Latif & Ferrara (2016) for a review. Here, we consider two representative scenarios, namely a light-seed (LS) model in which seeds are provided by the remnants of Population III stars forming in high-redshift, low-metallicity environments (Madau & Rees 2001), and a heavy-seed (HS) model where the seeds are instead formed by the collapse of protogalactic disks following the onset of bar instabilities (Volonteri et al. 2008).

In the LS scenario, we account for the mass losses during stellar collapse by assuming that seeds have mass $\sim 2/3$ of the mass of the initial Population III star, which we draw randomly from a log-normal distribution function centered on $300M_{\odot}$, with a standard deviation of 0.2 dex and a gap between 140 and 260 $M_{\odot}$ (to account for the fact that pair-instability SNe are not believed to leave a black hole remnant; see Heger & Woosley 2002). Moreover, following Volonteri et al. (2003), in this scenario we only seed the most massive halos, forming from the collapse of the $3.5\sigma$ peaks of the matter density field at high redshifts $z \gtrsim 15$. To ease the possible discrepancy between LS models and the high-redshift luminosity function (Madau et al. 2014), we allow for mildly super-Eddington MBH accretion ($A_{\text{Edd}} \approx 2$) in the LS scenario.

For the HS scenario, we adopt the model by Volonteri et al. (2008), where protogalactic disks are assumed to become unstable when their Toomre parameter gets lower than a critical threshold $2 \lesssim Q_c \lesssim 3$. The resulting seeds have mass related to the properties of the host halos (see Volonteri et al. 2008 for details), but typically of the order of $\sim 10^5 M_{\odot}$. Following Volonteri et al. (2008), also in this scenario we only seed halos at $z \gtrsim 15$, but we set $A_{\text{Edd}} = 1$ (unlike for LSs). Note also that in this paper, we fix $Q_c = 3$, unlike Bonetti et al. (2019), where $Q_c = 2.5$ was used. A higher value of $Q_c$ increases the number of seeds; therefore, in this paper, we explore a case where HSs form abundantly. This has to be taken into account when comparing this paper’s results to previous studies using the same semianalytic model and also to other semianalytic models (Sesana et al. 2007; Ricarte & Natarajan 2018a; Dayal et al. 2019), where HSs are normally rare. Note that Klein et al. (2016) assume $Q_c = 3$, as in this paper. However, the model for the evolution of MBH pairs in Klein et al. (2016) was simplified with respect to this work, especially when it comes to triple/quadruple MBH interactions (see Bonetti et al. 2019) and to the model for the intermediate-scale dynamical “delays” between galaxy/halo mergers and MBH coalescences (see below).

When it comes to the delays between galaxy and MBH mergers, Barausse (2012), Sesana et al. (2014), Antonini et al. (2015a), and Bonetti et al. (2019) employed a rather sophisticated model for the halo merger times, as well as for the evolution of MBH binaries on scales below a parsec (where stellar hardening, gas-driven migration, formation of triple/

<footnote n="6">⁶ Note that these star-formation prescriptions are slightly different from those used in previous works based on Barausse (2012)—such as Sesana et al. (2014), Antonini et al. (2015a), and Bonetti et al. (2018b, 2019)—but our results are robust against these changes.</footnote>
quadruple MBH systems, and gravitational wave emission were all accounted for), but lacked a description for the dynamics of MBH pairs on scales from several kiloparsecs down to a few parsecs (i.e., the time between the onset of the galaxy interaction and the formation of a bound MBH binary).

In more detail, the model of Barausse (2012), Sesana et al. (2014), Antonini et al. (2015a), and Bonetti et al. (2019) identifies the “nodes” of the underlying dark matter merger tree with the epoch at which the smaller (“satellite”) halo first enters the host halo. The satellite halo is then assumed to survive within the host as a subhalo while losing mass because of tidal stresses/evaporation and sinking toward the center under the effect of dynamical friction. We model the mass loss due to tidal effects by following Taffoni et al. (2003), while for the dynamical-friction time, we adopt the fit (to N-body simulations) of Boylan-Kolchin et al. (2008). Only once the dynamical friction time has elapsed is the subhalo assumed to lose its individual identity. At that point, Barausse (2012), Sesana et al. (2014), Antonini et al. (2015a), and Bonetti et al. (2019) assume that the satellite and host galaxies merge and that their MBHs, if present, form a binary.

The evolution of MBH binaries was then followed by accounting (via semianalytic prescriptions) for the three-body interactions with stars (Quinlan 1996; Sesana & Khan 2015), which extract energy and angular momentum from the binary. This results in a slow secular “hardening” on timescales of ~gigayears from separations comparable to the hardening radius $a_h = Gm_2/(4\sigma^2)$ (where $m_2$ is the secondary MBH’s mass; see Quinlan 1996) down to the separation (typically $10^{-2} - 10^{-3}$ pc) at which gravitational-wave emission is sufficient to drive the MBHs to merger in less than a Hubble time. In gas-rich galactic nuclei, where the nuclear gas mass contained within the influence radius of the binary exceeds the total binary mass, the model assumes instead that the binary’s evolution is driven by the gas on its viscous timescale ~$10^7 - 10^8$ yr, again down to the separation at which gravitational-wave emission becomes dominant. Finally, first in Antonini et al. (2015a) and Klein et al. (2016) and then, in a more rigorous fashion, in Bonetti et al. (2018b, 2019), the model accounts for the possible formation of triple/quadruple MBH systems, which can form when two galaxies merge before the hosted MBH binary (or binaries) has had time to coalesce. These MBH triple/quadruple systems can trigger the merger of at least two MBHs via hierarchical Kozai–Lidov interactions (Kozai 1962; Lidov 1962) or chaotic three-body interactions (Bonetti et al. 2018a).

## 2.1. Delays due to the Dynamical Evolution of Merging Galaxies

The timescale of Boylan-Kolchin et al. (2008), which as mentioned above we use to model the dynamical friction from the satellite subhalo within its host, correctly accounts for the evolution of the dark matter in a galaxy merger, but may be too short to account also for the subsequent dynamics of the baryonic components (and particularly the MBHs). We therefore augment it by several additive terms, accounting for the dynamical friction between galaxies, and later for the dynamical friction exerted on the individual MBHs (naked or still surrounded by a core of stars) by the stellar background, down to the hardening radius.

In more detail, on the scale of interacting galaxy pairs, we describe the dynamical friction timescale of two galaxies by following Binney & Tremaine (2008):

$$T_{\text{df}} = \frac{2.7 \text{ Gyr}}{\ln\Lambda} \frac{R_i}{30 \text{ kpc}} \left( \frac{\sigma_h}{200 \text{ km s}^{-1}} \right)^2 \left( \frac{100 \text{ km s}^{-1}}{\sigma_s} \right)^3, \quad (1)$$

$$\Lambda^{-1} = \max \left[ \frac{\sigma_s}{2^{3/2}\sigma_h}, \sqrt{2} \left( \frac{\sigma_s}{\sigma_h} \right)^3 \right], \quad (2)$$

where $\sigma_s$ and $\sigma_h$ are the velocity dispersions of the satellite and host galaxies, while $R_i$ is the initial separation, which we set to the half-light radius of the galaxy that will form from the merger (see Section 2.2.2 of Barausse 2012). Note that this expression accounts for the progressive tidal stripping of the secondary galaxy along the evolution (Binney & Tremaine 2008, at least as long as both galaxies are modeled simplistically as isothermal spheres). However, when considering the earlier dynamical-friction-driven evolution of the satellite subhalo within its host, if the satellite galaxy has already been completely disrupted by tidal stripping and evaporation (which we describe by the model of Taffoni et al. 2003), then applying Equations (1)–(2) is not physically meaningful. For those systems, we therefore replace Equations (1)–(2) by the dynamical friction timescale of a “naked” MBH (Binney & Tremaine 2008):

$$T_{\text{df}} = \frac{19 \text{ Gyr}}{\ln[R_i \sigma_h^2/(GM_{\text{bh,sat}})]} \left(\frac{R_i}{5 \text{ kpc}}\right)^2 \left(\frac{\sigma_h}{200 \text{ km s}^{-1}}\right) \left(\frac{10^8 M_{\odot}}{M_{\text{bh,sat}}}\right), \quad (3)$$

where again $R_i$ is set to the half-light radius of the host galaxy (because the satellite galaxy has been completely destroyed).

## 2.2. Delays due to Dynamical Evolution of MBH Pairs on Kiloparsec Scales

To account for the possibility that simplified prescriptions such as Equations (1)–(3) may not be sufficiently realistic on smaller scales, we also define an additional timescale $T_{\text{Romulus}}$ designed to fit the times spent by MBH binaries at separations at hundreds of parsecs (or larger) in the ROMULUS simulation of Tremmel et al. (2017). In more detail, to compute $T_{\text{Romulus}}$ for each galaxy merger, we sample a probability distribution function

$$\frac{dp}{d\log T_{\text{Romulus}}} = \frac{F}{\sqrt{2\pi}\sigma} e^{-\frac{(\log T_{\text{Romulus}} - \log \mu)^2}{2\sigma^2}} + (1-F)\delta(\log T_{\text{Romulus}} - \log t_0), \quad (4)$$

which is the linear combination of a log-normal distribution centered on $\mu$ and with standard deviation $\sigma$, and a Dirac delta peaked at the age of the universe $t_0 = 13$ Gyr. This bimodal distribution models the fact that only a fraction $F < 1$ of galaxy mergers results in a close MBH pair in the simulation of Tremmel et al. (2018), as shown in their Figure 1 (from which we extract $F$ as function of the stellar mass of the primary galaxy). The times spent by close MBH pairs at separations of
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Number of total (thin lines) and detected (thick lines) MBH mergers per unit redshift during a 4 yr LISA mission in the HS and LS models that include SN feedback (upper and lower panels, respectively; see Table 1 for a summary of the models). The difference between total and detected mergers for HS models is very small and not visible in this plot.</figcaption></figure>

hundreds of parsecs are instead shown in Figure 6 of Tremmel et al. (2018) as the differences between the orange and blue points. We fit these differences with a log-normal distribution, which yields $\mu = 1.3$ Gyr and $\sigma = 0.7$ (dex).

In order to avoid double counting the galaxy merger timescale naturally present in the cosmological simulations, we assume that the time spent by MBH binaries at separations of 100 s–1000 s pc to be max($T_{\text{df}}$, $T_{\text{Romulus}}$), where $T_{\text{df}}$ is the galaxy dynamical friction timescale discussed in the previous section.<fnref n="7" />

## 2.3. Dynamical Friction on Bound MBH Pairs (below 100 pc Scales)

At separations comparable with the influence radius of the primary MBH, the MBHs become bound in a binary, which changes the details of the dynamical friction (Dosopoulou & Antonini 2017). The dynamical friction timescale for the binary’s evolution from the primary’s influence radius $r_{\text{infl}}$ down to a smaller separation $\chi r_{\text{infl}}$ is given by (Dosopoulou & Antonini 2017)

$$ T_{\text{bare}}^{\text{df,infl}} = 1.5 \times 10^7 \frac{[\ln\Lambda'\alpha + \beta + \delta]^{-1}}{(3/2 - \gamma)(3 - \gamma)} (\chi^{\gamma-3/2} - 1) \left(\frac{M_1}{3 \times 10^9 M_{\odot}}\right)^{1/2} \left(\frac{M_2}{10^8 M_{\odot}}\right)^{-1} \left(\frac{r_{\text{infl}}}{300 \text{ pc}}\right)^{3/2} \text{yr, (5)} $$

where $M_1$ and $M_2$ are the primary and secondary MBH masses, and we assume for simplicity $\Lambda' \approx r_{\text{infl}}\sigma^2/(GM_2) \approx M_1/M_2$. The coefficients $\alpha$, $\beta$, and $\delta$ are functions of the power-law exponent $\gamma$ of the stellar density near $r_{\text{infl}}$, i.e., $\rho_* \propto (r/r_{\text{infl}})^{-\gamma}$, and are given by Equations (21)–(23) of Dosopoulou & Antonini (2017). Here, we assume $\gamma = 1$.

Because Equation (5) does not account for the fact that some stellar mass from the satellite galaxy (if it has not been completely disrupted by tidal effects earlier on) can remain bound to the secondary MBH even within $r_{\text{infl}}$, we follow again Dosopoulou & Antonini (2017) and model this effect by the timescale

$$ T_{\text{dressed}}^{\text{df,infl}} = 1.2 \times 10^7 \frac{[\ln\Lambda\alpha + \beta + \delta]^{-1}}{(3 - \gamma)^2} (\chi^{\gamma-3} - 1) \left(\frac{M_1}{3 \times 10^9 M_{\odot}}\right) \left(\frac{100 \text{ km s}^{-1}}{\sigma_s}\right)^3 \text{yr, (6)} $$

where $\Lambda$ is given by Equation (2). We therefore assume that the decay timescale from the influence radius is given by

$$ T^{\text{df,infl}} = \min(T_{\text{bare}}^{\text{df,infl}}, T_{\text{dressed}}^{\text{df,infl}}). \quad (7) $$

## 2.4. Effect of SN Feedback on MBH Growth

Another ingredient that we add to our semianalytic model is the possibility that SN feedback may stunt the growth of MBH in low-mass galaxies. Indeed, while the effect of SN feedback on star formation is already included in the model of Barausse (2012), its effect on the MBH growth had not been included yet. In fact, for MBH growth, the relevant quantity is not the overall gas fraction in the galaxy, but the physical state (density, temperature) of the gas near the MBH and how this gas is distributed (Dubois et al. 2015).

The exact way SN explosions affect gas, and therefore whether they are able to evacuate the gas near MBHs, strongly depends on the details of the process, i.e., on how the energy released in the explosion couples to the gas. For instance, Habouzit et al. (2017) find that weak thermal and kinetic SN feedback do not have a dramatic effect on the MBH growth, while in SN feedback models where gas cooling is delayed because of the (unresolved) shocks, accretion onto MBHs is suppressed in low-mass systems. Importantly, the observational properties of galaxies and MBHs are well reproduced in the simulations of Habouzit et al. (2017) only with the latter implementation of SN feedback. To account for the possible suppression of MBH growth in low-mass galaxies caused by SNe, we assume that the growth rate $\dot{M}_{\text{nucl}}$ of the nuclear gas reservoir from which the MBH accretes is quenched in galaxies where the escape velocity from the spheroidal bulge is lower than 270 km s⁻¹. The latter is indeed the speed of SN winds in the delayed cooling simulations of Habouzit et al. (2017). This suppresses MBH growth and lengthens the viscous timescale for the evolution of MBH binaries in circumbinary disks.

## 3. Catalogs of Merging MBH Binaries

Based on the output of our semianalytic model, we compute the expected detection rate of merging MBH binaries by the LISA mission. In more detail, we consider the models in Table 1, and for each of them, we produce synthetic catalogs of merging MBH binaries, including all relevant information on each binary, such as the component masses, merger redshift

<footnote n="7">⁷ The reader may think that because the ROMULUS simulation self-consistently describes the dynamics of MBH pairs, one may simply use $T_{\text{Romulus}}$ in place of $T_{\text{df}}$. However, this would not be a robust prescription, as our modeling of the ROMULUS results via Equation (4) is purely stochastic. For instance, an MBH whose host galaxy has been tidally disrupted and for which a long $T_{\text{df}}$, given by Equation (3), would be in order, may be assigned a short $T_{\text{Romulus}}$ by sampling Equation (4). Therefore, to remain conservative on our merger rates, we use the maximum of these values.</footnote>
<table><caption>**Table 1**
Summary of All Models Explored in This Work</caption><thead><tr><th>Model</th><th>SN Feedback on MBH growth<br>(Habouzit et al. 2017)</th><th>Galaxy/BH Dynamical Friction<br>(Binney & Tremaine 2008)</th><th>Kiloparsec-scale Delays<br>(Tremmel et al. 2018)</th><th>≤100 pc Delays<br>(Dosopoulou & Antonini 2017)</th></tr></thead><tbody><tr><td>SN-delays</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>SN-delays-medium</td><td>✓</td><td>✓</td><td>×</td><td>✓</td></tr><tr><td>SN-delays-short</td><td>✓</td><td>×</td><td>×</td><td>✓</td></tr><tr><td>SN-nodelays</td><td>✓</td><td>×</td><td>×</td><td>×</td></tr><tr><td>noSN-delays</td><td>×</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>noSN-nodelays</td><td>×</td><td>×</td><td>×</td><td>×</td></tr></tbody></table>

and spins, as well as the properties of the host merging galaxies (e.g., their merger redshift, their masses, etc.). The merger rate per unit redshift is calculated by summing the contributions within each redshift bin $\Delta z$:

$$ \frac{d^2N}{dzdt} = \frac{4\pi c}{\Delta z} \sum_{N \in \Delta z} W(z) \left[ \frac{d_L(z)}{1+z} \right]^2, \quad (8) $$

where $d_L$ is the luminosity distance, and $W(z)$ is the comoving number density of the binary (obtained from the comoving number density of its host galaxy).

For each MBH binary, we then calculate the signal-to-noise ratio (S/N) $\rho$, averaged over the polarization, inclination, and sky position, i.e. (see, e.g., Cornish & Robson 2019):

$$ \rho^2 = 4 \int_{f_{\min}}^{f_{\max}} \frac{|h(f)|^2}{S_n(f)} df, \quad (9) $$

where $h(f)$ is the gravitational-wave strain amplitude in Fourier space, and $S_n(f)$ is the LISA sensitivity. We calculate the gravitational-wave strain using the PhenomC inspiral-merger- ringdown model (Santamaría et al. 2010), and we use the LISA sensitivity curve from Amaro-Seoane et al. (2017), including the contribution from the foreground from Galactic binaries (Cornish & Robson 2017, 2019). The time to merger is drawn from a uniform distribution between 0 and the nominal mission duration (4 yr). We use $\rho = 8$ as the detection threshold. The number of sources detectable in a 4 yr LISA mission is shown in Table 2 and will be discussed below.

## 4. LISA Detection Rates: Delay Mechanisms and SN Feedback

In this section, we discuss in detail the results of our various models for the predicted merger rates of MBHs and their potential detection with LISA. Our models are summarized in Table 1. Each model incorporates a different combination of prescriptions pertaining to the intermediate-scale dynamical evolution of MBHs (as described in Sections 2.1–2.3), as well as SN-regulated MBH growth (Section 2.4). Each evolutionary model is applied starting from either an LS or HS high-redshift population, as described in Section 2.

In Table 2, we report the total number of MBH mergers, as well as the number of detections expected in 4 yr of observations with LISA, for each model. For the noSN-delays model, we find the same general trend found in previous work

<table><caption>**Table 2**
Total Number of Sources and Detections Expected in 4 yr of Observation with LISA for All the Models Explored Here (See Table 1 for a Summary)</caption><thead><tr><th rowspan="2">Model</th><th colspan="2">LS</th><th colspan="2">HS</th></tr><tr><th>Total</th><th>Detected</th><th>Total</th><th>Detected</th></tr></thead><tbody><tr><td colspan="5" style="text-align: center;">SN feedback</td></tr><tr><td>SN-delays</td><td>48</td><td>16</td><td>25</td><td>25</td></tr><tr><td>SN-delays-medium</td><td>157</td><td>38</td><td>89</td><td>88</td></tr><tr><td>SN-delays-short</td><td>169</td><td>33</td><td>74</td><td>73</td></tr><tr><td>SN-nodelays</td><td>178</td><td>36</td><td>1269</td><td>1269</td></tr><tr><td colspan="5" style="text-align: center;">No SN feedback</td></tr><tr><td>noSN-delays</td><td>192</td><td>146</td><td>10</td><td>10</td></tr><tr><td>noSN-nodelays</td><td>1159</td><td>307</td><td>1288</td><td>1288</td></tr></tbody></table>

(Sesana et al. 2007; Klein et al. 2016; Bonetti et al. 2019), i.e., that common, low-mass seeds produce significantly more detected mergers than high-mass seeds. Interestingly, this trend is reversed for models accounting for the effect of SN feedback on MBH growth. Low-mass black holes are those whose merger and detection prospects are most affected by SN feedback, because their arrested growth yields not only longer binary formation timescales (see Equations (3) and (5)), but also lower S/Ns, as low-mass binaries tend to coalesce at the high-frequency end of the LISA sensitivity curve.

The impact of SN explosions on the HS models is more subtle. Without any intermediate-scale delays, SN feedback has little effect, as the MBHs are seeded at masses large enough to form binaries in a timely manner and in a mass range readily detectable with LISA. However, when delay times to MBH binary formation are included, SN feedback increases the merger rate. Without SN regulation, massive seeds experience more early growth, generally increasing their spins (Barausse 2012; Sesana et al. 2014). Merging MBHs therefore experience stronger recoil kicks, which tend to remove the merger remnant from the center of the host galaxy. This obviously prevents the remnant from coalescing with other MBHs brought in by future galaxy mergers. Without the delay mechanisms in place for MBH binary formation, MBHs merge before they have a chance to grow and substantially spin up, irrespective of whether SN feedback is present or not.
High-mass seed models are much more affected by the inclusion of MBH dynamical evolution. The inclusion of even one of the delay timescales discussed in Section 2 results in an order of magnitude decrease in both the total and detected MBH merger rate. More detailed choices of additional delay mechanisms have a reduced effect, but can still be important. In particular, the inclusion of kiloparsec-scale dynamical evolution of MBHs within galaxy merger remnants, as predicted from cosmological simulations (Tremmel et al. 2018), results in a factor of ~3 decrease in MBH mergers for HS models. LS models are much less sensitive to these intermediate-scale delays, but also see a decrease of a factor of ~2 when including these galaxy-scale delays.

In the following sections, we describe in more detail the effect of SN feedback and intermediate-scale dynamical delay timescales on the predicted distribution of MBH merger properties and their astrophysical implications. We note that even in our most pessimistic models, we predict several detectable mergers within the LISA nominal mission duration of 4 yr, similar to previous works (Klein et al. 2016; Bonetti et al. 2019; Dayal et al. 2019). Our prediction is higher than those derived from the Illustris cosmological simulation (Katz et al. 2019), though this is likely due to their limited resolution and lack of MBHs in low-mass dwarf galaxies (Volonteri et al. 2020).

## 4.1. Delay Mechanisms

We first focus on the effect of MBH binary formation delay timescales. Beyond standard prescriptions for dark matter halo sinking timescales, for MBH binary hardening/gas-driven migration, and for triple/quadruple MBH interactions, which our model already includes (Barausse 2012; Sesana et al. 2014; Antonini et al. 2015a; Bonetti et al. 2019), we incorporate three different intermediate-scale delay mechanisms in this work (see Section 2). These are to be added to the aforementioned delays already included in the model. As outlined in Table 1, the “delays-short” models account for the dynamical friction on bound MBH pairs beginning at the primary MBH’s sphere of influence (Dosopoulou & Antonini 2017). The “delays-medium” models additionally account for dynamical friction acting on merging galaxies (Binney & Tremaine 2008). Finally, the “delays” models further incorporate a third delay timescale associated with the dynamical evolution of MBHs within galaxies, which has been calibrated to cosmological simulations (Tremmel et al. 2018). In this section, all of the models include SN-regulated MBH growth. Note also that the “nodelays” model, which we include as a baseline, may be physically relevant to describe a scenario where stellar hardening is efficient at driving the evolution of MBH binaries already at separations of hundred of parsecs (Sesana et al. 2006; Khan et al. 2012, 2016).

Figure 1 shows the redshift distribution of the MBH mergers detectable during LISA’s 4 yr of nominal mission duration. The top panel plots the results for HS models. Without intermediate-scale delays, there are a large number of MBH mergers predicted, the majority of which occur at high redshift ($z > 5$). However, the inclusion of even one of these intermediate-scale delays (i.e., the SN-delays-short model; green line) shows a substantial decrease in mergers at $z > 5$. The addition of galaxy merger dynamical friction timescales (SN-delays-medium; orange line) has a small effect of shifting high-$z$ MBH mergers to lower redshift, though the difference between

<figure><img src="image_3.png" /><figcaption>Figure 2. Number of detected mergers as a function of total MBH–binary mass in the HS and LS models that include SN feedback (upper and lower panels, respectively). The mass distribution of LS models is barely unaffected by time delays. In HS models, low-mass (and high-redshift; see Figure 1) binaries acquire large delays due to dynamical friction.</figcaption></figure>

SN-delays-medium and SN-delays-short is much less significant that that between SN-nodelays and SN-delays-short. A more significant difference is seen when the delays from MBH dynamical evolution on kiloparsec scales is included (SN-delays; blue line), which results in an additional decrease in the merger rates, particularly at $z > 2$.

For low-mass seeds (bottom panel of Figure 1), only when kiloparsec-scale delays are included (SN-delays; blue line) is there an appreciable effect on the MBH merger rate. Regardless of the binary formation delay timescale, many of the MBH mergers from low-mass seeds will be missed by LISA (as can be seen by comparing to the total number of MBH mergers, shown by the dashed lines for two of the LS models). This is because, unlike the HS models, LSs result in many MBH mergers with total mass low enough to give S/Ns below our detection threshold ($\rho = 8$).

Figure 2 shows the total mass distribution of detected MBH mergers for both HS (top) and LS (bottom) models. The inclusion of intermediate-scale delays mostly affects the number of low-mass mergers, which typically take place as the result of low-mass, high-$z$ galaxy mergers. For HSs, this means that longer delay timescales result in fewer mergers of mass $10^4-10^6 M_{\odot}$. For LSs, the effect is less pronounced, but the decrease in the overall merger rate seen in Table 2 and Figure 1 is also due to an overall decrease of mergers with mass below $10^4 M_{\odot}$.

Figure 3 shows the cumulative distribution of MBH merger mass ratios $q = M_2/M_1 \le 1$ for events detectable with LISA. The distribution is largely unaffected by the MBH binary delay timescale for HSs, except for a slight steepening trend toward higher mass ratios when delays are included. The merger mass ratio for the LS models is much more affected by binary formation delays, but interestingly, our model incorporating the longest binary formation timescales (SN-delay) is more similar to the model without delays (SN-nodelays) than to the short and medium delay models. This is caused by the implementation of
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Fraction of detected mergers with mass ratio below *q* in the HS and LS models with SN feedback (upper and lower panels, respectively). The distribution in HS models slightly shifts toward larger *q* (more equal mass ratio) in models with large delays (note that about 10% of the mergers in the HS models have *q* ≲ 0.1). In contrast, including all time delays in LS models produces a distribution with a pronounced tail toward low *q*.</figcaption></figure>

dynamical friction on bound MBH pairs, which has an explicit
dependence on $M_2^{-1}$ (see Equation (5)): this will preferentially
delay low-mass-ratio MBH mergers, which have already longer
sinking timescales (see, e.g., Ogiya et al. 2020), resulting in
fewer mergers with low mass ratios. Dynamical friction between
galaxies delays galaxy mergers and allows the MBHs (and
particularly the primary) to grow further prior to binary
formation, resulting in more low-mass-ratio mergers. The
additional delay associated with kiloparsec-scale dynamical
evolution increases this effect until the distribution is nearly the
same as that without any delays.

4.2. SN Feedback

Next, we explore the effect of SN feedback (Habouzit et al. 2017) on the detection rates and properties of MBH binaries. We compare this effect to that of binary formation delay times by including both “delays” and “nodelays” models, each with and without SN feedback. In Figure 4, we show the redshift distribution of the predicted 4 yr LISA detections. As also evident from Table 2, the HS model is not strongly affected by SN feedback, while the low-mass seeds are significantly impacted (compare the pink and red curves in Figure 4).

The reason for this sensitivity of LS models to SN feedback
can be understood from Figure 5, which shows that the number
of detected mergers in the 10⁵–10⁷*M*☉ mass range, to which
LISA is most sensitive, is drastically reduced relative to models
without SN feedback. This deficit is inherited from the fact that
the total number of these binaries (i.e., before the S/N cutoff) is
reduced in LS models with SN feedback, because SN winds
tend to expel nuclear gas, which results in suppression of seed
growth (due to lack of accretion) and in longer timescales for
binary migration in circumbinary disks.

Next, in Figure 6, we compare the shape of the distribution
(normalized to total detections) of S/N values when the four
antipodal model prescriptions, delay/no-delay and SN/no-SN,

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Number of detected mergers per unit redshift during a 4 yr LISA mission in the HS and LS models (upper and lower panels, respectively; see Table 1 for a summary of the models). Note that the HS models with no delays (both with and without SN feedback) are almost indistinguishable.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Number of detected mergers as a function of total MBH mass in the HS and LS models (upper and lower panels, respectively). The mass distribution of LS models is practically unaffected by time delays. In HS models, low-mass (and high-redshift; see Figure 1) binaries are ejected from their host halos in models that include time delays.</figcaption></figure>

are considered. As in the previous figures, the upper panel
shows results from the HS scenario, while the bottom panel
concerns the LS case. From the figure, it can be seen that the
S/N distribution reflects the mass distribution of the sources
(Figure 5). As expected, the S/N distribution in LS models
peaks at lower S/N values due to the low source mass. The
shape of the mass distribution for LSs is not greatly affected by
dynamical delays (see Section 4.1) but is significantly affected
by SN feedback. As a result, when SN explosions are included,
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Distribution (normalized to the total number of detected events) of S/N in the HS and LS models (upper and lower panels, respectively). The distribution of S/N in the LS models is very mildly affected by dynamical effects that cause time delays, whereas in the HS case, dynamical effects on 100 pc scales affect mostly low-mass systems (see Figure 2) and/or systems with unequal mass ratios (see Figure 3). As a result, the “surviving” binaries have higher S/N in HS models with large time delays. Note also that the S/N in LS models without SN feedback shift toward slightly higher values relative to the case with SN feedback.</figcaption></figure>

the number of high-S/N sources in the LS scenario declines. For HS models, the S/N distribution is more affected by MBH dynamics, with longer MBH binary formation timescales resulting in more numerous high-S/N systems than low-S/N ones. This is due to the decrease in low-mass MBH mergers, to the delays allowing for the MBHs (and especially the primary) to grow prior to binary formation and merger, and to the mergers being delayed to lower redshift, where their S/N is naturally higher (see Figure 1). Note that the more massive MBHs in the HS model are much less affected by SN winds because of their high mass, which is already in the range to which LISA is most sensitive.

Given the steeply growing distribution toward lower S/Ns in the LS models, an important conclusion from Figure 6 is that, if LSs are indeed realized in the universe, the total number of LISA detections may be extremely sensitive to the S/N threshold or, equivalently, to the LISA noise budget (see also Bonetti et al. 2019). Therefore, in view of the small number of expected detections per year in the LS scenario (in which the inclusion of SN feedback is crucial), particular care should be put into optimizing LISA's sensitivity at high frequencies, as any degradation can translate into significant event losses. For instance, changing the single link optical measurement system noise from 10pm/√Hz, adopted in Amaro-Seoane et al. (2017), to 15pm/√Hz, in order to account for the margin on this noise contribution inserted in the LISA Science Requirements Document (ESA 2018), does not affect the detection rates for the HS models, but LS detection rates are reduced by a factor of ~2. To improve the sensitivity to low-mass systems, joint observation campaigns with space-based missions that are scheduled at the same time as LISA—e.g., TianQin (Luo et al. 2016)—and which are also sensitive to MBH binaries (Wang et al. 2019; Shi et al. 2019) may be particularly useful.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Number of detected MBH mergers per unit stellar mass of the host galaxy in the HS and LS models (upper and lower panels, respectively). MBH merger rates are suppressed in low-mass halos by SN winds and are everywhere suppressed by delays due to MBH dynamics on tens to thousands of parsec scales. Regardless of the model’s details, a consistent prediction is that the majority of MBH mergers detectable by LISA occur in dwarf galaxies of mass $M_* \lesssim 10^{9.5} M_\odot$.</figcaption></figure>

### 4.3. Astrophysical Implications

In this section, we examine our model’s predictions for the MBH binary population in terms of the properties of their host galaxies and their galaxy mergers. These predictions are crucial for interpreting LISA detections in the context of galaxy formation and evolution, so as to use LISA data to constrain models for the synergic coevolution of MBHs with galaxies (Sesana et al. 2007; Berti & Volonteri 2008; Volonteri & Natarajan 2009; Klein et al. 2016; Ricarte & Natarajan 2018a; Bonetti et al. 2019) and to attempt multiwavelength follow-up observations of LISA sources (Tamanini et al. 2016).

#### 4.3.1. Host Galaxies of MBH Mergers

In Figure 7, we examine how both delays and SN feedback affect the distribution of the host galaxy stellar mass. Once again, HS models are not strongly affected by SN feedback but are sensitive to the dynamics of MBH binary formation. The inclusion of delays to MBH formation decreases the number of detected mergers predicted by the model, without greatly affecting the shape of the host galaxy distribution. Conversely, the LS models are affected by both delays and SN feedback. The inclusion of delays (in particular the kiloparsec-scale dynamical evolution of MBH pairs) decreases the number of MBH mergers in higher mass galaxies, while SN feedback decreases the number of mergers at low masses ($M_* \lesssim 10^8 M_\odot$). The result is that models that include SN feedback have flatter galaxy mass distributions at the low-mass end. This is due to the fact that SN winds are efficient at removing gas from the shallow potential wells of low-mass galaxies, thus effectively shutting off accretion on LSs—which cannot reach the masses at which LISA is most sensitive (~$10^5-10^7 M_\odot$)—
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Distribution of delay times between galaxy and binary merger for the detected binaries, in the LS and HS models.</figcaption></figure>

and resulting also in longer binary migration timescales in circumbinary disks.

It is important to note that all models, despite their large differences in merger rates and mass distributions, predict that mergers in low-mass galaxies ($M_* \lesssim 10^{9.5} M_\odot$) dominate the MBH merger population that will be detected by LISA. This is in agreement with previous cosmological simulations (Volonteri et al. 2020). It is therefore critical that models that hope to make predictions relevant for LISA fully resolve high-redshift mergers between low-mass galaxies, a particularly difficult challenge for large-scale cosmological simulations.

### 4.3.2. *The (Dis)Connection between Galaxy and MBH Mergers*

Figure 8 shows the distribution of total delay times between the coalescence of detected MBH binaries and the merger of their progenitor galaxies. Unsurprisingly, the inclusion of delay times results in significantly longer timescales between galaxy and MBH mergers, similar to what has been predicted from cosmological simulations (Tremmel et al. 2018; Volonteri et al. 2020). This fact has critical implications for multimessenger astrophysics, as MBHs are likely to coalesce long after their host galaxies have merged. It is unlikely that by observing a clearly disturbed or an actively merging galaxy pair, one could identify the cradle of a LISA MBH binary (unless the MBH binary formed as a result of an even earlier galaxy merger event). Electromagnetic counterparts to MBH binaries may therefore be commonly found in relatively undisturbed galaxies, billions of years following the progenitor galaxy merger event.

The timescale for LS models is also affected by the presence of SN feedback. As discussed previously, SN winds remove gas in low-mass, high-redshift galaxies and therefore curtail MBH growth and gas-driven migration in those systems. This has a clear impact on Figure 8 (though not as dramatic as the inclusion of intermediate-scale delays). Indeed, by ejecting gas from the nuclear region, SN explosions result in the

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Bolometric quasar luminosity function at $z=6$ for HS and LS models (colored lines, upper and lower panels, respectively). Black line: observed bolometric luminosity function (Shen et al. 2020). Black arrows: upper limits on the X-ray luminosity function (Vito et al. 2016), transformed into bolometric luminosity.</figcaption></figure>

suppression of circumbinary disk migration, thereby increasing the typical delays between galaxy and MBH coalescences.

## 5. Quasar Luminosity Functions

While we have discussed the predictions and implications of the model as they relate to gravitational wave detection, we can also examine the model results against electromagnetic observations, namely the quasar luminosity function. We adopt the most recent estimate of the bolometric luminosity function (Shen et al. 2020) and supplement it with an upper limit to the faint end at $z=6$. The latter is determined in X-rays (Vito et al. 2016) but is transformed into bolometric luminosity, after applying a correction for Compton-thick AGNs (Ueda et al. 2014), by using the same bolometric correction as in Shen et al. (2020).

All of our models (both LS and HS, with and without delays, with and without SN feedback) reproduce equally well the quasar luminosity function at $z < 3$, with differences appearing only at higher redshift. We therefore focus our comparison at $z=6$ (Figure 9) and discuss lower redshifts only briefly.

Dynamical delays do not have a noticeable effect on the luminosity function. Inclusion of the effect of SN feedback on MBH growth has the strongest impact on LS models. The specific implementation used here for SN feedback causes a clear underestimate of the luminosity function at $z=6$ for all LS models that include this effect.<fnref n="8" /> LS models without SN feedback also struggle somewhat to produce sufficiently massive MBHs and bright quasars, but overestimate the upper limit at the faint end. In summary, LS models without SN feedback overestimate the number of faint quasars, while LS

<footnote n="8">⁸ Note that Habouzit et al. (2017) find relatively good agreement with the same upper limit starting with MBH seeds with mass $\sim 10^2-10^3 M_\odot$. This confirms again that the specific implementation of how SN feedback affects MBH growth has a significant bearing on the results.</footnote>
models with SN feedback underestimate the number of bright AGNs.

All HS models fare well with bright AGNs (note that the semianalytic model does not include halos with mass >10$^{12}$–10$^{13}$ at z = 6, where quasars above the knee of the luminosity function are expected to reside). HS models without SN feedback overestimate the faint end of the luminosity function, being above the upper limit derived by Vito et al. (2016). HS models with SN feedback are compatible with the upper limit at the faint end and sit nicely on the observed portion of the luminosity function.

An interesting point is that models with very different merger rates produce very similar luminosity functions. Table 2 shows that the HS models SN-delays and SN-nodelays have a merger rate that differs by almost two orders of magnitude (25 versus 1269), but the luminosity functions are very similar. This comes about because models featuring the same (large) number of MBH seeds at high redshift, but with different merger rates due to dynamical delays, can produce luminosity functions compatible with observations if SN feedback efficiently suppresses MBH accretion.

The results are similar for the various models at z = 5, while by z = 4 they start to converge to the same luminosity function, with LS and HS models faring equally well for bright quasars, and with only LS models without SN feedback slightly overpredicting the faint end down to z = 3. At even lower redshift, all models produce indistinguishable luminosity functions.

## 6. Stochastic Background in the Pulsar-timing Array Band

Observations of the unresolved stochastic background of gravitational waves by pulsar-timing array experiments are scientifically complementary to LISA operations. These experiments include the European Pulsar Timing Array (EPTA) collaboration (Desvignes et al. 2016), the Australian Parkes Pulsar Timing Array (PPTA) experiment (Reardon et al. 2016), and the American NANOGrav collaboration (The NANOGrav Collaboration et al. 2015). These experiments also share their data under the patronage of the International Pulsar Timing Array (IPTA) collaboration (Manchester & IPTA 2013; Verbiest et al. 2016; Perera et al. 2019).

Pulsar-timing arrays attempt to detect gravitational waves by cross-correlating the timing residuals of ms pulsars (Sazhin 1978; Foster & Backer 1990). The presence of a gravitational wave stochastic background at frequencies of ~nanohertz would produce a potentially detectable quadrupolar correlation between the residuals of pulsars at different sky locations (Hellings & Downs 1983). The same technique can also detect individual gravitational wave signals, if those are strong enough to be resolved above the stochastic background (Sesana et al. 2009; Babak et al. 2016; Arzoumanian et al. 2018).

The stochastic background at nanohertz frequencies is expected to be produced mainly by inspiraling binaries of MBHs with total masses between 10$^8$ and 10$^{10}$ M$_{\odot}$ at redshifts $\lesssim 2$ (Rajagopal & Romani 1995; Jaffe & Backer 2003; Wyithe & Loeb 2003; Sesana et al. 2008, 2009; Ravi et al. 2012; Sesana 2013; McWilliams et al. 2014; Kulier et al. 2015; Ravi et al. 2015; Sesana et al. 2016; Kelley et al. 2017; Bonetti et al. 2018b). While these systems are heavier than those targeted by LISA, their hierarchical formation mechanism is the same (halo/ galaxy mergers followed by the formation of MBH pairs/ binaries). Therefore, pulsar-timing arrays provide a distinct,

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Characteristic strain of the stochastic background from MBH binaries in the band of pulsar-timing array experiments. The red and blue shaded regions encompass the LS and HS models, respectively. Unlike for the merger rate in the LISA band, the predictions for the pulsar-timing array signal are quite robust and show only a mild dependence on the model. Also shown are the sensitivity curves of ongoing pulsar-timing array experiments (Shannon et al. 2015; Desvignes et al. 2016; Arzoumanian et al. 2018).</figcaption></figure>

complementary way of exploring astrophysics similar to LISA's, but at larger scales and in different host environments.

Interestingly, the aforementioned pulsar-timing array experiments are already ongoing and analyzing data, and have put strong upper bounds on the stochastic background in the nanohertz band (Lentati et al. 2015; Arzoumanian et al. 2016; Shannon et al. 2015; Verbiest et al. 2016; Perera et al. 2019). The most robust constraint to date comes from NANOGrav's 11 yr data set (Arzoumanian et al. 2018). Stronger bounds have been put forward by PPTA (Shannon et al. 2015), but it is unclear if those robustly account for uncertainties in the position of the solar system barycenter. Overall, these upper bounds have already put constraints on models of MBH mergers (Rajagopal & Romani 1995; Jaffe & Backer 2003; Wyithe & Loeb 2003; Sesana et al. 2008, 2009; Ravi et al. 2012; Sesana 2013; McWilliams et al. 2014; Kulier et al. 2015; Ravi et al. 2015; Sesana et al. 2016; Kelley et al. 2017; Bonetti et al. 2018b; Middleton et al. 2018), even ruling out the most extreme ones (McWilliams et al. 2014) in which such mergers are very abundant.

In Figure 10, we show the predictions of the models presented in this paper for the background's characteristic strain $h_c$, as a function of gravitational wave frequency and compared to the upper bounds from EPTA, PPTA, and NANOGrav (see, e.g., Phinney 2001; Sesana et al. 2008; Dvorkin & Barausse 2017 for the definition and computation of $h_c$). The pink and purple shaded areas denote the envelope of the predictions of all models (i.e., the minimum area enclosing the predictions of all models) for our LS and HS scenarios, respectively. The slope of the background follows from the assumption that circular binaries lose energy only through the emission of gravitational waves. Indeed, most of the background signal comes from MBH binaries in their early inspiral phase, where GW emission is well described by the quadrupole formula, which gives $h_c \propto f^{-2/3}$. Note that the predictions from all models are very similar, as the signal is mostly emitted by binaries involving MBHs with masses above 10$^8$ M$_{\odot}$ at low redshift. For these systems, the impact of the different seeding and delay prescriptions is minor.
Reassuringly, the scatter of the predictions of the different models is very small, unlike what happens for the LISA detection rate predictions shown above. This was to be expected (see, e.g., Bonetti et al. 2018b), because the dependence on the seeding mechanism fades out when MBHs evolve to very large masses, and because the bulk of the pulsar-timing array signal comes from comparable-mass MBH binaries, for which the delays between galaxy and MBH merger are generally shorter (Dvorkin & Barausse 2017). By comparing to the results of Bonetti et al. (2018b), Figure 10 suggests that the stochastic background may be detected by pulsar-timing arrays in ≈15–20 yr of data collection, assuming a putative array of 50 ms pulsars monitored at the 100 ns level of precision. A significantly earlier detection would be achievable with the Square Kilometer Array (SKA) telescope (Dvorkin & Barausse 2017).

## 7. Conclusions

In this paper, we explore how the expected event rates for MBH mergers and their corresponding gravitational wave signals depend on the physical processes delaying the evolution of MBHs at intermediate separations (~tens to thousands of parsec) and on processes like SN feedback, which can regulate their growth by accretion and affect their gas-driven migration. To this purpose, we perform semianalytic simulations of the coevolution of galaxies and MBHs, starting from either low- or high-mass seeds for the MBH population at high redshift. We find that, regardless of the MBH seed mass model, the predicted rates of MBH mergers are heavily dependent on the inclusion of both delayed binary formation and SN-regulated MBH growth.

The effect of including even moderate delay timescales (i.e., the inclusion of dynamical friction timescales on scales below 100 pc as in our delays-short models) results in significantly fewer MBH mergers for models incorporating high-mass MBH seeds. Merger rates from low-mass seed models are less sensitive to such delay timescales when SN-regulated growth is included. The predicted MBH merger rates are also relatively insensitive to the detailed model choices for binary formation timescales. However, the inclusion of galactic-scale dynamical evolution (i.e., from several kiloparsecs to ~100 pc separations), a phase of evolution often disregarded in semianalytic models, results in an additional factor of ~2–4 decrease in MBH merger and detection rates for both low- and high-mass seed models. The inclusion of these delay timescales mostly affects the abundance of high-redshift, low-mass ($M_{tot} < 10^6 M_{\odot}$) mergers in the case of high-mass MBH seeds. For models assuming low-mass MBH seeds, binary formation delay timescales also affect the tail of the mass ratio distribution of MBH mergers.

Feedback from SN winds primarily affects the merger rate of low-mass MBH seeds. In these models, SN feedback is able to expel gas from the center of low-mass galaxies/dark matter halos, thereby regulating the gas reservoir available to MBHs in the early universe. For low-mass seed models, on the one hand, this decreases the intrinsic merger rate of MBHs, as the depletion of gas results in longer binary hardening timescales. On the other hand, SN winds also suppress seed growth via accretion, preventing MBHs from growing to sufficiently high masses for their mergers to be detectable with LISA.

Models with high-mass MBH seeds are generally less affected by SN feedback due to their already substantial

masses, which typically are already detectable by LISA. However, when binary formation delay times are included, SN feedback results in slightly more numerous MBH mergers. This happens because rapid, early mass growth spins MBHs up if SN feedback is unaccounted for. As a result, when MBHs merge at high redshift, the MBH remnant forming from the merger experiences a strong gravitational recoil kick. The kick can even be high enough to remove the MBH remnant from the host galaxy, depleting the number of central MBHs available to experience future mergers.

Another important consequence of SN-regulated MBH growth is on the predicted luminosity function at high redshift. Both low- and high-mass seed models predict too many low-luminosity MBHs without SN feedback. The effect of SN feedback is less dramatic for high-mass MBH seeds, but still significant. In summary, only models that include SN-regulated MBH growth are consistent with high-z quasar observations.

The most pessimistic model that we use in this work, which includes the longest binary formation delay times as well as SN-regulated growth, predicts that LISA should be able to detect several MBH mergers during its nominal mission duration of 4 yr. This work highlights how LISA will be a critical tool for constraining and discriminating models of MBH growth and dynamical evolution. Conversely, we predict that pulsar-timing array detection of the gravitational wave stochastic background should be relatively insensitive to model variations, including, unsurprisingly, the MBH seed masses. Rather, our results robustly show that the stochastic background should be detectable in the near future as the sensitivity of pulsar-timing arrays improves (Bonetti et al. 2018b; Dvorkin & Barausse 2017) and that this remains true even for dramatically different MBH evolutionary models.

We acknowledge the financial support provided under the European Union’s H2020 ERC Consolidator Grant “GRavity from Astrophysical to Microscopic Scales” grant agreement No. GRAMS-815673. I.D. thanks SISSA for hospitality during the early stages of this work. This work has made use of the Horizon Cluster, hosted by the Institut d’Astrophysique de Paris. We thank Stephane Rouberol for smoothly running this cluster for us.

## ORCID iDs

Enrico Barausse https://orcid.org/0000-0001-6499-6263

Michael Tremmel https://orcid.org/0000-0002-4353-0306

Marta Volonteri https://orcid.org/0000-0002-3216-1322

Matteo Bonetti https://orcid.org/0000-0001-7889-6810

## References

* Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, arXiv:1702.00786

* Antonini, F., Barausse, E., & Silk, J. 2015a, ApJ, 812, 72

* Antonini, F., Barausse, E., & Silk, J. 2015b, ApJL, 806, L8

* Arzoumanian, Z., Baker, P. T., Brazier, A., et al. 2018, ApJ, 859, 47

* Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2016, ApJ, 821, 13

* Babak, S., Petiteau, A., Sesana, A., et al. 2016, MNRAS, 455, 1665

* Baldassare, V. F., Geha, M., & Greene, J. 2020, ApJ, 896, 10

* Barausse, E. 2012, MNRAS, 423, 2533

* Barausse, E., Shankar, F., Bernardi, M., Dubois, Y., & Sheth, R. K. 2017, MNRAS, 468, 4782

* Berti, E., & Volonteri, M. 2008, ApJ, 684, 822

* Bieri, R., Dubois, Y., Rosdahl, J., et al. 2017, MNRAS, 464, 1854

* Binney, J., & Tremaine, S. 2008, Galactic Dynamics (2nd ed.; Princeton: Princeton Univ. Press)

* Blecha, L., Sijacki, D., Kelley, L. Z., et al. 2016, MNRAS, 456, 961

* Bonetti, M., Haardt, F., Sesana, A., & Barausse, E. 2018a, MNRAS, 477, 3910
Bonetti, M., Sesana, A., Barausse, E., & Haardt, F. 2018b, MNRAS, 477, 2599

Bonetti, M., Sesana, A., Haardt, F., Barausse, E., & Colpi, M. 2019, MNRAS, 486, 4044

Bortolas, E., Capelo, P. R., Zana, T., et al. 2020, MNRAS, 498, 3601

Boylan-Kolchin, M., Ma, C.-P., & Quataert, E. 2008, MNRAS, 383, 93

Cattaneo, A., Dekel, A., Devriendt, J., Guiderdoni, B., & Blaizot, J. 2006, MNRAS, 370, 1651

Colpi, M., Holley-Bockelmann, K., Bogdanović, T., et al. 2019, BAAS, 51, 432

Cornish, N., & Robson, T. 2017, JPhCS, 840, 012024

Cornish, N., & Robson, T. 2019, CQGra, 36, 105011

Correa, C. A., Schaye, J., Wyithe, J. S. B., et al. 2018, MNRAS, 473, 538

Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11

Cuadra, J., Armitage, P. J., Alexander, R. D., & Begelman, M. C. 2009, MNRAS, 393, 1423

Dayal, P., Rossi, E. M., Shiralilou, B., et al. 2019, MNRAS, 486, 2336

Dekel, A., & Birnboim, Y. 2006, MNRAS, 368, 2

Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Natur, 457, 451

Desvignes, G., Caballero, R. N., Lentati, L., et al. 2016, MNRAS, 458, 3341

Di Matteo, T., Springel, V., & Hernquist, L. 2005, Natur, 433, 604

Dickey, C., Geha, M., Wetzel, A., & El-Badry, K. 2019, ApJ, 884, 180

Dosopoulou, F., & Antonini, F. 2017, ApJ, 840, 31

Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2012, MNRAS, 420, 2662

Dubois, Y., Volonteri, M., Silk, J., et al. 2015, MNRAS, 452, 1502

Duffel, P. C., D'Orazio, D., Derdzinski, A., et al. 2020, ApJ, 901, 25

Dutton, A. A., & van den Bosch, F. C. 2009, MNRAS, 396, 141

Dvorkin, I., & Barausse, E. 2017, MNRAS, 470, 4547

ESA 2018, LISA Science Requirements Document, https://dms.cosmos.esa.int/COSMOS/doc_fetch.php?id=3752747

Foster, R. S., & Backer, D. C. 1990, ApJ, 361, 300

Gehren, T., Fried, J., Wehinger, P. A., & Wyckoff, S. 1984, ApJ, 278, 11

Granato, G. L., de Zotti, G., Silva, L., Bressan, A., & Danese, L. 2004, ApJ, 600, 580

Greene, J. E., Trader, J., & Ho, L. C. 2019, arXiv:1911.09678

Habouzit, M., Volonteri, M., & Dubois, Y. 2017, MNRAS, 468, 3935

Heger, A., & Woosley, S. E. 2002, ApJ, 567, 532

Hellings, R. W., & Downs, G. S. 1983, ApJL, 265, L39

Hirschmann, M., Dolag, K., Saro, A., et al. 2014, MNRAS, 442, 2304

Holley-Bockelmann, K., & Khan, F. M. 2015, ApJ, 810, 139

Jaffe, A. H., & Backer, D. C. 2003, ApJ, 583, 616

Katz, M. L., Kelley, L. Z., Dosopoulou, F., et al. 2019, MNRAS, 491, 2301

Kelley, L. Z., Blecha, L., Hernquist, L., Sesana, A., & Taylor, S. R. 2017, MNRAS, 471, 4508

Khan, F. M., Berentzen, I., Berczik, P., et al. 2012, ApJ, 756, 30

Khan, F. M., Fiacconi, D., Mayer, L., Berczik, P., & Just, A. 2016, ApJ, 828, 73

Khan, F. M., Holley-Bockelmann, K., Berczik, P., & Just, A. 2013, ApJ, 773, 100

Klein, A., Barausse, E., Sesana, A., et al. 2016, PhRvD, 93, 024003

Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511

Kormendy, J., & Richstone, D. 1995, ARA&A, 33, 581

Kozai, Y. 1962, AJ, 67, 591

Kulier, A., Ostriker, J. P., Natarajan, P., Lackner, C. N., & Cen, R. 2015, ApJ, 799, 178

Lapi, A., Raimundo, S., Aversa, R., et al. 2014, ApJ, 782, 69

Latif, M. A., & Ferrara, A. 2016, PASA, 33, e051

Lentati, L., Taylor, S. R., Mingarelli, C. M. F., et al. 2015, MNRAS, 453, 2576

Lidov, M. L. 1962, P&SS, 9, 719

Lodato, G., Nayakshin, S., King, A. R., & Pringle, J. E. 2009, MNRAS, 398, 1392

Luo, J., Chen, L.-S., Duan, H.-Z., et al. 2016, CQGra, 33, 035010

MacFadyen, A. I., & Milosavljević, M. 2008, ApJ, 672, 83

Madau, P., Haardt, F., & Dotti, M. 2014, ApJL, 784, L38

Madau, P., & Rees, M. J. 2001, ApJL, 551, L27

Manchester, R. N., & IPTA 2013, CQGra, 30, 224010

McConnell, N. J., & Ma, C.-P. 2013, ApJ, 764, 184

McWilliams, S. T., Ostriker, J. P., & Pretorius, F. 2014, ApJ, 789, 156

Middleton, H., Chen, S., del Pozzo, W., Sesana, A., & Vecchio, A. 2018, NatCo, 9, 573

Muñoz, D. J., Miranda, R., & Lai, D. 2019, ApJ, 871, 84

Nelson, D., Pillepich, A., Springel, V., et al. 2019, MNRAS, 490, 3234

Nixon, C. J., Cossins, P. J., King, A. R., & Pringle, J. E. 2011, MNRAS, 412, 1591

Ogiya, G., Hahn, O., Mingarelli, C. M. F., & Volonteri, M. 2020, MNRAS, 493, 3676

Parkinson, H., Cole, S., & Helly, J. 2008, MNRAS, 383, 557

Perera, B. B. P., DeCesar, M. E., Demorest, P. B., et al. 2019, MNRAS, 490, 4666

Pfister, H., Volonteri, M., Dubois, Y., Dotti, M., & Colpi, M. 2019, MNRAS, 486, 101

Phinney, E. S. 2001, arXiv:astro-ph/0108028

Plowman, J. E., Hellings, R. W., & Tsuruta, S. 2011, MNRAS, 415, 333

Pontzen, A., Tremmel, M., Roth, N., et al. 2017, MNRAS, 465, 547

Press, W. H., & Schechter, P. 1974, ApJ, 187, 425

Quinlan, G. D. 1996, NewA, 1, 35

Rajagopal, M., & Romani, R. W. 1995, ApJ, 446, 543

Ravi, V., Wyithe, J. S. B., Hobbs, G., et al. 2012, ApJ, 761, 84

Ravi, V., Wyithe, J. S. B., Shannon, R. M., & Hobbs, G. 2015, MNRAS, 447, 2772

Reardon, D. J., Hobbs, G., Coles, W., et al. 2016, MNRAS, 455, 1751

Reines, A. E., Greene, J. E., & Geha, M. 2013, ApJ, 775, 116

Reines, A. E., Sivakoff, G. R., Johnson, K. E., & Brogan, C. L. 2011, Natur, 470, 66

Ricarte, A., & Natarajan, P. 2018a, MNRAS, 481, 3278

Ricarte, A., & Natarajan, P. 2018b, MNRAS, 474, 1995

Ricarte, A., Tremmel, M., Natarajan, P., & Quinn, T. 2019, MNRAS, 489, 802

Roedig, C., Dotti, M., Sesana, A., Cuadra, J., & Colpi, M. 2011, MNRAS, 415, 3033

Santamaría, L., Ohme, F., Ajith, P., et al. 2010, PhRvD, 82, 064016

Sazhin, M. V. 1978, SvA, 22, 36

Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521

Schramm, M., & Silverman, J. D. 2013, ApJ, 767, 13

Sesana, A. 2013, MNRAS, 433, L1

Sesana, A., Barausse, E., Dotti, M., & Rossi, E. M. 2014, ApJ, 794, 104

Sesana, A., Gair, J., Berti, E., & Volonteri, M. 2011, PhRvD, 83, 044036

Sesana, A., Haardt, F., & Madau, P. 2006, ApJ, 651, 392

Sesana, A., Haardt, F., Madau, P., & Volonteri, M. 2004, ApJ, 611, 623

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

Sesana, A., Vecchio, A., & Volonteri, M. 2009, MNRAS, 394, 2255

Sesana, A., Volonteri, M., & Haardt, F. 2007, MNRAS, 377, 1711

Shankar, F., Bernardi, M., Sheth, R. K., et al. 2016, MNRAS, 460, 3119

Shannon, R. M., Ravi, V., Lentati, L. T., et al. 2015, Sci, 349, 1522

Sharma, R., Brooks, A., Somerville, R. S., et al. 2020, ApJ, 897, 103

Shen, X., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2020, MNRAS, 495, 3252

Shi, C., Bao, J., Wang, H.-T., et al. 2019, PhRvD, 100, 044036

Somerville, R. S., Hopkins, P. F., Cox, T. J., Robertson, B. E., & Hernquist, L. 2008, MNRAS, 391, 481

Somerville, R. S., & Primack, J. R. 1999, MNRAS, 310, 1087

Taffoni, G., Mayer, L., Colpi, M., & Governato, F. 2003, MNRAS, 341, 434

Tamanini, N., Caprini, C., Barausse, E., et al. 2016, JCAP, 4, 002

Tamburello, V., Capelo, P. R., Mayer, L., Bellovary, J. M., & Wadsley, J. W. 2017, MNRAS, 464, 2952

The NANOGrav Collaboration, Arzoumanian, Z., Brazier, A., et al. 2015, ApJ, 813, 65

Tremmel, M., Governato, F., Volonteri, M., Pontzen, A., & Quinn, T. R. 2018, ApJL, 857, L22

Tremmel, M., Governato, F., Volonteri, M., & Quinn, T. R. 2015, MNRAS, 451, 1868

Tremmel, M., Governato, F., Volonteri, M., Quinn, T. R., & Pontzen, A. 2018, MNRAS, 475, 4967

Tremmel, M., Karcher, M., Governato, F., et al. 2017, MNRAS, 470, 1121

Ueda, Y., Akiyama, M., Hasinger, G., Miyaji, T., & Watson, M. G. 2014, ApJ, 786, 104

Vasiliev, E., Antonini, F., & Merritt, D. 2015, ApJ, 810, 49

Verbiest, J. P. W., Lentati, L., Hobbs, G., et al. 2016, MNRAS, 458, 1267

Vito, F., Gilli, R., Vignali, C., et al. 2016, MNRAS, 463, 348

Vogelsberger, M., Genel, S., Springel, V., et al. 2014, Natur, 509, 177

Volonteria, M., Dubois, Y., Pichon, C., & Devriendt, J. 2016, MNRAS, 460, 2979

Volonteri, M., Haardt, F., & Madau, P. 2003, ApJ, 582, 559

Volonteri, M., Lodato, G., & Natarajan, P. 2008, MNRAS, 383, 1079

Volonteri, M., & Natarajan, P. 2009, MNRAS, 400, 1911

Volonteri, M., Pfister, H., Beckman, R. S., et al. 2020, MNRAS, 498, 2219

Volonteri, M., & Reines, A. E. 2016, ApJL, 820, L6

Wang, H.-T., Jiang, Z., Sesana, A., et al. 2019, PhRvD, 100, 043003

Wyithe, J. S. B., & Loeb, A. 2003, ApJ, 590, 691